${L}_{\mathrm{CO}}^{\prime }/{L}_{\mathrm{FIR}}$ RELATIONS WITH CO ROTATIONAL LADDERS OF GALAXIES ACROSS THE HERSCHEL SPIRE ARCHIVE

All spectra used in the sample were reprocessed with HIPE⁷ developer's version 13.0.3849 and spire_cal_13_0, obtained from Rosalind Hopwood on 2014 September 30. This calibration corrects for rapidly changing telescope temperatures near the beginning of observation cycles, which has the largest effect on faint sources (Swinyard et al. 2014). Overall the calibration errors, even from earlier calibration sets, are within 6% for point sources and 7% for extended sources. The majority of the observations were performed in sparse sampling mode, for which we took the spectra from the central pixels (SLWC3, SSWD4). For those mapping observations in the intermediate or fully sampled modes, we extracted the spectrum from the pixel corresponding to the central coordinates of the map (those in Table 1). All SPIRE photometry was downloaded in 2014 September from the HSA (SPG v11.1.0) and not reprocessed.

The majority of the sources in our band are quite point-like compared to the SPIRE FTS beam, which varies from ∼45'' to 17''. Because the beam size is discontinuous between the upper frequency end of the SLW band and the lower frequency end of the SSW band, galaxies which are not point-like will show a notable discontinuity at this juncture. Even galaxies that are relatively small compared to the beam, but still not perfectly point-like, will show this discontinuity and require a correction to properly compare the emission across the SPIRE bandpass. This is necessary because we only use the central FTS detectors. An example is shown in Appendix.

We perform the same source/beam correction as described in full detail in Section 2.2 of K14. Briefly, we use the SPIRE Photometer Short Wave (PSW, 250 μm) maps, observed with 1932 beams, convolve them to larger beam sizes (Ω_b), and measure the new peak flux density. We compare this flux density to that of a b = 435 beam, which corresponds to the beam at the CO J = 4−3 transition. The ratio of the two flux densities, η_b,43.5, as a function of beam size, is between that expected for a point source (1) and a fully extended source (Ω_b/Ω_43.5). For each galaxy's unique distribution, for any beam size, we have a value of η_b,43.5 to refer the emission to a 435 beam. We divide the SPIRE spectrum by η_b,43.5 to refer the flux density at all wavelengths to that observed by a 435 beam. We also use these values to refer CO integrated flux values measured from other facilities with smaller beam sizes to the 435 beam (Section 2.4).

We apply an additional correction (also used in K14) to match the total flux density of the spectrometer with the photometer flux density. At the high frequency end, we match the total SSW flux density integrated over the photometer PSW bandpass, $\hat{F}(\mathrm{PSW})$, to that of the PSW photometer-integrated flux density at 435, F''(PSW) by multiplying the spectrum by X_SSW = F''(PSW)/$\hat{F}$(PSW). There are two photometer bands (PMW and PLW) which overlap with the SLW band, so we define a line that connects those two ratios, F''(PMW)/$\hat{F}$(PMW) and F''(PLW)/$\hat{F}$(PLW), and multiply the spectrum by that ratio as a function of wavelength, X_SLW (ν). This photometry correction step is often most significant in the SSW, which can overestimate the measured flux compared to the photometry. Some spectra had somewhat over-subtracted telescope emission, giving slightly negative flux densities, in particular at the lowest-frequency end. In these cases, no correction was applied to the SLW band to match the PLW photometry, which would use negative ratios. In a few cases, we also did not correct the SLW band if the ratios derived from the PLW and PMW were significantly discrepant (i.e., would produce a nonsensical SLW continuum). The spectra that were not corrected are marked in Table 1.

