THE SCALE DEPENDENCE OF THE MOLECULAR GAS DEPLETION TIME IN M33

The observed correlation between gas and star formation rate (SFR) surface densities (the "star formation law") is one of the most widely used scaling relations in extragalactic astronomy (e.g., Schmidt 1959; Kennicutt 1998). However, its connection to the fundamental units of star formation, molecular clouds, and young stellar clusters remains poorly understood. On the one hand, averaged over substantial areas of a galaxy, the surface density of gas correlates well with the amount of recently formed stars (e.g., Kennicutt 1998). On the other hand, in the Milky Way, giant molecular clouds (GMCs), the birthplace of most stars, and H ii regions, the ionized interstellar medium (ISM) around (young) massive stars, are observed to be distinct objects. While they are often found near one another, the radiation fields, stellar winds, and ultimately supernovae make H ii regions and young clusters hostile to their parent clouds on small (∼10 pc) scales. Thus, while correlated on galactic scales, young stars and molecular gas are in fact anti-correlated on very small scales. The details of the transition between these two regimes remain largely unexplored (though see Evans et al. 2009).

Recent observations of nearby galaxies have identified a particularly tight correlation between the distributions of molecular gas (H₂) and recent star formation on ∼kpc scales (Murgia et al. 2002; Wong & Blitz 2002; Kennicutt et al. 2007; Bigiel et al. 2008; Leroy et al. 2008; Wilson et al. 2008). While the exact details of the relation are still somewhat uncertain, in the disks of spiral galaxies the parameterization seems to be a power law, Σ_SFR = aΣ^b_H2, with power law index b ≈ 1–1.7 and coefficient a corresponding to H₂ depletion times of ∼2 Gyr in normal spirals.

Both parts of this relation, the surface densities of H₂ and recent star formation, resolve into discrete objects: GMCs and H ii regions, young associations, and clusters. In this paper, we investigate whether the ∼kpc H₂–SFR relation is a property of these individual regions or a consequence of averaging over large portions of a galactic disk (and the accompanying range of evolutionary states and physical properties). To do so, we compare CO and extinction-corrected Hα at high spatial resolution in the nearby spiral galaxy M33. We examine how the ratio of CO-to-Hα changes as a function of region targeted and spatial scale. M33 is a natural target for this experiment: it has favorable orientation and is close enough that peaks in the CO and Hα maps approximately correspond to individual massive GMCs (Rosolowsky et al. 2007) and H ii regions (Hodge et al. 2002).

Perhaps not surprisingly, we find that the ratio of CO luminosity (a measure of the molecular gas mass) to extinction-corrected Hα flux (a measure of the star formation rate) depends on the choice of aperture and spatial scale of the observations. After describing how we estimate H₂ masses and the recent star formation rate (Section 2) and outlining our methodology (Section 3), we show the dependence of the depletion times on spatial scale and region targeted (Section 4). We then explore physical explanations for these results (Section 5).

We require the distributions of H₂ and recently formed stars which we trace via CO emission and a combination of Hα and IR emission, respectively.

2.1. Molecular Gas from CO Data

Star-forming clouds consist mainly of H₂, which cannot be directly observed under typical conditions. Instead, H₂ is usually traced via emission from the second most common molecule, CO. We follow this approach, estimating H₂ masses from the CO J = 1 − 0 data of Rosolowsky et al. (2007), which combines the BIMA (interferometric) data of Engargiola et al. (2003) and the FCRAO 14 m (single-dish) data of Heyer et al. (2004). The resolution of the merged data cube is 13'' × 2.03 km s⁻¹ with a median 1σ noise of 240 mK (∼2.1 M_☉ pc⁻² for our adopted X_CO). Rosolowsky et al. (2007) showed that this combined cube recovers the flux of the Heyer et al. (2004) FCRAO data.

We convert integrated CO intensities into molecular gas surface densities assuming X_CO = 2.0 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹. This is approximately the Milky Way conversion factor and agrees well with work on M33 by Rosolowsky et al. (2003). For this X_CO,

