[C i](1–0) and [C i](2–1) in Resolved Local Galaxies^*

The Beyond the Peak (BtP) project used the SPIRE spectrometer (Griffin et al. 2010) on the Herschel Space Observatory (Pilbratt et al. 2010) to target 22 local star-forming galaxies taken from the KINGFISH project sample (see Table 1). Galaxies were chosen to represent a range of gas conditions, star-formation rates and galaxy masses while maintaining a central infrared surface brightness of 9 × 10⁻⁶ W m⁻² sr⁻¹ or brighter in a 40''aperture. Table 1 documents the physical resolution achieved for each galaxy and whether each galaxy's [C i] emission is not detected (below a S/N of 3), detected only in the center or detected in multiple resolution elements (resolved). Note that the physical resolution is computed based upon 40''extraction apertures; the actual spatial resolution is either 37'' for the longer-wavelength SLW bolometer data and 19'' for the shorter-wavelength SSW bolometer data. Additional properties of these galaxies such as distance, metallicity, stellar mass and star-formation rate are documented in Table 1 of Kennicutt et al. (2011).

The Herschel SPIRE FTS data were obtained and reduced as explained in detail in Pellegrini et al. (2013). To more reliably extract the faint, extended emission, many custom modifications were made to the Herschel pipeline. A semi-extended source correction was applied to Level 1 data, after which spectra were extracted from 40'' diameter regions. We choose to extract on a rectilinear grid with beam centers separated by 40'' so that the regions are statistically independent. We then fitted a modified blackbody with a scalar offset to each SPIRE/FTS spectrum (omitting ranges with line emission) and removed this offset. Synthetic photometry in the three SPIRE photometric bands was then computed from these offset-corrected spectra. We flux-calibrated the spectra using a scale factor to match this synthetic photometry to that measured by the SPIRE photometer in identical apertures. After this flux calibration was applied, line fluxes were extracted by fitting Gauss-sinc profiles to the FTS spectra. This profile is appropriate because even at the highest spectral resolution at short wavelengths (300 km s⁻¹) nearly all the galaxies had unresolved spectral lines. The one exception was NGC 7331, which required convolution with a Gaussian profile due to physical spectral broadening from an inclined star-forming ring in the beam. Spectral lines were considered detected if they were above 3σ. The line intensities for [C i](1–0) and [C i](2–1) are tabulated in Table 2 for all 40'' diameter regions. The CO line intensities will be presented in a subsequent paper.

Ground-based CO data were compiled from multiple sources. CO(1–0) data came from the BIMA-SONG survey (only those galaxies with short-spacings; Helfer et al. 2003) or the Nobeyama CO Atlas of Nearby Spiral Galaxies (Kuno et al. 2007). CO(2–1) data came from the HERACLES survey (Leroy et al. 2013) and CO(3–2) data came from the JCMT Nearby Galaxies Legacy Survey and archival data (Wilson et al. 2012). The data cubes from these surveys were convolved to match our spatial resolution; line fluxes were extracted based upon velocity ranges from the high signal-to-noise CO(2–1) data. α_CO values are from Sandstrom et al. (2013), where available.

We fit single-component modified blackbody (MBB) functions to the FIR photometry for regions matching our spectroscopic apertures:

$$\begin{eqnarray}&&{I}_{\nu }=A{\nu }^{\beta }{B}_{\nu }(T),\end{eqnarray} \tag{ 1 }$$

where A is an amplitude, β is the emissivity, and B_ν(T) is the Planck function at temperature T. We use far-IR photometry at 70, 160, 250, 350, and 500 μm from the Herschel PACS and SPIRE maps obtained under the KINGFISH project (Kennicutt et al. 2011). The fits were performed using curve_fit function from the SciPy library; two example fits are shown in Figure 1. The fit parameters give an estimate of the cold dust temperature and the β parameter modifying the blackbody shape. We choose to use a single MBB component and drop the 24 μm photometry available from Spitzer, because the addition of this one more point, but three more free parameters (associated with a second MBB component) led to more degeneracies.

Zoom In Zoom Out Reset image size

