Metallicity Dependence of Molecular Cloud Hierarchical Structure at Early Evolutionary Stages

To investigate the development of thermal instability and the formation of multiphase ISM from the WNM, we calculate supersonic WNM converging flows by solving the following basic equations:

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {v}_{i})}{\partial t}+{{\rm{\nabla }}}_{j}({T}_{{ij}}+\rho {v}_{j}{v}_{i})=-\rho {{\rm{\nabla }}}_{i}{\rm{\Phi }},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{T}_{{ij}}=\left(P+\displaystyle \frac{{B}^{2}}{8\pi }\right){\delta }_{{ij}}-\displaystyle \frac{{B}_{i}{B}_{j}}{4\pi },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial e}{\partial t}+{{\rm{\nabla }}}_{i}\left[(e{\delta }_{{ij}}+{T}_{{ij}}){v}_{j}\right]\\ \qquad =\,{{\rm{\nabla }}}_{i}\left[\kappa (T){{\rm{\nabla }}}_{i}T\right]-\rho {v}_{i}{{\rm{\nabla }}}_{i}{\rm{\Phi }}-\rho { \mathcal L }(T,Z),\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{i}}{\partial t}+{{\rm{\nabla }}}_{j}\left({v}_{j}{B}_{i}-{v}_{i}{B}_{j}\right)=0,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}=4\pi {\rm{G}}\rho ,\end{eqnarray} \tag{ 6 }$$

where ρ is the mass density, v represents the velocity, P represents the thermal pressure, T without any subscript is the temperature, Φ is the gravitational potential, G is the gravitational constant, and ∇_i = ∂/∂x_i , where x_i spans x, y, and z. δ_ij is the identity matrix. We calculate the total energy density, e, as e = P/(γ − 1) + ρ v²/2 + B²/8π where the ratio of the specific heat is γ = 5/3. We introduce the thermal conductivity, κ, as κ(T) = 2.5 × 10³ T^0.5 erg cm⁻¹ s⁻¹ K⁻¹, which considers collision between hydrogen atoms (Parker 1953).

${ \mathcal L }(T,Z)$ is the net cooling rate. We employ the functional form of ${ \mathcal L }(T,1.0\,{Z}_{\odot })$ from Kobayashi et al. (2020) in the case of 1.0 Z_⊙ condition. This functional form combines the results from Koyama & Inutsuka (2002) in T ≤ 14,577 K and Cox & Tucker (1969) and Dalgarno & McCray (1972) in T > 14, 577 K by considering the cooling rates due to Lyα, C_II, He, C, O, N, Ne, Si, Fe, and Mg lines with the photoelectric heating. The shock heating and compression destabilize the WNM to join the UNM (defined as ${(\partial ({ \mathcal L }/T)/\partial T)}_{P}\lt 0;$ see, e.g., Balbus 1986, 1995), leading to the formation of the CNM clumps.

To investigate the lower-metallicity cases in this article, we apply three modifications to ${ \mathcal L }(T,1.0\,{Z}_{\odot })$ to prepare ${ \mathcal L }(T,Z)$. First, the cooling rate due to the metal lines is set to be linearly scaled with the metallicity, which is a good approximation in >10⁻⁴ Z_⊙ as long as the metal lines dominate the cooling (Inoue & Omukai 2015). Second, we set the photoelectric heating rate to be linearly scaled with the metallicity by assuming that the dust abundance is also scaled with the metallicity. We, therefore, use ${{\rm{\Gamma }}}_{\mathrm{pe}}=2.0\times {10}^{-26}\left(Z/{Z}_{\odot }\right)$ erg s⁻¹. Third, we implement the X-ray and cosmic-ray heating rates as Γ_X = 2.0 × 10⁻²⁷ erg s⁻¹ and Γ_CR = 8.0 × 10⁻²⁸ erg s⁻¹, respectively (Koyama & Inutsuka 2000). The X-ray and cosmic-ray heating processes are subdominant in the metallicity range of this study; for example, X-rays (cosmic rays) contribute 10% (30%) of the total heating rates in the 0.2 Z_⊙ environment. Their relative importance decreases even more when the metallicity increases because the photoelectric heating rate increases with metallicity. We show the detailed functional form of this revised ${ \mathcal L }(T)$ in Appendix A. Figure 1 shows the thermal equilibrium state, i.e., ${ \mathcal L }(T,Z)=0$. This shows that, in lower-metallicity environments, the combination of the inefficient cooling and metallicity-independent X-ray and cosmic-ray heatings allows for the existence of the WNM phase until higher density in the range of 1–10 cm⁻³. In each calculation, we assume a uniform metallicity distribution in space, as representation of low-metallicity environments.

Zoom In Zoom Out Reset image size

We use the publicly available MHD simulation code Athena++ (Stone et al. 2020) to solve the basic equations, where we employ the HLLD MHD Riemann solver (Miyoshi & Kusano 2005) and the constrained transport method to integrate the magnetic fields (Evans & Hawley 1988; Gardiner & Stone 2005, 2008). The self-gravity is calculated with the full multigrid method (Tomida et al. 2023).

Our simulation domain has sizes of L_x,y,z = 20, 10, and 10 pc on each side in Cartesian coordinates. We continuously inject supersonic WNM flows from the two x-boundaries so that the two flows collide head-on. We employ the periodic boundary condition on the y- and z-boundaries. The collision forms a shock-compressed layer sandwiched by two shock fronts, in which the thermal instability converts the WNM to CNM. We employ V_inflow = 20 km s⁻¹ as the flow velocity; this choice corresponds to a representation of the late phase of supernova remnant expansion or normal component of galactic spiral shocks. The initial velocity field is set as ${v}_{x}(x)\,={V}_{\mathrm{inflow}}\tanh (-x/0.78\,\mathrm{pc})$ and v_y,z = 0 km s⁻¹.

