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

We investigate the formation of molecular clouds from atomic gas by using three-dimensional magnetohydrodynamical simulations,including non-equilibrium chemical reactions, heating/cooling processes, and self-gravity by changing the collision speed $V_0$ and the angle $\theta$ between the magnetic field and colliding flow. We found that the efficiency of the dense gas formation depends on $\theta$. For small $\theta$, anisotropic super-Alfv\'enic turbulence delays the formation of gravitationally unstable clumps. An increase in $\theta$ develops shock-amplified magnetic fields along which the gas is accumulated, making prominent filamentary structures. We further investigate the statistical properties of dense clumps identified with different density thresholds. The statistical properties of the dense clumps with lower densities depend on $V_0$ and $\theta$ because their properties are inherited from the global turbulence structure of the molecular clouds. By contrast, denser clumps appear to have asymptotic universal statistical properties, which do not depend on the properties of the colliding flow significantly. The internal velocity dispersions approach subsonic and plasma $\beta$ becomes order of unity. We develop an analytic formula of the virial parameter which reproduces the simulation results reasonably well. This property may be one of the reasons for the universality of the initial mass function of stars.


INTRODUCTION
The formation and evolution of molecular clouds (MCs) are crucial to set the initial condition of star formation. The formation of MCs arises from accumulation of the diffuse atomic gas, which consists of the cold neutral medium (CNM) and warm neutral medium (WNM) (Field et al. 1969).
The formation of MCs has been intensively studied by numerical simulations of atomic colliding flows (e.g., Heiles & Troland 2005;Vázquez-Semadeni et al. 2006;Hennebelle et al. 2008;Inoue & Inutsuka 2012). The colliding flows model accumulation of the atomic gas owing to supernova explosions (e.g., McKee & Ostriker 1977), expansion of super-bubbles (e.g., McCray & Kafatos 1987), galactic spiral waves (e.g., Wada et al. 2011), and large-scale colliding flows. The shock compression induces the thermal instability (Field 1965) to form cold clouds (Hennebelle & Pérault 2000;Koyama & Inutsuka 2000). At the same time, a part of the kinetic energy of the colliding flows is converted into the turbulent energy of the cold clouds and intercloud gas, generating supersonic translational velocity dispersions of the cold clouds (Koyama & Inutsuka 2002). MCs form if a sufficient amount of the gas is accumulated into the cold clouds.
Magnetic fields play a crucial role on the physical properties of MCs. Iwasaki et al. (2019) (hereafter Paper I) investigated the early stage of the molecular cloud formation through compression of the two-phase atomic gas with various compression conditions (also see Inoue & Inutsuka 2009;Heitsch et al. 2009;Körtgen & Banerjee 2015;Dobbs & Wurster 2021). Paper I found that the physical properties of the post-shock layers strongly depend on the angle θ between the magnetic field and collision velocity. For a sufficiently small θ, accumulation of the two-phase atomic gas drives super-Alfvénic turbulence (Inoue & Inutsuka 2012). As θ increases, the Alfvén mach number of the turbulence decreases, and the shock-amplified magnetic field plays a major role in the dynamics of the post-shock layer.
In this paper, we investigate how the θ-dependence of the turbulent properties affects the subsequent star formation by using the numerical simulations including self-gravity. We focus on the physical properties of dense clumps that are often discussed by using the virial theorem both observationally (Bertoldi & McKee 1992) and theoretically (McKee & Zweibel 1992). In most colliding-flow simulations, the virial theorem of dense clumps have not been discussed quantitatively. For instance, Inoue & Inutsuka (2012) investigated the physi-cal properties of dense clumps and discuss their stability using the virial theorem, but their simulations did not include self-gravity and they did not provide quantitative discussions. We will derive an analytic estimate of each term in the virial theorem on the basis of the results of the numerical simulations.
This paper is structured as follows: In Section 2, we describe the methods and models adopted in this paper. We present our results including the global physical properties of the post-shock layers and the statistical properties of dense clumps in Section 3. Discussions are made in Section 5. Finally, we summarize our results in Section 6.

