Bimodal Behavior and Convergence Requirement in Macroscopic Properties of the Multiphase Interstellar Medium Formed by Atomic Converging Flows

We utilize the hydrodynamics part from the magnetohydrodynamics code by Inoue & Inutsuka (2008). This solves equations of hydrodynamics based on the second-order Godunov scheme (van Leer 1979). Heating and cooling are explicitly time integrated and have second-order accuracy. We solve the following basic equations:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{{\rm{\nabla }}}_{\mu }(\rho {v}_{\mu })=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {v}_{\mu })}{\partial t}+{{\rm{\nabla }}}_{\nu }(P{\delta }_{\mu \nu }+\rho {v}_{\nu }{v}_{\mu })=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+{{\rm{\nabla }}}_{\mu }\left((e+P){v}_{\mu }\right)={{\rm{\nabla }}}_{\mu }(\kappa (T){{\rm{\nabla }}}_{\mu }T)-\rho { \mathcal L }.\end{eqnarray} \tag{ 3 }$$

Here, ρ represents the mass density, v represents the velocity, P is the thermal pressure, T is the temperature, and ${{\rm{\nabla }}}_{\mu }=\partial /\partial \mu$, where μ spans x, y, and z. The total energy density, e, is given as e = P/(γ − 1) + ρv²/2, where γ means the ratio of the specific heat (=5/3). We implement the thermal conductivity, κ, as κ(T) = 2.5 × 10³ T^0.5 erg cm⁻¹ s⁻¹ K⁻¹ by considering collisions between hydrogen atoms (Parker 1953). The net cooling rate per mass ${ \mathcal L }$ is the integration of heating and cooling processes (see Figure 1). We utilize the following net cooling function,

$$\begin{eqnarray}&&\begin{array}{l}\rho { \mathcal L }=-\left(\displaystyle \frac{\rho }{{m}_{\mathrm{gas}}}\right){\rm{\Gamma }}+{\left(\displaystyle \frac{\rho }{{m}_{\mathrm{gas}}}\right)}^{2}{\rm{\Lambda }},\\ {\rm{\Gamma }}=2\times {10}^{-26}\,\mathrm{erg}\,{{\rm{s}}}^{-1},\\ \displaystyle \frac{{\rm{\Lambda }}(T)}{{\rm{\Gamma }}}\\ \ =\ \left\{\begin{array}{l}{10}^{7}\exp \left(\displaystyle \frac{-118400}{T+1000}\right)\\ \ +\ 1.4\times {10}^{-2}\sqrt{T}\exp \left(\displaystyle \frac{-92}{T}\right)\,{\mathrm{cm}}^{3}\,(\mathrm{for}\,T\leqslant 14,577\,{\rm{K}}),\\ 5\times {10}^{3}+1.4\times {10}^{-2}\sqrt{T}\exp \left(\displaystyle \frac{-92}{T}\right)\,{\mathrm{cm}}^{3}\\ (\mathrm{for}\,14,577\,{\rm{K}}\lt T\leqslant 19,449\,{\rm{K}}),\\ 3.75\times {10}^{4}\left(1-\tanh \left(\displaystyle \frac{T-2\times {10}^{5}}{2\times {10}^{5}}\right)\right)\\ \ \exp \left(\displaystyle \frac{-5\times {10}^{4}}{T}\right)\\ \ +\ {10}^{3}\exp \left(\displaystyle \frac{-5\times {10}^{4}}{T}\right)\,{\mathrm{cm}}^{3}\,(\mathrm{for}\,T\gt 19,449\,{\rm{K}}),\end{array}\right.\end{array}\end{eqnarray} \tag{ 4 }$$

where Γ is the heating rate, m_gas is the mean particle mass, and all T represent the temperature in units of Kelvin. The heating process, Γ, comes from photoelectric heating by polycyclic aromatic hydrocarbons. The cooling process, Λ(T), is proposed in Koyama & Inutsuka (2002), which is based on the detailed calculation of heating and cooling rates in optically thin ISM from Koyama & Inutsuka (2000). We use this formula in T ≤ 14,577 K, which mainly consists of two terms: the Lyα cooling (the first term) and the C_II cooling (the second term). We modify this function at a higher temperature regime as shown in the equations for T > 14,577 K by considering He, C, O, N, Ne, Si, Fe, and Mg (Cox & Tucker 1969; Dalgarno & McCray 1972; Figure 1).⁶ The typical cooling timescale of the injected WNM, ${\tau }_{\mathrm{cool}}$, is 1.3 Myr under this cooling function. This timescale becomes shorter once the injected WNM is compressed by the shock.⁷

Zoom In Zoom Out Reset image size

We define the thermally unstable neutral medium (UNM) as ${(\partial ({ \mathcal L }/T)/\partial T)}_{P}\lt 0$ (e.g., Balbus 1986, 1995), and the gas state warmer (colder) than the UNM is classified as the WNM (CNM). We will use this definition hereafter when we analyze the mass fraction of each phase. Although the boundary between WNM, UNM, and CNM depends on both the density and temperature, the three phases roughly correspond to the WNM as T ≥ 5000 K, the UNM as 100 ≤ T < 5000 K, and the CNM as T < 100 K based on Equation (4).

Figure 2 shows a schematic view of our simulations. The shock-compressed layer forms at the box center sandwiched by two shock fronts and becomes thicker while the flow continues as the two shocks propagate outwards. Our three-dimensional simulation box has size L_x,y,z = 20, 10, 10 pc, which is a typical size of giant molecular clouds in the Milky Way galaxy. Here, x is defined as the flow direction. The left and right parts of the WNM have velocities in the opposite direction with $\left|{v}_{x}\right|=20$ km s⁻¹, colliding at the box center. The relative velocity between these two flows is thus 40 $\mathrm{km}\,{{\rm{s}}}^{-1}$ in this 20 $\mathrm{km}\,{{\rm{s}}}^{-1}$ + 20 $\mathrm{km}\,{{\rm{s}}}^{-1}$ collision. This velocity is likely even faster in galaxy mergers (e.g., >100 $\mathrm{km}\,{{\rm{s}}}^{-1}$), but by employing 20 km s⁻¹, we are opting to focus on more common situations in galactic disks (e.g., the late phase of supernova remnant expansions, H ii region expansions, and normal shock due to galactic spirals⁸ ). For the initial velocity field, we may simply flip the sign of the velocity at x = L_x/2. However, as a more conservative approach to avoid any artifacts that may arise from such a step-function collision, we apply $\tanh$ smoothing as

$$\begin{eqnarray}&&{v}_{x}={V}_{\mathrm{in}}\tanh \left(\displaystyle \frac{x-{L}_{x}/2}{{n}_{\mathrm{smx}}{\rm{\Delta }}x}\right)\end{eqnarray} \tag{ 5 }$$

so that the initial collision occurs smoothly. Here, ${V}_{\mathrm{in}}=20$ km s⁻¹ and Δx is the mesh size. ${n}_{\mathrm{smx}}$ is chosen such that the physical scale of this smoothing is constant, ${n}_{\mathrm{smx}}{\rm{\Delta }}x=0.78$ pc between different resolution runs.

Zoom In Zoom Out Reset image size

The WNM is continuously injected through the two x boundaries at x = 0 and 20 pc until the calculation ends at 3 Myr, whereas the y and z boundaries have the periodic boundary condition. The WNM flow is in a thermally stable phase having the mean number density ${n}_{0}=0.57\,{\mathrm{cm}}^{-3}$ and pressure ${P}_{0}/{k}_{{\rm{B}}}=3500\,{\rm{K}}\,{\mathrm{cm}}^{-3}$. The corresponding mean temperature, sound speed, and dynamical pressure are 6141 K, ${C}_{{\rm{s}}}=8.16\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and 2.6 × 10⁴ K cm⁻³, respectively, where we use ρ₀ = n₀μ_Mm_p with μ_M = 1.27 as the mean molecular weight (cf. Inoue & Inutsuka 2008, 2012).

The mass density of the WNM flow has fluctuation ρ (x, y, z) = ρ₀ + δρ₀(x, y, z):

$$\begin{eqnarray}\begin{array}{rcl}\delta {\rho }_{0}(x,y,z) & = & \displaystyle \sum _{{k}_{x},{k}_{y},{k}_{z}}A(k)\\ & & \times \ \sin ({k}_{x}x+{k}_{y}y+{k}_{z}z+{\alpha }_{{k}_{x},{k}_{y},{k}_{z}}).\end{array}\end{eqnarray} \tag{ 6 }$$

k_x, k_y and k_z are the wavenumbers in the x, y and z directions: k_x = 2πl_x/(L_x/2), k_y = 2πl_y/L_y, and k_z = 2πl_z/L_z while the integers l_x, l_y, and l_z span from −32 to 32. A(k) is set such that the density power spectrum follows the Kolmogorov spectrum P_ρ(k) ∝ k^−11/3, where $k=\sqrt{{k}_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}}$ (Kolmogorov 1941; Armstrong et al. 1995). As indicated in the definition of k_x, we generate this fluctuation over a 10 pc × 10 pc × 10 pc volume so that the initial left half (x = 0–10 pc) and right half (x = 10–20 pc) of the WNM have the same density distribution (as the temperature distribution does to achieve the initial pressure equilibrium). The WNM flow injected from the boundaries also follow this density distribution as the flows move inward, and the WNM always accretes onto the shock-compressed layer with this density fluctuation.

