© ESO, 2017

1. Introduction

The cosmic star formation rate density (ρ_SFR) has declined by a factor of 20 since an observed peak at z ~ 2.5 (Hopkins & Beacom 2006; Madau & Dickinson 2014), and it remains unknown whether this is due to the exhaustion of the galactic interstellar medium (ISM), a reduction in the accretion of the pristine intergalactic medium, or a decline in the efficiency in the conversion of gas to stars. One approach towards disentangling the physical processes contributing to the decline in the overall ρ_SFR requires the characterisation of the content and physical conditions of interstellar gas in galaxies at all redshifts.

Probing the physical conditions of the ISM requires observations of emission lines, such as the far-infrared (FIR) fine-structure lines of [Cii] 158 μm, [Nii] 122, 205 μm, [Oi] 63, 145 μm, and [Oiii] 88 μm lines, which play a crucial role in the thermal balance of the gas, and the rotational transitions of carbon monoxide (CO). In particular, the [Cii] 158 μm line (ν_rest = 1900.54 GHz), which originates from the ²P_3/2 → ²P_1/2 transition of the ground state of singly ionised carbon, typically has a luminosity of 0.1–1% of the far-infrared luminosity in normal star-forming galaxies, thus making it one of the dominant cooling lines (e.g. Dalgarno & McCray 1972; Crawford et al. 1985; Stacey et al. 1991a). The [Cii] line emission comes from both neutral and ionised gas, as the low ionization potential of atomic carbon means C⁺ can be produced from far-ultraviolet (FUV) photons with energies greater than just 11.26 eV (cf. hydrogen’s ionisation potential of 13.6 eV). In star-forming galaxies, the majority of the line emission (~70%; e.g. Stacey et al. 1991a, 2010a) is shown to arise from photodissociation regions (PDRs), with the remaining fraction coming from star-forming Hii regions, lower density warm gas and diffuse Hi clouds, X-ray dominated regions (XDRs), and cosmic ray dominated regions (CRDRs). Deeper within the PDR regions, C⁺ becomes converted to CO, which is a standard tracer of the molecular gas content. The [Cii]/CO(1–0) ratio is mostly dependent on the C⁺ column density and surface temperature, which both decrease with stronger gas shielding (Kaufman et al. 1999). These lines are thus useful for constraining the ISM conditions in galaxies (see e.g. Tielens & Hollenbach 1985; Wolfire et al. 1990).

Recent advancements in space- and ground-based facilities for observing these emission lines at FIR and submillimetre wavelengths are rapidly expanding our knowledge of the ISM in both nearby and high redshift galaxies (see e.g. Solomon & Vanden Bout 2005; Carilli & Walter 2013). The emission of [Cii] and other fine-structure lines has been observed in low-z galaxy samples with airborne- or space-based observatories, such as the Kuiper Airborne Observatory (KAO; e.g. Stacey et al. 1991a; Madden et al. 1993) and the Infrared Space Observatory (ISO; e.g. Hunter et al. 2001; Malhotra et al. 2001; Brauher et al. 2008). The Herschel Space Observatory (Pilbratt et al. 2010) with the PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) instruments was capable of observing both the FIR cooling lines and FIR/submm spectral energy distribution at unprecedented resolution, enabling the study of gas heating and cooling (via the [Cii]/L_TIR or ([Cii]+[Oi]63)/L_TIR ratios) on galactic and spatially resolved, sub-kiloparsec scales (see e.g. Croxall et al. 2012; Lebouteiller et al. 2012; Parkin et al. 2013; Hughes et al. 2015).

At higher redshifts, studies of the ISM primarily rely on observations of the [Cii] line and the rotational transitions of CO as diagnostics of the physical conditions. However, sources at z > 1 avoid strong atmospheric absorption at submillimetre wavelengths and can be observed with ground-based instrumentation, such as the Northern Extended Millimetre Array (NOEMA), the Atacama Pathfinder EXperiment (APEX), and the Atacama Large Millimeter/submillimeter Array (ALMA). Numerous studies over the past decade (see e.g. Gullberg et al. 2015, and references therein) report the detection of [Cii] emission in high-z sources that are classified as possible active galactic nuclei (AGN) hosts (e.g. Maiolino et al. 2005; Stacey et al. 2010a; Wang et al. 2013) or starburst galaxies (see e.g. Ivison et al. 2010; Hailey-Dunsheath et al. 2010; Stacey et al. 2010a; Cox et al. 2011; Swinbank et al. 2012; Riechers et al. 2013; Magdis et al. 2014; Brisbin et al. 2015). A significant fraction of these detections also has CO(1–0) observations, such as the sample of gravitationally lensed, dusty star-forming galaxies in the redshift range z ~2.1–5.7 (Gullberg et al. 2015; Aravena et al. 2016) discovered by the South Pole Telescope (SPT; Carlstrom et al. 2011). The increasing number of observations available for systems over a wide redshift range means we can begin to characterise the ISM physical conditions at various epochs.

The physical properties of the gaseous components of the ISM may be determined by comparing the observed ratio of the [Cii] 158 μm and CO(1–0) line emission to the predictions of a PDR model. There are numerous PDR models available for determining the gas density, temperature and strength of the FUV radiation field (see Röllig et al. 2007, for a discussion, and references within). One of the most commonly used PDR models is that of Tielens & Hollenbach (1985), which characterises the physical conditions in a semi-infinite, plane-parallel slab PDR by two free variables: the hydrogen nucleus density, n, and the strength of the FUV (6 < hν < 13.6 eV) radiation field, G₀, which is in units of the Habing Field (Habing 1968). The model has since been updated by Wolfire et al. (1990), Hollenbach et al. (1991) and Kaufman et al. (1999, 2006). Predictions from PDR models have been compared to Herschel observations of both Galactic PDRs and nearby galaxies. For example, Croxall et al. (2012) studied a late-type spiral, NGC 1097, and a Seyfert 1 galaxy, NGC 4559, finding 50 ≤ G₀ ≤ 1000 varying with 10^2.5 cm^-3 ≤ n ≤ 10³ cm^-3 across both discs. Most recently, Parkin et al. (2013) examined the n and G₀ in various regions of M 51; the hydrogen density and FUV radiation peak in the nucleus and similarly decline in both the spiral arm and interarm regions, suggesting similar physical conditions in clouds in these environments (see also Parkin et al. 2014; Hughes et al. 2015). Stacey et al. (2010a) posit that the observed L_[Cii]/L_FIR and L_[Cii]/L_CO luminosity ratios suggest that the gas density and FUV radiation field in their sample galaxies at z ~ 1−2 are similar to the physical conditions in local starburst systems, which is also supported by the observations of the SPT sample (Gullberg et al. 2015).

In this paper, we study the physical conditions in the interstellar gas for a sample of 27 dusty galaxies at 0.03 < z < 0.2, selected from the Valparaíso ALMA Line Emission Survey (VALES; Villanueva et al., in prep.). These are characterised as normal star-forming galaxies (see Fig. 1 of Villanueva et al., in prep.) that generally lie on or slightly above the main sequence (see e.g. Elbaz et al. 2011). We use new ALMA Band 3 CO(1–0) observations combined with Herschel-PACS [Cii] 158 μm line emission data and FIR luminosities determined with photometry from Herschel-SPIRE and other facilities (Ibar et al. 2015), to investigate the physical properties of the interstellar gas in these galaxies by using the PDR model of Kaufman et al. (1999, 2006). We compare the physical conditions in our VALES sample to similar studies at low- and high-z. Our paper is structured as follows: in , we describe our sample, observations and data reduction methodology. In and , we describe the characteristics of the gas and compare our observations to theoretical PDR models. Finally, and present our discussion and conclusions. Throughout this paper, we adopt a ΛCDM cosmology with H₀ = 70 km s^-1 Mpc^-1, Ω_M = 0.27 and Ω_Λ = 0.73.

2. The sample and data

