SILCC-Zoom: The dynamical and chemical evolution of molecular clouds

TL;DR

This paper presents high-resolution 3D simulations of molecular cloud formation, analyzing their chemical and dynamical evolution, and examining the roles of H2 and CO as tracers within these clouds.

Contribution

It introduces detailed adaptive mesh refinement simulations that resolve cloud formation processes and chemical evolution, highlighting the importance of resolution and cloud definition.

Findings

CO traces dense gas with n ≥ 300 cm⁻³

H₂ appears at lower densities due to turbulence

CO-to-H₂ ratio increases over 2 Myr

Abstract

We present 3D "zoom-in" simulations of the formation of two molecular clouds out of the galactic interstellar medium. We model the clouds - identified from the SILCC simulations - with a resolution of up to 0.06 pc using adaptive mesh refinement in combination with a chemical network to follow heating, cooling, and the formation of H$_2$ and CO including (self-) shielding. The two clouds are assembled within a few million years with mass growth rates of up to $\sim$ 10$^{-2}$ M$_\mathrm{sun}$ yr$^{-1}$ and final masses of $\sim$ 50 000 M$_\mathrm{sun}$. A spatial resolution of $\lesssim$ 0.1 pc is required for convergence with respect to the mass, velocity dispersion, and chemical abundances of the clouds, although these properties also depend on the cloud definition such as based on density thresholds, H$_2$ or CO mass fraction. To avoid grid artefacts, the progressive increase of resolution has to occur within the free-fall time of the densest structures (1 - 1.5 Myr) and $\gtrsim$ 200 time steps should be spent on each refinement level before the resolution is progressively increased further. This avoids the formation of spurious, large-scale, rotating clumps from unresolved turbulent flows. While CO is a good tracer for the evolution of dense gas with number densities $n \geq$ 300 cm$^{-3}$, H$_2$ is also found for $n \lesssim 30$ cm$^{-3}$ due to turbulent mixing and becomes dominant at column densities around 30 - 50 M$_\mathrm{sun}$ pc$^{-2}$. The CO-to-H$_2$ ratio steadily increases within the first 2 Myr whereas $X_\mathrm{CO}\simeq$ 1 - 4 $\times$ 10$^{20}$ cm$^{-2}$ (K km s$^{-1}$)$^{-1}$ is approximately constant since the CO(1-0) line quickly becomes optically thick.

Tables

Table 1. Table 1: Centre (second column) and extent (third column) of the zoom-in regions for runs MC1 and MC2. The starting time t 0 subscript 𝑡 0 t_{0} of the zoom-in elapsed from the start of the SILCC simulation is given in column 4.

run	centre	volume	$t_0$
	[pc]	[pc3]	[Myr]
MC1	(146,126,0)	88 $\times$ 78 $\times$ 71	11.9
MC2	(42,188,0)	97 $\times$ 81 $\times$ 58	11.9

Table 2. Table 2: Overview of the progressive increase of the spatial resolution d x curr subscript 𝑥 curr x_{\mathrm{curr}} for four different total refinement times τ 𝜏 \tau . The times given below refer to the start of the zoom-in simulation at t 0 subscript 𝑡 0 t_{0} .

$\tau$ =	4.5 Myr	2.25 Myr	1.5 Myr	1.0 Myr
d$x_{\mathrm{curr}}$	time	time	time	time
[pc]	[Myr]	[Myr]	[Myr]	[Myr]
2.0	0 - 1.25	0 - 0.625	0 - 0.5	0 - 0.4
1.0	1.25 - 2.5	0.625 - 1.25	0.5 - 1.0	0.4 - 0.65
0.5	2.5 - 3.5	1.25 - 1.75	1.0 - 1.25	0.65 - 0.85
0.24	3.5 - 4.5	1.75 - 2.25	1.25 - 1.5	0.85 - 1.0
0.12	$>$ 4.5	$>$ 2.25	$>$ 1.5	$>$ 1.0

Table 3. Table 4: Effective spatial resolution d x 𝑥 x , corresponding particle mass m SPH subscript 𝑚 SPH m_{\mathrm{SPH}} in an SPH simulation (obtained via Eq.  7 ), and a selection of recent numerical works with that effective resolution. An italic font-style implies that no chemical network is included in the simulations.

d$x$	$m_{\mathrm{SPH}}$	references
d$x$ $\ge$ 3.9 pc	$m_{\mathrm{SPH}}\ge$ 5.4 M☉	Dobbs et al. (2008); Dobbs & Pringle (2013); Duarte-Cabral et al. (2015);
		Hopkins et al. (2013)a; Pettitt et al. (2014); Richings & Schaye (2016a, b)
3.9 pc $>$ d$x$ $\ge$ 1 pc	5.4 M${}_{☉}{}^{}>m_{\mathrm{SPH}}$ $\ge$ 1.4 M☉	Dobbs (2015); Duarte-Cabral & Dobbs (2016); Hennebelle & Iffrig (2014)a;
		Kim & Ostriker (2015)a, Hu et al. (2016)
1 pc $>$ d$x$ $\ge$ 0.5 pc	1.4 M${}_{☉}{}^{}>m_{\mathrm{SPH}}$ $\ge$ 0.79 M☉	Ibáñez-Mejía et al. (2016, 2017)a
0.5 pc $>$ d$x$ $\ge$ 0.24 pc	0.79 M${}_{☉}{}^{}>m_{\mathrm{SPH}}$ $\ge$ 0.33 M☉	Clark et al. (2012b); Smith et al. (2014)
0.24 pc $>$ d$x$	0.33 M${}_{☉}{}^{}>m_{\mathrm{SPH}}$	Renaud et al. (2013, 0.05 pc resolution)a

Equations

$$\rho (z)={\rho }_0\times \text{exp}\left[-\frac{1}{2}{\left(\frac{z}{h_z}\right)}^2\right]$$

$${\tau }_{\mathrm{diss}}={\left(3.3\times 10^{-11}{\mathrm{s}}^{-1}{\mathrm{G}}_0{f}_{\mathrm{shield},{\mathrm{H}}_2}{e}^{-3.5A_{\mathrm{V}}}\right)}^{-1},$$

$$\sigma =\sqrt{\frac{\Sigma (\text{d}m\times (\vec{v}-{\vec{v}}_{\mathrm{COM}}{)}^2)}{\Sigma \text{d}m}},$$

$$10\mathrm{d}x=L_{\mathrm{J}}=c_{\mathrm{s}}\sqrt{\frac{\pi }{G\rho }},$$

$$N_{\mathrm{neigh}}\times m_{\mathrm{SPH}}=M_{\mathrm{J}}=\frac{4\pi }{3}\rho {\left(\frac{L_{\mathrm{J}}}{2}\right)}^3=\frac{{\pi }^{5/2}}{6}\frac{c_s^3}{G^{3/2}{\rho }^{1/2}}.$$

$$m_{\mathrm{SPH}}=10\mathrm{d}x\frac{{\pi }^2}{6}\frac{c_s^2}{N_{\mathrm{neigh}}G}.$$

$$\frac{m_{\mathrm{SPH}}}{1{\mathrm{M}}_{\text{sun}}}\simeq 1.4\left(\frac{\mathrm{d}x}{1pc}\right).$$

$$f_{{\mathrm{H}}_2}=\frac{M_{{\mathrm{H}}_2,\mathrm{tot}}}{M_{\mathrm{H},\mathrm{tot}}}\text{and}{f}_{\mathrm{CO}}=\frac{\frac{12}{28}{M}_{\mathrm{CO},\mathrm{tot}}}{M_{\mathrm{C},\mathrm{tot}}}.$$

$${\tau }_{\mathrm{CO}}=\frac{1}{(5\times 10^{-16}{\mathrm{cm}}^3{\mathrm{s}}^{-1}\times 0.1n)}\simeq 600\mathrm{Myr}{\left(\frac{n}{1{\mathrm{cm}}^{-3}}\right)}^{-1},$$

$${\tau }_{\mathrm{CO},1-0}=\frac{3h{c}^3{A}_{10}}{16\pi k_{\mathrm{B}}{\nu }_{10}^2{T}_{\mathrm{ex}}\Delta v}{N}_{\mathrm{CO}}\left(1-e^{-h{\nu }_{10}/k_{\mathrm{B}}{T}_{\mathrm{ex}}}\right),$$

$${\tau }_{\mathrm{CO},1-0}\simeq 79\left(\frac{\text{km s−1}}{\Delta v}\right)\left(\frac{10\text{K}}{T_{\mathrm{ex}}}\right)\left(\frac{mag}{A_{\mathrm{V}}}\right)\left(1-e^{-5.53/T_{\mathrm{ex}}}\right).$$

Full text

††pagerange: SILCC-Zoom: The dynamical and chemical evolution of molecular clouds–A††pubyear: 2015

1 Introduction

Molecular clouds (MCs) condense out of the warm interstellar medium (ISM) on scales of several 100 pc and host filamentary substructures on sub-pc scales (see e.g. the reviews by André et al., 2014; Dobbs et al., 2014; Klessen & Glover, 2016). They consist of molecular hydrogen H2, which can only be traced indirectly with observations, mostly by means of CO line emission, dust emission, or dust extinction. Since they are observed to have large non-thermal line widths (e.g. Larson, 1981; Solomon et al., 1987), the cloud sub-structure appears to be shaped by supersonic turbulence. The crux, however, is that supersonic turbulence is expected to decay in a crossing time (e.g. Stone et al., 1998; Mac Low et al., 1998; Elmegreen & Scalo, 2004) unless it can be sustained by some physical process, e.g. stellar feedback (see e.g. Elmegreen & Scalo, 2004; Mac Low & Klessen, 2004; Ballesteros-Paredes et al., 2007, for reviews), gravity-driven turbulence (e.g. Vázquez-Semadeni et al., 2008; Klessen & Hennebelle, 2010; Goldbaum et al., 2011; Ballesteros-Paredes et al., 2011; Krumholz & Burkhart, 2016), or magnetorotational instabilities (Kim et al., 2003). With a turbulent dynamical state being observed, the idea of short-lived, dynamically evolving MCs is widely accepted (Mac Low & Klessen, 2004). It is thus likely that the turbulent, internal sub-structure of MCs is imprinted already in the earliest phase of their formation process.