The initial WNM flow, and also the injected WNM flow, are thermally stable, with a mean number density n₀ = 0.57 cm⁻³. The corresponding pressure, temperature, and sound speed at each metallicity are listed on Table 1, where we use ρ₀ = n₀ μ_M m_p with the mean molecular weight μ_M of 1.27 (Inoue & Inutsuka 2012). The ram pressure of the converging flow is ${P}_{\mathrm{ram}}/{k}_{{\rm{B}}}^{}=3.5\times {10}^{4}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$. The WNM flows have density fluctuation following the Kolmogorov spectrum P_ρ (k) ∝ k^−11/3 (Kolmogorov 1941; Armstrong et al. 1995). We impose a random phase in each k up to k/2π = 3.2 pc⁻¹, which corresponds to a wavelength of 0.32 pc. The mean dispersion of the density fluctuation is chosen as $\sqrt{\langle \delta {n}_{0}^{2}\rangle }/{n}_{0}=0.5$. The injected flows have periodic distributions smoothly connected to the initial condition, so that we impose ρ(t, x = − 10 pc, y, z) = ρ(t = 0, x = 10 pc − V_inflow t, y, z) and ρ(t, x = 10 pc, y, z) = ρ(t = 0, x = − 10 pc + V_inflow t, y, z) on the x-boundaries where t represents the time. The interaction between this density inhomogeneity and shock fronts generates turbulence (e.g., Koyama & Inutsuka 2002; Inoue & Inutsuka 2012; Carroll-Nellenback et al. 2014; Kobayashi et al. 2022).

Table 1. Parameters at Different Metallicities

	${P}_{0}/{k}_{{\rm{B}}}^{}$	T0	Cs,0	tcool	tfinal
	[K cm−3]	[K]	[km s−1]	[Myr]	[Myr]
1 Z⊙	3.6 × 103	6.4 × 103	6.4	1.0	3.0
0.5 Z⊙	3.5 × 103	6.2 × 103	6.3	2.0	6.0
0.2 Z⊙	3.4 × 103	6.1 × 103	6.2	5.1	15.0

Note. n₀ = 0.57 cm⁻³ and B₀ = 1 μG Are Common at All Metallicities.

Download table as:  ASCIITypeset image

The magnetic field is initially threaded in the x-direction with 1 μG strength. Previous studies show that successful molecular cloud formation occurs in such a configuration where the flow and the mean magnetic field are close to parallel (Hennebelle & Pérault 2000; Inoue & Inutsuka 2008; Iwasaki et al. 2019). The typical mean field strength in H i gas varies from 1–10 μG in the Milky Way (Crutcher 2012) and 0.5–5 μG in the Magellanic Clouds (Gaensler et al. 2005; Beck 2013; Livingston et al. 2022). We, therefore, opt to choose 1 μG as representative strength in investigating the metallicity dependence. The detailed dependence on the field strength in low-metallicity environments is left for future studies at this moment.

We employ a uniform spatial resolution of 0.02 pc to resolve the typical cooling length of the UNM. This resolution is required to have convergence in the CNM mass fraction after 1.0 t_cool (where t_cool is the cooling time defined as $P/(\gamma -1)/\rho { \mathcal L }$; Kobayashi et al. 2020). We identify the shock front position by P > 1.3P₀ to define the volume of the shock-compressed layer.

Figure 2 shows the three-dimensional view of the initial density field with a uniform 1 μG magnetic field.

Zoom In Zoom Out Reset image size

Figure 1 shows that the CNM thermal state at n > 10² cm⁻³ is similar within the 1.0–0.2 Z_⊙ range. In this metallicity range, the cooling rate is proportional to the metallicity (see Inoue & Omukai 2015). These suggest that the t_cool is longer in low-metallicity environments as ∝ Z⁻¹ but CNM properties may become similar between 1.0–0.2 Z_⊙ at the same time measured in units of the cooling time.

The typical t_cool of the injected WNM is ∼1 Myr at 1 Z_⊙, ∼2 Myr at 0.5 Z_⊙, and ∼5 Myr at 0.2 Z_⊙ (see Appendix B.1). Kobayashi et al. (2020) investigated the converging flow at 1 Z_⊙ and suggested that the shock-compressed layer achieves a quasi-steady state at ∼3t_cool. Therefore, to make a comparison between different metallicities both at the same physical time and at the same time measured in units of the cooling time, we integrate until 3 t_cool in each metallicity. Table 1 lists these parameters.

Figure 3 shows examples of the three-dimensional view of our simulation results. Panels (a), (b), and (c) compare the results at 3 Myr from the 1.0 Z_⊙, 0.5 Z_⊙, and 0.2 Z_⊙ runs. Panels (d) and (e) compare the results at 3t_cool from the 0.5 Z_⊙ and 0.2 Z_⊙ runs. These panels show that CNM clumpy/filamentary structures form more slowly in lower-metallicity environments. The development of such CNM structures similar to those in 1 Z_⊙ environment requires similar times measured in units of the cooling time, for example, at 3t_cool. Compared to the same physical time of 3 Myr, the geometry of the shock-compressed layer in a lower-metallicity environment is less disturbed (e.g., panel (c)), close to a plane-parallel configuration because inefficient cooling keeps the layer more adiabatic.

Zoom In Zoom Out Reset image size

Figure 4 compares the thermal states of the different metallicities at 3 Myr and at 3 t_cool. This shows that, at 3 Myr, the shock-compressed layer is still dominated more by the shock-heated WNM/UNM in the lower-metallicity environment. Their thermal pressure is still close to the flow ram pressure. Therefore, the temperature just starts to decrease in the shock-compressed layer of 0.2 Z_⊙ at 3 Myr (0.6 t_cool), so that the plane-parallel geometry of the shock-compressed layer seen in Figure 3 is close to a simple one-dimensional shock compression (see also Section 3.3 and Appendix C for its impact on the turbulent velocity).

Zoom In Zoom Out Reset image size

The net magnetic flux in our simulation does not increase in time because the mean magnetic field is completely parallel to the WNM flow. Nevertheless as we see in Figure 3, the pre-shock density fluctuation induces a number of oblique shocks at the shock front, which locally fold the field lines. The magnetic fields in the shock-compressed layer are further twisted and stretched due to the turbulence, which introduces a larger scatter in the local field strength. We investigate the relation between the field strength and the number density at 3 t_cool. Figure 5 shows the phase histogram on the plane of the magnetic field strength and the number density in the shock-compressed layer. This figure shows that the field strength is initially amplified toward the shock-heated WNM volume through the shock compression almost following B ∝ n¹ (as indicated with the translucent gray line). The strength further varies due to the turbulence at n < 10 cm⁻³. Note that the initial mean magnetic field is completely parallel to the flow velocity, and not all of the volume experiences a perfect one-dimensional shock compression. Therefore, the most frequent field strength is slightly weaker than that the relation B ∝ n¹ (see the upper panels of Figure 5).

Zoom In Zoom Out Reset image size

