STAR FORMATION IN DISK GALAXIES. I. FORMATION AND EVOLUTION OF GIANT MOLECULAR CLOUDS VIA GRAVITATIONAL INSTABILITY AND CLOUD COLLISIONS

As we move to smaller scales (∼100 pc) we see, based mostly on surveys of CO emission in the Milky Way (Solomon et al. 1987; Dame et al. 2001; Jackson et al. 2006) that most of the dense, cold, potentially star-forming gas is organized into giant molecular clouds (GMCs) with masses ∼10⁶ M_☉ and radii ∼30 pc (with these particular values of mass and size corresponding to volume-averaged number densities of H nuclei of n_H = 260 cm^-3 for a spherical cloud). Given that most of the Galactic molecular gas is organized in GMCs and that GMCs have approximately equal mass atomic envelopes (Blitz 1990), the results of Wolfire et al. (2003) for the radial distribution of Galactic molecular and atomic gas indicate that a large fraction, ∼1/3, of the total gas mass in the Milky Way inside the solar circle is associated with these structures. Many questions concerning GMCs are still debated, including whether they are long or short-lived compared to their free-fall timescale,

$$\begin{equation} t_{\rm ff} = \left(\frac{3\pi }{32G\rho }\right)^{1/2} = 4.35 \times 10^6 \left(\frac{n_{\rm H}}{100\,{\rm cm^{-3}}}\right)^{-1/2}{\rm yr}, \end{equation} \tag{ 5 }$$

and whether they are gravitationally bound (McKee & Ostriker 2007, and references therein).

Observations of GMCs are most complete in the Milky Way interior to the solar circle (Solomon et al. 1987; Williams & McKee 1997; Heyer et al. 2008). Except for the vicinity of the Galactic center, there is a dearth of clouds in the inner ∼2 kpc. GMCs occupy a very thin vertical distribution in the Galaxy: Bronfman et al. (2000) derived a half-intensity height of z_1/2 ≃ 60 pc; Stark & Lee (2005) derived a vertical scale height of ≲35 pc.

Stark & Brand (1989) derived a one-dimensional cloud to cloud rms velocity dispersion of 7.8 ± 0.5 km s^-1 from a study of GMCs within 3 kpc of the Sun, although this estimate includes small-scale streaming motions.

Based on the ¹²CO survey data of Solomon et al. (1987), Williams & McKee (1997) derived a cloud mass function of the form

$$\begin{equation} \frac{d{\cal N}_c}{d \ln M_c} = {\cal N}_{\rm cu} \left(\frac{M_c}{M_u} \right)^{-\alpha _{c}} \end{equation} \tag{ 6 }$$

for M_c ⩽ M_u, with $d {\cal N}_c(M_c)$ being the number of clouds with masses in the range [M_c, M_c(1 + dln M_c)]. Based on clouds in the mass range of a few ×10⁵ M_☉ to a few ×10⁶ M_☉, Williams & McKee (1997) estimated that the population of GMCs inside the solar circle is described by α_c = 0.6, ${\cal N}_{\rm cu}=63$, and M_u = 6 × 10⁶ M_☉. If the observational surveys have higher degrees of incompleteness at lower masses, Williams and McKee estimated α_c = 0.85, ${\cal N}_{\rm cu}=25$, and M_u = 6 × 10⁶ M_☉. Note that these mass estimates did not include any atomic gas that might be associated with the GMCs, which could be a substantial fraction (∼1/2) of the total mass bound to the cloud (Blitz 1990).

The mean mass surface density, Σ_c, of Galactic GMCs was derived by Solomon et al. (1987) to have a median value of about 200 M_☉ pc⁻². Recently, Heyer et al. (2008) derived a median value of 42 M_☉ pc^-2, based on an LTE analysis of ¹³CO survey data. They estimate that the true values are larger because of subthermal excitation and abundance variations and are in the range 80–120 M_☉ pc^-2. Again these estimates do not include any atomic gas associated with the clouds.

Bertoldi & McKee (1992) define the virial parameter,

$$\begin{equation} \alpha _{\rm vir} \equiv \frac{5 \sigma _{c}^2 R_{c,A}}{G M_c}, \end{equation} \tag{ 7 }$$

where σ_c is the mass-averaged one-dimensional velocity dispersion of the cloud, which includes thermal and nonthermal contributions via σ_c ≡ (c²_s + σ²_nt,c)^1/2 where σ_nt,c is the one-dimensional rms velocity dispersion about the cloud's center-of-mass velocity, and R_c,A is a measure of cloud radius based on its projected area. The virial parameter is a measure of the ratio of the kinetic to gravitational energy of a cloud. For a uniform, spherical cloud, α_vir = 1 implies the total kinetic energy of the cloud is half the magnitude of the gravitational energy. Nonspherical and nonuniform density distributions typically make only modest effects, parameterized via the dimensionless factor a in the gravitational energy equation W = −(3/5)aGM²_c/R_c,A: Bertoldi & McKee (1992) estimate that cloud major to minor axis ratios of a factor of 3 only cause a to deviate from unity by ≲10%, while a power-law density profile ρ ∝ r^−3/2 results in a = 1.25. Observationally, Heyer et al. (2008) find a median value of α_vir of about unity.

Phillips (1999) considered the rotational properties of Galactic GMCs, finding a significant spread in the directions of the angular momentum vectors (as derived from radial velocity gradients and as projected on the plane of the sky), indicating that a substantial fraction of GMCs rotate in a retrograde direction with respect to Galactic rotation.

Observations of GMCs in other galaxies tend to find they have similar properties as Galactic GMCs, including their velocity dispersions and virial parameters (Bolatto et al. 2008). Rosolowsky et al. (2003) measured the rotational properties of large GMCs in M33, finding relatively small values of specific angular momentum, and a substantial population (∼40%) of clouds with rotation directions retrograde with respect to that of the galaxy. Fukui et al. (2008) studied 164 GMCs in the Large Magellanic Cloud, concluding that these clouds were also close to virial equilibrium.

