THE VIRUS-P EXPLORATION OF NEARBY GALAXIES (VENGA): THE X_CO GRADIENT IN NGC 628

All the measurements of nebular emission lines used to estimate the metallicity, ionization parameter, dust extinction through the Balmer decrement, and the Hα SFR, are made on the Mitchell Spectrograph (formerly VIRUS-P) integral field unit (IFU) data cube of NGC 628 produced by the VENGA survey (Blanc et al. 2010). The data-cube samples a rectangular area of 52 × 17 centered in the nucleus of the galaxy, and has a spatial resolution of 56 FWHM. The spectra cover the 3555–6790 Å wavelength range with an instrumental spectral resolution of ≃ 5 Å FWHM. At the assumed distance of NGC 628 the VENGA spatial resolution corresponds to ∼235 pc, a few times larger than the typical sizes of large GMCs (≲ 60 pc; Heyer et al. 2001) and individual H ii regions ionized by single clusters (≲ 120 pc; Gutiérrez et al. 2011). The observations, data reduction, and calibrations will be described in an upcoming publication.

To ensure that our results are robust against the choice of the SFR tracer used, along with the VENGA Hα derived SFRs, we use a linear combination of far-ultraviolet (FUV) and mid-infrared (mid-IR) 24 μm emission. The former traces unobscured star formation, while the latter recovers dust obscured star formation by means of reprocessed UV radiation reemitted as thermal IR emission from heated interstellar dust grains. A number of calibrations have been proposed to use linear combinations of obscured and unobscured tracers to estimate the SFR (Calzetti et al. 2007; Leroy et al. 2008, 2012; Kennicutt et al. 2009; Hao et al. 2011; Murphy et al. 2011).

We use the Galaxy Evolution Explorer (GALEX) FUV image of NGC 628 from the GALEX Nearby Galaxy Survey (Gil de Paz et al. 2007) which samples the 1350–1750 Å wavelength range, has a point-spread function (PSF) FWHM of 45, and is deep enough for us to measure the FUV flux at high signal-to-noise ratio over the whole area of interest. To trace the IR dust emission we use the MIPS 24 μm data of NGC 628, taken as part of the Spitzer Infrared Nearby Galaxy Survey (SINGS; Kennicutt et al. 2003) and the Local Volume Legacy survey (Dale et al. 2009). The map has a PSF FWHM of 6'' and like the FUV map, it is deep enough for us to measure the 24 μm flux over the whole area of interest in this paper. The data processing of both the FUV and 24 μm maps, including the suppression of backgrounds and the masking of foreground stars is detailed in Leroy et al. (2012).

We use three independent data sets to measure I(CO) across the disk of NGC 628. The Heterodyne Receiver Array CO Line Extragalactic Survey (HERACLES; Leroy et al. 2009) CO(2–1) map, obtained at the IRAM 30 m single-dish telescope, has a beam size of 136 and rms noise of ∼22 mK per 2.6 km s⁻¹ channel, which translates into a 1σ limit on the molecular gas surface density of $\Sigma _{{\rm H_2}}\sim 3$ M_☉ pc⁻². Following Schruba et al. (2011) we assume a constant CO(2–1)/CO(1–0) line ratio of 0.7 to estimate I(CO) from the HERACLES data.

The Berkeley–Illinois–Maryland Array (BIMA) Survey of Nearby Galaxies (BIMA-SONG; Helfer et al. 2003) CO(1–0) map combines zero spacing single-dish data from the NRAO 12 m telescope and interferometric BIMA C and D array data, resulting in a map with a robust beam size of 62, and rms noise of 51 mJy beam⁻¹ in a 10 km s⁻¹ channel or $\Sigma _{{\rm H_2}} \sim 13$ M_☉ pc⁻².

Finally, the CARMA CO(1–0) map of NGC 628 (Rahman et al. 2012) has a spatial resolution of 36 and rms noise of ∼20 mJy beam⁻¹ in a 10 km s⁻¹ channel ($\Sigma _{{\rm H_2}}\sim 5$M_☉ pc⁻²). We convolve the BIMA-SONG and CARMA maps with a Gaussian kernel to match the 136 PSF of the HERACLES map, which corresponds to ∼570 pc at the assumed distance to NGC 0628. The three PSF-matched CO maps are shown in Figure 1.

Zoom In Zoom Out Reset image size

In order to account for systematic errors in the absolute flux calibration of the data cubes, and in the assumed CO(2–1) to CO(1–0) conversion factor, we scale the maps to match the total CO(1–0) luminosity of the HERACLES map in the inner 5 kpc of the galaxy. We find scaling factors of 1.36 and 1.36 for the BIMA-SONG and CARMA maps, respectively. Figure 2 presents both the scaled and original flux measurements in each of the radial bins over which we will measure X_CO (see Section 3.2) for the three maps. The ∼30% difference implied by these scaling factors is consistent with the systematic uncertainties in the absolute flux calibration of the three data sets, and the observed scatter in the CO(2–1) to CO(1–0) ratio across the HERACLES galaxies (Leroy et al. 2009). Assuming that the CO(2–1) to CO(1–0) ratio has no significant radial dependence, these corrections only introduce a scaling in the derived X_CO factors. Since we are interested in measuring relative changes in X_CO within the galaxy, a scaling of this type does not affect our results.

Zoom In Zoom Out Reset image size

The radial dependence in the CO(2–1) to CO(1–0) ratio across the disk of spiral galaxies has been studied by Leroy et al. (2009) by comparing the CO(2–1) emission in the HERACLES maps to the Nobeyama 45 m single-dish CO(1–0) maps of (Kuno et al. 2007; see Figure 34 in Leroy et al. 2009). Although the scatter seen in the line ratio is large, no evidence for a radial dependence is seen outside the very central parts of the galaxies (R < 0.05 R₂₅). Furthermore, as will be discussed in Section 4, the relatively good consistency between the different data sets (see also Figure 2), and the agreement between our measured X_CO radial profiles and independent measurements based on dust modeling indicate that any potential radial trend in the CO(2–1) to CO(1–0) ratio is at a significantly smaller level than the radial trends seen in X_CO.