To identify the shock-compressed layer, we search cells along the x direction at every given (y, z) inward from both x boundaries and label the first cells with P > 1.3P₀. We calculate the arithmetic mean of the distances from the box center to these positions as the representative shock-front position:

$$\begin{eqnarray}&&\langle {x}_{\mathrm{shock}}\rangle =\displaystyle \frac{{\displaystyle \sum }_{y,z}\left({x}_{{\rm{R}}}(y,z)-{x}_{{\rm{L}}}(y,z)\right)}{2{N}_{{yz}}},\end{eqnarray} \tag{ 7 }$$

where the subscripts R and L denote the two shock-front positions (the right and left sides of the layer), and N_yz denotes the total number of cells on the y–z plane. The factor of 2 in the denominator means that there are two shocks, and the overall mean position of the shock front is half of x_R − x_L at every given (y, z). We then measure the mean density within the shock-compressed layer $\langle n\rangle$ accordingly, similar to the gas density averaged over the entire volume of the shock-compressed layer defined above.

For a systematic study, we opt to vary two properties: the amplitude of the density fluctuation $\langle A\rangle$ and the spatial resolution Δx. We vary $\langle A\rangle$ such that the mean dispersion of the density ${\rm{\Delta }}{\rho }_{0}=\sqrt{\langle \delta {\rho }_{0}^{2}\rangle }/{\rho }_{0}$ spans 100%, 31.6%, 10%, 3.16%, 1%, 0.361%, and 0.1%. The 100% is motivated by the fact that density inhomogeneity exists on all spatial scales in the diffuse ISM (e.g., Armstrong et al. 1995; Lazarian & Pogosyan 2000; Chepurnov & Lazarian 2010), and the column density of H i gas likely varies with the order of unity (e.g., Burkhart et al. 2015; Fukui et al. 2018). We also investigate the conditions with extremely low levels of fluctuation (as low as 0.1%) to understand the dependence of ${\rm{\Delta }}{\rho }_{0}$. We choose the spatial resolution of Δx = 7.8 × 10⁻², 3.9 × 10⁻², 2.0 × 10⁻², and 9.8 × 10⁻³ pc by splitting the calculation domain from 256 × 256 × 128 cells at the lowest resolution up to 2048 × 1024 × 1024 cells at the highest resolution (hereafter noted as Δx = 0.08 pc, 0.04 pc, 0.02 pc, and 0.01 pc for simplicity in the text and figures). In addition, we repeat simulations with the same ${\rm{\Delta }}{\rho }_{0}$ and Δx but with three different random phases α₁, α₂, and α₃ to investigate the possible range over which the average properties vary due to such randomness even under the same statistical condition. Note that the highest resolution is computationally expensive and is limited to study three representative cases only (${\rm{\Delta }}{\rho }_{0}=31.6$, 10, and 3.16% with phase α₁).

The combination of these 3(+1) resolutions, 7 fluctuation amplitudes, and 3 phases corresponds to 66 cases. We additionally perform three controlled runs as a reference, where we investigate head-on collisions with the upstream WNM density completely uniform (${\rm{\Delta }}{\rho }_{0}=0 \%$). We thus investigate 69 cases in total. Table 1 summarizes the parameter sets, where we list 45 cases that we will present in this article, out of the total 69 cases. The calculation results are sampled every 0.1 Myr in each run.

Table 1.  Studied Parameters

Runs	Parameters
	Adiabatic (A)	Resolution	Δρ0	Phase	γ	Figures
	or Heating	Δx	$=\sqrt{\langle \delta {\rho }_{0}^{2}\rangle }{/\rho }_{0}$	α
	& Cooling (C)	(pc )	(% )
	or Effective index (E)
A2D0000	A	0.02	0	N/A	γ0 = 5/3	4
C2D0000	C	0.02	0	N/A	γ0 = 5/3	4 and 3
C2D1000	C	0.02	100	α1, α2, α3	γ0 = 5/3	3, 7, 8, 9, 10, 11, 13
C8D0316	C	0.08	31.6	α1, α2, α3	γ0 = 5/3	5, 6, 7
C4D0316	C	0.04	31.6	α1, α2, α3	γ0 = 5/3	5, 6, 7
C2D0316	C	0.02	31.6	α1, α2, α3	γ0 = 5/3	3, 4, 5, 6, 7, 8, 10, 11, 13
C1D0316	C	0.01	31.6	α1	γ0 = 5/3	5, 6, 7
C8D0100	C	0.08	10	α1, α2, α3	γ0 = 5/3	7
C4D0100	C	0.04	10	α1, α2, α3	γ0 = 5/3	7
C2D0100	C	0.02	10	α1, α2, α3	γ0 = 5/3	3, 7, 8, 10, 11, 13
C1D0100	C	0.01	10	α1	γ0 = 5/3	7
C8D0031	C	0.08	3.16	α1, α2, α3	γ0 = 5/3	5, 6, 7
C4D0031	C	0.04	3.16	α1, α2, α3	γ0 = 5/3	5, 6, 7
C2D0031	C	0.02	3.16	α1, α2, α3	γ0 = 5/3	3, 5, 6, 7, 8, 11
C1D0031	C	0.01	3.16	α1	γ0 = 5/3	5, 6, 7
C2D0010	C	0.02	1	α1, α2, α3	γ0 = 5/3	3, 7, 8, 9, 10, 11, 13
C2D0003	C	0.02	0.316	α1, α2, α3	γ0 = 5/3	3, 8, 11
C2D0001	C	0.02	0.1	α1, α2, α3	γ0 = 5/3	3, 7, 8, 11
E2D0000	E	0.02	0	N/A	γ0 = 5/3
					& ${\gamma }_{\mathrm{eff}}=0.863$	13

Note. Forty-five cases (out of total 69 cases) that we will present in this article. Details of individual parameters are described in Section 2. The run names are in the first column. Each column is as follows: the second column indicates whether we calculate an adiabatic fluid (indicated by "A") or include heating and cooling processes based on Equation (4) (indicated by "C"), or use the effective index ${\gamma }_{\mathrm{eff}}$ (indicated by "E"). Resolution shows the spatial resolution. Phases α₁, α₂, and α₃ represent three different random phases for initial density fluctuation (${\alpha }_{{k}_{y},{k}_{z}}$ in Equation (6)), where N/A represents no fluctuation (i.e., $\sqrt{\langle \delta {\rho }_{0}^{2}\rangle }=0.0$). γ lists the original specific heat γ₀ and the effective index ${\gamma }_{\mathrm{eff}}$. Figures list the corresponding figure numbers.

Download table as:  ASCIITypeset image