Accompanying the density enhancement toward n ∼ 10³ cm⁻³ by the thermal instability, the field strength gradually increases but with a limited level. The black solid curve in Figure 5 shows the volume-weighted average of the magnetic field strength as a function of the number density 〈B〉(n), where $B=\left|{\bf{B}}.\right|$. ¹⁰ This shows that the amplification roughly scales with 〈B〉 ∝ n^1/5 in the n > 1 cm⁻³ range. This slow evolution occurs because condensing motion by the thermal instability is confined along the magnetic field orientations. Similar results have also been obtained by previous numerical studies of the solar metallicity environment (e.g., Hennebelle et al. 2008; Heitsch et al. 2009; Inoue & Inutsuka 2012). Our results show that this gradual amplification of the magnetic field occurs also in lower-metallicity environments.

Figure 5 indicates that there is a maximum magnetic field strength attained. We can estimate this maximum strength by assuming there is balance between the postshock magnetic pressure and the ram pressure of the injected WNM as ${B}_{\mathrm{eq}}^{2}/8\pi ={\rho }_{0}{V}_{\mathrm{inflow}}^{2}$. This gives a typical value of

$$\begin{eqnarray}&&{B}_{\mathrm{eq}}=11\,\mu {\rm{G}}{\left(\displaystyle \frac{{n}_{0}}{0.57\,{\mathrm{cm}}^{-3}}\right)}^{1/2}\left(\displaystyle \frac{{V}_{\mathrm{inflow}}}{20\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right).\end{eqnarray} \tag{ 7 }$$

Although this is an extreme case where the magnetic energy dominates the shock-compressed layer, Equation (7) should give a good approximation even when we consider local enhancement by the turbulence and the condensing motion by the thermal instability, because the inflow ram pressure determines the typical maximum pressure of the shock-compressed layer. B_eq, shown in Figure 5 (the horizontal black dashed line), indeed outlines the typical maximum strength of the magnetic fields.

Note that the field strength can increase beyond B_eq at the densest volume where the self-gravity plays a role. This occurs in n ≳ 2.6 × 10³ cm⁻³ in our simulation setup. Such a critical density of n ∼ 2.6 × 10³ cm⁻³ can be estimated as the density at which the CNM plasma beta with B = B_eq becomes unity (Iwasaki & Tomida 2022). It is difficult to clearly confirm this amplification in our current simulations due to our limited volume and due to our diffuse WNM-only initial condition. Nevertheless, we expect that the field strength also increases in low-metallicity environments because the CNM thermal states are similar between 1.0–0.2 Z_⊙, as are the critical densities at which the self-gravity starts to dominate. Previous simulations of a converging flow at 1.0 Z_⊙ show 〈B〉 ∝ n^1/2 at n > 10³ cm⁻³ either when they integrated much longer times to accumulate mass or when they started with two-phase atomic flows with their mean number density already n ≳ 5 cm⁻³ (e.g., Hennebelle et al. 2008; Inoue & Inutsuka 2012; Iwasaki & Tomida 2022).

Our results show that the development of CNM clumpy/filamentary structure and the magnetic field strength is similar between metallicities if they are measured at the same time in units of the cooling time in each metallicity (i.e., the evolution slows down linearly with the metallicity in terms of the physical time). The field strength is relatively constant from a few microgauss to 10 μG up to n < 10³ cm⁻³ as 〈B〉 ∝ n^1/5, and possibly starts to amplify further by self-gravity. Note that this field strength in n < 10³ cm⁻³ is consistent with Zeeman measurements of low (column) density regions of molecular clouds in the Milky Way (e.g., Crutcher 2012) and Faraday rotation measurements toward the ISM in the LMC and SMC (Gaensler et al. 2005; Beck 2013; Livingston et al. 2022).

In Kobayashi et al. (2022), we showed that the coevolution of the turbulence and the thermal evolution by the thermal instability plays a significant role to determine the multiphase density structure and its density probability distribution function at the molecular cloud formation stage. Kobayashi et al. (2022), however, neglected the magnetic fields and studied simple hydrodynamic converging flows. In this section, we will confirm those findings even in the current MHD simulations and investigate how the turbulent structure differs between/resembles the cases with different metallicities.

Panels (a), (b), and (c) of Figure 6 show the evolution of the velocity dispersion, the mean density, and the turbulent pressure, respectively, as a function of the time measured in units of the cooling time (i.e., t/t_cool). In the early stage at t < 0.5 t_cool, the shock-compressed layer has plane-parallel geometry in the low-metallicity environment (see Section 3.2). Due to the resultant efficient deceleration of the inflow and the inefficient cooling, the mass is dominated by the shock-heated WNM state with slow turbulence, with negligible volume/mass of the CNM. Therefore, the volume-weighted velocity dispersion $\sqrt{\langle \delta {v}^{2}\rangle }$ is close to the density-weighted velocity dispersion $\sqrt{\langle \delta {v}^{2}{\rangle }_{\rho }}$ of ∼2 km s⁻¹ in 0.2 Z_⊙ (panel (a) of Figure 6; see also Appendix C).

Zoom In Zoom Out Reset image size

At t ≳ 0.5 t_cool in the lower-metallicity environment, the phase transition from the shock-heated WNM to the CNM tends to occur at a pressure closer to the inflow ram pressure (see panel (c) Figure 4). The mean number density is higher in lower-metallicity environments accordingly (panel (b) of Figure 6). During this evolution, the turbulent pressure gradually grows until it balances against the inflow ram pressure. As a result, the turbulent pressure is almost the same between the three metallicities at 3 t_cool (panel (c) of Figure 6).

These results suggest that the turbulent pressure is almost the same at 3 t_cool. The difference exists due to efficient compression in the lower-metallicity environments, because the inflow and the shock front incidents with an almost 90° angle, which results in a denser postshock WNM with slower velocity. However, this difference is limited by a factor of 2 (4) in the velocity dispersion (in the mean density, respectively). The value of $\sqrt{\langle \delta {v}^{2}\rangle }\sim 6$–9 km s⁻¹, close to the sound speed of the WNM, suggests that turbulence on the molecular cloud scale is powered by the WNM super-Alfvénic turbulence not only in the Milky Way galaxy but also in the LMC and SMC. ¹¹

To understand the turbulence structure, we perform a Fourier analysis of the turbulence by using a fast Fourier transform in the west (3.3; Frigo & Johnson 2005), and decompose the turbulence into the solenoidal and compressive modes, which are defined, respectively, as

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{v}}}}_{\mathrm{sol}}(\hat{k})=\left(\hat{k}\times \tilde{{\boldsymbol{v}}}\right)\times \hat{k},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{v}}}}_{\mathrm{comp}}(\hat{k})=\left(\hat{k}\cdot \tilde{{\boldsymbol{v}}}\right)\hat{k}.\end{eqnarray} \tag{ 9 }$$