The consistency between the different measurements shown in Figure 2 is reassuring. While the BIMA-SONG map includes both BIMA interferometric data and NOAO 12 m single-dish data, the CARMA map is constructed without the addition of zero spacing information, and therefore lacks sensitivity on large scales. The agreement seen in Figure 2 implies that even in the outermost radial bin for which we can measure the flux in the CARMA map, the lost large-scale extended emission does not contribute significantly to the total CO surface brightness.

We produce maps of Σ_SFR using two independent SFR indicators. First we use the Hα emission line flux, corrected for dust extinction using the Balmer decrement. The methods used to measure dust-corrected Hα fluxes and to correct the emission line map for the contribution from diffuse ionized gas are analog to the ones described in Blanc et al. (2009). An updated description of the methods used, including slight differences in the methodology used for flux calibration and construction of the IFU data cubes will be described in an upcoming publication. We use the Hα SFR calibration of Murphy et al. (2011) and Hao et al. (2011), taken from the compilation in Kennicutt & Evans (2012), which assumes a Kroupa & Weidner (2003) initial mass function (IMF). Note that this differs from the Salpeter (1955) IMF used in Blanc et al. (2009).

The second SFR indicator used is a linear combination of FUV and 24 μm flux. We use the calibration of Leroy et al. (2012), which assumes a Chabrier (2003) IMF. Calibrations based on the Kroupa and Chabrier IMFs yield nearly identical SFRs (Chomiuk & Povich 2011), which are typically ∼30% lower than those obtained assuming a Salpeter IMF. Since Hα and FUV photons trace star formation over different timescales (roughly 10 and 100 Myr, respectively), systematic differences between the two independent estimates might arise as a consequence of differences in star formation history. Analogously to what was done with the CO maps, we attempt to remove any systematic differences in the flux calibration of the different data sets and the SFR calibrations used by matching the total SFR within the central 5 kpc of the galaxy. This translates to us multiplying the Σ_{SFR, FUV + 24 μm} map by a factor of 0.8 to match the Σ_{SFR, Hα} map. Again, this correction does not affect the relative differences in X_CO which we are trying to measure, but only introduces a scaling of the derived values across the whole galaxy. Both maps are convolved with a Gaussian kernel in order to match the common 136 FWHM PSF, and are presented in Figure 1.

We use the PSF-matched SFR maps to measure Σ_{SFR, FUV + 24 μm}, Σ_{SFR, Hα}, and I(CO) in elliptical annuli of constant galactocentric radius. We chose an annuli width of 40'' which roughly corresponds to three times the spatial FWHM of the maps, and a physical scale of 1.68 kpc at the adopted distance to NGC 628 (see Figure 1). The sensitivity of the BIMA-SONG and CARMA maps only allows the measurement of I(CO) out to 5 kpc (∼120'') which corresponds to the inner three annuli, while the HERACLES map is able to provide a reliable measurement in a fourth bin, out to a galactocentric distance of 7 kpc (∼170''). While Schruba et al. (2011) was able to measure CO emission using the HERACLES map of NGC 628 out to larger radii (∼10 kpc) using a stacking technique, here we are also limited by the coverage of the VENGA IFU observations which only reach out to ∼7.5 kpc. We use the extended stacking measurements of Schruba et al. (2011) separately in the next section to extend our measurements to larger radii. Having measured these three quantities (Σ_{SFR, FUV + 24 μm}, Σ_{SFR, Hα}, and I(CO)), we use Equation (2) to obtain Σ_H2, and we transform it to units of column density. Finally we calculate X_CO for each radial bin by applying Equation (1). The central radius and measured X_CO values for each annuli are reported in Table 1.

The uncertainty in the measured X_CO values is dominated by the intrinsic scatter in the molecular SFL. A series of studies have shown that this quantity depends on the physical scales over which the observed surface densities are integrated (Verley et al. 2010; Onodera et al. 2010; Schruba et al. 2011; Liu et al. 2011). For the large physical scales considered in this work (∼4 kpc given the mean area of our radial bins), we assume a value of ∼0.15 dex for the scatter, which we propagate into the uncertainty in X_CO. Thanks to the high signal-to-noise ratio of the data sets and the large areas over which we are integrating, the contribution from random errors in Σ_SFR, which includes photometric errors in the Hα and Hβ (the latter entering in the Balmer decrement dust extinction correction), and in the FUV and 24 μm fluxes, as well as from noise in the CO maps, is negligible (a factor of ∼5 smaller) compared to the scatter in the molecular SFL.

In Figure 3 we present the measured X_CO factor as a function of galactocentric radius for the six combinations of data sets (two SFR indicators and three CO maps from different telescopes), and two different assumptions for the molecular SFL slope (N = 1.0 and 1.5). At all radii, all combinations of data sets are consistent with each other within the 1σ formal error bars, and show the same systematic trend of increasing X_CO with radius. This implies that our results are robust against potential systematics associated with the adopted SFR indicator, and the nature of the CO maps utilized (single-dish versus interferometer, CO(1–0) versus CO(2–1)). The consistency between different data sets also confirms that photometric errors are not a significant source of uncertainty, when compared to the systematic errors associated with the scatter in the SFL.

Zoom In Zoom Out Reset image size

As can be seen in Figure 1, the agreement between the Hα and FUV+24 μm SFR maps is good. After performing the scaling described in Section 3.1, the ratio of the two different SFRs measured in the inner four radial bins shows an rms scatter of 13%. When analyzing the morphology of these two SFR tracers we note a strong similarity. All individual star-forming regions picked up by one method are also recovered by the other. This is important as it implies that the dust-corrected Hα method is not missing a significant fraction of completely obscured star formation which should be recovered by dust emission at 24 μm.