We separate the disk emission from the central apertures in order to investigate the properties of [C i] emission in normal galaxy disks as well as in isolated galaxy centers. We plot the disk and central [C i](2–1)/[C i](1–0), CO(4–3)/[C i](1–0) and CO(7–6)/[C i](2–1) ratios against dust temperature, IRAC 3.6 μm intensity (for the disks; a proxy for stellar surface density), and central X-ray luminosity (for the centers) in Figure 2. All three ratios are chosen such that higher values reflect higher excitation conditions (temperature or density; see Table 3). The motivation for pairing these specific CO lines with the [C i] lines is that they are the closest in frequency and thus may be conveniently observed together, with similar angular resolutions. With 40'' diameter apertures separated by 40'' between aperture centers, we consider all non-central apertures to be dominated by disk emission.

Zoom In Zoom Out Reset image size

We look for correlations between both the disk and central [C i] ratios with dust temperature, as a proxy of local ISM excitation. Using a Kendall's tau statistic appropriate for censored data as implemented in the R package NADA (Nondetects And Data Analysis; Lee 2017), we can reject the null hypothesis of no correlation for the [C i](2–1)/[C i](1–0) and CO(4–3)/[C i](1–0) ratios for both disks and galaxy centers and the CO(7–6)/[C i](2–1) ratio for galaxy centers. The p-values of these correlations are shown in the upper right corner of the appropriate subplot in Figure 2. In all of these cases, the data show the expected correlation of more highly excited [C i] ratios with increasing dust temperature.

There are very few detections of the CO(7–6)/[C i](2–1) ratio within galaxy disks, but surprisingly, the detections tend to decrease in ratio as dust temperature increases in contrast to what is seen for the galaxy centers and reported elsewhere for LIRGs (Lu et al. 2017). We note that all of the high-temperature low CO(7–6)/[C i](2–1) ratio points are contributed by the same galaxy, NGC 4736. This galaxy's central near-IR spectrum lacks Brγ emission (Walker et al. 1988) and optical spectroscopy diagnoses its nuclear region as a 1 Gyr post-starburst with very little current star formation (Taniguchi et al. 1996). NGC 4736's dust is likely heated by this high-density post-starburst population instead of ongoing star formation. Furthermore, NGC 4736 has the strongest central [C ii] suppression (Smith et al. 2017), which is similarly interpreted to be due to a high-intensity, but softer (than ongoing star formation) stellar radiation field. Thus, NGC 4736's very low CO(7–6)/[C i](2–1) ratios at relatively high dust temperature may be explained by intense starlight with a relative lack of recent star formation.

Within galaxy disks, the IRAC 3.6 μm intensity is a proxy for the density of the radiation field contributed by old stars, which has been found to correlate with the [C ii] deficit (Smith et al. 2017). The relations in Figure 2 do not show any statistically significant correlations between these three [C i] ratios and the IRAC 3.6 μm intensity. Note that galaxy centers are excluded from this plot, because a low-luminosity AGN or central starburst can make even the 3.6 μm emission a less reliable tracer of old stars.

The galactic centers of the BtP galaxies are a mix of low-luminosity AGNs, nuclear starbursts and quiescent nuclei (Moustakas et al. 2010) in otherwise normal galaxies. Due to the resolution of our observations, the central extraction measures not only the true nucleus, but extended circumnuclear regions from a 0.9 to 6 kpc scale (see Table 1). Line ratios for these central regions are plotted against their nuclear (within 23) X-ray luminosity in Figure 2 in order to look for correlations with AGNs or star-formation power. The majority of these X-ray luminosities are from the compilation by Grier et al. (2011) of Chandra measurements, with additional Chandra measurements for NGC 1266 (Alatalo et al. 2014) and NGC 4536 (McAlpine et al. 2011). Note that these X-ray luminosities include both hard and soft emission (0.3–8.0 keV) and thus may have contributions from an AGN and/or from central star formation. No statistically significant correlations are observed between X-ray luminosity and any of these central [C i] ratios.

Both [C i] lines are indisputably correlated with all observed transitions of CO, as shown in the intensity versus intensity plots of Figure 3 (first and third rows of each column). In these plots, gray circles indicate detections, green left arrows indicate upper limits in CO, and blue down arrows indicate upper limits in [C i]. For clarity, data that are non-detections in both lines are not shown. We fit a power-law relation between the intensities (indicated with a dashed line) and the residual scatter around this relation (second and fourth row of each column).