Our sample of galaxies is drawn from a Herschel programme, the Herschel Astrophysical Terahertz Large Area Survey (Eales et al. 2010; Valiante et al. 2016; Bourne et al. 2016), which is capable of providing a sufficient number of far-IR bright galaxies over ~600 deg². In addition to a wealth of high-quality ancillary data, a significant sample have both Herschel-PACS [Cii] 158 μm (Ibar et al. 2015) and our follow-up ALMA CO(1–0) emission line observations from VALES (Villanueva et al., in prep.). From the three equatorial fields (totaling ~160 deg²) covered by H-ATLAS, galaxies were selected based on the following criteria: (1) a flux limit of S_{160 μm} > 150 mJy, i.e. near the peak of the SED of a typical, local star-forming galaxy; (2) no neighbours with S_{160 μm} > 160 mJy (3σ) within 2 arcmin from their centroids; (3) an unambiguous identification (reliability > 0.8, Bourne et al. 2016) in the Sloan Digital Sky Survey (SDSS DR7; Abazajian et al. 2009); (4) a Petrosian SDSS r-band radius <15″, i.e. smaller than the PACS spectroscopic field of view; (5) high-quality spectroscopic redshifts (zqual > 3) from the Galaxy and Mass Assembly survey (GAMA; Driver et al. 2009, 2011; Liske et al. 2015); and (6) a redshift between 0.02 < z < 0.2 (median of 0.05), beyond which the [Cii] emission becomes redshifted to the edge of the PACS spectrometer 160 μm band. After applying these criteria, 324 galaxies remain to comprise a statistically significant sample spanning a wide range of optical morphological types and IR luminosities. Of these, 27 objects have observations of both the [Cii] and CO(1–0) emission lines that form the focus of our study.

2.1. Herschel-PACS [CII] 158 micron line observations

The [Cii] 158 μm line observations of our sample are presented in Ibar et al. (2015). The end of the Herschel mission meant only 28 galaxies, which are all located in the GAMA 9h field, of the parent sample could be observed during our Herschel-PACS [Cii] spectroscopic campaign. Their selection arises purely on the basis of scheduling efficiency, so this sample is representative of the original sample albeit smaller in number (see Fig. 1 of Ibar et al. 2015). We first briefly summarise their properties.

The redshifted [Cii] 158 μm line emission was observed with the PACS first order (r1) filter covering ~2μm of bandwidth in a 47 arcsec × 47 arcsec field of view in a single pointing (pointed mode). The central spaxel of the 5 × 5 spaxal array captured the majority of the line emission for each target. Data cubes comprising 5 × 5 spaxels × n_chan rebinned spectral channels were generated from calibrated PACS level-2 data products processed with SPG v12.1.0 with an effective spectral resolution of ~190–240 km s^-1. The [Cii] 158 μm line flux density was determined from the weighted sum (aided by the instrumental noise cube) of the central 3 × 3 spaxels via the simultaneous fitting of a linear background slope and a Gaussian to the spectra. Uncertainties at the 1σ level were derived for the line parameters with a Monte Carlo realisation (1000 ×), randomly varying the signal per spectral channel using the instrumental error cube. These [Cii] line flux measurements are presented in ; Ibar et al. (2015) provide more details.

2.2. ALMA CO(1–0) line observations

Despite the fact that ALMA is not an ideal telescope for creating galaxy surveys, we use a novel observational approach to target large numbers of galaxies at low and intermediate redshifts to create a reference sample for interpreting observations of the high redshift Universe. Since our sample is drawn from the GAMA fields that are each ~4° × 14° in size and therefore provide large numbers of galaxies at similar redshifts, we can minimise the number of spectral tunings needed to observe all sources independently by setting the source redshifts to zero and fixing the spectral windows (SPW) manually to cover the widest possible spectral (thus redshift) range. The central frequency position of each SPW is then optimised to cover the redshifted CO(1–0) line for the maximum number of sources, therefore significantly minimising overheads. Using this observing strategy, we obtained observations targeting the CO(1–0) line for 67 galaxies during cycle 1 and 2 (project 2013.1.00530.S; P.I.: E. Ibar). Whereas Villanueva (in prep.) present the observations, data reduction, and analysis in detail for the complete sample, in this paper we focus solely on galaxies with the cycle 2 CO(1–0) and Herschel [Cii] data. Here, we briefly summarise these observations and data reduction steps.

The cycle-2 observations taken in Band 3 covered the ¹²CO(1–0) emission line down to 2 mJy beam^-1 at 30 km s^-1 for a sample of 27 sources at 0.03 <z< 0.2. All observations were reduced homogeneously within the Common Astronomy Software Applications¹ (CASA; McMullin et al. 2007) using a common pipeline, developed from standard pipelines, for calibration, concatenation and imaging. The bandpass calibrators for our sources were J0750+1231, J0739+0137 and J0909+0121, flux calibrators were Callisto, Mars, Ceres and Titan, and the phase calibrators were J0909+0121 and J0901-0037. The optimal image resolution was that which provides the highest number of non-cleaned point-like >5σ sources when varying the cube spectral resolution from 20 to 100 km s^-1 in steps of 10 km s^-1 with the clean task; cubes without detections were set at 100 km s^-1 channel width. The final cubes, of 256 × 256 pixels in spatial size, were corrected for the primary beam, manually cleaned to a threshold of 3σ, and created using natural weighting. The barycentre velocity reference was set as the optical spectroscopic redshift (z_spec) of each source.

For measuring the CO(1–0) line emission, we first perform a Gaussian fit to the spectra to identify the central frequency, ν_obs, and ν_FWHM of the line. The velocity integrated CO(1–0) flux densities (S_COΔv in units of Jy km s^-1) were then obtained by collapsing the data cubes between ν_obs−ν_FWHM and ν_obs + ν_FWHM, and fitting these cubes with a Gaussian. We detect >95% (26 of 27) of the targets with a >5σ peak line detection. We note that our selection criteria includes targets with SDSS r-band radii smaller than 15″ to obtain reliable Herschel-PACS observations (Ibar et al. 2015). Considering that the maximum recoverable angular scales of ALMA in Band 3 are approximately 25″ (with a 60″ primary beam), these detections do not suffer from any missing flux. For the upper limit of source G09.DR1.328, we collapse the 100 km s^-1 spectral resolution whilst setting ν_FWHM = 250 km s^-1 and adopt the limit as 5 × the measured RMS. Our CO(1–0) line flux measurements are presented in ; Villanueva (in prep.) provide more details. We calculate the CO line luminosity, L_CO, in units of L_⊙ following ${\mathit{L}}_{\mathrm{CO}}\mathrm{=}\mathrm{1.04}\mathrm{\times }{\mathrm{10}}^{-3}{\mathit{S}}_{\mathrm{CO}}\mathrm{\Delta }\mathit{v}\hspace{0.17em}{\mathit{\nu }}_{\mathrm{rest}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-1}{\mathit{D}}_{\mathit{L}}^{\mathrm{2}}\mathit{,}$(1)where ν_rest is the rest frequency in GHz and D_L is the luminosity distance in Mpc (from Eq. (1) of Solomon & Vanden Bout 2005).

2.3. Ancillary data for far-infrared luminosities

We use previously published values of the FIR luminosity, L_FIR, which were measured for each galaxy by fitting the SED constructed from Herschel PACS and SPIRE, WISE-22 μm, and IRAS photometry (see Ibar et al. 2013, 2015). In brief, each rest-frame SED is fit with a modified black body that is forced to follow a power law at the high-frequency end of the spectrum and the flux of the best-fitting SED is integrated between 8 μm and 1000 μm, i.e. ${\mathit{L}}_{\mathrm{FIR}}\mathrm{(}\mathrm{8}\mathrm{-}\mathrm{1000}\hspace{0.17em}\mathit{\mu }\mathrm{m}\mathrm{)}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{D}}_{\mathit{L}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\hspace{0.17em}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathrm{1}}}^{{\mathit{\nu }}_{\mathrm{2}}}{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }$(2)to obtain the dust temperature (T_d), dust emissivity index (β), mid-IR slope (α_{mid - IR}), and normalisation, which provides the total FIR luminosity (rest-frame 8–1000 μm). These L_FIR estimates and their uncertainties, which were obtained from randomly varying the broadband photometry within their uncertainties in a Monte Carlo simulation (100 ×), are listed in , whereas the complete set of derived parameters may be found in Table 3 of Ibar et al. (2015).