In this section, we provide a brief overview of our simulation results, showing how the geometrical structure of the shock-compressed layer varies with ${\rm{\Delta }}{\rho }_{0}$ and the cooling process. Figure 3 shows the gallery of the density slice from runs with ${\rm{\Delta }}{\rho }_{0}$ = 100% (run C2D1000), 31.6% (run C2D0316), 10% (run C2D0100), 3.16% (run C2D0031), 1% (run C2D0010), and 0% (run C2D0000) with phase α₁ and Δx = 0.02 pc. The shock-compressed layer is wider and wound more significantly with larger ${\rm{\Delta }}{\rho }_{0}$ due to the larger density fluctuation, whereas the layer is narrower and tends to be less deformed with smaller ${\rm{\Delta }}{\rho }_{0}$. Nonzero ${\rm{\Delta }}{\rho }_{0}$ provides the interaction between shocks and density inhomogeneity, which drives turbulence in the shock-compressed layer. In contrast, ${\rm{\Delta }}{\rho }_{0}=0 \%$ is an extreme condition where the flow becomes one dimensional (i.e., the layer maintains a completely straight geometry), and all the mass of the WNM flow cools down and eventually accretes onto the thin CNM sheet formed at the center.

Zoom In Zoom Out Reset image size

Figure 4 shows the time evolution of $\langle {x}_{\mathrm{shock}}\rangle$ from three runs: an adiabatic flow with ${\rm{\Delta }}{\rho }_{0}=0 \%$ (run A2D0000, "1D adiabatic"), a flow with cooling with ${\rm{\Delta }}{\rho }_{0}=0 \%$ (run C2D0000, "1D with cooling"), and a flow with cooling with ${\rm{\Delta }}{\rho }_{0}=31.6 \%$ with phase α₁ (run C2D0316, "3D with cooling"). As seen from the difference between runs A2D0000 and C2D0000, the cooling process transforms the WNM to the CNM and the shock-compressed layer becomes denser and narrower. As shown in run C2D0316, a nonzero ${\rm{\Delta }}{\rho }_{0}$ converts a fraction of the WNM to the CNM and drives some turbulence, which makes the shock-compressed layer less dense and wider compared with the uniform case (run C2D0000), but still significantly denser and narrower than the adiabatic case (run A2D0000). All of the shock-compressed layer with nonzero ${\rm{\Delta }}{\rho }_{0}$ in our simulations evolves somewhere between runs A2D0000 and C2D0000.

Zoom In Zoom Out Reset image size

We first investigate the convergence with respect to the spatial resolution Δx by varying it from 0.08 pc to 0.01 pc in all ${\rm{\Delta }}{\rho }_{0}$ cases as shown in Table 1. We find that the mean properties of the shock-compressed layer is converged with 0.02 pc spatial resolution in any ${\rm{\Delta }}{\rho }_{0}$ case, but the trend of convergence differs between the cases with larger ${\rm{\Delta }}{\rho }_{0}\gt 10 \%$ and with smaller ${\rm{\Delta }}{\rho }_{0}\leqslant 10 \%$. Therefore, in this section, we show the results mainly from ${\rm{\Delta }}{\rho }_{0}=31.6 \%$ and 3.16% as examples of large and small ${\rm{\Delta }}{\rho }_{0}$ and investigate how the physical properties of the shock-compressed layer converge.

Panels (a) and (b) in Figure 5 show the mass histogram on the P–n diagram at 3 Myr from ${\rm{\Delta }}{\rho }_{0}=31.6 \%$ (run C2D0316, panel (a)) and 3.16% (run C2D0031, panel (b)). Most of the mass resides in two regions: the shock-heated component (WNM and UNM at $n\simeq 2\,{\mathrm{cm}}^{-3}$) and the cooled component (CNM at $n\gt 100\,{\mathrm{cm}}^{-3}$). As shown here, the WNM density is typically 100 times more diffuse than that of CNM, and the volume of the shock-compressed layer is dominated by the WNM. Therefore, resolving the typical transition scale from the WNM into the CNM is crucial for the convergence in the mean properties of the shock-compressed layer, while resolving the CNM cooling length is important when we investigate the detailed structures of CNM clumps. We thus investigate the effective cooling length λ_cool,eff in the shock-compressed layer at 3 Myr. Panels (c) and (d) in Figure 5 show its histogram. Here, we calculate λ_cool,eff as ${C}_{{\rm{s}}}P/(\gamma -1)/({n}^{2}{\rm{\Lambda }}-n{\rm{\Gamma }}+P{{\rm{\nabla }}}_{\mu }{v}_{\mu })$ in each cell, where ${C}_{{\rm{s}}}$ and P are the sound speed and thermal pressure in each cell. The panels show that 0.02 pc corresponds to the typical spatial scale on which the dominant components in the cooling length transits from the UNM to CNM. We also find that the spatial resolution of 0.02 pc fully resolves the WNM cooling length, and resolves the cooling length of more than 99% (92%) of the UNM in volume (mass) when ${\rm{\Delta }}{\rho }_{0}=31.6 \%$ and more than 98% (80%) when ${\rm{\Delta }}{\rho }_{0}=3.16 \%$. We thus expect that the mean properties converge when Δx = 0.02 pc or higher, with which we can follow the dynamical condensation from WNM and UNM to CNM due to cooling.

Zoom In Zoom Out Reset image size

This characteristic scale of 0.02 pc is also consistent with the resolution requirement empirically suggested by previous numerical simulations. For example, Inoue & Omukai (2015) performed converging-flow calculations similar to our studies⁹ and suggested that more than 60 cells on the most frequent cooling length (defined as ${\lambda }_{\mathrm{cool}}={C}_{{\rm{s}}}P/(\gamma -1)/({n}^{2}{\rm{\Lambda }})$) are required to have convergence in the mass probability distribution function and CNM clump mass function. In our simulations, the mean cooling scale $\langle {\lambda }_{\mathrm{cool}}\rangle$ is peaked at the WNM and UNM components with 3.4 pc in ${\rm{\Delta }}{\rho }_{0}=31.6 \%$ and 2.2 pc in 3.16%. The requirement of more than 60 cells over one cooling scale corresponds to Δx < 0.037 pc, accordingly.

Panels (a)–(d) in Figure 6 show such convergence in the mean properties; the time evolution of $\langle {x}_{\mathrm{shock}}\rangle$ and $\langle n\rangle$ as a function of Δx with phase α₁, with panels (a) and (c) for ${\rm{\Delta }}{\rho }_{0}=31.6 \%$ and (b) and (d) for 3.16%. In both ${\rm{\Delta }}{\rho }_{0}$ cases, the linear growth of $\langle {x}_{\mathrm{shock}}\rangle$ (i.e., the constant speed of the shock propagation) corresponds to the quasi-steady $\langle n\rangle$. As already seen in Figure 3, the larger (smaller) ${\rm{\Delta }}{\rho }_{0}$ results in a wider (narrower) layer and smaller (denser) $\langle n\rangle$. The measured $\langle {x}_{\mathrm{shock}}\rangle$ and $\langle n\rangle$ at 3 Myr are listed in Table 2. These results suggest that $\langle {x}_{\mathrm{shock}}\rangle$ and $\langle n\rangle$ indeed show a convergence with Δx = 0.02 pc, with less than 17% difference between the results of Δx = 0.02 pc and 0.01 pc.

Zoom In Zoom Out Reset image size

Table 2.  Physical Properties at 3 Myr from ${\rm{\Delta }}{\rho }_{0}=31.6 \%$ (Runs C*D0316) and 3.16% (Runs C*D0031) with Phase α₁, with Respect to the Spatial Resolution Δx

	$\langle {x}_{\mathrm{shock}}\rangle$ [pc], $\langle n\rangle$ $[{\mathrm{cm}}^{-3}$], fWNM [%], fUNM [%], fCNM [%], nCNMpeak $({\mathrm{cm}}^{-3}$)
Δx (pc)	${\rm{\Delta }}{\rho }_{0}=31.6 \%$	${\rm{\Delta }}{\rho }_{0}=3.16 \%$
0.08	5.74, 6.61, 9.7, 52.0, 38.3, 100	3.02, 12.0, 8.4, 33.5, 58.1, 200
0.04	4.73, 7.93, 8.6, 44.2, 47.2, 158	2.40, 15.0, 5.9, 31.3, 62.8, 251
0.02	4.71, 7.99, 7.9, 45.3, 46.8, 158	1.69, 21.2, 4.2, 24.6, 71.2, 316
0.01	5.14, 7.38, 7.6, 49.9, 42.5, 158	1.45, 24.7, 3.3, 24.4, 72.3, 400

Download table as:  ASCIITypeset image