Here, $\hat{k}$ represents a unit wavevector and $\tilde{{\boldsymbol{v}}}$ represents the Fourier component of the velocity field. We here denote $k=\left|{\boldsymbol{k}}\right|$. Figure 7 shows the power spectrum at 3 t_cool. The solenoidal mode dominates the turbulence power on all scales. Table 2 summarizes the total fraction of each turbulent mode at 3 t_cool. The solenoidal mode fraction amounts to 80%–90%. In calculating this, we employ powers between k/2π = 0.1 and 2.56 pc⁻¹ to avoid numerical diffusion effects on a small scale (we refer the reader to the caption of Figure 7). ¹²

Zoom In Zoom Out Reset image size

Figure 8 shows the time evolution of this mode fraction as a function of t/t_cool. The solenoidal (compressive) mode fraction evolves quasi-steadily with ∼80% (20%, respectively) after 1 t_cool.

Zoom In Zoom Out Reset image size

Such solenoidal-mode-dominated turbulence is a natural consequence of the accreting system followed with the thermal instability, as we discuss in our previous converging flow simulations at 1.0 Z_⊙ (Kobayashi et al. 2022). At early stages, the shock fronts deform by interaction with the inflow density inhomogeneity. These curved shock fronts introduce inhomogeneity in the postshock entropy distribution. Such entropy inhomogeneity generates turbulent vorticities accounting for a small fraction of the solenoidal mode (see Kida & Orszag 1990). At later stages after 1.0 t_cool, the interaction between the formed CNM clumps and the shock fronts significantly deforms the shock fronts. This deformation creates a number of oblique shocks, whose maximum size is comparable to the size of the shock-compressed layer. This induces strong shear motion into the shock-compressed layer, so that the solenoidal mode dominates the turbulence. Our results show that the solenoidal-mode-dominated turbulence also emerges in the low-metallicity environment.

The results presented in this Section suggest that in the range of 1.0–0.2 Z_⊙, molecular clouds forming in the shock-compressed layer have a turbulent pressure close to the inflow ram pressure. Note that ≳80% of the turbulence power is in the solenoidal mode at all of the metallicities. This indicates that even in a galactic-scale converging region, forming molecular clouds are always solenoidal-mode dominated. Therefore, galactic-scale compressive motion is important in the formation of molecular clouds, but it does not immediately imply enhancement of star formation efficiency by enhancing compressive motion in molecular clouds.

We identify the CNM structures to further investigate their properties. In this section, we define the CNM as the volume with T < 200 K and n > 20 cm⁻³.

First, we perform the Friends-of-Friends algorithm to identify CNM clumps as groups of connected CNM cells. Each clump has more than 64 member cells to avoid numerical noise on small scales. Figure 9 compares the size and mass functions of the identified CNM clumps at 3 t_cool. We define the size l of each CNM clump as $l=\sqrt{{I}_{\max }/M}$ where ${I}_{\max }$ and M are the maximum eigenvalue of its inertia matrix and the mass of each clump. The size distribution peaks at ∼0.1 pc at all of the metallicities. ¹³ The mass functions follow a power-law distribution of ${dn}/{dm}\propto {m}^{-1.7}$ up to ∼10² M_⊙. Some CNM mass tends to stagnate in the central region of the shock-compressed layer due to the converging flow configuration (see Figure 3 and Appendix C). Large CNM clumps obtain more mass or coagulate with other CNM gas so that the most massive clumps deviate from the power-law distribution.

Zoom In Zoom Out Reset image size

An index of −1.7 has often been reported in the 1.0 Z_⊙ environment by previous converging WNM flow simulations in two dimensions (e.g., Heitsch et al. 2005; Hennebelle & Audit 2007; Hennebelle et al. 2008) and in three dimensions (e.g., Inoue & Inutsuka 2012). This power-law function is explained via the statistical growth of the thermal instability with a given density fluctuation spectrum. A functional form of  $dn/dm\propto {m}^{(\alpha -3)/3-2}$ is expected when the seed density fluctuation has a three-dimensional averaged spectral index of α (we refer the readers to Hennebelle & Audit (2007), for a detailed explanation and the derivation of the index (α -3)/3-2)). When α = 11/3, i.e., the Kolmogorov fluctuation, this gives ${dn}/{dm}\propto {m}^{-1.7}$, which is indeed consistent with the numerical results of previous papers in 1.0 Z_⊙ environments. Our results suggest that, at the same t/t_cool, the same CNM mass spectrum can also appear in lower-metallicity environments due to the thermal instability with the Kolmogorov turbulence background.

Second, we investigate the turbulence structure of the CNM volume alone by measuring the two-point correlation of the CNM velocity field. We measure the second-order velocity structure function

$$\begin{eqnarray}&&S(r)=\langle {\left|v({\boldsymbol{r}}+{\boldsymbol{x}})-v({\boldsymbol{x}})\right|}^{2}\rangle ,\end{eqnarray} \tag{ 10 }$$

where $r=\left|{\boldsymbol{r}}\right|$, and x is the three-dimensional position of the CNM cells. We select 27 subvolumes in the shock-compressed layer, where each subvolume is a (3.3 pc)³ cube, and their volume-centered position is (x, y, z) = (0 ± 3.3 pc, 0 ± 3.3 pc, 0 ± 3.3 pc). Figure 10 shows $\sqrt{S(r)}$ from each metallicity at 3 t_cool. The CNM volume overall shows the transition toward the Larson-type scale-dependence as $\sqrt{S(r)}\propto {r}^{1/2}$. This relation extends from small scales within individual CNM clumps to large scales beyond those clumps. This indicates that the commonly observed Larson's law in molecular clouds is inherited from the CNM phase and is also ubiquitous in low-metallicity environments down to 0.2 Z_⊙.

Zoom In Zoom Out Reset image size

Note that the dynamic range is still limited if we consider only the scales without the numerical diffusion (i.e., the scales larger than the vertical dotted line). S(r) in this range is between the Kolmogorov and Larson relations. Further high-resolution simulations are required to finally conclude that the CNM velocity structure function follows the Larson-type relation on all scales. Nevertheless, the strongest turbulence power within the CNM clumps is in eddies whose scale is comparable to the typical CNM clump size ∼0.1 pc. Our results, therefore, suggest that the internal velocity dispersion within the CNM clumps remains ≲1 km s⁻¹ while the clump-to-clump relative velocity is 3–5 km s⁻¹.