While both maps trace the same structures, a more detailed inspection shows that the FUV+24 μm SFR map does seems to have a more diffuse morphology than the Hα SFR map. In particular more diffuse FUV+24 μm emission is seen in the inter-arm regions of the galaxy, while the Hα emission seems to be more concentrated in the arms. This is consistent with the typical star formation timescales traced by these two methods. Hα traces young stars (⩽10 Myr), which are seen close to their birthplaces, while both FUV and 24 μm emission receive a significant contribution from older stars (⩽100 Myr) which are more homogeneously distributed across the galactic disk, and not necessarily associated with the spiral arms.

The morphological agreement between the different CO maps is also very good. After applying the scalings described in Section 2.3 the CO luminosities in the inner three radial bins show a typical rms scatter of 20% (see Figure 2). The agreement justifies our adoption of a single CO (2–1)/(1–0) line ratio across the galaxy. It is interesting to note that while both the CO emission and the SFR trace similar morphological structures across the galaxy, that is, mainly the grand design spiral arms of NGC 628, a detailed comparison shows that there is not a one-to-one correlation between individual CO bright regions and regions which are bright in Hα or FUV+24 μm. This reflects the temporal dependence of the star formation process, in which recently formed star clusters disrupt and dissociate their molecular gas birth clouds by injecting both radiative and mechanical energy into their surrounding ISM. As this process takes place, an individual star-forming region evolves from a CO bright molecular complex into a 24 μm bright, highly obscured, star-forming region, which eventually dissociates and removes the bulk of the molecular gas and dust surrounding it, becoming an Hα and FUV bright region in the process. This temporal evolution translates into a spatial offset between SFR and molecular gas tracers, and is the main cause behind the need to average over large kiloparsec scales in order to do a measurement like the one presented in this work (Schruba et al. 2010; Madore 2010; Feldmann et al. 2012b; Kennicutt & Evans 2012).

In order to quantify the previous statements regarding the robustness of our measurement, we calculate the standard deviation of all X_CO measurements in each radial bin, and find it to range from a few percent to ∼25% of the measured values in the worst case (outermost bin). The mean rms across all four bins corresponds to ∼15% of the measured values. This dispersion is caused by a combination of different factors including systematic errors in the adopted CO(2–1)/CO(1–0) ratio, a possible systematic lack of sensitivity to low surface brightness emission in the interferometric CO maps, and potential variations in the ratio between the Hα and FUV+24 μm SFRs which might arise as the consequence of changes in the stellar populations and the properties of dust across the disk of the galaxy. The magnitude of these systematic deviations is smaller than the formal 0.15 dex (∼35%) uncertainty adopted. Nevertheless, we add this dispersion in quadrature to our formal errors in order to take into account the effects of theses systematics in our error budget.

The dotted horizontal line in Figure 3 marks the canonical MW X_CO factor of 2.3 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹ from Pineda et al. (2010), which is adopted as the preferred value in Kennicutt & Evans (2012). While at radii >2 kpc the X_CO factor in NGC 628 is consistent with the MW value, there is an evident drop in the central 2 kpc of the galaxy, where the measured X_CO is a factor of two lower than the MW value.

The drop in the radial profile of X_CO toward the central regions is observed independently of the assumed SFL slope. A nonlinear SFL with N > 1 implies lower H₂ surface densities, and therefore lower derived X_CO values for high SFR surface density regions than a linear relation. Since the SFR in NGC 628, and in most disk galaxies, decreases with radius, a nonlinear slope implies an even larger drop in X_CO toward the central regions. In any case, the magnitude of the changes in X_CO between the N = 1.0 and 1.5 cases is smaller than the changes in X_CO observed as a function of radius across the galaxy. This implies that potential changes in the SFL slope within the galaxy would not affect our results significantly. While a sub-linear (N < 1) slope could reduce the magnitude of the observed drop in X_CO toward the central regions, this scenario is not currently supported by observations (Calzetti et al. 2012; Kennicutt & Evans 2012, and references within).

This dependence with the SFL slope is important, as it implies that the observed radial trend of X_CO is not a consequence of changes in the SFE with radius. Adopting a super-linear slope for the SFL is effectively testing the case in which the efficiency of star formation is enhanced in higher density regions, and as stated above, this effect only boosts the decrease of X_CO at small radii, although not significantly. We further discuss this subject in the next section.

Also shown as open circles in Figure 3 are the measurements of X_CO in NGC 628 from dust spectral energy distribution (SED) modeling (Sandstrom et al. 2012). We only plot their data points with an uncertainty ⩽0.4 dex in X_CO. This study uses the method of Leroy et al. (2011) to simultaneously model the dust mass surface density and the dust-to-gas ratio across the disks of a sample of nearby spiral galaxies with far-IR SED measurements from the SINGS and KINGFISH (Kennicutt et al. 2011) surveys. The dust mass surface density, in combination with the dust-to-gas ratio, provides an estimate of the total gas surface density. Subtraction of the atomic component (measured from H i 21 cm maps) yields the surface density of the molecular component. This technique is completely independent of the method adopted in this work. The data points from the dust SED modeling have a typical uncertainty of 0.2 dex at R < 7 kpc, hence they agree with our measurements across the whole range in radii sampled by our data, and significantly show the same trend of decreasing X_CO toward the central regions of NGC 628. No scaling has been applied to the measurements in Figure 3, so the agreement is not only in the shape of the X_CO radial profile but also in its absolute value.

A similar analysis was conducted by Aniano et al. (2012) for NGC 628 and NGC 6946. While a drop in X_CO toward the center of the galaxy is evidently seen in their data for NGC 6946, the authors do not claim the detection of a changing X_CO profile in NGC 628. Inspecting their Figure 4, it is evident that the results of their measurements of X_CO in NGC 628 are subject to strong systematic uncertainties, and show large changes (from an increasing to a flat, or even decreasing, radial trend) depending on which subset of the Spitzer+Herschel data is used and the spatial resolution at which the modeling is conducted. This is not the case for NGC 6946 where different data sets at different resolutions yield consistent results. Given the large systematic uncertainties in their measurement of X_CO for NGC 628, we refrain from conducting a detailed comparison to their results.