$$\begin{equation} \Sigma _{\rm H2} [ {M}_\odot \;{\rm pc}^{-2}] = 4.4 I_{\rm CO} [ {\rm K}\; {\rm km \;s}^{-1}], \end{equation} \tag{ 1 }$$

where I_CO is the integrated CO intensity over the line of sight and Σ_H2 is the mass surface density of molecular gas, including helium.

The data cover a wide bandpass, only a small portion of which contains the CO line. As a result, direct integration of the cube over all velocities produces an unnecessarily noisy map. Therefore, we "mask" the data, identifying the velocity range likely to contain the CO line along each line of sight. We integrate over all channels with ±25 km s⁻¹ of the local mean H i velocity (using the data from Deul & van der Hulst 1987). To ensure that this does not miss any significant emission, we also convolve the original CO cube to 30'' resolution and then identify all regions above 3σ in 2 consecutive channels. Any region within or near such a region is also included in the mask. We blank all parts of the data cube that do not meet either criteria and then integrate along the velocity axis to produce an integrated CO intensity map. Figure 1 shows this map at full resolution (middle left) and smoothed to ∼45'' resolution (top left) to increase the signal-to-noise ratio (S/R) and highlight extended emission. The noise in the integrated intensity map varies with position but typical 1σ values are 8–10 M_☉ pc⁻²; the dynamic range (peak S/R) is ∼20. The 3σ mass sensitivity in an individual resolution element is ∼10⁵ M_☉.

Zoom In Zoom Out Reset image size

2.2. Recent Star Formation from Hα and IR Data

We trace the distribution of recent star formation using Hα emission, which is driven by ionizing photons produced almost exclusively in very young (massive) stars. We account for extinction by combining Hα and infrared (24 μm) emission, a powerful technique demonstrated by Calzetti et al. (2007) and Kennicutt et al. (2007, 2009). Assuming continuous star formation over the past 100 Myr and studying a set of extragalactic star-forming regions, Calzetti et al. (2007) found the recent SFR to be

$$\begin{equation} {\rm SFR} \left[ {M}_\odot \;{\rm yr}^{-1} \right] = 5.3 \times 10^{-42} \left[ L \left({\rm H}\alpha \right) + 0.031\ L \left(24 \;\mu {\rm m}\right) \right], \end{equation} \tag{ 2 }$$

where L(Hα) and L(24 μm) = νL_ν(24 μm) are the luminosities of a region in Hα emission and at 24 μm, measured in erg s⁻¹.

The assumption of continuous star formation is certainly inapplicable to individual regions, which are better described by instantaneous bursts (e.g., Relaño & Kennicutt 2009). We only report averages of a large set (∼150) of regions, which together constitute a large part of M33's total Hα, and argue that this justifies the application of Equation (2) (see Section 3.2). In any case, SFR units allow ready comparison to previous work.

2.2.1. Hα Data

We use the narrowband Hα image obtained by Greenawalt (1998) with the KPNO 0.6 m telescope. The reduction, continuum subtraction, and other details of these data are described by Hoopes & Walterbos (2000). Before combining with the IR map, we correct the Hα map for Galactic extinction using a reddening of E(B − V) = 0.042 (Schlegel et al. 1998) and a ratio of Hα narrowband extinction to reddening of R(Hα) = 2.33.