The star formation activity within GMCs appears highly clustered. Most stars are born from star-forming clumps that turn into star clusters (Lada & Lada 2003), with initial radii ∼1 pc. Here, the local overall star formation efficiency is relatively high, _* ≡ M_*/M_gas ∼ 0.1–0.5, where M_* is the total mass of stars formed from the mass of gas M_gas that occupies the same volume as forming star cluster. Conversely, most of the volume of GMCs is not actively forming stars, perhaps because of magnetic field support (Crutcher 2005).

Theoretical work on the formation of GMCs has led to different schools of thought (McKee & Ostriker 2007, and references therein). "Top-down" formation mechanisms suggest that GMCs form via large-scale gravitational or magnetic disk instabilities (e.g., Kim et al. 2003; Shetty & Ostriker 2006; Glover & Mac Low 2007a, 2007b), whereas "bottom–up" processes have GMCs forming via agglomeration from inelastic collisions of clouds (e.g., Kwan 1979) or from turbulent or colliding flows (e.g., Bergin et al. 2004; Heitsch et al. 2008). It is possible that both processes may be important depending on the galactic environment (Dobbs 2008). If GMCs are relatively long-lived and contain a large fraction of the total ISM, then the formation and destruction of GMCs to and from the atomic phase may be less important than the interactions between those already existing (Tan 2000): in other words, the nature of a particular cloud may change more drastically via merging collisions with other GMCs than via flows into or out of the atomic phase.

In this paper, as a first step toward understanding galactic SFRs, we will examine the formation and evolution of GMCs in the context of the global dynamics of a flat rotational curve disk galaxy.

Numerous theories have been proposed to explain the observed galaxy-scale, Kennicutt star formation relations. One group of theories is based on the growth rate of gravitational perturbations in a disk. The timescale for perturbation growth can be expressed as t_grow ∝ (Gρ_g)^−0.5 (e.g., Larson 1988; Elmegreen 1994; Wang & Silk 1994), and so ρ_sfr ∝ ρ_g/t_grow ∝ ρ^1.5_g. Assuming a constant disk scale height, we obtain Equation (1) with α_sfr = 1.5 for local mass surface densities. However, disk-averaged quantities will depend on the radial gas distribution. We can also express t_grow ∼ σ_g/(πGΣ_g) ∼ Q/κ. Perturbation growth via swing amplification in a differentially rotating disk occurs at a similar rate (e.g., Larson 1988). By assuming that star formation self-regulates and keeps Q constant, Larson (1988) and Wang & Silk (1994) predicted Σ_sfr ∝ Σ_g/t_grow ∝ Σ_gΩ, since κ ∝ Ω, for disks with flat rotation curves. Li et al. (2006) presented isothermal smooth particle hydrodynamic (SPH) simulations of disk galaxies, from which they concluded that the rate of the nonlinear development of gravitational instability determines the local and global Kennicutt relations.

However, these theories that involve the SFR being set by the growth rate of large-scale gravitational instabilities leading to GMC formation, have difficulty in explaining why a large mass fraction of the gas is already organized into GMCs. If the GMCs are forming rapidly in ∼1 dynamical timescale the mass flux into (and out of) GMCs would be ∼100 times greater than that from GMCs into star clusters (Zuckerman & Evans 1974). To understand global SFRs, it seems more reasonable to look for processes that create the star-forming parsec-scale clumps within GMCs, rather than the processes that create the GMCs themselves.

The spatial correlation of star formation with large-scale spiral structure in some disk galaxies motivates theories for the triggering of star formation during the passage of gas through density waves. Wyse (1986) and Wyse & Silk (1989) propose a SFR law of the form

$$\begin{equation} \Sigma _{\rm sfr}\propto \Sigma _{g}^{\alpha _{\rm sfr}}(\Omega -\Omega _{p}), \end{equation} \tag{ 8 }$$

where Ω_p is the pattern frequency of the spiral density wave. In the limit of small Ω_p and for α_sfr = 1 we recover Equation (2). Increased cloud collision rates and increased perturbation growth rates in the arms, where Q is locally lowered, have been suggested as the star formation triggering mechanism (e.g., Dobbs 2008). Kim et al. (2008) examined the development of thermal instabilities in galactic spiral shocks.

One prediction of these theories is a correlation of SFR with the density wave amplitude. However, this is not observed (Elmegreen & Elmegreen 1986; McCall & Schmidt 1986; Kennicutt 1989). Furthermore, such theories have difficulty explaining star formation in galaxies where there is a lack of organized star formation features, as in flocculent spirals (e.g., Thornley & Mundy 1997a, 1997b; Grosbol & Patsis 1998). GMCs are present and SFRs are similar to those systems where star formation is organized into spiral patterns. This suggests stellar disk instabilities, which create spiral density waves, and gas instabilities, which lead to GMCs and large-scale star formation, are decoupled (Kennicutt 1989; Seiden & Schulman 1990).

Tan (2000) proposed a model in which the majority of star formation in disk galaxies occurs in the pressurized regions triggered by GMC–GMC collisions. GMCs are observed to occupy a very thin vertical distribution in the Galaxy (Stark & Lee 2005), which is similar to the actual sizes of the clouds. In this effectively two-dimensional system, collisions are set by galactic shear at impact parameters of about the tidal radius of the clouds (Gammie et al. 1991), and for a Q ∼ 1 disk with a relatively large mass fraction, ∼1/2, in GMCs and associated gas, collisions are expected to occur approximately every 20% of an orbital time. In this way, for flat rotation curve galaxies, Equation (2) is recovered, and the cloud collision mechanism is the link between the global galactic dynamics and the parsec-scale star-forming clumps of GMCs. One prediction of this model is a reduction in the SFR per orbital time for galaxies with lower rates of shear, i.e., those with rising rotation curves that are closer to solid body rotation. This model also requires GMCs to be relatively long-lived, i.e., several tens of Myr.