Zoom In Zoom Out Reset image size

The power law we are interested in is written as follows:

$$\begin{eqnarray}&&I([C\,{\rm{I}}])\propto I{\left(\mathrm{CO}\right)}^{\gamma }.\end{eqnarray} \tag{ 2 }$$

Values of γ close to one imply that it is reasonable to convert directly from one line intensity to the other. Sub- or super-linear correlations imply a dependence on a physical effect that correlates with the line intensity of a region. Fitting in log–log space means the slope of a linear fit will be the desired power-law exponent, γ. However, the task of fitting such a line is complicated by the multitude of upper limits within these data sets. We follow advice from Feigelson & Babu (2012) and use the doubly censored Theil–Sen estimator given in Akritas et al. (1995) for fitting lines with doubly left-censored data as implemented in the R package NADA. We determine best slopes using this algorithm for each combination of CO line and [C i] line, which are shown with dashed lines in Figure 3. It should be noted that this algorithm does not return a best estimate for the y-intercept. Thus, we used least-squares to find a y-intercept estimate with a fixed slope at the NADA value, and ignored the censored data.

The calculated values of the power-law exponent are shown in Figure 4 for both [C i](1–0) and [C i](2–1) against the various J CO lines. Uncertainties on these slopes are determined by bootstrap resampling of the data. For [C i](2–1), the relationship with CO(4–3) is very nearly linear. For [C i](1–0), the relationship between [C i] and CO is always sublinear, which means that higher-intensity CO regions have higher CO/[C i] ratios, with a larger effect for lines farther from CO(4–3). For both [C i](1–0) and [C i](2–1), we expect the following physical reasons are at work. At high intensities of low-J CO, the majority of gas is cold, and unable to excite the transitions of [C i] with their higher upper level temperatures. Thus these regions have lower [C i]/CO ratios and result in sublinear correlations. At high intensities of high-J CO, there is more dense, warm gas, which is more effective at exciting the high critical density high-J CO lines. Again, this results in lower [C i]/CO ratios at the high intensity end and thus, again, a sublinear correlation.

Zoom In Zoom Out Reset image size

The intrinsic scatter in the relation can be determined using the statistical approach of Cappellari et al. (2013; their Equation (6)). The basic idea is to compute the reduced χ² statistic including a term for the intrinsic scatter that may be increased until the reduced χ² equals its expected value of 1. Using this approach, the intrinsic scatter appears minimal for the [C i](2–1) to CO(4–3) relation as shown in the bottom panels of Figure 4. This low scatter is a further indication of a strong relationship between [C i](2–1) and CO(4–3). Relations between [C i](1–0) and CO(2–1), CO(3–2) and CO(4–3) all have similar intrinsic scatter. Thus, we cannot make a strong statement about which CO line [C i](1–0) is best correlated with.

We first consider what the [C i] lines alone can tell us about the physical state of the gas. Unfortunately, there are only two [C i] lines, which means we must make some assumptions if we are to derive gas properties. A common assumption is that the neutral carbon is in LTE and the two [C i] lines are optically thin (e.g., Kamenetzky et al. 2014). In this case, the excitation temperature is simply calculated using the Boltzmann distribution for the [C i] level populations and Einstein A coefficients to determine the line fluxes from these level populations. Figure 5 documents these optically thin LTE temperatures for all galaxies where both [C i](1–0) and [C i](2–1) are detected. These temperatures are generally lower than those found for the more active galaxies (starbursts, AGNs, ultraluminous infrared galaxies (ULIRGs)) in the Kamenetzky et al. (2014) sample. Within their sample, only 6 out of 18 studied regions have LTE [C i] temperatures below 24 K, the highest temperature seen in the BtP sample. These optically thin LTE temperatures are consistent with the kinetic temperatures found for the [C i]-emitting gas for simulated molecular clouds at both a solar neighborhood interstellar radiation field and an ISRF 10 times stronger (Clark et al. 2019).

Zoom In Zoom Out Reset image size

