Since stars form out of the cold, dense regions of the interstellar medium (ISM) where conditions favor the presence of molecules, the low-level rotational transitions of 12 CO are commonly used to study star-formation properties in galaxies (e.g., Solomon & Sage 1988;Fukui & Kawamura 2010;Kennicutt & Evans 2012). CO has been extensively targeted because it is the most abundant molecule after molecular hydrogen (H 2 ) and cold H 2 cannot be observed directly in emission. CO starts to form at visual extinctions of 1-3 mag corresponding to column densities of 1 -3 × 10 21 cm -2 at solar metallicity, while hydrogen becomes mostly molecular for column densities greater than a few 10 20 cm -2 . Under normal metallicity and moderate radiation field conditions, most of the cold, dense ISM is not expected to be dark in CO, and CO and its rarer isotopologues are expected to trace well the H 2 column density of clouds (Tielens & Hollenbach 1985;van Dishoeck & Black 1988;Sternberg et al. 2014).
Over the past decades, the 12 CO(J=1-0) emission has been calibrated to provide a measure of the total mass of molecular hydrogen via the CO-to-H 2 conversion factor X CO or α CO (e.g., Bolatto et al. 2013). In external galaxies, the main calibration techniques employed are based on the virial method, dust emission, optically thin molecular tracers, or radiative transfer of multiple molecules/transitions. These techniques often rely on making strong assumptions regarding, e.g., abundances, grain properties, filling factors, or the virialisation of molecular clouds. In addition, dependencies of the X CO factor on physical conditions within clouds, such as density, temperature, turbulence or metallicity are expected (e.g., Wolfire et al. 2010;Shetty et al. 2011;Narayanan et al. 2011). As a consequence, the X CO factor varies, as seen from galaxy to galaxy and across individual galaxies in observations (e.g., Arimoto et al. 1996;Downes & Solomon 1998;Leroy et al. 2011;Papadopoulos et al. 2012;Sandstrom et al. 2013;Cormier et al. 2014;Kamenetzky et al. 2014). In normal, star-forming disc galaxies, the amplitude of those variations is up to an order of magnitude (Bolatto et al. 2013;Sandstrom et al. 2013).
In addition to environmental dependencies of a luminosity-mass conversion factor, the high abundance and densities of 12 CO makes the J = (1 -0) transition optically thick in most molecular clouds, which complicates interpretation and can hamper accurate determination of cloud properties. Rarer isotopologues of the most abundant molecules containing carbon, such as 13 CO, are on the one hand less abundant and hence more difficult to observe in galaxies. On the other hand, they have the advantage of being more optically thin, allowing us to access the full column density of the material they arise from. Beyond our Galaxy (see Heyer & Dame 2015, for a review), 13 CO has mainly been observed in the centres or in small, targeted regions of nearby galaxies or integrated over entire, bright galaxies (e.g., Encrenaz et al. 1979;Young & Sanders 1986;Casoli et al. 1992;Aalto et al. 1995;Wilson et al. 1997;Paglione et al. 2001;Krips et al. 2010;Tan et al. 2011;Danielson et al. 2013;Alatalo et al. 2015;Vila-Vilaro et al. 2015;Sliwa et al. 2017a,b). With the "EMIR Multiline Probe of the ISM Regulating Galaxy Evolution" survey (EM-PIRE; Bigiel et al. 2016), we have obtained complete and high signal-to-noise maps of 13 CO as well as the main dense molecular gas tracers (HCN, HCO + , HNC) in the J = 1 → 0 transition across the discs of nine nearby spiral galaxies with the IRAM 30-m telescope (Bigiel et al. 2016;Jiménez-Donaire et al. 2017b,c;Gallagher et al. 2017). We also obtained full maps of 12 CO(1-0) emission for those nine galaxies in follow-up programs. These galaxies are drawn from the HERACLES 12 CO(2-1) survey (Leroy et al. 2009). They are selected to have diverse structural properties (barred/unbarred, flocculent/grand-design spiral arms) and to reside in different environments (field/Virgo Cluster galaxy) in order to test whether these parameters influence their observable ISM properties.
In this paper, we present an analysis of the 13 CO(1-0) and 12 CO(1-0) observations from EMPIRE. Since 13 CO has not been observed or mapped as extensively as 12 CO in galaxies, our goal is to investigate variations in the 12 CO(1-0)-to-13 CO(1-0) ratio among and within those galaxies, and to understand if/under which conditions 13 CO may be a better tracer of the molecular gas mass than 12 CO. Throughout this paper, the 12 CO(1-0)-to-13 CO(1-0) integrated intensity ratio is denoted ℜ. Section 2 describes the observations. Section 3 presents an analysis of ℜ with radial profiles and correlation diagrams, as well as a qualitative comparison to models and a derivation of column densities. The physical origin of ℜ variations and the ability of the CO lines to trace the molecular gas mass are discussed in section 4. Finally, we summarize our conclusions in section 5.
2.1 EMPIRE observations of 13 CO(1-0)
Observations of the IRAM 30-m large program EMPIRE (PI Bigiel) were carried out in 2012 (pilot program for NGC 5194) and between December 2014 and December 2016 (program 206-14, for the other galaxies). We mapped the full discs of 9 nearby spiral galaxies with the EMIR E0 receiver in the on-the-fly mapping mode. The half-power beam width (HPBW) at 110 GHz is 22 arcsec and the adopted spectral resolution is 4 km s -1 .
The data were reduced with our in-house pipeline. The main steps include: baseline subtraction with a polynomial function of order 2, rejection of spectra above 3 times the theoretical noise, conversion to main beam temperature assuming main beam and forward efficiencies of 0.78 and 0.94, projection of the spectra onto grids of pixel size 4 arcsec. After gridding, the full width half-maximum (FWHM) of the 13 CO(1-0) data is 27 arcsec. This corresponds to a linear resolution of ∼1.5 kpc for our sample of galaxies. We refer to Jiménez-Donaire et al. (2017c,a) for a detailed description of the data acquisition and reduction. Line calibrators were observed during each run of the campaign and their intensities vary by about 5 per cent only. Table 1 provides the list of our targets and noise levels achieved at the frequency of the 13 CO(1-0) line (110.20 GHz).
We used the 12 CO(2-1) data from HERACLES (Leroy et al. 2009) as a guide to create integrated intensity maps for the 13 CO(1-0) line. The data were retrieved from the HER-ACLES repository 1 , convolved to a common resolution of 27 arcsec using Gaussian kernels, and put on the same spatial grid as the EMPIRE data using the IDL procedure hastrom. At each position in the map, we fitted the 12 CO(2-1) line with a single Gaussian. For pixels below a signal-to-noise ratio of 5 for the velocity-integrated intensity, we interpolated central velocities and line widths from well-detected neighboring pixels by fitting a plane to the maps. In regions where most pixels are not detected in 12 CO(2-1), the interpolated central velocities are not allowed to take values lower (higher) than the minimum (maximum) velocity measured in the well-detected pixels, and the interpolated line widths are set to the average line width measured in the well-detected pixels.
The central velocities and line widths of the 12 CO(2-1) line were used as initial guesses for those of the 13 CO(1-0) line. We created intensity maps for 13 CO(1-0) both by fitting and integrating a single Gaussian and by integrating directly the signal in specific velocity windows. We defined the windows as ∼ 3× FWHM of the 12 CO(2-1) line, thus the windows vary for each line of sight. The two methods yield differences in integrated intensities that are typically less than 7 per cent. Since the CO line profiles are not always Gaussian, especially in the galaxy centres, we prefer to use the direct integration maps and not the line-fitted maps. We also produced error maps. For each pixel, the error on the integrated intensity is calculated as the standard deviation in the line-free parts of each spectrum, multiplied by the square root of the number of (4 km s -1 wide) channels inside the FWHM given by the Gaussian fit. Figure A1 in the Appendix shows final integrated intensity maps for each galaxy.
In order to make homogeneous and matched-quality measurements of ℜ, we performed new observations of 8 galaxies of the EMPIRE survey (all except NGC 5194) in the 12 CO(1-0) line, with the IRAM 30-m telescope. Maps of the entire discs were obtained as part of the programs 061-15, 059-16 (PI Jiménez-Donaire) and . The data were reduced with the same pipeline as for EM-PIRE and we produced final cubes at a spatial resolution of ≃25 arcsec and a spectral resolution of 4 km s -1 . We reached sensitivities of 17-30 mK (T mb ) per 4 km s -1 channel. For NGC 5194, the 12 CO(1-0) and 13 CO(1-0) data are taken from the PAWS survey (30-m cubes 2 ; Schinnerer et al. 2013;Pety et al. 2013).
All datasets are convolved to a common resolution of 27 arcsec using Gaussian kernels and put on the same spatial grid as the EMPIRE data. Integrated intensity maps are 1 http://www.iram-institute.org/EN/content-page-242-7-158-240-242-0.html 2 http://www.mpia.de/PAWS/PAWS/Data.html created following the same steps as for the 13 CO(1-0) data, described in section 2.1.2.
All of the EMPIRE targets were observed with GALEX as well as with the MIPS instrument onboard Spitzer as part of the programs LVL (Dale et al. 2009) and SINGS (Kennicutt et al. 2003). Herschel photometry exists for all of our targets except NGC 2903. We use user-provided products from the key programs VNGS and KINGFISH (data release 3) (Bendo et al. 2012;Kennicutt et al. 2011). No PACS 100 µm observations are available for NGC 5194. The reduction and map-making were done in HIPE versions 9 and 8 and Scanamorphos versions 21 and 16.9, respectively. Convolutions are done using the kernels from Aniano et al. (2011).
From the photometry and spectroscopy, we derive physical quantities such as star-formation rate surface densities (Σ SFR ), molecular gas surface densities (Σ mol ), and H 2 column densities (N(H 2 )). All surface densities are corrected for inclination.
Our SFR estimates are based on IR data. For all galaxies overlapping with KINGFISH, we use the TIR maps calculated with dust models from Galametz et al. (2013). For NGC 5194 and NGC 2903, we compute TIR surface brightness using the generic calibration from Galametz et al. (2013). We combine the bands MIPS 24 µm, PACS 70, 160 µm and SPIRE 250 µm for NGC 5194 and the bands MIPS 24 and 70 µm for NGC 2903. For the galaxies in KING-FISH, the comparison of TIR maps calculated from SED models and from the generic calibration yield differences of about 10 per cent for the first combination (MIPS, PACS and SPIRE) and 20 per cent for the second combination (MIPS only). The generic calibration tends to systematically overpredict slightly the TIR flux in the brightest regions. Uncertainty maps are generated by adding in quadrature errors on the fluxes and errors on the calibration coefficients and from the choice of method as quoted above in the case of the TIR calibration. The TIR maps are then converted to starformation rate surface density maps using the calibration from Murphy et al. (2011).
Molecular gas mass surface densities and H 2 column densities are commonly derived using the 12 CO(1-0) line as a proxy for H 2 and a CO-to-H 2 conversion factor. For Σ mol , the standard conversion factor is α CO of 4.4 M ⊙ pc -2 (K km s -1 ) -1 which includes helium (Bolatto et al. 2013). For N(H 2 ), the standard conversion factor is X CO of 2 × 10 20 cm -2 (K km s -1 ) -1 which does not account for helium. These are considered as reference, Milky-Way values. Variations of the conversion factors are discussed in section 4.2.
The 12 C/ 13 C isotope abundance ratio results from stellar processing and evolution (e.g., Wilson & Rood 1994). The 12 CO/ 13 CO abundance ratio further depends on chemical processes within molecular clouds. While we discuss possible sources of variations in those abundances in section 4, in our analysis we make no difference between the isotope ( 12 C/ 13 C) abundance ratio and the isotopologue ( 12 CO/ 13 CO) abundance ratio. The abundance ratio is largely unconstrained for external galaxies but has been observed to vary within galaxies and from galaxy to galaxy.
In the Milky Way, it increases from 25 in the centre to 100 in the outer disc, with a value of ∼70 in the solar neighborhood (Langer & Penzias 1990;Wilson & Rood 1994;Milam et al. 2005). In the Large Magellanic Cloud, it is ∼50 (Wang et al. 2009). In starburst galaxies, it is found to be >40, and in some ULIRGs it is even > 100 (e.g., Martín et al. 2010;Henkel et al. 2014;Sliwa et al. 2017b). Here, we adopt a fiducial value for the isotope abundance ratio of 60. In the Milky Way, this value corresponds to a distance of ≃ 6.7 kpc or 0.6 r 25 .
The 12 CO/H 2 abundance ratio is also sensitive to radiative processes and gas chemistry. It can vary by orders of magnitudes from translucent to dense lines-of-sight, with a typical scatter of 0.5 dex at any given H 2 column density (van Dishoeck et al. 1992;Sheffer et al. 2008). In our large telescope beam, such different lines of sights are mixed and beam-to-beam variations of the abundance are expected to be smaller. In the following, we take the canonical value of 10 -4 for the 12 CO/H 2 abundance ratio, which corresponds to a 13 CO/H 2 abundance ratio of 1.7 × 10 -6 (e.g., Dickman 1978) and could be uncertain by a factor of a few.
3.1 13 CO(1-0) intensities across galaxies Figure 1 shows the 13 CO(1-0) and scaled 12 CO(1-0) intensities as a function of the PACS 70 µm intensity, which can be seen as a proxy of the star-formation rate surface density, as well as histograms of the 13 CO(1-0) intensities for each galaxy. 13 CO(1-0) is detected throughout our maps with signal-to-noise ratio > 10 in the brightest regions, and the majority of the pixels in our maps are detected at a > 5σ level. 13 CO(1-0) peaks on the galaxy centres, except in NGC 3627 where it peaks on the star-forming knots at the end of the bar (though this may be somewhat an artifact of resolution). In the inter-arm/outer-disc regions, the 13 CO(1-0) intensity has values <1 K km s -1 (<0.3 K km s -1 for the two faintest galaxies NGC 0628 and NGC 3184).
We quantify correlations with the Spearman's rank correlation coefficient and its significance using the IDL procedure r_correlate.pro. The correlation coefficient is computed for data points with signal-to-noise ratio of the 13 CO(1-0) intensity greater than 5. With this condition, the signal-to-noise of the ancillary data sets described above is generally not a concern. The significance corresponds to the p-value or probability of null hypothesis. It is reported in parenthesis in the figures. The uncertainty (±) on the coefficient is estimated with a Monte-Carlo simulation. We find that 13 CO(1-0) and 12 CO(1-0) are generally well correlated with the dust emission from PACS 70 µm. Both CO lines show very similar distributions and scatter, with departure from each other in some galaxy centres (e.g., NGC 5055). At the resolution of our data, their line profiles and line widths are also similar (see Fig. A3).
Given the high quality of the data, we perform a lineof-sight based analysis. Only in the following section 3.2, we perform a stacking analysis to extract representative spectra in the different environments of our sample of galaxies.
ℜ is the ratio of the 12 CO(1-0) intensity and the 13 CO(1-0) intensity, which have units of K km s -1 . We measure ℜ in the galaxy centres, arm, inter-arm regions and total emission in our maps by means of stacking. In the stacking step, spectra corresponding to a given region are aligned in velocity using the 12 CO(2-1) data as reference and averaged. Jiménez-Donaire et al. (2017b) describes the method in detail and Figure A3 shows stacked spectra of 13 CO(1-0) and 12 CO(1-0) for the galaxy centres and entire galaxies. CO intensities are then measured by direct integration of the stacked spectra over a velocity window large enough to encompass all the signal. Table 2 reports the ℜ values ob-Figure 1. Top panels: EMPIRE observations of the CO line intensity (in K km s -1 ) on the y-axis following well the broadband 70 µm intensity (in MJy sr -1 ) which is a proxy for the star-formation rate surface density on the x-axis. 13 CO(1-0) data are in colour and 12 CO(1-0) data, scaled by a factor of 10, are in grey. For both CO lines, smaller circles correspond to pixels below a signal-to-noise ratio of 5 in the respective maps. The Spearman's rank correlation coefficients between each CO line intensity and the 70 µm intensity, their uncertainty, and their significance (in parenthesis), are indicated. Those are measured on data with signal-to-noise ratio above 5. Bottom panels: Histogram of the 13 CO(1-0) intensities. Only few pixels are below the 1σ (5σ) noise level as indicated by the hashed black (beige) parts of the histograms.
tained this way for each galaxy as well as averages and dispersions over the sample. Figure A1 in the Appendix shows maps of ℜ with contours delineating the centre and arm regions. The regions were defined using cuts in intensity of the 12 CO(1-0) line (at our working resolution of 27 arcsec or ≃1.5 kpc), and the cut levels were chosen by eye. 'Centre' refers to the galaxy centres (distance to the centre of the galaxy <16 arcsec or inner 0.8 kpc, typically) where the 12 CO(1-0) emission is brightest (first cut at CO intensity levels of: 5.5, 35, 7, 30, 32, 22, 30, 32, 40 K km s -1 , for NGC 0628, 2903, 3184, 3627, 4254, 4321, 5055, 5194, 6946, respectively). 'Arm' refers to the galaxy discs where the 12 CO(1-0) emission is bright (second cut at CO intensity levels of : 3, 18, 3.5, 11, 11, 8, 10, 10, 10 K km s -1 , for NGC 0628, 2903, 3184, 3627, 4254, 4321, 5055, 5194, 6946, respectively). 'Inter-arm' refers to the inter-arm and outer parts of the maps where the 12 CO(1-0) emission is fainter but detected at a signal-to-noise ratio > 3. Finally, 'total' refers to the entire map where 12 CO(1-0) is detected. We note that we also defined contours by hand, identifying centres and arm regions based on 24 µm and 70 µm continuum images. Values of ℜ obtained by stacking with these hand contours vary by at most 10 per cent compared to values reported in Table 2. This 10 per cent discrepancy can be seen as a methodology uncertainty.
Thanks to our sensitive, full maps, we can achieve much better galaxy-integrated and environment-specific measurements than previous investigations of ℜ in nearby galaxies. In all galaxies and all environments, the noise of the data indicates that we could have measured ratios up to 150 in the stacks (but the observed ratios are much lower). The mean global (i.e. full-galaxy) ℜ value that we measure over our sample is 11 with a standard deviation of 1.5. When the galaxies are divided into different environments (centre, arm, inter-arm regions), we obtain similar average ℜ values, with standard deviation ∼2. We notice that the dispersion in ℜ is highest in the galaxy centres, with values varying between 7 and 15 (see Table 2). The dispersion is also high in the inter-arm regions, but with larger error bars on individual measurements. Studying, also at a kpc-scale, centres of ∼10 nearby galaxies with AGN activity, signs of a recent merger, or an intense central starburst, Israel (2009a,b) found values of ℜ in the range 8 -16.
In the EMPIRE survey, we find that the galaxies with bright, starburst-dominated nuclei (NGC 2903, NGC 3627, NGC 4321, and NGC 6946) have the largest central ℜ values. Two galaxies in our sample have strong bars (NGC 2903 and NGC 3627) and three other galaxies have weaker bars (NGC 3184,NGC 4321,NGC 6946,and possibly NGC 5194).
Those barred galaxies also show generally higher central ℜ values than the non-barred galaxies. Both properties of having a bar and a starburst-dominated nucleus might be related as bars may help to funnel gas towards the galaxy center -leading to high dense gas fractions (Gallagher et al. 2017) -and to fuel star formation (e.g., Ho et al. 1997). Finally, two galaxies with AGN activity have low central ℜ values (NGC 5055, NGC 5194) but the trend is not systematic (e.g., NGC 3627 also has AGN activity).
We generate radial profiles for the main tracers discussed in this paper: 13 CO(1-0), 12 CO(1-0), and 12 CO(2-1). We choose a step size in radius of 13.5 arcsec, corresponding to half of the spatial resolution. At a given radius, we measure the average of the intensities within a tilted ring. The width of the rings is taken as the chosen step size (defined along the minor axis), so that the rings are not overlapping, but they are correlated because the beam is oversampled. Table 1 gives the assumed inclinations and position angles. Error bars on those measurements are calculated as the root-mean-squared error within each ring of our error maps and multiplied by the square root of the oversampling factor N s , defined as N s = 1.13×(map resolution/pixel size) 2 . All intensity profiles are multiplied by cos(i) to correct for inclination. Profiles of intensity ratios, such as ℜ, are built by dividing the radial profile of the quantity at the numerator with the radial profile of the quantity in the denominator. Figure A2 shows radial profiles of the individual CO line intensities for each galaxy and Figure 2 shows radial profiles of ℜ for all galaxies.
Radial profiles of 13 CO(1-0), 12 CO(1-0), and 12 CO(2-1) globally follow each other very well for all galaxies (Fig. A2).
Intensities peak in the centre and decrease by an order of magnitude at r ≃ 0.6 r 25 , and even more for NGC 6946 which has a very prominent centre. Differences in the radial behavior of 13 CO(1-0) and 12 CO(1-0) are made more obvious by inspecting ℜ (Fig. 2). We have reliable measurements of ℜ out to r ≃ 8 kpc or r ≃ 0.7 r 25 .
• The profiles of ℜ in NGC 2903, NGC 3627, and NGC6946 decrease by a factor 1.2 -1.4 from centre to disc (r = 3 kpc) and stay mostly flat at larger radii. For NGC 3627, we notice a clear suppression of the 13 CO(1-0) peak intensities in the galaxy centre (see Figures A1 andA2 in the Appendix).
• The profiles of NGC 3184 and NGC 4321 are flat at all radii.
• The profiles of NGC 4254, NGC 5194, and NGC 5055 increase steadily by a factor of 1.4 from centre to outer disc (r ≃ 7 kpc).
• The profile of NGC 0628 is peculiar. It increases by a factor of 1.8 from centre to r = 2 kpc and stays mostly flat at larger radii.
Centres aside, we find that the profiles increase slightly as a function of radius on average. The increase is mild compared to that observed in the Milky Way (Roman-Duval et al. 2016). At large radii (r = 7 -8 kpc), our sample of galaxies and the Milky Way have similar ℜ values, but at lower radii, the Milky Way shows systematically lower values, between 5 and 8. ℜ increases by a factor of ∼ 2 from the inner disc to the outer disc of the Milky Way. Such measurements, along with observations of rarer isotopologues of CO, have been used to infer abundance gradients in the Milky Way (Milam et al. 2005, shown as the black dotted line in Fig. 2). ℜ profiles are different in our sample of galaxies, implying that the Milky Way as a massive, rather quiescent galaxy (e.g., Chomiuk & Povich 2011) that may be compact for its stellar mass, might represent a different physical regime than the galaxies in our survey (e.g., different abundance patterns), or that there could be issues (e.g., geometric) in comparing Milky Way and extragalactic work. For example, the study of Roman-Duval et al. (2016) was limited to regions less than 50 pc away from the Galactic plane and it could be missing substantial 12 CO emission at high latitude (Dame et al. 2001) that we do capture in our sample.
Paglione et al. ( 2001) observed 13 CO(1-0) and 12 CO(1-0) along the major axes of NGC 3184, NGC 3627, NGC 5055, NGC 5194 and NGC 6946 with the FCRAO 14-m telescope (beam size ∼47 arcsec). Their central and outer values globally agree with our observations within errors, except for NGC 3184 for which Paglione et al. (2001) find ℜ values two times lower but with marginal detections. For NGC 3627, we find a lower average in the disc (12 instead of 17), but they have much lower statistics. For NGC 6946, we find similar central values but larger values in the disc (13-14 instead of 2016) and of the 12 C/ 13 C abundance gradient (black dotted line, divided by a factor of 7 for display) from Milam et al. (2005). Right panel: same as the left panel with radii normalized to r 25 (x-axis) and the profiles normalized to their value at a radius of r = 0.2 r 25 . We assume r 25 = 11.5 kpc for the Milky Way (de Vaucouleurs & Pence 1978). 10). At high resolution (5 arcsec or ∼150 pc), Meier & Turner (2004) find a range of values between 7-20 in the nucleus of NGC 6946, which is compatible with our central average of 15. Muraoka et al. (2016) mapped NGC 2903 with the NRO 45-m telescope (beam size ∼14 arcsec). They find values of ℜ that are around 10, which is globally consistent with our results. They divide the galaxy in ten distinct regions and find that ℜ varies by a factor of about two in the different environments. Although our resolution is coarser, we do not find as large spatial variations within NGC 2903 as they do. However, we achieve much better signal-to-noise ratios for individual 13 CO(1-0) measurements. Garcia-Burillo et al. (1993) 2015) also observed the centres of NGC 0628, NGC 2903, NGC 4254, and NGC 5055 with the ARO KP 12-m telescope (beam size ∼56 arcsec). We find similar values as theirs, except for the centre of NGC 4254 (8 instead of 11.5). The discrepancy could be attributed to beam size differences, as we find a global average for NGC 4254 closer to 11, or to calibration uncertainties. Our value in the centre of NGC 4254 is closer to that reported by Li et al. (2015).
In high-resolution mapping observations of 23 nearby disc galaxies selected to lie on the blue sequence and to be actively star-forming, IR-bright galaxies (CARMA STING surveyfoot_1 ), probing scales of 300-500 pc, Cao et al. (2017) find similar ℜ ratios and flat ℜ profiles. They report that ℜ varies mostly from galaxy to galaxy, with values between 5 and 15. This behavior resembles what we find in the centres of our sample of galaxies, but our profiles in the discs show less scatter. While they are limited by sensitivity and focus on the inner bright, molecular gas-rich regions, we detect fainter and more diffuse regions with EMPIRE, enabling us to probe a wider range of environments (centre, arm, interarm regions). gal Sakamoto et al. (1997) performed stripscan observations of the edge-on spiral NGC 891 at a resolution of 14 arcsec. They find a high value of ℜ ≃ 15 in the centre and a general increase of ℜ, from 4 at r = 3 kpc to 20 at r = 10 kpc. The increase of ℜ with radius is less pronounced in our sample of galaxies.
Overall, there is generally good agreement between our ℜ measurements and those reported in the literature for the EMPIRE galaxies or for similar types of galaxies and at similar scales. When there are discrepancies, these tend to be due to noise in the data. The uniqueness of the EMPIRE survey lies in the large, homogeneous spatial coverage (out to r ≃ 0.5 r 25 ) and high signal-to-noise ratios achieved.
In this section, we investigate how ℜ correlates with physical parameters such as the line ratio CO(2-1)/CO(1-0), the Figure 3. Correlation between the 12 CO(1-0)-to-13 CO(1-0) intensity ratio ℜ and several quantities: the 12 CO(2-1)/ 12 CO(1-0) ratio probing the gas physical conditions, the PACS70/PACS160 flux density ratio probing the dust temperature, the star-formation rate surface density Σ SFR , and the Σ SFR /CO(1-0) ratio which is a proxy for the star-formation efficiency. Quantities are averaged over entire galaxies and do not show significant correlations. Uncertainties on the galaxy averages are plotted but smaller than the symbol sizes. We also indicate the Spearman's rank correlation coefficients, their uncertainty, and their significance (in parenthesis). The grey data correspond to all kpc-size data points with smaller symbol size for 13 CO(1-0) data with signal-to-noise ratio below 5. The binned values, in colour, consider all data points in the bin. We also indicate Spearman's rank correlation coefficients, their uncertainty, and their significance (in parenthesis), which are measured on data with signal-to-noise ratio above 5.
IR colour PACS70/PACS160, the star-formation rate surface density (Σ SFR ), and the Σ SFR /CO(1-0) ratio (∝ starformation efficiency). We use the CO(2-1)/CO(1-0) ratio as a probe of gas conditions (temperature, density, opacity), and the IR colour PACS70/160 as an indicator of dust temperature that could be coupled to the gas temperature in the case where densities are high enough that collisions between dust grains and the gas particles lead to equal gas and dust temperatures.
Figure 3 shows values averaged over entire galaxies while Figure 4 shows the individual ∼kpc-sized regions as well as binned values for each galaxy. To produce these measurements, we first sample our maps with a hexagonal grid of spacing 13.5 arcsec (half of the spatial resolution). For total averages (Fig. 3), we measure the average of each observed quantity (CO intensity, PACS flux density, Σ SFR , etc.) by considering all sampling points, and for ratios, we divide these averages with each other. For the bins, we construct irregular bins of the physical quantities such that each bin contains 20 sampling points. Similarly (Fig. 4), we measure the intensity of each CO line and its error, either by considering the individual sampling points (for the kpc-size regions) or by averaging intensities of the sampling points falling in the bin (binned averages). We then calculate ℜ by taking the ratio of these intensities. We do not mask sampling points with low signal-to-noise when doing the binned averages. Error bars are calculated as the dispersion in ℜ measurements obtained with a Monte-Carlo simulation.
Inspecting Fig. 3, we find that the global averaged values of ℜ are not correlated with any quantity. Looking at integrated galaxies, Aalto et al. (1995) find that the warmest galaxies (F60/F100≥0.7), which are absent in our EMPIRE sample, display systematically high ℜ values (> 10), while ℜ and F60/F100 do not correlate for cold-intermediate IR colours (F60/F100≤0.6), which is in line with our results. Moreover, Davis (2014) report a correlation between ℜ and Σ SFR for integrated early-type galaxies. The Σ SFR range that they investigate is much larger than ours (5 orders of magnitude as opposed to a dynamic range in the average Σ SFR of less than 1 dex in our case). This may be why we do not find a significant trend with Σ SFR . We discuss Figure 4 in detail next.
CO(2-1)/CO(1-0) intensity ratio (Fig. 4a): We find that ℜ systematically decreases with increasing CO(2-1)/CO(1-0) ratio within galaxies. Moreover, for a given value of CO(2-1)/CO(1-0), there are significant offsets in ℜ from galaxy to galaxy. The 13 CO(1-0), 12 CO(1-0) and 12 CO(2-1) lines have upper energy levels E/k of 5.3 K, 5.5 K and 16.6 K, and critical densities of 2 × 10 2 cm -3 , 1 × 10 3 cm -3 and 6.7 × 10 3 cm -3 (for T kin = 20 K and τ line ≃ 1), respec-tively. We note that the critical density is a function of optical depth when line trapping effects are important (e.g., Scoville & Solomon 1974;Shirley 2015). In the case of 12 CO, the 12 CO(2-1)-to-12 CO(1-0) ratio would not only be sensitive to the temperature and density of the gas but also to the optical depth of 12 CO which may play a role in dictating where sub-thermal excitation happens (Peñaloza et al. 2017). On kpc-scales, our sample of galaxies span a range of CO(2-1)/CO(1-0) ratios between 0.3 and 2 and show a weak-to-moderate anti-correlation between ℜ and CO(2-1)/CO(1-0), indicating that ℜ decreases for increasing temperature/density/opacity (see Fig. 7 for model predictions). Both ℜ and CO(2-1)/CO(1-0), i.e. the y-axis and x-axis in Fig. 4a, are correlated via 12 CO(1-0). We investigate if the correlations are real or not by performing two Monte-Carlo tests that are detailed in Appendix B. Those tests indicate that the observed scatter in the line ratios is physical, it cannot be explained purely by noise. They also indicate that the correlation coefficients that we measure are robust and not driven by the correlated axes. Hence both the variations in the line ratios and the observed correlations are real. PACS70/160 IR colour (Fig. 4b): All galaxies in our sample span a similar range of PACS70/PACS160 values (0.2 -0.8). This IR colour can be used as a proxy for the dust temperature. We observe a weak-to-moderate anti-correlation of ℜ with PACS70/PACS160 for 5 galaxies: NGC 0628, NGC 4254, NGC 5055, NGC 5194, and NGC 6946. The high PACS70/PACS160 and low ℜ values are found in the centre of those galaxies, except for NGC 6946, where they are found in more diffuse regions of the disc. Comparing ratio behaviors between Fig. 4a and Fig. 4b, temperature effects may be at play in those five galaxies, but it is probably not the dominant/unique condition affecting ℜ.
SFR surface density (Fig. 4c): All galaxies in our sample span a similar range of Σ SFR values, with NGC 0628 and NGC 3184 being a bit less active. We observe a weak correlation between ℜ and Σ SFR for NGC 3184 and a moderate anti-correlation for NGC 4254 and NGC 5055. Those two galaxies have ℜ profiles steadily increasing with radius and also show an anti-correlation with the dust temperature. For NGC 3627, the behavior of ℜ with Σ SFR is not immediately apparent from the radial profiles (Fig. 2), probably because the star-forming knots at the end of the bar and the centre have high Σ SFR but different ℜ values. For the CARMA STING survey, Cao et al. ( 2017) also find no trend on global scales and moderately decreasing ℜ with increasing Σ SFR for some galaxies, their galaxies probing the range
Σ SFR /CO(1-0) (Fig. 4d): All galaxies in our sample span a similar range of Σ SFR /CO(1-0) values (0.3 -2 Gyr -1 ), which we have normalized such that they have unit of starformation efficiency (SFE mol ) under a constant, Galactic α CO value. We find that ℜ is moderately anti-correlated with Σ SFR /CO(1-0) within all galaxies, i.e. that ℜ is lower at high efficiencies, though the two quantities are correlated by construction. As for the CO(2-1)/CO(1-0) intensity ratio, we test those correlations with Monte-Carlo simulations that are detailed in Appendix B. The test indicates that the scatter in the Σ SFR /CO(1-0) values is physical and cannot be purely explained by noise. The correlation coefficients are also robust and highest for NGC 5055 and NGC 6946. Previous works have suggested that variations in the starformation efficiencies of massive galaxies could be linked to, e.g., a change in the relative fractions of diffuse and dense molecular gas (Saintonge et al. 2012;Shetty et al. 2014), or to enhanced CO excitation in galaxy centres (Leroy et al. 2013). Both effects have implications on the gas opacity and, though moderate, our trends indicate that the optical depth of CO, through ℜ, could indeed account for some changes in star-formation efficiencies.
We also explored how ℜ correlates with UV/TIR (not shown). This ratio can be viewed as an indicator of the visible/obscured star-formation activity and, for a given geometry, of the ISM photodissociation/shielding ability. The EM-PIRE galaxies span a wide range of values for the UV/TIR luminosity ratio. However, we find no clear correlation between ℜ and UV/TIR within galaxies.
Overall, most of the trends or lack of trends of ℜ with the physical parameters (Fig. 3 and Fig. 4) resemble the trends with galactocentric radius (Fig. 2). None of the parameters is revealing a strong positive or negative correlation, indicating that the kpc-scale resolution may be too coarse to isolate and identify the local effects affecting ℜ.
3.5 Correlation between the X CO factor and ℜ Since the physical quantities described above may also influence the conversion of 12 CO(1-0) intensity to molecular gas mass, we investigate empirically and with models how the variations in ℜ and in the X CO factor, calibrated on the dust reference, are linked. The X CO factor can be determined indirectly from dust emission. Modeling of the dust emission provides a dust mass that is converted to a total gas mass with a dust-togas ratio and to a molecular mass by subtracting the mass of atomic gas. This molecular gas mass or surface density (denoted below Σ(mol, dust)) is then divided by the CO intensity to determine X CO or α CO . Leroy et al. (2011) first applied this by solving for the dust-to-gas ratio and α CO simultaneously on spatially-resolved scales of a few kpc for the Local Group. That method was further employed in Sandstrom et al. (2013) for the HERACLES galaxies (including 8 out of the 9 galaxies of our sample) and is similarly applied to NGC 5194 in Leroy et al. (2017b) and Groves et al. (in prep.). In Sandstrom et al. (2013), α CO is measured for hexagonal pixels of size 37.5 arcsec, using molecular surface densities from dust emission and HERACLES 12 CO(2-1) observations. They formally measure α CO(2-1) that they express as α CO(1-0) because 12 CO(1-0) is more commonly used. However, they adopt a fixed CO line ratio of 0.7. Here, with measurements of the 12 CO(2-1)/ 12 CO(1-0) ratio for each location, we update those α CO values. As we Figure 6. Correlation within galaxies of ℜ and the conversion factor α CO , which is calibrated on dust emission. The trends are quite different from galaxy to galaxy. The colour-code corresponds to distance to the galaxy centre and it is the same for all panels. The dashed lines are shown to guide the eye. We also indicate the Spearman's rank correlation coefficients and their significance in parenthesis.
find 12 CO(2-1)/ 12 CO(1-0) ratios spanning a range of values between 0.3 and 2 on a kpc-scale in our sample of galaxies, we expect the α CO values to change a bit. We calculate the average CO intensities and ℜ value in a circular aperture (roughly matched radius of 34.1 arcsec) centered on each of the pixels. The new conversion factor α CO is equal to:
(1)
We show a pixel-to-pixel comparison of the old and new α CO values in Figure 5, and we report averages of the new values, of ratios with the old values and dispersions in Table 3. The new α CO values are about 15 per cent lower than the old values. There is a significant offset for NGC 3627 because its 12 CO(2-1)/ 12 CO(1-0) ratio is the lowest, with an average value of 0.5 as suggested by our new IRAM 30-m observations as opposed to the canonical value of 0.7 assumed by Sandstrom et al. (2013).
Figure 6 shows ℜ as a function of α CO , normalized to the Milky Way value. ℜ and α CO are positively correlated in NGC 5055 and NGC 3184, and more weakly in NGC 0628. ℜ and α CO are anti-correlated in NGC 3627, NGC4321, and rather weakly in NGC 6946, and the lowest α CO values (≤ 0.1 α CO,MW ) are reached in those galaxies. In galaxy centres, one possible explanation for the low α CO values observed is a change in the 12 CO optical depth (e.g., Sandstrom et al. 2013). Those optical depth effects seem to be at work in NGC 3627, NGC 6946, and to a lesser extent, NGC 4321.
To visualize how changes in physical conditions, specifically optical depth, density, and temperature, can affect ℜ quantitatively, we compare our observations to models from the non-LTE code RADEX (van der Tak et al. 2007). We note that those physical conditions, especially density and optical depth, could be correlated in real molecular clouds.
We use the grid of models built by Leroy et al. (2017a) to predict ℜ, the 12 CO(2-1)/ 12 CO(1-0) intensity ratio, and X CO . These models consider a distribution of volume densities (as an approximation for sub-beam density variations) and take into account sub-thermal excitation, but they assume fixed abundances and the spatial coexistence of 12 CO and 13 CO that fill the beam and are at a common kinetic temperature. In that sense, the models are representative of a single-phase medium, not of a multi-phase medium. Model predictions are scaled to our adopted 13 CO and 12 CO abundances (section 2.3.2). X CO is computed as the inverse of the 12 CO(1-0) emissivity, ǫ 12 , where ǫ 12 = I 12 /N 12 , I 12 is the integrated intensity of 12 CO(1-0) and N 12 is the column density of 12 CO. We vary the optical depth of the 12 CO(1-0) Figure 7. Correlation between ℜ and the 12 CO(2-1)/ 12 CO(1-0) ratio (top) and between ℜ and the X 12CO conversion factor (bottom) predicted by the RADEX grid of models from Leroy et al. (2017a). Observations are plotted in grey. The models show the expected trend, in LTE with filling factors of unity and fixed abundances, of X 12CO being inversely proportional to ℜ (see section 3.7). Optical depths τ 12 are varied from 3 (dark-color symbols) to 10 (light-color symbols) in steps of 0.1 dex, while the ratio τ 13 /τ 12 is fixed to 0.03. Left panels: model predictions for kinetic temperatures of T kin = 15 K (blue diamonds), T kin = 25 K (purple asterisks), and T kin = 35 K (red triangles). The mean volume density is set to 10 3 cm -3 . Right panels: model predictions for mean volume densities of n H = 10 2 cm -3 (blue diamonds), n H = 10 3 cm -3 (purple asterisks), and n H = 10 4 cm -3 (red triangles). The kinetic temperature is set to 25 K. line, τ 12 , from 3 to 10 in steps of 0.1 dex. The ratio of the 13 CO(1-0) and 12 CO(1-0) line optical depths, τ 13 /τ 12 , is set to 0.03 and the effects of varying this value are discussed below. We assume a lognormal distribution of H 2 volume densities with a width σ = 0.8 dex. Several values of mean volume density (n H = 10 2 , 10 3 , 10 4 cm -3 ) and temperature (T kin = 15, 25, 35 K) are tested. The predicted quantities that we investigate and the trends with density are somewhat sensitive to the selected width of the distribution, but not enough to change the results reported below.
Figure 7 shows predictions of ℜ, 12 CO(2-1)/ 12 CO(1-0), and X CO . We discuss the effects of optical depth and density separately, but these could be correlated. As expected, the 12 CO(2-1)/ 12 CO(1-0) line ratio (top panels) increases for increasing temperature and density, and it is rather constant for increasing τ 12 (same symbols, different shade). Concerning X CO (bottom panels), ℜ decreases and X CO increases with increasing τ 12 because the 12 CO emissivity goes down.
ℜ is inversely proportional to X CO (see also equation 6 in section 3.7), as long as τ 13 remains small. The models predict an increase of X CO or decrease of 12 CO emissivity -for a fixed τ -with increasing temperature. Increasing the temperature reduces the opacity per unit mass of 12 CO. But to keep τ constant, the increase in 12 CO intensity is not as large as the increase in gas column density because 12 CO is optically thick, hence its emissivity goes down. The models predict only mild variations of ℜ and X CO with density: ℜ decreases most for a density of 10 3 cm -3 and X CO increases linearly for the densities considered. This can be understood as the emissivity of 12 CO peaks at densities of 10 2 cm -3 (high ℜ and low X CO ) while that of 13 CO peaks for densities 10 3 cm -3 (Leroy et al. 2017a).
Overall, we find that the models can only partially reproduce the observations. Density and temperature have a weaker impact on ℜ predictions than optical depth. The range of ℜ values found in the observations is covered by the models, but the range of 12 CO(2-1)/ 12 CO(1-0) ratios and X CO values is not. While one can imagine that a broader range of temperatures or densities (or a different opacity ratio) than tested could explain the range of observed 12 CO(2-1)/ 12 CO(1-0) ratios, the trends of ℜ with X CO require additional modifications. Focusing on the isotopic abundance ratio (which is proportional to the optical depth ratio in LTE), here we discuss qualitatively which parameter space in the models is compatible with our observations:
(i) Galaxies with values in the disc of ℜ ≃ 10 and α CO ≃ α CO,MW (NGC 0628, NGC 3184, NGC 4254, NGC 5194), require the models in Fig. 7 to shift by a factor of 2 either to the right, which can be accomplished by decreasing the 12 CO/H 2 abundance (keeping the isotopic ratio constant), or up, which can be accomplished by increasing the isotopic abundance ratio (keeping the 12 CO/H 2 abundance constant).
(ii) The trend of higher ℜ values and lower α CO , noted above for NGC 3627, NGC 6946, and to a lesser extent, NGC 4321, can naturally be explained with the models having low τ 12 values (≃ 4). In their centres, the lowest α CO values still require the models to shift by a factor of 2 -3 either to the left, which can be accomplished by increasing the 12 CO/H 2 abundance (keeping isotopic ratio constant), or down for very low τ 12 values ( 12 CO becoming optically thin), which can be accomplished by lowering the isotopic abundance ratio (to ≃ 20).
(iii) The low ℜ and low α CO values in the centre of NGC 5055 can be accounted for with the same modifications as for (ii), but for higher τ 12 values.
In summary, for the discs of most of galaxies in our sample, the range of ℜ and X CO values observed in can be explained by variations in optical depths (for a fixed optical depth ratio or isotopic abundance ratio). However some values, especially in galaxy centres, require a change by a factor of 2 -3 in the optical depth ratio (or isotopic abundance ratio), which seems reasonable. We should keep in mind that physical conditions may also be changing within the model grid (i.e. a multi-phase model would be more representative).
Since 13 CO(1-0) emission remains optically thin over larger parts of molecular clouds, it is often advocated as a more accurate tracer of the molecular gas mass than 12 CO(1-0) in the intermediate-density regime (n H ≃ 10 3 cm -3 ). In this section, we present the optical depth of the 13 CO(1-0) line (τ 13 ), 13 CO column densities (N 13 ), and H 2 column densities (N(H 2 ) 13CO ) assuming local thermodynamic equilibrium (LTE). We compare the H 2 column densities to those obtained from 12 CO(1-0) emission. We opt for a simple framework as an alternative to the non-LTE models such as those presented in the previous section to have more flexibility on the choice of conditions/parameters and, in particular, explore the possibility of 12 CO and 13 CO having different filling factors. Ultimately, more lines and transitions would be needed for a complete, multi-phase modeling.
Given LTE, the equations of radiative transfer give the general expression for the observed brightness temperature of a line (T obs ):
where η bf is the beam filling factor, J ν is the line intensity given by the Planck function, T ex is the excitation temperature, T bg is the background temperature equal to 2.7 K, and τ is the optical depth of the line. To reduce this equation, we make the following assumptions for the 13 CO(1-0) and 12 CO(1-0) lines:
τ 12 > 1 while τ 13 ≤ 1, η bf,12 = 1 ( 12 CO fills the beam) while η bf,13 is let free.
The optical depth τ 13 and the column density N 13 are given by:
J ν,12 (T ex,12 ) -J ν,12 (T bg ) J ν,13 (T ex,13 ) -J ν,13 (T bg ) × I 13 I 12
(3)
where I is the integrated CO line intensity (in K km s -1 ), obtained by integrating the line profiles which have similar shape and width for both CO lines. Equation 2 corresponds to equation (15.36) from Wilson et al. (2012). At the adopted T ex values (≥ 20 K), we are in the Rayleigh-Jeans regime. The temperature term is an approximation of the partition function of CO with all energy levels populated in LTE. This simplification is valid for T ex ≥ 8 K. If τ 13 is always small, N 13 is an increasing function of T ex because upper states get preferentially populated. We also assume that the continuum is weak (no IR pumping). We refer to Jiménez-Donaire et al. (2017b) for details of the calculations. In case the beam filling factor of 13 CO is lower than unity (hence lower than that of 12 CO), the quantities τ 13 and N 13 represent local values as opposed to beam-averaged values. To be able to compare H 2 column densities derived from 12 CO(1-0) and from 13 CO(1-0), we shall then consider beam-averaged quantities. The beam-averaged column density of H 2 can be deduced from 13 CO(1-0) with:
From these formulae, for τ 13 < 1, the 12 CO-to-H 2 conversion factor will depend on ℜ, which captures the optical depth of the 12 CO line. In this case: To start with (Case 1), we assume thermalization of the lines and make the most simplifying assumptions, i.e. that both lines have the same excitation temperature, a beam filling factor equal to unity, abundances equal to the fiducial values, and X CO equal to the Milky Way value. The assumed parameters as well as the results for the 13 CO optical depths (τ 13 ) and column densities (N 13 and N(H 2 )) are reported in Table 4. We give average values of τ 13 and N 13 and their dispersion within the maps, for entire galaxies, their centres, and their discs. Since our sample of galaxies display quite a narrow range of ℜ values, the average values for τ 13 are very similar in all galaxies and around 0.07 -0.11. Those are beamaveraged values and correspond to τ 12 ≃ 6 (for an isotope abundance ratio of 60) and to a 12 CO effective critical density of n crit = n crit,thin /τ ≃ 2 × 10 2 cm -3 . For N 13 , we find values in the range 0.4-2× 10 15 cm -2 with a factor of two dispersion in the maps. Those correspond to 12 CO column densities of 0.2 -1.2 × 10 17 cm -2 and are averages over large beams, not local quantities. The implied optically thin 13 CO-to-H 2 conversion factor is 8 × 10 20 cm -2 (K km s -1 ) -1 .
We find that values of N(H 2 ) 13CO are systematically lower by a factor of ∼3 compared to N(H 2 ) 12CO values. If our assumptions hold, either 13 CO traces poorly the total H 2 column densities, or 12 CO overpredicts the total H 2 column densities. Similar offsets in column densities were found by Meier et al. (2001); Meier & Turner (2004); Meier et al. (2008). Taking the 12 CO-based estimates as reference, we can speculate about mechanisms to increase the 13 CO-based estimate. To increase N(H 2 ) 13CO by a factor of 3, we would need either an average excitation temperature of 60 K, an isotopic abundance ratio of 180, or a difference in beam filling factors of about 3. The maximum values allowed for each of those parameters and for each galaxy are reported in Table 4. We note that the beam filling factor can reflect variations in both abundances and emissivities of the lines but here we aim to separate the two effects and we consider that beam filling factors are dominated by emissivity variations. Although uncertainties associated to the adopted abundances are large, and temperatures in the galaxy centres may approach the required 60 K value (e.g., for NGC 6946, Walsh et al. 2002;Meier & Turner 2004), these requirements seem unlikely to be valid throughout the discs and in all discs of our sample. Moreover, if the temperature is indeed higher than the typical temperature of Galactic molecular clouds, then a Galactic conversion factor would also no longer apply for 12 CO. Hence, the discrepancies between H 2 column densities derived from 13 CO(1-0) and those derived from 12 CO(1-0) are most likely due to differences in beam filling factor of the two lines. Presumably, 12 CO(1-0) emission fills the beam and traces a more diffuse phase (n < 10 3 cm -3 ) while 13 CO(1-0) emission is confined to a denser phase. In our sample of galaxies, the H 2 column densities can be explained by a filling factor fraction of diffuse versus dense gas that varies between 2.2 and 3.5, with a possible trend of higher fraction in galaxies with more clumpy star formation (e.g., NGC 0628, NGC 6946) and lower fraction in galaxies with strong spiral modes (e.g., NGC5055, NGC 5194).
The discrepancy between the H 2 column densities derived from 13 CO(1-0) and those derived from 12 CO(1-0) may be due to oversimplifying assumptions, probably because, in our large beam size, we are probing clouds with a range of properties (densities, temperatures, opti-cal depths; Szűcs et al. 2016;Leroy et al. 2017a). For example, abundance variations dependent on column density can change N(H 2 ) 13CO by a factor of 2 -3 in Galactic molecular clouds (Goldsmith et al. 2008), though not sufficiently enough to reconcile LTE masses with the higher virial masses (Heyer et al. 2009). In the limit of co-existing 12 CO and 13 CO emission, and for a given set of physical conditions as chosen for Case 1, considering a sub-beam density distribution (section 3.6) would not be sufficient to reconcile column densities (models predict X CO below the Milky-Way value). In numerical simulations of realistic molecular clouds, Szűcs et al. (2016) investigate how standard methods ( 13 CO and LTE, 12 CO and X CO , the virial method) perform in recovering the true molecular mass/column density. They find that the 13 CO method is the worst predictor and systematically underpredicts the true mass by a factor of 2 -3 because of chemical and optical depth issues. Throughout our maps and with a beam size larger than the typical size of a molecular cloud, τ 13 remains small, below 0.2, and the variations from beam to beam are not significant enough to change N(H 2 ) 13CO . It is possible that our N(H 2 ) 13CO values are slightly underestimated because 13 CO is not completely optically thin and hides dense gas within the beam, although this explanation is unlikely. Indeed, in Galactic molecular clouds, τ 13 can be locally enhanced on local, subparsec scales (e.g., Kramer et al. 1999Kramer et al. , 2004;;Jakob et al. 2007), but the impact on the global H 2 column densities derived is marginal (e.g., Wong et al. 2008). Our low values for τ 13 are indicative of dense gas mixed with large amounts of diffuse gas within our (kpc-scale) beam.
12 CO(1-0) being optically thick leads to line trapping with the effect of lowering the critical density of the 12 CO line. 12 CO can emit strongly in low-density gas that makes up a large part of molecular clouds, and where 13 CO emission is weak because it is sub-thermally excited (e.g. Goldsmith et al. 2008;Leroy et al. 2017a). In that case, the fact that the bulk of the 12 CO and 13 CO emission does not trace the same gas can naturally explain the lower N(H 2 ) 13CO values compared to the 12 CO-based estimates. To reasonable approximation, for optically thin emission the emissivity per molecule scales with the density for n < n crit , and is constant for n > n crit . The CO critical density in the optically thin case is ∼ 2, 000 cm -3 , so a factor of 3 correction in column densities would be expected if the volume density of the gas dominating the 13 CO emission were ∼ 700 cm -3 . Conversely, one would only expect that the 13 CO would produce LTElike levels of emission if the density were above ∼ 2, 000 cm -3 , which is clearly a large density for the bulk of molecular clouds. The 12 CO(1-0) emission is presumably more extended than the 13 CO(1-0) emission, as seen in resolved (pcscales) studies (e.g., Pety et al. 2013). In the Milky Way, Roman-Duval et al. (2016) quantify the fraction of diffuse gas (gas detected in 12 CO(1-0) but not in 13 CO(1-0)) and dense gas (gas detected in both 12 CO(1-0) and 13 CO(1-0)). In terms of luminosity, they find that half of the gas is diffuse and half is dense in the outer disc, while in the inner disc, most of the gas (≃ 80 per cent) is dense.
Non-LTE considerations, such as different excitation temperatures, abundances, and filling factors for the two lines (mimicking a two-phase model) are explored in the following (Case 2).
T l ex,13 Notes. Column 1: galaxy name. Columns 2-3: assumed excitation temperatures of the 12 CO(1-0) and 13 CO(1-0) lines. Level population of the two levels of interest corresponding to the level population as predicted by LTE at the specific temperature. Column 4: assumed ratio of beam filling factors. Column 5: assumed isotope abundance ratio. Column 6: assumed value of the 12 CO(1-0)-to-H 2 conversion factor, normalized to the Milky Way value. The values are based on the results from Sandstrom et al. (2013) when available and otherwise (as indicated by the symbol * ) motivated by Bolatto et al. (2013). Column 7: mean optical depth of the 13 CO(1-0) line and dispersion in the map in parenthesis. Column 8: column density of the 13 CO molecule and dispersion in the map in parenthesis.
Column 9: ratio of H 2 column densities obtained from 13 CO(1-0) emission with the conditions given in columns 3-5 and from 12 CO(1-0) emission using the conversion factor given in column 6. Columns 10-11-12: maximum values allowed for the physical conditions such that column 9 equals to unity.
Results for τ 13 , N 13 , N(H 2 ): motivated choice of conditions (Case 2)
In the second case, we allow for the lines to be nonthermalized and we adopt reasonable physical parameters for the centre and discs of our sample based on Galactic studies and on results from section 3.5. In the centres, the temperatures and filling factor of dense gas are assumed higher than in the discs, and the isotopic abundance is assumed lower. The assumed parameters and resulting quantities are listed in Table 4. By allowing the excitation temperatures and the beam filling factors of the two CO lines to differ, the resulting τ 13 values generally increase with respect to Case 1. τ 13 is found around 0.14 -0.34 in the centres and around 0.26 -0.41 in the discs. N 13 has similar values than in Case 1 for the discs, but larger values for the centres, around 1 -11 × 10 15 cm -2 , because of the higher temperatures used.
When we adopt more realistic conditions for the centres and discs of our sample, the values of N(H 2 ) 13CO get closer to the values of N(H 2 ) 12CO than in Case 1. The results are also very sensitive to the adopted X CO factors. The X CO values found by Sandstrom et al. (2013), based on dust emission, are generally lower than the standard Milky Way value. For the discs, the need for a higher filling factor of diffuse versus dense gas (by a factor of about 2) still persists in most galaxies. In the centres, the effect of a low X CO value is somewhat compensated by adopting a lower isotope abundance ratio, without necessarily requiring a change in the beam filling factors. However, in NGC 3627, NGC 4321, and NGC 5055, the X CO values are so low that this leads to 13 CO predicting twice more H 2 than 12 CO. In those galaxies, reconciling N(H 2 ) 13CO with N(H 2 ) 12CO requires either lower temperatures than assumed, or lower abundance ratios, or larger filling factors for 13 CO than for 12 CO. Both conditions on the temperatures and filling factors seem unrealistic and the need for changing abundances is in line with results from section 3.5. It is also possible that uncertainties on X CO in those galaxies dominate here.
The main trends of ℜ as a function of galactocentric radius that we find are: (1) a slight increase from centre to outer disc (by a factor of 2 at most); (2) offsets between galaxies (by a factor of 2 at most); (3) a high value in the centres of NGC 3627 and NGC 6946. Causes for ℜ variations are essentially linked to isotopic (or isotopolgue) abundances or gas physical conditions. They are discussed extensively for external galaxies in, e.g., Paglione et al. (2001); Tan et al. (2011); Danielson et al. (2013); Davis (2014) and briefly examined here:
Changes in isotope abundance due to stellar nucleosynthesis (e.g., Henkel & Mauersberger 1993;Casoli et al. 1992). After a recent burst, one expects the 12 C and 18 O abundances to be enhanced relative to the 13 C abundance (e.g., Meier et al. 2014;Sliwa et al. 2017b) and ℜ to increase, although the actual abundances are sensitive to the star-formation history and chemical effects within a galaxy. In the Milky Way, observations show an increase of the 12 C/ 13 C and 16 O/ 18 O abundance ratios with galactocentric radius (e.g., Milam et al. 2005) that can be reproduced by time-dependent models (e.g., Romano et al. 2017). A mild radial increase of the isotope abundance ratio (similar to or shallower than the gradient observed in the Milky Way) would naturally explain the steadily increasing ℜ profiles of some of the galaxies in our sample.
Changes in abundance due to selective photodissociation (e.g., Bally & Langer 1982;Visser et al. 2009). Since 13 CO is less abundant than 12 CO and their abundances are coupled, one expects 13 CO to be preferentially photo-dissociated under hard radiation fields. The effect is less clear for C 18 O as it can be formed in a separate way (Bron et al. 2017). Simulating a suite of molecular clouds, Szűcs et al. (2014) showed that the effect of selective photo-dissociation on the 12 CO/ 13 CO abundance ratio is minimal. Jiménez-Donaire et al. (2017c) also argue that shielding from dust and H 2 dominates over self-shielding of CO molecules. Therefore, we do not consider selective photo-dissociation as a dominant effect.
Changes in abundance due to chemical fractionation (e.g., Liszt 2007). This process can enhance the abundance of 13 CO by a factor of ∼2-3 at low temperature and low optical depth in molecular clouds (Szűcs et al. 2014). It is believed to have a non-dominant effect on large (kpc) scales in the discs of galaxies (Paglione et al. 2001), though our C 18 O(1-0) observations do not rule out fractionation as an important effect (Jiménez-Donaire et al. 2017c). The lowest ℜ ratios are found in the centres of a few galaxies and ℜ anti-correlates with dust temperature and Σ SFR for those galaxies, the inverse of what is expected from chemical fractionation. Hence we consider chemical fractionation as an unlikely explanation for the lowest ℜ values.
Changes in physical conditions such as gas density, temperature, and opacity (e.g., Pineda et al. 2008;Wong et al. 2008). Theoretically, we expect 13 CO(1-0) emission to trace denser and cooler gas than 12 CO(1-0). As explored in section 3.4, temperature/excitation may drive some of the lower ℜ values observed but it is not the only condition affecting ℜ. The presence of diffuse emission or increased turbulence (for example, due to a stellar bar) would lower the optical depth of 12 CO(1-0) and boost its emission relative to 13 CO(1-0) (high ℜ). Using several J transitions of 12 CO and 13 CO, Israel (2009a,b) model the centres of galaxies that have active nuclei with two ISM components, one hot and tenuous component with low optical depth, and one cooler and denser component. A different mixing of such two components (with different physical conditions or in different proportions) may explain the range of ℜ values observed in our sample of galaxies. The centre of NGC 3627 displays the largest 12 CO(1-0) line widths (noting that this includes resolution and inclination effects, and its strong stellar bar could add significant non-circular motions as well). In the centre of NGC 6946, Wu et al. (2017) also find larger velocity dispersions by analysing cloud-scale 12 CO(2-1) observations, although Meier & Turner (2004) do not find a correlation between line width and ℜ.
We conclude that changes in the isotopic abundance due to nucleosynthesis or/and changes in the gas physical conditions, such as a different mixture of dense, cold gas and diffuse, warm gas that impacts on the mean opacity, could account for the trends with radius and offsets between galaxies. In the galaxy centres, turbulence/bars could account for the highest ℜ values that we observe, but their low α CO values still argue for changes in isotopic abundances (section 3.5). To disentangle the effects of abundance and physical conditions on ℜ, at least another 13 CO transition would be needed for a future, non-LTE modeling analysis.
4.2 12 CO, 13 CO, and the X CO conversion factor
In this section, we aim to understand under which circumstances 13 CO(1-0) can be used to improve estimates of the molecular gas mass.
Galaxy centres are particularly interesting because they show some of the strongest variations in observations (including in ℜ) and in physical properties in spiral galaxies. Several studies have found that the X CO factor that converts 12 CO(1-0) intensity to molecular gas surface density is lower in the centres of galaxies (e.g., Oka et al. 2001;Zhu et al. 2009;Sandstrom et al. 2013). This could be due to factors such as an excess of 12 CO(1-0) emission from increased temperatures and large velocity dispersions lowering the optical depth of 12 CO (e.g., Bolatto et al. 2013). In our sample, high 12 CO(2-1)/ 12 CO(1-0) ratios, tracing higher temperatures/densities, are indeed observed in some galaxy centres (Leroy et al. 2009(Leroy et al. , 2013)). The lowest X CO factors are also We also report the Spearman's rank correlation coefficients, their uncertainty, and their significance (in parenthesis), which are measured for a given CO line on all data points above a signal-to-noise ratio of 5 and measured on galaxy centres only.
found for galaxies with the highest ℜ (see section 3.5). In those cases, knowledge on both 12 CO(1-0) and 13 CO(1-0) emission can help reduce uncertainties on X CO if the later is unknown. However, we do not find a systematic increase of 12 CO(1-0) emission relative to 13 CO(1-0) in the galaxy centres. This implies that, at kpc-scales in normal star-forming disc galaxies, the issues encountered when using 12 CO(1-0) as a tracer of the molecular gas mass are not systematically solved by using 13 CO(1-0). Figure 8 shows the correlation between Σ SFR and the 13 CO(1-0) and 12 CO(1-0) intensities of all EMPIRE galaxies. Density of data points are plotted as contours and data points corresponding to galaxy centres are overplotted (filled circles, centres are defined by a cut in CO intensity and radial distance, see section 3.2). The strengths of the correlations as well as the contour shapes are similar for both lines. The correlation between 12 CO(1-0) and Σ SFR is marginally tighter than with 13 CO(1-0) when considering all data points and it is the same when considering centre points only. In the brightest regions (Σ SFR ≥ 0.1 M ⊙ yr -1 kpc -2 ; corresponding to the centres of all galaxies except NGC 0628 and NGC 3184), the 13 CO(1-0) data points show slightly less scatter with Σ SFR than 12 CO(1-0), and the star-formation efficiency from 13 CO(1-0) also has less scatter. In those bright regions, the distribution of 13 CO(1-0) data points tend to follow a superlinear relation with Σ SFR while the distribution of 12 CO(1-0) data points appears to be bimodal with points (corresponding to NGC 6946) following a linear relation and points (corresponding to NGC 5194) following a super-linear relation. Hence 13 CO does not appear to be a more stable tracer of the SFR than 12 CO in our sample of disc galaxies, except at the high Σ SFR end.
The analysis of H 2 column densities under the LTE as- X CO,13-total X CO,13-centre X CO,13-disc average all 10.0 (0.29) 3.6 (0.25) 10.8 (0.29) average nuc. (a) 4.9 (0.10) 2.2 5.3 (0.11) Average values, in units of 10 20 cm -2 (K km s -1 ) -1 . The standard deviation is given in parenthesis and is in dex. The averages are measured as the mean of the logarithmic values of X CO,13 . (a) Galaxies with nuclear starbursts (NGC 2903, NGC 3627, NGC 4321, and NGC 6946).
sumption in section 3.7 has shown that 13 CO(1-0) systematically underpredicts masses of H 2 , both in all pixels of our maps and if we consider integrated measurements. The behavior of the central pixels is obscured by the behavior of the more numerous disc pixels when looking at galaxies in their entirety. We have attributed the low H 2 column densities derived from 13 CO(1-0) in the discs to, mainly, the presence of diffuse gas that is traced by 12 CO(1-0) or dust emission but not by 13 CO. General methods to convert the dust and 12 CO(1-0) emission of galaxies to mass, however, are not calibrated on those diffuse phases and may overestimate the molecular gas mass. In turn, the presence of COdark gas could compensate somewhat for this overestimation (see Liszt et al. 2010;Liszt & Pety 2012). While 13 CO fails in retrieving the bulk molecular gas mass, its emission may vary less than 12 CO in the galaxy centres where the fraction of dense gas is larger and 13 CO traces better the intermediate-density regime (∼ 1, 000 cm -3 ).
The X CO factor from 13 CO(1-0) emission
We can derive an empirical X CO factor from 13 CO(1-0) emission by dividing dust-based H 2 column densities (see section 3.7.3) by the 13 CO(1-0) intensity: X CO,13 = N(H 2 , dust)/I 13 or X CO,13 = X CO,12 × ℜ. Figure 9 shows values of X CO,13 (absolute and normalized to the average in each galaxy) as a function of distance to the centre for all galaxies. Averaged values in different parts of the galaxies are reported in Table 5. For entire galaxies, we find an average value of X CO,13 = 1.0 × 10 21 cm -2 (K km s -1 ) -1 over our sample. If one were to consider a galactic conversion factor appropriate, then X CO,13 can just be obtained by scaling X CO with ℜ, and the average X CO,13 factor to use would be 2.2 × 10 21 cm -2 (K km s -1 ) -1 . As for X CO,12 (Sandstrom et al. 2013), X CO,13 has lower values in the galaxy centres, by a factor of ∼ 3, and increases with distance. The scatter is large for absolute X CO,13 values and reduces to less than a factor of 2 for X CO,13 values normalized to each galaxy average. Such trends are also observed but less pronounced with Σ SFR and the stellar surface density decreasing with increasing X CO,13 because those quantities have higher values in galaxy centres. Sandstrom et al. (2013) found that the 12 CO(1-0)-to-H 2 conversion factor, X CO,12 , varies more within the galaxies dominated by a central starburst: it is about one order of magnitude lower than the Galactic value in the centres and three times lower on average. Such low con-Figure 9. Top panel: Empirical X CO factor showing values lower than average in the galaxy centres and increasing with distance from the centre. X CO,13 is derived from dust-based H 2 column densities and 13 CO(1-0) emission, in units of cm -2 (K km s -1 ) -1 . The black open squares and curve show the mean values and dispersion of all data points within bins of width 0.1 r 25 . The horizontal dotted line indicates the typical value of the 12 CO-to-H 2 conversion factor in the Milky Way. We also report Spearman's rank correlation coefficients and their significance in parenthesis. Middle panel: Same as above, normalized to the mean X CO,13 value in each galaxy. Bottom panel: Same as the top panel for 12 CO(1-0) emission.
version factors can be found for extreme starbursts like LIRGs (Downes & Solomon 1998;Kamenetzky et al. 2014;Sliwa et al. 2017a). The galaxies dominated by a central starburst in Sandstrom et al. (2013) that are overlapping with our EMPIRE sample are NGC 3627, NGC 4321, and NGC 6946 (bottom panel of Figure 9). Those galaxies have higher ℜ values in their centres and on average. The top panel of Figure 9 shows that they have lower X CO,13 than the other galaxies on average, and slightly less variation from centre to outer disc than with 12 CO(1-0) because of their declining ℜ profiles. Within galaxies, X CO,13 and X CO,12 show a similar dispersion of ∼0.2 dex, but there are still large variations from galaxy to galaxy. Taking all galaxies together, the dispersion per radial bin is about 0.30 dex for X CO,13 and 0.35 dex for X CO,12 .
Overall, our results suggest that the observed changes in X CO,12 cannot be explained purely by 12 CO optical depth effects. In that case we would expect X CO,13 to stay approximately constant, while X CO,12 changes. In contrast, our observations favor variations in X CO,13 as well, with lower values in the centres. This suggests that some combination of varying parameters (abundance, optical depth, etc.) is also affecting X CO,13 . We note that these conclusions are all tied to H 2 column estimates from the dust-based method, so they are subject to the same systematic uncertainties inherent in that technique.
We present new observations of the 13 CO(1-0) emission from the EMPIRE survey (Bigiel et al. 2016;Jiménez-Donaire et al. 2017a) and of the 12 CO(1-0) emission from follow-up programs obtained with the IRAM 30-m in 9 nearby spiral galaxies (NGC 0628, NGC 2903, NGC 3184, NGC 4254, NGC 4321, NGC 5055, NGC 6946;and NGC 5194 from PAWS, Pety et al. 2013). 13 CO(1-0) is detected at high signal-to-noise across the entire molecular disc in those galaxies. We summarize our results as follows:
-The integrated intensity of 13 CO(1-0) is on average 11 times fainter than that of 12 CO(1-0). The 12 CO(1-0)-to-13 CO(1-0) intensity ratio (ℜ) does not vary significantly within or amongst galaxies in our sample (at a resolution of ∼1.5 kpc), by a factor of ∼2 at most, and the spread in values is largest in the galaxy centres. On those spatial scales, the 12 CO(2-1)-to-12 CO(1-0) ratio varies between 0.3 and 2.
-We correlate ℜ with several physical quantities on global and resolved scales. We find no strong trend on global scales. On resolved scales, ℜ anti-correlates weakly to moderately with the 12 CO(2-1)/ 12 CO(1-0) ratio, the IR colour, Σ SFR , and the SFE (Σ SFR /Σ H 2 ). Anti-correlations are more visible for galaxies with lowest ℜ values in their centres (NGC 4254, NGC 5055, NGC 5194), which we interpret as a local temperature/excitation effect. However, other conditions such as density/optical depths must be at play to explain all ℜ behaviors.
-We find that galaxies with starburst-dominated nuclei (which are also barred galaxies) have higher ℜ values in their centres than the other galaxies in our sample. In those galaxies, the anti-correlation of ℜ with the α CO conversion factor from dust is compatible with optical depth effects. Changes in the isotope abundance ratio from stellar nucleosynthesis due to the recent central bursts are also possible.
-Assuming LTE and fixed abundances, we compute optical depths and column densities for the 13 CO(1-0) line. Beam-averaged optical depths are on the order of 0.1. We find that the H 2 column densities derived from 13 CO(1-0) are systematically lower by a factor of 2 -3 than those derived from 12 CO(1-0). Those discrepancies can be mainly explained by non-LTE effects (significant diffuse phase in the discs where 13 CO is sub-thermally excited) or departure of abundances from nominal galactic values.
-Assuming H 2 column densities from dust emission, we calculate an empirical 13 CO(1-0)-to-H 2 conversion factor. The average value found in our sample is 1.0 × 10 21 cm -2 (K km s -1 ) -1 . It increases from galaxy centre to outer disc with similar scatter (∼0.3 dex) than the 12 CO(1-0)-to-H 2 conversion factor (based on 12 CO(2-1) data; Sandstrom et al. 2013). Except in the centres, 13 CO does not appear as a more stable tracer of the molecular gas mass or star-formation rate than 12 CO(1-0) in normal star-forming disc galaxies.
Despite being optically thin, the use of 13 CO(1-0) emission as a tracer of total mass seems to be limited in normal star-forming disc galaxies, due to the presence of a significant diffuse phase and variations in physical conditions of the gas and abundances. To constrain those, systematic observations of at least one more 13 CO transitions will be an important next step. less well peaked, probably because the dynamic range in the ratios is not very large. Σ(SFR)/CO(1-0) and 1/ℜ: The same two tests as described above are performed for Σ(SFR)/CO(1-0) instead of CO(2-1)/CO(1-0). Here, we find for all galaxies that the peak of the distribution of the standard deviations in the simulated Σ(SFR)/CO(1-0) values is systematically lower by a factor of ∼2 compared to the observed standard deviation. Hence the scatter in the observed Σ(SFR)/CO(1-0) values is physical.
The second test indicates that the correlation coefficients and slope values between ℜ and Σ(SFR)/CO(1-0) reported in Fig. 4d are also robust. The distribution of the slopes is more peaked than with CO(2-1)/CO(1-0).
MNRAS 000, 1-25 (2018)
http://www.astro.umd.edu/~bolatto/STING/ MNRAS 000, 1-25 (2018)
This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 000, 1-25(2018)