In order to compare the CO emission to L_FIR, we must also correct L_FIR to properly represent the same amount of emission as the CO within our beam. Similar to the procedure above, we convolve the SPIRE photometer maps at wavelength λ to the beam size of 435 and find the ratio of the peak flux density in Jy measured with such a beam (F_beam,λ) to the total integrated emission in the map (F_total,λ). Assuming L_FIR,beam = L_FIR,total × F_{beam,250 μm}/F_{total,250 μm}, and likewise for the 350 and 500 μm maps, we can determine the proper L_FIR,beam for comparison to the CO emission. The three photometers agreed well, and so we use the average of the F_beam,λ/F_total,λ ratios. Of the 232 observations with known redshift and available spectra, 118 have ratios of <0.8, and 42 have ratios of <0.5. We propagate the errors from the total measured integrated flux density through to the final measurement of L_FIR,beam.

The CO J = 4−3 to J = 13−12 lines, both [C i] lines, and the [N ii] 205 μm line are the brightest lines in the FTS spectra. To fit these, we start with the FTFitter code from the University of Lethbridge.⁸ Treating each detector (SLW and SSW) separately, the code fits a polynomial to the baseline, and then simultaneously fits unresolved lines at the expected frequencies of the lines listed above (given known redshifts). We place the lower limit of the total area of the line profile to be above 0; we do not expect any of these lines to be in absorption. We limit the line center to within ±500 km s⁻¹ of that expected from the redshift to allow for uncertainty in the velocity scale and physical differences in the gas kinematics. In wavenumbers, this is about 0.025–0.084 cm⁻¹ over the band, compared to the FWHM of the line profile of 0.048 cm⁻¹.

We manually inspected the resulting fits to determine if any lines were clearly resolved. This is most likely to be the case for the [N ii] line as its velocity resolution is the highest at the higher frequencies, and it is much brighter than the CO lines at similarly high frequencies which may be undetected. Resolved lines do not show the same characteristic ringing of the sinc function; the ringing is significantly lower, if not imperceptible, smeared out by the effective convolution of the emission line profile and instrumental profile. We refit the lines that meet this criteria as a Gaussian convolved with the instrumental line profile. In this case, the lines are barely resolved, thus Gaussian profiles are perfectly adequate (no more detailed velocity profiles can be determined from the FTS).

The fact that SPIRE utilizes a FTS introduces a special problem in the treatment of line fitting. The true measured quantity is the interferogram, or the interference pattern at the focus as a function of optical path difference (OPD) as the mirror of the interferometer moves linearly. The spectrum itself is the Fourier transform (FT) of this interferogram, which leads to two important consequences: (1) the wavelength bins are not truly independent, which many fitting routines assume, and (2) the resulting noise pattern closely resembles the FWHM = 0.048 cm⁻¹ sinc function line profile. The result is that the errors output by a least-squares fitter, like the FTFitter and the built-in HIPE Spectrum Fitter routines, may not be an accurate representation of the line flux uncertainty. Moreover, it can do an excellent job of fitting a "ripple" in the spectrum which, to the observer's eye, may not be particularly distinguishable from any other ripple nearby, other than that we expect the, e.g., CO line to correspond to the fitted ripple's wavelength, and no similarly strong lines to be adjacent in the spectrum.

Although the ideal situation would be to the fit the interferogram itself,F this is not a user-accessible option for SPIRE data considering the many calibration steps that occur in processing after the FT. Instead, we created a Bayesian analysis method to determine the probability distribution function (PDF) of the true line flux given the observed line flux, $P({f}_{\mathrm{true}}| {f}_{\mathrm{obs}})$, which is heavily influenced by the correlated noise pattern in the spectrum. The noise itself is difficult to accurately characterize, varying from observation to observation, and across the bandpass of a given observation. Therefore, instead of attempting to describe the correlated noise for our entire sample, we focus on the area around each individual (unresolved) line.

We describe the procedure briefly here, but show a more in-depth example with illustrative figures for NGC 4388 in the Appendix. This procedure is not used for lines that were manually identified as resolved, which already have a high signal-to-noise ratio (S/N). For each line, we input sinc profile lines of varying amplitudes f_true over the region ±2 cm⁻¹ from the line center (excluding the area immediately around any CO, [C i], or [N ii] lines) and then refit the spectrum. We compare the measured integrated fluxes, f_obs to the known input values, f_true. The PDF for our CO line is the distribution of the input fluxes that produced that particular measured flux value, a slice of the $P({f}_{\mathrm{obs}}| {f}_{\mathrm{true}})$ two-dimensional distribution that we created.