Previous studies use differing definitions of the total infrared band (TIR, e.g. 3–1100 μm; Parkin et al. 2013, 2014) and also the FIR bands that are commonly used as proxies of the total infrared emission, such as from 42 μm to 122 μm (e.g. Helou et al. 1988; Dale et al. 2001; Stacey et al. 2010a), or 40 μm to 500 μm (e.g. Graciá-Carpio et al. 2008; Gullberg et al. 2015). These quantities are related by F_TIR(3–1100 μm) ≈2 × F_FIR (42–122 μm) in Dale et al. (2001) and F_TIR (8–1000 μm) ≈1.3 × F_FIR(40–500 μm) in Graciá-Carpio et al. (2008). The difference in the IR luminosities obtained over 3–1100 μm and 8–1000 μm is typically less than 15% (e.g. Rosario et al. 2016), thus comprising only a small fraction of the TIR emission. Therefore, the FIR luminosity, as defined in this present work (8–1000 μm), is a good proxy for the TIR emission and is approximately equivalent to the bolometric FIR flux of the PDR model (see e.g. Parkin et al. 2013, 2014), which we discuss shortly.

3. The [CII] – CO luminosity correlation

Before deriving the physical conditions via PDR modelling, we first perform a brief sanity check on the observations of our VALES sample by examining the well-studied [Cii]–CO luminosity correlation (see e.g. Crawford et al. 1985). Previous studies have shown that a variety of sources, from normal star-forming galaxies to starbursts and AGN hosts, follow the correlation up to z ~ 6 with an observed L_[Cii]/L_CO ratio ranging from 1300 to 6300 with a typical median of 4400 (Crawford et al. 1985; Stacey et al. 1991a, 2010a; Swinbank et al. 2012). In , we plot the [Cii] luminosity versus the CO(1–0) luminosity and superimpose our H-ATLAS-based sample at 0.03 <z< 0.2 onto the observations of samples at higher and lower redshifts presented by Gullberg et al. (2015, see their Fig. 7). We find the average L_[Cii]/L_CO ratio for our sample is 3500 (median of ~2400) with a standard deviation of 1200; this value is slightly higher than the average value (1300 ± 440) for the low-z, normal star-forming galaxies, but much lower than that of the high-z SPT sample (5200 ± 1800; Gullberg et al. 2015).

Following the reasoning and methodology of Gullberg et al. (2015), we use the L_[Cii]/L_CO ratio to constrain the optical depths and excitation temperatures of the lines from a comparison of the source functions. The L_[Cii]/L_CO ratio is then given by $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\frac{{\mathit{L}}_{\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}}{{\mathit{L}}_{\mathrm{CO}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}}}\mathrm{=}{\left(\frac{{\mathit{\nu }}_{\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}}{{\mathit{\nu }}_{\mathrm{CO}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}}}\right)}^{\mathrm{3}}\mathrm{\times }\left(\frac{\mathrm{\Delta }{\mathit{\nu }}_{\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}}{\mathrm{\Delta }{\mathit{\nu }}_{\mathrm{CO}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}}}\right)\mathrm{\times }\frac{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}}}{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{CO}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}}}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{\times }\frac{{\mathrm{e}}^{\mathit{h}{\mathit{\nu }}_{\mathrm{CO}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}}\mathit{/}\mathit{k}{\mathit{T}}_{\mathrm{ex}\mathit{,}\mathrm{CO}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}}}\mathrm{-}\mathrm{1}}{{\mathrm{e}}^{\mathit{h}{\mathit{\nu }}_{\mathrm{\left[}\mathit{C}{ii}\mathrm{\right]}}\mathit{/}\mathit{k}{\mathit{T}}_{\mathrm{ex}\mathit{,}\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}}\mathrm{-}\mathrm{1}}\hspace{0.17em},\end{array}& & \end{array}$(3)under the assumption that the [Cii] and CO(1–0) filling factors are equal. Knowing the L_[Cii]/L_CO ratio leaves the two excitation temperatures (T_{ex, [Cii]}, T_ex,CO(1−0)) and two opacities (τ_[Cii], τ_CO(1−0)) as free parameters. We can then vary these free parameters to match the observations.

We first consider the scenario in which the excitation temperatures are equal, T_{ex, [Cii]} = T_ex,CO(1−0). In this case, underpredicts the observed L_[Cii]/L_CO ratio by an order of magnitude when the [Cii] is optically thin (τ_[Cii] = 0.1) and CO is optically thick (τ_CO(1−0) = 1), but overpredicts the ratio by an order of magnitude when [Cii] is optically thick (τ_[Cii] = 1) and CO is optically thin (τ_CO(1−0) = 0.1). The observed ratio is only reproducible when the optical depths are the same (τ_[Cii] = τ_CO(1−0)) but requires equal excitation temperatures of >50 K and, given that the CO is usually optically thick (e.g. from the ¹²CO/¹³CO line ratios; see Hughes et al., in prep.), this implies that the [Cii] line would also be approaching the optically thick regime. However, such a scenario, in which both [Cii] and CO lines are optically thick with equal excitation temperatures over 50 K, is not a viable solution given observational evidence to the contrary. Firstly, [Cii] optical depths are consistently of the order unity or less for even very bright Galactic star formation regions (Stacey et al. 1991b; Boreiko & Betz 1996; Graf et al. 2012; Ossenkopf et al. 2013). Secondly, although the observed [Cii]/[Oi]line ratios (see e.g. Stacey et al. 1983a; Lord et al. 1996; Brauher et al. 2008) and measurements using velocity-resolved peak [Cii] antenna temperatures (Graf et al. 2012; Ossenkopf et al. 2013) yield [Cii] excitation temperatures that are far in excess of 50 K, CO excitation temperatures on galactic scales are typically less than 50 K (e.g. Gullberg et al. 2015).

Fig. 1 [Cii] luminosity vs. the CO(1–0) luminosity for the VALES sample at 0.03 <z< 0.2 superimposed on the low- and high-z samples taken from Fig. 7 of Gullberg et al. (2015), with symbols as specified in the legend. Downward triangles represent 5σ upper limits. The grey shaded region the 1σ spread in the mean observed L[Cii]/LCO ratio (3500 ± 1200; blue dashed line) of our sample, which we compare to the mean ratios of the low-z sample (1300 ± 440; black dotted line).

The alternative scenario to consider is when the two excitation temperatures are different. We examine the relationship between T_{ex, [Cii]} and T_ex,CO(1−0) for our observed L_[Cii]/L_CO ratio considering various optical depths of the two lines (see ), finding that T_{ex, [Cii]}>T_ex,CO(1−0) in all cases and that we cannot exclude that both [Cii] and CO could be optically thick (τ_[Cii] = 1, τ_CO(1−0)> 1). In addition, we can fix the CO excitation temperature, such as in Weiß et al. (2013) and Gullberg et al. (2015), by making the assumption that the CO traces molecular gas within which dust is thermalised, such that the dust temperatures obtained for our VALES sample (20 K <T_d< 55 K; Ibar et al. 2015) are representative of T_ex,CO(1−0). From , the fixed range of T_ex,CO(1−0) implies T_{ex, [Cii]} ≈ 35–90 K for τ_[Cii] = τ_CO(1−0) and T_{ex, [Cii]} ≈ 45–110 K for τ_[Cii] = 1 with τ_CO(1−0) = 4. However, we also predict that our L_[Cii]/L_CO ratio can arise from gas with T_{ex, [Cii]} ≈ 100–300 K with optically thin [Cii] (τ_[Cii] = 0.1) and optically thick CO (τ_CO(1−0) = 1). When the gas density is greater than the [Cii] critical density, the excitation temperature becomes equivalent to the gas temperature and, as we present in the following section (), the PDR models support this latter scenario.