Note that the amplitude of $\sqrt{S(r)}$ at 0.2 Z_⊙ is smaller by a factor of 2 than that at 1.0 Z_⊙, as we see in $\sqrt{\langle \delta {v}^{2}{\rangle }_{\rho }}$ (see Figure 6). This is consistent with the comparison of the CO line width between clouds in the LMC/SMC and the Milky Way (Bolatto et al. 2008; Fukui & Kawamura 2010; Ohno et al. 2023), except for those in extreme star-forming systems, such as 30 Dor R136 cluster (Indebetouw et al. 2013). Also note that a variation at a given spatial displacement r is larger in the lower-metallicity environment because the shock-compressed layer is thicker, and thus the total number of cells in this analysis increases toward lower metallicity.

The results presented in this Section show that CNM mass spectrum is ${dn}/{dm}\propto {m}^{-1.7}$, common in the 1.0–0.2 Z_⊙ range with a Kolmogorov turbulence background. The second-order structure function of the CNM volume alone indicates the transition toward the Larson-type turbulence scale-dependence ubiquitously at all metallicities with its amplitude smaller by a factor of 2 at 0.2 Z_⊙ compared with 1.0 Z_⊙.

Combining all of the results in Section 3, Figure 11 schematically summarizes the hierarchical thermal/turbulent structure of the multiphase medium in this metallicity range.

Zoom In Zoom Out Reset image size

On one hand in the 1.0 Z_⊙ context, previous authors often investigated cloud formation by supersonic flows in supernova remnants, galactic spirals, etc. (Kim et al. 2020; Chevance et al. 2022). For example, a single supernova remnant expands to 50 pc size in 0.3 Myr (∼0.3 t_cool) and accumulates 10⁴ M_⊙ H i mass (e.g., Kim & Ostriker 2015). The compilation of multiple supernovae events is able to create even more-massive molecular clouds (Inutsuka et al. 2015). On the other hand, in a low-metallicity environment, the molecular cloud formation out of the WNM alone requires longer physical time. This is the case even in the configuration where the mean magnetic field is completely parallel to the WNM inflow, as we have studied, which forms molecular clouds most efficiently compared with the magnetic fields inclined against the flow (Inoue & Inutsuka 2008; Iwasaki et al. 2019).

To estimate how this timescale impacts molecular cloud formation in low-metallicity environments, let us assume that all of the CNM volume eventually evolves to the molecular gas. Based on our simulation results, ¹⁴ the expected mass of the cloud, M_cloud, is approximately