This conception has several important implications. First, the formation and evolution of MCs cannot be modelled in isolation but should be simulated together with the surrounding, multi-phase ISM. Since the turbulent motions within the cloud seem to be inherited from large scale flows (e.g. Brunt et al., 2009), the physical processes driving these flows need to be modelled self-consistently. This includes e.g. thermal instability to seed cold clouds (Vázquez-Semadeni et al., 2000), self-gravity leading to the collapse of cold clouds, gas accretion onto the forming MCs, as well as nearby supernova explosions, which affect the ISM on scales of several 100 pc and larger (Klessen & Glover, 2016).

Second, the chemical abundances of the gas evolve on timescales comparable to the dynamical time of cloud formation and therefore, cloud formation and molecule formation need to be modelled simultaneously and coupled together. The formation timescale of molecular hydrogen roughly scales as $\tau_{\rm H_{2}}\approx 1\;{\rm Gyr}\left(n/1\;{\rm cm}^{-3}\right)^{-1}$ (Hollenbach & McKee, 1979; Glover & Mac Low, 2007a, b; Glover et al., 2010) leading to $\tau_{\rm H_{2}}\sim$ 10 Myr at a mean number density of $n=100\;{\rm cm}^{-3}$. The free fall timescale at this density is $\tau_{\rm ff}=\sqrt{3\pi/(32G\rho)}\approx 5$ Myr, and thus the chemical state of the gas cannot be assumed to be in equilibrium. Glover & Mac Low (2007b) have shown that in turbulent environments H2 can form more rapidly than $\tau_{\rm H_{2}}$ at the respective mean gas density (see also Micic et al., 2012), possibly also enhanced to due turbulent mixing (Valdivia et al., 2016).

Simulating the formation and evolution of a MC in its large-scale environment with a non-equilibrium chemical network is a numerically challenging task due to the large dynamic range of the problem. Moreover, an on-the-fly treatment of molecule formation and of approximate radiative transfer to compute (self-)shielding are challenging and time-consuming, but are becoming feasible (Dobbs et al., 2008; Clark et al., 2012b; Inoue & Inutsuka, 2012; Pettitt et al., 2014; Smith et al., 2014; Seifried & Walch, 2016; Walch et al., 2015; Richings & Schaye, 2016a, b; Hu et al., 2016; Szűcs et al., 2016; Valdivia et al., 2016; Hu et al., 2017). It is also evident that the detailed modelling of the dense, molecular phase is important for the understanding of galaxy formation on larger scales (Naab & Ostriker, 2016).

Smith et al. (2014) study MC formation from kpc scales down to 0.3 pc scales. They investigate the amount of CO-dark molecular gas in a Milky Way-like galaxy and report extremely elongated filamentary clouds with length $>$ 100 pc as observed by, e.g., Goodman et al. (2014) and Ragan et al. (2014). However, these simulations do not account for the effects of gas self-gravity or stellar feedback. When studying synthetic observations of H2, CO, and HI of a Milky Way-like disc simulated by Dobbs & Pringle (2013), Duarte-Cabral et al. (2015) show that models without feedback and self-gravity fail to reproduce the observed properties of the Milky Way. Furthermore, by re-simulating a MC formed in a full galactic disc simulation including chemistry, Dobbs (2015) demonstrate that clouds tend to only approximately follow the linewidth - size relation found by Larson (1981) with a significant scatter. Moreover, synthetic emission maps of this cloud show that CO picks up only small sub-complexes within giant MCs, Duarte-Cabral & Dobbs (2016, but see also ). CO emission estimates thus miss a significant fraction of the molecular H2 content, which agrees with observational results by e.g. Tielens & Hollenbach (1985); van Dishoeck & Black (1988); Xu et al. (2016).

Richings & Schaye (2016a, b) study the chemical evolution of low-mass disc galaxies and their embedded MCs showing that the CO-to-H2 conversion factor can have a significant scatter and that up to 86% of the H2 might be not traceable by CO observations. However, with a mass resolution of 750 M*☉*, they are unable to resolve the detailed substructure of the clouds as well as the accurate formation of H2 and CO. At significantly higher resolution, Hu et al. (2016, 2017) have investigated the formation of H2 in dwarf galaxies and the impact of interstellar radiation.

On somewhat smaller scales ($\sim$ 500 pc), Ibáñez-Mejía et al. (2016) present a zoom-in study of MCs (down to scales of $\sim$ 0.5 pc) that form out of the multi-phase ISM in supernova-driven stratified galactic discs, but do not include a chemical network. They find that with self-gravity turned off initially, the velocity dispersions of the clouds are too low and move onto the Larson relation only once self-gravity is switched on.

Furthermore, using a setup with periodic boundary conditions, Padoan et al. (2016a, b) investigate MC and star formation in supernova-driven periodic boxes (see also Gatto et al., 2015; Li et al., 2015). Like Dobbs (2015) the authors argue that the MCs tend to follow the Larson relation with a significant scatter of about one order of magnitude. However, as pointed out by Rey-Raposo et al. (2015), starting simulations with an already pre-existing cloud can be problematic since it does not match the full complexity of galactic clouds and the potential impact of star formation in the earliest stage.

This long-standing problem of realistic initial conditions for simulations of MC formation and evolution is the motivation for the present study. Although, as discussed before, a wealth of simulations have tackled this problem, they usually lack one or more important aspects: the larger-scale environment, the inclusion of a chemical network, self-gravity, or high ($\sim$ 0.1 pc) spatial resolution. In particular for the latter aspect, so far a detailed resolution study related to the formation of MCs as well as their chemical evolution is lacking in the literature.

In this work we try to remedy the aforementioned shortcomings by making use of the SILCC111SImulating the LifeCycle of molecular Clouds; www.astro.uni-koeln.de/$\sim$silcc simulations (Walch et al., 2015; Girichidis et al., 2016), which model the evolution of the multi-phase ISM in a part of a stratified galactic disc (500 pc $\times$ 500 pc $\times$ $\pm$ 5 kpc) including a chemistry network describing the formation of H2 and CO. The simulations thus provide detailed thermodynamical and chemical conditions, which allow a self-consistent modelling of the formation of MCs out of the diffuse ISM. We make use of a zoom-in strategy which allows us – for the first time – to simultaneously model MCs with

a high, effective spatial resolution of up to 0.06 pc, 2. 2.

their large-scale environment at lower resolution, and 3. 3.

their detailed chemical evolution.

Zoom-in simulations are becoming more prominent in the context of MC formation (e.g. Renaud et al., 2013; Dobbs, 2015; Ibáñez-Mejía et al., 2016, 2017; Kuffmeier et al., 2016) and we address important questions related to the zoom-in strategy in this paper. We will also investigate which resolution is required to accurately model the dynamical and chemical evolution of clouds. Furthermore, we will show that similar to ISM simulations on larger scales (Walch et al., 2015; Girichidis et al., 2016), a simple density threshold criterion does not recover the same cloud structure that we find when we examine the chemical state of the gas. Moreover, we demonstrate that the H2 abundance is enhanced in low-density gas ($n$ $<$ 30 cm*-3*) due to turbulent mixing as previously suggested by Glover & Mac Low (2007b) and Valdivia et al. (2016).

The structure of the paper is as follows: We first present the zoom-in strategy used in this work (Section 2). In Section 3 we discuss MC properties obtained by different selection criteria and the evolution of the H2 content. Next, we investigate which resolution is required to properly model the chemical and dynamical evolution of the MCs in Section 4.2. In Section 4.3, we further show that the time, over which the resolution is increased, can change the structure of the resulting cloud. We show a first application in the context of synthetic observations in Section 5 and conclude in Section 6.

2 Numerical methods and simulation setup

The zoom-in simulations discussed in this work are based on the SILCC simulations presented in detail in Walch et al. (2015) and Girichidis et al. (2016). In the following we briefly describe the numerical methods (Section 2.1) and initial conditions (Section 2.2). The zoom-in technique is introduced in Section 2.3.

2.1 Numerical methods

We use the adaptive mesh refinement (AMR), finite-volume code FLASH version 4 (Fryxell et al., 2000; Dubey et al., 2008) to solve the (magneto-)hydrodynamical equations. The ’Bouchut 5-wave solver’ guarantees positive entropy and density (Bouchut et al., 2007, 2010; Waagan, 2009; Waagan et al., 2011). The resolution on the base-grid is 3.9 pc.

We model the chemistry of the ISM using a simplified chemical network for hydrogen and carbon chemistry (Nelson & Langer, 1997; Glover & Mac Low, 2007a, b; Glover et al., 2010; Glover & Clark, 2012). The network is designed to follow the chemical abundances of H+, H, H2, C+, and CO as well as free electrons and atomic oxygen. We assume solar metallicity, with fixed elemental abundances of carbon and oxygen given by $x_{\mathrm{C}}=1.41\times 10^{-4}$ and $x_{\mathrm{O}}=3.16\times 10^{-4}$, respectively (Sembach et al., 2000). Initially, all carbon is in the form of C+ and all hydrogen is in the form of atomic hydrogen.

The chemical network also describes the thermal evolution of the gas including detailed heating and cooling processes using the chemical abundances provided by the network (Glover et al., 2010; Glover & Clark, 2012; Walch et al., 2015). For the heating via the photoelectric effect as well as photo-dissociation reactions, the simulations include a uniform interstellar radiation field (ISRF) with a strength of G0 = 1.7 (Draine, 1978) in units of the Habing field (Habing, 1968). The ISRF is attenuated due to shielding by the surrounding gas and dust. The necessary column densities of H2, CO, and the total gas are calculated via the TreeCol algorithm (Clark et al., 2012a). The cosmic ray ionisation rate is set to a constant value of 1.3 $\times$ 10*-17* s*-1* in the entire simulation domain. More details of the implementation of the chemical network in FLASH and its coupling to the TreeCol algorithm is explained in Walch et al. (2015) and Wünsch et al. (2017).

Both self-gravity of the gas as well as a background potential, which represents the old stellar component of the galactic disc, are taken into account. For self-gravity of the gas we solve the Poisson equation with a tree based method (Wünsch et al., 2017), which is part of the official FLASH release. The external gravitational acceleration due to the pre-existing stellar component in the galactic disc is modelled with an isothermal sheet with $\Sigma_{\mathrm{star}}$ = 30 M*☉* pc*-2* and a scale height of 100 pc, originally proposed by Spitzer (1942).