To summarise, we conclude that the best scenario to explain the L_[Cii]/L_CO ratios observed in our VALES sample is where the [Cii] and CO emission originates from gas with T_{ex, [Cii]}>T_ex,CO(1−0), where [Cii] is optically thin and CO is optically thick.

Fig. 2 CO(1–0) vs. [Cii] excitation temperatures determined via the ratio of the source functions ()) assuming various line opacities, as specified in the legend, and adopting the observed L[Cii]/LCO ratio of 3500. The light and dark grey shaded regions denote denotes the 1σ spread in the mean observed ratio, whereas the red shaded region represents the CO excitation temperatures of 20 K <Tex,CO(1−0)< 55 K. Combined with the results of PDR modelling and previous observational evidence in the literature (), the best scenario to explain the observed L[Cii]/LCO ratios is gas where Tex, [Cii]>Tex,CO(1−0) and with [Cii] is optically thin and the CO is optically thick (τ[Cii] = 0.1, τCO(1−0) = 1; black dotted line).

Fig. 3 Diagnostic diagram of the observed L[Cii]/LFIR ratio vs. LCO/LFIR ratio for our sample. The observations are superimposed onto the grid of constant hydrogen nuclei density, log n (black solid lines), and FUV radiation field strength, log G0 (black dotted lines), determined from the PDR model of Kaufman et al. (1999, 2006). The unadjusted observations (blue solid circles) are compared to the adjusted observations, as described in , where the fraction of the [Cii] emission arising from neutral gas is fixed at 70% (red semi-open circles) and 50% (red open circles). Downward triangles represent 5σ upper limits. Two sources lie above the log n = 2 contour, i.e. outside the model parameter space, and the blue arrow points towards G09.DR1.328, which remains outside the plot regardless of adjustments. The grey shaded region denotes the 1σ spread in the mean observed L[Cii]/LCO ratio (blue dashed line). This figure is adapted from Stacey et al. (2010a) and Hailey-Dunsheath et al. (2010).

4. Results from PDR modelling

We now compare our observations to the PDR model of Kaufman et al. (1999, 2006), which is an updated version of the Tielens & Hollenbach (1985) model. The model treats PDR regions as homogeneous infinite plane slabs of hydrogen with physical conditions characterised by the hydrogen nuclei density, n, and the strength of the incident FUV radiation field, G₀, which is normalised to the Habing Field in units of 1.6 × 10^-3 erg cm^-2 s^-1 (Habing 1968). The gas is collisionally heated via the ejection of photoelectrons from dust grains and PAH molecules by FUV photons, and gas cooling from line emission is predicted by simultaneously solving the chemical and energy equilibrium in the slab. For a given a set of observations of spectral line intensities, the corresponding G₀ and n values predicted by the PDR model are available online² via the “Photo Dissociation Region Toolbox” (PDRT, Pound & Wolfire 2008), where the models cover a density range of 10¹ ≤ n ≤ 10⁷ cm^-3 and a FUV radiation field strength range of 10^-0.5 ≤ G₀ ≤ 10^6.5. In the following, the [Cii] and CO(1–0) observations are compared to the PDR model grid lines from the Kaufman et al. (1999) diagnostic plots, for which we must assume that each emission component – the [Cii] line emission, CO(1–0) emission, and the FIR continuum – originates from a single PDR component in our sources.

In , we superimpose the observed (i.e. unadjusted) L_[Cii] versus the L_CO line luminosities for our sample on the PDR model grid lines of constant log (n/ cm^-3) and log G₀ constructed from the Kaufman et al. (1999, 2006) diagnostic plots, adapted from figures in Stacey et al. (2010a) and Hailey-Dunsheath et al. (2010). The majority of our observations lie within the parameter space covered by the PDR model and exhibit moderate FUV radiation field strengths (2.0 < log G₀ < 3.0) and moderate hydrogen densities (2 < log n/ cm^-3 < 4.5). We stress, however, that there is much uncertainty in the lower limit of the latter parameter, owing to the degeneracy in the parameter space between densities of log n/ cm³ = 2 and 3. In addition, three of our galaxies fall outside of the L_[Cii]/L_FIR versus the L_CO/L_FIR parameter space defined by the PDR model; the closest contours of constant log n and log G₀ are those of the lowest density and weakest field strength, respectively. We remove galaxy G09.DR1.328 from the remainder of the analysis, since the estimated L_CO/L_FIR of 2 × 10^-8 (from 5 × the RMS) places the galaxy far outside the parameter space explored in . Before continuing further, several adjustments to our observations are necessary to facilitate a proper comparison with the PDR model.

4.1. Adjustments to observed quantities

In order to draw a direct comparison between our observations and the PDR model of Kaufman et al. (1999, 2006), we must first make several adjustments to the observed [Cii] and CO line emission.

Firstly, the [Cii] emission originates from both ionised and neutral gas, owing to the lower ionisation potential of carbon (11.26 eV) with respect to hydrogen (13.6 eV). Because the PDR models consider that the [Cii] emission originates purely from the neutral gas, we must therefore take the fraction of the [Cii] emission arising from the ionised gas into account. A direct method to determine this fraction is by comparing ratios of the [Cii] 158 μm and the [Nii] 122, 205 μm fine-structure emission lines, particularly the [Cii]/[Nii]205 and [Nii]122/[Nii]205 ratios. The [Nii]122/[Nii]205 ratio is a sensitive probe of the ionised gas density in Hii regions, since the ionisation potential of nitrogen (14.5 eV) is greater than that of neutral hydrogen (13.6 eV). Because the [Cii] and [Nii] 205 μm lines have very similar critical densities for collisional excitation by electrons (46 and 44 cm^-3 at T_e = 8000 K, respectively), the [Cii]/[Nii]205 line ratios are mainly dependent on the relative abundances of C and N in the Hii regions. The ionised gas density can thus be inferred from the theoretical [Nii]122/[Nii]205 ratio, and used to predict the theoretical [Cii]/[Nii]122 ratio arising from the ionised gas and, subsequently, estimate the neutral gas contribution to the [Cii] emission (for details see Oberst et al. 2006, 2011). However, the lack of observations targeting these [Nii] transitions for our sample means we cannot exploit this method in this work. Instead, we must adopt the correction factors obtained from similar previous studies.

Using the direct method, Oberst et al. (2006) found that ~73% of the observed [Cii] line emission of the star-forming Carina nebula in the Galaxy arises from neutral gas in PDRs. Modelling of Hii and PDR regions of starburst galaxies NGC 253 and M 82 (Carral et al. 1994; Lord et al. 1996; Colbert et al. 1999) has also shown that PDRs account for ~70% of the [Cii] emission, with similar results (70–85%) reported by Kramer et al. (2005) in their study of M 51 and M 83. Based on this evidence, numerous studies of higher redshift (z ~ 1–6) IR-bright galaxies thus adopt a 0.7 correction factor to account for the ionised gas contribution to the [Cii] emission (see e.g. Stacey et al. 2010a; Hailey-Dunsheath et al. 2010; Gullberg et al. 2015). Since our sample comprises luminous (L_FIR ~ 10¹⁰–10¹¹L_⊙), actively star-forming (SFR ≈ 40 M_⊙ yr^-1) main-sequence galaxies (see Villanueva et al., in prep.), we also primarily assume ~70% of the observed [Cii] line emission (hereafter ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$) arises from PDRs. We note, however, the broader range of values found in the literature. The survey of Malhotra et al. (2001), for example, found that about 50% of the observed [Cii] emission in their sample of galaxies originated in PDRs when using the [Cii]/[Nii]205 emission combined with the [Nii]122/[Nii]205 ratio of our Galaxy. Spatially resolved studies with the Herschel Very Nearby Galaxies Survey have also demonstrated the varying distributions of neutral and ionised gas within local galaxies using the method of Oberst et al. (2006, 2011) with both direct (e.g. Parkin et al. 2013, 2014; Hughes et al. 2015) and inferred (Hughes et al. 2015) measurements of the [Nii] fine-structure lines. Within M 51, the average ionised gas contribution to the [Cii] emission is ~50% in the spiral arm and interarm regions and reaches up to 80% in the nucleus (Parkin et al. 2013). In NGC 891, although the diffuse neutral component dominates the [Cii] emission in extraplanar regions, regions within the disk have up to 50% of the [Cii] emission originating from ionised gas (Hughes et al. 2015). Thus, at present it is still unclear what fraction of the [Cii] line emission arises from the ionised medium in galaxies, and a constant global correction is clearly not appropriate in all sources. We therefore also choose to examine the results of the PDR modelling when we assume that ~50% of the observed [Cii] line emission (hereafter ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{50}\mathrm{\%}}$) comes from PDRs for all galaxies.