The trend of convergence, however, differs between large and small ${\rm{\Delta }}{\rho }_{0}$. For example, the $\langle {x}_{\mathrm{shock}}\rangle$ and $\langle n\rangle$ of runs C*D0316 do not vary monotonically along with Δx whereas those of runs C*D0031 do. In addition, the difference in $\langle {x}_{\mathrm{shock}}\rangle$ and $\langle n\rangle$ at 3 Myr between Δx = 0.08 pc and 0.01 pc is limited to a factor of 0.11 in runs C*D0316 but to a factor of 2.09 in runs C*D0031. Thus, overall, convergence with large ${\rm{\Delta }}{\rho }_{0}$ is nonmonotonic, and calculations with coarse resolutions of ≥0.02 pc show similar values in $\langle {x}_{\mathrm{shock}}\rangle$ and $\langle n\rangle$, whereas the convergence with small ${\rm{\Delta }}{\rho }_{0}$ is monotonic and stringently requires 0.02 pc as expected from the cooling length analysis above.

To understand these trends, we investigate the mass frequency as a function of density within the shock-compressed layer, which is shown in panels (a) and (b) of Figure 7 with $\langle n\rangle$ as vertical lines. We also measure the mass fraction of WNM (f_WNM), UNM (f_UNM), and CNM (f_CNM), and the peak density of the CNM n_CNMpeak, which are summarized on Table 2. The bimodality commonly appears in both ${\rm{\Delta }}{\rho }_{0}$ cases, which shows the shock-heated component peaked at n ∼ 2–4 ${\mathrm{cm}}^{-3}$, mostly corresponding to the WNM and UNM, and the cooled component peaked at n ∼ 100–400 ${\mathrm{cm}}^{-3}$, mostly corresponding to the CNM. The relative mass fraction, however, depends on ${\rm{\Delta }}{\rho }_{0}$. f_CNM remains ∼45% when ${\rm{\Delta }}{\rho }_{0}=31.6 \%$, and both the WNM+UNM and CNM contribute to $\langle n\rangle$, whereas f_CNM is ∼70% when ${\rm{\Delta }}{\rho }_{0}=3.16 \%$ and $\langle n\rangle$ depends more on the CNM component.¹⁰ For example, $\langle n\rangle$ with ${\rm{\Delta }}{\rho }_{0}=31.6 \%$ shifts nonmonotonically in accordance with the nonmonotonic change in the WNM fraction with Δx, resulting in a similar $\langle n\rangle$ between different Δx. In contrast, $\langle n\rangle$ with ${\rm{\Delta }}{\rho }_{0}=3.16 \%$ shifts monotonically by some factor as n_CNM,peak monotonically shifts with Δx (e.g., from Δx = 0.08 pc to 0.01 pc, $\langle n\rangle$ shifts from 12.0 to 24.7 ${\mathrm{cm}}^{-3}$ with n_CNM,peak from 200 to 400 ${\mathrm{cm}}^{-3}$).

Zoom In Zoom Out Reset image size

Two distinct behaviors depending on ${\rm{\Delta }}{\rho }_{0}$ also appear with different random phases α₁, α₂, and α₃ in the upstream density fluctuation. Panel (c) of Figure 7 shows the mass frequency variation due to those three random phases in each ${\rm{\Delta }}{\rho }_{0}$ under a fixed power spectrum of P_ρ(k) ∝ k^−11/3 and a fixed spatial resolution of Δx = 0.02 pc. This shows that the mass frequency has a larger variation by the random phases α in the larger ${\rm{\Delta }}{\rho }_{0}$ cases. To compare this variation due to different α with the variation due to Δx, panel (d) of Figure 7 shows the ratio of $\langle n\rangle ({\rm{\Delta }}x,\alpha )$ to $\langle n\rangle (0.01\,\mathrm{pc}$, α₁) in each ${\rm{\Delta }}{\rho }_{0}$. In the case of ${\rm{\Delta }}{\rho }_{0}=31.6 \%$, we found that the variation due to different α and that due to different Δx are comparable (i.e., the width of the gray shades is comparable to the scatter of points around unity). This suggests that resolving the density fluctuation with different ${\rm{\Delta }}x$ also intrinsically produces the variation in $\langle n\rangle$ comparable to the variation produced by different α. Therefore, $\langle n\rangle$ converges nonmonotonically and calculations with coarse resolutions of Δx > 0.02 pc practically provide similar values in $\langle {x}_{\mathrm{shock}}\rangle$ and $\langle n\rangle$ even though the physical convergence still requires Δx = 0.02 pc. The physical origin of such a variation due to different Δx is explained by the nonlinear nature of the system, which we will discuss in Section 4.2. In contrast, in the case of ${\rm{\Delta }}{\rho }_{0}\leqslant 10 \%$, we found that the variation due to different α is limited, and the effect of changing Δx directly appears as a monotonic change of $\langle n\rangle$, showing the convergence as predicted by the cooling length analysis.

Therefore, in the following sections, we will focus on the results based on Δx = 0.02 pc and highlight how the physical properties of the shock-compressed layer varies with ${\rm{\Delta }}{\rho }_{0}$ and the random phases.

We found that two distinct behaviors depending on ${\rm{\Delta }}{\rho }_{0}$ also appear in the geometry of the shock-compressed layer corresponding to large/small ${\rm{\Delta }}{\rho }_{0}$. The left panel of Figure 8 summarizes the time evolution of the average shock-front position $\langle {x}_{\mathrm{shock}}\rangle$ for each ${\rm{\Delta }}{\rho }_{0}$ with Δx = 0.02 pc. The shades in each ${\rm{\Delta }}{\rho }_{0}$ show the range of the maximum and minimum due to the three random phases α. This shows that the shock-compressed layer becomes wider (with larger variation with different α) when ${\rm{\Delta }}{\rho }_{0}$ is larger, which one can visually confirm also in the density snapshots already shown in Figure 3. In addition, we investigate the dispersion of $\langle {x}_{\mathrm{shock}}\rangle$, defined as

$$\begin{eqnarray}&&\sqrt{\langle \delta {x}_{\mathrm{shock}}^{2}\rangle }=\displaystyle \frac{{\displaystyle \sum }_{y,z}\sqrt{{\left({x}_{{\rm{R}}}(y,z)-{\bar{x}}_{{\rm{R}}}\right)}^{2}+{\left({x}_{{\rm{L}}}(y,z)-{\bar{x}}_{{\rm{L}}}\right)}^{2}}}{2{N}_{{yz}}},\end{eqnarray} \tag{ 8 }$$

where ${\bar{x}}_{{\rm{R}}}$ and ${\bar{x}}_{{\rm{L}}}$ are the mean position of the right and left shocks, respectively, and the results are summarized in the right panel of Figure 8. $\sqrt{\langle \delta {x}_{\mathrm{shock}}^{2}\rangle }$ gradually grows in time in all ${\rm{\Delta }}{\rho }_{0}$ cases, especially for density inhomogeneity with a larger ${\rm{\Delta }}{\rho }_{0}$ which more significantly deforms shock fronts whose amplitude reaches a fraction of 10 pc (the size of the cross section of the computational domain). Given that the Kolmogorov spectrum in the upstream density fluctuation has a larger power on larger scales (Equation (6)), we expect that the shock deformation is apparent on larger scales. The shock-front geometry of ${\rm{\Delta }}{\rho }_{0}\,=100 \%$, for example, indeed indicates the impact from the largest-scale mode of $| k| =2\pi /{L}_{y}$ (e.g., panel (a) of Figure 3). A smaller ${\rm{\Delta }}{\rho }_{0}$ introduces less deformation, and the shock-front geometry becomes closer to the completely straight front observed in ${\rm{\Delta }}{\rho }_{0}=0 \%$ (e.g., panels (d) and (e) of Figure 3).

Zoom In Zoom Out Reset image size