$$\begin{eqnarray}&&\begin{array}{l}{M}_{\mathrm{cloud}}\simeq 1.1\times {10}^{3}\,{M}_{\odot }\\ \left(\displaystyle \frac{{n}_{0}}{0.57\,{\mathrm{cm}}^{-3}}\right)\left(\displaystyle \frac{{V}_{\mathrm{inflow}}}{20\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right){\left(\displaystyle \frac{L}{10\,\mathrm{pc}}\right)}^{2}\left(\displaystyle \frac{t}{3{t}_{\mathrm{cool}}}\right)\end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\propto {\left(\displaystyle \frac{{n}_{0}}{0.57\,{\mathrm{cm}}^{-3}}\right)}^{2}\left(\displaystyle \frac{{V}_{\mathrm{inflow}}}{20\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right){\left(\displaystyle \frac{L}{10\,\mathrm{pc}}\right)}^{2}\left(\displaystyle \frac{Z}{0.2{Z}_{\odot }}\right)\left(\displaystyle \frac{t}{15\,\mathrm{Myr}}\right).\end{eqnarray} \tag{ 12 }$$

Here L² is the inflow cross section in our calculation domain, and thus the first three factors in Equation (11) come from the inflow mass flux. The exact dependence of t_cool on n₀ and V_inflow is still left for future studies, which is involved in the conversion from Equation (11) to Equation (12). We here employ ${t}_{\mathrm{cool}}\propto {n}_{0}^{-1}{Z}^{-1}$ as the fiducial dependence, which is consistent with our definition of t_cool as calculated in Table 1 (see Appendix B.2 and Equation (B7)).

If we suppose the case of a single flow event with n₀ = 0.57 cm⁻³ and V_inflow = 20 km s⁻¹ creates a typical maximum cloud of a few 10⁵ M_⊙ in the LMC/SMC, Equation (12) indicates that such a WNM flow should continue coherently on scales of L ≃ 100 pc and 15 Myr. Note that 15 Myr is 1 order of magnitude longer compared with the typical expansion timescale of a single supernova remnant. A superbubble rather than a single supernova event is more likely to keep such a coherent one-directional flow over a 15 Myr timescale. However, the prevalent existence of molecular clouds that are not associated with superbubbles in the LMC/SMC indicates that faster flows and/or the preexistence of CNM structure in the WNM are important for cloud formation and evolution in a lower-metallicity environment (see Dawson et al. 2013). ¹⁵

One possibility is a fast H i gas flow, which is observed as the tidal interaction between the LMC and the SMC (Fukui et al. 2017; Tsuge et al. 2019). Its velocity is as high as 50–100 km s⁻¹, comparable to the escape velocity of the LMC/SMC. We can straightforwardly expect the coherent continuation of such a flow over a few 10 Myr timescale because it travels on a 10 kpc scale between the LMC and the SMC. However, such a fast flow induces strong turbulence in molecular clouds, which can also destroy dense structures (see Klein et al. 1994; Heitsch et al. 2006a). In addition, even for such efficient mass accumulation, thermal evolution from the WNM to molecular clouds still requires cooling. For example, albeit on an L = 100 pc scale with a coarser spatial resolution of 0.2 pc, Maeda et al. (2021) performed convergence of WNM flow simulations with V_inflow = 100 km s⁻¹ but without any CNM structure preexisting in the WNM flow. In the 0.2 Z_⊙ environment with an initial WNM density of n₀ = 0.75 cm⁻³, they showed that the molecular cloud formation takes 23 Myr. The cloud mass with n > 10⁴ cm⁻³ is a few 10⁵ M_⊙ at 23 Myr, which is consistent with Equation (12). ¹⁶

Therefore, even in the fast flow environment, the preexistence of CNM structure in the WNM is important (see Pringle et al. 2001; Inoue & Inutsuka 2009). This leads to another question of, "how does the ISM form the CNM in low-metallicity environments?" As can be seen in our simulations, low-mass CNM clumps form if the mean magnetic field is parallel to compressive events, and this inflow can be as slow as V_inflow = 20 km s⁻¹. Therefore, some previous generations of flows with a few tens of kilometers per second along with the mean magnetic field (such as a single supernova) are responsible for introducing CNM clumps in the WNM. Such CNM preexistence and subsequent compression due to shocks presumably determine the formation site, the mass, and the structure of molecular clouds. Indeed, based on the recent ALMA observations toward N159E/W regions in the LMC, Tokuda et al. (2022) indicated that a fast H i gas flow with a multiple-scale substructure potentially explains the formation of filamentary molecular clouds on various size scales. Such a hierarchical structure may also determine the fractal nature of the young stellar structures in the LMC, suggested by the spatial clustering analysis of star clusters (e.g., Miller et al. 2022). In addition, the similarity between the power-law spectrum of the molecular cloud mass function (e.g., ${dn}/{dm}\propto {m}^{-\alpha }$ where α = 1.7–1.9 in the LMC; Fukui et al. 2001; Ohno et al. 2023) and that of CNM clump mass function in our simulation also indicates that molecular cloud formation is predetermined by CNM structure (Figure 9). This implication requires further study to be confirmed in the future.

The H₂ is an important coolant in the context of the primordial gas cooling. Chemical paths to H₂ formation exist even in the primordial gas, especially that via the gas-phase interaction between H and H⁻, where the electron fraction impacts the abundance of H⁻. Susa et al. (1998) calculated the electron fraction in the postshock layer of supersonic ISM flows, and Susa & Umemura (2004) showed that, with V_inflow = 30 km s⁻¹, H₂ cooling is dominant over the metal line cooling in the postshock region even in 0.1 Z_⊙ without any UV background that dissociates the H₂. This is not the case in our calculation where we employ G₀ = 1 for all of the metallicity by assuming ongoing active star formation as in the LMC/SMC. The metal line cooling is always dominant in this setup.

The impact of H₂ cooling on thermal instability is limited typically to G₀ ≤ 10⁻³ environments (Inoue & Omukai 2015; see also Glover & Clark 2014). It is left for future studies to investigate the dynamical formation and evolution of the multiphase ISM driven by supersonic flows in lower-metallicity environments with a UV radiation field of G₀ < 1.

In this section, we introduce several previous studies to list how the adopted assumptions impact our results.

First, many previous studies also assume that the dust-to-gas ratio is linearly scaled with the gas-phase metallicity in their fiducial models. Meanwhile, based on the extragalactic observational indications (e.g., Rémy-Ruyer et al. 2014), there are also studies testing a broken power-law dependence. For example, ${ \mathcal D }\propto {Z}^{2}$ in Glover & Clark (2014) and Z³ in Bialy & Sternberg (2019), where ${ \mathcal D }$ is the dust-to-gas ratio. Such a superlinear decrease of the dust abundance at ≲0.2 Z_⊙ reduces the photoelectric heating rate relatively to the C ii cooling rate. The thermally stable states of the UNM and the CNM are cooler in this superlinear case ¹⁷ . However, in the range of ≥0.2 Z_⊙, they show that this difference has a limited impact on the CNM thermal equilibrium state (see also Figure 5 of Kim et al. 2023b). In addition, this superlinear relation of the dust-to-gas ratio does not significantly impact the overall dynamics in our simulations because C ii cooling at the shock-heated WNM right after the shock compression provides the typical evolutionary timescale (as discussed in Appendix B.1).

Second, the gas-phase electron fraction changes the photoelectric heating efficiency. The photoelectric heating rate by FUV radiation can be estimated as Γ_phot(1.0 Z_⊙)= 1.3 × 10⁻²⁴ G₀ erg s⁻¹ where is the photoelectric heating efficiency as

$$\begin{eqnarray}\begin{array}{rcl}\epsilon & = & \displaystyle \frac{4.9\times {10}^{-2}}{1+4.0\times {10}^{-3}{({G}_{0}{T}^{1/2}/{n}_{e}{\phi }_{\mathrm{PAH}})}^{0.73}}\\ & & +\,\displaystyle \frac{3.7\times {10}^{-2}{(T/{10}^{4})}^{0.7}}{1+2.0\times {10}^{-4}({G}_{0}{T}^{1/2}/{n}_{e}{\phi }_{\mathrm{PAH}})}\,\end{array}\end{eqnarray} \tag{ 13 }$$

(see Equation (43) in Bakes & Tielens 1994). Here n_e is the electron number density so that the electron fraction can be defined as x_e = n_e /n_H, and the polycyclic aromatic hydrocarbon (PAH) reaction rate is ϕ_PAH = 0.5 (Wolfire et al. 2003). The G₀ T^1/2/n_e dependence in the denominator describes the relative importance of the ionization against the recombination on the dust surface. ¹⁸ Our fiducial value, Γ_phot(1.0 Z_⊙) = 2.0 × 10⁻²⁶ erg s⁻¹ in Equation (A2), corresponds to the typical condition of the thermally stable WNM at 1.0 cm⁻³ with x_e ∼ 0.02 determined by cosmic-ray ionization (see Equation (12) of Bialy & Sternberg 2019). The photoelectric heating efficiency may deviate from this fiducial value right after the shock compression where the collisional ionization increases the electron fraction. Koyama & Inutsuka (2000) performed a high-resolution simulation (albeit one-dimensional) of shock propagation into the WNM. The postshock shock-heated WNM volume reaches n ∼ 5 cm⁻³, x_e ∼ 0.1, and T ∼ 6400 K, so that Equation (13) predicts Γ_phot(1.0 Z_⊙) = 6 × 10⁻²⁶ erg s⁻¹. The factor of 3 difference from our fiducial value can expand the parameter space for the thermally stable WNM. However, note that a factor of a few uncertainty also exists in the assumption on the PAH size distribution and shape (Kim et al. 2023b). Also note that such high x_e volume is limited to a 10⁻¹ pc scale from the shock front by the recombination (i.e., Figure 3 of Koyama & Inutsuka 2000). Therefore, we opt to use the fiducial constant photoelectric heating efficiency (and its linear dependence on the metallicity) for our comparison between metallicities in this article.

Lastly, the interstellar radiation field varies across space. It is ideal to numerically investigate the 10 pc-scale local radiation field, e.g., zoomed-in approach consistently coupled with galactic-scale star formation (e.g., Peters et al. 2017; Kim et al. 2023a, 2023b), and this is left for future studies.

Our results indicate that the solenoidal-mode-dominated turbulence is generally realized in molecular clouds undergoing supersonic compressions in the metallicity range 1.0–0.2 Z_⊙. The density inhomogeneity (and/or velocity fluctuation) in the converging flows can be larger (can exist) in reality, especially if the flows already have multiple phases over the one-phase WNM (Inoue & Inutsuka 2012; Iwasaki et al. 2019). Such a strong density/velocity inhomogeneity leads to stronger shock deformation and induces stronger shear motion into the shock-compressed layer. The stronger inhomogeneity in the inflow density/velocity induces stronger turbulence (e.g., the systematic survey by Kobayashi et al. 2020), but the turbulent speed is expected to be on the order of the WNM sound speed (see Appendix C).

Expansion of multiple supernova remnants also induces cloud compressions as often observed in galactic-scale numerical simulations (Girichidis et al. 2016; Kim et al. 2023a) and may explain the bubble features in nearby galaxies observed in JWST (Watkins et al. 2023). Each compression event occurs at diverse angles at various times, and the generation of solenoidal mode turbulence locally occurs due to the deformed shock fronts. It is likely that the compressive mode of turbulence does not become dominant unless multiple compressions from various angles occur simultaneously with relatively uniform density or unless the cloud itself becomes massive enough to gravitationally collapse. Previous analytic studies also show that the solenoidal mode of turbulence can grow faster than the compressive mode even in a gravitationally contracting background (e.g., Higashi et al. 2021).

Our simulations show that physical properties of molecular clouds resemble one another if compared at the same time in units of the cooling time. What about those in lower-metallicity ranges, which are important in understanding the overall cosmic star formation history beyond the Magellanic Clouds? The metallicity dependence that we have shown comes primarily from the line coolings of [O i] (63.2 μm) and [C ii] (157.7 μm). This metallicity dependence in cooling rate is expected to continue down to Z ∼ 10⁻³ Z_⊙, until which [O i] and [C ii] are the dominant coolants to induce thermal instability. In contrast, the main heating process changes from the photoelectric heating to the X-ray and cosmic-ray heating at ∼0.1 Z_⊙. The heating rate in <0.1 Z_⊙ becomes independent of the metallicity if we employ the same constant X-ray and cosmic-ray heating rate. Under such conditions in <0.1 Z_⊙, Bialy & Sternberg (2019) showed that the density of the thermally stable CNM scales roughly as n_CNM ∝ Z⁻¹ (see their Equation (18)). This results in a roughly constant cooling time in <0.1 Z_⊙. In an even lower metallicity range of Z ≲ 10⁻³ Z_⊙, the UV field strength impacts the growth of thermal instability (due to H₂ dissociation), and thus the comprehensive investigations remain by also changing the UV field strength for future studies. Overall, we assume that, in the range of ≥0.1 Z_⊙, physical properties of molecular clouds are similar between metallicities at the same time measured in units of the cooling time; we require further investigation for the range of <0.1 Z_⊙.

Meanwhile, our results also indicate that a difference in cloud formation conditions (e.g., different n₀ and V_inflow) is important in the formation of molecular clouds with different properties, which is required to comprehensively understand cosmic star formation history. For example, recent ALMA observations started to spatially resolve the galactic disk of star-forming galaxies at Cosmic Noon at the redshifts of 2–4. They showed the existence of molecular clouds in gravitationally unstable galactic disks and that clouds tend to have higher column density/higher star formation rate density (e.g., Σ_gas > 100 M_⊙ pc⁻² and Σ_SFR > 1 M_⊙ yr⁻¹ kpc⁻²) than those in the Milky Way (e.g., Tadaki et al. 2018). Such properties resemble those observed in nearby luminous infrared galaxies (Kennicutt & De Los Reyes 2021), who possibly experience galaxy mergers. These indicate a possibility that some drastic mechanisms, such as galaxy mergers with high velocity, form high column density molecular clouds that host active star formation at Cosmic Noon.

In addition, even starting with the same cloud formation condition to form a similar molecular cloud across the metallicity, the variation of the stellar initial mass function with the metallicity, if any, may arise from subsequent fragmentation in the collapse of clouds and circumstellar disks. We remark that some previous authors investigated this phenomena with numerical simulations using the same initial cloud properties for all metallicities (e.g., Bate 2019; Chon et al. 2021).

For the time being, the above discussions are merely speculations. Further simulations with various initial density, inflow velocity, etc., are needed and left for future studies.

The cooling time, $P/(\gamma -1)/\rho { \mathcal L }$, varies locally with space and time, depending on the local density/temperature variation. Nevertheless, the typical cooling time of this converging flow system can be estimated based on the typical state of the shock-heated WNM as follows.

To determine the representative thermal state, let us start with a simple estimation in a perfect one-dimensional shock compression. Once the inflow WNM passes the shock front, the density in the downstream increases to n ≃ 4n₀ = 2.3 cm⁻³. Meanwhile, the temperature also increases to T ≃ P_ram/k_B/(4n₀) = 1.5 × 10⁴ K, but it rapidly decreases to T ≃ T₀ because of the efficient coolings by Lyα and heavy metals' lines. For example, the temperature decreases on the timescale of 2.0 × 10⁻⁴ Myr, 4.1 × 10⁻⁴ Myr, and 1.0 × 10⁻³ Myr at 1.0, 0.5, and 0.2 Z_⊙ based on the T > 14,577 K regime of Equation (A3) with n = 4n₀ and T = P_ram/k_B/(4n₀). These timescales are shorter than the typical sound crossing time on a single numerical cell of 0.02 pc size with T = P_ram/k_B/(4n₀), which is 1.2 × 10⁻³ Myr. Therefore, we need to perform simulations with a higher spatial resolution to confirm whether this shock-heated WNM evolves isochorically or isobarically during the return to T≃T₀. In the following, let us assume that this rapid cooling leads to rather isochorical evolution and employ n = 4n₀ and T = T₀ as the representative initial density and temperature to estimate the typical t_cool (see also the discussions in Appendix B.2). The corresponding pressure, 4n₀ k_B T₀, is smaller than the ram pressure by roughly a factor of 2 (P/k_B ≃ P_ram/2k_B = 1.75 × 10⁴ erg cm⁻³). This is consistent with most of the volumes in our simulations as shown in Figure 4. Such a thermal pressure smaller than the ram pressure is also achieved because the turbulent pressure almost balances with the inflow ram pressure, as shown in Figure 6.

Starting with n = 4n₀ and T = T₀, the net cooling rate during the subsequent thermal evolution is mostly determined by the C_II cooling given by the first equation of Equation (A3):

$$\begin{eqnarray}\begin{array}{rcl}\rho { \mathcal L } & \simeq & {n}^{2}{\rm{\Lambda }}(T,Z)\\ & = & {n}^{2}\times 2.8\times {10}^{-28}(Z/{Z}_{\odot })\sqrt{T}\exp (-92/T)\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-3}.\end{array}\end{eqnarray} \tag{ B1 }$$

Therefore, the cooling time can be estimated in general as

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}=\displaystyle \frac{P/(\gamma -1)}{{n}^{2}{\rm{\Lambda }}(T,Z)}\end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray}&&\begin{array}{l}\simeq \,2.99\,\mathrm{Myr}\left(\displaystyle \frac{P/{k}_{{\rm{B}}}}{{10}^{4}\,{\rm{K}}\,{\mathrm{cm}}^{-3}}\right){\left(\displaystyle \frac{n}{1\,{\mathrm{cm}}^{-3}}\right)}^{-2}\\ \ \times \,{\left(\displaystyle \frac{Z}{{Z}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{T}{6400\,{\rm{K}}}\right)}^{-1/2}\exp \left(1-\displaystyle \frac{T}{6400\,{\rm{K}}}\right)\end{array}\end{eqnarray} \tag{ B3 }$$

with n = 4n₀, T = T₀, and P/k_B = P_ram/2k_B; the typical cooling time is

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{cool}} & = & 2.99\,\mathrm{Myr}\left(\displaystyle \frac{{P}_{\mathrm{ram}}/2{k}_{{\rm{B}}}}{{10}^{4}\,{\rm{K}}\,{\mathrm{cm}}^{-3}}\right){\left(\displaystyle \frac{4{n}_{0}}{1\,{\mathrm{cm}}^{-3}}\right)}^{-2}\\ & & \times \,{\left(\displaystyle \frac{Z}{{Z}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{{T}_{0}}{6400\,{\rm{K}}}\right)}^{-1/2}\exp \left(1-\displaystyle \frac{{T}_{0}}{6400\,{\rm{K}}}\right)\end{array}\end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray}&&\simeq 1.01\,\mathrm{Myr}{\left(\displaystyle \frac{Z}{{Z}_{\odot }}\right)}^{-1}.\end{eqnarray} \tag{ B5 }$$