Krumholz & McKee (2005) proposed that the SFR is regulated by turbulence in a galaxy's molecular gas, so that a fixed fraction of the turbulent, molecular gas is converted to stars per free-fall time of the GMC. To recover Equation (2), one has to assume that GMC dynamical timescales are a fixed fraction of galactic dynamical timescales and that GMC mass fractions are constant in different galactic systems (see also Wada & Norman 2007).

The large range of spatial scales involved from GMCs to global galactic properties (a range of ∼10⁴ from ≲ 10 parsec-scale clouds to the galactic scale) has made it extremely difficult for simulations to encompass GMC formation on global scales. Work that does perform this is either limited to two dimensions (Shetty & Ostriker 2008) or has to assume a fixed two-phase medium for the ISM (Dobbs 2008). While the clouds do largely reside in the plane of the disk (which is in essence a two-dimensional system), cloud collisions and feedback can eject gas from the surface, a process that controls the pressure of the ISM from which the clouds are forming (McKee & Ostriker 1977; Cox 2005). Similarly, a fixed phase ISM results in clouds being created from gas with a predetermined structure, from which is it harder to discern the main physical effects controlling their formation and evolution. Three-dimensional models that contained a self-consistent multiphase atomic ISM on global scales were performed by Tasker & Bryan (2006, 2008) and Wada & Norman (2007), but either did not reach the resolution needed to probe below the most massive GMCs (e.g., the simulations of Tasker & Bryan (2008) had a limiting resolution of 25–50 pc), or only considered the inner galaxy (e.g., the simulations of Wada & Norman (2007) were restricted to r < 2.56 kpc). Robertson & Kravtsov (2008) presented a global SPH simulation, with similar effective resolution to the simulation of Tasker & Bryan (2008) (A. V. Kravtsov 2008, private communication) of star formation in disk galaxies including a multiphase ISM, and subgrid models for various forms of stellar feedback and molecular hydrogen formation and destruction (see also Gnedin et al. 2009).

Local models can resolve down to smaller scales and, if set up in a shearing box, can approximate the effects of galactic shear on cloud formation and evolution (e.g., Kim & Ostriker 2001, 2006, 2007). Other examples of local models have focused on cloud formation from imposed colliding flows (e.g., Heitsch et al. 2008), including the effects of magnetic fields (Heitsch et al. 2009). Other groups have studied cloud formation in the context of the local interplay between supernova feedback and a turbulent, multiphase ISM (e.g., Slyz et al. 2005). Simulations including the nonequilibrium formation and destruction of H₂ molecules have been carried out by Glover & Mac Low (2007a, 2007b).

While these local models are very important for studying many aspects of GMC formation and evolution, they cannot explore the global evolution of GMCs as they orbit through the disk or measure how GMC (and star formation) properties are related to the more readily observed global galactic properties such as mean gas mass surface density. Trends with galactocentric radius are also easier to study in global models.

In this paper, we present results from a three-dimensional global galaxy (∼32 kpc box containing a gravitationally unstable galactic disk with diameter 20 kpc) simulation that tracks the formation and evolution of clouds in a self-consistent multiphase atomic ISM. Our main simulation has a limiting resolution of 7.8 pc. As discussed below, we do not at this stage include the physics of star formation and feedback. We also do not include magnetic fields, cooling below 300 K, or any treatment of the formation and destruction of molecules. Our definition of "GMCs," which we also refer to interchangeably as "clouds," is discussed in detail below, but basically involves a density threshold of n_H ⩾ 100 cm^-3. In real galaxies, the structures that correspond to this definition will typically be composed of a mixture of atomic and molecular gas, but we argue that the bulk dynamical properties of the clouds can still be reasonably well captured by our simulations.

