On the resolution requirements for modelling molecular gas formation in solar neighbourhood conditions

TL;DR

This study determines the spatial resolution needed in simulations to accurately model molecular hydrogen and CO formation in the solar neighborhood, emphasizing the importance of resolving specific density and timescale conditions.

Contribution

It provides critical resolution criteria and formulas to ensure convergence of molecular gas formation modeling in hydrodynamical simulations.

Findings

H$_2$ converges at resolutions below 0.2 pc.

CO requires resolutions below 0.04 pc for convergence.

Derived conditions help verify molecule formation convergence in simulations.

Abstract

The formation of molecular hydrogen (H$_2$) and carbon monoxide (CO) is sensitive to the volume and column density distribution of the turbulent interstellar medium. In this paper, we study H$_2$ and CO formation in a large set of hydrodynamical simulations of periodic boxes with driven supersonic turbulence, as well as in colliding flows with the \textsc{Flash} code. The simulations include a non-equilibrium chemistry network, gas self-gravity, and diffuse radiative transfer. We investigate the spatial resolution required to obtain a converged H$_2$ and CO mass fraction and formation history. From the numerical tests we find that H$_2$ converges at a spatial resolution of $\lesssim0.2$~pc, while the required resolution for CO convergence is $\lesssim 0.04$~pc in gas with solar metallicity which is subject to a solar neighbourhood interstellar radiation field. We derive two critical conditions from our numerical results: the simulation has to at least resolve the densities at which (1) the molecule formation time in each cell in the computational domain is equal to the dissociation time, and (2) the formation time is equal to the the typical cell crossing time. For both H$_2$ and CO, the second criterion is more restrictive. The formulae we derive can be used to check whether molecule formation is converged in any given simulation.

Tables

Table 1. Table 1: Initial conditions for the different TB and CF simulations. Initial gas density ρ 0 subscript 𝜌 0 \rho_{0} (2 nd column), number density n 0 subscript 𝑛 0 n_{0} of H nuclei (3 rd columns), temperature T 0 subscript 𝑇 0 T_{0} (4 th column), and the sound speed c s subscript 𝑐 𝑠 c_{s} (5 th column). The initial relative abundances of x He = 0.1 subscript 𝑥 He 0.1 x_{\text{He}}=0.1 , x C = 1.4 × 10 − 4 subscript 𝑥 C 1.4 superscript 10 4 x_{\text{C}}=1.4\times 10^{-4} , x O = 3.2 × 10 − 4 subscript 𝑥 O 3.2 superscript 10 4 x_{\text{O}}=3.2\times 10^{-4} , and x M = 1 × 10 − 7 subscript 𝑥 M 1 superscript 10 7 x_{\text{M}}=1\times 10^{-7} are common for all the TB and CF simulations. The initially atomic gas with μ = 1.27 𝜇 1.27 \mu=1.27 contains all of the hydrogen as H atoms and all of the carbon as C + ions.

Setup	${\rho }_0$ [g cm-3 ]	$n_0$ [cm-3 ]	$T_0$ [K]	$c_s$ [km s-1 ]
TB	$7\times {10}^{-24}$	3	880	2.4
TB	$7\times {10}^{-23}$	30	122	0.9
TB	$7\times {10}^{-22}$	300	88	0.7
CF	$1.7\times {10}^{-24}$	0.75	4082	5.1

Table 2. Table 2: List of simulations performed for this study. TB indicates the Turbulent Box setup and CF indicates the Colliding Flows setup. R[4-10] denotes the refinement level, where R4 corresponds to a cell size of Δ ​ x = 0.500 Δ 𝑥 0.500 \Delta x=0.500  pc and R10 corresponds to an effective resolution of Δ ​ x = 0.008 Δ 𝑥 0.008 \Delta x=0.008  pc. There are two sets of TB runs, one in a (32 pc) 3 box with v rms = 10 subscript 𝑣 rms 10 v_{\rm rms}=10  km s -1  and one in a (8 pc) 3 box with different v rms subscript 𝑣 rms v_{\rm rms} . The latter has a “-v[x]" tag and is used in Section  5.4 to verify our resolution requirement. In the TB runs, the initial density of the H nuclei in the homogeneous medium is denoted by the ending tag (e.g. n003 corresponds to n 0 = 3 subscript 𝑛 0 3 n_{0}=3  cm -3  , see column 3). The (effective) resolution of the simulations is denoted by Δ ​ x Δ 𝑥 \Delta x . In the TB runs, the value in parentheses denotes the number of finite volume cells used. In the CF runs, the “AMR" in parentheses denotes that the effective resolution is obtained via adaptive mesh refinement.

Run	$\Delta x$	$n_0$	$v_{\text{rms}}$
		[cm-3 ]	[km s-1 ]
Turbulent box: box size (32 pc)3
TB-R4-n003	0.500 pc (643)	3	10
TB-R4-n030	0.500 pc	30	10
TB-R4-n300	0.500 pc	300	10
TB-R5-n003	0.250 pc (1283)	3	10
TB-R5-n030	0.250 pc	30	10
TB-R5-n300	0.250 pc	300	10
TB-R6-n003	0.125 pc (2563)	3	10
TB-R6-n030	0.125 pc	30	10
TB-R6-n300	0.125 pc	300	10
TB-R7-n003	0.063 pc (5123)	3	10
TB-R7-n030	0.063 pc	30	10
TB-R7-n300	0.063 pc	300	10
Small turbulent box: box size (8 pc)3
TB-R5-n030-v1	0.250 pc (323)	30	1
TB-R5-n030-v2	0.250 pc	30	2
TB-R5-n030-v3	0.250 pc	30	3
TB-R5-n030-v6	0.250 pc	30	6
TB-R6-n030-v1	0.125 pc (643)	30	1
TB-R6-n030-v2	0.125 pc	30	2
TB-R6-n030-v3	0.125 pc	30	3
TB-R6-n030-v6	0.125 pc	30	6
TB-R7-n030-v1	0.063 pc (1283)	30	1
TB-R7-n030-v2	0.063 pc	30	2
TB-R7-n030-v3	0.063 pc	30	3
TB-R7-n030-v6	0.063 pc	30	6
Colliding flow
CF-R4	0.500 pc	0.75	–
CF-R5	0.250 pc	0.75	–
CF-R6	0.125 pc (AMR)	0.75	–
CF-R7	0.063 pc (AMR)	0.75	–
CF-R8	0.032 pc (AMR)	0.75	–
CF-R9	0.016 pc (AMR)	0.75	–
CF-R10	0.008 pc (AMR)	0.75	–

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0,$$

$$\frac{\partial (\rho \mathbf{v})}{\partial t}+\nabla \cdot \left[\rho \mathbf{v}{\mathbf{v}}^{\text{T}}+\left(P+\frac{{\mathbf{B}}^2}{8\pi }\mathbf{I}\right)-\frac{\mathbf{B}{\mathbf{B}}^{\text{T}}}{4\pi }\right]=\rho (\mathbf{g}+\mathbf{F}),$$

$$\frac{\partial E}{\partial t}+\nabla \cdot \left[\left(E+\frac{{\mathbf{B}}^2}{8\pi }+\frac{P}{\rho }\right)\mathbf{v}-\frac{(\mathbf{B}\cdot \mathbf{v})\mathbf{B}}{4\pi }\right]=\rho \mathbf{v}\cdot (\mathbf{g}+\mathbf{F})+{\dot{u}}_{\text{chem}},$$

$$\frac{\partial \mathbf{B}}{\partial t}-\nabla \times (\mathbf{v}\times \mathbf{B})=0,$$

$$E=u+\frac{\rho {\mathbf{v}}^2}{2}+\frac{{\mathbf{B}}^2}{8\pi },$$

$$\Delta \Phi =4\pi G\rho ,$$

$$\text{d}\tilde{\mathbf{F}}(\mathbf{k},t)=\widetilde{F_0}(\mathbf{k},t){\mathcal{P}}^{\zeta }(\mathbf{k})\frac{\text{d}\mathbf{W}(t)}{T_{\text{ac}}}-\tilde{\mathbf{F}}(\mathbf{k},t)\frac{\text{d}t}{T_{\text{ac}}},$$

$$c_{\text{s}}=\sqrt{\frac{k_{\text{B}}T}{\mu m_{\text{p}}}},$$

$$x^{\prime }=A\left[\cos (2-yz)\cos (k_{y}y+\varphi )+\cos (0.5-yz)\sin (k_{z}z)\right]$$

$${\lambda }_{\text{J}}=\sqrt{\frac{\pi c_{\text{s}}^2}{\rho G}},$$

$$\Delta x=\frac{32\text{ pc}}{2^{R-1}\times 8}=\frac{32\text{ pc}}{2^{\text{R}+2}}=2^{3-\text{R}}\text{ pc}.$$

$${\alpha }_{\text{C}}=\frac{1}{x_{\text{C,tot}}\times {\mu }_{\text{C}}}=595,$$

$$f_{{\text{H}}_2}=\frac{M_{{\text{H}}_2}}{M_{\text{H,tot}}}\equiv \frac{2N_{{\text{H}}_2}}{N_{\text{H,tot}}},$$

$${\tau }_{{\text{H}}_2,\text{form}}=\frac{M_{\text{H,tot}}}{{\dot{M}}_{{\text{H}}_2}},$$

$${\eta }_{{\text{H}}_2}=\frac{f_{{\text{H}}_2}}{f_{{\text{H}}_2,{\text{R}}_{\text{max}}}}.$$

$${\tau }_{{\text{H}}_2,\text{eq}}\sim 10^3\text{ Myr}{\left(\frac{\bar{n}}{1{\text{ cm}}^{-3}}\right)}^{-1},$$

$$f_{\text{CO}}=\frac{\frac{12}{28}{M}_{\text{CO}}}{M_{\text{C,tot}}}\equiv \frac{N_{\text{CO}}}{N_{\text{C,tot}}},$$

$${\tau }_{\text{CO},\text{form}}=\frac{M_{\text{C,tot}}}{\frac{12}{28}{\dot{M}}_{\text{CO}}},$$

$${\eta }_{\text{CO}}=\frac{f_{\text{CO}}}{f_{\text{CO},{\text{R}}_{\text{max}}}}.$$