Equation (B4) gives the typical value of t_cool as 1.01 Myr at 1.0 Z_⊙, 2.03 Myr at 0.5 Z_⊙, and 5.11 Myr at 0.2 Z_⊙, as summarized in Table 1. Equation (B5) gives a rough metallicity dependence.

In the current article, we focus on the metallicity dependence starting with fixed initial conditions. Therefore, it is left for future studies to investigate the exact dependence of t_cool on the initial conditions other than the metallicity. This requires parameter surveys by changing the mean density, the inflow velocity, the magnetic field strength, the inclination of the mean magnetic fields, the UV background, etc. The resultant impacts by all of these changes are coupled to each other because they modify the ram pressure, the turbulent pressure, and the following resultant shock-heated WNM state.

Nevertheless, as a first step, let us list several possible dependences of the two parameters, n₀ and V_inflow. For simplicity, we here assume that the metallicity dependence comes only from the cooling/heating rate and is always independent from the choice of n₀ and V_inflow. The other conditions are the same as our simulation setup in this article (e.g., the WNM inflow is parallel to the orientation of the mean magnetic field with B = 1 μG, the injected WNM is on the thermally stable state, etc.). These simplistic assumptions have to be investigated by further simulations, which is left for future studies.

Suppose that the postshock pressure, density, and temperature have dependences of $P\propto {n}_{0}^{\alpha }{V}_{\mathrm{inflow}}^{\beta }$, $n\propto {n}_{0}^{\gamma }{V}_{\mathrm{inflow}}^{\delta }$, and $T=P/{{nk}}_{{\rm{B}}}\propto {n}_{0}^{\alpha /\gamma }{V}_{\mathrm{inflow}}^{\beta /\delta }$, respectively, then

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}\propto {n}_{0}^{\tfrac{\alpha -3\beta }{2}}\,{V}_{\mathrm{inflow}}^{\tfrac{\beta -3\delta }{2}}\,{Z}^{-1}.\end{eqnarray} \tag{ B6 }$$