Methods
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 ρ is the gas density, v i is the gas velocity, B i is the magnetic field, E = ρv 2 /2+P/(γ −1)+B 2 /8π is the total energy density, φ is the gravitational potential, T is the gas temperature, κ is the thermal conductivity, and L is a net cooling rate per unit volume. The thermal conductivity for neutral hydrogen is given by κ(T ) = 2.5 × 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 × 1024 2 , which is two times higher than that in Paper I. The corresponding mesh size is ∆x ∼ 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 λ J ≡ c s π/Gρ is shorter than 4∆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 max ≡ 1.4 × 10 5 cm −3 (N x /2048) 2 (c s /0.3 km s −1 ) −2 above which the gas is artificially treated as the adiabatic fluid with γ = 5/3.
The maximum density n 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 km s −1 and V 0 = 10 km s −1 . The cold phase of the atomic gas is as dense as n cold ∼ 50 cm −3 . The maximum densities obtained by the shock compression are estimated to be n cold (V 0 /0.3 km s −1 ) 2 ∼ 10 4 cm −3 for V 0 = 5 km s −1 and ∼ 4 × 10 4 cm −3 for V 0 = 10 km s −1 , both of which are smaller than n max .
Paper I took into account 9 chemical species (H, H 2 , 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 2 , 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 farultraviolet (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 ξ H varies significanly depending on regions. In this paper, we adopt the recently updated one in diffuse molecular clouds ξ H = 2×10 −16 s −1 (Indriolo et al. 2007(Indriolo et al. , 2015Neufeld & 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 ξ H decreases with A V in proportional to A −1 V for A V > 0.5. Our adopted value of ξ H = 2 × 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 4thorder Runge-Kutta scheme (Press et al. 1986). The time width ∆t 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.

A Brief Review of Paper I
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 ( n 0 , V 0 , B 0 , θ), where n 0 is the mean density of the pre-shock atomic gas, V 0 is the collision velocity, and B 0 is the field strength, and θ is the angle between the collision velocity and magnetic field of the atomic gas. At a given parameter set of ( n 0 , B 0 , V 0 ), they found that there is a critical angle θ cr which characterizes the physical properties of the post-shock layers. They obtained an analytic formula of the critical angle as follows: sin θ cr = 0.1 n 0 10 cm −3 • For an angle smaller than θ 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 θ cr , the shockamplified magnetic field weakens the post-shock turbulence, making the post-shock layer denser. The mean post-shock density has a peak around θ ∼ θ cr as a function of θ.
• For an angle larger than θ cr , the magnetic pressure plays an important role in the post-shock layer. Their mean densities decrease with θ and are determined by the balance between the post-shock magnetic pressure and pre-shock ram pressure (see also Inoue & Inutsuka 2012).

Models
We consider colliding flows of the two-phase atomic gas whose mean density is n 0 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 θ. The colliding flows move in the ±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 twophase 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 n 0 and a uniform magnetic field of B 0 in a cubic simulation box of −10 pc ≤ x, y, z ≤ 10 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 km s −1 (L/ 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 ∼ 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 pc ≤ x ≤ 20pc, −10 pc ≤ y, z ≤ 10 pc. In the y-and z-directions, periodic boundary conditions are imposed. The x-boundary conditions at x = ±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, n 0 = 10 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 ∼ 10 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 µ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 km s −1 . In order to investigate the dependence of the results on the collision velocity, we perform additional simulations with V 0 = 10 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 θ cr = 0.1 for V 0 = 5 km s −1 and 0.2 for V 0 = 10 km s −1 (Equation (7)).
For each V 0 , we consider cases with θ = 0.25θ cr and θ = 2θ cr as representative models of the turbulencedominated and magnetic-pressure-dominated layers, respectively (Section 2.2). We name a model by attaching  Table 1. List of the model parameters. In the rightmost column, t first indicates the epoch when the first unstable clump is formed. In all the models, n0 = 10 cm −3 and B0 = 3 µG are adopted.
"Θ cr " in front of a value of θ/θ cr followed by "V" in front of a value of V 0 /1 km s −1 (Table 1). It is useful to estimate global mass-to-flux ratios µ 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 layer = 2 ρ 0 V 0 L 2 t, where L = 20 pc is the box size along the y-and z−axes. The magnetic flux held in the postshock layer depends on direction. In the plane perpendicular to the x-axis (y-axis), the magnetic fluxes The global mass-to-flux ratio is estimated as µ layer = min(M layer /Φ x , M layer /Φ y ). It increases with time as M layer /Φ x ∝ t, but does not exceed M layer /Φ y = ρ 0 L/(B 0 sin θ). The global gravitational instability of the post-shock layers occurs if µ layer exceeds a critical value of µ cr = (2π √ G) −1 . The ratio is given by where n 0 1 = n 0 /10 cm −3 , B 3 = B 0 /3 µG, V 5 = V 0 /5 km s −1 , L 20 = L/20 pc, and t 10 = t/10 Myr. The mass-to-flux ratios of the post-shock layers increase in proportion to time and exceed unity at t = 2.5V −1

5
Myr. They continue to increase at least until t = 10 Myr for the θ = 0.25θ cr models while they reach the maximum values of 2 at t = 5 Myr for model Θ cr 2V5 and 4 at t = 2.5 Myr for model Θ cr 2V10. The star formation will be suppressed for largely-inclined compressions with sin θ > 0.8 n 0 10 B 3 L 20 because M layer /Φ y < µ cr unless magnetic dissipation processes, such as ambipolar diffusion and/or turbulent reconnection work.

RESULTS
In this section, we show the simulation results up to t = t first + 0.5 Myr, where t first is the formation time of the first unstable clump, which will be defined in Section 3.2. In our simulations, the resolution is not high enough to resolve dense cores (> 10 5 cm −3 ) and the thermal processes are switched-off in the regions where the Jeans condition is violated (Truelove et al. 1997) instead of inserting sink particles. Therefore, we focus on the early evolution of the clumps rather than performing longterm simulations in this paper. Column density maps integrated along the zaxis at t = 5 Myr for models (a) Θcr0.25V5, (b) Θcr2V5, (c) Θcr0.25V10, and (d) Θcr2V10. The black arrows indicate the direction of the density-weighted mean magnetic field along the z-axis.

Global Properties of post-shock Layers
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 θ = 0.25θ 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 kin , E mag , E grv of the post-shock layer 1 for θ = 0.25θ cr . We first focus on the results with V 0 = 5 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. (Top panels) Time evolutions of (red) kinetic energy, (green) magnetic energy, and (blue) gravitational energy densities of the post-shock layers are compared between the models with (a) θ = 0.25θcr and (b) θ = 2θcr. All the energy densities are divided by V 2 0 . The gray line correspond to Egrv,ana shown in Equation (10). The horizontal black lines indicate the kinetic energy density flowing into the postshock layer from the atomic gas, ρ0 V 2 0 /2. (Bottom panels) Time evolutions of the density-weighted velocity dispersions in the (red) x-, (green) y-, and (blue) z-directions for (c) θ = 0.25θcr and (d) θ = 2θcr. In Panel (d), the two black lines correspond to the predictions from δvx/V0 = 0.45τ −1/2 , where τ = t n0 /5 cm −3 V0/20 km s −1 /(1 Myr) (Paper I) with V0 = 5 km s −1 and 10 km s −1 . In all the panels, the solid and dashed lines show the results with V0 = 5 km s −1 and 10 km s −1 , respectively.
An increase in V 0 from 5 km s −1 to 10 km s −1 does not change the physical properties of the post-shock layers for θ = 0.25θ cr . Figure 2a clearly shows that all the energy densities are roughly proportional to V 2 0 , 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 2 0 . Figure 2c shows that δ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 δv x / (δv 2 y + δv 2 z )/2, is as large as ∼ 3 for V 0 = 5 km s −1 and ∼ 4 for V 0 = 10 km s −1 at t = 5 Myr. In addition, δv z is almost identical to δv y for each model. This suggests that there is no preferential directions in the (y, z)-plane. When θ is increased from 0.25θ cr to 2θ cr for V 0 = 5 km s −1 (model Θ cr 2V5), 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 θ = 2θ cr because the Lorentz force owing to the shock-amplified magnetic field decelerates the CNM clumps efficiently for θ = 2θ 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 δv x for model Θ cr 2V5 is as large as that for model Θ cr 0.25V5 because the magnetic field has not been amplified yet. In the subsequent evolution, δv x decreases in proportion to t −1/2 (Paper I), and all the three components δv x,y,z become of similar magnitudes around t ∼ 5 Myr.
Unlike model Θ cr 0.25V5, there is anisotropy in the velocity dispersion in the (y, z) plane for model Θ cr 2V5. The velocity dispersion along the y-axis increases with time while δv z is constant with time. The continuous increase in δv y corresponds to the formation of filamentary structure by self-gravity, which will be shown in Figure 4.
Unlike the models with θ = 0.25θ cr , the time evolution of the physical quantities depends on the collision speed for the θ = 2θ cr models. The relative importance of E mag to the total energy is more significant for V 0 = 10 km s −1 than for V 0 = 5 km s −1 . After the early evolution (t < 1 Myr) in which the velocity dispersions δv x,y,z increase in proportion to V 0 , δv x begins to decrease earlier for V 0 = 10 km s −1 than for V 0 = 5 km s −1 . This behavior is consistent with the results of Paper I where they found that δv The formula suggests that δv x /V 0 is independent of ρ 0 and V 0 if δv x is measured at the time of reaching the same column density. The effect of the rapid decrease in δv x can be seen in the fact that the shock fronts are flatter for V 0 = 10 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 km s −1 . Interestingly, the epochs (t ∼ 5 Myr) when δv x becomes equal to δv y ap-pear to be independent of V 0 .
There are two main mechanisms to drive the postshock 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 n 0 V 2 0 , 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 ∼ 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 δn 0 / n 0 is less than 10%, where δn 0 is the upstream density perturbation. In such situations, the energy conversion efficiency from the upstream kinetic energy ρ 0 V 2 0 /2 to the post-shock turbulence energy is as low as a few percents. When δn 0 /n 0 > 10 %, the interaction between shock fronts and density inhomogeneity drives strong turbulence where is as high as 10%, which is roughly consistent with the values obtained in this paper, which is ∼ (δv x /V 0 ) 2 ∼ 16% (Paper I).
3.2. Structure Formation due to Self-gravity where the integration is done over the post-shock layer 3 . The continuous accumulation of the atomic gas deepens the gravitational potential well, and −E 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 Σ 0 (t) = 2 ρ 0 V 0 t is the mean column density of the post-shock layer. Figures 2a and 2b show that the time evolution of E grv is consistent with that of E grv,ana for all the models, indicating that the gravitational energy densities increase in proportion to V 2 0 . The deviations of E grv from E grv,ana come from the non-uniform density distributions. Time evolution of the mass fraction of dense gas with n > nmax for models (solid red) Θcr0.25V5, (solid blue) Θcr2V5, (dashed red) Θcr0.25V10, and (dashed blue) Θcr2V10.
As an indicator of star formation, we measure the mass fraction of the dense gas whose density is larger than n 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 max does not depend on V 0 significantly for each θ ( 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 n 0 accelerates the formation of the dense gas with n > n max because the gravitational energy density is proportional to n 0 2 while the ram pressure of the preshock gas is proportional to n 0 .
As is expected in Paper I, the formation of the dense gas with n > n max is delayed by compressions with smaller θ because anisotropic super-Alfénic turbulence is driven for θ less than θ cr .
Another interesting feature is that the difference in V 0 affects the formation of the dense gas with n > n max in a different way at different θ. As V 0 increases, the time evolution of the dense gas fraction does not change significantly for θ = 2θ cr while the formation of the dense gas is suppressed slightly for θ = 0.25θ cr . This is due to the fact that the θ-dependence of the velocity field differs depending on the collision velocity. For θ = 0.25θ cr , the anisotropy of the velocity dispersion becomes more significant with increasing V 0 ( Figure 2c).
We define t 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 max becomes negative, where the clumps are identified by a friends-of-friends method. Table 1 lists t first , which are almost identical to the epochs when the mass with n > n max increases rapidly as shown in Figure 3. Color maps of the column densities integrated along the x-axis for models (a)Θcr0.25V5, (b)Θcr2V5, (c)Θcr0.25V10, and (d)Θcr2V10 at t first + 0.5 Myr. In each panel, the black arrows indicate the direction of the magnetic fields.
The column density distributions along the collision direction strongly depend on θ as shown in Figure 4. For the models with θ = 2θ 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 selfgravity preferentially along with the magnetic fields. By contrast, the column density maps for the models with θ = 0.25θ cr do not show filamentary structures clearly. The projected magnetic fields are randomized. The detailed structure of filamentary clouds is investigated in forthcoming papers.

Field strength-density Relation
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 first + 0.5 Myr for the models with V 0 = 5 km s −1 and V 0 = 10 km s −1 , respectively.
We first consider the models with V 0 = 5 km s −1 (Figure 5a). In low densities less than 10 3 cm −3 , model Θ cr 2V5 shows that the mean field strength is determined by the shock compression that amplifies the field strength up to (11) which is derived by balance between the pre-shock ram pressure and post-shock magnetic pressure. The independence of B ∼ B sh on n for model Θ cr 2V5 reflects gas condensation along the magnetic field, corresponding to the formation of filaments. By contrast, for model Θ cr 0.25V5, the mean field strength (∼ 6 µ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 is 20%, the following estimated field strength B tur is consistent with the simulation results, (12) Once the gas density becomes larger than ∼ 10 3 cm −3 , the θ-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 β=1 , which is given by β = 1 at T = 20 K that is a typical temperature of the dense regions, where β 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 km s −1 behave similarly to those of the models with V 0 = 5 km s −1 . At each θ/θ 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 km s −1 , in the high density regions, the field strengths follow B β=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 km s −1 to 10 km s −1 .
For the models with θ = 2θ cr , the mean field strengths take a constant value of B sh for lower densities and it in-creases following B β=1 ∝ n 1/2 for higher densities (Figure 5b). From these facts, a critical density n grv can be derived by equating B sh and B β=1 as follows: When the density exceeds n 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 grv and B sh , respectively. The B-n relation is characterized by n grv and B sh reasonably well for the models with θ = 2θ cr . The mean field strength for the models is well approximated by where β d is a typical plasma β 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.

STATISTICAL PROPERTIES OF DENSE CLUMPS
In this section, we investigate the physical properties of dense clumps. The clumps are identified by connecting adjacent cells whose densities are larger than a threshold number density of n th . We extract dense clumps consisting of more than 100 cells, whose internal structure is sufficiently well resolved at least with 6 cells along the diameter. Dib et al. (2007) estimated errors for the terms in the virial theorem by using a uniform sphere, and found that the errors are not larger than ∼ 15 % if the diameter of the sphere is resolved by four cells. Kobayashi et al. (2022) also found that the internal velocity dispersions of dense clumps are numerically converged if the clump size is resolved at least by 5 cells. This is because the largest scale of a clump, which corresponds to the clump size, has the largest power of the internal turbulence spectrum.
We choose a minimum threshold density of 10 3 cm −3 because we could not identify dense clumps as an isolated structure at lower densities. We here consider the dense clumps identified with four different threshold densities (n th = 10 3 cm −3 , 10 3.5 cm −3 , 10 4 cm −3 , and 10 4.5 cm −3 ). The spatial distributions of the dense clumps identified at t = t first + 0.5 Myr are displayed in Figure 6. The filamentary structures seen in the left panels of Figures 4 correspond to the dense clumps with n th ≥ 10 3.5 cm −3 .
We define the mass M cl , the position of the center of Spatial distributions of the dense clumps identified with four different threshold densities, n th = (blue)10 3 cm −3 , (orange)10 3.5 cm −3 , (green)10 4 cm −3 , and (red)10 4.5 cm −3 for models (a)Θcr0.25V5, (b)Θcr2V5, (c)Θcr0.25V10, and (d)Θcr2V10 at t = t first + 0.5 Myr. mass x cl , and density-weighted mean velocity v cl of an identified clump as follows: where V denotes the volume of each identified clump. In this paper, when considering the virial theorem, the surface terms are neglected although they can be important to judge the stability of dense clumps (Dib et al. 2007). In a future paper, we will investigate the stability of dense clumps by taking into account the surface terms as well as the time evolution of filamentary structures.
The virial theorem without the surface terms is 1 2 where I is the moment of inertia, E th and E kin are the thermal and kinetic energies and E mag is the magnetic energy, where u = (v − v cl ). In Equation (19), the gravitational energy is denoted by where r = (x − x cl ). E grv is not exactly the same as the gravitational energy of the clump, and includes the tidal effect from the external gas. In this paper, E grv is called the gravitational energy. It can be positive when the tidal force exerted by the external gas destroys the clump.
In preparation for interpreting the results, we express the energies in terms of the mean values, where c s is the mean sound speed, δu cl,1D is the onedimensional internal velocity dispersion, B cl is the mean field strength, V cl is the volume of dense clumps, ρ cl = M cl /V cl is the mean density. In derivation of Equation (26), we assume a uniform spherical cloud with a radius of R cl . In this paper, as a size of clumps, we use All the statistical quantities shown in this section are estimated using the dense clumps identified in t first − 0.5 Myr ≤ t ≤ t first + 0.5 Myr to increase the number of samples.

The Virial Parameter
From the virial theorem without the surface terms (Equation (19)), a stability criterion can be derived bÿ I < 0, which is rewritten as The virial parameter α 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 α tot on n th , M cl , V 0 , and θ, we investigate the velocity structure of the dense clumps that are related to E kin in Sections 4.1.1 and 4.1.2.

Bulk Speeds
In Figure 2, we found that the velocity dispersions of the post-shock layers and their anisotropy in the postshock layer depend not on V 0 significantly but on θ/θ cr . How are these velocity structures imprinted the motion and internal turbulence of the dense clumps? Figures 7 show the bulk speeds of the clumps as a function of the clump size at three different threshold densities (n th = 10 3.5 , 10 4 , and 10 4.5 cm −3 ). Figures  7a and 7b show that for all the models most clumps have supersonic bulk speeds because the sound speed is as low as 0.2 − 0.3 km s −1 (Koyama & Inutsuka 2002). The bulk speeds decrease with the clump size on average although there are large dispersions.
In order to investigate how the bulk velocities of the dense clumps depend on θ in more detail, we measure the one-dimensional translational velocity dispersions of the bulk motion of the dense clumps with R cl ≤ 0.2 pc that are defined as where the average is taken over all the dense clumps identified with each threshold density. Figure 8a shows ∆v cl as a function of the mean clump density n cl . All the model shows that ∆v cl decreases at a similar rate as n cl increases. The difference in θ does not have significant effect on ∆v cl . By contrast, ∆v cl increases roughly in proportion to V 0 . This comes from that the global velocity dispersions are in proportion to V 0 , but are roughly independent of θ (Figures 2c and  2d). At t = t first , the global velocity dispersions are ∼ 1.7 km s −1 for both the models with V 0 = 5 km s −1 , are ∼ 3.4 km s −1 for model Θ cr 0.25V10 and ∼ 2.4 km s −1 for model Θ cr 2V10.
We further examine how the anisotropy of the velocity dispersions of the post-shock layers shown in Figure 2b carries over the anisotropy of v cl . To characterize it, f aniso ≡ ∆v x,cl /∆v ⊥,cl is defined as a representative value of the degree of anisotropy, where ∆v 2 ⊥,cl = (∆v 2 y,cl + ∆v 2 z,cl )/2. (a) Velocity dispersions of the bulk motion of the dense clumps with R cl ≤ 0.2 pc for models (solid red)Θcr0.25V5, (solid blue)Θcr2V5, (dashed red)Θcr0.25V10, and (dashed blue)Θcr2V10 as a function of the mean clump density at four different threshold densities, n th = 10 3 , 10 3.5 , 10 4 , and 10 4.5 cm −3 . In each panel, the horizontal solid and dashed lines correspond to the collision speed and sound speed in the dense regions. Figure 8b shows that unlike ∆v cl , f aniso depends more strongly on θ than on V 0 . For both the models with θ = 2θ cr , f aniso is almost constant with n cl and roughly equal to unity. This is consistent with Figure 2b. By contrast, the anisotropy in the models with θ = 0.25θ cr are more significant toward lower n cl . At n th = 10 3 cm −3 , f iso is as high as ∼ 2.8, which is consistent with δv x / (δv 2 y + δv 2 z )/2 ∼ 2.8 for V 0 = 5 km s −1 and ∼ 3.8 for V 0 = 10 km s −1 at t = t first (Figure 2b).
For θ = 0.25θ cr , the larger the collision speed, the higher f aniso remains until the gas density becomes higher. This indicates that the ballistic motion of the dense clumps is maintained without much deceleration up to higher densities for larger V 0 (Figure 8).

Internal Velocity Dispersions
We here investigate the size-dependence of the internal velocity dispersion of each dense clump. The onedimensional internal velocity dispersions of the clumps δu cl,1D are plotted in Figures 9. They increase with the clump size in a power-law manner although the range of the clump sizes is small. We fit them using a power-law function of The best fit parameters are listed in Table 2. In the fitting, the clumps with n th = 10 4.5 cm −3 are removed because a increase of δu cl,1D due to self-gravity is found for R cl > 0.1 pc in Figure 9. The power-law indices shown in Table 2 are consistent with 0.4 − 0.5 for all the models, which is close to those of shock-dominated turbulence (a = 0.5, Elsässer & Schamel 1976) than incompressible Kolmogrov tur-bulence (a = 1/3). They are consistent with observed values (e.g., Larson 1981; Solomon et al. 1987).
In order to investigate the parameter dependence of δu cl,1D in more detail, we plot δu cl,1D mass-weighted averaged over the clumps with R cl ≤ 0.2 pc as a function of n cl in Figure 10. Like the bulk velocity dispersions ∆v cl , δu cl,1D depends more on V 0 than θ. For V 0 = 5 km s −1 , the internal velocity dispersions are as low as ∼ 0.2 − 0.3 km s −1 , regardless of n cl for both θ = 0.25θ cr and 2θ cr . The internal velocity dispersions at n th = 10 3 cm −3 increases roughly in proportion to V 0 . Their scatters in δu cl,1D also increase with V 0 . Unlike for V 0 = 5 km s −1 , δu cl,1D decreases with n cl for V 0 = 10 km s −1 . As a result, the difference in δu cl,1D between the models with different V 0 decreases as n cl . This feature is consistent with the results of Audit & Hennebelle (2010). Mass-weighted average of the internal velocity dispersions as a function of the mean clump density at n th = 10 3 , 10 3.5 , 10 4 , and 10 4.5 cm −3 for models (solid red)Θcr0.25V5, (solid blue)Θcr2V5, (dashed red)Θcr0.25V10, and (dashed blue)Θcr2V10. The average is performed over the dense clumps with R cl < 0.2 pc. The horizontal and vertical error bars indicate the dispersions of the mean density and internal velocity dispersion at each threshold density, respectively.

Virial Parameters
The virial parameter (Equation (28)) is divided into three contributions, thermal, kinetic, and magnetic virial parameters as follows:  Table 2. List of the best-fit formula of the internal velocity dispersions. From the second to fifth columns, the numbers in the parentheses represent δv0 in km s −1 . Figure 11. (first column) Thermal, (second column) kinetic, and (third column) magnetic virial parameters as a function of the clump mass at (red) n th = 10 3.5 , (green) 10 4 , and (blue) 10 4.5 cm −3 for models (first row) Θcr0.25V5, (second row) Θ2V5, (third row) Θcr0.25V10, and (forth row) Θcr2V10. The forth column shows the total virial parameter αtot = α th + α kin + αmag.
In each panel, the three dashed lines correspond to the predictions from the analytic formula (Equations (32), (34), (36), and (37)). Figure 11 illustrates the three virial parameters (α th , α kin , and α mag ) as a function of the clump mass at three density thresholds. One can see that the three virial parameters have different M cl and n th -dependence. Focusing on the M cl -dependence, as M cl increases, α th and α mag both decrease at a rate of ∝ M −0.6 cl while α kin decreases more slowly following ∝ M −0.3 cl . All the virial parameters decrease with n th , keeping their trend unchanged while the decreasing rate depends on θ and V 0 .
In order to explain the simulation results, we here derive analytic formula of α th , α kin , and α mag . First we consider the thermal virial parameter. Using Equations (23) and (26), Equation (31) is rewritten as where c 0.3 = c s /0.3 km s −1 , ρ 4 = ρ cl /(10 4 µm H cm −3 ), and M 1 = M cl /1 M . Equation (32) is essentially the same as that derived by Bertoldi & McKee (1992) and Lada et al. (2008). The thermal virial parameters are consistent with the predictions from Equations (32) for all the models (the first column of Figure 11). The kinetic virial parameter is estimated as where the internal velocity dispersion δu cl,1D = δu 0 (R cl /R 0 ) 0.5 is used (Section 4.1.2). Equation ( where we use δu 1 = δu 0 /1 km s −1 and R 0 = 1 pc for all the models because δu cl,1D roughly converges to ∼ 0.6 km s −1 (R cl /1 pc) 0.5 toward higher clump densities, regardless of the parameters V 0 and θ (Section 4.1.2). The mass-dependence of α kin is weaker than that of α th , indicating that α kin dominates over α th for more massive clumps. Comparison between α kin and the predictions from Equation (34) shows that Equation (34) reproduces the simulation results reasonably well in higher densities although deviations are found in lower clump densities for V 0 = 10 km s −1 as expected.
The magnetic virial parameter is expressed as The field strength B cl is estimated from the B-n relation shown in Figure 5. For the models with θ = 2θ cr , the density-dependence of the field strength is characterized by n grv . Substituting Equation (15) This equation shows that the magnetic virial parameter decreases in proportion to M −2/3 1 . This dependence is the same as that of α th as shown in Figure 11. The predictions from Equation (36) are consistent with the simulation results (the third column of Figure 11). Note that α mag is closely related with the mass-to-flux ratio that is estimated observationally by the Zeemann measurement (Crutcher 1999;Heiles & Troland 2005;Crutcher et al. 2010). We will discuss an asymptotic property of the mass-to-flux ratio of dense clumps in Section 5.1.
Combining Equation (32), (34), and (36), one obtains the analytic expression of the total α parameter, where the first term in the right-hand side corresponds to α th + α mag and the remaining term corresponds to α kin . It is worth noting the importance of the magnetic support in α tot,ana (Equation (37)). For clumps with ρ cl > ρ grv , the field strength is roughly determined by the condition β ∼ O(1). The third term in the parentheses of the first term in the right-hand side of Equation (37) is negligible compared with the other two terms. In such a dense clump, the magnetic field is only of secondary importance in α tot,ana because α th,ana is three times larger than α mag,ana when β d = 1 (α th /α mag = 2E th /E mag ∼ 3β). This finding suggests that the masses of clumps and cores are determined almost independently of the magnetic field strength once ρ cl exceeds ρ grv .
The diffuse clumps with n cl < n grv 3 = 10 3 cm −3 n 0 10 cm −3 V 0 5 km s −1 2 c s 0.3 km s −1 −2 (38) are expected to be mainly supported by the shockamplified magnetic fields if the contribution of the kinetic energy is excluded. This is because the third term dominates over the other two terms in the parentheses of the first term in the right-hand side of Equation (37). In our setup, the condition is met when n cl < 10 3 cm −3 for V 0 = 5 km s −1 and n cl < 4 × 10 3 cm −3 for V 0 = 10 km s −1 . This is why magnetically-supported dense clumps with n th ≥ 10 4 cm −3 are not seen in our results.

Mass-to-Flux Ratio
The mass-to-flux ratio µ is often used to evaluate the importance of magnetic fields in cloud stability. If the mass-to-flux ratio is larger than the critical value of µ cr = 1/2π √ G, magnetic fields are not strong enough to support clumps against self-gravity A factor of 1/2π slightly changes depending on the clump shape (Strittmatter 1966;Mouschovias & Spitzer 1976;Nakano & Nakamura 1978;Tomisaka et al. 1988;Tomisaka 2014). Observationally µ is often estimated by taking the ratio between the column density and the field strength measured by the Zeeman effect (Crutcher 1999;Heiles & Troland 2005;Crutcher et al. 2010).
In terms of α mag , the mass-to-flux ratio is expressed as From Equation (36), µ/µ cr is expected to be proportional to M 1/3 cl , regardless of θ. This mass dependence was confirmed in the third column of Figure 11. In simulations of cloud-cloud collision, Sakre et al. (2021) found that as the clump mass increases, the magnetic energy increases more slowly than the gravitational energy, which is consistent with the positive M cl dependence of the mass-to-flux ratio. Considering compressions along the magnetic field, Banerjee et al. (2009) and Inoue & Inutsuka (2012) found a similar mass-dependence of the mass-to-flux ratios which follow ∝ M 0.4 cl . Iffrig & Hennebelle (2017) also found a similar relation in semiglobal galactic-disk simulations, and explained the massdependence by a similar argument as in this paper. We plot the mass-to-flux ratios averaged over 1 M ≤ M cl ≤ 5 M as a function of the mean clump density in Figure 12. The θ = 2θ cr models show that the mass-to-flux ratios are consistent with the predictions from the analytic formula (Equation (36)). Equation (36) for n cl > n grv and µ sh µ cr = 1.1 ρ 0 for n cl < n grv . The mass-to-flux ratios increase in proportion to n 2/3 cl until n cl reaches n grv . When n cl exceeds n grv , the mass-to-flux ratios approach the asymptotic line given by the condition β ∼ 1. The difference in V 0 appears in how fast µ approaches µ β .
The difference in µ/µ cr between the models disappear when n cl exceeds n grv because the field strength is determined by the condition β ∼ O(1) (Figure 5). Comparison between the simulation results and the predictions from Equation (36) with M cl = 3 M shows that they are consistent for both the models.
For n cl > n grv , the mass-to-flux ratios are insensitive to n cl (Figure 12) since µ β ∝ n 1/6 cl . This means that the mass-to-flux ratio does not change significantly once n cl becomes larger than n grv . We thus define the characteristic mass-to-flux ratio µ ch for which n cl is equal to n grv , µ ch µ cr = 0.5β Mass-to-flux ratios of the dense clumps with n th = 10 4.5 cm −3 as a function of the clump mass for models (red)Θcr0.25V5, (green)Θcr2V5, (blue)Θcr0.25V10, and (magenta)Θcr2V10. The solid and dashed lines correspond to µ ch /µcr with V0 = 5 km s −1 and V0 = 10 km s −1 , respectively. Figure 13 shows the mass-to-flux ratios of the dense clumps with n th = 10 4.5 cm −3 , which are almost independent of the parameters V 0 and θ. This is consistent with what we found in Figure 12. We compare the massto-flux ratios of the dense clumps with n th = 10 4.5 cm −3 with the predictions from Equation (42) in Figure 12. One can see that they are consistent. The reason why the measured mass-to-flux ratios are slightly larger than the predictions is that the clump density is larger than n grv .
Interestingly, µ ch /µ cr is roughly equal to unity, and it is insensitive to the collision parameters ( n 0 , V 0 ). The sound speed is expected to decrease to 0.2 km s −1 when the gas density increases further, µ ch /µ cr increases by a factor of two (Equation (42)). We thus expect that dense clumps with ∼ 1 M become super-critical if their mean densities exceed n grv in a wide range of the collision parameters.

Internal Alfén Mach number
We here show the ratio between the kinetic energy and magnetic energy, or Alfvén mach number, Using Equations (15) and (30) for n < n grv . If the ram pressure of the colliding flow ρ 0 V 2 0 is extremely large, the Alfvén Mach number can be small in the low density regions with n < n grv . Once the number density exceeds n grv , Equation (44) shows that the Alfvén Mach number converges to order of unity. This result is robust because it is extremely insensitive to both the clump density and mass. This feature is consistent with observations by Myers & Goodman (1988) and Crutcher (1999).

Structures of Filamentary Clouds
In our models, prominent filamentary structures form only for θ > θ cr (Figure 4) because super-Alfvénic turbulence prevents such a coherent structure from forming for θ < θ cr . Our results are consistent with the observational results (e.g., André et al. 2010) because compressions with θ < θ cr is extremely rare if compressions are isotropic (Paper I). In reality, the magnetic fields are not uniform but have different orientations, depending on the location. Compressions with θ < θ cr are unlikely to occur over a wide area.
In this section, the internal structures of the filaments are briefly shown for model Θ cr 2V5 although they will be investigated in detail in forthcoming papers. We extract one of the filaments from the simulation results and its column density is shown in the top-left panel of Figure 14. Overall, the projected magnetic fields are perpendicular to the filament, which is consistent with observational results (e.g., Alves et al. 2008;Sugitani et al. 2010;Chapman et al. 2011).
The three dimensional structure of filaments is not a cylinder but more like a "ribbon". The top-right and bottom panels of Figure 14 correspond to the density slices at z = z a , z b , and z c , whose positions are shown in the top-left panel of Figure 14. The cross sections of the filament are flattened along the local magnetic field direction. Ribbon-like flattened structures are a natural configuration of strongly-magnetized filaments because the anisotropic nature of the Lorentz force, whose direction is perpendicular to the local magnetic field (Tomisaka 2014;Auddy et al. 2016).
Interestingly, the elongated direction of the flattened density structures varies significantly along the filament axis because the directions of the magnetic fields significantly vary as shown in the density-slice plots of Figure  14. For instance, at z = z b , the minor axis of the flattened structure found around (y, x) ∼ (7 pc, 0 pc) is parallel to the line-of-sight direction (the x-axis). That is why the widths of the filament around z = z b appear to be larger than those of the other regions (the top panel of Figure 14).

A Universal B-n Relation
The B-n relation for high densities reflects how the magnetic field is amplified during gas contraction. The gravitational collapse of a spherical cloud with extremely weak magnetic fields leads to the scaling law B ∝ n 2/3 (Mestel 1965). By contrast, a cloud flattened by the existence of the magnetic field undergoes gravitational contraction with B ∝ n 1/2 , indicating that the plasma beta remains constant (Mouschovias 1976;Tomisaka et al. 1988).
Our analytic formulae derived in Section 4 rely on the speculation that the field strength is determined by a condition of β ∼ O(1) in the high density limit, regardless of θ.
We first present a simple argument of why the central plasma β of gravitationally contracting flattened clouds should be O(1). We then compare the B-n relations of our work with those of the previous studies in Section 5.4.2, and those inferred from observational studies in Section 5.4.3.

Flattened Clouds
In Section 5.3, we found that the dense clumps have ribbon-like structures flattened along the local magnetic The top-left panel shows the column density map in 4 pc ≤ y, z ≤ 9 pc for model Θcr2V5. The white arrows indicate the directions of the projected magnetic fields density-weighted-averaged along the x-axis. The remaining three panels correspond to the density slices at z = za = 7.6 pc, z b = 6.6 pc, and zc = 4.6 pc. The red lines indicate the contours of n = 10 3.5 cm −3 . The black lines show the stream lines of the vector field (By, Bx). field for model Θ cr 2V5. Also for the other filaments of models Θ cr 2V5 and Θ cr 2V10, similar ribbon-like structures are found. Thus, at least for θ = 2θ cr , the dense clumps form through gas accumulation along the magnetic fields.
We here investigate how magnetic fields are amplified by self-gravity using a simple ribbon model. As shown in Figure 15, a flattened ribbon-like filament with the central density ρ c , the line mass M L , and the width d 0 is considered. The uniform magnetic field, whose strength is B 0 , is perpendicular to the ribbon surface. The thickness of the ribbon is determined by the thermal The relation between ρ c and column density Σ = M L /d 0 is The flux-freezing condition gives We here consider the situation that the ribbon grows Figure 15. Schematic picture of the ribbon with the width d0 and line mass ML. The ribbon extends in the direction parallel to the sheet. The central density is assumed to be ρc, and the density is assumed to be independent of the width. The uniform magnetic field (field strength B0) is parallel to the minor axis of the ribbon.
by gas accretion along the magnetic field. Initially, the disk is not massive enough to contract radially due to the magnetic pressure. In order for the disk to contract radially, the mass-to-flux ratio be larger than the critical value, (Nakano & Nakamura 1978). From this, we define a critical column density of Once the column density of the ribbon reaches this critical value, the ribbon begins to contract radially. The line mass at this epoch is denoted by M L,cr , This is essentially the same as the critical line mass at the large magnetic flux limit derived by Tomisaka (2014). Using Equation (47) and (48) This result leads to the important conclusion that the plasma β of dense clumps is order of unity when their gravitational contraction starts. A similar argument can be done in an oblate cloud flattened along the magnetic field. The corresponding expression for the central plasma β is given by where M is the cloud mass and where R 0 is the radius of the oblate. An important conclusion given from Equations (52) and (53) is that the plasma β should be larger than unity in contracting flattened clouds.
We should note that the above analytic argument cannot be applied to the models with θ = 0.25θ cr . Despite the fact that the dense clumps are not necessarily flattened along the magnetic field, the B-n relations of the models with θ = 0.25θ cr are consistent with the line β ∼ O(1) (Figure 5). More general physics may be hidden. Figure 16 shows the B-n relations of the simulation studies with various setups. The references prefixed by (C) (Hennebelle et al. 2008;Banerjee et al. 2009) investigate the colliding flows of atomic gases as in our study. All the B-n relations of the colliding-flow simulations follow the lines of β = 1 for high densities.

B-n Relations in Previous Theoretical Studies
Next, our B-n relations are compared with those of the colliding-flow of MCs that is a promising mechanism to form massive cores which will evolve into massive stars (Inoue & Fukui 2013;Inoue et al. 2017). Their field strengths approach the lines of β = 1 as the density increases, which is consistent with our finding. For densities higher than ∼ 10 8 cm −3 , the density-dependence of the field strength becomes weaker than B ∝ n 1/2 . This may be explained if the mass increases during the gas contraction (see Equation (52)).
We should note that the magnetic field for lower densities is stronger than ours because their ram pressure is an order of magnitude larger than ours. The strong magnetic field prevents the gas from becoming gravitationally unstable before massive cores form, and the dense pre-shock gas provides a sufficient amount of mass to form massive cores.
We also compare our B-n relations with those of supernova-regulated ISM MHD simulations (Padoan et al. 2016;Hennebelle 2018). Padoan et al. (2016) carry out a simulation with periodic boundary conditions in a simulation box of (250 pc) 3 . Their results are consistent with our results although their simulation setup is quite different from ours. Hennebelle (2018) consider the stratification of the galactic disk, and their simulations resolve the gas dynamics from the kpc scale to 0.004 pc by using a zoom-in technique. Their results show that the field strength follows B ∝ n 1/2 for higher densities, but it is slightly stronger than that predicted from β = 1 even when the gas density is as high as 10 6 cm −3 if the gas temperature is as low as 10 K.  (2013) and Inoue et al. (2017). The two gray lines correspond to β10K = 1 and β20K = 1. The observational prediction in Crutcher et al. (2010), B ∼ max(10 µG(n/300 cm −3 ) 0.65 , 10 µG), is shown by the dashed line. Mocz et al. (2017) performed isothermal simulations of driven turbulence with different initial field strengths, and demonstrated that the outer regions (∼ 10 4 au) of the cores have β ∼ 1, regardless of the initial field strength. The regions of B ∝ n 2/3 , where gravitationally collapse isotropically, emerge on small scales (r < 10 4 au) only when the initial Alfvén Mach number is less than unity. In our simulations, the phase of B ∝ n 2/3 cannot be resolved because of the lack of resolution. Li et al. (2015) found that the mean field strength of clumps does not follow n 1/2 but it is proportional to n 0.7 , which is comparable to those found in observations of the Zeeman measurement (see Section 5.4.3). The mean field strength of their results is B ∼ 20 µG n/10 4 cm −3 0.7 , regardless of the initial field strength. It is higher than those obtained in our simulations. An possible reason is the difference of the simulation setups. In our simulations, dense clumps evolve from the stable stage to the unstable stage through coalescence of smaller clumps and the gas accumulation of the atomic gas. By contrast, Li et al. (2015) prepared a gravitationally unstable gas as the initial condition.
Rapid super-Alfvénic contraction of the dense clumps may lead to a strong density-dependence of the field strength (B ∝ n 0.7 ).
In summary, the plasma β in dense regions is order of unity in previous studies conducting colliding flow simulations (Hennebelle et al. 2008;Banerjee et al. 2009;Inoue & Fukui 2013;Inoue et al. 2017). However, in other setups, some simulations show consistent results (Padoan et al. 2016;Girichidis et al. 2018), but some studies are not consistent with ours (Hennebelle 2018;Li et al. 2015). Further investigations are needed to understand what makes this difference in the high density limit of the B-n relation

B-n Relations in Observational Studies
Analysing the results of Zeeman surveys of Hi, OH, and CN, Crutcher et al. (2010) derived an expected field strength as a function of density that is given by B obs ∼ 10 µG(n/300 cm −3 ) 0.65 for n > 300 cm −3 . It is about three times stronger than our results at n = 10 3 cm −3 , and the difference increases with the density, indicating that the plasma β of the observed clouds is much lower than unity. Crutcher (1999) found that the plasma β is as low as 0.04 +0.03 −0.02 , which appears to be contradict with our results.
A reason of this discrepancy may come from the fact that the Zeeman surveys have mainly been conducted in massive star forming regions. Suppose that a spherical cloud with a density of ρ and mass of M has a magnetic field with a strength of B, the mass-to-flux ratio is given by Solving equation (55) The equation indicates that in order for a cloud with β ∼ 0.04 to be gravitationally unstable, it should be massive. This is consistent with the fact that most data in high density range analyzed in Crutcher et al. (2010) are associated with massive star formation. Inoue & Fukui (2013) and Inoue et al. (2017) demonstrated that magnetically-supported dense cores formed by cloudcloud collisions are possible sites of massive star formation ( Figure 16). An averaged plasma β of such dense cores is expected to be much lower than unity as seen in the Zeeman observations although the plasma beta of the central densest region is order of unity ( Figure 16). We focuses on low and intermediate mass star formation. For model Θ cr 2V10, a typical plasma β and den-sity of the post-shock gas are ∼ 0.025 and 200 cm −3 , respectively, if we assume that the pre-and post-shock gas are uniform and isothermal with a sound speed of 0.3 km s −1 . Only forty-times density enhancement is enough for the plasma β to become unity. Thus, our results are not contradict with the observations measuring the field strength through the Zeeman effect.

Relative Orientations between Density Structures and Magnetic Fields
To characterize the relation between the density and magnetic fields, Soler et al. (2013) define that the angle ϕ between the magnetic field and gradient of the density using the following expression: where ∇n is calculated using the Gaussian derivative kernel with the standard deviation ∆x/ √ 2 to reduce numerical fluctuations.
The two-dimensional version of Equation (57), where n is replaced by observed surface densities and the direction of B is estimated from polarization observations, is often used to link the simulation and observation studies (e.g., Planck Collaboration et al. 2016a). Three-dimensional HROs obtained from the snapshots at t = t first + 0.5 Myr for models (a) Θcr0.25V5, (b) Θcr2V5, (c) Θcr0.25V10, and (d) Θcr2V10. In each panel, the lines correspond to the HROs at the five different density bins. Figure 17 shows the histogram of cos ϕ called Histogram of Relative Orientations (HROs), which was proposed by Soler et al. (2013). For the lowest density bin (10 2.5 cm −3 ≤ n < 10 3 cm −3 ), the HROs have sharp peaks around cos ϕ ∼ 0 for all the models, indicating that the density structure tends to be elongated along the magnetic field. The directions of the density elongation are different between the models with θ = 0.25θ cr and 2θ cr . For θ = 0.25θ cr , the strong shear motion, which is generated by the anisotropic super-Alfvénic turbulence (Figure 2c), stretches both the density structure and magnetic fields in the collision direction. By contrast, for θ = 2θ cr , the magnetic fields amplified by the shock compression are perpendicular to the collision direction, and the gas is elongated along the magnetic field.
As the density increases, the peaks around cos ϕ ∼ 0 become less prominent while the number of counts around cos ϕ ∼ ±1 increases for all the models. The highest density bin (10 4.5 cm −3 ≤ n ≤ 10 5 cm −3 ) corresponds to the densest clumps identified in this paper. Figure 17 shows that the HROs at the highest density bin depend on θ/θ cr . For θ = 0.25θ cr , the HROs become almost flat, indicating that any angles between ∇n and B are equally realized. For θ = 2θ cr , the magnetic fields tend to be perpendicular to the direction of the elongation of the density structures for n > 10 4 cm −3 . This clearly reflects that the filamentary structures are elongated perpendicular to the magnetic field (the right panels of Figure 4).
These properties of the HROs are consistent with other simulation studies where the flip from the regime of ∇n ⊥ B to that of ∇n B occurs only in high magnetization cases with β 1 (Soler et al. 2013;Barreto-Mota et al. 2021) although what determines the critical densities at which the flip occurs is not fully understood theoretically (Pattle et al. 2022). For the models with θ = 2θ cr , the critical density (∼ 10 4 cm −3 ) roughly agrees with the mean density at which the dense clumps with > 1 M change from sub-critical to supercritical ( Figure 12). This is consistent with the results of Seifried et al. (2020).
Planck Collaboration et al. (2016a,b) found that the relative angle ϕ systematically changes with increasing surface density in a consistent way as the simulation results. With the clump mass ∼ 10 M around which α tot becomes less than unity (the bottom panels of Figure 11), the transition density n ∼ 10 4 cm −3 corresponds to the surface density N ∼ n × (4π/3 × (10 M )/(1.4nm H cm −3 )) 1/3 ∼ 10 22 cm −2 , which appears to be consistent with the value N H ∼ 10 21.7 cm −2 obtained from the Planck observations (Planck Collaboration et al. 2016b). More rigorous comparisons will be made in forthcoming papers by applying synthetic observations to our results.

Density Dependence of Velocity Dispersions
Figures 9 and 10 show that the δu cl,1D -R relations depends on the mean clump density although the dynamic range of the clump sizes is less than one decade. The density dependence is more significant for larger V 0 . Observationally, the density-dependence of the δu cl,1D -R relation can be investigated by using multiple lines of different molecules and isotopes (Goodman et al. 1998).
In the L1551 cloud in Taurus, which is an isolated low-mass star-forming region, Yoshida et al. (2010) do not find any clear deferences in the three line width-size relations (σ NT ∼ 0.2 km s −1 (L/1 pc) 0.5 ) obtained from the 12 CO (J = 1 − 0), 13 CO (J = 1 − 0), and C 18 O (J = 1 − 0) lines. The obtained velocity dispersions are transonic ∼ 0.2 km s −1 even at L ∼ 1 pc, and appear to be smaller than a typical line width-size relation (e.g., Larson 1981; Solomon et al. 1987). This observational result may be qualitatively consistent with the models with V 0 = 5 km s −1 .
Recently, in the filamentary infrared dark cloud G034.43+00.24, Liu et al. (2022) found that the velocity dispersion-size relation in smaller scales (∼ 0.09 pc) gives systematically lower velocity dispersions than that in larger scales (∼ 0.6 pc). This may be qualitatively consistent with the models with V 0 = 10 km s −1 , which show the negative dependence of δu cl,1D on the mean clump density (Figure 10).
On the other hand, Liu et al. (2022) found that the turbulence spectrum toward a filamentary infrared dark cloud G034.43+00.24 observed in the dense-gas tracer H 13 CO + shows a systematic deviation from the Larson's law and the slope is steeper. This may be qualitatively consistent with the models with V 0 = 10 km s −1 . We attribute this deviation from the Larson's law to the gravitational collapse in the dense region.

Core Formation
We found that for n > n grv , µ/µ cr becomes order of unity. This indicates that the magnetic energy is comparable to the gravitational energy. Since plasma beta is order of unity for n > n grv , the thermal energy is also comparable to the magnetic energy. In Section 5.2, we found that the Alfvén mach number is also order of unity for n > n grv . Combining these results suggests that all the energies (thermal, kinetic, magnetic, and gravitational energies) become equipartition when n > n grv , regardless of θ. This property has been found in previous studies. Lee & Hennebelle (2019) for instance found that E mag , E th , and E grv reach equipartition in the highdensity ends of their runs. Chen & Ostriker (2014) found that core masses do not depend on the upstream magnetic field direction in their colliding-flow simulations. They estimated the core mass by using a theoretical argument considering gas accumulation along the magnetic field. In terms of our notation, their estimated core mass corresponds to the Jean mass with n = n grv = n 0 (V 0 /c s ) 2 ∼ 10 5 cm −3 .
From Equation (42), the cores with n ∼ n grv are expected to be super-critical when the core mass is around 1 M . Our results suggest that the contribution of the magnetic field to the virial theorem is negligible when n > n grv /3 = 3 × 10 4 cm −3 (Equation (38)). This is consistent with their results, which show that a typical core mass is independent of magnetic fields. Chen & Ostriker (2015) investigated the dependence of the core properties on the collision speed and field strength. The plasma β of the cores takes values roughly between 1 and 10, consistent with our results. They found that the core collapse time is proportional to V −1/2 0 (also see Iwasaki & Tsuribe 2008;Gong & Ostriker 2011). This appears to be contradict with our results, which show that the clump formation time is insensitive to V 0 (Table 1 and Figure 3). This discrepancy may come from the difference in the disturbances of the pre-shock gas. Relatively weak turbulence (δv = 0.14 km s −1 in a box of (1 pc) 3 ) was considered in Chen & Ostriker (2015) while the two-phase structure of the atomic gas is considered in this paper. In our setting, higher collision speeds do not make the clump formation time shorter because stonger post-shock turbulence is driven.

Long-term Evolutions
In this paper, we focus on the early evolution of the dense clumps since our simulations are terminated at 0.5 Myr after the formation of the first unstable clump because a star formation treatment is not included and the resolution is not high enough to resolve very dense cores. In this section, we here discuss whether our results are applicable to the long-term evolution.
In the long-term evolutions, the clumps will grow through gas accretion and denser cores will form. The analytic formulae of the virial parameters shown in Equations (32), (34), and (36) are all applicable in the long-term evolutions. α th parameter (Equation (32)) is unlikely to change significantly unless the gas is heated up locally by the star formation. The magnetic energies of the dense clumps are strongly related with the B-n relation. If the magnetic field amplification is caused by the gravitational contraction also in the long-term evolutions, Equation (36) is also unlikely to change significantly. It is possible that E kin increases as the gravitational potential becomes deeper (Arzoumanian et al. 2013), or δu 1 in Equation (34) increases. The longterm evolutions with higher resolution are investigated in forthcoming papers.

SUMMARY
We investigated the formation and evolution of MCs by colliding flows of the atomic gas with a mean density of n 0 = 10 cm −3 and field strength of B 0 = 3 µG. We additionally include self-gravity in the simulations done in Paper I and investigate the global properties of post-shock layers and the statistical properties of dense clumps.
As parameters, we consider the collision speed and the angle between the magnetic field and upstream flow. As examples of the turbulent-dominated and magnetic field-dominated layers, we adopt θ = 0.25θ cr and 2θ cr , respectively. The critical θ cr is found in Paper I and defined in Equation (7). We consider two collision speeds of 5 km s −1 and 10 km s −1 .
Our findings are as follows: 1. The θ-dependence of the physical properties of the post-shock layers affect on the clump formation.
• For θ = 0.25θ cr , anisotropic super-Alfvénic turbulence driven for θ = 0.25θ cr suppresses the formation of self-gravitating clumps. The post-shock layers are significanly disturbed by turbulence, and no prominent filamentary structures are found.
• For θ = 2θ cr , the gas accumulates along the shock amplified magnetic field, and prominent filamentary structures form.
Self-gravitating clumps are formed more efficiently for θ = 2θ cr than for θ = 0.25θ cr .
2. The difference of V 0 does not influence the epoch of the formation of the first unstable clump t first significantly although the mass accumulation rate increases in proportion to V 0 . This is because all the energies are in proportional to V 2 0 , and their ratios are independent of V 0 .
3. The parameter dependence of the field strengthdensity relation appears only for low densities. The mean field strength approaches a univeral relation that determined by β ∼ O(1) as the density increases, regardless of values of θ and V 0 (Figure 5). The critical density which divides the two asymptotic behaviors is expressed as n grv = n 0 (V 0 /c s ) 2 (Equation (14)). We confirm that the B-n relations derived in previous studies with different setups approach β ∼ O(1) in the highdensity limit, consistent with our results ( Figure  16).
4. The physical properties of the bulk velocities of dense clumps are inherited from the turbulence structure of the post-shock layers. The models with θ = 0.25θ cr show that the bulk velocities are biased in the collision direction for less dense clumps. As the clump density increases, the anisotropy of the bulk velocities becomes insignificant. Anisotropy in the bulk velocities is not found in the models with θ = 2θ cr .
5. The internal velocity dispersions of dense clumps roughly obey ∝ R 0.5 cl , where R cl is a typical clump radius, which is defined as R cl = (3V cl /(4π)) 1/3 , where V cl is the volume of dense clumps. The internal velocity dispersions δu cl,1D depend not on θ but on V 0 . For diffuse clumps with ∼ 10 3 cm −3 , δu cl,1D is roughly in proportion to V 0 . The larger the collision speed, the faster δu cl,1D decreases with the clump density. As a result, the parameter-dependence of δu cl,1D becomes insignifiant for denser clumps.
6. The virial parameter α tot is divided into three contributions from the thermal pressure α th , Raynolds stress α kin , and Maxwell stress α mag in Equation (28). The thermal and magnetic virial parameters are proportional to M −2/3 cl and kinetic virial parameter is proportional to M −1/3 cl , where M cl is the clump mass. We develop an analytic formulae for α th in Equation (32), for α kin in Equation (34), and for α mag in Equation (36). They explain our results reasonably well. 7. We found that the mass-to-flux ratios of dense clumps with ∼ 1 M are expected to be order of unity if the clump density is higher than n grv in a wide range of the collision parameters which reproduces the observed relations.