Studies of M33 and other nearby galaxies find typically ∼40% of the Hα emission to come from "diffuse ionized gas" (DIG; Hoopes & Walterbos 2000; Thilker et al. 2002; Hoopes & Walterbos 2003; Thilker et al. 2005). The origin of this emission is still debated; it may be powered by leaked photons from bright H ii regions or it may arise "in situ" from isolated massive stars. We choose to remove this diffuse component from the Hα map and any discussion of Hα emission in the following analysis refers to this DIG-subtracted map (we assess the impact of this step in the Appendix). We do so using the following method from Greenawalt (1998). We begin by median filtering the Hα map with a 900 pc kernel. We then identify H ii regions as areas in the original map that exceed the median-filtered map by an emission measure of 50 pc cm⁻⁶ (outlined by black contours in the top right panel of Figure 1). We blank these regions in the original map and smooth to get an estimate of DIG emission toward the H ii regions. Our working map consists of only emission from the H ii regions after the DIG foreground has been subtracted. The integrated Hα flux (inside R_gal = 4.5 kpc) allocated to the diffuse map is 1.05 × 10⁴⁰ erg s⁻¹ (44%) while the part allocated to the DIG-subtracted, H ii-region map is 1.35 × 10⁴⁰ erg s⁻¹ (56%), in good agreement with previous results on M33 and other nearby galaxies.

2.2.2. IR Data

We measure IR intensities from 24 μm maps obtained by the Spitzer Space Telescope (PI: Gehrz et al. 2005, see also Verley et al. 2007). The data were reduced by K. Gordon (2009, private communication) following Gordon et al. (2005). Spitzer's point spread function at 24 μm is ∼6'', well below our smallest aperture size (∼18'') and so is not a large concern.

As with the Hα image, the 24 μm map includes a substantial fraction of diffuse emission—infrared cirrus heated by an older population, emission from low-mass star-forming regions, and dust heated leakage from nearby H ii regions, with a minor contribution from photospheric emission of old stars. Verley et al. (2007) argue that this diffuse emission accounts for ∼2/3 of all 24 μm emission in M33. To isolate 24 μm emission originating directly from H ii regions, we follow a similar approach to that used to remove DIG from the Hα map. The key difference is that instead of trying to identify all 24 μm bright sources by filtering and applying a cutoff to the 24 μm emission, we use the existing locations of H ii regions to isolate any local 24 μm excess associated with H ii regions. We extinction-correct the DIG-subtracted Hα emission using only this local excess in 24 μm emission. The total integrated flux at 24 μm (R_gal ⩽ 4.5 kpc) is 3.92 × 10⁴¹ erg s⁻¹, the fraction of DIG-subtracted 24 μm inside the H ii region mask is 1.63 × 10⁴¹ erg s⁻¹ (42%). The 24 μm correction implies Hα extinctions of A_Hα ∼ 0.3–0.4 mag.

To quantify the scale dependence of the molecular star formation law, we measure the H₂ depletion time,⁵τ_dep = Σ_H2/Σ_SFR, for apertures centered on bright CO and Hα peaks. We treat the two types of peaks separately and vary the sizes of the apertures used. In this way, we simulate a continuum of observations ranging from nearly an entire galaxy (>1 kpc apertures) to studies of (almost) individual GMCs or H ii regions (75 pc apertures). The CO data limit this analysis to galactocentric radii smaller than 4.5 kpc (∼0.6 r₂₅).

3.1. Identifying CO and Hα Peaks

3.2. Measuring Depletion Times

3.3. Uncertainties

Figure 2 shows a well-known result for M33. There is a strong correlation (rank correlation coefficient of r ≈ 0.8) between the surface densities of SFR and H₂ at 1200 pc scales. Power law fits to the different samples (types of peaks) and Monte Carlo iterations yield H₂ depletion times, τ_dep = M_H2/SFR, of ∼1 Gyr (at H₂ surface densities of 3 M_☉ pc⁻²) and power law indices of ∼1.1–1.5. These results (modulo some renormalization due to different assumptions) match those of Heyer et al. (2004) and Verley et al. (2010) in their studies of the star formation law in M33. The important point here is that there is good evidence for an internal H₂–SFR surface density relation in M33.

Zoom In Zoom Out Reset image size