In Section 3, we discussed the inherit degeneracy in our method between the SFL normalization constant (or equivalently the gas depletion time when N = 1), and the measured values of X_CO. A basic assumption in this work is that of a constant SFL across the disk of NGC 628. In order to explain the data in Figure 3 solely as a variation in the normalization of the SFL would require the depletion timescale to be longer, or the SFE to be lower, by a factor of two in the central (denser) regions than in the outer disk. This is not in line with expectations from measurements of the SFE across different environments in galaxies in the local universe. A series of studies have shown that luminous infrared galaxies, ULIRGS, and SMGs, in which the molecular gas surface density is one to two orders of magnitude larger than in the disks of normal spirals, show typical depletion timescales which are factors of 3–4 shorter than normal spiral galaxies. The difference becomes even larger (factors of 4–10 shorter) if differences in the X_CO factor in these systems are taken into account (Daddi et al. 2010; Genzel et al. 2010; García-Burillo et al. 2012). Although the central region (<2 kpc) of NGC 628 is only factors of a few denser than the outer disk, it would be surprising to find a reversal from the global trend seen in the SFE toward denser environments. Furthermore, the consistency between the results from this work and the independent dust SED modeling method, which is not subject to such a degeneracy between X_CO and the depletion timescale, is very encouraging, and gives us confidence that our working assumption is valid within the range of current uncertainties.

If we model the X_CO profile in the inner 7 kpc of the galaxy as a simple gradient of the form:

$$\begin{equation} {\rm log}(X_{\rm CO})={\rm log}(X_{{\rm CO},0})+\Delta {\rm log}(X_{\rm CO})\times R \end{equation} \tag{ 3 }$$

a fit to the average of the six measurements in each radial bin (two measurements in the outer bin), which are shown as stars in Figure 3, yields X_CO gradients of Δlog(X_CO) = 0.06 ± 0.02 dex kpc⁻¹ for N = 1, and Δlog(X_CO) = 0.08 ± 0.02 dex kpc⁻¹ for N = 1.5.

The best-fit central X_CO values are X_{CO, 0} = 1.3 ± 0.2 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹ for N = 1 and X_{CO, 0} = 1.1 ± 0.2 × 10²⁰ cm⁻¹ (K km s⁻¹)⁻¹ for the N = 1.5 case. Our data imply that X_CO in the galactic center of NGC 628 is about a factor of two higher than the typical values of X_{CO, 0} = 0.1–0.5 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹ measured in the MW Galactic center (e.g., Sodroski et al. 1995; Strong et al. 2004; Oka et al. 1998).

Inspecting Figure 3, particularly for the N = 1 case, it is not evident if the radial profile of X_CO follows a linear gradient, or if it is fairly flat across the disk of the galaxy, and only drops in the central regions. Our measurements are limited at R > 7 kpc by the sensitivity of the HERACLES CO map, and the coverage of the VENGA IFU data cube. The dust SED modeling data points go out to ∼10 kpc, but they show a very large scatter, and at such large radii, the uncertainty in their measurements increases sharply, so it is not easy to draw conclusions from these data points regarding the behavior of X_CO at large radii.

In an attempt to extend our measurements into the outer disk of the galaxy, we use our method to estimate X_CO on the stacked CO and FUV+24 μm radial profiles of NGC 628 presented in Schruba et al. (2011). In this work, the authors use the 21 cm H i velocity field to predict the exact wavelength of the CO line at every position in the HERACLES data cube. This information is then used to register and stack the spectra of all spatial resolution elements in 15'' wide radial bins across the galaxy. Using this technique, Schruba et al. (2011) is able to significantly measure I(CO) out to ∼10 kpc. They also measure the FUV+24 μm SFR profile using the same maps used in this work. We show these measurements as filled black circles in Figure 3. Error bars are not shown for clarity, but they are also dominated by the intrinsic scatter in the SFL, and are of the same order of magnitude as the error bars for the measurements in the inner regions of the galaxy (color filled circles and squares). The results of this exercise are inconclusive. While in the N = 1 case the value of X_CO seems to flatten around the canonical MW value, if we assume N = 1.5 X_CO such hints for a flattening in the outer regions are less obvious but still present.

As will be further discussed in Section 5, NGC 628 shows a metallicity gradient described by a fairly constant slope out to ∼1 R₂₅ (∼13 kpc; Rosales-Ortega et al. 2011). If X_CO scales as some power of the metal abundance (Z), we would expect it to keep on increasing out to a similar radius. On the other hand there are reasons why a flattening in the radial profile might be expected. For example, if X_CO increases with a decreasing molecular cloud surface density, as proposed by Narayanan et al. (2012), and the GMCs in the outer disk have uniform properties but become denser toward the central regions, a certain level of flattening in the X_CO profile might occur, although the metallicity dependence discussed above should still be present. We will explore the role of the molecular gas surface density and the UV radiation field in Section 6, and will discuss this possibility in more detail.

For all the results presented in the following sections, we adopt the N = 1.0 case. We have checked that changing the slope to N = 1.5 does not significantly affect the rest of the results presented below.

The VENGA IFU data cube provides maps of many strong nebular emission lines typically used for metallicity diagnostics. As discussed in detail in Kewley & Ellison (2008) and Moustakas et al. (2010), different strong-line methods used to measure metallicity can show large discrepancies in the derived abundances. Particularly, methods calibrated against theoretical photoionization models typically yield higher metallicities than methods calibrated against samples of H ii regions with direct electron temperature measurements (a.k.a. empirically calibrated methods). Differences as large as 0.6 dex are seen between these two families of methods. On the other hand, relative differences in metallicity from one system to another do not change drastically if different methods are applied, allowing meaningful conclusions to be drawn from comparing different samples, as long as a single method is used, or measurements using different methods are properly transformed to a common metallicity scale (Kewley & Ellison 2008). In order to study the metallicity dependence of X_CO across NGC 628, and compare it to previous measurements in the literature, we must not only choose a reliable strong-line abundance indicator, but we must also put ours and the literature measurements on a common metallicity scale.