The simulations performed in this paper were run using Enzo: a three-dimensional adaptive mesh refinement (AMR) hydrodynamics code (Bryan & Norman 1997; Bryan 1999; O'Shea et al. 2004). Enzo has been used previously in galactic disk simulations (e.g., Tasker & Bryan 2006, 2008), where it successfully produced a self-consistent multiphase atomic ISM, consisting of a wide range of densities and temperatures. Grid codes are particularly adept at modeling multiphase gases since the grid cells form natural boundaries allowing the gas to evolve with a range of temperatures, densities, and pressures. Particle-based techniques also struggle to resolve fluid instabilities (Tasker et al. 2008) and can suffer from overmixing problems unless specific algorithmic steps are taken.

Enzo evolved the gas using a three-dimensional version of the Zeus hydrodynamics algorithm (Stone & Norman 1992). This scheme uses an artificial viscosity term to model shocks and the variable associated with this, the quadratic artificial viscosity, was set to 2.0 (the default) for all simulations. Radiative gas cooling followed the solar metallicity cooling curve of (Sarazin & White 1987) down to temperatures of T = 10⁴ K and extends down to T = 300 K using rates from Rosen & Bregman (1995). These temperatures take us to the upper end of the atomic cold neutral medium (Wolfire et al. 2003). In this study, we do not include the formation and destruction of molecules or any cooling processes below 300 K. By a combination of dust and molecular cooling, gas in real GMCs reaches temperatures of ∼10 K, more than an order of magnitude below our minimum radiative cooling temperature. However, since we are not able to resolve the detailed internal structure of clouds and their internal turbulence (a typical cloud with a diameter of ∼100 pc would only have ∼13 cells across each linear dimension in our highest resolution simulation) and since we are not including magnetic pressure support, we view our temperature floor of 300 K as being equivalent to imposing a minimum one-dimensional signal speed equal to the sound speed c_s = (γP/ρ)^1/2 = (γkT/μ)^1/2 = 1.80(T/300 K)^1/2 km s⁻¹, where γ = 5/3, μ = 1.273m_p (for an assumed n_He = 0.1n_H). This signal speed is somewhat smaller than observed internal velocity dispersions of GMCs. In fact we will see that, even though our simulation does not include feedback, our simulated GMCs typically attain internal nonthermal (i.e., turbulent) velocity dispersions much greater than this minimum.

The galaxy was modeled in a three-dimensional simulation box of side 32 kpc with isolated gravitational boundary conditions and outflow fluid boundaries. For our main, high-resolution simulation, the root grid was 256³ with an additional four levels of refinement, producing a minimum cell size of 7.8 pc. To examine how simulation results depend on resolution, we also carried out medium and low-resolution runs with a minimum cell sizes of 15.6 pc and 31.2 pc, respectively. Refinement of a cell occurred when the Jeans length drops below four cell widths, in accordance with the criteria suggested by Truelove et al. (1997) for resolving gravitational instabilities. Resolution of the fragmentation is discussed further in Section 2.4. The disk was allowed to evolve for 324 Myr, about 1.3 orbital periods at r = 8.0 kpc (see Section 2.2) and substantially longer than the local free-fall time of the initial mid-plane gas at this location, which was about 60 Myr.

In this paper, the first of a series, we examine the formation and evolution of GMCs in a flat rotation curve galactic disk without the presence of star formation and feedback mechanisms, which we defer to future papers. Since the amount of mass removed by star formation is expected to be only a few percent per local free-fall time (Zuckerman & Evans 1974; Krumholz & Tan 2007), the gas depletion over the course of the simulation time would be relatively modest. The structure of the ISM and the GMC population we derive can be regarded as that resulting from the limiting case when there is very little coupling of stellar feedback (including FUV heating) to the ISM. Gravitational scattering of bound clouds in a shearing disk and dissipative cloud collisions are the processes that will act in our simulations to extract orbital energy and regulate the ISM structure.

To mimic the present-day Milky Way, the simulated galaxy was set up as an isolated disk of gas orbiting in a static background potential which represented both dark matter and a stellar disk component. This minimized the effects from the evolution of the galaxy on the interstellar structure, allowing us to investigate the formation and evolution of GMCs in a steady-state environment. The background potential was chosen to produce a flat rotation curve at r ≫ r_c, where the core radius was set to be r_c = 0.5 kpc. The form of the potential is (Binney & Tremaine 1987)

$$\begin{equation} \Phi = \frac{1}{2}v_{c,0}^2\ln \left[\frac{1}{r_c^2}\left(r_c^2 + r^2 + \frac{z^2}{q_\phi ^2}\right)\right], \end{equation} \tag{ 9 }$$

where v_c,0 is the constant circular velocity in the limit of large radii, here set equal to 200 km s^-1, r and z are the radial and vertical coordinates, respectively, and the axial ratio of the potential field is q_ϕ = 0.7. The form of the circular velocity of the disk, v_c, is then given by

$$\begin{equation} v_c=\frac{v_{c,0}r}{\sqrt{r_c^2 + r^2}}. \end{equation} \tag{ 10 }$$

The initial radial profile of the gas mass surface density in the disk was chosen to give a constant value of the Toomre Q parameter for gravitational instability (Toomre 1964):

$$\begin{equation} \Sigma _g = \frac{\kappa \sigma _g}{\pi G Q}, \end{equation} \tag{ 11 }$$

where σ_g, the one-dimensional velocity dispersion of the gas, is equivalent to the sound speed c_s for the case of an razor thin disk with only thermal pressure. The particular choices of Q at different radii are discussed below. For our simulated gas disk we define σ_g ≡ (σ²_nt + c²_s)^1/2, averaged over the mass in particular regions of the disk, where σ_nt is the one-dimensional velocity dispersion of the gas motions in the plane of the disk after the subtraction of the circular velocity, i.e., representing nonthermal motions. Our disks are initialized with σ_nt = 0.

The initial vertical profile of the gas was set proportional to sech²(z/z_h), where z_h is the vertical scale height, which was assumed to vary, i.e., increase, with galactocentric radius based on observations of the H i in the Milky Way presented in Binney & Merrifield (1998). At the solar radius of 8 kpc, z_h = 290 pc. Then we have Σ_g = ∫^∞_−∞ρ₀ sech²(z/z_h)dz = 2ρ₀z_h, where ρ₀ is the midplane (z = 0) density, so that the gas distribution becomes

$$\begin{equation} \rho (r,z) = \frac{\kappa \sigma _g}{2 \pi G Q z_h} {\rm sech}^2{\left(\frac{z}{z_h}\right)}. \end{equation} \tag{ 12 }$$

The initial disk profile is divided radially into three sections. In our main region of interest, between radii of r = 2–10 kpc (that is, encompassing the part of the Galaxy inside the solar circle with radius of 8 kpc) Σ_g is set so that Q = 1 if σ_g = c_s were equal to 6 km s⁻¹, similar to the observed velocity dispersion of the ISM (Kennicutt 1998; Stark & Brand 1989). The actual initial velocity dispersion (i.e., sound speed) is set to c_s = 9.0 km s⁻¹, so that T_init ∼ 7450 K and Q = 1.5. The threshold of gravitational instability is crossed by allowing the gas to cool. The other regions of the galaxy, from 0 to 2 kpc and from 10 to 12 kpc, are initialized in a similar way, but are designed to be gravitationally stable with Q = 20 for a flat rotation curve if c_s were to equal 6 km s⁻¹ (although note the rotation curve is not flat in the center). Their initial temperatures are set to the same value as in the main disk region. Even after cooling to the temperature floor of 300 K, for which c_s = 1.8 km s⁻¹, these regions remain gravitationally stable, except for a small amount of gas that collects at the galaxy center. Beyond 12 kpc, the disk is surrounded by a static, very low density medium, that has negligible influence. Only the main disk region from 2 to 10 kpc is analyzed in this paper, and in fact, due to boundary effects, we typically restrict analysis to radii from 2.5 to 8.5 kpc.

In total, the gas mass was 6 × 10⁹ M_☉. Note, this is a factor of 10 smaller than the disks presented in Tasker & Bryan (2008).

The clouds in the galactic disk were located using a number density of H nuclei threshold of n_H,c = 100 cm⁻³, similar to the mean (volume-averaged) densities of typical Galactic GMCs. Also recall that the formation of molecules is not being followed in this simulation. When we refer to "clouds" or "GMCs" we are describing the gas that has achieved densities of n_H ⩾ n_H,c in the atomic phase. We expect that most of the gas above this density would form molecules via surface reactions on dust grains, although there could still be substantial (∼ equal mass) atomic components present, as are observed in and around Galactic GMCs (Wannier et al. 1991; Blitz 1990).

The process of locating and tracking the clouds in the disk at a given time in the simulation is outlined below.

• 1.  
  Clouds were identified as peaks in the baryon density field with n_H ⩾ n_H,c = 100 cm⁻³. If two peaks were ⩽4 minimum cell widths apart, i.e., about 30 pc in our high-resolution run, then only the higher density peak was retained for cloud definition.
• 2.  
  Nonpeak cells with n_H ⩾ n_H,c were assigned to the nearest peak that was connected to them by cells with n_H ⩾ n_H,c, i.e., clouds are continuous structures.

Once the cloud had been defined, properties including the cloud's center of mass, velocity, and angular momentum were calculated. Multiple clouds could exist in the same continuous density structure if it contained more than one well-separated peak.

To follow the evolution of the clouds, simulation outputs were analyzed every 1 Myr and the clouds mapped between outputs with a tag number assigned to each cloud to follow its life through the simulation. To map a cloud between times t₀ and t₁, the code performed the following steps.

• 1.  
  Assuming linear motion, a predicted position of each cloud's center of mass found at time t₀ is calculated for the cloud at t₁.
• 2.  
  A volume of radius 50 pc, larger than the expected deviations from linear motion due to typical accelerations (e.g., over 1 Myr these are about 20 pc at r = 2 kpc due accelerations in the galactic potential), centered on the predicted position of each cloud is searched for clouds present at time t₁. If multiple clouds are found in this region at t₁, the nearest one is chosen to be associated with the cloud from t₀. In cases where two or more t₀ clouds are associated with the same t₁ cloud, the nearest one is matched and the volume around the predicted positions of the other t₀ clouds is searched for alternative candidates.
• 3.  
  If no clouds at t₁ are associated with a t₀ cloud, then a volume with radius equal to 3× the average radius of the t₀ cloud is searched. This allows for large, extended clouds whose radius may be ≳50 pc and whose centers may have shifted due to external perturbations.
• 4.  
  Any t₀ clouds remaining unassigned may have merged with neighboring t₀ clouds. A volume of radius 2× the average radius of each t₁ cloud is searched for unassigned t₀ clouds. This value was chosen to be lower than for the previous step since a recently merged cloud is likely to be fairly extended. If found, a merger is declared between the t₀ cloud previously associated with the t₁ cloud and the unassigned t₀ cloud. The t₁ cloud inherits the tag number from the more massive of the two merged clouds. If there were multiple t₁ clouds close to an unassigned t₀ cloud, the closest was chosen. Multiple mergers involving more than two clouds were possible, though were typically rare.
• 5.  
  Any t₀ clouds remaining unassigned after these steps are declared to have been destroyed by nonmerger processes.

Note that this method assumes that any cloud that has been destroyed in close proximity to another cloud has suffered a merger. This assumption works well in the current simulation, which has no stellar feedback, but it remains to be tested once feedback processes are operating.

Figure 1 illustrates two possible examples of the cloud mapping. The blue circles show clouds at t = t₀ and the red circles show clouds located 1 Myr later at t = t₁. In panel (A), the predicted position of cloud 1 at t = t₁ is marked by the black dotted open circle. The surrounding region with radius 50 pc is searched for clouds present at t₁ and two are found. The closest of these is identified to be cloud 1 at t₁ whereas the other cloud is either a newly formed cloud or a different cloud present at t₀. Panel (B) shows how the algorithm deals with a more complex situation involving a merger. Here, a cloud at t₁ is associated with both t₀ clouds 2 and 3. The code checks there are no unassigned clouds at t₁ that are a viable match for cloud 2 or 3 but finding none, it declares this a merger event. Note that the cloud detection algorithm does not allow distinct clouds to be found within four minimum cell widths of each other.

Zoom In Zoom Out Reset image size

Figure 2 shows how the cloud identification and tracking work in the simulation. The two images show the gas density in a one-cell thick slice of a square patch of the galaxy disk midplane 2 kpc across and 5 kpc from the galactic center, taken 1 Myr apart. Contour lines mark the density threshold of n_H,c ⩾ 100 cm⁻³ and the squares and triangles (prograde and retrograde rotators, see Section 4.2) show the center of mass of each of the clouds identified, with their tag number written below. Areas where contours or a peak appear to be visible without a cloud are due to the cloud's center of mass being above or below the slice shown. The cloud tracking algorithm correctly maps clouds between the two outputs and shows a merger of clouds 9282 and 13580 on the right-hand side.

Zoom In Zoom Out Reset image size

Finite computational resources impose a limit on the number of levels of refinement we are able to employ in the simulation. Cells at the highest level of refinement are themselves not refined further, which means they can grow in mass and density to the point that the Truelove et al. (1997) condition to avoid artificial fragmentation is no longer satisfied, i.e., that the Jeans length, λ_J, associated with gas in each cell should be at least four times larger than the size of the cell, Δx. In Figure 3, we show the mean density of cells as a function of their Jeans number, N_J = Δx/λ_J. This figure indicates that for the high-resolution simulation the Truelove et al. criterion is satisfied for densities up to n_H ≃ 70 cm^-3, which is close to the threshold density used to define our GMCs. Thus, we expect that the gravitational fragmentation process that leads to the formation of GMCs is reasonably well resolved in our simulation. It is not the goal of this paper to follow the fragmentation inside GMCs, and for this reason a minimum effective sound speed of about 1.8 km s⁻¹ is imposed (see above). It is possible that large GMCs may suffer spurious fragmentation in the simulation, however it should be noted that these GMCs typically have internal velocity dispersions that are much larger than the sound speed (see Section 4.2). We examine the effect of simulation resolution on the GMC mass function and other properties in Section 4.2.

Zoom In Zoom Out Reset image size

As the disk becomes gravitationally unstable, local patches increase in density, many of which reach our threshold value of n_H,c = 100 cm⁻³, at which point we recognize them as a GMC. The formation rate of these clouds over the course of the simulation is shown in Figure 8. There is a burst of cloud formation associated with the initial fragmentation of the disk, with the rate peaking at about 400 Myr⁻¹ at a simulation time of t = 100 Myr. By t = 140 Myr, cloud formation settles down to a slower and nearly constant rate of ∼200 Myr⁻¹. This is partly due to there being a smaller noncloud diffuse gas reservoir and partly due to heating of the disk by gravitational interactions of the cloud population. We are most interested in the cloud and ISM properties after t = 140 Myr, since before this time their properties are influenced by the artificially smooth initial conditions and their fragmentation via the Toomre ring instability.

Zoom In Zoom Out Reset image size

Figure 8 also shows the evolution of the total number of clouds and the merger rate of clouds. This information is used to calculate the mean merger time of clouds, t_merger, relative to their orbital time, as a function of galactocentric radius (Figure 9). This calculation assumes mergers are all binary mergers (the non-binary merger fraction is indeed negligibly small). From Figure 9 we see that the average merger time settles to be a small fraction, ∼0.2 of the orbital time, with only a modest dependence on galactocentric radius and simulation time for t ⩾ 140 Myr. This confirms the results of Tan (2000), who estimated t_merger/t_orbit = 0.2 based on the simplified calculation of cloud orbits and binary interactions by Gammie et al. (1991). Note that the orbital time for v_c = v_c,0 = 200 km s^-1 is

$$\begin{equation} t_{\rm orbit} = 123 (r/4\;{\rm kpc}) (v_c/200\;{\rm km\;s^{-1}})^{-1} {\rm Myr}, \end{equation} \tag{ 13 }$$

so that the typical merger time is only ∼25 Myr, which is relatively short compared to traditional estimates of cloud collision timescales of hundreds of Myr that do not allow for the essentially two-dimensional geometry of the GMC population, gravitational focusing and that cloud interaction velocities are set by galactic shear and are larger (by about a factor of 2) than the local cloud velocity dispersion (see discussion in McKee & Ostriker 2007). Our derived collision timescales are also shorter than many estimates of GMC lifetimes (e.g., McKee & Williams 1997; Matzner 2002), which suggests that collisions are likely to be important even when feedback mechanisms are taken into account, i.e., an individual GMC is just as likely to have its properties dramatically altered by a merger than by a destructive mechanism such as supernova or ionization feedback.

Zoom In Zoom Out Reset image size

The nearly constant ratio of t_merger/t_orbit means that the cloud collision rate is tied to the global galactic dynamical timescale and, if cloud collisions are the trigger for the majority of star formation, then this provides a natural mechanism to explain the global SFR–gas content correlations observed by Kennicutt (1998) in which the overall star formation efficiency per galactic dynamical time is small, while having star formation occur mostly on small scales at relatively high efficiency in a highly clustered mode (Tan 2000). The relation of star formation to cloud collisions in global galaxy simulations will be investigated in Paper II.

The following caveats should be considered in the above estimates of the merger timescale. The simulated merger rate depends on the mass function and gravitational boundedness of the GMCs. More massive and more gravitationally bound GMCs will experience higher merger rates, although for typical cloud mass functions most of their collisions will be with lower-mass clouds. In Figure 9, we also show the merger timescale for clouds with M_c ⩾ 10⁶ M_☉, which is about a factor of 2 shorter than that of the average cloud. We will see below that the simulated GMC population mass function and gravitational boundedness are similar to observed values. The total number of clouds is not fully resolved in these simulations (the most common cloud mass is several×10⁵ M_☉ in the high resolution run at late times). Higher resolution simulations are likely to resolve larger number of clouds, which would increase the merger rate, although not necessarily the total gas mass that has been compressed by collisional shocks. The inclusion of stellar feedback would likely lead to a reduced total mass fraction in clouds and a reduced gravitational boundedness of those clouds that do form, and these effects would reduce the merger rate.

The cloud mass function is shown in Figure 10(a). Since there is no star formation feedback in these models, it is difficult to destroy gravitationally bound gas clouds, except via their merger into more massive clouds or via disruptive collisions. Unbound and weakly bound clouds can also be destroyed by shearing tidal forces due to grazing and close cloud interactions. As a result, we expect the mean cloud mass to grow over the course of the simulation, and this behavior is seen in Figure 10(a). Note, however, that the peak of the cloud mass function is at 6 × 10⁵ M_☉ (for uniform binning in logM_c) and this does not vary significantly from t=100 to 300 Myr. Note the initial Toomre mass in a Q = 1 disk is M_T ≡ Σ_gλ²_T = 8.9 × 10⁶(σ_g/6 km s^-1)³(r/4 kpc)(v_c/200 km s^-1)⁻¹M_☉, about an order of magnitude larger than the peak of the simulated GMC mass function. The Jeans mass in the disk is M_J = πΣ_gλ²_Jσ⁴_g/(G²Σ_g) → 2.2 × 10⁶(σ_g/6) km s^-1)³(r/4 kpc)(v_c/200 km s^-1)⁻¹(Q/1) M_☉, where λ_J ≡ σ²_g/(GΣ_g) is the shortest wavelength permitting gravitational instability in a thin, nonrotating disk (Kim & Ostriker 2001). The gas tends to cool to the temperature floor of 300 K, for which σ_g ≃ 2 km s^-1 and thus Q ≃ 1/3, resulting in M_J ≃ 3 × 10⁴ M_☉. The effect of simulation resolution on the cloud mass function at t = 250 Myr is shown in Figure 11(a). The number of GMCs appears reasonably well resolved down to masses of about 2 × 10⁶ M_☉.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The simulated cloud mass function can be compared to the observed mass spectrum of GMCs in the Milky Way and M33, which can be fitted by a power law $d{\cal N}_c /(d \ln M_c) \propto M_c^{-\alpha _c}$. In the Milky Way, α_c is observed to be ≃0.6–0.8 (Williams & McKee 1997), whereas a steeper index of ≃1.6 appears to hold in M33 (Rosolowsky et al. 2003). After collision and agglomeration processes have had enough time to operate, i.e. by t = 200 Myr, the shape of the simulated cloud mass function (in the mass range ∼(0.5–10) × 10⁶ M_☉ relevant to observations) approaches that seen in the Milky Way GMC population, indicating that most of the total cloud mass is contained in the more massive clouds. It remains to be determined whether these processes would still be as important for shaping the cloud mass function if stellar feedback processes were also operating, which would tend to reduce cloud lifetimes.