Secondly, we apply a correction to the CO(1–0) emission to account for the likely case that this line becomes optically thick in dense star-forming regions much faster than the [Cii] line or the total infrared flux (see e.g. Stacey et al. 1983b; Tielens & Hollenbach 1985). The PDR infinite plane slab experiences an incident radiation field from one side, whereas each cloud in the ensemble of clouds in extragalactic sources falling within the ALMA beam is not orientated with the irradiated side facing towards our line of sight. We may therefore observe all the optically thin [Cii] and FIR emission but overlook some of the optically thick CO(1–0) line emission that escapes from those clouds with their irradiated sides orientated in the opposite direction from us. Previous studies make the assumption that we only observe about half of the total CO(1–0) emission from all PDRs and thus multiply their observed CO(1–0) emission by a factor of two. In the similar case of comparing the [Cii] 158 and [Oi] 63 μm, Stacey et al. (2010b) reason that such geometric issues should be accounted for by as much as a factor of four for high optical depth in a spherical cloud geometry. Here, we increase our observed CO(1–0) emission by a factor of two for all galaxies, but keep in mind this may be a conservative correction. We refer to the combination of ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$ and 2 × CO as standard adjustments throughout the following text.

Finally, some previous studies also include a correction to the total infrared flux from extragalactic sources. These works argue that because the PDR model assumes the L_FIR emission originates purely from the front side of the cloud, the observations must be reduced by a factor of two to account for the optically thin infrared continuum flux emitting not just towards the observer but from both sides of the PDR slab (see e.g. Parkin et al. 2013; Hughes et al. 2015) and, thereby, make the L_FIR emission directly equivalent to the bolometric far-infrared flux of the PDR model (see Kaufman et al. 1999). We do not apply this correction here, choosing instead to follow the methodology of Stacey et al. (2010a) and Gullberg et al. (2015) to facilitate a comparison with these high-z studies. As such, the possible contamination from an unknown fraction of L_FIR that arises from ionised gas remains the main uncertainty in the FIR emission.

4.2. Insights from the L_[CII] – L_CO diagnostic diagram

Following these adjustments, we now return to to examine their effects on the observations in the L_[Cii] versus L_CO parameter space. Focussing first on our ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$-based correction (, semi-filled red circles), we see that the comparison of these ratios to the contours on the diagnostic diagram indicates the majority of the galaxies have hydrogen nuclei densities in the range of 4 < log n/ cm³ < 5, and typically experience incident FUV radiation fields with strengths varying between log G₀ ≈ 2.0 up to 3.0, although this increases towards log G₀ = 3.5 for some objects. Those galaxies previously falling outside of the L_[Cii]/L_FIR versus L_CO/L_FIR parameter space defined by the PDR model now lie within these same broad parameter ranges, although they appear to contain the least dense gas (log n/ cm³ ~ 3.5) in comparison to the rest of the sample. When we consider the effect of the ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{50}\mathrm{\%}}$-based correction (, open red circles), we find that the expected overall shift of the data points to lower values of L_[Cii]/L_FIR translates into a slight increase in hydrogen nuclei densities, with an approximate range of 4.25 < log n/ cm³< 5.5. In this region of the model parameter space, the FUV radiation field strengths depend primarily on the L_CO/L_FIR ratio and so there is no significant difference in our G₀ values when adopting either ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$ or ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{50}\mathrm{\%}}$ as the diagnostic quantity.

Despite the high quality of our ALMA CO(1–0) and Herschel-PACS [Cii] line observations, it is still difficult to accurately constrain the PDR model parameters using just two lines – one ratio – as diagnostics of the gas conditions due to degeneracies between the line ratio and model parameter space. With access to more observations of far-infrared fine-structure lines (e.g. [Oi] 63 or 145 μm) or lines arising from higher J transitions of the CO molecule, it would be possible to determine more accurately the corresponding best-fit n and G₀ values from the model by fitting the line ratios (using e.g. the PDRT online). Lacking this option, we nevertheless attempt to quantify the mean physical conditions in each galaxy by comparing the observations to the model parameter space. We use a KDTree search coupled with a χ² minimisation technique to find the closest contour of the models to the [Cii]/CO(1–0) observations and assign the corresponding n and G₀ as the best-fitting values. We estimate the errors on the best-fitting model parameters via a bootstrap technique. After the best-fit n and G₀ values are determined, 1000 new sets of data points are created by generating random flux values from within the error bars of the observed or adjusted [Cii], CO(1–0) and FIR flux measurements. Each new dataset is then refit to find the alternative best-fitting parameter values. We then calculate the 68% interval in both the upper and lower parameter distributions and set these intervals as the new upper and lower limits. The differences in the parameter values of the original best-fit solution and the extreme values generated from the bootstrap technique are taken as the uncertainties in the best-fit PDR model parameters. This technique is satisfactory for assigning approximate PDR parameter values and their uncertainties based on the diagnostic diagram in , yielding typical errors of ~0.1–0.3 dex in log n and ~0.1 dex in log G₀. However, a robust quantitative analysis remains limited by uncertainties in the coefficients used to adjust the observations (that are not considered with the bootstrap technique), the degeneracies in the model parameter space, and the fact we are treating each complete galaxy as a single PDR component. We discuss each of these issues in detail later.

The results of our model fitting considering both the original and adjusted observations are presented in . For each case, the parameter ranges of course remain as discussed above. The unadjusted observations yield a mean incident field strength of 2.63 ± 0.32 and density of 4.00 ± 0.95 in logarithmic space. Using these best-fit G₀ and n values, we predict the temperatures of the hot surface layer of the PDR atomic gas that faces towards the FUV source (the surface temperature; see Fig. 1 in Kaufman et al. 1999), which range from approximately 100 to 260 K with a mean temperature of 181 K. In reality, the majority of the PDR gas will be much cooler than this surface temperature. Equating the gas pressure to the product of the density and temperature, i.e. P = nT with units of K cm^-3, as in Malhotra et al. (2001), yields a range of 2 × 10⁴ < P < 24 × 10⁶ K cm^-3 with a mean of 5 × 10⁶ K cm^-3. On the other hand, the physical conditions of the gas traced with ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$ and CO adjusted upward by a factor of two are on average denser (log n = 4.52 ± 0.49), warmer (T = 207 ± 68 K), and more pressurised (47 × 10⁶ K cm^-3) than those results from the unadjusted observations, despite experiencing radiation field strengths of the same order of magnitude (log G₀ = 2.25 ± 0.25). Similar results are found when using ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{50}\mathrm{\%}}$ as the PDR diagnostic (cf. ).

5. Discussion

To summarise the results from a comparison of our adjusted observations to the Kaufman et al. (1999, 2006) PDR model predictions, we find that, when using standard adjustments to observable quantities, the majority of the galaxies in our sample have hydrogen nuclei with densities ranging from 4 < log n/ cm³ < 5.5 with a mean of log n/ cm³ ~ 4.9, and experience an incident FUV radiation field with a strength between 2 < log G₀ < 3 normalised to the Habing (1968) Field (see ). These correspond to mean gas temperatures of 207 K with pressures of 47 × 10⁶ K cm^-3. However, these results rely heavily on certain assumptions made in our analysis. Before comparing our results to those of previous studies in the literature, we first discuss the effects of these assumptions on our results.

5.1. Assumptions and uncertainties