Our preferred metallicity calibration is the N2O2 method of Kewley & Dopita (2002), which estimates the oxygen abundance from the [N i]λ6584/[O ii]λ3727 ratio, and is calibrated against photoionization models. We have a series of reasons to favor this calibration. Mainly, the N2O2 indicator is highly independent of the ionization parameter, is single valued at all metallicities, and the secondary production of nitrogen makes this indicator very sensitive to metallicity changes at high metallicities like the ones we expect in massive spiral galaxies like NGC 628 (see discussion in Kewley & Dopita 2002).

On the other hand, the metallicity compilations for the systems with measured X_CO values reported in Bolatto et al. (2008) and Leroy et al. (2011) are largely based on direct electron temperature measurements. The metallicities reported in Schruba et al. (2012) are on an intermediate scale, since the authors averaged the values from the two indicators used by Moustakas et al. (2010): the photoionization model calibrated R₂₃ metallicities of Kobulnicky & Kewley (2004, hereafter R₂₃ − KK04) and the empirically calibrated R₂₃ method of Pilyugin & Thuan (2005, hereafter R₂₃ − PT05).

We decide to apply a simple offset in metallicity to both our N2O2 measured values in NGC 628 and the reported metallicity values in Schruba et al. (2012), to put them on a common scale with the direct method metallicities reported in Bolatto et al. (2008) and Leroy et al. (2011). To do so, we measured the metallicity in the NGC 628 data cube using the three methods mentioned above (N2O2, R₂₃ − KK04, and R₂₃ − PT05), and measure the mean offset in the recovered metallicity values between different methods for all spaxels above an Hα surface brightness cut chosen to trace H ii regions.¹⁰ As expected, both the N2O2 and the average between the R₂₃ − KK04 and R₂₃ − PT05 methods (i.e., the method used in Schruba et al. 2012) yield higher metallicities than R₂₃ − PT05 alone, which is calibrated against direct method measurements on H ii regions and should fall in a similar scale to the values reported in Bolatto et al. (2008) and Leroy et al. (2011). We find and apply mean offsets of −0.66 dex and −0.35 dex to the N2O2 and Schruba et al. (2012) metallicities before doing any comparisons. These values are in good agreement with the 0.6 dex offset found between R₂₃ − KK04 and R₂₃ − PT05 by Moustakas et al. (2010).

Note that we decided to apply a single offset in metallicity instead of using the conversion relations of Kewley & Ellison (2008). This is done deliberately, as the nonlinear nature of these conversions imply that the measured metallicity gradient would change its value depending on the adopted scale. Also, these conversions have been calibrated against integrated measurements of galaxies, and it is not clear how valid it would be to apply them to spatially resolved regions within a single object. This issue will be the subject of a future publication in which we will analyze the impact of different methods on the observed metallicity distributions in the VENGA galaxies. For the purpose of this work, we prefer to conserve the shape of the gradient as measured by the N2O2 method which we consider more robust than other formulations. In any case, using the Kewley & Ellison (2008) conversions, or using R₂₃ − PT05 directly on our data instead of N2O2 does not significantly change our results.

Figure 4 presents the spatial distribution in metallicity of NGC 628. The top panel shows a two-dimensional map of 12+log(O/H), where thick contours mark regions above the Hα surface brightness cut where the nebular emission is dominated by H ii regions. The bottom panel shows the distribution of metallicity as a function of galactocentric radius for data-cube spaxels inside these regions. Black error bars show the median measurement errors in each radial bin, red error bars mark the standard deviation for all spaxels in each bin, and green error bars show the 0.1 dex intrinsic scatter associated with the N2O2 calibration (Kewley & Ellison 2008). We measure a gradient of Δlog(O/H) = 0.036 ± 0.002 dex kpc⁻¹, in good agreement with the IFU measurements of Rosales-Ortega et al. (2011).

Zoom In Zoom Out Reset image size

In Figure 5, red stars show the X_CO values measured in each radial annuli in NGC 628, as a function of the metallicity at the radius corresponding to the center of each annuli (evaluated from the metallicity gradient measured in Section 5.1). Also shown are the data of Bolatto et al. (2008), Leroy et al. (2011), and Schruba et al. (2012),¹¹ as well as predictions from the theoretical models of Narayanan et al. (2012) and Feldmann et al. (2012a).

Zoom In Zoom Out Reset image size

Across the disk of NGC 628, X_CO increases toward lower metallicity regions. This is expected from the positive and negative gradients measured for X_CO and the oxygen abundance, respectively (Sections 4 and 5.1). This behavior is consistent with integrated measurements of galaxies in the local universe. Both Leroy et al. (2011) and Schruba et al. (2012) see an increase in X_CO toward low metallicities by conducting integrated measurements across galaxies using dust mass modeling in the first case, and a similar method to the one used here in the later. The X_CO values inferred from GMC virial masses in Bolatto et al. (2008) do not show an obvious trend with metallicity, but are consistent with other studies given the large scatter in the observed relation. Also, as discussed in the original Bolatto et al. (2008) paper, and in Schruba et al. (2012), estimates of X_CO using virial mass measurements from CO observations most likely miss the "CO dark" envelopes of molecular clouds in low metallicity galaxies, as they can only trace mass enclosed in the CO emitting region. This could translate in an underestimation of X_CO at low metallicity in the Bolatto et al. (2008) data. The same effect could explain the constancy of X_CO seen by Rosolowsky et al. (2003) across 45 GMCs in M33, even in the presence of a 0.8 dex change in metallicity within their sample.