Turbulence is generated in the simulations by supernova (SN) feedback. For a single SN explosion, an energy of 1051 erg is injected in the form of thermal energy if the mean density in the surrounding spherical region (i.e. the injection region with a minimum radius of 4 cells) is low enough to resolve the Sedov-Taylor phase of the SN. If the Sedov-Taylor phase is unresolved, we heat the gas within the injection region to $10^{4}$ K and inject the momentum, which the swept-up shell has gained at the beginning of the momentum-driven snowplough phase (see Gatto et al., 2015, for details, but also Walch & Naab 2015; Haid et al. 2016).

The rate at which SNe are injected is based on the Kennicutt-Schmidt relation (Schmidt, 1959; Kennicutt, 1998), which relates the disc’s gas surface density $\Sigma_{\mathrm{gas}}=10\;{\rm M}_{\odot}\;{\rm pc}^{-2}$ (see Section 2.2) to a typical star formation rate surface density $\Sigma_{\mathrm{SFR}}$, i.e. $\Sigma_{\mathrm{SFR}}\propto\Sigma_{\mathrm{gas}}^{1.4}$. In a next step, $\Sigma_{\mathrm{SFR}}$ is translated into a SN rate by assuming a standard initial mass function (Chabrier, 2001). This defines an SN rate with which the gas is heated and stirred.

In the original SILCC simulations (Walch et al., 2015; Girichidis et al., 2016), we investigate how a different positioning of the SNe relative to the dense gas influences the structure of the ISM. We find that mixed driving results in an ISM with solar neighbourhood properties. Mixed driving means that 50% of the SNe are injected at random positions (modulo a weighting of the vertical positions with a Gaussian distribution with a scale height of 50 pc) and 50% are injected at local density peaks. In this way we mimic the explosions of massive stars in either their parental cloud or diffuse regions due to runaway O/B-stars. The simulation with mixed driving is the starting point of the work presented here.

We switch off the SN driving when we start the zoom-in simulations and discuss the impact of nearby SN explosions on the evolution of the forming MC in a subsequent paper.

2.2 Initial conditions of the SILCC simulation

Our initial conditions represent a small section of a galactic disc at low redshift with solar neighbourhood properties and solar metallicity. The simulation box has a size of 500 pc $\times$ 500 pc $\times$ $\pm$ 5 kpc with periodic boundary conditions applied along the $x$- and $y$-direction and outflow conditions along the $z$-direction. The gas surface density is $\Sigma_{\mathrm{gas}}$ = 10 M*☉* pc*-2* and the initial vertical distribution of the gas is modelled with a Gaussian profile

[TABLE]

with $h_{z}$ = 30 pc and $\rho_{0}=9\times 10^{-24}$ g cm*-3*. The resulting total gas mass is $M_{\mathrm{disc}}=2.5\times 10^{6}$ M*☉. The gas near the disc midplane has an initial temperature of 4500 K and consists of atomic hydrogen and C+. At large scale heights the density is set to a constant value of $\rho=10^{-28}$ g cm-3*, representing hot, ionised halo gas with T = 4 $\times$ 108 K.

The gas is initially at rest. From the start of the simulation until the time when we start to zoom-in on the forming molecular clouds, we drive SNe at a constant rate of 15 Myr*-1* using mixed driving (see Section 2.1). In addition to the SN driving, the heating and cooling of the gas plus the gravitational forces self-consistently generate structures and motions within the disc gas.

2.3 Initial conditions of the zoom-in simulations

A number of dense and cold MCs develop from the diffuse ISM in this SILCC simulation (see Walch et al., 2015, for a discussion). These MCs have different histories. They continuously accrete gas and sometimes merge with each other. In addition, they are pushed from different sides by nearby SN explosions. We select some of these clouds by eye and re-simulate them at a high spatial and temporal resolution to study the evolution of their internal sub-structure and chemistry. We call this a ”zoom-in” simulation.

To properly follow the evolution of the clouds’ internal structure, we rewind the simulation and start to zoom-in at a time when the typical number density within the zoom-in region does not exceed a few times 10 cm*-3*. Such low densities are a crucial requirement to properly model the initial formation process of the cloud. This limits the number of possible cloud candidates since in particular at later times the clouds will already be significantly denser. Moreover, in order to allow for a clean convergence study, we choose MCs which form in (relative) isolation and are not affected by any significant merger. In particular at late times the MCs start to merge, which leaves us with very few compact clouds after $\sim$ 100 Myr (see e.g. Figs. 6 and 7 in Walch et al., 2015).

In this paper, we choose two molecular clouds, MC1 and MC2, which form at the same time but at different positions in the galactic midplane. Initially, the precursors of the two MCs are formed from gas swept-up by SN explosions (Koyama & Inutsuka, 2000; Inoue & Inutsuka, 2009; Inutsuka et al., 2015). As discussed later on, the cloud properties are representative for typical Galactic MCs. While the surrounding multi-phase ISM is modelled at a fixed resolution of 3.9 pc, we gradually increase the resolution within the zoom-in region (see below, Section 2.4). The starting time for both zoom-in runs is $t_{0}=11.9$ Myr and the extent of the two zoom-in regions are specified in Table 1. Since the MCs do not have any significant bulk motion with respect to the grid, we do not have to shift the zoom-in region with time.

2.4 Refinement strategy

Flash is an adaptive mesh refinement (AMR) code based on the PARAMESH library (MacNeice et al., 2000) with which the resolution can be increased in a spatially and temporally adaptive way (Fryxell et al., 2000). The user has to define refinement criteria upon which a decision whether the resolution should be increased or decreased is made. Further, in Flash the spatial resolution can only differ by a factor of 2 for neighbouring cells.

The minimum refinement level, with which the SILCC simulation is run before $t_{0}$, corresponds to a base grid resolution of 3.9 pc. To increase the resolution further, we use two refinement criteria. Up to a spatial resolution of 0.5 pc (3 levels above the base grid), we refine on variations in the gas density. This estimate is based on the second derivative of the density field normalised by the average of the gradient (following Lohner, 1987) and allows to capture low-density fluctuations. This criterion is quite sensitive such that the entire region is refined to the maximum possible resolution. Furthermore, we also refine on the local Jeans length, $L_{\rm Jeans}$, which is computed for each cell. We require that $L_{\rm Jeans}$ is resolved with at least 16 cells in each dimension (Truelove et al., 1997; Federrath et al., 2011).

For refinement above 0.5 pc, we only use the Jeans refinement criterion. To carry out a computationally feasible simulation222In collapsing cold cores the minimum Jeans length decreases with increasing density up to $\rho\approx 10^{-13}\;{\rm g\;cm}^{-3}$, where the gas becomes optically thick to infrared radiation. It is unfortunately impossible to resolve such a high density in our simulations., we have to set a global maximum refinement level, $l_{\rm max,tot}$ with a corresponding effective spatial resolution d$x$. In our fiducial simulation, $l_{\rm max,tot}$ corresponds to d$x$ = 0.12 pc (or 5 levels above the base grid). Density structures that form on these scales will not be resolved and we cannot make a prediction on their further fragmentation. However, we show that an effective resolution of d$x$ = 0.12 pc is sufficient to reach convergence in terms of H2 and CO mass fractions within the forming molecular clouds (see Section 4.2).

A stepwise refinement is applied within the zoom-in region. Therefore, we define a refinement time scale, $\tau$, over which the current maximum refinement level, $l_{\rm max,curr}$ (with the current maximum resolution d$x_{\mathrm{curr}}$) is increased from the base grid level at $t_{0}$ to $l_{\rm max,tot}$ at time $t_{0}+\tau$. We test different values of $\tau=$4.5, 2.25, 1.5, 1.0, and 0.0 Myr, where $\tau=0.0$ Myr implies that the code is enabled to instantaneously increase the resolution to $l_{\rm max,tot}$. For $\tau>0$ we spend a certain time on each intermediate maximum refinement level $l_{\rm max,curr}$ and thus, slowly increase $l_{\rm max,curr}$ to $l_{\rm max,tot}$. In Table 3 we list the time spent on each refinement level for different $\tau$. The time spent corresponds to a typical number of time steps for which the code runs with an intermediate $l_{\rm max,curr}$. For $\tau=1.0$ Myr, we typically spend 200 time steps on each level before advancing to the next higher level. According to the CFL condition, a typical time step decreases with decreasing d$x_{\mathrm{curr}}$ and therefore the time spent on an intermediate $l_{\rm max,curr}$ decreases with increasing $l_{\rm max,curr}$. For other values of $\tau$, the times on each intermediate level are scaled up.

We note that so far the choice of $\tau$ is somewhat heuristic. However, in Section 4.3, we show that a refinement time of $\tau$ = 1.5 Myr yields good results, whereas for $\tau$ $<$ 1 Myr – due to the low number of time steps per intermediate $l_{\rm max}$ – numerical artefacts appear. We also perform runs with different maximum refinement levels and determine the resolution needed to reach convergence in cloud mass and chemical abundances in Section 4.2. All runs performed in this work are listed in Table 3.

3 Fiducial runs

In order to study the impact of resolution and refinement time on the simulation results (see Section 4), we have to define certain quantities that are used to compare the simulations. In this section, we therefore discuss the results for runs MC1 and MC2 (see Section 2.3) with the fiducial parameters $\tau$ = 1.5 Myr and d$x$ = 0.12 pc (see Section 2.4). These are called MC1_tau-1.5_dx-0.12 and MC2_tau-1.5_dx-0.12 (see Table 3). Starting from $t_{0}$, every simulation is evolved for 5 Myr, so the final time is $t_{\rm end}=t_{0}+5$ Myr. In all figures shown in this paper, the time refers to the time elapsed since $t_{0}$.

3.1 Overview

In Fig. 1 we show the edge-on view of the total gas column density $\Sigma$ of run MC1_tau-1.5_dx-0.12 at time $t_{\rm end}$ on different spatial scales of (from left to right) 500 pc, 100 pc, and 16 pc. The figure demonstrates the power of the zoom-in technique, which allows us to cover spatial scales over about 5 orders of magnitude to resolve the filamentary structure on pc scales and densities over about 8 orders of magnitude from the diffuse ISM to the densest structures. Further, in Fig. 2 we show the time evolution of the column density of run MC1_tau-1.5_dx-0.12 (top) and run MC2_tau-1.5_dx-0.12 (bottom). Neither of the clouds has a simple spherical shape. MC1 has an almost “T”-shaped structure with significant filamentary and clumpy substructure developing over time. While MC1 is dominated by one massive central clump, MC2 fragments into several clumps. Furthermore, we find extended low-(column-)density regions, which are created by nearby SNe that have exploded before $t_{0}$.