For high S/N line detections, the procedure replicates a Gaussian distribution of similar median and error (σ) as the parameters estimated by the FTFitter. This is because a very high amplitude line input, added anywhere in the spectrum, will return the same integrated flux value as we input (f_obs = f_true). However, a line with smaller amplitude may add constructively or destructively to the underlying ripple pattern of the spectrum, returning a higher or lower flux than input. A local variation from the detector's average baseline may also influence the final fitted value, which may shift the median value of the PDF of the line flux.

The procedure makes the most difference for the high-J CO lines; 60% of the 3σ detections of CO J = 13−12 from least-squares fitting were shown to have >3σ uncertainty. For all the CO lines in the SSW band (J_upper ≥ 9), this number is 44%. For the CO lines in the SLW band, it is only 15%. The numbers are the lowest for [N ii] (4%) because it is so bright, and for [C i] J = 2−1 and CO J = 7−6 (5%, 8%) because they lie in the lowest noise part of the spectrum and are relatively bright.

The results from this line fitting procedure are shown in Table 2. The median, −1σ, and +1σ values are derived from the points at which the cumulative distribution functions (CDFs) equal 0.5, 0.159, and 0.841, respectively. If −1σ/median is less than 3, a value for a 3σ upper limit is also shown (where CDF = 0.997).

The bandpass of the Herschel FTS starts around the CO J = 4−3 line, but the majority of the molecular mass in galaxies is cool and populates the lower rotational levels. We complement the line fluxes derived from the FTS with the CO J = 1−0, J = 2−1, and J = 3−2 lines available from ground-based observatories. Many of these galaxies have already been studied in the literature, particularly in large CO surveys.

For some galaxies, we also performed single-dish measurements using the ARO. Measurements of the CO J = 1−0 line were conducted with the 12 m dish on Kitt Peak in 2015 May, and those of CO J = 2−1 and J = 3−2 were conducted with the Submillimeter Telescope (SMT) located on Mt. Graham from 2014 November to 2015 February. At the 12 m, we used the ALMA Type 3 mm receiver with two 2 MHz backends in series, yielding a 2.6 km s⁻¹ channel resolution and an about 670 km s⁻¹ bandwidth. At the SMT, we used the ALMA Type 1.3 mm sideband separating receiver (for CO J = 2−1) and the 0.8 mm double sideband receiver (for CO J = 3−2) with the 1 MHz filterbanks in 2IF mode. Most observations were conducted with beam switching, except for highly extended sources which required position switching. Pointing and focus were checked periodically on planets or bright continuum sources.

Beam efficiency measurements were conducted with Venus and Jupiter⁹ using the procedure described in Schlingman et al. (2011). For CO J = 1−0, J = 2−1, and J = 3−2 we found η_mb of 0.86–0.89, 0.57–0.62, and 0.58–0.65, respectively. The beam sizes are approximately 55'', 32'', and 22'', respectively. The spectra were reduced, baseline subtracted, and converted to the T_mb scale using η_mb in CLASS. Finalized spectra were smoothed to approximately 20 km s⁻¹ bins from which the total integrated fluxes were derived.

All low-J lines utilized in this work, including the ones from ARO, are available in Table 3. As Table 3 shows, the low-J lines we use come from a variety of telescopes with different beam sizes. For subsequent comparison to Herschel CO lines, all line fluxes are referenced to the 435 beam size using the same ratios (η_b,43.5) as described in Section 2.2. Table 3 lists both the literature reported values (third column, in their original beam sizes and units) and the 435 referenced fluxes in Jy km s⁻¹ we use in our analysis (ninth column).