The NGC 628 data points are consistent with previous observations, especially considering the large scatter seen in the X_CO–metallicity relation (∼0.5 dex). The limited dynamic range in metallicity that we can sample across the inner 7 kpc of a single galaxy prevents us from being able to put strong constraints in the slope of the X_CO–metallicity relation. Modeling this relation as a power law of the form X_CO∝Z^α, Schruba et al. (2012) measures a slope of α = −2.0 ± 0.4 for the full HERACLES sample after rejecting starburst galaxies (green solid line in Figure 5), while Israel (1997) finds a slope of α = −2.7  ±  0.3 in their study of dust emission in local dwarfs. A linear fit to our data points yields a slope of α = −1.1 ± 1.6, shallower than previous measurements but still consistent with both Israel (1997) and Schruba et al. (2012). Again, given the limited dynamic range in metallicity of our data, and the large scatter in the X_CO–metallicity relation, detailed comparisons are difficult.

The slope of the observed relation in NGC 628 is also in agreement with expectations from theoretical models. At fixed molecular gas surface density, the simulations in Narayanan et al. (2012) predict a linear dependence of X_CO with metallicity (i.e., α = −1.0),¹² and the study by Feldmann et al. (2012a) predicts slopes in the −0.5 to −0.8 range, depending on the treatment given to the sub-grid gas dynamics in their simulations. Regarding the normalization of the observed relation, our data show a remarkably good agreement with the Narayanan et al. (2012) model for an assumed molecular surface density of Σ_H2 = 100 M_☉ pc⁻². On the other hand, the Feldmann et al. (2012a) model, in which X_CO is independent of the molecular gas surface density, predicts values that are a factor of ∼4 higher than our observations. In the following section we discuss the characteristics of these models in more detail, paying particular attention to the discrepancies which arise from the different treatments given to the role of the local gas surface density at setting the X_CO conversion factor.

In the theoretical model of Feldmann et al. (2012a), X_CO is a strong function of the molecular gas surface density of a given GMC, with the conversion factor growing away from its minimum toward both low and high densities. This behavior is a consequence of two different effects taking place in the low- and high-density regimes. At high column densities the numerator in Equation (1) increases linearly, while the denominator (i.e., the CO intensity) is either constant or increases sub-linearly (as in the virialized cloud approximation for which $I({\rm CO})\propto \Delta v \propto \sqrt{N({\rm H_2})}$). This effect is supported by spatially resolved observations of GMCs in the MW (Heiderman et al. 2010). Therefore X_CO increases toward high column densities. On the other hand, toward lower N(H2), the CO abundance drops due to CO dissociation as the dust extinction (which is proportional to the column density) decreases, giving rise to large "CO dark" envelopes and also implying an increase of X_CO. Both the location of the minimum, and the value of X_CO at the minimum are strong functions of the metallicity in the model.

When attempting to study these trends on kiloparsec scales within or across galaxies, we must take into account the fact that we measure averages for ensembles of molecular clouds with a distribution in their physical properties. In the Feldmann et al. (2012a) model, this average washes away the dependence on column density seen on small scales, and the authors predict almost no variation in X_CO as a function of H₂ surface density. It is important to point out that in this model the molecular gas temperature is fixed to 10 K in the calculation of I(CO), and the authors either impose a constant CO line width, or assume virialized clouds (i.e., $\Delta v \propto \sqrt{N({\rm H_2})}$).

As discussed in the previous section, the models of Narayanan et al. (2012) and Feldmann et al. (2012a) predict similar behaviors for X_CO as a function of metallicity. On the other hand, the simulations conducted by Narayanan et al. (2012) do not assume a fixed temperature for the molecular component of the ISM, and instead include relevant heating and cooling calculations, which allow the authors to trace the molecular gas temperature. Furthermore, for molecular clouds which are larger than the spatial resolution of the simulation, they use the simulations themselves to follow the dynamical state of the molecular gas (i.e., the gas velocity dispersion sets the CO line width). For sub-resolution clouds they assume virialization.

The ability to follow the temperature and dynamical state of the gas allows Narayanan et al. (2012) to study the impact of turbulence and heating from both ongoing star formation, and large-scale gas dynamics, on the X_CO factor under high-density conditions like the ones present in starbursts, mergers, and the central regions of galaxies. While gas density itself should have little direct impact on X_CO when averaging over kiloparsec scales, the existence of a correlation like the SFL implies that higher molecular gas density is accompanied by stronger star formation activity. Star formation is associated with feedback from protostellar jets, stellar winds, and supernova (SN) explosions, as well as heating (both photoelectric and by cosmic rays produced in SNe), which rises both the temperature and turbulence of the gas. Furthermore, high-density regions are also commonly associated with dynamically violent environments, in which for example, merger-induced shocks and cloud–cloud collisions can induce higher levels of turbulence in the gas. Given the optically thick nature of the CO(1–0) transition, any turbulence-induced line broadening translates into a brighter I(CO) and a lower X_CO. Similarly, a higher gas temperature translates into a higher brightness temperature for the CO line. By construction, these effects are not recovered in the Feldmann et al. (2012a) simulations.

Comparing the observed dependence of X_CO with gas surface density, or with CO surface brightness, to theoretical numerical model predictions is challenging, and a series of caveats arise which must be kept in mind when interpreting this type of comparison. These caveats arise on one side from the limitations of the observations, where the large physical spatial resolution typically achieved in extragalactic studies (∼600 pc in NGC 628 for the HERACLES beam size) dilutes the actual surface brightness (or mass surface density) of GMCs in the beam, by some factor which depends on both the intrinsic surface brightness distribution of the clouds, and the covering fraction of clouds across the area over which the emission is being integrated. In denser regions, overlap of optically thick clouds can also affect the measured surface brightness. On the other hand, in numerical simulations, while physical numerical quantities are readily accessible, the definition of a "cloud," and the exact method used to measure its surface brightness or density can be non-trivial. This can be especially problematic when the resolution of the simulation is larger, or of the same order of magnitude as the typical sizes of clouds.