In Fig. 3 we show the column density maps of atomic hydrogen, H2, and CO at the end of both runs. We find that atomic hydrogen envelopes the molecular component, whereas CO is mainly present in the densest parts of the clouds (Smith et al., 2014; Duarte-Cabral & Dobbs, 2016; Xu et al., 2016). We note that the chemical model used here, might slightly overestimate the production of CO (Glover & Clark, 2012). However, we do not expect a more accurate chemical model to qualitatively change the aforementioned findings.

Also, we depict the evolution of the total mass and the mass in the five species which are traced with our chemical network in Fig. 4. The total mass in the zoom-in region shows only minor variations over time, which indicates that we are only marginally affected by the flows through the boundaries of the zoom-in region. The H2 mass in both MCs is 3 – 4 $\times$ 104 M*☉* towards the end. The mass in H+ is negligible compared that of H and H2. For both, hydrogen and carbon, the molecular component starts to dominate after about 3 Myr. However, the growth rates of the CO and H2 masses are different. At $t_{0}$ approximately 104 M*☉* of H2 are already present while the mass of CO – which forms at higher visual extinction, $A_{\mathrm{V}}$, than H2 (e.g. Röllig et al., 2007; Glover et al., 2010) – is still negligible. In the following 2 Myr the H2 mass shows only a moderate increase by a factor of 2 – 3, whereas the CO mass increases by a factor of 100 once the high density regions start to form. At later stages the evolution of both species seems to be similar. We discuss the resulting time-variable CO-to-H2 ratio in more detail in Section 4.2.2 and investigate its effect on the conversion factor between measured CO, $J$ = 1 – 0 line intensity and H2 mass, $X_{\mathrm{CO}}$ (see e.g. Shetty et al., 2011a, b; Bolatto et al., 2013) in Section 5.

3.2 Molecular cloud parameters

We consider the time evolution of global properties like mass, mean density, and velocity dispersion of the simulated MCs. In this context, a basic question is how to define a MC. In numerical works this has often been done by using a volume density (e.g. Tasker & Tan, 2009; Ibáñez-Mejía et al., 2016) or column density threshold (e.g. Dobbs, 2008; Ward et al., 2014; Bertram et al., 2016). Motivated by observations (e.g. Goldsmith et al., 2008; Heyer et al., 2009), recent studies which include the chemical evolution of MCs also use intensity thresholds for CO line emission (Clark et al., 2012b; Smith et al., 2014; Padoan et al., 2016a; Richings & Schaye, 2016a). In this paper, we use five simple definitions based on the 3D structure and chemistry of the cloud:

• •

three “volume density thresholds”, where all gas above a threshold value is counted as part of the MC. We consider the three density thresholds $\rho_{\rm thres}=$ 1.15 $\times$ 10*-22*, 3.84 $\times$ 10*-22*, and 1.15 $\times$ 10*-21* g cm*-3*. Under the assumption of a mean molecular weight of 2.3, this corresponds to particle densities of $n_{\rm thres}=$ 30, 100, and 300 cm*-3*. In all figure labels we refer to $n_{\rm thres}$.

• •

an ”H2 threshold”, where we sum up the total gas mass in cells with a mass fraction of H2 – with respect to the total hydrogen mass in that cell – that is greater or equal to 50%, and

• •

a ”CO threshold”, where we sum up the total gas mass in cells with a mass fraction of carbon locked up in CO – with respect to the total mass in carbon in that cell – that is greater or equal to 50%.

In Fig. 5, we plot the masses (top row), mean number densities $\bar{n}$ (middle), and velocity dispersions $\sigma$ (bottom row) for MC1 (run MC1_tau-1.5_dx-0.12, left column) and MC2 (run MC2_tau-1.5_dx-0.12, right column) as a function of time. The time refers to the start of the zoom-in simulation at $t_{0}$. Each of these quantities is derived for the 5 different criteria. We refer to each criterion via a subscript, e.g. the masses for the three $n_{\rm thres}$ thresholds are $M_{30}$, $M_{100}$, $M_{300}$, and for H2 and CO criterion $M_{\mathrm{H_{2}}}$, and $M_{\mathrm{CO}}$.

3.3 Mass assembly and evolution of the mean density

The MCs condense out of the atomic ISM within $\sim$ 2 Myr from the start of the zoom-in simulation at $t_{0}$. Therefore, the masses above the given density thresholds, $M_{30}$, $M_{100}$, and $M_{300}$ grow within this period of time. Depending on $n_{\rm thres}$, the final masses of MC1 are between $M_{30}\approx 7.0\times 10^{4}\;{\rm M}_{\odot}$ and $M_{300}\approx 4.3\times 10^{4}\;{\rm M}_{\odot}$. For MC2 we reach masses between $M_{30}\approx 4.0\times 10^{4}\;{\rm M}_{\odot}$ and $M_{300}\approx 2.8\times 10^{4}\;{\rm M}_{\odot}$. Overall, $M_{300}$ accounts for about 40% of the total mass within the zoom-in region in both cases. The derived masses and sizes (see Fig. 2) indicate that these clouds are representative of typical MCs in the Milky Way (e.g. Solomon et al., 1987; Elmegreen & Falgarone, 1996; Heyer et al., 2001; Roman-Duval et al., 2010; Miville-Deschênes et al., 2017).

Using the H2 and CO thresholds, $M_{\mathrm{H_{2}}}$ and $M_{\mathrm{CO}}$ start to increase later than $M_{300}$ and show a different time dependence. Again, we remind the reader that $M_{\mathrm{H_{2}}}$ and $M_{\mathrm{CO}}$ are not the total mass in H2 and CO (as shown in Fig. 4), but those obtained by the criteria discussed before. $M_{\mathrm{CO}}$ grows more slowly than $M_{300}$, but steadily increases with a mean rate of $5-7\times 10^{-3}$ M*☉* yr*-1*. Therefore, $M_{\mathrm{CO}}$ approaches $M_{300}$ towards the end of the simulations. We note that initially $M_{\mathrm{H_{2}}}$ is close to zero, although about 7 $\times$ 103 M*☉* of H2 are present within both zoom-in regions at that time (see Fig. 4). This is due to the fact that at this stage hydrogen is predominantly in atomic form, and thus the diffuse H2 component is not taken into account in $M_{\mathrm{H_{2}}}$. Once the cloud starts to become molecular, $M_{\mathrm{H_{2}}}$ shows a higher growth rate than $M_{\mathrm{CO}}$ with $\dot{M}_{\mathrm{H_{2}}}\sim 1\times 10^{-2}$ for run MC1_tau-1.5_dx-0.12 (and $\dot{M}_{\mathrm{H_{2}}}\sim 7\times 10^{-3}$ M*☉* yr*-1* for MC2_tau-1.5_dx-0.12) and thus, overtakes $M_{300}$ after $\sim$ 2.5 Myr as well as $M_{100}$ after $\sim$ 3.5 Myr. These growth rates are similar to recent observations (e.g. Fukui et al., 2009; Kawamura et al., 2009) but the different time dependence of $M_{\mathrm{H_{2}}}$ and $M_{\mathrm{CO}}$ supports the idea of a time-dependent X-factor (see Fig. 4 and Section 4.2.2 and e.g. Richings & Schaye, 2016a, b).

The mean, volume-weighted densities $\bar{n}$ (middle row of Fig. 5) inferred using the $n_{\rm thres}$ criteria show an increase over time. Given that the masses, $M_{30}$, $M_{100}$, and $M_{300}$, differ only by a factor of $\sim$ 2, the differences in $\bar{n}$ by more than one order of magnitude reflect that the occupied volume, i.e. the volume of all cells which fulfil the applied criterion, increases significantly with decreasing density threshold.

For the CO threshold, we find that $\bar{n}_{\rm CO}$ is similar to, but slightly higher than, $\bar{n}_{300}$. On the other hand, $\bar{n}_{\rm H_{2}}$ decreases as a function of time and even drops below $\bar{n}_{30}$ after $\sim$ 4 - 4.5 Myr. The decrease of $\bar{n}_{\mathrm{H_{2}}}$ is due to the progressive extension of H2 into the low-density regime. This is in agreement with our analysis so far, which implies that gas with $n<100\;{\rm cm}^{-3}$ can become predominantly molecular and shows that the volume occupied by this gas significantly increases over time.

3.4 Evolution of H2

In Fig. 6 we show the mass fraction of H2 with respect to the total amount of hydrogen in the dense ($n$ $>$ 100 cm*-3*) gas of MC1 and MC2. Initially, the gas is predominantly atomic, with an H2 mass fraction of $\sim$ 0.2 – 0.3. After about 1.5 Myr it increases more rapidly reaching an almost fully molecular state with an H2 mass fraction of about 0.9 at $t$ = $t_{0}$ + 5 Myr. We note that for density thresholds of 30(300) cm*-3* the evolution is qualitatively very similar with 5 – 10% lower(higher) H2 fractions. In Fig. 7 we plot the ratio of the H2 to atomic hydrogen column density, $\Sigma_{\mathrm{H_{2}}}/\Sigma_{\mathrm{HI}}$, against the total column density for MC1 at $t$ = $t_{0}$ + 5 Myr. The transition to H2 dominated gas occurs around 30 – 50 M*☉* pc*-2* ($\Sigma_{\mathrm{H_{2}}}/\Sigma_{\mathrm{HI}}$ $>$ 1, indicated by the blue line representing the mean value at each $\Sigma$). This is somewhat higher than the value of 10 M*☉* pc*-2* found in the Perseus cloud (Lee et al., 2012, 2015) and that of theoretical predictions for solar metallicity by Krumholz et al. (2009, their Equation 39, green line in Fig. 7). It agrees, however, well with recent observations of W43 (Bihr et al., 2015). We note that for other times and MC2 we find a very similar behaviour.