The two modes in these geometries, combined with the WNM:CNM mass fraction in Section 3.2, suggest that large-scale oblique shocks induced by significant shock deformation with larger ${\rm{\Delta }}{\rho }_{0}$ suppresses the energy dissipation at the shock fronts, driving stronger turbulence and forming less CNM in the shock-compressed layer, and vice versa for smaller ${\rm{\Delta }}{\rho }_{0}$ cases. We thus investigate the velocity and vorticity structure of the shock-compressed layer, and the results are shown in Figure 9. We here choose ${\rm{\Delta }}{\rho }_{0}=100 \%$ and 1% to clearly highlight the difference between large and small ${\rm{\Delta }}{\rho }_{0}$. There are indeed fast WNM flows that continue into the shock-compressed layer without significant vorticity generation when ${\rm{\Delta }}{\rho }_{0}$ is large (e.g., 3 pc ≲ y ≲ 7 pc in panels (a) and (c)), whereas most of the flow is well decelerated when ${\rm{\Delta }}{\rho }_{0}$ is small (in panels (b) and (d)). We also measure the density-weighted velocity dispersion $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ as a function of n and ${\rm{\Delta }}{\rho }_{0}$, which is shown in Figure 10. The bimodal behavior depending on ${\rm{\Delta }}{\rho }_{0}$ also appears in this $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ as we expected; the less-decelerated fast flow of the WNM in larger ${\rm{\Delta }}{\rho }_{0}$ cases drives a larger $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ with significant variation due to different α.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Note that the high-density structures predominantly contribute to the turbulent energy density $\rho \langle \delta {v}_{\mathrm{dw}}^{2}\rangle$ in any ${\rm{\Delta }}{\rho }_{0}$ cases because the overall n dependence of $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ is $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }\propto {n}^{-m}$ with m ≤ 0.5 in the range of n = 1–1000 ${\mathrm{cm}}^{-3}$. In addition, $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ is dominated by the x component because the converging flow has a directionality in x. We found that the x component of $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ amounts to >70% of the total $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ in all density ranges in any ${\rm{\Delta }}{\rho }_{0}$ case. Such anisotropic turbulence is reported even in magnetized converging-flow simulations when the magnetic fields are mostly parallel to the gas flow (e.g., Vázquez-Semadeni et al. 2007; Inoue & Inutsuka 2012; Iwasaki et al. 2019).

Also note that the overall decreasing trend of $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ with n is consistent with other simulations and observations. For example, the velocity dispersion ≥2 $\mathrm{km}\,{{\rm{s}}}^{-1}$ at n < 1 ${\mathrm{cm}}^{-3}$ is reported from recent observations of the WNM absorption feature (e.g., Patra et al. 2018). The velocity dispersion ≤2 $\mathrm{km}\,{{\rm{s}}}^{-1}$ at n > 10 ${\mathrm{cm}}^{-3}$ is reported in previous converging WNM flow studies without magnetic fields (e.g., Heitsch et al. 2006). Fukui et al. (2018) also report this decreasing trend of $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ with n, where they perform synthetic observations of H i line profiles in magnetohydrodynamics simulations and compare them with emission-absorption measurements along quasar lines of sight.¹¹

Finally, we investigate the energy partition in the shock-compressed layer and its dependence on ${\rm{\Delta }}{\rho }_{0}$. The left panel of Figure 11 shows the time evolution of $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}{\rangle }_{\mathrm{tot}}}$ as a function of ${\rm{\Delta }}{\rho }_{0}$. The subscript "tot" means that we take the average both over the entire volume of the shock-compressed layer and over the three random phases. This again clearly shows the bimodal behavior with respect to ${\rm{\Delta }}{\rho }_{0}$, as expected from Figure 10. In addition, we found that a small ${\rm{\Delta }}{\rho }_{0}\leqslant 10 \%$ creates the velocity dispersion of 2.0–3.2 $\mathrm{km}\,{{\rm{s}}}^{-1}$ even down to ${\rm{\Delta }}{\rho }_{0}=0.1 \%$. The value of ∼2 $\mathrm{km}\,{{\rm{s}}}^{-1}$ is consistent with the dispersion driven by the thermal instability alone (e.g., Koyama & Inutsuka 2002, 2006). This result indicates that shocks are always able to drive the velocity dispersion of ∼2 $\mathrm{km}\,{{\rm{s}}}^{-1}$ through the thermal instability even when they propagate though the ISM with almost uniform density.

Zoom In Zoom Out Reset image size

Based on this $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}{\rangle }_{\mathrm{tot}}}$, we measure the conversion rate of the upstream kinetic energy into the postshock turbulent and thermal energy as

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{turb}}(t)=\displaystyle \frac{{M}_{\mathrm{total}}(t)\langle \delta {v}_{\mathrm{dw}}^{2}(t){\rangle }_{\mathrm{tot}}/2}{\int {\dot{M}}_{\mathrm{total}}(t){V}_{\mathrm{in}}^{2}/2\,{dt}}=\displaystyle \frac{\langle \delta {v}_{\mathrm{dw}}^{2}(t){\rangle }_{\mathrm{tot}}}{{V}_{\mathrm{in}}^{2}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{th}}(t)=\displaystyle \frac{\int P({\boldsymbol{x}},t)d{\boldsymbol{x}}/(\gamma -1)}{\int {\dot{M}}_{\mathrm{total}}(t){V}_{\mathrm{in}}^{2}/2\,{dt}}=\displaystyle \frac{\langle P(t){\rangle }_{\mathrm{tot}}/(\gamma -1)}{\langle \rho (t){\rangle }_{\mathrm{tot}}{V}_{\mathrm{in}}^{2}/2}.\end{eqnarray} \tag{ 10 }$$

Here, ${\dot{M}}_{\mathrm{total}}(t)$ denotes the mass accretion rate onto the shock-compressed layer at time t, which can be approximated as ${\dot{M}}_{\mathrm{total}}(t)\simeq 2{\rho }_{0}{L}_{y}{L}_{z}({V}_{\mathrm{in}}+{V}_{\mathrm{shock}}(t))$. Note that ${V}_{\mathrm{in}}\,=20\,\mathrm{km}\,{{\rm{s}}}^{-1}=\mathrm{Const}.$ and we can assume ${V}_{\mathrm{shock}}=\mathrm{Const}.$ as a zeroth-order estimation based on the time evolution of $\langle {x}_{\mathrm{shock}}\rangle$. Therefore, we use ${M}_{\mathrm{total}}(t)=\int {\dot{M}}_{\mathrm{total}}(t){\rm{d}}t$ to obtain the second equality both in Equations (9) and (10).¹²

The right panel of Figure 11 shows _turb and _th as a function of ${\rm{\Delta }}{\rho }_{0}$ and time. The bimodal behavior depending on ${\rm{\Delta }}{\rho }_{0}$ again appears in _turb and _th individually, and appears also in the relative importance of _turb and _th; for example, at 3 Myr, the turbulent and thermal pressures equally support the shock-compressed layer when ${\rm{\Delta }}{\rho }_{0}\leqslant 10 \%$ whereas turbulence dominates when ${\rm{\Delta }}{\rho }_{0}\gt 10 \%$.

In small ${\rm{\Delta }}{\rho }_{0}$ cases, the shock-compressed layer is initially supported by the thermal pressure with limited turbulence. As CNM formation proceeds and turbulence develops (as seen in the left panel of Figure 11), the turbulent and thermal pressures start to equally support the shock-compressed layer. The sum of _turb and _th is limited to $\lt 10 \%$ of the injected kinetic energy due to the radiation energy loss (implemented as the source term $\rho { \mathcal L }$ in Equation (3)). In larger ${\rm{\Delta }}{\rho }_{0}$ cases, strong turbulence prevents the dynamical condensation by cooling and the following CNM formation. The shock-compressed layer is occupied by the low-density high-temperature WNM and UNM (see the P–n diagram of Figure 5), which keeps _th higher than that in smaller ${\rm{\Delta }}{\rho }_{0}$ cases. The turbulence dominantly supports the shock-compressed layer at 3 Myr and _turb reaches 12.4% at maximum. Such a high conversion rate is consistent with the one driven by the interaction between shocks and density inhomogeneity (e.g., Inoue et al. 2012, 2013; Iwasaki et al. 2019). Inoue et al. (2013) showed that the growth velocity of the Richtmyer–Meshkov instability (Richtmyer 1960; Nishihara et al. 2010) is able to account for the velocity dispersion of the turbulence.

These results indicate that the observed supersonic turbulence within molecular clouds originated not only from a weak turbulence of a few $\mathrm{km}\,{{\rm{s}}}^{-1}$ owing to the thermal instability, but also (or even overridden) by the interaction between shocks and density inhomogeneity because the ISM in reality has a density fluctuation close to ${\rm{\Delta }}{\rho }_{0}=100 \%$ (see also Section 4.1). This is also consistent with the conclusion of Inoue & Inutsuka (2012).