Williams & McKee (1997) also derived the normalization of the Galactic GMC mass function and noted there appears to be a truncation above a molecular mass of M_u = 6 × 10⁶ M_☉ (see Equation (6)). Allowing for an equal amount of atomic gas associated with GMCs (Blitz 1990), this truncation would occur at 1.2 ×  10⁷ M_☉. We see that at late simulation times, the cloud mass function has grown beyond this value, though there are very few clouds with M_c > 10⁷ M_☉ (see also Figure 11(a)). The presence of these very massive clouds is likely due to the lack of stellar feedback cloud destruction mechanisms in these current simulations, especially the ability of stars to destroy their natal clouds. Williams & McKee (1997) estimate there are ∼100–200 inner Milky Way GMCs with M_c > 10⁶ M_☉, and about 1000 with M_c > 10⁵ M_☉. If one allows for an equal mass of atomic gas associated with the observed molecular gas of these GMCs, then in comparison we find about 400 clouds with M_c > 2 × 10⁶ M_☉, a factor of 2–4 higher than the Williams and McKee estimate. We do not adequately resolve the number of 2 × 10⁵ M_☉ clouds (Figure 11(a)) to make a useful comparison at that mass scale. These results may indicate that the number of massive GMCs in our simulation is a factor of a few larger than in the Milky Way. Indeed the mass fraction of gas in GMCs in the simulation is about 0.69 (Figure 6), somewhat higher than the mass fraction implied by the analysis of Wolfire et al. (2003). Lack of feedback processes in the simulation is an obvious potential cause of this discrepancy. Another potential contributing factor is our somewhat arbitrary choice of n_H,c = 100 cm^-3 for the cloud threshold density. Note that M33, being a much smaller galaxy than the Milky Way, has fewer massive GMCs.