To further analyse the properties of the H2 gas, we show a density slice through the centre of MC1 at $t$ = 5 Myr in Fig. 8. The contours show the volume occupied by gas above $n_{\rm thres}=$ 30, 100, and 300 cm*-3* (blue, green, and cyan lines) and above the H2 threshold (red line). It can be seen that the H2 dominated gas extends into regions with gas densities below 30 cm*-3*. The presence of H2 in such low-density gas is difficult to explain with in-situ formation of the molecules in this regions: The time scale $\tau_{\mathrm{form}}$ for H2 formation in gas with $n=30$ cm*-3* is about 1 Gyr $\times\left(\frac{30\,\mathrm{cm}^{-3}}{1\,\mathrm{cm}^{-3}}\right)^{-1}\approx 33$ Myr (Hollenbach & McKee, 1979; Glover & Mac Low, 2007a, b; Glover et al., 2010), and thus significantly longer than the entire simulated time.

We rather argue that H2 forms in dense regions and is efficiently re-distributed into the lower density environment by turbulent mixing. The turbulent mixing timescale $\tau_{\mathrm{mix}}$, on which this is expected to happen, can be estimated by assuming that the typical distance $l_{\mathrm{mix}}$ of H2 in the low-density environment from the high density clumps is a few pc (Fig. 8) and the typical velocity $v_{\mathrm{mix}}$ a few km s*-1* (see bottom row of Fig. 5). Hence, we obtain $\tau_{\mathrm{mix}}=l_{\mathrm{mix}}/v_{\mathrm{mix}}\simeq$ 1 Myr $<<$ $\tau_{\mathrm{form}}$. Furthermore, the dissociation timescale of H2 due to UV radiation, which is the dominant dissociation process, is given by

[TABLE]

where $f_{\mathrm{shield,H_{2}}}$ is the UV shielding due to H2 molecules (van Dishoeck & Black, 1988; Glover et al., 2010). In the regions with $n$ $\gtrsim$ 30 cm*-3* we typically find $A_{\mathrm{V}}$ $>$ 0.3 and $f_{\mathrm{shield,H_{2}}}$ $<$ 3 $\times$ 10*-4*, which results in $\tau_{\mathrm{diss}}$ $\gtrsim$ 5 Myr $>$ $\tau_{\mathrm{mix}}$. Therefore, the relatively short time scale of a few Myr on which H2 starts to dominate in gas with $n\lesssim 30\;{\rm cm}^{-3}$ implies that the mixing of H2 into the low density gas is significantly more important than the in-situ formation of the molecules. This result is similar to the enhanced H2 formation process reported in more idealised simulations of molecular cloud formation using turbulent periodic boxes (Glover & Mac Low, 2007b) and simulations of isolated clouds (Valdivia et al., 2016). We note that a part of the H2 found in low-density regions might also be caused by re-expansion of dense, molecular regions regions on timescales shorter than $\tau_{\mathrm{diss}}$.

Since turbulent mixing seems to play an essential role for the formation of H2, the determination of chemical abundances in a post-processing step can be inaccurate, unless the full time history of the fluid element is known (e.g. by the usage of tracer particles). A simple approach based on the current state of the cell, i.e. using $\tau_{\mathrm{form}}$ as an approximation (e.g. Padoan et al., 2016a), will generally underestimate the amount of H2. Moreover, this also directly affects abundance of CO: first, H2 contributes to the shielding of radiation which dissociates CO (van Dishoeck & Black, 1988; Glover et al., 2010). Second, H2 is involved in initiating a number of chemical pathways that lead to the formation of CO. Hence, underestimating H2 will also result in a underestimation of the CO abundance. Our findings also demonstrate that for analytical estimates of atomic and molecular hydrogen column densities (e.g. Krumholz et al., 2009; Bialy et al., 2017) turbulent mixing should be taken into account.

3.5 Velocity dispersion

Next, we consider the velocity dispersion $\sigma$ of the two MCs (bottom row of Fig. 5). We define $\sigma$ with respect to their centre-of-mass (COM) velocity, $\vec{v}_{\mathrm{COM}}$,

[TABLE]

where d$m$ and $\vec{v}$ are the mass and velocity of an individual cell, respectively, and we sum over all cells identified with the corresponding definition. Note that $\sigma$ contains contributions from all the individual fragments (see Fig. 2) as well as their relative motions with respect to the COM. A more detailed investigation of what dominates the value of $\sigma$ (i.e. turbulence, interclump motions, infall, or outflow motions) is postponed to a subsequent paper.

The velocity dispersion in both clouds reaches final values around 4 – 8 km s*-1* in good agreement with observations (e.g. Solomon et al., 1987; Elmegreen & Falgarone, 1996; Heyer et al., 2001; Roman-Duval et al., 2010; Miville-Deschênes et al., 2017). With typical gas temperatures around 10 K, i.e. a sound speed $c_{\mathrm{s}}$ around 0.2 km s*-1*, this corresponds to highly supersonic motions with Mach numbers above 10. This strengthens our statement in Section 3.4 that turbulent mixing is the main reason for the presence of H2 in regions with $n$ $\simeq$ 30 cm*-3*.

For MC1, $\sigma$ is steadily increasing with time. This could be attributed to two effects: First, turbulent motions on small scales become resolved with increasing resolution. However, since only a small fraction of energy is stored on smaller scales (Kolmogorov, 1941), we rather attribute the increase of $\sigma$ for MC1 to a continuous gravitational collapse (see Fig. 2). The resulting infall velocities might be converted to turbulent motions due to the complex structure of the cloud (Klessen & Hennebelle, 2010; Goldbaum et al., 2011; Matzner & Jumper, 2015), and thus counteract the decay of turbulent, kinetic energy, which is believed to happen within one crossing time (Stone et al., 1998; Mac Low et al., 1998; Mac Low & Klessen, 2004; Elmegreen & Scalo, 2004).

This picture is also supported by the saturation of $\sigma$ for MC2 just above 4 km s*-1* after $\sim$ 3 Myr: MC2 has fragmented into a number of spatially separated, small objects (see Fig. 2), which might be less gravitationally unstable. We emphasize, however, that each of the fragments is collapsing itself, which is reflected by the increase in column density and the continuously increasing value of $M_{\mathrm{300}}$, $M_{\mathrm{H_{2}}}$, and $M_{\mathrm{CO}}$ (top right panel of Fig. 5). On the other hand, MC1 seems to collapse towards a common centre (with higher infall velocities), which is reflected by the larger amount of dense gas (Fig. 5). Furthermore, we find that at later times, for both MCs, $\sigma$ increases slightly with increasing density threshold. This further corroborates the picture that gravity/infall motions contribute to the overall velocity dispersion, in particular the observed increase. This interpretation is in line with a number of recent works on the origin of $\sigma$ in MCs (e.g. Heyer et al., 2009; Ballesteros-Paredes et al., 2011; Kritsuk et al., 2013; Ibáñez-Mejía et al., 2016, 2017), although for the initial level of turbulence ($\sigma$ = 2 – 3 km s*-1*, Fig. 5) most likely also SN feedback is important.

Finally, we note that the velocity dispersion derived for the various cloud criteria differs only little, which we attribute to the fact that it is a bulk property (i.e. its magnitude is independent of the size of the system).

4 Convergence study

In recent years a number of authors have studied the formation of MCs in galactic environments. However, little work has been done to investigate which resolution is required to accurately model the chemical, dynamical, and structural properties of the forming MCs. Typical maximum resolutions of a number of galactic disc simulations modelling the formation of MCs are of the order of a few pc or a few M*☉* in case of smoothed-particle hydrodynamics (SPH) (see Table 4). Here we present a comprehensive study to identify resolution requirements for modelling the formation of MCs in galactic discs.

In the following, we test the convergence of our results for increasing spatial resolution (see section 4.2) as well as different refinement times $\tau$ (see section 4.3). All runs are listed in Table 3.

4.1 Comparison with other work

In order to compare our work with SPH simulations, we convert the spatial resolution to a mass resolution using the Jeans length and mass. In a grid code the Jeans length $L_{\mathrm{J}}$ typically has to be resolved with about 10 grid cells (Truelove et al., 1997; Federrath et al., 2011), i.e.

[TABLE]

while in an SPH simulation the Jeans mass $M_{\mathrm{J}}$ has to be resolved with a certain number $N_{\mathrm{neigh}}$ of particles with mass $m_{\mathrm{SPH}}$, i.e

[TABLE]

Next, Eq. 4 can be used to eliminate $\rho$ on the right-hand side of Eq. 5, which allows us to solve for $m_{\mathrm{SPH}}$:

[TABLE]

Here, we adopt $N_{\mathrm{neigh}}=100$ typical for 3D simulations (Bate & Burkert, 1997) and consider isothermal ($\gamma$ = 1), molecular gas with $T$ = 10 K and $\mu$ = 2.3 to obtain a simple equation to convert between the effective spatial resolution of a grid code and the corresponding particle mass in an SPH code,

[TABLE]

In Table 4 we list a number of recent works, which study the formation of MCs in their galactic environment. We group the papers by their spatial/mass resolution to compare them with the different spatial resolutions used in this work. Most of the works listed neither reach our highest spatial resolution nor include a chemical network. Our work is thus the first which combines high spatial resolution as well as an on-the-fly modelling of molecule formation in an extensive resolution study.

4.2 Spatial resolution study

Our fiducial runs, which we discussed so far, are MC1_tau-1.5_dx-0.12 and MC2_tau-1.5_dx-0.12 with a maximum spatial resolution of d$x=0.12$ pc and a refinement time of $\tau=1.5$ Myr.

We repeat the simulations for both MC1 and MC2 using a different maximum refinement level. We perform runs with a lower effective spatial resolution of d$x$ = 0.24, 0.5, 1.0, and 3.9 pc, i.e. we stop the refinement procedure after $\tau=$1.25, 1.0, 0.5, and 0.0 Myr. These additional runs are called MC1_tau-1.25_dx-0.24, MC1_tau-1.0_dx-0.5, MC1_tau-0.5_dx-1.0, and MC1_dx-3.9, as well as MC2_tau-1.25_dx-0.24, MC2_tau-1.0_dx-0.5, MC2_tau-0.5_dx-1.0, and MC2_dx-3.9. In addition, we carry out two runs with a higher effective spatial resolution of d$x$ = 0.06 pc with $\tau=1.65$ Myr. These runs are called MC1_tau-1.65_dx-0.06 and MC2_tau-1.65_dx-0.06333Since these two simulations become numerically unstable, we were not able to run them for more than 2.7 Myr.. All runs are listed in Table 3.

4.2.1 Cloud structure, mass, and velocity dispersion