We plot the median τ_dep as a function of scale (aperture size) in Figure 3, giving results for apertures centered on CO (red circles) and Hα (blue stars) peaks separately. For the largest scales, we find a similar τ_dep for both sets of apertures (as was evident from Figure 2). Going to smaller aperture sizes, τ_dep becomes a strong function of scale and type of peak targeted. Small apertures centered on CO peaks have very long τ_dep (up to 10 Gyr). Small apertures targeted toward Hα peaks have very short τ_dep (0.3 Gyr). This may not be surprising, given the expectations that we outlined in Section 1 and the distinctness of the bright Hα and CO distributions seen in the lower left panel of Figure 1, but the dramatic difference as one goes from ∼kpc to ∼100 pc scales is nonetheless striking.

Zoom In Zoom Out Reset image size

A few caveats apply to Figure 3. First, in subtracting the diffuse emission (DIG) from the Hα map, we removed ∼40% of the flux. This could easily include faint regions associated with CO peaks, which instead show up as zeros in our map. Perhaps more importantly, we use the 24 μm map only to correct the DIG-subtracted Hα map for extinction. Any completely embedded star formation will therefore be missed. For both of these reasons, the SFR associated with the red points, while it represents our best guess, may be biased somewhat low and certainly reflects emission from relatively evolved regions—those regions that have Hα fluxes above our DIG-cutoff value. There is no similar effect for the CO map.

Figure 3 implies that there is substantial movement of points in the star formation law parameter space as we zoom in to higher resolution on one set of peaks or another. Figure 4 shows this behavior, plotting the median Σ_SFR and median Σ_H2 for each set of apertures (N.B., the ratio of median Σ_H2 to median Σ_SFR does not have to be identical to the median τ_dep; the difference is usually ≲30%). We plot only medians because individual data are extremely uncertain, include many upper limits, and because we are primarily interested in the systematic effects of resolution on data in this parameter space.

Zoom In Zoom Out Reset image size

Apertures centered on CO peaks (red points) have approximately constant Σ_SFR, regardless of resolution. This can be explained if emission in the Hα map is homogeneously distributed as compared to the position of CO peaks. Meanwhile, there is a strong change in Σ_H2 for decreasing aperture sizes on the same peaks; Σ_H2 goes up as the bright peak centered on fills more and more of the aperture. A similar effect can be seen for the Hα (blue stars), though there is more evolution in Σ_H2 with increasing resolution because most bright Hα peaks also show some excess in CO emission.

Figure 3 shows that by zooming in on an individual star-forming region, one loses the ability to recover the star formation law observed on large scales. For apertures ≲300 pc in size, the relative amounts of CO emission and Hα intensity vary systematically as a function of scale and what type of region one focuses on. Another simple way to put this, demonstrated in Figure 4, is that scatter orthogonal to the SFR–H₂ relation increases with increasing resolution. Eventually this washes out the scaling seen on large scales and the star formation law may be said to "break down."

What is the origin of this scale dependence? In principle, one can imagine at least six sources of scale dependence in the star formation law.

• 1.  
  Statistical fluctuations due to noise in the maps.
• 2.  
  Feedback effects of stars on their parent clouds.
• 3.  
  Drift of young stars from their parent clouds.
• 4.  
  Region-to-region variations in the efficiency of star formation.
• 5.  
  Time evolution of individual regions.
• 6.  
  Region-to-region variations in how observables map to physical quantities.

Our observations are unlikely to be driven by any of the first three effects. In principle, statistical fluctuations could drive the identification of Hα and CO peaks leading to a signal similar to Figure 3 purely from noise. However, our Monte Carlo calculations, the overall S/R in the maps, and the match to previous region identifications make it clear that this is not the case.

Photoionization by young stars can produce CO shells around H ii regions inside a larger cloud or complex. This is a clear case of a small-scale offset between Hα and CO. However, the physical scales in Figure 3 are too large for this effect to have much impact; it should occur on scales more like ∼10 pc.

Similarly, the scales over which τ_dep diverges between CO and Hα peaks (75–300 pc) are probably too large to be produced by drift between young stars and their parent cloud. A typical internal GMC velocity dispersion in M33 is a few km s⁻¹ (1σ; Rosolowsky et al. 2003). Over an average cloud lifetime (∼30 Myr; Blitz et al. 2007; Kawamura et al. 2009), this implies a drift of at most 100 pc. This extreme case is just large enough to register in our plot but unlikely to drive the signal we see at scales of 150–300 pc. See Engargiola et al. (2003) for a similar consideration of GMCs and H i filaments.