However, the conditions for LTE may not be met, for instance, if subthermal excitation is important. Section 3.3 of Glover et al. (2015) discusses the problems they see when assuming LTE for the [C i] emission of their simulated cloud. If lines are (partially) optically thick, we also cannot simply derive the excitation temperature, as the observed flux will be lower than the emitted flux. We will later discuss how important these effects are for the BtP sample, if we use the gas properties derived from our two-component fit.

Next, we use models to help connect the [C i] lines to a physical gas phase. As CO and [C i] emission have been observed to be cospatial on parsec scales within Milky Way clouds (e.g., Figures 3 and 4 of Shimajiri et al. 2013) as well as on hundreds of parsec scales within resolved galaxies (e.g., Krips et al. 2016), we use the emission from both species to constrain the CO- and [C i]-emitting gas parameters.

LVG models are commonly used to describe the emission from molecular clouds because they provide an easily computable non-LTE approximation. To model the line intensities of a single molecular species, the parameters required are the kinetic temperature, the volume density of dominant collision partners, and the column density of the molecule of interest per unit velocity spread. When such models are applied to regions within galaxies, we assume those regions are made up of clouds with identical values of these parameters. An additional normalization parameter allows for the intensity of the region to be lower than if it were uniformly filled by such clouds. LVG models may also include multiple species, as long as either an abundance ratio is assumed, or there are enough lines to allow the abundance ratio to be fit.

Two-component LVG models are required to fit the [C i] and CO emission of these galaxies. Israel et al. (2015) noted that their lowest-excitation galaxies cannot be well fit by a single component. Schirm et al. (2017) document in their Appendix A why a one-component LVG model is not realistic for regions in M51. Similarly, we found unrealistically high gas kinetic temperatures and low densities when the BtP data were fit with a single LVG component.

A full two-component model requires 10 free parameters: density, kinetic temperature, column density per velocity interval, filling factor, and C to CO abundance for both components. In the very best cases within the BtP sample, only 9 lines were detected, therefore a more constrained model was necessary. We first tried a combination of a low-excitation template based upon the Milky Way and fit a high-excitation LVG model. However, this led to high reduced χ² values for many regions due to poor fits at the low-excitation end of the spectral line energy distribution. Thus, we instead use an empirically derived template for the high-excitation emission and an LVG model to account for the lower-excitation emission. The high-excitation template we use is based on the detailed empirical observations of an actively star-forming region in the LMC, N159W (Lee et al. 2016). Lee et al. (2016) analyze the sources powering the [C i] and CO emission in this region, acknowledging some is from star-formation-powered PDRs, while other emission must be powered by star-formation feedback related shocks. Instead of modeling these components separately, we simply use their integrated observed emission from the region as a template.

For the lower-excitation component, we use the LVG code RADEX to provide models of the emission (van der Tak et al. 2007). We constrained the models to n(H₂) = 10³ cm⁻³ because at higher densities the shape simply resembles the high-excitation template and we lose diagnostic power. Kinetic temperature is allowed to vary from 10 to 35 K in 5 K increments. N(CO)/Δv = 0.25, 0.5, 1, 2, 4 × 10¹⁷ cm⁻² km⁻¹ s are options for the column density per unit velocity. In addition, we allow the C⁰:CO abundance ratio to be 2:1, 1:1, 1:2, 1:4, 1:8, or 1:16. These choices are all based upon the "quiescent" and "very quiescent" models for [C i] and CO from Papadopoulos & Greve (2004), which summarize observations of the Milky Way, M31 and M33. They also cover the parameters found to fit the cold component of M51 (Schirm et al. 2017). The scale factor for each component is defined as the average number of Δv = 1 km s⁻¹ clouds along the line of sight and is hereafter denoted ${f}_{{ \mathcal L }}^{\star }$. While this appears to be an awkward definition, it prevents us from either having to fix the typical cloud width within our galaxies or from assuming the full velocity line width of the galaxy is the relevant line width for the LVG approximation. Based upon this definition, the N159W high-excitation region has ${f}_{{ \mathcal H }}^{\star }=1$ because it is a 10 km s⁻¹ wide cloud with a beam filling factor of 0.1 (Lee et al. 2016). Thus, the scale factor we apply to the high-excitation N159W template for a given region gives us that region's ${f}_{{ \mathcal H }}^{\star }$ scale factor.