Figure 10(b) shows the average radius of the clouds defined as R_c,A ≡ (A_c/π)^1/2, where A_c is the projected area of a cloud in the y–z plane (i.e., as it would be observed by an observer embedded in the plane of the galaxy). The most common radius for a cloud is around 20–30 pc (for uniform linear binning in R_c,A). By 300 Myr, agglomeration processes have created a population of larger clouds. It is clear from Figures 10(a) and (b) that the numbers of these larger, more massive clouds are steadily increasing in time, and in this respect a steady state has not been reached. Such a steady state likely requires feedback processes. Nevertheless, the typical sizes of the cloud population are similar to observed sizes of GMCs in the Milky Way and other galaxies, such as M33 (Rosolowsky et al. 2003). The effect of simulation resolution on the distribution of cloud sizes is shown in Figure 11(b).

Figure 10(c) shows the distribution of mass surface density of the clouds, which is defined as Σ_c ≡ M_c/A_c. The peak of the distribution (for uniform binning in logM_c) is at about 300 M_☉ pc^-2, and does not change significantly during the simulation. Although the number of more massive clouds is increasing, these are also larger, and so have similar values of Σ_c. Most clouds are within a factor of 3 of this peak value. This surface density is similar to the mean value of ∼200 M_☉ pc^-2 derived by Solomon et al. (1987) from a ¹²CO survey of the Galaxy. Heyer et al. (2008) have derived smaller values ∼100 M_☉ pc^-2 based on better sampled ¹³CO surveys. To compare these observed values, which are based solely on molecular line (CO) emission, with our simulated clouds, one should also include their associated atomic gas, which may contribute at about the factor of 2 level. We conclude that, even without disruptive effects of star formation feedback, our simulated GMCs attain values of Σ_c that are similar to observed values. This may indicate that the dominant source of turbulent pressure support in GMCs is injected via cloud collisions and interactions (i.e., turbulent converging flows) rather than via star formation feedback. Note that Joung & Mac Low (2006) found that supernova-driven turbulence was insufficient to explain the observed ISM velocity dispersions. The effect of simulation resolution on the distribution of Σ_c is shown in Figure 11(c): the peak of the distribution is relatively insensitive to resolution.