Firstly, there are uncertainties in the PDR parameters associated with our choice of PDR model. Röllig et al. (2007) performed a detailed comparison of PDR models to examine the effects on the physical properties and chemical structures of the model clouds when using different PDR model codes in the literature. An important feature of a PDR model is the adopted geometry of the PDR region; the plane-parallel geometry of the Kaufman et al. (1999, 2006) model is a first order approximation and it is likely that a spherical model could be a more realistic approximation of the complex geometry of the PDR regions. Whilst the benchmarking exercise demonstrated that resulting trends in physical parameters are consistent between the participating codes, they warn that discrepancies remain between the observables computed with different codes (e.g. the atomic fine-structure line intensities) and that these uncertainties should be kept in mind when comparing PDR model results to observations in order to constrain physical parameters, such as density, temperature and radiation field strength.

In applying the PDR model of Kaufman et al. (1999, 2006), we adjusted the observations by, first, correcting the [Cii] 158 μm line emission to remove the contribution to the emission arising from diffuse ionised gas and, second, increasing the CO(1–0) emission by a factor of two to account for the likely case that this line becomes optically thick in dense star-forming regions. We assume these adjustments are correct to the first order for facilitating a proper comparison of our observations and the model. As we previously mention (), without an observational constraint of the fraction of the [Cii] line emission arising from the ionised medium in each of our galaxies, we must use evidence in the literature for global (e.g. Carral et al. 1994; Lord et al. 1996; Colbert et al. 1999; Malhotra et al. 2001) and resolved observations of galaxies (e.g. Parkin et al. 2013, 2014; Hughes et al. 2015) and assume that ~50–70% of the observed [Cii] line emission comes from PDRs for all galaxies. Although this is a standard assumption (see e.g. Stacey et al. 2010a; Hailey-Dunsheath et al. 2010; Gullberg et al. 2015), applying a constant global correction will not be appropriate for all galaxies and L_[Cii] calculated from ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{50}\mathrm{\%}}$ should be considered a lower limit (thus an upper limit in gas density). In the case of the correction to the CO(1–0) emission, there remains a possibility that the line intensity should be further corrected by a factor of two (i.e. a factor of four in total) for high optical depth in a spherical cloud geometry, in contrast to the infinite slab geometry of the PDR model (Stacey et al. 2010a). This correction would shift our galaxies to very low FUV radiation field strengths (log G₀< 1.0) not typically seen in nearby or higher redshift galaxies. A similar effect would be seen if optical depth effects become increasingly important in the handful of high-inclination systems in our sample (see e.g. NGC 891; Hughes et al. 2015).

Another possible correction would be to reduce the FIR emission by a factor of two, in order to account for the optically thin continuum flux emitting not just towards the observer but from both sides of the PDR slab (see ), since the PDR model assumes the L_FIR emission originates purely from the front side of the cloud. In addition, it may also be necessary to reduce the TIR emission further to account for continuum emission from other non-PDR sources, such as Hii regions and ionised gas (see e.g. Croxall et al. 2012, who assume 70% of the TIR emission originates from cool diffuse gas). Whilst we do not make such adjustments in this work, any reduction to L_FIR would shift our ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$- and ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{50}\mathrm{\%}}$-based adjusted observations upward and to the right in the L_[Cii]/L_FIR versus L_CO/L_FIR parameter space in the diagnostic diagram (and also in ), shifting n and G₀ to lower densities (log n/ cm³ ≈ 4.0) and lower FUV field strengths (log G₀ ≈ 1.5).

We further caution that additional shifts in the diagnostic diagrams may arise because of uncertainties with the ALMA and Herschel flux calibration, and small, unknown offsets in the spatial coincidence of the [Cii], CO(1–0) and dust emission. At present, we expect such random shifts to be much smaller, secondary effects compared to the significant uncertainties in the observable quantities and model predictions discussed above. Although there are thus many uncertainties that apply to galaxies on an individual basis, we expect that the systematic adjustments we make are correct on average such that the dervied G₀ and n values are representative. All of these limitations should be kept in mind throughout the following discussion.

Fig. 4 Diagnostic diagram of the L[Cii]/LFIR vs. LCO/LFIR ratios (see also caption). We superimpose our adjusted observations, assuming 70% of the total [Cii] emission arises from PDR regions (red semi-open circles), onto the grid of constant hydrogen nuclei density, log n (black solid lines), and FUV radiation field strength, log G0 (black dotted lines), determined from the PDR model of Kaufman et al. (1999, 2006). Downward triangles represent upper limits. The comparative galaxy samples of Hailey-Dunsheath et al. (2010), where this figure is adapted from, are also shown with the corresponding symbols in the legend. The high-z submillimetre galaxies (blue stars) are from Gullberg et al. (2015). The grey shaded region denotes the 1σ spread in the mean observed L[Cii]/LCO ratio (blue dashed line).

5.2. Comparison with previous studies

As a first step in comparing the physical conditions derived for our sample to those in the literature, in we present the L_[Cii]/L_FIR versus L_CO/L_FIR diagnostic diagram with the contours of constant hydrogen nuclei density and FUV radiation field strength predicted by the Kaufman et al. (1999, 2006) PDR model (see ). We superimpose our adjusted observations, assuming 70% of the total [Cii] emission arises from PDR regions, together with the galaxy samples compared in Fig. 3 of Hailey-Dunsheath et al. (2010), from where this figure is taken. To the first order, our galaxies fall in the region of the parameter space (i.e. 2 ≤ log G₀ ≤ 3 and 3 ≤ log n/ cm³ ≤ 5) predominantly populated by normal galaxies and non-starburst nuclei, in contrast to the Galactic star-forming regions and starburst nuclei that can exhibit higher radiation field strengths (G₀ > 10³) and slightly higher densities (n > 10⁵ cm^-3).

The gas properties we find in our sample are also consistent with other previous surveys of global, integrated observations not present in the figure of Hailey-Dunsheath et al. (2010), although the hydrogen nuclei density tends to be higher on average than the literature values for nearby galaxies. For example, the Malhotra et al. (2001) ISO survey uncovered densities of 2 ≤ log n/ cm³ ≤ 4.5, radiation field strengths of 2 ≤ log G₀ ≤ 4.5, and surface temperatures of ~270–900 K. Similar results have since been uncovered in numerous other samples (cf. Table 9 in Parkin et al. 2013; see also e.g. Lebouteiller et al. 2012). The detailed view afforded by spatially resolved studies of nearby objects also indicates that our sample consists of normal star-forming galaxies. In particular, the Very Nearby Galaxy Survey (VNGS; P. I.: C. D. Wilson; see e.g. Parkin et al. 2013, 2014; Schirm et al. 2014; Hughes et al. 2014, 2016), a Herschel Guaranteed Time Key Project that aims to study the gas and dust in the ISM of a diverse sample of 13 nearby (D_L < 90 Mpc) galaxies, has provided several spatially resolved studies. In the case of M 51 (Parkin et al. 2013), the spiral arm and interarm regions exhibit hydrogen densities and FUV radiation field strengths of 2.75 ≤ log n ≤ 3 cm^-3 and 2.25 ≤ log G₀ ≤ 2.5, respectively, and typically reach higher ranges of n and G₀ for both the central (3 ≤ log n/ cm³ ≤ 3.5, 2.75 ≤ log G₀ ≤ 3) and nuclear (3.75 ≤ log n/ cm³ ≤ 4, 3.25 ≤ log G₀ ≤ 3.75) regions. Comparable distributions in the physical conditions are also seen in the disks of Centaurus A (Parkin et al. 2014) and NGC 891 (Hughes et al. 2015), although the high inclination of the latter system presents an additional challenge in their interpretation. The Kingfish team also report that the PDR gas in NGC 1097 and NGC 4559 has densities between 10^2.5 and 10³ cm^-3 whereas G₀ ranges from 50 to 1000 when adopting conservative corrections for ionised gas (see Fig. 10 of Croxall et al. 2012).