For robust fits of the five parameters (T, N(CO)/Δv, C:CO, ${f}_{{ \mathcal L }}^{\star }$, ${f}_{{ \mathcal H }}^{\star }$) we decided to require detections in both [C i] lines and at least four CO lines, including one of CO(1–0) or CO(2–1). This cut results in 64 independent regions that are fit. The best fit is determined by minimizing the χ². We show the data and best fits for six example regions in Figure 6; histograms documenting the fit results for all regions are shown in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the adequate decomposition provided by the two-component fit, we indicate distributions of the low-excitation fraction of each line in Figure 8. For [C i], nearly all of the [C i](1–0) line is contributed by the low-excitation component. For [C i](2–1), the fraction varies significantly from region to region, although the median is close to 70%. For CO, there is a steady decrease in the contribution of the low-excitation component as J increases. We note that our range for the CO(2–1) high-excitation fraction corresponds nearly identically to the range Israel et al. (2015) report for their LIRG/ULIRGs CO(2–1) dense fraction (rightmost panel of their Figure 8). As expected, very little of the highest-J CO line intensities are contributed by the low-excitation component.

Zoom In Zoom Out Reset image size

These distributions represent the fraction of line emission from a low-excitation component for the brightest regions within a sample of normal, local spiral galaxies. By this analysis, [C i](1–0) and CO(1–0) are both reliable tracers of the low-excitation component for these regions. [C i](2–1) is less reliable, as even in these normal spiral galaxies there are cases where half of its emission is from a high-excitation component. The CO J_up = 2−5 lines have low-excitation fraction ranges of typically 30%–70%. So for these regions of spiral galaxies they are not dominated by one component or the other and indeed show a wide variation. While this limited two-component analysis is only a slight improvement on trying to account for the physical properties with a single LVG model, already we can see that considering multiple components is very important for the [C i]- and CO-emitting gas.

One of the goals of this modeling work is to explain the strong link found between the CO(4–3) and [C i](2–1) emission. In order to do this, we consider the range of line ratios produced by these two-component models. For various abundances and temperatures (shown on the x-axis; temperature increases from 10 to 35 K left to right in 5 K steps for each abundance), Figure 9 shows the line ratios produced by our two-component model. The different color circles represent different ${f}_{{ \mathcal H }}^{\star }/{f}_{{ \mathcal L }}^{\star }$ filling factor ratios. We depict the median ratio in orange, the ratio at 90% of the observed distribution in green and the ratio at 10% of the observed distribution in blue. We only show the N(CO)/Δv = 10¹⁷ cm⁻² km⁻¹ s models to avoid clutter. The number in the bottom right corner shows the multiplicative spread from the minimum to maximum line ratio for the models.

Zoom In Zoom Out Reset image size

For [C i](2–1), the models collapse to a small range at the CO(4–3)/[C i](2–1) ratio. Note that this also occurs for the four other column densities individually and all five column densities as a set (not shown). This collapse fits precisely with the linear correlation with minimal scatter observed for CO(4–3) versus [C i](2–1) in Section 3.2. While Figure 9 graphically shows the limited range in ratios, we may understand the narrow range physically as being due to the very similar upper level temperatures of 55 K versus 62 K for CO(4–3) and [C i](2–1), respectively. While their critical densities differ by a factor of about 20 and this is surely responsible for some of the differences seen in these ratios, the similar upper level temperatures are more important in these cases where density is not too far below critical.

However, the observed ratios exhibit an even narrower range than that populated by these "reasonable" models. The full range of observed ratios is shown by the gray band in all of these plots. So even if the models alone suggest a narrow range for the CO(4–3)/[C i](2–1) ratio, we find that real regions within galaxies have an even narrower range in this ratio, suggesting that all of our model parameter space is not explored by real galaxies.

For [C i](1–0), with its 24 K upper level temperature and lower critical density, the narrowest range of ratios should occur for CO(1–0)/[C i](1–0) based purely on the reasonable range models. But at least for the regions with [C i](1–0) detections, the analysis summarized by Figure 5 showed minimal scatter for the relationships between [C i](1–0) and CO(2–1), CO(3–2) or CO(4–3). This may be due to our detections in [C i](1–0) being biased toward the brightest (likely higher-excitation) regions of these spiral galaxies or the lower quality of the ground-based CO(1–0) data.