Figure 10(d) shows the distribution of GMC virial parameters (see Equation (7)), α_vir ≡ 5σ²_cR_c,A/(GM_c), where σ_c is the mass-averaged one-dimensional velocity dispersion of the cloud, i.e., σ_c ≡ (c²_s + σ²_nt,c)^1/2, where σ_nt,c is the one-dimensional rms velocity dispersion about the cloud's center-of-mass velocity (Bertoldi & McKee 1992). A virial parameter of unity implies a spherical, uniform cloud with negligible surface pressure and magnetic fields is virialized, so that its total kinetic energy is half the magnitude of the gravitational energy. Surface pressure terms and a centrally concentrated density distribution tend to raise the value of α_vir corresponding to virial equilibrium (see Equation (7) and the discussion above). We see that the distribution of cloud virial parameters peaks at α_vir < 1, and decreases slightly as the simulation progresses: the peak of the virial parameter distribution at 300 Myr is about 0.6. This indicates that the clouds in the simulation are typically gravitationally bound, perhaps somewhat more strongly bound than observed GMCs, though there is still a significant population of unbound clouds. The effect of increasing simulation resolution on the distribution of α_vir is shown in Figure 11(d): the peak of the distribution moves to larger values and the fraction of unbound clouds grows.

Figure 10(e) shows the distribution of the vertical (z) component of specific angular momentum, j_z, of the clouds. The fraction of retrograde (i.e., j_z < 0) clouds grows at late times as clouds have had more opportunities to suffer mergers and close interactions (see also Dobbs 2008). The vertical line indicates the value of j_z of a spherical (≃110 pc radius) region of the initial conditions at galactocentric radius r = 4 kpc containing 10⁶ M_☉. (Kim et al. 2003 derive an analytic expression for j_z ∝ ΩR²_c,A for two-dimensional patches of galactic disks.) We see that GMCs typically have much smaller amounts of angular momentum compared to similar masses of diffuse gas that would be spread out over larger scales and be more susceptible to galactic differential rotation. The effect of simulation resolution on the distribution of j_z is shown in Figure 11(e): at higher resolution it becomes more sharply peaked toward small value of j_z.