In Fig. 9 we show the total gas column density at $t$ = $t_{0}$ + 5 Myr for runs MC1_tau-0.5_dx-1.0, MC1_tau-1.0_dx-0.5, MC1_tau-1.25_dx-0.24, and MC1_tau-1.5_dx-0.12 (from left to right). While for the runs with d$x$ = 0.12 and 0.24 pc, the distribution looks relatively similar with a prominent filamentary structure, the filamentary structure seems to be only marginally resolved for d$x=$ 0.5 pc and the structure is rather clumpy for d$x$ = 1 pc. A similar change in cloud structure can be seen for the runs of MC2 (not shown here).

In Fig. 10 we show the time evolution of the masses (top panels) and velocity dispersions (bottom panels) of the two clouds as defined via the three number density thresholds $n_{\rm thres}=$ 30, 100, and 300 cm*-3*. We show all runs with different spatial resolution, ranging from d$x$ = 3.9 pc (base grid resolution) to d$x$ = 0.06 pc. As expected, higher densities are increasingly well resolved with higher spatial resolution. Therefore $M_{30}$ is already well converged at d$x\leq$0.5 pc, while $M_{100}$ requires $x\leq$ 0.24 pc, and $M_{300}$ requires d$x\leq$ 0.12 pc to meet the same time evolution. For lower resolutions, the masses are systematically underestimated. For example, for $M_{300}$, the cloud masses in the simulations with d$x=$ 0.24 pc are $\sim$ 5 – 10% below the results of the runs with d$x\leq$ 0.12 pc. For low resolutions, d$x=$ 1 pc and d$x=$ 3.9 pc, $M_{300}$ is reduced by 30% and 60%, respectively. Since gas with number densities greater than 100 cm*-3* can quickly become molecular, the insufficient resolution of high density gas is accompanied by an artificially reduced molecular gas fraction in the simulations (see 4.2.2).

The velocity dispersion (bottom row of Fig. 10) systematically increases with $n_{\rm thres}$ and with resolution. However, the differences in velocity dispersion between runs with different spatial resolution hardly depend on $n_{\rm thres}$. Again we find that runs with d$x=$ 1 pc and d$x=$ 3.9 pc highly underestimate $\sigma$. For d$x=$ 0.24 pc, the velocity dispersion is a few 0.1 km s*-1* ($\sim$ 10%) smaller than for d$x=$ 0.12 pc. Since the cloud dynamics at late times is dictated by self-gravity and turbulence, we do not expect $\sigma$ to be fully converging at d$x\leq$ 0.12 pc since small scales (small turbulent eddies and/or small, collapsing structures) still contribute to the kinetic energy content of the cloud.

4.2.2 Chemical properties

In Fig. 11 we plot the mass fractions of H2 and CO, $f_{{\rm H}_{2}}$ (left column) and $f_{{\rm CO}}$ (middle column), in both zoom-in regions, with

[TABLE]

Here, $M_{\mathrm{H,tot}}$ and $M_{\mathrm{C,tot}}$ are the total mass of hydrogen and carbon in the considered volume, $M_{\mathrm{H_{2},tot}}$ and $M_{\mathrm{CO,tot}}$ the mass of all H2 and CO molecules, respectively, and the factor $\frac{12}{28}$ corrects for the mass of oxygen. We refer to Fig. 4 for the time evolution of the total masses $M_{\mathrm{H_{2},tot}}$ and $M_{\mathrm{CO,tot}}$ in the zoom-in regions.

Interestingly, although MC1 and MC2 show a different evolution in the total masses (Figs. 3 and 5), the mass fractions are remarkably similar. In general, $f_{\mathrm{CO}}$ is always below $f_{\mathrm{H_{2}}}$ because CO forms at higher visual extinction ($A_{\mathrm{V}}$ $\gtrsim$ 1) than H2 ($A_{\mathrm{V}}$ $\gtrsim$ 0.3) (e.g. Röllig et al., 2007; Glover et al., 2010). In agreement with the convergence of the mass evolution discussed in Fig. 10, we find that $f_{\mathrm{CO}}$ and $f_{\mathrm{H_{2}}}$ are grossly underestimated for low resolution runs and only start to converge for d$x\leq$ 0.24 pc (the difference in $f_{\mathrm{CO}}$ and $f_{\mathrm{H_{2}}}$ between runs with d$x$ = 0.12 and 0.06 pc is smaller than 5%). Although the low resolution runs would also form significant amounts of molecular gas towards the end of the runs, the corresponding timescale would be significantly longer. This resolution dependence also affects the transition to H2 dominated gas (Fig. 7, blue line). With decreasing resolution, the value of 30 – 50 M*☉* pc*-2* for d$x$ = 0.12 pc, where $\Sigma_{\mathrm{H_{2}}}/\Sigma_{\mathrm{HI}}$ $>$ 1 , increases up to $\sim$ 100 M*☉* pc*-2* for d$x$ = 1 pc. For d$x$ = 3.9 pc, up to this time ($t$ = $t_{0}$ + 5 Myr) the gas does not become H2 dominated at all. In addition, the molecular gas sits in blobby configurations at low resolution, rather than in thin filamentary structures (see Fig. 9).

Furthermore, we plot the ratio of the total number of CO and H2 molecules within the zoom-in region, $N_{\rm CO}/N_{{\rm H}_{2}}$, in the right column of Fig. 11. The molecular clouds mostly form within the first $\sim$ 2 Myr of the zoom-in simulations and thus, $N_{\rm CO}/N_{{\rm H}_{2}}$ is rapidly increasing in this period. From this one could naively assume a highly variable $X_{\rm CO}$ during cloud formation (see Section 5). At later times the ratio saturates and we find that $N_{\rm CO}/N_{{\rm H}_{2}}\sim 1.8-1.9\times 10^{-4}$ at $t=t_{0}+5$ Myr for both molecular clouds and for all but the lowest resolution run. The fractional abundance of carbon with respect to hydrogen nuclei used in the simulations is 1.4 $\times$ 10*-4*. Therefore the final value of $N_{\rm CO}/N_{{\rm H}_{2}}$ is somewhat below the ratio of 2.8 $\times$ 10*-4* for fully molecular gas. Our resolution study indicates that d$x\lesssim$ 0.24 pc should be used to determine $X_{\rm CO}$ in simulations of molecular cloud formation.

The initial increase of $N_{\rm CO}/N_{{\rm H}_{2}}$ is a consequence of the rapid formation of CO in this time period. The rate-limiting reaction is the formation of hydrocarbon radicals like CH or CH2 out of the (already present) reservoir of H2 (left panel of Fig. 11) with a rate coefficient of 5 $\times$ 10*-16* cm3 s*-1* (Nelson & Langer, 1997; Glover et al., 2010). Assuming a fractional abundance of H2 of 10%, the typical timescale for CO formation is thus

[TABLE]

which is slightly faster than the canonical formation time of H2, $\tau_{\rm H_{2}}\approx 1\;{\rm Gyr}\left(n/1\;{\rm cm}^{-3}\right)^{-1}$, thus explaining the initial increase of $N_{\rm CO}/N_{{\rm H}_{2}}$.

The formation of H2 and CO occurs at high column densities, which shield the forming molecules from the ambient interstellar radiation field (ISRF). A different way to analyse at which resolution we expect to converge with respect to the H2 and CO mass fractions is to look at the mass-weighted probability-density-function (PDF) of visual extinction $A_{\mathrm{V}}$ vs. gas density. Here, $A_{\mathrm{V}}$ is the effective visual extinction of radiation in each cell in the simulation domain obtained via the TreeCol algorithm (Clark et al., 2012a). In Fig. 12, we plot the resulting PDFs of MC1 at $t$ = $t_{0}$ + 5 Myr for different resolutions (from left to right: d$x$ = 1.0, 0.5, 0.24, and 0.12 pc). We expect most hydrogen to be in the form of H2 for $A_{\mathrm{V}}\gtrsim 0.3$ (see e.g. Röllig et al., 2007, for benchmark tests for photodissociation regions, but also Glover et al. 2010) and most carbon to be in the form of CO for $A_{\mathrm{V}}\gtrsim 1$. The distribution of $A_{\mathrm{V}}$ is broad due to the different amount of shielding of different regions in the cloud, i.e. at a given density cells near the cloud surface are found at low $A_{\mathrm{V}}$ and cells within the deeply embedded interior are found at high $A_{\mathrm{V}}$. Most of the volume is occupied by lower density, low $A_{\mathrm{V}}$ gas. An $A_{\mathrm{V}}$ of 0.3 is reached for number densities between $\sim$ 1 – 100 cm*-3*, while an $A_{\mathrm{V}}$ of 1.0 typically occurs at number densities of $\sim$ 300 cm*-3*, which is barely resolved in the lower resolution simulations.

4.2.3 Physical interpretation

Combining the results from the Sections 4.2.1 and 4.2.2, we conclude that an effective spatial resolution of d$x$ $\lesssim 0.12$ pc is required to obtain a reasonably well converged result for structural, dynamical and chemical properties.

For CO this agrees well with the fact that the mass in gas above $n_{\rm thres}=300\;{\rm cm}^{-3}$ is well converged at d$x\lesssim$ 0.12 pc. For H2, we could expect convergence at d$x\lesssim$ 0.5 pc (from the mass evolution and the $A_{\mathrm{V}}$ distribution). However, the evolution of $f_{{\rm H}_{2}}$ indicates that d$x\lesssim$ 0.12 pc is also required to properly resolve the formation of molecular hydrogen. This is due to the fact that H2 forms quickly in dense regions and is then efficiently mixed into the lower density environment, an effect which is not well resolved at lower resolution.

We note that the resolution requirement of $\sim$ 0.1 pc also agrees with the sonic scale of interstellar turbulence (e.g. Goodman et al., 1998; Vázquez-Semadeni et al., 2003; Federrath et al., 2010, and references therein). On these scales velocity-coherent structures are believed to form, which fits with our findings that the structural properties converge around this resolution (see Fig. 9.)

This finding has implications for a number of recent studies (see Table 4), in which a lower spatial resolution is used. Our results suggest that there are limitations for the accuracy of the quantities derived from these simulation data. We also speculate that there are similar limitations for associated synthetic observations produced for comparison with actual observations, e.g. the linewidth-size relation.

4.3 Impact of the refinement time