Instead of drifts or offsets, what we observe is simply a lack of direct correspondence between the CO and Hα luminosities of individual star-forming regions. The brightest 150 CO peaks are simply not identical to the brightest 150 Hα peaks. The bottom right panel of Figure 1 shows this clearly; about a third of the peaks in M33 are nearer to another peak of their own type (i.e., CO to CO or Hα to Hα) than to a peak of the other type. Thus, Figure 3 shows that the ratio of CO to Hα emission varies dramatically among star-forming regions. In this case, the size scale on the x-axis in Figure 3 is actually a proxy for the number of regions inside the aperture. In M33, apertures of 75 pc diameter usually contain a single peak. At 150 pc, this is still the case ∼70% of the time, and at 300 pc only a few regions are included in each aperture.

Why does the ratio of CO-to-Hα vary so strongly from region to region? The efficiency with which gas forms stars may vary systematically from region to region (with high Hα peaks being high-efficiency regions), star-forming regions may undergo dramatic changes in their properties as they evolve (with Hα peaks being evolved regions), or the mapping of observables to physical quantities (Equations (1) and (2)) may vary from region to region.

It is difficult to rule out region-to-region efficiency variations, but there is also no strong evidence for them. Leroy et al. (2008) looked for systematic variations in τ_dep as a function of a number of environmental factors and found little evidence for any systematic trends. Krumholz & McKee (2005) and Krumholz & Tan (2007) suggested that the cloud free-fall time determines τ_dep to first order, but based on Rosolowsky et al. (2003), the dynamic range in free-fall times for M33 clouds is low. On the other hand, Gardan et al. (2007) found unusually low values of τ_dep in the outer disk of M33.

There is strong evidence for evolution of star-forming regions. Fukui et al. (1999), Blitz et al. (2007), Kawamura et al. (2009), and Chen et al. (2010) showed that in the LMC, the amount of Hα and young stars associated with a GMC evolves significantly across its lifetime. In our opinion this is the most likely explanation for the behavior in Figure 3. Star-forming regions undergo a very strong evolution from quiescent cloud, to cloud being destroyed by H ii region, to exposed cluster or association. When an aperture contains only a few regions, the τ_dep for that aperture will be set by the evolutionary state of the regions inside it. That state will in turn determine whether the aperture is identified as a CO peak or an Hα peak. CO peaks will preferentially select sites of heavily embedded or future star formation while Hα peaks are relatively old regions that formed massive stars a few Myr ago.

Region-to-region variations in the mapping of observables (CO and Hα) to physical quantities (H₂ mass and SFR) are expected. Let us assume for the moment that the ratio of H₂ to SFR is constant and independent of scale. Then to explain the strong scale dependence of the ratio of CO to Hα in Figure 3 there would need to be much more H₂ per unit CO near Hα peaks and many more recently formed stars per ionizing photon near the CO peaks. At least some of these effects have been claimed: e.g., Israel (1997) find a strong dependence of X_CO on radiation field and Verley et al. (2010) suggest that incomplete sampling of the initial mass function in regions with low SFRs drive the differences they observe between star formation tracers. However, both claims are controversial and it seems very contrived to invoke a scenario where only this effect drives the breakdown in Figure 3. It seems more plausible that the mapping of observables to physical quantities represents a secondary source of scatter correlated with the evolutionary state of a region (e.g., the age of the stellar population).

5.1. Comparison to a Simple Model

5.2. τ_dep Versus Scale at Different Radii