As mentioned in Section 2.2, all fluxes including L_FIR are referred to the emission within a 435 beam. To determine the coefficients of the relation log(L_FIR) = a log(${L}_{{\rm{CO}}}^{\prime }$) + b, we used the python module lnr.bces¹¹ , which utilizes the Bivariate Correlated Errors and intrinsic Scatter method of Akritas & Bershady (1996). This is important because we have errors on both variables (we introduced a non-negligible error into the L_FIR variable through our source/beam correction). As stated above, we chose the examination of L_FIR as a function of ${L}_{{\rm{CO}}}^{\prime }$ to match the most recent literature and note that the solution to the inverse problem (${L}_{{\rm{CO}}}^{\prime }$ as a function of L_FIR) does not simply produce best-fit slopes that are the inverse of those presented here. For the case of low-J lines collected from the literature and our ARO follow-up, multiple measurements for the same galaxies (positioned at the location matching the FTS coordinates) are treated independently.

When including all spectra in our sample, we find slopes starting at 1.3 for CO J = 1−0, and lowering to about ∼1 for the mid- to high-J CO lines, with no discernible trend with increasing J. The results are shown in Table 4. There is no significant difference whether we include or exclude FTS lines with S/N < 3. We also separated our samples into a few categories that are somewhat overlapping, and summarize the differences here:

• 1.  
  We separated our galaxies into known AGN (categorized in Hyperleda¹² as a quasar or any type of Seyfert galaxy) or not. This classification (using the agnclass category) is taken from the Véron-Cetty & Véron (2006) catalog. Within those classified by Hyperleda as AGNs (78 of the 232 galaxies), no information is provided regarding the relative SF versus AGN contributions to the total L_FIR, so this is a somewhat crude division of the galaxy sample. Looking only at AGNs (Table 5), compared to the whole sample we find higher slopes for the low-J lines (1.5 ± 0.1, 1.2 ± 0.1, 1.3 ± 0.1, 1.1 ± 0.1, 1.15 ± 0.04 for J = 1−0 through J = 5−4), but the slope error bars overlap by CO J = 6−5 and continue to follow the same trend as the whole sample. Looking at the sample that completely excludes AGNs, we find a lower slope than the whole sample for CO J = 1−0 (1.13 ± 0.06), but at subsequently higher-J the slopes are not distinguishable from the combined sample. In summary, comparing the AGN to the non-AGN sample, the most significant difference is in the CO J = 1−0 line. There are also differences, at less significance, up to the J = 5−4 line. The low-J CO emission is not expected to be affected by the AGN; the L_FIR and high-J CO are likely being influenced. We likely do not see any difference because our AGN-designated galaxies (which are only about one third of the full sample) may not be entirely dominated in their molecular excitation from the AGN; as mentioned, we do not separate the relative SF versus AGN contributions to total L_FIR. Better quantification of the AGN influence, and higher spatial resolution, may result in differences in the slope.
• 2.  
  Astronomers often separate galaxies into (U)LIRGs or lower luminosity galaxies. We separate our galaxies at L_FIR = 6 × 10¹⁰L_⊙ which we found is approximately equal to L_IR(8–1000 μm) = 1 × 10¹¹, based on the luminosities listed in Greve et al. (2014) and K14. (The exact cutoff value does not change the conclusions that follow.) The CO J = 1−0 line has been known to not be fit by one slope among (U)LIRGs and lower luminosity galaxies; the superlinear slope that results from a single power line fit is due to higher dense gas fractions in (U)LIRGs (Gao & Solomon 2004b; Greve et al. 2014), essentially creating a higher intercept (but the same slope). We do find lower slopes in J = 1−0 among these two populations fit separately than their combined slope of 1.3 ± 0.4. We find that the mid-J lines of CO are well fit by a single slope across many orders of magnitude; the slopes in all three cases (full sample; just (U)LIRGs; just galaxies with lower luminosities than (U)LIRGs) are not statistically distinguishable given the error bars. Our three highest-J CO lines (J = 11−10, J = 12−11, J = 13−12), however, do show a measurable difference in the best-fit slope. Focusing primarily on (U)LIRGs changes the slopes at higher-J, decreasing to about 0.83 ± 0.03 (weighted average of the three highest-J CO lines in Table 6), a highly significant difference from 1. The results for these two populations are shown in Tables 6 and 7.