If we employ n = 4n₀ and T = T₀ as the shock-heated WNM state right behind the shock, as we did in Appendix B.1, α = γ = 1 and β = δ = 0. The cooling time is

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}\propto {n}_{0}^{-1}\,{Z}^{-1}.\end{eqnarray} \tag{ B7 }$$

We employ Equation (B7) as our fiducial dependence to derive Equation (12).

Alternatively, most of the shock-heated WNM volume evolves close to T = T₀ even beyond n > 4n₀ (Figure 4), and thus we may consider that an isothermal shock approximation describes the overall postshock state. In this case, α = γ = 1 and β = δ = 2. The cooling time is

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}\propto {n}_{0}^{-1}{V}_{\mathrm{inflow}}^{-2}\,{Z}^{-1},\end{eqnarray} \tag{ B8 }$$

accordingly.

We may also consider a simpler condition with P = P_ram and T = T₀, and n = P_ram/k_B/T. This gives α = β = δ = 2 and γ = 1, and the cooling time is

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}\propto {n}_{0}^{-1/2}{V}_{\mathrm{inflow}}^{-2}\,{Z}^{-1},\end{eqnarray} \tag{ B9 }$$

However as shown in our simulations, such a perfect isothermal compression is limited to the central volume during the very early stage when the shock compression is still close to a plane-parallel geometry without significant deformation (see panel (c) of Figure 3 and panel (c) of Figure 4).

Lastly, we would like to note that the above arguments assume the maximum pressure of the thermally stable state is independent of the metallicity. The flow ram pressure in the <0.1 Z_⊙ environments must be higher (e.g., with a higher n₀, a faster V_inflow) to successfully overcome this maximum pressure to enter the thermally unstable regime, because the maximum pressure depends on the metallicity roughly as ∝ 1/Z (see Figure 6 of Bialy & Sternberg 2019).