In the case of higher redshift studies, Gullberg et al. (2015) present [Cii] observations of gravitationally lensed, dusty star-forming galaxies in the redshift range z ~ 2.1–5.7 discovered by the SPT (Carlstrom et al. 2011). A subsample also has CO(1–0) observations (Aravena et al. 2016). Using the L_[Cii]/L_FIR versus L_CO/L_FIR diagnostic diagram (see their Fig. 12), the SPT sample observations suggest the radiation field strength and average gas density to be in the range 100 <G₀ < 8000 and 10² <n< 10⁵ cm^-3 for unadjusted observations. Applying the adjustments to obtain 2 × CO and ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$ for these galaxies increases the gas density range by an order of magnitude. In addition to the SPT galaxies, Gullberg et al. (2015) also assemble a high redshift (2.33 < z < 6.42) sample of galaxies in the literature with available measurements of the [Cii] 158 μm and CO(1–0) line emission and FIR luminosities (see their Table B1 and accompanying references therein). On the L_[Cii]/L_FIR versus L_CO/L_FIR diagnostic diagram, these systems have observed ratios corresponding to field strengths of 2 < log G₀< 4 and densities of 2 < log n/ cm³< 6, which translate to much denser PDR regions (~10⁵ cm^-3) that experience similar field strengths if we apply the adjustments of 2 × CO and ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$, as detailed above. In comparison, our sample galaxies therefore typically have lower density gas encountering weaker FUV fields than these higher redshift SPT systems.

Fig. 5 Average G0 (left panel) and n (right panel) plotted as a function of galaxy redshift. We compare the values derived from our unadjusted observations (blue solid circles) and observations including the adjustments applied to the [Cii] and CO(1–0) emission (see ), assuming either 70% (red semi-open circles) or 50% (red open circles) of the total [Cii] emission arises from PDR regions. Downward triangles represent upper limits. Comparative galaxy samples of the literature are shown with the symbols specified in the legends, together with the ranges seen in spatially resolved studies of M 51 (blue line; Parkin et al. 2013), Cen A (red line; Parkin et al. 2014) and NGC 891 (green line; Hughes et al. 2015). To facilitate a comparison, we apply the same standard adjustments (i.e. adopting ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$ and 2CO) to the high redshift (z> 1) systems as those made for galaxies at lower redshift.

5.3. Is there a cosmic evolution of the physical conditions?

As hinted in the previous section, there appear to be some offsets between our galaxies and those at low and high redshift, which may suggest an evolution in the average gas density of the PDR gas. An evolution of the mean gas density would be an interesting result, providing an explanation for the evolution in the star formation rate density of the Universe, ρ_SFR (see e.g. Madau & Dickinson 2014, and references therein), whereby on average galaxies contain denser neutral gas at higher redshifts that can more readily be converted into stars than in their low redshift counterparts. Given that the star formation activity in the Universe was significantly higher in the past and the FUV radiation field strength within galaxies is determined by the star formation history and the initial mass function, with a small contribution from the cosmic FUV background, we expect to observe an evolution of G₀ with redshift. Additional constraints may also be gained from the evolution of the physical conditions; for example, following standard assumptions that the galaxy mass-to-light ratio (M/L) and surface brightness (L/R²) are constant, that n is a constant fraction of the total galaxy mass M such that n ∝ M/R³, and the FUV radiation field strength scales as G₀ ∝ L, we would expect that the quantity G₀n² is also constant as a function of redshift. A robust refute of this would therefore have important implications for the Tully-Fisher relation (Tully & Fisher 1977), for example, since recent studies show the relation exhibits weak/no evolution with redshift (see e.g. Miller et al. 2012; Molina et al., in prep.).

In , we present the G₀ and n values for each galaxy plotted as a function of galaxy redshift. Whilst we primarily focus on our adjusted observations with adjustments based on the standard assumptions of ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$ and optically thick CO emission (i.e. increasing CO by a factor two), we also plot our original, unadjusted observations together with the adjusted observations via ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{50}\mathrm{\%}}$ to demonstrate the possible ranges of our observations. We compare our data to as many literature values as possible, under the condition that they are derived using the same PDR model, in order to avoid any complications arising from the aforementioned offsets between different PDR codes (e.g. Röllig et al. 2007), and with similar adjustments to account for the fraction of [Cii] from ionised gas. Comparative galaxy samples of the literature are shown with varying symbols in the legends, together with the ranges seen in spatially resolved studies of M 51 (blue line; Parkin et al. 2013), Cen A (red line; Parkin et al. 2014), and NGC 891 (green line; Hughes et al. 2015).

Interestingly, the FUV radiation field strength does not appear to show any apparent overall evolution with redshift; the low redshift (z< 0.05) sample of global integrated observations exhibit an average FUV radiation field strength of ~10³, with most galaxies lying within an order of magnitude above and below this value, and the high redshift (z> 1) systems show a similar mean and distribution around G₀ ~ 10³. Our sample typically lies below ~10³ regardless of whether or not we apply any adjustments to our observations. It is worth noting that to match the mean G₀ of our sample to the low and high redshift values would require a significant reduction (~50%) in both the [Cii] and CO(1–0) emission, whereas most of the uncertainties in our methodology discussed in would in fact lead to an increase of the emission. In contrast, the hydrogen nuclei density appears to show a significant redshift evolution. The mean density of both the integrated and spatially resolved observations of galaxies at low redshift is n ~ 10³. Our galaxies, however, exhibit a higher mean density between n ~ 10⁴–10⁵, depending on the adopted correction (see the right panel of ), that is roughly equivalent to the density distribution of the high-z systems, which suggests a strong evolution of the density up to z = 0.2. The value of G₀n² is not constant, contrary to our prediction, as we observe a similar variation with redshift that is driven by the variation in the gas density.

However, it is not possible at this stage to conclude firmly that this observed evolution is in fact physical and we must be cautious not to confuse clear observational biases with an evolution of the physical conditions. One bias may be introduced from differences in the samples at low and high redshift. The SPT sample comprises submillimetre-selected dusty star-forming galaxies, which are typically massive in both their stellar content (M_∗ ~ 10¹¹M_⊙) and gas content (M_gas ~ 5 × 10¹⁰M_⊙) and have star formation rates of the order of 1000M_⊙ yr^-1 (see Gullberg et al. 2015, and references therein). In contrast, the corresponding key properties of our H-ATLAS-based sample suggest our galaxies are normal star-forming systems (see Villanueva et al., in prep.) and so may not be directly comparable to these high redshift systems that generally lie above the main sequence (see e.g. Elbaz et al. 2011). The potential for misclassification of high-z sources and contamination of the sample (e.g. from AGN hosts) only adds to these difficulties. Nevertheless, we observe an offset between the VALES sample and the low redshift systems that share similar properties and appear comparable (e.g. Malhotra et al. 2001; Negishi et al. 2001; Parkin et al. 2013), but this potential evidence for an evolution in the gas density is also undermined by a bias that may arise from using different suites of emission lines as PDR diagnostics. The results at low redshift are based on a whole suite of FIR fine structure lines (e.g. the [Cii] 158 μm, [Oi] 63 and 145 μm, [Oiii] 88 μm and [Nii] 122 and 205 μm lines) and derivative diagnostic ratios (e.g. [Cii]/[Oi]63 or ([Cii]+[Oi]63)/L_TIR) that may be probing a different density regime of the PDR regions compared to the higher redshift observations that rely on the [Cii] and CO(1–0) line emission alone. Furthermore, we note that very nearby galaxy studies tend to be more hand-picked and less statistical in their selection, thus may not be appropriate for comparison to IR-selected high-z samples. The difference between the nearby galaxies and our sample is most likely owing to differences in methodology and/or sample selection and, for these reasons, we recommend our VALES sample as a better local reference for more distant galaxies when using the [Cii] and CO(1–0) emission lines.