So far we have only considered simulations with a refinement time $\tau=1.5$ Myr at d$x=0.12$ pc. However, one might ask why an instantaneous refinement up to the maximum refinement level, i.e. $\tau$ = 0 should not be allowed. We discuss the results of simulations with $\tau$ = 0 in section 4.3.1 and discuss why it does not work well. Furthermore, we could decrease/increase $\tau$ without going to the extreme case of instantaneous refinement. The results for runs with $\tau=$ 1.0, 2.25, and 4.5 Myr are shown in section 4.3.2, where the impact of $\tau$ on the cloud mass evolution, molecular mass fraction evolution, velocity dispersion, and cloud substructure is discussed.

4.3.1 Instantaneous refinement $\tau=0$

We consider both molecular clouds with an instantaneous refinement, i.e. the runs MC1_tau-0_dx-0.12 and MC2_tau-0_dx-0.12 (see Table 3). In Fig. 13 we show the column density of run MC1_tau-0_dx-0.12 at different times. Soon after we start the zoom-in calculation (up to 1 – 2 Myr), the instantaneous refinement causes strong grid artefacts. The ”blocky” structures are formed if a coarse (parent) cell is refined multiple times within just a few time steps and thus, all of its sub-cells still contain almost equal hydrodynamical quantities since the small-scale velocity and density fluctuations had no time to developed. This problem cannot be easily solved even by a first- or second-order prolongation scheme from coarse to fine cells since the conservation of mass and momentum has to be guaranteed. As a consequence, at the boundaries of some parent cells, gas will be pushed together and compressed, thus leading to artificial filamentary structures which were unresolved on the coarse grid.

At later times, the grid artefacts are washed out, but the filamentary structures they seeded can survive and develop into dense filaments that cannot be distinguished from the ”real” filamentary structures which form in turbulent molecular clouds. Comparing the $\tau=0$ Myr and the fiducial $\tau$ = 1.5 Myr runs (compare with Fig. 2, in particular the southern part of the cloud at $t\geq$ 3 Myr), it can be seen that the cloud substructure in run MC1_tau-0_dx-0.12 is different. Similar results are also found for MC2.

We see such grid artefacts for any $\tau<$ 1 Myr. Therefore, we use $\tau=$ 1.0 Myr, which corresponds to about 200 time steps on each refinement level, as a lower limit and do not consider runs with $\tau<1$ Myr in the following section.

4.3.2 Influence of different $\tau>0$

As discussed before, a refinement time $\tau$ $<$ 1 Myr results in undesired grid artefacts, which imposes a lower limit on $\tau$. In order to show the impact of different $\tau$, we repeat the runs for both molecular clouds with d$x=$ 0.12 pc and different $\tau$. Consequently, the time (number of time steps) which we spend on each refinement level is adapted (see Table 3). In particular, we use one refinement time that is shorter than our fiducial value of 1.5 Myr, that is $\tau=$ 1.0 Myr, and two larger values, $\tau=$ 2.25 and 4.5 Myr. The runs are called MC1_tau-1.0_dx-0.12, MC1_tau-2.25_dx-0.12, and MC1_tau-4.5_dx-0.12 for MC1, and MC2_tau-1.0_dx-0.12, MC2_tau-2.25_dx-0.12, and MC2_tau-4.5_dx-0.12 for MC2 (see Table 3).

In Fig. 14 we show the column densities of MC1 at $t=t_{0}+1.5$ Myr (top row) and $t=t_{0}+5$ Myr (bottom row) for the four different $\tau$ used (from left to right). At $t=t_{0}+1.5$ Myr the runs with $\tau>$ 1.5 Myr have not yet reached their highest resolution. This affects the cloud structure at later times: As $\tau$ increases, less sub-structure is visible at $t=t_{0}+5$ Myr. While the differences between runs MC1_tau-1.0_dx-0.12 and MC1_tau-1.5_dx-0.12 are small, we see significantly less fragmentation in MC1_tau-4.5_dx-0.12. Again, the results for MC2 are qualitatively similar.

In Fig.  15 we compare the total mass fractions of H2 and CO within the zoom-in region as well as the cloud masses and velocity dispersions obtained for three cloud definitions: a density threshold of $n_{\rm thres}=100\;{\rm cm}^{-3}$, and the H2 and CO threshold. The runs with long refinement times (in particular $\tau=$ 4.5 Myr) are not able to properly capture the formation of both CO and H2 at early times. In particular for H2 (top panels) the mass fractions remain lower than for shorter $\tau$ throughout the simulation. We attribute this to two reasons: first, the turbulent formation mechanism of H2 (see Section 3.4 but also Glover & Mac Low, 2007b) is more efficient for higher spatial resolution and thus underestimated for long $\tau$. Second, due to the lower resolution for longer $\tau$ the mass at high densities ($M_{100}$, third row of Fig. 15), in which H2 is formed, is lower. For this reason, also $M_{\mathrm{H_{2}}}$ and $M_{\mathrm{CO}}$ grow more slowly for longer refinement times, in agreement with the H2 and CO mass fractions. For all shorter values of $\tau\leq$ 2.25 Myr, we see only small differences within the first 2-3 Myr of evolution. In particular for $\tau$ = 1.0 Myr and 1.5 Myr the changes are minor.

Once each simulation reaches the highest refinement level, both masses and velocity dispersions seem to converge for different $\tau$, however, there are several problems which become apparent if the runs are inspected more closely.

First, the velocity dispersions do not evolve smoothly for $\tau$ = 4.5 Myr, but jump to a higher value as soon as a higher refinement level is introduced. This clearly shows that the dynamics of the cloud is not well captured if $\tau$ is too long. Second, the cloud substructure changes irrevocably if $\tau=$ 4.5 Myr, which we show for MC1444This particularly bad example occurred in a simulation of MC1 including SN explosions close to the zoom-in region between $t_{0}$ and $t_{0}$ + 5 Myr. The effect of nearby SNe will be discussed in detail in a subsequent paper. However, we note that the changes between the simulations presented here and the simulations with nearby SNe are overall minor in this particular run. in Fig. 16. A too long refinement time leads to a rotationally supported structure of size $\sim$ 5 – 10 pc in the centre of the cloud (top panels). This disc-like cloud is broken up into smaller filaments and blobs if the refinement time is shorter (see $\tau$ = 1.5 Myr in the bottom panels). As turbulent motions in the molecular cloud locally create converging flows which carry some amount of net angular momentum with respect to the collision point, rotating structures on the grid scale can be formed when the flows converge under the influence of self-gravity. These structures might be unphysical if smaller-scale, turbulent motions are unresolved. Since for the run with $\tau$ = 4.5 Myr the resolution is lower than for $\tau$ = 1.5 Myr at the moment when the flows collide ($t\simeq$ $t_{0}$ + 3.5 Myr), this results in a rather large, rotating disc which does not disappear, even after the resolution is increased further. A numerical artefact like this could lead to the frequent appearance of rotating molecular clouds in larger-scale galaxy-scale simulations (e.g. Renaud et al., 2013; Bournaud et al., 2014), while observations of such large-scale rotating molecular cloud structures are very rare (e.g. Li et al., 2017, for one example).

4.3.3 Physical interpretation

Overall, we have shown that there is a lower limit ($\sim$ 1 Myr) as well as an upper limit ($\sim$ 1.5 Myr) for $\tau$. The lower limit on $\tau$ is set by a minimum number of time steps ($\gtrsim$ 200) which should be spent on each newly introduced refinement level, and is thus most likely applicable to any kind of grid code, although the exact numbers may vary.

The upper limit of $\tau$ can be explained by means of the free-fall time of the clouds. At $t$ = $t_{0}$, there is gas with number densities between 100 and 300 cm*-3* but no gas with $n$ $>$ 300 cm*-3* in both zoom-in regions555In the Appendix we discuss the effect when starting the entire refinement procedure earlier and show that overall we find a good agreement between the runs with different zoom-in starting times.. The corresponding free-fall times are $\tau_{\rm ff}$ = 3.40 and 1.96 Myr, respectively. Moreover, since the clouds seem to assemble within about 2 Myr (see upper row of Fig. 5), $\tau$ $<$ 2 Myr seems to be advisable. This is also indicated by the moderate differences occurring for the runs with $\tau$ = 2.25 Myr, the runs with $\tau$ = 4.5 Myr $>\tau_{\rm ff}$ lead to a clear delayed collapse of gas into dense regions (see Fig. 15). We emphasize that $\tau_{\rm ff}$ – and thus the upper limit of $\tau$ – naturally depends on the mean density of the simulated object.

4.4 Impact of magnetic fields

In the simulations presented here we have neglected magnetic fields, which are able to significantly delay the collapse time of dense objects by exerting an additional magnetic pressure, which counteracts gravity (e.g. Heitsch et al., 2001; Vázquez-Semadeni et al., 2011; Collins et al., 2012; Federrath & Klessen, 2012; Walch et al., 2015; Girichidis et al., 2016; Mocz et al., 2017). This can alleviate the tight constraints on the refinement time $\tau$ discussed in Section 4.3: Since in the presence of magnetic fields the actual collapse time can become longer – in particular longer than $\tau_{\mathrm{ff}}$, which we take as the rough upper limit for $\tau$ – the upper limit of $\tau$ could be increased. Furthermore, due to the slower collapse, large and well-shielded structures appear later, which in turn delays the formation of H2 and CO (see Fig. 10 in Walch et al., 2015). Hence, the initial formation time (see Fig. 4), which has to be resolved, becomes longer, which makes it easier to accurately model the chemical evolution. We note that we plan to study the impact of magnetic fields in a follow-up paper.

5 $X_{\mathrm{CO}}$-factor

As demonstrated in Fig. 11, the ratio of CO to H2 molecules increases for the first 2 Myr and then remains roughly constant. This might have some implications for the conversion factor $X_{\mathrm{CO}}$, which is usually found to have values around 2 $\times$ 1020 cm*-2* (K km s*-1*)-1 (Bolatto et al., 2013). In order to test this, we produce 12CO, $J$ = 1 – 0 line emission maps of our two MCs at $t$ = $t_{0}$ + 1.0, 1.5, 2.0, and 5.0 Myr with RADMC-3D (Dullemond, 2012) using the Large Velocity Gradient (LVG) method. The molecular data, e.g. the Einstein coefficients, are taken from the Leiden Atomic and Molecular data base (Schöier et al., 2005). We use 201 velocity channels with a width of 100 m s*-1* to guarantee that all emission is captured properly. We produced emission maps for line-of-sight (LOS) directions along the $x$, $y$, and $z$ axis. We note that we do not include any observational effects like noise or resolution limitations.