• 3.  
  We restricted the fit to only those galaxies that are not particularly well resolved, those with L_FIR,beam ≥ 0.8 L_FIR,total based on the photometry maps, to see if our infrared luminosity correction could be influencing the y-axis values. We find lower slopes in this case, 0.88 ± 0.06 on average for high-J lines, but this may be due to the fact that 2/3 of this population are (U)LIRGs, not because large L_FIR corrections are necessarily inaccurate.

Certainly, these slopes are masking a considerable amount of scatter in the data, and different trends can be discerned when fitting different populations (also seen in Liu et al. 2015). The remaining figures of this paper examine the survey populations themselves with comparisons to these derived, "global" relations.

Fitting a linear relation in log space requires creating error bars that are symmetric in log space, which ours are not. The exaggeration of the linear errors when viewed in log space is highest for those line measurements with the lowest S/N. To test whether the conversion to these symmetric error bars (which decreases the size of the lower error bar) introduces systematic bias in deriving the slopes, we chose fixed values of the slopes, set the nominal y-values to match the chosen slope, and draw randomly from each data point's x and y error bars in linear space. Fitting such scattered distributions many times results in distributions centered upon the chosen slopes and with scatter comparable to our error bars in Table 4. We find no evidence of systematic bias due to this effect.

We also investigated whether our error bars are responsible for the slightly lower value of the CO J = 1−0 line's slope at 1.3 instead of 1.4–1.6. This value is not due to the error bars in either the x- or y-directions; we still find a slope of about 1.3 if we exclude one, the other, or both in the line fitting. We do find a slope of 1.44 if we restrict the fit to those data points with log (${L}_{{\rm{CO}}}^{\prime }$) > 9, indicating that the lower luminosity points are bringing down the slope overall.

Greve et al. (2014) conducted a comparison of CO emission to L_FIR for the non-extended galaxies of the HerCULES sample (van der Werf et al. 2010) and found sharply decreasing slopes between log(L_FIR) and log(${L}_{{\rm{CO}}}^{\prime }$), starting around 1 at CO J = 1−0 and decreasing as one moves to higher-J starting around CO J = 6−5 (0.93 ± 0.05, down to 0.47 ±0.20 by the J = 13−12 line).

They interpret this sublinear slope as an indication that far-UV (FUV) radiation fields are not responsible for the the dense, warm gas emitting in the high-J lines. Mechanical heating and shocks are the likely explanation for the high-J excitation, as was shown in many galaxies of the sample of K14 (and references therein). We are in agreement with Greve et al. (2014) that mechanical feedback can be related to SF, such as stellar outflows, winds, and supernova remnant expansion; not being related to the FUV radiation from stars does not mean the excitation is not related to SF at all.

The fluxes used in Greve et al. (2014) are reported in Rosenberg et al. (2015). Their sample is included in ours, but our measured fluxes do not match in all cases. For lines up to J = 10−9, 2/3 of their 18 non-extended galaxies match our flux values within 30%. At higher-J, this number is more like 1/3. When discrepant, their values are often ∼50% higher, or even factors of a few for IRAS F05189-2524, for which we do not use the same spectra. The major difference for other galaxies is likely the treatment of the source/beam correction. Importantly, Greve et al. (2014) only included (U)LIRGs (log(L_FIR) ≥ 11). As discussed above, we find lower slopes if we restrict ourselves to this population, but not as low as the slopes reported for Figure 1 by Greve et al. (2014), due to the differences in our measured fluxes.