At present, we therefore cannot draw strong conclusions on whether the gas density evolves with redshift. Placing a more robust constraint on the evolution in the average G₀ and n values would require additional gas diagnostics via observations of other far-infrared fine-structure lines in higher redshift systems. For example, the ratio of the two submillimetre C⁰ transitions, [Ci] 370 μm/[Ci] 609 μm, is particularly useful for tracing the gas density for the following reasons (see e.g. Kaufman et al. 1999): (1) a weak temperature dependence due to low upper-state energies for both transitions (E₆₀₉/k ~ 24 K, E₃₇₀/k ~ 63 K); (2) an insensitivity to the FUV radiation field strength in temperatures warm enough to excite both transitions; (3) similar critical densities (n_crit([Ci] 370 μm) ~ 3 × 10² cm^-3 and n_crit([Ci] 609 μm) ~ 2 × 10³ cm^-3); and (4) an insensitivity to the C⁰ abundance. For high values of the FUV radiation field strength per hydrogen nuclei density (G₀/n), the FUV heating of the gas occurs at A_V< 1 and the [Ci] lines are generally weaker than the [Cii] and [Oi] lines, whereas C⁰ lines tend to dominate gas cooling at small column densities in clouds with low values of G₀/n, i.e. 3 × 10^-3. The [Ci] lines are therefore useful for determining the details of the chemistry and penetration of the incident FUV in photodissociation regions. Of course, the [Nii] 122 and 205 μm fine-structure emission lines for constraining the contribution to the [Cii] emission from ionised gas would also be valuable.

6. Conclusions

We have used new ALMA CO(1–0) observations to study the physical conditions in the interstellar gas of a sample of 27 dusty main-sequence star-forming galaxies at 0.03 <z< 0.2, bridging the gap between local and high-z galaxy samples. Our sample was drawn from far-IR bright galaxies over ~160 deg² in the Herschel Astrophysical Terahertz Large Area Survey. We combined these observations with high-quality ancillary data, including Herschel-PACS [Cii] 158 μm spectroscopy and far-infrared luminosities determined with photometry from Herschel and other facilities. Our main results and conclusions are:

• 1.

  The average line luminosity ratio of the [Cii] andCO(1–0) detections is 3500 ± 1200 for our sample, in agreement with previous works. The observed line emission is also consistent with that of both local and high-z galaxies. The observed L_[Cii]/ L_CO ratios appear to arise from gas with T_{ex, [Cii]} > T_ex,CO(1−0), where the [Cii] is optically thin and the CO is optically thick.

• 2.

  Our sample covers the same distribution in the L_[Cii]/L_FIR and L_CO/L_FIR ratio as normal star-forming galaxies and non-starburst nuclei.

• 3.

  A comparison with the L_[Cii]/L_FIR versus L_CO/L_FIR diagnostic diagram in provides a first order estimate of the radiation field strength and average gas density found to be in the range of 2 < log G₀ < 3 and 2 < log (n/ cm³) < 4.5. Assuming standard adjustments increases the gas density by an order of magnitude. These values are consistent with previous surveys employing either integrated or spatially resolved observations, although the hydrogen nuclei density tends to be higher on average than the literature values for observations of nearby galaxies.

• 4.

  The average FUV radiation field strength appears constant up to redshift z ~ 6.4, yet we find an increase of a factor of ~100 in the neutral gas density as a function of redshift that persists regardless of various adjustments to our observable quantities. The apparent evolution most likely arises from a combination of observational biases rather than physical processes and we highlight the need to take a consistent methodology between different studies.

Of the many uncertainties in this work, the major limitation is the use of only two emission lines – one ratio – as diagnostics of the physical conditions of the gas. Future studies with ALMA could observe these galaxies using other fine-structure lines, particularly the two submillimeter C⁰ transitions, [Ci] 370 and 609 μm, the ratio of which is useful for tracing the gas density and determining the details of the chemistry and penetration of the incident FUV in photodissociation regions. The capability of ALMA to resolve spatially the [Cii] emission and other fine structure lines in galaxies at higher redshifts will provide in the future a more detailed picture of whether or not the gas density evolves with cosmic time.

Acknowledgments

We thank the anonymous referee for their clear and constructive report. T.M.H. and E.I. acknowledge CONICYT/ALMA funding Program in Astronomy/PCI Project N°:31140020. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00530.S. M.A. acknowledges partial support from FONDECYT through grant 1140099. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The Herschel-ATLAS is a project with Herschel, which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. The H-ATLAS website is http://www.h-atlas.org/. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the NASA (USA).

References

All Tables

All Figures

Fig. 1 [Cii] luminosity vs. the CO(1–0) luminosity for the VALES sample at 0.03 <z< 0.2 superimposed on the low- and high-z samples taken from Fig. 7 of Gullberg et al. (2015), with symbols as specified in the legend. Downward triangles represent 5σ upper limits. The grey shaded region the 1σ spread in the mean observed L[Cii]/LCO ratio (3500 ± 1200; blue dashed line) of our sample, which we compare to the mean ratios of the low-z sample (1300 ± 440; black dotted line).

In the text

Fig. 2 CO(1–0) vs. [Cii] excitation temperatures determined via the ratio of the source functions ()) assuming various line opacities, as specified in the legend, and adopting the observed L[Cii]/LCO ratio of 3500. The light and dark grey shaded regions denote denotes the 1σ spread in the mean observed ratio, whereas the red shaded region represents the CO excitation temperatures of 20 K <Tex,CO(1−0)< 55 K. Combined with the results of PDR modelling and previous observational evidence in the literature (), the best scenario to explain the observed L[Cii]/LCO ratios is gas where Tex, [Cii]>Tex,CO(1−0) and with [Cii] is optically thin and the CO is optically thick (τ[Cii] = 0.1, τCO(1−0) = 1; black dotted line).

In the text

Fig. 3 Diagnostic diagram of the observed L[Cii]/LFIR ratio vs. LCO/LFIR ratio for our sample. The observations are superimposed onto the grid of constant hydrogen nuclei density, log n (black solid lines), and FUV radiation field strength, log G0 (black dotted lines), determined from the PDR model of Kaufman et al. (1999, 2006). The unadjusted observations (blue solid circles) are compared to the adjusted observations, as described in , where the fraction of the [Cii] emission arising from neutral gas is fixed at 70% (red semi-open circles) and 50% (red open circles). Downward triangles represent 5σ upper limits. Two sources lie above the log n = 2 contour, i.e. outside the model parameter space, and the blue arrow points towards G09.DR1.328, which remains outside the plot regardless of adjustments. The grey shaded region denotes the 1σ spread in the mean observed L[Cii]/LCO ratio (blue dashed line). This figure is adapted from Stacey et al. (2010a) and Hailey-Dunsheath et al. (2010).

In the text

Fig. 4 Diagnostic diagram of the L[Cii]/LFIR vs. LCO/LFIR ratios (see also caption). We superimpose our adjusted observations, assuming 70% of the total [Cii] emission arises from PDR regions (red semi-open circles), onto the grid of constant hydrogen nuclei density, log n (black solid lines), and FUV radiation field strength, log G0 (black dotted lines), determined from the PDR model of Kaufman et al. (1999, 2006). Downward triangles represent upper limits. The comparative galaxy samples of Hailey-Dunsheath et al. (2010), where this figure is adapted from, are also shown with the corresponding symbols in the legend. The high-z submillimetre galaxies (blue stars) are from Gullberg et al. (2015). The grey shaded region denotes the 1σ spread in the mean observed L[Cii]/LCO ratio (blue dashed line).

In the text

Fig. 5 Average G0 (left panel) and n (right panel) plotted as a function of galaxy redshift. We compare the values derived from our unadjusted observations (blue solid circles) and observations including the adjustments applied to the [Cii] and CO(1–0) emission (see ), assuming either 70% (red semi-open circles) or 50% (red open circles) of the total [Cii] emission arises from PDR regions. Downward triangles represent upper limits. Comparative galaxy samples of the literature are shown with the symbols specified in the legends, together with the ranges seen in spatially resolved studies of M 51 (blue line; Parkin et al. 2013), Cen A (red line; Parkin et al. 2014) and NGC 891 (green line; Hughes et al. 2015). To facilitate a comparison, we apply the same standard adjustments (i.e. adopting ${\mathrm{\left[}\mathrm{C}{ii}\mathrm{\right]}}_{{PDR}}^{\mathrm{70}\mathrm{\%}}$ and 2CO) to the high redshift (z> 1) systems as those made for galaxies at lower redshift.

In the text