As already seen in Figure 7 and Table 2, the CNM mass fraction varies from 45% to 70% as ${\rm{\Delta }}{\rho }_{0}=0.1 \%$ to 100% at 3 Myr. There is also a time evolution such that f_CNM is initially 0 because we inject the WNM flows and gradually increase to the aforementioned values. Typically, after the typical cooling time of the injected WNM flow of ${\tau }_{\mathrm{cool}}=1.3\,\mathrm{Myr}$, $\langle {x}_{\mathrm{shock}}\rangle$ evolves almost linearly with time where the CNM mass fraction is constant (i.e., the mass conversion rate into the WNM and CNM per injected WNM is constant, and therefore the layer widens quasi-steadily). The time evolution of the mean density $\langle n\rangle$ reflects such evolution as shown in Figure 6.

We suggest that the mass fraction measured at 3 Myr is a typical value achieved in the ISM. Any volume of the ISM is typically swept up by at least one supersonic shock in every 1 Myr (see McKee & Ostriker 1977) due to supernovae. Thus, once the ISM reaches the quasi-steady state as observed at 3 Myr in our simulations, such frequent shock events presumably sustain that state. We also suggest that the density structure with ${\rm{\Delta }}{\rho }_{0}=100 \%$ is more common in the ISM because the density probability distribution function in solar neighborhood star-forming regions has a wide density range well fitted by a log-normal function (e.g., measured from dust, Schneider et al. 2013, 2016). Therefore, our results suggest that the typical CNM mass fraction is ∼50% in the typical ISM. Large-scale simulations overestimate the dense gas mass available for star formation if their subgrid models immediately convert 100% of the diffuse WNM into CNM and molecular gas within 3 Myr. It is left for future studies to numerically simulate another shock passage through the multiphase ISM whose CNM mass fraction is already ∼50% (see Inoue & Inutsuka 2012; Iwasaki et al. 2019).

The CNM formation in any converging-flow simulation has a chaotic behavior due to the nonlinear nature of the system, in a sense that a slight difference in the shock deformation later results in the formation of CNM clumps with slightly different masses and sizes, and this impacts subsequent CNM formation by altering the density/turbulent structures. For example, as shown in panels (d) and (e) of Figure 3, even a limited amplitude difference between ${\rm{\Delta }}{\rho }_{0}=3.16 \%$ and 1% already resulted in a different density distribution in the shock-compressed layer by 1.8 Myr. Such a chaotic behavior always occurs also within a given ${\rm{\Delta }}{\rho }_{0}$ between various α under a fixed Δx and between various Δx under a fixed α, which is able to introduce variation in the mean properties of the layer. Our results in Sections 3.2, 3.3, and 3.4 show that, when ${\rm{\Delta }}{\rho }_{0}$ is small, this nonlinear nature is less prominent if we average over the shock-compressed layer as $\langle {x}_{\mathrm{shock}}\rangle$ and $\langle n\rangle$. This is because >90% of the WNM flow kinetic energy is dissipated, and most of the injected mass condense into CNM with limited turbulence. Therefore, the monotonic convergence with Δx in the mean properties directly reflects whether or not we resolve the condensation process well enough.

In contrast, when ${\rm{\Delta }}{\rho }_{0}$ is large, such a chaotic behavior appears in the mean properties, which we found is primarily due to the interaction between shocks and CNM clumps. In a large ${\rm{\Delta }}{\rho }_{0}$ case, the layer becomes strongly turbulent where CNM clumps have a large velocity dispersion; for example, $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ at n = 100–1000 ${\mathrm{cm}}^{-3}$ reaches 2.7–4.0 $\mathrm{km}\,{{\rm{s}}}^{-1}$ when ${\rm{\Delta }}{\rho }_{0}=100 \%$, which is >5 times faster than the CNM sound speed (0.67 $\mathrm{km}\,{{\rm{s}}}^{-1}$ at the thermally balanced state of n = 100 ${\mathrm{cm}}^{-3}$ and T = 41 K). Once those fast CNM clumps form, they sometimes push/penetrate shock fronts (e.g., as seen at y ≃ 1 pc and 5 pc in panel (b) and y ≃ 1 pc in panel (c) of Figure 3). This process provides an additional deformation of shock fronts on top of the deformation driven by the upstream density inhomogeneity. This additional deformation impacts the following CNM formation and turbulence in the shock-compressed layer, which introduces another shock deformation subsequently. The mean properties (e.g., the CNM mass fraction) also vary accordingly. This process always prevents a perfect convergence, and our results in Sections 3.2, 3.3, and 3.4 suggest that such variation due to different Δx becomes comparable to the variation due to different α when ${\rm{\Delta }}{\rho }_{0}$ is large.

Up to now, we have ignored self-gravity in this article. Self-gravity decelerates the expansion of the shock-compressed layer, keeps CNM clumps in the layer, and pulls CNM clumps back to the layer even when they penetrate the shock fronts. Let us estimate two timescales determined by self-gravity and show the time range where the absence of self-gravity is valid.

First, we define a timescale after which the shock-compressed layer becomes self-gravitating rather than ram pressure confined. Let us label this timescale as ${t}_{\mathrm{sg}}$. ${t}_{\mathrm{sg}}$ can be estimated as the force balance between the self-gravity and the ram pressure of the converging flows as $\pi G{{\rm{\Sigma }}}^{2}(t)/2\gt {\rho }_{0}{V}_{\mathrm{in}}^{2}$, where G is the gravitational constant, and Σ(t) is the column density of the shock-compressed layer at time t. Given that ${\rm{\Sigma }}(t)\simeq 2{\rho }_{0}{V}_{\mathrm{in}}t$, ${t}_{\mathrm{sg}}$ can be estimated as

$$\begin{eqnarray}&&{t}_{\mathrm{sg}}\simeq \sqrt{\displaystyle \frac{1}{2\pi G{\rho }_{0}}}=44\,\mathrm{Myr}{\left(\displaystyle \frac{{n}_{0}}{0.57{\mathrm{cm}}^{-3}}\right)}^{-1/2}.\end{eqnarray} \tag{ 11 }$$

Second, there is a typical timescale over which the self-gravity of the shock-compressed layer pulls back CNM clumps that push/penetrate the shock fronts (see Appendix C in Iwasaki et al. 2019). Let us label this timescale as ${t}_{\mathrm{stop}}$, which is ${t}_{\mathrm{stop}}\,\simeq {\left({v}_{\mathrm{ej}}/4\pi G{\rho }_{0}{V}_{\mathrm{in}}\right)}^{1/2}$, where ${v}_{\mathrm{ej}}$ is the CNM clump's ejection velocity when they penetrate the shock fronts. ${t}_{\mathrm{stop}}$ can be estimated as