$$M_i=\sum _{j=1}^{{\text{N}}_{\text{cells},i}}{\rho }_j{V}_j,$$

$$\log \rho \in \left[\log {\rho }_i,\log {\rho }_i+d(\log {\rho }_i)\right[.$$

$$P_{\rho ,i}=\frac{M_i}{M_{\text{tot}}},$$

$${\tau }_{\text{mol},\text{form}}\le {\tau }_{\text{mol},\text{dissoc}},$$

$${\tau }_{\text{mol},\text{form}}\le t_{\text{cell,cross}}.$$

$$\text{H}+\text{H(s)}⟶{\text{H}}_2,$$

$${\text{H}}_2+\gamma ⟶\text{H}+\text{H},$$

$$R_{{\text{pd,H}}_2}=3.3\times 10^{-11}{G}_0{f}_{{\text{dust,H}}_2}{f}_{{\text{shield,H}}_2}{\text{ s}}^{-1},$$

$$n_{{\text{H}}_2,\text{physical}}\gtrsim 10\text{ cm−3}.$$

$${\tau }_{{\text{H}}_2,\text{form}}=10^{3.1}\text{ Myr}{\left(\frac{n}{\text{1 cm−3}}\right)}^{-1},$$

$${\text{C}}^{+}+{\text{H}}_2⟶{\text{CH}}_2^{+}+\gamma ,$$

$$\text{C}+{\text{H}}_3^{+}⟶{\text{CH}}_2^{+}+\text{H},$$

$${\text{H}}_2+\text{cosmic-ray}⟶{\text{H}}_2^{+}+{\text{e}}^{-},$$

$${\text{H}}_3^{+}+{\text{e}}^{-}⟶{\text{H}}_2+\text{H},$$

$$R_{{\text{CH}}_2^{+}}=k_{{\text{CH}}_2^{+},1}{n}_{{\text{C}}^{+}}{n}_{{\text{H}}_2}+2{\zeta }_{\text{H}}{n}_{{\text{H}}_2}\frac{k_{{\text{CH}}_2^{+},2}}{k_{\text{dr}}}\frac{n_{\text{C}}}{n_{{\text{e}}^{-}}}$$

Full text

††pubyear: 2018††pagerange: On the resolution requirements for modelling molecular gas formation in solar neighbourhood conditions–B

1 Introduction

Molecular clouds (MCs) are turbulent structures in the interstellar medium (Scalo & Elmegreen, 2004; Mac Low & Klessen, 2004), rich in chemical species that are not in chemical equilibrium (Leung et al., 1984). The overall density structure of the MC is known to significantly affect the chemical composition since it influences the shielding experienced by various regions within the cloud (e.g. Stutzki & Guesten, 1990; Safranek-Shrader et al., 2017). Glover et al. (2010) suggest that the MC composition is dependent on the history of the gas since the dynamical time-scale and the molecule formation time-scale are comparable in MC conditions. Thus, simulations of MCs should couple the dynamical and chemical evolution to obtain correct approximations of the MC formation in nature. In doing so, the chemical abundances have to be evaluated at every hydrodynamical time-step; thus, only simple chemical networks are affordable. Koyama & Inutsuka (2000, 2002) and Bergin et al. (2004) implemented the non-equilibrium chemistry of H2 and CO in their one- and two-dimensional hydrodynamic simulations because CO is the most important tracer of H2 and an essential coolant in the local interstellar medium (ISM). With the increase in computational capabilities, chemistry networks with more species and reactions were used in three-dimensional magneto-hydrodynamic (MHD) simulations (e.g. Glover & Mac Low, 2007a, b). Later, Glover et al. (2010) presented a chemistry network with 32 species linked by 218 reactions, which allows for an accurate treatment of H2 and CO formation. A simpler, computationally cheaper network, based on Glover & Mac Low (2007a, b) and Nelson & Langer (1997) has been applied in many recent high-resolution (M)HD simulations (e.g. Walch et al., 2011; Micic et al., 2012; Walch et al., 2015; Girichidis et al., 2016; Gatto et al., 2017; Pardi et al., 2017; Peters et al., 2017; Seifried et al., 2017). Similarly, a few works (e.g. Glover & Clark, 2012; Clark et al., 2012b; Szűcs et al., 2014) have used a more complex network (termed “NL99") that combines the CO chemistry of Nelson & Langer (1999) with the hydrogen chemistry of Glover & Mac Low (2007a, b). The NL99 network is used in the work presented here. Alternative, simplified networks for H2 and/or CO formation have also been presented by Keto & Caselli (2008), Valdivia et al. (2016), and Gong et al. (2017). In simulations of star-forming filaments, Seifried & Walch (2016) presented the use of a heavier network containing 287 reactions between 37 chemical species, implemented as part of the Krome package (Grassi et al., 2014). Recently, Mackey et al. (2018) presented simulations with an updated NL99 chemistry network that incorporates updated reaction rates suggested by Gong et al. (2017).

Many of the simulations mentioned above have been performed using one of the two standard models of MC formation, namely the turbulent periodic box and the colliding flow simulations. The Turbulent box (TB) simulations became popular following the understanding that the energy dissipation in supersonic turbulence observed in MC conditions is comparable to the typical free-fall time (Mac Low et al., 1998; Stone et al., 1998), while the MCs themselves are a few free-fall times old; thus, energy has to be continuously injected into the MC to maintain such turbulence. In the last two decades, the nature of MCs as dynamical objects dominated by turbulence has been explored via simulations of driven supersonic turbulence (e.g. Padoan et al., 1999; Mac Low & Klessen, 2004; Kritsuk et al., 2011; Federrath, 2013). Recently, Padoan et al. (2016) suggested that such TB models are close approximations to the supernova-driven turbulence in the ISM. By now, multiple TB simulations using coupled chemistry have been performed (e.g. Glover & Mac Low, 2007a; Glover et al., 2010; Walch et al., 2011; Micic et al., 2012; Glover & Clark, 2012; Padoan et al., 2017).

The TB models describe random flows, and to form a dense cloud, some sort of convergent flow is required. The Colliding flow (CF) simulation is a simple scenario to model such gas accumulation, where two flows of warm neutral medium (WNM) converge to quickly form clumps of cold molecular gas in the resulting dense and turbulent collision interface. An extensive set of CF simulations have been performed (e.g. Ballesteros-Paredes et al., 1999; Heitsch et al., 2006; Vázquez-Semadeni et al., 2007; Hennebelle & Audit, 2007; Hennebelle et al., 2007; Hennebelle et al., 2008; Banerjee et al., 2009; Clark et al., 2012b; Körtgen & Banerjee, 2015; Fogerty et al., 2016) in which various dynamical and thermal instabilities generate turbulence at the collision interface and the inflowing WNM from both sides injects energy into the system. These CF simulations are also able to reproduce the typical velocity dispersion and column densities observed in the MCs (e.g. Klessen & Hennebelle, 2010; Heitsch et al., 2011; Valdivia et al., 2016).

Some of these studies have discussed the influence of the numerical resolution required to model the structure and composition of MCs. Federrath (2013) perform very high resolution isothermal TB simulations without gravity and concluded that a numerical resolution of 10243 or higher is needed to get a converged density distribution with $<20\%$ difference from the “infinite-resolution limit". Gong et al. (2018) conclude that a resolution of $\sim 2$ pc is sufficient to obtain a converged CO-to-H2 conversion factor in the chemical post-processing of their simulation of a region of a galactic disk; however, the convergence of the total CO abundance is not achieved for a spatial resolution of up to $\sim 1$ pc. In the TB simulations by Glover et al. (2010) and Micic et al. (2012), the CF simulations by Valdivia et al. (2016), and the zoom-in simulations of MC evolution from the turbulent ISM by Seifried et al. (2017) – all of which have coupled chemistry networks – a spatial resolution of $\sim 0.1$ pc seems to be necessary and also sufficient to follow the H2 abundance accurately. CO has a comparatively complex chemistry, and its abundance is mostly dependent on the local density and shielding properties. Although many studies with a coupled chemistry network have done some convergence tests, the spatial resolution needed to model the different chemical species in such 3D simulations has not been consistently determined. Since the idea of dynamical modelling with coupled chemistry is flourishing, it is necessary to have a thorough look into the resolution requirements to model the chemical evolution of turbulent MCs. The multiple TB and CF simulations in this paper are designed to investigate these resolution requirements and will serve to indicate whether such requirements have been met in existing studies.

The paper is organized as follows: the numerical methods implemented in the simulations are described in Section 2 and the setups of the TB and CF simulations along with the complete list of simulation runs are described in Section 3. The results of the simulations are presented in Section 4 and the resolution requirements for converged H2 and CO production are then derived in Section 5. Finally, the conclusions are drawn in Section 6 followed by some complementary information in the Appendix.

2 Numerical Methods

The three-dimensional (3D) MHD adaptive mesh refinement (AMR) code Flash4 (Fryxell et al., 2000; Dubey et al., 2008) is used to perform the TB and CF simulations. The included physical processes are introduced below (see Walch et al., 2015, for more details).

2.1 Magnetohydrodynamics

We solve the ideal MHD equations using the five-wave Bouchut MHD solver (Waagan et al., 2011). Additional source terms for self-gravity, turbulent driving, and radiative heating and cooling are handled using operator splitting. Since this project studies the chemical evolution in hydrodynamic simulations, $\mathbf{B}=0$ throughout this paper. The impact of magnetic fields will be discussed in a subsequent paper (Weis et al. in prep). The set of equations reads:

[TABLE]

where $t$ is time, $\rho$ is the density, $\mathbf{v}$ is the gas velocity, $P=(\gamma-1)u$ is the thermal pressure ($u$ is the internal energy density, $\gamma=5/3$ is the adiabatic index of the ideal gas), $\mathbf{B}$ is the magnetic field, $\mathbf{I}$ is the identity matrix, $\mathbf{g}$ is the acceleration due to self-gravity, $\mathbf{F}$ is the random, turbulent acceleration introduced for the TB simulations ($\mathbf{F}=0$ for the CF simulations), $E$ is the total energy density and $\dot{u}_{\text{chem}}$ is the net rate of change of internal energy due to the heating from the attenuated uniform background interstellar radiation field (ISRF) as well as due to radiative cooling.

2.2 Self-gravity

The acceleration due to gas self-gravity $(\mathbf{g}=-\nabla\Phi)$ is obtained by solving the Poisson equation for the gas

[TABLE]

where $\Phi$ and $G$ refer to the gravitational potential and the gravitational constant, respectively. The Gravity module by Wünsch et al. (2018) calculates $\mathbf{g}$ with an OctTree-based method. We use the geometric Multipole Acceptance Criterion (MAC) introduced by Barnes & Hut (1986) with a limiting angle of $\theta_{\text{lim}}=0.5$. In this project, self-gravity is activated after 10 Myr for the TB simulations, i.e. after the turbulence has developed for $\sim 3$ crossing times. The CF simulations, on the other hand, incorporate self-gravity from the very beginning.

2.3 Dust shielding and molecular (self-)shielding

The attenuation of the ISRF due to dust, H2 and CO (self-)shielding and the shielding of CO by H2 in every single cell in the computational domain is incorporated using the OpticalDepth module (Wünsch et al., 2018). The algorithm is similar to the Treecol algorithm (Clark et al., 2012a) and has been used in Gatto et al. (2015); Walch et al. (2015); Girichidis et al. (2016); Gatto et al. (2017); Peters et al. (2017) and Seifried et al. (2017). When constructing the column density and extinction maps for each cell, we use the recommended number of Healpix pixels of $N_{\text{PIX}}=48$ corresponding to $\theta_{\text{lim}}=0.5$ (see Wünsch et al., 2018). Furthermore, only the contributions from gas within a sphere with radius 16 pc from the cell of interest are included. This radius corresponds to half the box side length in the TB simulations and half of the y– and z– side length for the CF simulations and avoids duplicate contributions from the same gaseous sub-structure across periodic boundary conditions (see Section 3).

2.4 Turbulence

In the TB simulations, turbulent velocities are driven in the gas using a random forcing field, which is evolved in time via the Ornstein-Uhlenbeck process (Konstandin et al., 2015). In Fourier space, the forcing acceleration ($\mathbf{F}$) in Equations 2 and 3 evolves as

[TABLE]

where $\mathbf{\tilde{F}}(\mathbf{k},t)$ is the Fourier transform of the forcing acceleration as a function of wave-vector $\mathbf{k}$ and time $t$, $\tilde{F_{0}}(\mathbf{k},t)$ is the forcing amplitude of the $\mathbf{k}$ mode, $\mathbf{W}(t)$ is the Wiener process that introduces Gaussian random perturbations in $\tilde{F_{0}}(\mathbf{k},t)$, $\mathcal{P}^{\zeta}(\mathbf{k})$ is the projection operator that projects the perturbed Fourier amplitudes to $\mathbf{k}$ in order to obtain a desired mixture of solenoidal to compressive modes (here we choose a thermal mix of modes), and $T_{\text{ac}}=L/v_{\text{rms}}$ is the auto-correlation time. The last term represents the exponential decay of the Fourier modes. Turbulence is introduced on the box scale such that $\tilde{F_{0}}(\mathbf{k},t)\neq 0$ only for $|\mathbf{k}|=1$, where $|\mathbf{k}|$ is measured in units of $2\pi/L$, such that the largest eddy is formed on the box scale with $L=32$ pc (see Section 3.1). The forcing mechanism is configured to maintain the mass-weighted 3D root-mean-square velocity $v_{\text{rms}}$ of the gas at 10 km s*-1* , and to re-normalise the bulk velocity of the simulation domain to zero in every forcing step throughout the simulation time. Therefore, we set $T_{\text{ac}}=L/v_{\text{rms}}=3.13$ Myr.

2.5 Chemical Model, Heating and Cooling

The chemical network used in this project is the “NL99" network (see Glover & Clark, 2012, for more details), which combines the hydrogen chemistry taken from Glover & Mac Low (2007a, b) and the CO network from Nelson & Langer (1999). The network tracks the abundances of 14 chemical species plus free electrons (non-equilibrium evolution of H2, H+, C+, CHx, OHx, CO, HCO+, He+, M+ while atomic H, He, C, O, M are calculated via conservation laws). Here, CHx and OHx refer to intermediate carbon-bearing or oxygen-bearing species such as CH, CH2, OH, H2O or OH*+, and M+ to ionized metals. This network considers multiple pathways of CO formation via CHx* and OHx species, and CO destruction via photodissociation as well as via reaction with ions such as H${}_{3}^{+}$ and He+ that are very sensitive to the cosmic ray ionization rate (Bisbas et al., 2015) and the X-ray energy density (Mackey et al., 2018).

Cooling and heating due to chemical reactions, the photoelectric effect on dust grains, UV radiation, X-ray, and cosmic rays via the diffuse ISRF with strength $G_{0}=1.7$ in units of the Habing field (Habing, 1968; Draine, 1978) are taken into account. The relative abundances (number of atoms relative to total H nuclei) of total He, C, O, and M are $x_{\text{He}}=0.1$, $x_{\text{C}}=1.4\times 10^{-4}$, $x_{\text{O}}=3.2\times 10^{-4}$, and $x_{\text{M}}=1.0\times 10^{-7}$, respectively (Sembach et al., 2000; Szűcs et al., 2014). The cosmic ray ionization rate of atomic hydrogen is $\zeta_{\text{H}}=3\times 10^{-17}$ s*-1* and the cosmic ray ionization rate for other chemical species are then scaled as $\zeta_{\text{H}_{2}}=2\times\zeta_{\text{H}}$, $\zeta_{\text{He}}=1.09\times\zeta_{\text{H}}$, and $\zeta_{\text{C}}=3.83\times\zeta_{\text{H}}$ using the factors from Liszt (2003) based on data from the UMIST111http://www.rate99.co.uk database.

3 Simulation Setup

3.1 Turbulent Box (TB) simulation

The TB simulations are performed in a (32 pc)3 box – shown in Figure 1 – with periodic boundary conditions for hydrodynamics and gravity. Typically, turbulent box simulations start with uniform atomic gas: Glover & Mac Low (2007a) performed TB simulations with three different initial number densities $n_{0}=10,\;30,$ and 100 cm*-3* , Glover et al. (2010) with $n_{0}=300$ cm*-3* , Micic et al. (2012) with $n_{0}=30$ and 300 cm*-3* , Glover & Clark (2012) with $n_{0}=100$ and 300 cm*-3* , and Padoan et al. (2017) simulated supernova-driven turbulence in a periodic box with $n_{0}=5$ cm*-3* . The high density medium naturally forms molecular gas much faster than the low density medium.

The initial densities of H nuclei in the TB simulations presented here span two orders of magnitude: $n_{0}=3$ cm*-3* (low density), $n_{0}=30$ cm*-3* (intermediate density), and $n_{0}=300$ cm*-3* (high density). Table 1 shows the initial conditions, namely the equilibrium values of temperature, sound speed, and the abundance of the chemical species relative to H nuclei in a quiescent medium (without gravity, shielding, and turbulence) for each initial density considered in the simulations. The equilibrium conditions are calculated self-consistently by running the chemical network on a one-zone model with density $n_{0}$ for 104 Myr and assuming no shielding of the gas. The gas remains atomic with all of the hydrogen in H and all of the carbon atoms as C+ ions.

The isothermal sound speed, $c_{\text{s}}$, in Table 1 is calculated via

[TABLE]

where $k_{B}$ is the Boltzmann constant, $T$ is the temperature of the gas in equilibrium, $\mu=1.27$ is the mean molecular weight of the gas, and m${}_{\text{p}}$ is the proton mass.

The simulations are named as TB-R[x]-n[ddd], where TB denotes the Turbulent Box model, R[x] denotes the refinement level x that corresponds to a certain uniform resolution and n[ddd] denotes the initial number density of H nuclei in the homogeneous gas. The TB simulations are evolved with four different uniform resolutions ranging from the cell size of $\Delta x=0.5$ pc ($64^{3}$; R4) to $\Delta x=0.063$ pc ($512^{3}$; R7). For example, TB-R7-n030 is a turbulent box simulation with $n_{0}=30$ cm*-3* and refinement level 7 (i.e. $512^{3}$). Note that the TB simulations do not use the AMR capabilities of Flash; however, the notation of the refinement level (e.g. R4) is used to be consistent with the nomenclature for the CF simulations. Supersonic turbulence is constantly driven in the simulation box to maintain $v_{\text{rms}}=10$ km s*-1* that results in the turbulent crossing time $T_{\text{ac}}=3.13$ Myr (see Section 2.4). TB simulations with lower $v_{\text{rms}}$ produce lower molecular content (Micic et al., 2012). To allow for a direct comparison of the different resolution runs, the same turbulence driving pattern is used in all the simulations with the (32 pc)3 box.

In addition, another set of test simulations with $n_{0}=30$ cm*-3* and $v_{\text{rms}}=1,2,3,6$ km s*-1* are performed in a smaller (8 pc)3 box. These simulations have the uniform spatial resolution corresponding to the TB-R[5,6,7]-n030 runs, and are designed to test the effect of $v_{\text{rms}}$ on the convergence of the H2 chemistry (only used in Section 5.4). These 12 special simulations have a “-v[x]" tag at the end, which indicates the driving velocity, and are run only up to 10 Myr.

All these hydrodynamic simulations exclude magnetic fields and self-gravity is activated for the simulations only after 10 Myr, when the turbulence is well-developed. The dense MC regions form and evolve in the TB simulations up to 20 Myr. The runs with $n_{0}=300$ cm*-3* are stopped shortly after 10 Myr because the medium is so dense that some regions collapse immediately due to gravity (further explanation in Section 4.1). The complete list of the TB simulations is given in Table 2.

3.2 Colliding Flow (CF) simulation

The CF simulations by Heitsch et al. (2006); Vázquez-Semadeni et al. (2007); Clark et al. (2012b); Körtgen & Banerjee (2015) and Valdivia et al. (2016) follow the collision of two streams of WNM with density $n_{0}=1$ cm*-3* , the flow velocity ranging from 5 km s*-1* to 31 km s*-1* , and the domain size ranging from 44 pc to 256 pc. The domain of the CF simulation presented here is a 128 pc $\times$ 32pc $\times$ 32 pc rectangular cuboid with inflow boundary conditions in the $x$–direction and periodic boundaries in the $y$– and $z$–directions (shown in Figure 2). Similarly, the gravity boundary condition is isolated in the $x$–direction and periodic along the $y$– and $z$–directions (see Wünsch et al., 2018, for the implementation of mixed gravity boundary conditions). The whole domain is initially filled with a warm, uniform density medium with $n_{0}=0.75$ cm*-3* . As in section 3.1, we use the chemical network to calculate the equilibrium temperature and chemical composition for this $n_{0}$. We summarize the initial conditions in Table 1.

An irregular interface is constructed at the centre of the simulation box to trigger nonlinear thin-shell instabilities. This is an extension of the one-dimensional sinusoidal interface implementation in Heitsch et al. (2006). The collision interface that separates the flows is created by perturbing the $x=0$ plane using the function

[TABLE]

where $x^{\prime}$ is the $x$–coordinate of the perturbed interface. The parameters that allow one to vary the structure of the collision interface are the amplitude of the perturbation $A$, the wave-numbers along the $y$– and $z$–axis, $k_{y}$ and $k_{z}$, and the phase $\phi$. While the velocities in $y$– and $z$–direction are initially zero, the $x$–velocity on either side of this interface is assigned with a value of $v_{x}=\pm 13.6$ km s*-1* , so that the collision occurs immediately upon the start of the simulation. Figure 3 shows the slices of multiple colliding interfaces at the $x=0$ plane (black regions correspond to positive $v_{x}$, white regions to negative $v_{x}$); Figure 4 shows the H2 and CO content developing in the CF simulations up to 15 Myr for 9 different interfaces. The interfaces with large, homogeneous sub-structures (I1, I2 and I3) do not provide a lot of surface area for the opposite flows to interact. As a result, gas is easily collected in homogeneous “pockets" and much more H2 and CO is produced due to more effective shielding. In all other interfaces that lead to more sub-structure, the molecular content is systematically smaller. Interfaces I5a, I5, and I5b are obtained by changing the amplitude ($A=$ 1 pc, 1.6 pc, and 2 pc, respectively) in Eq. 9, while all other parameters are fixed; the slice of the collision interface at the $x=0$ plane in Figure 3 is identical for all three cases. Figure 4 shows that the molecular content changes slightly with the change in $A$. For the presented CF simulations, the interface I5 with $A=$ 1.6 pc, $k_{y}=2$, $k_{z}=1$, and $\phi=0$, that produces an intermediate amount of molecular gas, is chosen. The system is evolved up to 20 Myr to allow for dense MC cores to form and evolve.

The complete list of CF simulations is given in Table 2. The simulations are named CF-R[x] where CF denotes the Colliding Flow model and R[x] denotes the refinement level x that corresponds to an effective resolution via the AMR technique. The low-resolution CF-R4 and CF-R5 runs have uniform resolution. For all higher resolution runs, the base-grid resolution is $\Delta x=0.25$ pc (corresponding to R5) and the effective resolution is achieved using AMR. For refinement level 6 ($\Delta x=0.125$ pc), we refine on a threshold on the second spatial derivative of the gas density (Paramesh refinement). For higher refinement levels, we use the Jeans criterion (Truelove et al., 1997), which requires that the Jeans length,

[TABLE]

is resolved with at least N${}_{\text{Jeans}}$ cells in one dimension. Here we use N${}_{\text{Jeans}}=8$.

The CF simulation domain is constructed using 4$\times$(32 pc)3 boxes (i.e. blocks) stacked along the $x$–direction, such that the first refinement level (R = 1) contains four blocks, each with $8^{3}$ cells. For each higher refinement level, the number of blocks along each dimension increases by 2; thus, the cell size obtained for refinement level R is

[TABLE]

4 Results

The focus of this paper is to study the changes in the chemical composition in both the TB and CF models with increasing resolution. At first, the basic differences between the results of the TB and CF runs are outlined. Then, the gas distribution as well as the formation of H2 and CO in models with increasing resolution is discussed. In all of the presented results of the CF runs, only the gas with temperature T $<$ 4000 K is taken into account so that the inflowing warm gas is excluded.

4.1 Turbulent Box and Colliding Flows : General differences

The setups of the TB and CF simulations are quite different and thus, the gas evolves differently in each model. The turbulence in the TB runs is continuously driven and energy is injected to maintain a constant velocity dispersion. The turbulence in the CF simulations is a result of various instabilities that have been described in Heitsch et al. (2006) and the energy is injected via the incoming flows. Since the evolution of the gas distribution differs for each initial condition in the different TB and CF runs, the chemical evolution should also be different.

To get a general overview of the different gas morphology produced by the TB and CF runs, their column density maps are compared in Figure 5. The first three maps correspond to the TB simulations with various mean initial densities and the last map corresponds to the CF simulation, all at the same (effective) resolution ($\Delta x=0.063$ pc, R7) and at $t=10$ Myr. The TB-R7-n003, TB-R7-n030 and TB-R7-n300 runs produce gas column densities as high as $10^{2},10^{3}$ and $10^{4}$ M⊙pc*-2*, respectively (note the corresponding scaling shown in the figure). Driven turbulence creates filamentary structures in all TB simulations, that evolve until $t=10$ Myr without self-gravity. On the other hand, the CF simulation with gravity produces comparatively larger, diffuse clumps with dense gas at their centre where the gas column density reaches $\Sigma\sim 10^{2}$ M⊙pc*-2*.

In Figure 6, the evolution of the mass-weighted mean number density ($\langle n\rangle_{\text{mass}}$, top panel) and 3D velocity dispersion ($\langle\sigma\rangle_{\text{mass}}$, bottom panel) of the gas is shown for these four runs (TB-R7-n003, TB-R7-n030, TB-R7-n300 and CF-R7). The vertical, dashed line at $t=10$ Myr denotes that gravity is switched on in the TB simulations after the turbulence is well-developed in the gas for $\sim 3$ crossing times. In all three TB runs, the jumps in $\langle n\rangle_{\text{mass}}$ after 1 Myr arise from local density enhancements when the homogeneous gas first experiences shocks due to the turbulent driving. The value of $\langle n\rangle_{\text{mass}}$ increases by an order of magnitude, and is maintained at this level until gas self-gravity is activated at 10 Myr. In the TB-R7-n003 run, the effect of gravity is not significant. In TB-R7-n030, gravity causes the fragmentation of gas and leads to high density regions. The results of the TB-R[4-7]-n300 runs with gravity are not used for analysis since, shortly after self-gravity is switched on at 10 Myr, most of the mass in these runs (above 30%) is contained in unresolved dense structures that undergo gravitational collapse. In the CF-R7 run, $\langle n\rangle_{\text{mass}}$ experiences a jump at $\sim$ 1 Myr when the colliding flow starts to cool down in the centre. After that, it smoothly increases and remains at $\langle n\rangle_{\text{mass}}\sim 30$ cm*-3* for the first 16 Myr, similar to the low density TB run. Self-gravity dominates the gas dynamics after $t\sim 16$ Myr and increases the mean gas density to $\langle n\rangle_{\text{mass}}\sim 2\times 10^{4}$ cm*-3* , similar to the intermediate density TB run.

The evolution of $\langle\sigma\rangle_{\text{mass}}$ in the bottom panel shows that all the TB runs maintain $\langle\sigma\rangle_{\text{mass}}\approx 10$ km s*-1* once the density inhomogeneities arise after 1 Myr. This is enforced by the turbulent driving described in Section 2.4. On the other hand, for the CF run, the initially high velocity dispersion due to the colliding and cooling gas is gradually lowered to $\langle\sigma\rangle_{\text{mass}}=6$ km s*-1* as the denser regions begin to dominate the gas dynamics (see also Appendix A).

4.2 Turbulent Box results

Figure 7 displays the column densities integrated along the $x$–axis of the total gas, H, H2, C+, C and CO (from top to bottom) for increasing resolution (from left to right) in the TB-R[4-7]-n030 simulations at $t=20$ Myr. The gas has already evolved with gravity for 10 Myr at this point. In order to make the low column density visible, the H column densities are scaled by $\alpha_{\text{H}}=10$. Similarly, the column densities of the carbon species are scaled by a factor

[TABLE]

where $x_{\text{C,tot}}=1.4\times 10^{-4}$ is the relative abundance of carbon atoms in the simulations and $\mu_{\text{C}}=m_{\text{C}}/m_{\text{H}}=12.011$ is the atomic mass of carbon.

Figure 7 shows that both the large scale and the small scale gas morphology changes with resolution. Numerous denser but thinner and more filamentary regions form at higher resolution. The H column density decreases as the resolution increases since more H is converted to H2 that resides either in the dense cores or in the filamentary structures. A low and diffuse C+ column density is found throughout the simulation domain. The column density of C is also low and mainly traces the outskirts of the CO rich regions, showing regions of active photodissociation. In regions of high CO column density, which are only found inside the H2-rich regions, the amount of both C+ and C is negligible. Thus, the distribution of all species greatly differ with increasing resolution.

4.2.1 Convergence study of H2

In this section, the evolution of H2 in the TB simulations with various initial conditions and resolutions are presented. Three different quantities are calculated to investigate the convergence of H2 in the simulations. The mass-fraction of H2 is defined as

[TABLE]

where $M_{\text{H}_{2}}$ is the total mass of H2 and $M_{\text{H,tot}}$ is the total mass of all hydrogen atoms. The second, equivalent expression is the ratio of the number of H atoms in H2 molecules ($2\,N_{\text{H}_{2}}$) and total hydrogen atoms ($N_{\text{H,tot}}$).

The time-scale of H2 formation in the simulation is defined as

[TABLE]

where $\dot{M}_{\text{H}_{2}}$ is the H2 formation rate that is averaged over 2 Myr for every point.

Finally, taking the simulation with the highest refinement level, R${}_{\text{max}}$, as a reference, the results for various resolutions are compared using the ratio

[TABLE]

Figure 8 shows the results of low, intermediate, and high density runs (denoted by n003, n030, and n300) in the first, second, and third column, respectively. Since the first density inhomogeneities due to the turbulence appear after 1 Myr (see Figure 6), the evolution of the molecular content starting from 2 Myr is shown. The vertical, dashed lines denote that gravity is switched on at $t=10$ Myr. Note that the TB-R[4-7]-n300 runs evolve only for a short time after 10 Myr due to the reasons mentioned in Section 4.1. The first row shows the evolution of $f_{\text{H}_{2}}$, with a maximum value of 1 when 100% of the H atoms are in the form of H2. In all three setups with different initial densities, the higher resolution runs generally have higher $f_{\text{H}_{2}}$ at any given time. One exception to this is the TB-R6-n003 run, that contains more H2 than the TB-R7-n003 run after 10 Myr. The runs with low and intermediate density do not show any significant gravitational effects on the evolution of $f_{\text{H}_{2}}$ for any resolution. In the high density runs, the values of $f_{\text{H}_{2}}$ for various resolutions look similar towards the end of the simulations. This is trivial because the molecular content in these high initial density runs saturates and the gas becomes fully molecular regardless of the resolution. The large variation in the evolution of $f_{\text{H}_{2}}$ for different resolutions derives directly from the changing gas morphology, as shown in Figure 7.

The second row shows the evolution of $\tau_{\text{H}_{2},\text{form}}$ for various resolutions. The dashed lines in each plot denote the mass-weighted mean value of the theoretical estimate for the equilibrium H2 formation time-scale $\langle\tau_{\text{H}_{2},\text{eq}}\rangle_{\text{mass}}$, which is given by Hollenbach et al. (1971); Hollenbach & McKee (1989); Goldsmith & Li (2005)

[TABLE]

where $\bar{n}$ is the mean gas number density in the region of interest. In all of the low, intermediate, and high density runs, the H2 formation at all resolutions is faster than $\langle\tau_{\text{H}_{2},\text{eq}}\rangle_{\text{mass}}$ before the saturation level is reached. This “faster-than-average" H2 formation in turbulent media has been shown by Glover & Mac Low (2007b) and recently noted by Micic et al. (2012), Valdivia et al. (2016), and Seifried et al. (2017). For the low density runs in the first column, $\tau_{\text{H}_{2},\text{form}}\rightarrow\infty$ when H2 is being destroyed. There is a common trend for the runs in all three columns: after the first density inhomogeneities set in and H2 formation starts, $\tau_{\text{H}_{2},\text{form}}$ becomes shorter with increasing resolution. Consequently, the saturation level of $f_{\text{H}_{2}}$ is reached faster for higher resolution runs. After $f_{\text{H}_{2}}$ has reached saturation in the higher resolution runs, $\tau_{\text{H}_{2},\text{form}}$ becomes longer for higher resolution. This is less obvious for the low density runs, but clearly visible in the intermediate and high density runs.

The third row shows the time evolution of $\eta_{\text{H}_{2}}$. In case of the TB runs, R${}_{\text{max}}=$R7 with $\Delta x=0.063$ pc. Note that $\eta_{\text{H}_{2}}$ only reflects the trend of the simulations with increasing resolution, and $\eta_{\text{H}_{2}}=1$ is not the expected converged state. In case of the low density runs, the TB-R4-n003 and TB-R5-n003 runs produce less H2 than the TB-R7-n003 run ($\eta_{\text{H}_{2}}<1$), while the TB-R6-n003 run has $\eta_{\text{H}_{2}}>1$ shortly after 10 Myr. In the intermediate and high density runs, $\eta_{\text{H}_{2}}<1$ for all the lower resolution runs and $\eta_{\text{H}_{2}}\rightarrow 1$ faster with increasing resolution. Hence, although $\eta_{\text{H}_{2}}\rightarrow 1$ with time in all cases, the evolution of total H2 is significantly different for the various resolutions and the formation of H2 in runs with resolutions $\Delta x\geq 0.125$ pc is not converged. Since these simulations are computationally expensive ($\sim 500,000$ CPU-hours for the TB-R7-n030 run), it is currently not possible to go to resolutions as high as achieved in the isothermal TB simulations by, for example, Federrath (2013).

4.2.2 Convergence study of CO

Similar to the study for H2, three quantities are calculated to study the formation of CO. The mass fraction of CO is defined as

[TABLE]

where $\frac{12}{28}M_{\text{CO}}$ is the total mass of carbon in CO molecules and $M_{\text{C,tot}}$ is the total mass of all the carbon atoms. The equivalent expression is the ratio of the number of carbon atoms in CO molecules, $N_{\text{CO}}$, to total carbon atoms, $N_{\text{C,tot}}$.

The CO formation time-scale in the simulation is defined as

[TABLE]

where $\dot{M}_{\text{CO}}$ is the CO formation rate, that is averaged over 2 Myr for each point and the factor 12/28 corrects for the mass of oxygen.

Finally, the values of $f_{\text{CO}}$ for various resolutions are compared using the ratio

[TABLE]

Figure 9 shows the evolution of $f_{\text{CO}}$ in the first row with a maximum value of 1 when 100% of the C atoms are in CO. In the low density runs (first column), there is very little mass contained in the box, the mean density in the turbulent medium is low and any kind of shielding is ineffective, leading to insignificant CO abundances. Nevertheless, the differences with increasing resolution can still be investigated. There is no clear tendency of the changes in CO formation with regard to the resolution. For the intermediate density runs (second column), $f_{\text{CO}}$ increases and its fluctuations with time decrease as the resolution increases. The TB-R5-n030, TB-R6-n030, and TB-R7-n030 runs have similar $f_{\text{CO}}$ towards the end of the simulation, although their gas distribution is quite different, as shown in Figure 7. For the high density runs (third column), CO is initially rapidly produced. With increasing resolution, $f_{\text{CO}}$ saturates faster and at higher values, similar to H2 in Figure 8.

The second row shows the evolution of $\tau_{\text{CO},\text{form}}$. Since CO is easily destroyed, multiple instances of $\tau_{\text{CO},\text{form}}\rightarrow\infty$ is seen. In the low density runs, the evolution of $\tau_{\text{CO},\text{form}}$ is chaotic since the environment suitable for CO formation is changing rapidly due to the turbulent stirring. In the intermediate density runs, the general trend of $\tau_{\text{CO},\text{form}}$ with increasing resolution is difficult to discern but it lies between 10 and 1000 Myr for the various resolution runs. In the high density runs, $1<\tau_{\text{CO},\text{form}}<10$ Myr during the early evolution at $t<4$ Myr; this corresponds to the rapid initial increase in $f_{\text{CO}}$ seen in the first row. Notably, $\tau_{\text{CO},\text{form}}<\tau_{\text{H}_{2},\text{form}}$ during this phase, due to the artificial setup of the cold neutral medium, in which the equilibrium CO formation is found to be shorter than the equilibrium H2 formation time-scale (Seifried et al., 2017; Gong et al., 2018, and see Section 5.1.2). Similar to $\tau_{\text{H}_{2},\text{form}}$, the value of $\tau_{\text{CO},\text{form}}$ increases for higher resolution once $f_{\text{CO}}$ begins to saturate.

The third row shows the time evolution of $\eta_{\text{CO}}$ calculated with respect to R${}_{\text{max}}=$R7. As expected, there are large, random variations in $\eta_{\text{CO}}$ for the low density runs. In the intermediate density runs, $\eta_{\text{CO}}\rightarrow 1$ with time for the TB-R5-n030 and TB-R6-n030 simulations, while the TB-R4-n030 run shows a large offset. The high density runs show that the large initial differences quickly settle down once $f_{\text{CO}}$ saturates and the $\eta_{\text{CO}}\approx 1$ at all resolutions towards the end of the runs. Again, similar to the H2 formation, the evolution of total CO content is not converged for all resolutions ($\Delta x\geq 0.125$ pc ; see Section 5 for explanation).

4.2.3 Density distribution

The effective formation of CO molecules occurs in the dense regions, reasonably well-shielded from external radiation (Glover et al., 2010). Thus, it is necessary to probe the mass-weighted density distribution at various resolutions, as shown in Figure 10. The mass contribution in the $i^{\text{th}}$ density bin is calculated as

[TABLE]

where $\rho_{j}$ and $V_{j}$ are the density and volume of the $j^{\text{th}}$ cell, respectively, and $\text{N}_{\text{cells},i}$ are all cells in the computational domain with density $\rho$ such that

[TABLE]

where $d(\log\rho_{i})$ is the bin size. Then, the mass-weighted density PDF, $P_{\rho}$, for every density bin $i$ is calculated as

[TABLE]

where $M_{\text{tot}}$ is the total mass in the simulation. The distributions are averaged from five simulation snapshots between $t=9.5$ Myr and 10 Myr, i.e. before gravity is activated. Hence, they display the density distributions in a medium with well-developed turbulence. The PDFs of low, intermediate, and high density runs are arranged from top to bottom. In all runs, the maximum density increases by a factor of 2 to 3 for each higher resolution. The low density tail is similar for all resolutions but the peak of the distribution as well as the high density tail of runs with different resolutions show significant differences. Thus, the mass-weighted distribution of the gas density is also not converged for the simulations with $\Delta x\geq 0.125$ pc (R7).

4.3 Colliding Flow results

Figure 11 shows the time evolution of the total gas mass with $T<4000$ K, H2 and CO in the CF-R10 simulation with an effective resolution of $\Delta x=0.008$ pc. The CO mass is scaled by $\alpha_{\text{C}}$ (see Equation 12) to visually compare with the mass of the total gas and H2. After about 1 Myr, the shocked gas at the collision interface starts to cool and rapid H2 formation sets in. At the end of the simulation, $\sim 5\times 10^{3}\,\text{M}_{\odot}$ of H2 has formed, which is about 30% of the total mass. When the rapid H2 formation slows down after $\sim$ 10% of the gas mass is in H2, the formation of CO starts. A significant amount of CO is produced only after 14 Myr. At the end of the simulation, about 3 $\text{M}_{\odot}$ of CO has formed, which accounts for $\sim 15\%$ of the total mass in carbon atoms.

Figure 12 shows the column density of the total gas, H, H2, C+, C, and CO (top to bottom) at $t=20$ Myr, integrated along the direction of the incoming flow ($x$–axis). Four selected resolutions (increasing from left to right) are displayed. The column densities of C+, C and CO are scaled by $\alpha_{\text{C}}$ to make the carbon content visible on the same scale. The total gas column density shows that there is a significant change in the structure of the dense gas as the resolution increases in the CF-R[4-6] runs. However, the CF-R6 run and the highest resolution run (CF-R10) have similar gas distributions except for the densest regions. The clumpy and spatially separated dense regions obtained in the low resolution runs (CF-R4, CF-R5) convert to a network of filaments with compact cores at higher resolutions. The H column density shows the diffuse gas distribution that surrounds the denser H2 regions. A diffuse C+ distribution is seen throughout the simulation domain. The low and extended column density of atomic carbon surrounds the CO-rich regions.

4.3.1 Convergence study of H2

The left column of Figure 13 shows the evolution of H2 in the CF runs. Similar to the results of the TB runs, the evolution of the molecular content is shown starting from 2 Myr. The top panel shows the evolution of $f_{\text{H}_{2}}$ (Equation 13). About 40% of H atoms are found in H2 towards the end of the simulations. Compared to the higher resolution runs, the CF-R4 and CF-R5 runs have lower $f_{\text{H}_{2}}$ before 10 Myr and then higher $f_{\text{H}_{2}}$ towards the end of the simulation. The values of $f_{\text{H}_{2}}$ in CF-R6 and higher resolution runs are similar.

The middle panel shows the evolution of $\tau_{\text{H}_{2},\text{form}}$ (Equation 14). The dashed line depicts the equilibrium H2 formation time scale $\langle\tau_{\text{H}_{2},\text{eq}}\rangle_{\text{mass}}$ given by Equation 16. After the rapid H2 formation phase (i.e. $t<6$ Myr), $\tau_{\text{H}_{2},\text{form}}$ is always shorter than $\langle\tau_{\text{H}_{2},\text{eq}}\rangle_{\text{mass}}$. As the mean density of the MC forming gas gradually increases with time, $\tau_{\text{H}_{2},\text{form}}$ decreases to $<20$ Myr towards the end of the simulation. The CF-R4 and CF-R5 runs have longer $\tau_{\text{H}_{2},\text{form}}$ than the other runs during the early evolutionary phase, but have shorter $\tau_{\text{H}_{2},\text{form}}$ after $t\sim 10$ Myr. There are only minor differences in the evolution of $\tau_{\text{H}_{2},\text{form}}$ for the other high resolution runs.

The bottom panel shows the time evolution of $\eta_{\text{H}_{2}}$ (Equation 15) in order to compare $f_{\text{H}_{2}}$ in the runs with various resolutions. Here, the reference simulation is CF-R10 with R${}_{\text{max}}=$ R10, corresponding to an effective resolution of $\Delta x=0.008$ pc. The CF-R4 and CF-R5 runs have $\eta_{\text{H}_{2}}<1$ at early times but eventually exceed the amount of H2 produced by the higher resolution runs ($\eta_{\text{H}_{2}}>1$) towards the end of the simulation. There is no significant difference in $f_{\text{H}_{2}}$ between the other higher resolution runs.

Overall, the H2 formation is not sensitive to the (Jeans-)refinement of the dense regions; the total H2 abundance and average H2 formation time-scale are converged for an effective resolution of $\Delta x=0.125$ pc (CF-R6) and higher. This is in agreement with the required resolution for H2 chemistry as noted by Glover et al. (2010); Micic et al. (2012); Valdivia et al. (2016) and Seifried et al. (2017) (see Section 5).

4.3.2 Convergence study of CO

The right column of Figure 13 shows the evolution of CO in the CF runs. The top panel shows the evolution of $f_{\text{CO}}$ (Equation 17). The CF-R4 run generally has lower $f_{\text{CO}}$ than the other higher resolution runs in the beginning. With time, $f_{\text{CO}}$ becomes significantly larger for the low resolution CF-R4 and CF-R5 runs. The spikes seen in the low resolution runs (e.g. in CF-R4 at $t=6$ Myr) are absent in higher resolution runs. Between $t=12$ and $16$ Myr, the CF-R6 and CF-R7 runs show some deviation in $f_{\text{CO}}$ from the higher resolution runs.

The middle panel shows the evolution of $\tau_{\text{CO, form}}$ (Equation 18). Only the CF-R4 run shows $\tau_{\text{CO, form}}\rightarrow\infty$ when a significant amount of CO is destroyed after $t\sim 6$ Myr. Before 10 Myr, the CO formation time-scale is extremely long. Afterwards, the CF-R4 and CF-R5 runs have shorter $\tau_{\text{CO,form}}$ than the higher resolution runs. With time, the differences in $\tau_{\text{CO,form}}$ at various resolutions become smaller and all the runs have $\tau_{\text{CO,form}}\sim$ 10 Myr towards the end of the simulations.

The bottom panel shows the time evolution of $\eta_{\text{CO}}$ (Equation 19, with R${}_{\text{max}}=$R10). Before 5 Myr, the CF-R4 and CF-R5 runs have $\eta_{\text{CO}}<1$, but they quickly produce as much CO as the CF-R10 run. After 10 Myr, $\eta_{\text{CO}}\gg 1$ for the CF-R4 and CF-R5 runs. With increasing resolution, the deviation from the reference value decreases. The runs CF-R6 and CF-R7 show fluctuations with $0.2<\eta_{\text{CO}}<2$ and only the CF-R8 and CF-R9 runs show $\eta_{\text{CO}}\approx 1$. Hence, the sub-structure of small, dense regions significantly affects the CO formation; the total CO abundance and the average CO formation time-scale are converged for an effective resolution of $\Delta x=0.032$ pc (CF-R8) and higher.

The formation of a higher amount of H2 and CO in lower resolution simulations is explained by the distribution of the visual extinction, $A_{V}$, experienced by the cells in the simulation, shown in Figure 14. For each cell, $A_{V}$ is the average value over the 48 Healpix pixels, computed by the OpticalDepth module (see Section 2.3). The top panel shows the $A_{V}$-PDF at $t=5$ Myr and the bottom panel at $t=20$ Myr, when both H2 and CO are being produced rapidly. It is clearly seen that the mass of well-shielded gas (around $A_{V}=1$) is similar at early times; however, at $t=20$ Myr, the shielding in lower resolution runs (CF-R4 and CF-R5) has become more effective than in the higher resolution runs. Indeed, the column densities of the CF runs in Figure 12 show that a coarser resolution produces large, clumpy structures that eventually lead to higher shielding of the gas within, resulting in a more effective molecule formation. As more molecules form, the non-linear effects due to (self–)shielding and density further increase the molecular content of the cold gas.

4.3.3 Density distribution

Figure 15 shows the mass-weighted density PDFs (Equation 22) of the CF simulations with various resolutions. The distributions are averaged from five simulation snapshots between $t=19.5$ Myr and 20 Myr. The spike at low density ($\rho\sim 10^{-24}$ g cm*-3* ) belongs to the incoming warm gas that begins to cool below 4000 K near the collision interface. The distribution of the gas below the density of $\rho\sim 10^{-22}$ g cm*-3* (peak of the distribution) is similar for all the resolutions. The high density tails in Figure 15 clearly show the effect of resolution on the distribution of gas. With increasing resolution, the density structures of interest – i.e. shocks or gravitating regions – are better resolved and higher densities are obtained. This occurs because the gas is collected in relatively large cells in a low resolution run, whereas the same gas is compressed into smaller sub-structures for higher resolution runs. It is also apparent that the density at which the PDFs start to deviate from one another increases with increasing resolution. This influences the convergence of the molecular gas mass in the presented runs and will be discussed in the next section.

5 The resolution criterion for the convergence of H2 and CO formation

Qualitatively, the required spatial resolution for converged H2 formation in multiple chemical networks has been shown to be $\Delta x_{\text{H}_{2},\text{form}}\sim 0.1$ pc by Glover et al. (2010), Micic et al. (2012), Valdivia et al. (2016), and Seifried et al. (2017); the CF runs agree with this requirement. A theoretical foundation, however, is still missing. Here, the resolution requirements for H2 and CO formation are discussed using the relevant time-scales and densities.

Three important time-scales should be considered to understand the chemical evolution of the gas, namely the molecule formation time-scale ($\tau_{\text{mol},\text{form}}$), the molecule dissociation time-scale ($\tau_{\text{mol},\text{dissoc}}$), and the local cell crossing time-scale of the gas ($t_{\text{cell,cross}}$). The relations between these time-scales result in two conditions relevant for the converged molecule formation in a simulation.

Condition 1 is a physical condition and is given by

[TABLE]

such that effective molecule formation is possible. If a simulation is not able to resolve the typical density at which Condition 1 is satisfied, then the qualitative correctness of the molecule formation is questionable; i.e. regions that should contain molecular gas are not likely to have any molecules (see Section 5.1).

Condition 2 is a dynamical condition and is given by

[TABLE]

If a simulation resolves the typical densities that are high enough to satisfy Condition 1 but too low to satisfy Condition 2, then the molecule formation may be qualitatively correct, but quantitatively incorrect; i.e. the simulation will correctly predict that the gas in a certain region should be molecular, but might get the molecular fraction or the time to reach the equilibrium wrong. Cells that fulfill the dynamical condition become fully molecular within one cell crossing time regardless of the physical properties of the gas. The dynamical condition is important because it addresses the fact that small dynamical changes in the environment easily change the history of molecule formation in a cell due to the non-linear nature of the turbulent gas evolution (see Section 5.2).

Both conditions are applied to our simulation results in Sections 5.3 and 5.4.

5.1 Condition 1: Physical condition

The physical condition (Eq. 23) equates the molecule formation and dissociation time scale, which can be investigated for any given molecule. It results in a density criterion, indicating which densities need to be resolved by simulation. Here we calculate the density criterion only for H2 and CO. For both molecules, we show that the resolution requirement following from Eq. 23 is weaker than the one following from Eq. 24, but this could be different for more complex molecules.

5.1.1 H2 molecule

The H2 formation occurs predominantly on the surface of dust grains (Hollenbach & McKee, 1989) via the reaction

[TABLE]

with the rate coefficient $k_{165}$ in Glover et al. (2010). Photo-dissociation is the main destruction mechanism of H2 molecules (e.g Draine & Bertoldi, 1996) via the reaction

[TABLE]

with the H2 photodissociation rate given by

[TABLE]

where $G_{0}=1.7$ in Habing units denotes the strength of the interstellar radiation field. The factors $f_{\text{dust,H}_{2}}$ and $f_{\text{shield,H}_{2}}$ account for the dust extinction and H2 self-shielding, respectively. If $f_{\text{dust,H}_{2}}=0$, the radiation incident on the H2 molecules is completely attenuated by dust. Similar holds for $f_{\text{shield,H}_{2}}$ which is a function of H2 column density and gas temperature (see Glover et al., 2010; Walch et al., 2015, for further details). Using these reaction rates, the distribution of the H2 formation time-scale ($\tau_{\text{H}_{2},\text{form}}$) and H2 dissociation time-scale ($\tau_{\text{H}_{2},\text{dissoc}}$) in a simulation can be obtained.

The 2D-histogram in Figure 16 shows how the two time-scales, $\tau_{\text{H}_{2},\text{form}}$ and $\tau_{\text{H}_{2},\text{dissoc}}$ are related with the gas number density in the CF-R10 run at $t=20$ Myr. The colour code shows the H2 mass distribution. The H2-mass-weighted mean of the each distribution is displayed by the black, dashed ($\tau_{\text{H}_{2},\text{form}}$) and dotted ($\tau_{\text{H}_{2},\text{dissoc}}$) lines. Cells with densities $n>10$ cm*-3* have $\tau_{\text{H}_{2},\text{form}}<\tau_{\text{H}_{2},\text{dissoc}}$ and thus, form H2 faster than it is dissociated222Since the photodissociation rate depends linearly on the incident radiation field, the required density will increase in environments with higher $G_{0}$, for instance near star-forming regions.. As a result, almost all of the H2 is found in cells that fulfill the physical condition

[TABLE]

The solid, black line in Figure 16 denotes the linear fit for $\log_{10}\tau_{\text{H}_{2},\text{form}}$($\log_{10}n$) for the cells that satisfy the physical condition. This leads to the power law relation

[TABLE]

which is similar to the equilibrium H2 formation time-scale ($\tau_{\text{H}_{2},\text{eq}}$, Equation 16).

5.1.2 CO molecule

The shielding of CO by H2 is important for CO formation. Therefore, effective CO formation sets in only when enough H2 is formed (Glover et al., 2010, see also Figure 11) and hence the formation of H2 on dust grains is a limiting reaction for CO formation. In addition, the formation of CH${}_{2}^{+}$ via the reactions

[TABLE]

have been identified as the limiting reactions for CO formation (see Glover & Clark, 2012, for further details). The net formation rate of CH${}_{2}^{+}$ is then given by Glover & Clark (2012) as

[TABLE]

where $k_{\text{CH}_{2}^{+},1}$ denotes the rate coefficient for reaction 30, $n_{x}$ denotes the number density of the chemical species $x$, $\zeta_{\text{H}}=3\times 10^{-17}\text{ s}^{-1}$ is the cosmic ray ionisation rate of atomic hydrogen, $k_{\text{CH}_{2}^{+},2}$ denotes the rate coefficient for reaction 31, and $k_{\text{dr}}$ denotes the rate coefficient for the dissociative recombination reaction 33.

When the cosmic ray (Bisbas et al., 2015) and X-ray (Mackey et al., 2018) ionisation rates are moderate, the CO molecules are mostly destroyed by photodissociation via the reaction

[TABLE]

with the reaction rate

[TABLE]

where the factor $f_{\text{dust,CO}}$ accounts for the shielding of CO by the dust and $f_{\text{shield,CO}}$ incorporates the shielding of the CO by H2 as well as CO self-shielding (see Glover et al., 2010; Walch et al., 2015, for further details). Using these reaction rates, the CO formation time-scale ($\tau_{\text{CO},\text{form}}$) and dissociation time-scale ($\tau_{\text{CO},\text{dissoc}}$) can now be determined.

As for H2, the 2D-histograms in Figure 17 show the distribution of the two time-scales, $\tau_{\text{CO},\text{form}}$ and $\tau_{\text{CO},\text{dissoc}}$, as a function of the gas number density in the CF-R10 run at $t=20$ Myr. The colour code denotes the mass of CO in each bin. The CO-mass-weighted mean of the each distribution is displayed by the black, dashed ($\tau_{\text{CO},\text{form}}$) and dotted ($\tau_{\text{CO},\text{dissoc}}$) lines. Given the significant amount of CO in cells with $n\gtrsim 10^{3}$ cm*-3* , it can be concluded that the physical condition for CO chemistry is satisfied by cells with

[TABLE]

The solid, black line in Figure 17 shows the linear fit of $\log_{10}\tau_{\text{CO},\text{form}}$($\log_{10}n$) for the cells that satisfy the physical condition. This leads to the power law relation

[TABLE]

Thus, the typical CO formation time-scale is slightly longer than the H2 formation time-scale (Equation 29).

5.2 Condition 2: Dynamical condition

It is evident that the physical condition in Equation 23 is a necessary but not a sufficient condition for converged molecule formation. With the expression for the typical H2 and CO formation time at hand (Equations 29 and 38), the dynamical condition, which equates the molecule formation and cell crossing time-scales (see Equation 24), can be examined once the typical cell crossing time of the dense gas is determined.

The local (M)HD time step is given by (Courant et al., 1967)

[TABLE]

where $\Delta x$ is the size of the grid cell, $v$ is the velocity of the gas in the cell, and the constant factor CFL $\leq 1$ (see e.g. Derigs et al., 2016). In the following discussion, CFL$=1$ is taken for simplicity. Furthermore, determining $\Delta t$ via the velocity of individual cells is a conservative estimate. The reason is that the chemical abundances are advected with the flow and therefore large advection velocities do not change the density distribution of the gas. Since molecule formation occurs in dense and turbulent gas, it is intuitive to take the mass-weighted mean velocity dispersion $\langle\sigma\rangle_{\text{mass}}$ of such dense gas as an estimate for the typical velocity in molecule forming gas. This results in the typical cell crossing time of the gas, defined as

[TABLE]

Substituting the time-scales from Equations 29 and 40 in Equation 24, the gas number density which fulfills condition 2 for the H2 molecule is

[TABLE]

Similarly, from Equations 38, 40 and 24, the gas number density which fulfills condition 2 for the CO molecule is found to be

[TABLE]

Therefore, the density required for the CO formation process to fulfill the dynamical condition is $\sim 4$ times higher than that for H2. The densities in Equations 41 and 42 are denoted by the intersection of the solid, black line and horizontal,dot-dashed lines in Figures 16 and 17, respectively.

5.3 Resolution requirements for the CF simulations

The density estimates such as in Equations 41 and 42 are sufficient for the theoretical understanding of the resolution criteria. However, it is practical to have an estimate of the spatial resolution required to simulate the formation of a molecular cloud with converged H2 and CO formation. For simulations with self-gravitating gas (e.g. the CF runs), the spatial resolution associated with each density in Equations 41 and 42 can be calculated via the Jeans density (Jeans, 1902)

[TABLE]

where $\lambda_{\text{J}}=N_{\text{J}}\,\Delta x$ is the Jeans length resolved with $N_{\text{J}}$ cells, $G$ is the gravitational constant, and the isothermal sound speed $c_{\text{s}}$ in the dense molecule-forming regions is given by Equation 8. Thus, the relation between the number density $n$ and the spatial resolution $\Delta x$ for a simulation is

[TABLE]

This leads to an expression for the required spatial resolution as a function of the gas number density, gas composition, and temperature

[TABLE]

For simplicity, $N_{\text{J}}=1$ is used for the conversion between the number density and spatial resolution.

5.3.1 Resolution requirement for H2 convergence

Substituting the density from Equation 28 in Equation 5.3, and taking $\mu=1.27$, $N_{\text{J}}=1$, and $T=100$ K for the H2 forming cold atomic gas, the required spatial resolution to satisfy the physical condition for H2 formation is

[TABLE]

Therefore, a low spatial resolution of $\sim 40$ pc is enough to obtain densities for which H2 formation from the warm neutral medium is faster than H2 dissociation.

Similarly, for the dynamical condition, substituting Equation 41 in Equation 5.3, the expression for H2 is

[TABLE]

The molecule formation is effective for the dense gas with $n>10$ cm*-3* (see Section 5.1.1), that has a mean velocity dispersion of $\langle\sigma\rangle_{\text{mass}}\sim 2$ km s*-1* in the CF runs (see Appendix A). Moreover, in the dense regions where the dynamical condition is important, $T=10$ K and $\mu\in[1.27,2.35]$ depending on whether the H2 forming region is locally predominantly atomic or molecular. This results in the required spatial resolution of

[TABLE]

5.3.2 Resolution requirement for CO convergence

Substituting the density from Equation 37 in Equation 5.3, and taking $\mu=2.35,\,N_{\text{J}}=1,\,\text{and~{}}T=10\text{~{}K}$ for the CO forming molecular (H2) gas, the required spatial resolution to satisfy the physical condition for CO formation is

[TABLE]

Therefore, a fairly high spatial resolution is required to obtain densities for which CO formation is faster than CO dissociation.

Similarly, for the dynamical condition, Equations 42 and 5.3 result in

[TABLE]

With $\langle\sigma\rangle_{\text{mass}}\sim 2$ km s*-1* and $\mu$, $N_{\text{J}},T$ as before, the required spatial resolution for converged CO formation is:

[TABLE]

For both H2 and CO, the resolution requirement from the dynamical condition is more strict. These two derived resolution requirements are in agreement with the simulation results of H2 and CO convergence in Section 4.3. For the typical velocity dispersion $\langle\sigma\rangle_{\text{mass}}\sim 2$ km s*-1* , the cell crossing time $t_{\text{cell,cross}}$ for two different resolution runs is shown by the dot-dashed lines in Figures 16 and 17 (cyan: CF-R6 and green: CF-R8). The coloured, vertical lines denote the densities corresponding to the spatial resolution of the two runs (Equation 44 with $\mu=2.35$, $N_{\text{J}}=1$, and $T=10$ K). Figure 16 shows that the CF-R6 run with $\Delta x=0.125$ pc is able to fulfill both the physical and dynamical condition for H2 formation. Similarly, Figure 17 shows that the CF-R8 run with $\Delta x=0.032$ pc is able to fulfill both conditions for CO formation. Note that the density corresponding to the CF-R8 run (vertical, green line) is also comparable to the density at which the PDF begins to deviate from that of the CF-R9 and CF-R10 runs in Figure 15.

5.4 Resolution requirements for the TB simulations

In contrast to the CF runs, the TB runs have uniform resolution and self-gravity is activated only after 10 Myr. Thus the conversion from the densities to the length-scale as for the CF runs is not straightforward. Here, we focus on the resolution requirement for converged H2 formation.

Figure 18 shows that the H2 formation time-scale ($\tau_{\text{H}_{2},\text{form}}$) in the TB-n030-R7 run at 20 Myr is similar to that seen in the CF runs (Figure 16). Therefore, the density expression in Equation 41 is valid for the TB runs as well. However, there is a notable difference in the density at which the physical condition (Equation 23) is satisfied because of the significantly different distribution of the dissociation time-scale ($\tau_{\text{H}_{2},\text{dissoc}}$). The TB-R7-n030 run starts with cold atomic gas ($n_{0}=30$ cm*-3* ) and Figure 18 shows that the physical condition is satisfied for $n>1$ cm*-3* . The coloured, horizontal, dot-dashed lines denote $t_{\text{cell,cross}}$ (Equation 40, with $\langle\sigma\rangle_{\text{mass}}=v_{\rm rms}=10$ km s*-1* by construction) corresponding to two different resolutions – cyan: R6 and green: R8. The coloured, vertical, solid lines refer to the density corresponding to the spatial resolution of R6 and R8, as in Figure 16. It is evident that the TB-R6-n030 run is no longer sufficient to fulfill the dynamical condition. Instead, the distribution indicates that refinement level 8 or higher would be required to fulfill both conditions for the H2 molecule.

For the gravitating gas in the TB-R*-n030 runs at $t=20$ Myr, Equation 47 can be employed to obtain the necessary spatial resolution for converged H2 convergence. For $\langle\sigma\rangle_{\text{mass}}=10$ km s*-1* , this results in

[TABLE]

Such high resolution TB runs are computationally very demanding. Equation 47 indicates that the convergence of H2 chemistry is possible at lower resolutions if the velocity dispersion is lowered (i.e. if $t_{\text{cell,cross}}$ is increased in Figure 18). This can be easily tested in the TB runs since the required $\langle\sigma\rangle_{\text{mass}}$ is a tuneable parameter of the setup.

The study of multiple velocity dispersion values is quite expensive if the usual (32 pc)3 box is used. Therefore, a smaller (8 pc)3 simulation box is taken to carry out the TB-n030 runs with various $\langle\sigma\rangle_{\text{mass}}$ values (see Table 2) and with resolutions corresponding to the TB-n030-R[5,6,7] runs. These test runs evolve for $t=10$ Myr without gravity. There is effectively less mass and less shielding capacity of the gas, and this affects the physical condition (Equation 23). However, the physical condition is still less restrictive than the dynamical condition. Thus, these runs are still valid to test the resolution requirement with various strengths of the supersonic turbulence in the gas.

Figure 19 shows the evolution of $\eta_{\text{H}_{2}}$ (Equation 15) in the test runs. It is evident that the low resolution simulations, that have converged evolution of H2 content at low $\langle\sigma\rangle_{\text{mass}}$, are not able to show convergence as the 3D velocity dispersion increases. The results can be summarized as follows:

• •

H2 convergence is possible even for the TB-n030-R5 run with $\Delta x=0.250$ pc if $\langle\sigma\rangle_{\text{mass}}=1$ km s*-1* .

• •

The TB-n030-R6 run with $\Delta x=0.125$ pc is almost converged, with $\eta\gtrsim 80$% for $\langle\sigma\rangle_{\text{mass}}=2$ km s*-1* or 3 km s*-1* .

• •

When $\langle\sigma\rangle_{\text{mass}}=6$ km s*-1* , $\eta_{\text{H}_{2}}$ is significantly lower than 1 for the TB-n030-R[5,6] runs.

Equation 47 requires $\Delta x_{\text{H}_{2},\text{dynamical}}\lesssim$ 0.4, 0.2, 0.12, and 0.06 pc for $\langle\sigma\rangle_{\text{mass}}=$ 1, 2, 3, and 6 km s*-1* , respectively. Therefore, the numerical result from the test runs is consistent with the spatial resolution given by Equation 47 for various $\langle\sigma\rangle_{\text{mass}}$ in the TB runs.

In addition, Figure 10 shows that the density PDFs in the pre-gravity phase are approximately log-normal, and higher resolution runs sample the high density regions more. According to Equation 41, simulations with higher $\langle\sigma\rangle_{\text{mass}}$ require higher densities to satisfy the dynamical condition. Therefore, if the spatial resolution is high enough to sufficiently sample the critical densities given by Equation 41 and 42, molecule formation converges. This holds even if Equation 5.3, which converts a density to a spatial resolution using the Jeans condition, is not applicable.

6 Conclusion

Recently, the inclusion of chemical networks has become more and more common in hydrodynamical simulations of the interstellar medium. However, a comprehensive study on the spatial resolution required to model different molecules is missing. We present 3D hydrodynamic simulations coupled with a non-equilibrium chemistry network carried out using the MHD code Flash4. The effects of gas self-gravity, turbulence, and diffuse radiative transfer are taken into account. The chemical evolution under solar neighbourhood conditions is studied in two prevalent molecular cloud formation scenarios, (a) a turbulent periodic box and (b) a colliding flow. A large set of simulations with spatial resolution ranging from $\Delta x=0.5$ pc to $\Delta x=0.063$ pc for the turbulent box and up to $\Delta x=0.008$ pc for the colliding flow simulations are performed to investigate the convergence of the evolving total H2 and CO content. For the turbulent boxes we also investigate different initial mean number densities of 3, 30, and 300 cm*-3*.

In all simulations, the morphology of the gas changes significantly with increasing resolution on both small and large scales. Typically, molecules form initially more slowly in low resolution runs. After some time, however, large clumps of gas form at low resolutions whereas numerous smaller, filamentary structures form at high resolutions. As the (self-)shielding in the large clumps becomes more effective, the low resolution runs eventually produce more H2 and CO than the higher resolution runs. The simulations suggest that the resolution required to obtain the convergence of simple quantities such as the total H2 and CO content is high, even in quiescent environments subject to moderate levels of turbulence. It follows that more complex simulations with e.g. magnetic fields, high radiation fields, or supernova feedback potentially demand an even higher spatial resolution for the convergence of various chemical species than the one derived here.

The resolution requirements we derive follow from the understanding that the simulation needs to resolve typical densities at which (1) the molecule formation and the dissociation time scale are equal, which we refer to as a physical condition; and (2) the molecule formation and the cell crossing time scale are equal, which we refer to as a dynamical condition. The dynamical condition is more restrictive than the physical condition for both molecules. It requires the simulation to resolve

[TABLE]

for converged H2 formation and

[TABLE]

for converged CO formation. We find that the required spatial resolution depends on the composition, temperature, and the strength of turbulence of the MC forming dense gas. In particular, we derive that a spatial resolution required in the colliding flow setup is

$\Delta x\lesssim 0.2\text{~{}pc}$ to model H2 formation 2. 2.

$\Delta x\lesssim 0.04\text{~{}pc}$ to model CO formation

This is consistent with the numerical results of the colliding flow simulations, which are converged in H2 and CO for these spatial resolutions.

In the driven turbulent box with a driving velocity of 10 km s*-1*, neither the total H2 content, nor the total CO content are converged when a spatial resolution up to $\Delta x=0.063$ pc is used. The driven turbulence in warm atomic gas severely affects the evolutionary history of the H2 and CO molecules because any dense region capable of sustaining the molecular gas is quickly disrupted. Similar to the colliding flow runs, the evolution of the CO content is more severely affected by the resolution. We can explain why the molecule formation is unresolved in these runs: condition 2 is not fulfilled because the required spatial resolution scales with one over the mass-weighted velocity dispersion, and therefore a factor of 5 higher resolution would be required to model the formation of H2 and CO.

Acknowledgements

The authors are grateful to the anonymous referee for the insightful comments that strengthened the paper. The authors are grateful to J. C. Ibáñez-Mejía and F. Dinnbier for fruitful discussions. PRJ, SW, and DS acknowledge the support from the DFG through SFB 956 sub-project C5. SW and SDC acknowledge funding by the European Research Council through the ERC Starting Grant “RADFEEDBACK: The radiative ISM" (project number 679852). SW and SCOG acknowledge the support of the DFG via the SPP 1573, “Physics of the Interstellar Medium” (grant number GL 668/2-1). SCOG acknowledges support from the DFG via SFB 881, “The Milky Way System” (sub-projects B1, B2 and B8). SCOG also acknowledges support from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007 - 2013) via the ERC Advanced Grant “STARLIGHT: Formation of the First Stars” (project number 339177). The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago333FLASH is available at http://www.flash.uchicago.edu/. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de) and thank the Regional Computing Center of the University of Cologne (RRZK) for providing computing time on the DFG-funded High Performance Computing (HPC) system CHEOPS as well as support. The authors are also thankful to VisIt (Childs et al., 2012) which allowed to produce the volume-rendered images.

Appendix A The velocity dispersion of the dense gas in Colliding Flow setup

Figure 20 shows the evolution of the mass-weighted mean velocity dispersion of the dense gas with number density $n>10$ cm*-3* for the CF-R10 run. This dense gas contains $\sim 90\%$ of the total mass of the gas with $T<4000$ K. After $t=5$ Myr, the velocity dispersion of $\sim 2$ km s*-1* is maintained during the cloud evolution as the dense regions become progressively more massive.

Appendix B The molecular content in Jeans unresolved regions

Figure 21 shows the evolution of (top) the percentage of the mass contained in the Jeans unresolved regions of the CF runs and (bottom) the amount of molecular gas (H2 and CO) in such unresolved regions. With increasing resolution, the mass of Jeans unresolved regions decrease significantly, i.e. $\sim 30$% in CF-R4 to $\sim 5$% in CF-R6 and higher resolution runs. For all resolutions, such dense regions hold $\sim 70$% of their mass in H2 and CO molecules. Therefore, most of the molecular gas investigated for the resolution criteria are contained in resolved density structures.