Liu et al. (2015) conducted a larger study, more comparable to ours, utilizing the full extent of the Herschel archive. With the inclusion of galaxies from log(L_FIR) ≥ 8, they found linear slopes between log(L_FIR) and log(${L}_{{\rm{CO}}}^{\prime }$) throughout the CO ladder starting at J = 4−3 (statistically consistent with our results). Also consistent with our results and Greve et al. (2014), restricting the fits to the HerCULES sample of (U)LIRGs yields sublinear slopes for the highest-J lines. Their work does not include a table of the galaxies included in their sample or the line fluxes used, so we cannot make a direct comparison of the fluxes used for our relations at this time. Some major differences may be important in the comparison of our work: first, they use SPIRE calibration version 12 and we use version 13, which may make the most difference for galaxies observed near the start of SPIRE cooling cycles (due to "cooler burp"). Second, each of their CO and L_FIR relations use a different beam size, that of the frequency of the CO line in the FTS. Third, for sparsely and intermediately sampled galaxies, they extract spectra from multiple bolometers in the detector array, meaning they are including multiple spatially Nyquist-sampled data points in their relationships (with their own matching L_FIR values) for such galaxies. Finally, their line fitting procedure uses sinc-convolved-Gaussian (SCG) lines in HIPE for sparse observations and sinc for intermediate or full sampling observations. They note that SCG-derived fluxes are systematically higher, but we found that few galaxies have CO linewidths large enough to be resolved by the FTS. Differences in the slopes of (U)LIRGs among Liu et al. (2015), Greve et al. (2014), and this work can also be attributed to the small dynamical range spanned by (U)LIRGs. Despite these differences, our results agree well in finding mid- to high-J CO slopes of 1 in a large sample of galaxies and less than 1 for (U)LIRGs.

Zoom In Zoom Out Reset image size

Lu et al. (2014) also examined the CO J = 4−3 through J = 13−12 emission of the 65 LIRGs in the Great Observatories All-Sky LIRG Survey, comparing to the IRAS 60–100 μm color, C(60/100), as a proxy for dust temperature. They demonstrated that L_IR is not the best predictor of SLED shapes; we find an overall trend in Figure 2 (the line luminosity ratios relative to J = 1−0 in average SLEDs by L_FIR bin increase with L_FIR), but also a significant amount of variation in Figure 3 (individual mid- to high-J CO/J = 1−0 luminosity ratios for each galaxy, unbinned), discussed in the next section. C(60/100), which is not presented here, is potentially a better predictor of the CO SLED shape (based on the line at which the SLED peaks in luminosity).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Within galaxy-wide log(L_FIR) ranges of approximately 0.5 dex, we compiled weighted-averaged (L_CO/L_FIR)/(L_CO(J = 1−0)/L_FIR) ratios, presented in Figure 2 (top; see caption for note about conversion to ${L}_{{\rm{CO}}}^{\prime }$). The SLEDs are normalized to the J = 1−0 line to show the relative excitation across the CO ladder. The CO J = 1−0 line measures emission from the same type of cold, ubiquitous molecular gas found throughout many types of galaxies. There are four main trends to examine in this plot: first, there is a trend with increasing L_FIR toward much more high-J CO luminosity compared to J = 1−0. This indicates a greater energy input relative to typical PDRs, to explain the high-J emission. Second, the location of the CO luminosity peak moves to higher-J with higher L_FIR. Third, the slope of the mid- to high-J CO emission SLED becomes flatter with increasing L_FIR. However, all the SLED slopes are relatively flat, none show an extreme drop-off. Fourth, the values of mid- to high-J CO relative to CO J = 1−0 in luminosity only range from about 10 to 100 across all lines.¹³ This is consistent with the CO SLED compilation shown in Figure 5 of K14, which showed high-J/J = 1−0 ratios from about 5 to 100 (with one outlier, NGC 6240, at 250 at its peak, see Meijerink et al. 2013). These last three trends are indicative of a higher average pressure (product of kinetic temperature and density) required to explain the shape of the CO SLED. We emphasize average because all of the CO emission is the sum of a gradient of conditions from different environments.

Figure 2 (top) also shows a difference between the bin-averaged SLEDs (with data points) and the SLEDs that would be predicted simply from our log(L_FIR) versus log(${L}_{{\rm{CO}}}^{\prime }$) slopes (light colored, no data points). The predicted SLEDs span a narrower range than the bin-averaged SLEDs. Describing each CO line by a single relationship over the entire sample averages over the very real differences in populations, such as those described in Section 3.1 (e.g., different L_CO/L_FIR slopes at higher luminosities). Also, the actual weighted averages are often influenced by high S/N galaxies that may not represent the whole bin, but still illustrate the overall trend. Either way of examining the data, which we do more in the following section, shows a relatively narrow range of high-J to CO J = 1−0 ratios compared to expectations from some theoretical models.