$$\begin{eqnarray}&&{t}_{\mathrm{stop}}\simeq 12\,\mathrm{Myr}{\left(\displaystyle \frac{{v}_{\mathrm{ej}}}{3\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{\tfrac{1}{2}}{\left(\displaystyle \frac{{n}_{0}}{0.57{\mathrm{cm}}^{-3}}\right)}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{{V}_{\mathrm{in}}}{20\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-\tfrac{1}{2}}.\end{eqnarray} \tag{ 12 }$$

Here we take 3 $\mathrm{km}\,{{\rm{s}}}^{-1}$ for ${v}_{\mathrm{ej}}$ based on the typical CNM $\sqrt{\langle \delta {v}_{\mathrm{dw}}^{2}\rangle }$ in our simulations with ${\rm{\Delta }}{\rho }_{0}=100 \%$ (see Figure 10).

Therefore, both timescales suggest that it is an acceptable assumption to ignore self-gravity when we focus on the early stages of multiphase ISM formation as we did in our simulation ≲3 Myr,

Thermal conduction plays an important role in determining the detailed structure of the CNM, especially the structure of the thin transition layer between the WNM and CNM (Field 1965). This is characterized by the Field length of the CNM. Under the cooling function and conduction rate (Equations (3) and (4)), the typical Field length of the CNM is as short as 10⁻³ to 10⁻⁴ pc (at T = 20 K and n = 150 cm⁻³; see Koyama & Inutsuka 2004), which we do not resolve in our current simulations. However, the thermal conduction does not appear to have a strong influence on the convergence in the macroscopic properties of the shock-compressed layer (Section 3.2). This is because the dynamics in our simulation is dominated more by the interaction between the shocks and upstream density inhomogeneity, and also by the dynamical condensation due to cooling, than by the thermal conduction alone. We describe our understanding of this situation below.

In a system where thermal conduction controls the dynamics, it is important to resolve the Field length. For example, starting from a thermal equilibrium UNM, Koyama & Inutsuka (2004) performed a one-dimensional numerical calculation to investigate the formation of the WNM and CNM and followed the long-term evolution of the motion between the two phases. They showed that thermal conduction drives motions of ∼0.1 $\mathrm{km}\,{{\rm{s}}}^{-1}$, and resolving the Field length is required to calculate this motion (by at least three cells—the Field condition). Iwasaki & Inutsuka (2014) investigated the two-dimensional cases, which also observed the ∼0.1 $\mathrm{km}\,{{\rm{s}}}^{-1}$ velocity dispersion driven by the thermal conduction and confirmed the requirement of the Field condition to achieve convergence in the velocity dispersion.

In contrast, in more dynamical systems like our simulations, supersonic shocks create the WNM and UNM with high pressure, which is far from thermal equilibrium. In this case, the cooling dominates the dynamical evolution of those phases (see, e.g., Appendix A of Iwasaki & Inutsuka 2012), and it is important to resolve the cooling length. Koyama & Inutsuka (2002) numerically investigated the evolution of the postshock medium and showed that the interaction between the shocks and upstream density inhomogeneity, as well as the following dynamical condensation of the WNM and UNM into CNM due to cooling, keep driving much faster turbulence of a few $\mathrm{km}\,{{\rm{s}}}^{-1}$. Hennebelle & Audit (2007) numerically investigated the effect of thermal conduction in their high-resolution, converging-flow simulations (with a 0.002 pc resolution albeit two dimensional) by changing the thermal conductivity, and they have shown that the total mass in the CNM does not significantly change with the thermal conductivity.

Therefore, thermal conduction does not play a critical role in our converging-flow calculations, and this seems to be the reason why we do not necessarily resolve the Field length of 10⁻³ pc in this case. Nevertheless, it is still important to resolve the cooling length of a few pc—10⁻² pc to achieve the convergence in the macroscopic properties of the shock-compressed layer (as shown in panels (c) and (d) of Figure 5 and discussions therein).

ISM models and star formation prescription below the spatial resolution is one of the challenges and uncertainties in large-scale simulations (e.g., evolution of the entire galactic disk and large-scale structure of the universe). To consistently calculate such sub-resolution-scale ISM evolution and star formation, semianalytical studies have been performed aiming to formulate an ISM effective equation of state (EoS) that describes the balance between the CNM formation by thermal instability and supernova feedback (e.g., Yepes et al. 1997; Springel & Hernquist 2003). There are also studies based on numerical simulations to model such an ISM effective EoS controlled by turbulence (Joung et al. 2009; Birnboim et al. 2015). In this section, we would like to propose a similar effective EoS on a ∼10 pc scale that approximates the multiphase ISM in the CNM formation epoch as a one-phase medium. This is in the form of $P\propto {\rho }^{{\gamma }_{\mathrm{eff}}}$, where we evaluate the effective index ${\gamma }_{\mathrm{eff}}$ based on the results of our converging-flow simulations.

The concept of our effective EoS is summarized in Figure 12. The ordinary Rankine–Hugoniot relations connect physical quantities across a shock front in a one-dimensional adiabatic flow. The density ratio of the post/preshock regions, r = ρ₂/ρ₁, is characterized as $r=(\gamma +1){{ \mathcal M }}^{2}/((\gamma -1){{ \mathcal M }}^{2}+2)$, where ${ \mathcal M }$ is the Mach number ${ \mathcal M }={V}_{\mathrm{in},\mathrm{pre}}/{C}_{{\rm{s}}}$ with ${V}_{\mathrm{in},\mathrm{pre}}$ as the flow speed in the shock-front rest frame, and γ represents the polytropic index of the fluid. This relation can be inverted as $\gamma =(2r-{{ \mathcal M }}^{2}(1+r))/({{ \mathcal M }}^{2}(1-r))$ to evaluate γ by measuring the density ratio r and ${ \mathcal M }$ of that fluid. As an analogy to such adiabatic shocks, we propose that measurement of r and ${ \mathcal M }$ of the multiphase ISM should also give an effective index ${\gamma }_{\mathrm{eff}}$, which approximates the multiphase ISM as a one-phase medium. Here, the one-phase approximation means that the adiabatic WNM with the EoS using ${\gamma }_{\mathrm{eff}}$ evolves while having its mean properties consistent with that of the multiphase ISM even without directly solving heating and cooling processes. Such an effective EoS based on our converging-flow simulations should be a relation of $P\propto {\rho }^{{\gamma }_{\mathrm{eff}}}$ that connects the initial state of the injected WNM and the simulated mean state of the shock-compressed layer (see Figure 12). Qualitatively speaking, when the mean density of the layer is low and the layer is geometrically wide, the multiphase ISM is stiff against a given ram pressure by the WNM inflow, and the corresponding ${\gamma }_{\mathrm{eff}}$ is expected to be large. When the mean density is high and the layer is geometrically narrow, the multiphase ISM is soft and the corresponding ${\gamma }_{\mathrm{eff}}$ is expected to be small.

Zoom In Zoom Out Reset image size

The converging-flow configuration is almost in the postshock rest frame whereas shock propagations in reality are mostly in the preshock rest frame. Therefore, we need to modify the original Rankine–Hugoniot relation as follows when evaluating ${\gamma }_{\mathrm{eff}}$ from our simulations:

$$\begin{eqnarray}&&{\gamma }_{\mathrm{eff}}=\displaystyle \frac{2{r}_{{\rm{d}}}-{{ \mathcal M }}_{\mathrm{shock}}^{2}(1+{r}_{{\rm{d}}})}{{{ \mathcal M }}_{\mathrm{shock}}^{2}(1-{r}_{{\rm{d}}})}.\end{eqnarray} \tag{ 13 }$$

Here, ${r}_{{\rm{d}}}$ is the effective density ratio, ${r}_{{\rm{d}}}=\langle \rho \rangle /{\rho }_{0}$, where $\langle \rho \rangle$ is the mean mass density in the shock-compressed layer. This $\langle \rho \rangle$ includes both the WNM and CNM because we aim to formulate an EoS representing the overall mean properties of the multiphase ISM. ${{ \mathcal M }}_{\mathrm{shock}}$ is the Mach number in the shock-front rest frame, and therefore, ${{ \mathcal M }}_{\mathrm{shock}}=({V}_{\mathrm{in}}+{V}_{\mathrm{shock}})/{C}_{{\rm{s}},\mathrm{WNM}}$, where ${V}_{\mathrm{shock}}$ is the shock propagation speed.

The left panel of Figure 13 shows the time evolution of ${\gamma }_{\mathrm{eff}}$. We measure the ${V}_{\mathrm{shock}}$ as (x_shock(t)−x_shock(t − 0.1 Myr))/0.1 Myr until ${\tau }_{\mathrm{cool}}$ (= 1.3 Myr) and perform a linear fit as $({x}_{\mathrm{shock}}(t)\,-{x}_{\mathrm{shock}}({\tau }_{\mathrm{cool}}))/(t-{\tau }_{\mathrm{cool}})$ after ${\tau }_{\mathrm{cool}}$. All ${\rm{\Delta }}{\rho }_{0}$ cases show a decreasing trend of ${\gamma }_{\mathrm{eff}}$ in time, which reflects the time evolution of $\langle n\rangle$ due to the CNM formation (i.e., the layer becomes softer by forming the CNM; see Figure 6). ${\gamma }_{\mathrm{eff}}$ tends to become constant after ${\tau }_{\mathrm{cool}}$ because the layer expansion becomes quasi-steady with an almost quasi-steady NM mass fraction (see the left panel of Figure 8). Larger/smaller ${\rm{\Delta }}{\rho }_{0}$ cases have lower/higher mean density, and ${\gamma }_{\mathrm{eff}}$ is stiffer/softer accordingly, as we expected.

Zoom In Zoom Out Reset image size

We found that ${\gamma }_{\mathrm{eff}}$ is softer than isothermal even in the stiffest case ${\rm{\Delta }}{\rho }_{0}=100 \%$, and that it becomes as soft as ${\gamma }_{\mathrm{eff}}\sim 0.7$ in smaller ${\rm{\Delta }}{\rho }_{0}$ conditions. Such a high compressibility indicates that the overall dynamics on <10 pc scales in large-scale simulations could be further improved by introducing this type of effective EoS, especially in regions where multiphase ISM formation is ongoing without any feedback, because most current simulations use ${\gamma }_{\mathrm{eff}}\gt 1$ as a subgrid model (c.f., Inoue & Yoshida 2019).

As the first step toward such an actual application, we perform a converging adiabatic WNM flow by employing ${\gamma }_{\mathrm{eff}}$ to demonstrate how well our ${\gamma }_{\mathrm{eff}}$ reproduces the properties of the multiphase ISM. In this demonstration, we directly update the thermal pressure based on the density evolution and ${\gamma }_{\mathrm{eff}}$ at each time step, instead of calculating the heating and cooling ($\rho { \mathcal L }$ in Equation (3)). For simplicity, we opt to solve equations only in the x direction to set this demonstration of ${\gamma }_{\mathrm{eff}}$ in the one-dimensional case and set a constant value of ${\gamma }_{\mathrm{eff}}=0.863$ without modeling any time evolution seen in the left panel of Figure 13. This value of 0.863 is based on the result of run C2D1000 with phase α₁ at 3 Myr. We also employ ${\rm{\Delta }}{\rho }_{0}=0$ to simulate a completely uniform head-on collision, just as a simple demonstration.

The right panel of Figure 13 shows the time evolution of $\langle {x}_{\mathrm{shock}}\rangle$. "3D with cooling" shows the result of our simulations with ${\rm{\Delta }}{\rho }_{0}=100 \%$ (run C2D1000 with phase α₁), whereas "1D with ${\gamma }_{\mathrm{eff}}$" shows the result of this the test demonstration using ${\gamma }_{\mathrm{eff}}=0.863$ (run E2D0000). Table 3 summarizes the measured properties at 3 Myr. Our constant ${\gamma }_{\mathrm{eff}}$ does not reproduce the detailed time evolution of "3D with cooling." Nevertheless, it successfully reproduces ${V}_{\mathrm{in}}$ + ${V}_{\mathrm{shock}}$, $\langle {x}_{\mathrm{shock}}\rangle$, and $\langle n\rangle$ at 3 Myr within an 11% difference, which is still smaller than the variation due to different α; for example, the variation of $\langle {x}_{\mathrm{shock}}\rangle$ due to phases α₁, α₂, and α₃ at 3 Myr is 1.763 pc in run C2D1000 (∼29% variation against $\langle {x}_{\mathrm{shock}}\rangle =6.183$ pc with phase α₁; see the shade of ${\rm{\Delta }}{\rho }_{0}=100 \%$ in the left panel of Figure 8).

Table 3.  Measured Properties at 3 Myr

Runs	${V}_{\mathrm{in}}$ + ${V}_{\mathrm{shock}}$	$\langle {x}_{\mathrm{shock}}\rangle$	$\langle n\rangle$
	(km s−1)	(pc)	(cm−3)
3D with cooling (C2D1000)	21.276	6.183	6.750
1D with ${\gamma }_{\mathrm{eff}}$ (E2D0000)	22.067	6.338	6.081

Note. ${V}_{\mathrm{in}}$ + ${V}_{\mathrm{shock}}$, $\langle {x}_{\mathrm{shock}}\rangle$, and $\langle n\rangle$ measured in runs C2D1000 and E2D0000.

Download table as:  ASCIITypeset image

Ideally, we would like to provide a time-evolving model of ${\gamma }_{\mathrm{eff}}(t)$ along with a time-evolving model of f_CNM, which, however, we reserve for future studies at this moment. In such studies, we should also consider the intrinsic variation of ${\gamma }_{\mathrm{eff}}$ due to random phases, for example, ${\gamma }_{\mathrm{eff}}=0.829$, 0.863, and 0.898 in run C2D1000 with phases α₃, α₁, and α₂ at 3 Myr (see the left panel of Figure 13).

Note that the implementation of this effective EoS is unfortunately not that straightforward. For example, to introduce a time-evolving ${\gamma }_{\mathrm{eff}}$, we would like to measure shock propagation speed and elapsed time since the last shock passage even in large-scale simulations. This is, however, computationally expensive and time consuming, similar to following the stellar population evolution below the spatial resolution to blow supernovae at the correct timing. Such difficulties have to be also discussed and left for future studies. Nevertheless, because the typical frequency of shock passages is as high as once per megayear (e.g., McKee & Ostriker 1977; see also Section 4.6), we may expect that the quasi-steady ${\gamma }_{\mathrm{eff}}$ at ∼3 Myr is still close to the typical time-averaged state of the actual ISM.

In this section, we briefly address some potential limitations/issues in converging-flow simulations (not only ours but also in general). Given that 1 Myr is the typical interval in the ISM between successive shock passages by multiple supernovae and/or H ii regions (McKee & Ostriker 1977; Inutsuka et al. 2015), supersonic flow in reality continue in a fixed direction only ≤1 Myr. Statistically speaking, successive flows essentially propagate from any direction, and they are incident on the shock-compressed layer at some angle. Almost all of the converging-flow studies therefore presumably keep the flow injection too long in a fixed direction. There have been previous studies introducing an inclination angle between flows to investigate the effect of magnetic diffusion and supercritical core formation (Körtgen & Banerjee 2015) and to investigate the reorientation of preexisting filaments (see Fogerty et al. 2016, 2017), but it is still left for future studies to reveal how the mean properties of the shock-compressed layer depend on such inclinations and multiple compressions by flows from various angles.

Similarly, the typical dynamical timescale of the shock-compressed layer is a few megayears (e.g., the crossing time of the WNM component over the shock-compressed layer is 10 pc/$10\,\mathrm{km}\,{{\rm{s}}}^{-1}$ in run C2D1000). It is thus also left for future studies to investigate how an already-created shock-compressed layer expands and/or shrinks once the inflow ceases and to measure whether or not the turbulence and f_CNM are maintained. Limiting the inflow mass is one of the techniques to study such a condition; Vázquez-Semadeni et al. (2007), for example, demonstrate that the global and local collapses of the shock-compressed layers occur once all of the gas finish accreting onto the shock-compressed layers.

The concept of a converging-flow setup is to easily perform calculations in the postshock rest frame. However, most of the postshock regions in reality are presumably not sandwiched by two shock fronts as in the converging flows, but instead by one shock front and one contact discontinuity. Converging flow is just an analog of such a shock-contact discontinuity system, and only one-half of the shock-compressed layer is meaningful. The two-shock-front configuration likely impacts the time evolution of the shock-compressed layer. For example, as shown in Figure 9, fast WNM flows continue deep into the layer and interact with each other, especially when ${\rm{\Delta }}{\rho }_{0}$ is large, which depends on the shock-front geometry on both sides. The turbulent properties and f_CNM may accordingly differ in a shock-contact discontinuity system. Simulation studies of a shock-contact discontinuity system (e.g., Koyama & Inutsuka 2002) are still limited, and careful comparison is left for future studies.

We also ignore magnetic fields for simplicity in this article, but they play a pivotal role; for example, magnetic field pressure supports the shock-compressed layer and turbulence decays, especially in the case where field lines are oriented perpendicular to the inflow (Heitsch et al. 2009; Vázquez-Semadeni et al. 2011; Inoue & Inutsuka 2012; Iwasaki et al. 2019). We expect that magnetic fields do not modify our proposed ${\gamma }_{\mathrm{eff}}$ significantly because the magnetized shock-compressed layer tends to be equipartitioned (see Iwasaki et al. 2019, for ${V}_{\mathrm{in}}\lt 20\,\mathrm{km}\,{{\rm{s}}}^{-1}$ cases), but this still has to be investigated with magnetized converging-flow simulations.