In Fig. 17 we show the integrated intensity as well as the $X_{\mathrm{CO}}$ maps of MC1 at $t$ = $t_{0}$ + 2.0 Myr for the LOS along the $x$-axis, i.e. the same direction as shown in Fig. 2. Overall, the $X_{\mathrm{CO}}$ factor shows a strong spatial variation of more than one order of magnitude in this region, which also holds for all other maps produced. Similar results were found in numerical (Clark & Glover, 2015; Glover & Clark, 2016; Szűcs et al., 2016) and observational works (e.g. Lee et al., 2014).

In order to quantify the time-dependence of $X_{\mathrm{CO}}$ for both clouds and all three LOS, we calculate the mean conversion factor $\langle X_{\mathrm{CO}}\rangle$ for all pixels in the maps. The values of $\langle X_{\mathrm{CO}}\rangle$ show only a weak time dependence (see Fig. 18) and scatter around 1 – 4 $\times$ 1020 cm*-2* (K km s*-1*)-1, which is in good agreement with observations (see e.g. Bolatto et al., 2013, for a recent review). From the CO-to-H2 ratio shown in the right panel of Fig. 11 one might expect a decreasing $\langle X_{\mathrm{CO}}\rangle$ factor with time. However, the total CO intensity increases slower than the number of CO molecules and is roughly balanced by the increasing amount of H2, which results in an almost constant $\langle X_{\mathrm{CO}}\rangle$ factor. Overall, our results thus indicate that the $X_{\mathrm{CO}}$ factor seems to be rather robust and independent of the evolutionary stage of the considered MC.

We attribute the slower increase of the emission of CO with time to the fact that it becomes optically thick in the densest regions. The optical depth of CO can be estimated via

[TABLE]

where $h$ and $k_{\mathrm{B}}$ are Planck and the Boltzmann constant, $A_{\mathrm{10}}$ and $\nu_{\mathrm{10}}$ the Einstein coefficient and frequency of the considered transition, $\Delta v$ the line-width, and $T_{\mathrm{ex}}$ excitation temperature. This can be rewritten as

[TABLE]

Hence, for $\Delta v$ $\simeq$ 4 km s*-1* (Fig. 5) and $T_{\mathrm{ex}}$ $\simeq$ 15 K, CO starts to become optically thick at $A_{\mathrm{V}}$ $\sim$ 1, i.e. there is a significant amount of optically thick gas in the clouds (see Fig. 12).

There appears to be a slight trend of an increase of $\langle X_{\mathrm{CO}}\rangle$ with time, which is in good agreement with the findings of Glover & Clark (2016, see their Fig. 7), but in contrast to the results reported by Richings & Schaye (2016a, b), who find a decreasing $\langle X_{\mathrm{CO}}\rangle$ factor over time. As we argue in Section 4.2, the result of Richings & Schaye (2016a, b) might be an artefact of the mass resolution in their SPH simulations, which is orders of magnitude lower than the corresponding spatial resolution used in our simulations, thus leading to a significant underestimation of CO.

6 Conclusions

We present simulations of the self-consistent formation of two molecular clouds (MCs) in a stratified galactic disc in the context of the SILCC project (Walch et al., 2015; Girichidis et al., 2016). The simulations include a chemical network to model the formation of H2 and CO. Furthermore, we make use of a zoom-in technique, which allows us to follow the large-scale, galactic environment and, at the same time, resolve the formation of filamentary MCs down to scales of $\sim$ 0.1 pc.

The simulated MCs typically grow from the low-density environment by accreting mass at a rate of a few times $10^{-3}$ to 10*-2* M*☉* yr*-1* for a couple of Myr and reach masses between 3 – 7 $\times$ 104 M*☉. We tentatively argue that gravity contributes a significant part to the velocity dispersion, which reaches 4 – 8 km s-1*. The transition to molecular hydrogen dominated gas occurs around 30 – 50 M*☉* pc*-2*. At later stages, the dense gas ($n\gtrsim 100$ cm*-3*) is almost fully molecular with H2 fractions around 0.9. We show that different cloud definitions (such as based on density thresholds, H2 or CO mass fraction) result in significantly different cloud properties. Gas which contains a significant fraction of H2 ($>$ 50% H2) grows differently from the dense and CO dominated gas and cannot be matched by a single density threshold. While CO follows the evolution of the denser gas with $n\geq$ 300 cm*-3*, H2 is found in gas with lower and lower number density (even in gas with $n\lesssim 30$ cm*-3*). From estimating the relative time scales, it seems that molecular hydrogen is efficiently mixed by turbulent motions into the lower density gas rather than being formed in-situ (see also Glover & Mac Low, 2007b; Valdivia et al., 2016). We speculate that calculating the chemical abundances from simulations without any chemical network in a post-processing step is significantly complicated by this mixing process.

The CO-to-H2 ratio changes as a function of time, in particular at early times when the mass accretion rate onto the cloud is high, later on it settles to an approximately constant value of $\sim 1.8\times 10^{-4}$. Since the CO(1-0) line, however, becomes quickly optically thick, the XCO conversion factor remains rather constant over time with values around 1 – 4 $\times$ 1020 cm*-2* (K km s*-1*)-1 and thus seems to be independent of the evolutionary state of a MC.

From the work presented here, we identify four requirements for numerically converged zoom-in simulations of molecular cloud formation:

An effective resolution of $\lesssim 0.1$ pc is required to obtain converged mass fractions of molecular hydrogen and CO as well as realistic morphological and dynamical properties. 2. 2.

The adaptive mesh refinement (AMR) should proceed stepwise, thus allowing the hydrodynamic quantities to relax on each refinement level. Therefore, in any grid code $\gtrsim$ 200 time steps should be spent on a level before a next higher level is introduced. Otherwise grid artefacts at the boundaries of coarse blocks might occur. 3. 3.

On the other hand, for simulations of self-gravitating objects, the total refinement time (after which the resolution reaches its maximum) should be shorter than its free-fall time to follow the actual collapse of the dense structures (here $<$ 2 Myr). This ensures that realistic, small-scale substructures can form and avoids the formation of numerical artefacts like rotationally supported objects. 4. 4.

The total refinement time should also be adjusted to the formation time of chemical species in turbulent environments in order to prevent the delay of H2 and CO formation.

If fulfilled, the forming molecular clouds have converged chemical compositions, mass distributions, and substructure, and grid artefacts are minimised.

Since the paper serves as a proof-of-concept, it does by no means cover all aspects of the simulations. For future work we plan to investigate in more detail important parameters such as the origin of the velocity dispersion, the filamentary structure, the chemical evolution, the virial state as well as column density PDFs of the clouds. Moreover, since we have detailed information about the chemical composition of the gas as well as the dust temperature, we plan to produce synthetic line and continuum emission maps for a detailed comparison with current observational data.

Acknowledgements

The authors like to thank the anonymous referee for the comments which helped to significantly improve the paper. DS and SW acknowledge the support of the Bonn-Cologne Graduate School, which is funded through the German Excellence Initiative. SW further acknowledges support via the ERC starting grant No. 679852 ”RADFEEDBACK”. PG, SW, TN, SCOG, RSK, and TP acknowledge support from the DFG Priority Program 1573 “Physics of the Interstellar Medium”. PG acknowledges funding from the European Research Council under ERC-CoG grant CRAGSMAN-646955. TN acknowledges support from the DFG cluster of excellence ”Origin and Structure of the Universe”. RW acknowledges support by the Czech Science Foundation grant 15-06012S and by the project RVO:67985815 of the Academy of Sciences of the Czech Republic. RSK and SCOG thank the DFG for funding via the SFB 881 “The MilkyWay System” (subprojects B1, B2, and B8). RSK furthermore acknowledges support from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013) via the ERC Advanced Grant STARLIGHT (project number 339177). PCC acknowledges support from the Science and Technology Facilities Coun- cil (under grant ST/N00706/1) and the European Com- munitys Horizon 2020 Programme H2020-COMPET-2015, through the StarFormMapper project (number 687528). The FLASH code used in this work was partly developed by the Flash Center for Computational Science at the University of Chicago. The simulations were performed at SuperMUC at the Leibniz-Rechenzentrum Garching within the projects pr94du and pr45si.

Appendix A Influence of the starting time

Here we present results of simulations of MC1 and MC2 when starting to zoom-in 0.9 Myr earlier than in our fiducial runs MC1_tau-1.5_dx-0.12 and MC2_tau-1.5_dx-0.12, i.e. at $t$ = $t_{0}$ - 0.9 Myr. In Fig. 19 we compare the evolution of the column density for run MC1 when zooming in earlier (bottom row) with the one obtained from the fiducial run (top row). In Fig. 20 we show the same for MC2. Qualitatively, the exact time at which we start to zoom-in seems to have only little effect for both clouds.

In order to test the dependence on the starting time more quantitatively, in Fig. 21 we compare the time evolution of the mass, mean density, and velocity dispersion obtained via the different criteria used in this paper for both MCs.

As expected, there are differences in the various quantities in particular in the beginning ($t\lesssim$ 2 Myr), i.e. during the refinement time. All considered quantities grow at earlier times when starting to zoom-in earlier, which reflects that in this case the gravitational collapse sets in earlier. However, towards later times the various values approach each other with typical deviations of a few 10% or less. The only exception is $\bar{n}_{\mathrm{H_{2}}}$, which is somewhat lower for the earlier zoom-in starting time. This is a consequence of the turbulent mixing of H2 into the low density regions, which increases the occupied volume (compare Section 3.4) and therefore apparently decreases the mean density of the H2 dominated gas. The total mass $M_{\mathrm{H_{2}}}$, however, is converge towards the end. All general trends of MC1 and MC2 are recovered, independent of the starting time.

Finally, we consider the dependence of the H2 and CO mass fractions on the starting time in Fig. 22. We find a similar behaviour to the mass evolution discussed above: Below $t_{0}$ + 2 Myr, i.e. during the zoom-in procedure, the molecular abundances are smaller for the fiducial runs where we start to zoom-in later. However, as time evolves, the differences become significantly smaller reaching relative deviations of a few percent towards the end of the runs.

In summary, for a complex, turbulent object like a molecular cloud, a perfect match between simulations with different zoom-in starting times cannot be expected. However, the overall good agreement (to within a few percent) of the structural, dynamical, and chemical properties makes us confident about the results presented in this paper.