Figure 10(f) shows the distribution of cloud center-of-mass vertical positions relative to the disk midplane, i.e., where the vertical acceleration due to the galactic potential is zero. The rms heights of the cloud population at 100, 200, and 300 Myr are 13, 25, and 51 pc, respectively, growing at late times because of cloud–cloud gravitational scattering. Note that Stark & Lee (2005) derived a vertical scale height of Milky Way GMCs of ≲35 pc. The effect of simulation resolution on the vertical distribution of clouds is shown in Figure 11(f).

Figure 12 shows the velocity dispersion versus size relation for clouds at t = 250 Myr. The effects of using the full velocity information to derive a mass-averaged one-dimensional internal velocity dispersion versus just using radial velocities at a particular viewing angle are examined. We also consider the effect of a mass cut, i.e., restricting the analysis only to the most massive clouds. Such effects are important in influencing the normalization and slope of the line-width size relation. Our simulated cloud population has a similar, though somewhat higher, velocity dispersion versus size relation to observed GMCs, such as those observed in M33 by Rosolowsky et al. (2003), in the Milky Way by Solomon et al. (1987), or the average of a number of extragalactic systems, $\sigma _c = \sigma _{\rm pc}(R/{\rm pc})^{\alpha _\sigma }$ with α_σ = 0.60 ± 0.10 and σ_pc = 0.44 ± 0.15 km s^-1, derived by Bolatto et al. (2008). The larger velocity dispersion of the simulated clouds at a given size scale may be due to the lack of star formation feedback causing them to be more tightly gravitationally bound than real clouds. A discrepancy may also result due to differences in how the effective radius of clouds is evaluated for the simulation clouds and the real clouds. The dispersion of the simulated clouds about the best-fit relation is about 2.2 km s^-1 (for the sample in the lower-right panel of Figure 12), corresponding to about 40% of σ_c. Solomon et al. (1987) found a dispersion of about 30% in their sample of Galactic GMCs (see also Heyer & Brunt 2004).

Zoom In Zoom Out Reset image size

Figure 13 shows the distribution of the angle, θ, between the cloud angular momentum vector and the galactic rotation axis at different times in the simulation. At early times, in the initial fragmentation phase clouds are all born with prograde rotation (0° < θ < 90°), i.e., in the same sense as the galactic rotation. At later times, cloud interactions lead to a build up of a retrograde (90° < θ < 180°) cloud population, which is about 30% of the total.

Zoom In Zoom Out Reset image size

To compare to observed systems, one needs to account for the viewing angle to the galaxy and the fact that only line-of-sight velocities are used to determine rotation. Figure 14 shows these effects by plotting the distribution of position angles, θ', of cloud rotation axes derived from radial velocities along sight lines at different inclination angles to the simulated galaxy at t = 250 Myr, when the retrograde fraction of θ is 0.28. These results indicate that the derived retrograde fraction depends on inclination angle: e.g., for the 0° (edge on) case the retrograde fraction inferred from θ' is about twice that as measured by the actual angular distribution of θ. The figure also shows the observational results of Rosolowsky et al. (2003), who find about 40% (18 out of 45) of the GMCs in M33 have apparently retrograde rotation. Their results are broadly consistent with our simulation results (note their slightly larger bin size), but we emphasize that we can only make a qualitative comparison at this stage, since our model galaxy was not set up to mimic the details of M33 (e.g., its rotation curve), nor does it yet include the effects of stellar feedback. Nevertheless cloud formation in a self-gravitating, shearing disk and cloud evolution via cloud–cloud interactions appear to be promising mechanisms to help explain the observed angular distributions of cloud angular momentum vector orientations.

Zoom In Zoom Out Reset image size

Our method of tracking clouds from one output time to the next allows us to study the evolution of cloud properties as a function of cloud age, t'. Figure 15, analogous to Figure 10, shows how the distribution of cloud masses, radii, mass surface densities, virial parameters, vertical component of specific angular momentum, and vertical positions changes with cloud age (t' = 0–1, 9–10, 49–50, 99–100 Myr) for clouds born after t = 140 Myr, i.e., in the fully fragmented phase of the simulation. Younger, especially newly formed, clouds tend to have lower masses, slightly smaller radii, lower mass surface densities, and larger virial parameters (i.e., are less gravitationally bound). However, young clouds have very similar distributions of the vertical component of specific angular momentum and of vertical position above and below the disk midplane.

Zoom In Zoom Out Reset image size

Figure 16 shows the distribution of the angles, θ, between cloud angular momentum vectors and the galactic rotation axis for different cloud ages for clouds born after t = 140 Myr in the fully fragmented stage. The fraction of retrograde clouds is about 25% of the total and nearly independent of cloud age, though it has risen to about 30% for 100 Myr old clouds. At times t > 140 Myr, clouds are being born from overdense gas that has already been disturbed by cloud interactions. At earlier times (t < 140 Myr) we find much smaller retrograde fractions, ∼10%. These results, especially for the later simulation times, may depend on the absence of feedback in the simulation, which would lead to a higher mass fraction in lower-density gas. Clouds forming from such gas are expected to form with a higher degree of prograde rotation, since the initial conditions of the clouds are more dispersed and thus influenced by galactic shear. The influence of feedback on these aspects of cloud formation will be investigated in a future paper.

Zoom In Zoom Out Reset image size