Our main conclusion is that the molecular star formation law observed in M33 at large scales (e.g., Heyer et al. 2004; Verley et al. 2010) shows substantial scale dependence if one focuses on either CO or Hα peaks. The median depletion time (or CO-to-Hα ratio) measured in a 75 pc diameter aperture (derived from averaging ∼150 such apertures) varies by more than an order of magnitude between CO and Hα peaks. At large (∼kpc) scales, this difference mostly vanishes. We argue that the scale for the breakdown is set by the spatial separation of high-mass star-forming regions, with the breakdown occurring when an aperture includes only a few such regions in specific evolutionary states (a scale that corresponds to ∼300 pc in M33).

In this case, the scaling relation between gas and star formation rate surface density observed at large scales does not have its direct origin in an instantaneous cloud–scale relation. Individual GMCs and H ii regions will exhibit a CO-to-Hα ratio that depends on their evolutionary state (likely with significant additional stochasticity) and as a result the ∼150 brightest objects at a given wavelength will be a function of the evolutionary state that observation probes. This divergence is consistent with recent results from the LMC (Kawamura et al. 2009 and Chen et al. 2010) indicating that individual GMCs exhibit a range of evolutionary states over their 20–30 Myr lifetime.

This does not mean that comparisons of tracers of recent and future star formation on small scales are useless. On the contrary, such observations contain critical information about the evolution of individual clouds as a function of time and location that is washed out at large scales (≳300 pc in M33). However, once one moves into the regime where a single object contributes heavily to each measurement, it is critical to interpret the results in light of the evolution of individual clouds.

We thank Mark Heyer and Edvige Corbelli for providing us with the CO data, Rene Walterbos for the Hα image, and Karl Gordon for the Spitzer images of M33. We thank the anonymous referee for a careful and constructive report that improved the paper. We acknowledge the use of the H ii region catalog from Hodge et al. A.K.L. thanks Michele Thornley and Jack Gallimore for helpful discussions. The work of A. S. was supported by the Deutsche Forschungsmeinschaft (DFG) Priority Program 1177. This research made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the JPL/Caltech, under contract with NASA as well as NASA's Astrophysical Data System (ADS). Support for A.K.L. was provided by NASA through Hubble Fellowship grant HST-HF-51258.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.

Here we test how our method of selecting peaks and the removal of a diffuse ionized component affect our results.

First, we repeated the analysis on the original maps without any DIG subtraction. We show the results in the left panel of Figure 5, along with the original measurements (in gray) from Figure 3. To first order, the scale dependence of τ_dep is unchanged, but the H₂ depletion times are offset; τ_dep derived from maps with no DIG subtraction is a factor of ∼2 smaller from the DIG-subtracted maps.

Zoom In Zoom Out Reset image size

As a second test, we assess the impact of the particular choice of aperture positions. In the main text, we used a peak-finding algorithm. Here we test the effect of using published positions of GMCs and H ii regions instead. The right panel in Figure 5 shows our original data in gray while the red and blue symbols are derived using the Rosolowsky et al. (2007) and Hodge et al. (2002) catalogs. The median τ_dep at large scales is unchanged from Figure 3 (gray values). However, for apertures centered on GMCs, τ_dep on small scales does change from our analysis. This difference originates in different numbers and locations of the positions that are studied. In our original analysis, we study 172 positions which have CO emission peaks above 3σ. The Rosolowsky catalog, on the other hand, consists of only 140 positions inside a galactocentric radius of 4.5 kpc. In addition, a subset of the two samples targets different regions in M33: first, the catalog positions tend to be more clustered than the "peak" positions which leads to a somewhat larger number (5%–15%) of objects in the smaller apertures and a smaller deviation in depletion times for CO and Hα centered apertures. Second, while the molecular gas surface densities at the positions of the two samples do not differ significantly, the star formation rate surface densities are a factor of ∼3 higher for the catalog positions as compared to the (more numerous) "peak" positions. This leads to shorter H₂ depletion times on small scales for the catalog sample.

Both tests show that the analyzed scale dependence of the star formation relation and the determination of its origin is not strongly dependent on the particular methodology chosen in this paper. While global shifts in the derived depletion times can arise due to the subtraction of diffuse Hα emission, we find only small variations in the scale dependence due to different selection of positions where we perform our study.