Narayanan et al. (2012) define X_CO as the ratio between the average column density, and the average CO intensity (luminosity divided by area) for each of their simulated galaxies, and provide a fitted formula to calculate X_CO from measurements of the metallicity (Z' = Z/Z_☉) and the luminosity-weighted CO intensity of GMCs (〈I(CO)_GMC〉) in a galaxy:

$$\begin{equation} X_{\rm CO}=\frac{{\rm min}[4,\;6.75\times \langle I({\rm CO})_{{\rm GMC}} \rangle ^{-0.32}]\times 10^{20}}{Z^{\prime 0.65}}. \end{equation} \tag{ 4 }$$

Unfortunately, 〈I(CO)_GMC〉 is not recoverable by CO observations at kiloparsec scales, and also note that it is not the same quantity used in their definition of X_CO. Instead, of the luminosity-weighted average of the CO intensity of individual clouds, we measure the average CO intensity over large areas of the galactic disk of NGC 628, which is equivalent to the average intensity going in the definition of X_CO in the model.

Because of averaging over large areas, with a lower than unity cloud filling factor, our measured I(CO) values should be significantly lower than the CO intensities one would measure for individual GMCs. In order to use the Narayanan et al. (2012) model to make meaningful predictions from our data, we adopt a "clumping factor" to scale our observed average CO intensities, in an attempt to recover the luminosity-weighted average intensity of individual GMC in the observed regions. A similar approach is adopted by Krumholz et al. (2009) to transform observed average H₂ surface densities into intrinsic GMC surface densities for input to their SFL model.

The choice of the clumping factor is somewhat arbitrary, but it only affects the absolute value of the predicted X_CO factors, and not any relative trends observed across the galaxy. We use a factor of 30 to scale up the observed I(CO) values for input into Equation (10) in Narayanan et al. (2012). The justification for adopting this value is as follows. For an MW X_CO factor (which is consistent with the average factor we measure in NGC 628), the average CO intensity of NGC 628 inside a 7 kpc radius implies a mean value of Σ_H2 which is nine times lower than the typical surface density of Σ_H2 ≃ 100 M_☉ pc⁻² measured for GMCs in the MW and the Local Group (Roman-Duval et al. 2010; Rosolowsky et al. 2003; Bolatto et al. 2008). On the other hand, comparison of Equations (6) and (10) in Narayanan et al. (2012) implies that in their model Σ_H2∝〈I(CO)_GMC〉^0.64. Therefore, a clumping factor of 30 ≃ 9^1/0.64 would imply that the luminosity-weighted average GMC surface density in the inner 7 kpc of NGC 628 is similar to the typical MW and Local Group value.

The predicted X_CO values for different radial bins in NGC 628, computed using the measured metallicity and observed CO intensity (average between the HERACLES, BIMA, and CARMA measurements) corrected by this clumping factor, are shown as blue stars in Figure 5. Predictions for the Feldmann et al. (2012a) model would fall right on top of one of the two cyan solid lines in Figure 5, as in their parameterization X_CO is independent of Σ_H2 and I(CO). As stated in Section 5, the Narayanan et al. (2012) model is in better agreement with our data than the model of Feldmann et al. (2012a).

Figure 6 presents a comparison between the average X_CO values for each radial bin in NGC 628 (red stars) and the values predicted by the Narayanan et al. (2012) model (blue stars). Both are presented as a function of I(CO) scaled by a clumpiness factor of 30, which, as stated above, should be roughly equivalent to the input 〈I(CO)_GMC〉 for the model. By inspecting Figures 5 and 6 we can see that the model reproduces very well the observed trends in X_CO with both metallicity and CO surface brightness for the inner three radial bins in NGC 628. The slopes of both trends are consistent, and the absolute values agree to 1σ.

Zoom In Zoom Out Reset image size

A caveat to the above approach is that adopting a single clumpiness parameter is equivalent to assuming a constant cloud filling factor over the observed area. A decrease in the cloud filling factor toward the outer disk is in fact likely, and could indeed explain both the flattening seen in X_CO for the N = 1 case when using the stacked Schruba et al. (2011) data in the outer disk of NGC 628, and the discrepancy observed between the data and the Narayanan et al. (2012) model for the outer radian bin. This latter discrepancy is reduced if one assumes N = 1.5 instead of N = 1 for the SFL.

At this point, it is important to remind the reader about the many scaling factors which affect the normalization of the X_CO values. These include the zero points in the flux calibration of the HERACLES CO and VENGA Hα maps, and the adopted Hα SFR calibration (we have scaled the other data sets to this reference), the adopted offsets used to match the metallicity scale of the different data sets and the models, and the adopted clumpiness parameter which is constrained only by demanding a typical GMC surface density of Σ_H2 ≃ 100 M_☉ pc⁻². Given the inherent systematic uncertainties associated with all these scaling factors, it is difficult to draw strong conclusions from the relative offset seen between the observed and predicted X_CO values. None of these scaling factors affect the relative trends seen across different radial bins in both the data and the model.

The main conclusion we can draw from the above comparison is the following. Our observations are in good agreement with the Narayanan et al. (2012) model. If X_CO follows the metallicity and CO brightness dependence proposed in this model, that is X_CO∝〈I(CO)_GMC〉^−0.32 × Z'^−0.65, then, the observed metallicity and CO surface brightness gradients imply that both quantities contribute to the formation of an X_CO gradient across the disk of NGC 628. As stated in Section 5, the measured metallicity gradient implies Z'∝10^−0.04R, and a linear fit to the CO surface brightness of the inner three radial bins (R < 5 kpc) implies 〈I(CO)_GMC〉∝10^−0.07R. The model implies that X_CO∝10^{(0.32 × 0.07 + 0.65 × 0.04)R}, or a gradient of Δlog(X_CO) = 0.05 dex kpc⁻¹, in excellent agreement with the measured gradient of Δlog(X_CO) = 0.06  ±  0.02 dex kpc⁻¹. From the above calculation it can be seen that both the metallicity and the CO surface brightness (or equivalently the molecular gas surface density) contribute to the observed X_CO gradient. The Narayanan et al. (2012) prescription, in which X_CO is inversely proportional to the metallicity, and inversely proportional to the square root of the molecular gas surface density, is consistent with our data.

Our results indicate that across the disk of NGC 628, both the metallicity and the molecular gas surface density (through its incidence in the temperature and turbulence of the gas) are important factors setting the CO to H₂ conversion factor. Future analysis of a larger subset of the VENGA sample will allows us to know if this result is general enough to be applicable to most massive spiral galaxies in the local universe.

In this section we use the ionization parameter, as measured from H ii region emission line ratios in the VENGA data, as a proxy for the intensity of the local interstellar UV radiation field S_UV across the disk of NGC 628. The goal is to study the potential effect that the UV radiation field has at setting the value of X_CO.

The strength of the local UV radiation field has been proposed as an important factor setting the CO to H₂ conversion factor. Analyzing a sample of individual molecular complexes in the Large Magellanic Cloud and the Small Magellanic Cloud, Israel (1997) finds a linear correlation between the X_CO and the ratio of the far-IR surface brightness to H i column density ($\sigma _{{\rm FIR}}/N({\rm H\,\scriptsize{I}})$) which the author considers a proxy for the local UV radiation field per hydrogen atom (reprocessed after absorption and reemission by dust).

The effect of the UV radiation field on X_CO is also discussed by Feldmann et al. (2012a) in the context of their numerical simulations and their model of CO emission. The authors propose that for a single molecular cloud, the CO abundance decreases in the presence of a stronger UV radiation field U_UV (which they parameterize in units of the local interstellar UV radiation field, so U_UV = S_1000 Å/S_MW, with S_MW = 10⁶ photons cm⁻² s⁻¹ eV⁻¹). This dependence is stronger at low A_V (i.e., low N(H₂), or low Z, or both), but an order of magnitude change in U_UV can still change the CO abundance by factors of a few at A_V = 1–3 mag. These extinction levels correspond to the outer envelopes CO emitting regions in GMCs, where CO(1–0) is starting to become optically thick. As the UV radiation field goes up, the CO abundance goes down, and X_CO increases for a single cloud (see their Figure 2).

When investigating this effect on galactic scales, Feldmann et al. (2012a) find that the strong dependence of X_CO with the UV radiation field present in their model for single GMCs at fixed N(H₂) is completely suppressed in their simulated galaxies. This is due to the fact that in the simulations, physical parameters are measured over ensembles of molecular clouds with a distribution in their properties, and the abundance of clouds with low A_V, and therefore high X_CO, decreases strongly for higher values of the UV radiation field. In the Feldmann et al. (2012a) simulations, increasing U_UV suppresses the low-density tail of the N(H2) distribution function, dissociating diffuse molecular clouds with low densities (≲ 10²¹ cm⁻²) which have large X_CO values. Therefore, when averaging over kiloparsec scales, the radiation field has little impact in the average X_CO value, which is dominated by the surviving higher density regions with N(H₂) ≃ 10²² cm⁻².

A competing effect, which might be important at setting the value of X_CO, is photoelectric heating working on large polycyclic aromatic hydrocarbon molecules and small dust grains in the outer envelopes of molecular clouds (Tielens 2005; Glover et al. 2010). The photoelectric heating rate is directly proportional to the intensity of the interstellar UV radiation field (Tielens 2005), therefore, the gas temperature and the CO transition brightness temperature can be increased in the presence of a stronger UV radiation field. By construction, this effect cannot be recovered by Feldmann et al. (2012a) as the authors impose a constant molecular gas temperature of 10 K in their simulations. This effect might be present in the Narayanan et al. (2012) simulations, which follow the main heating and cooling processes in the gas, and are able to trace changes in the gas temperature, though the authors do not attempt to separate it from dynamical effects on the brightness of the CO line. As mentioned above, part of the dependence of X_CO with molecular gas surface density might indeed be associated with this mechanism, as higher gas density translates into higher star formation activity and hence, a stronger interstellar UV radiation field.

In Figure 7 we present the radial distribution of the ionization parameter

$$\begin{equation} q=\frac{S_{H^0}}{n}, \end{equation} \tag{ 5 }$$

where $S_{H^0}$ is the ionizing photon flux per unit area and n is the number density of hydrogen atoms. We have estimated q from the [O iii] λ5007/[O ii] λ3727 line ratio, as measured from the VENGA data cube of NGC 628, following the iterative procedures described in Kobulnicky & Kewley (2004). The measurements are limited to spaxels dominated by H ii region emission (black dots) as in the case of the metallicity measurements described in Section 5.1. Red circles in Figure 7 show median values in at low densities (≲ 10²¹ cm⁻²) 0.5 kpc wide radial bins, and red and green error bars show the standard deviation in each bin and the systematic uncertainty in the calibration, respectively.

Zoom In Zoom Out Reset image size

In NGC 628, the ionization parameter falls toward larger radii. A linear fit implies a gradient of Δlog(q) = 0.061 ± 0.003 cm s⁻¹ kpc⁻¹. A drop in the ionizing photon flux per hydrogen atom is consistent with the results of Aniano et al. (2012) who find a decreasing level of dust heating with radius due to a decrease in starlight intensity by modeling the Spitzer+Herschel spatially resolved dust SED of NGC 628. Considering the X_CO radial profile measured in Section 4, we find no evidence of an increase in X_CO in regions where the UV radiation field is enhanced, but we rather observe the opposite trend. The increase in X_CO for higher U_UV expected on clouds scales in the model of Feldmann et al. (2012a), and observed on molecular complexes scales by Israel (1997), is not seen in kiloparsec scales across the disk of NGC 628. This could be cause in part by the effects of averaging over large ensembles of clouds, discussed in Feldmann et al. (2012a) and described above, but this could only suppress this dependence.

The fact that we observe the opposite behavior is consistent with the photoelectric heating effect described above. It is difficult to decouple the impact that the UV radiation field might have on the gas temperature and hence, on the value of X_CO, from surface density dependence discussed in Section 6.1. All we can conclude, is that the observed radial distribution in the ionization parameter, implies that thermal effects might be important at setting X_CO. Further study will be necessary to properly model these thermal effects, and decouple them from dynamical effects which also impact the radiative transfer of CO photons.