Deriving a multivariate CO-to-H$_2$ conversion function using the [CII]/CO(1-0) ratio and its application to molecular gas scaling relations

TL;DR

This study develops a multivariate CO-to-H2 conversion function based on the [CII]/CO(1-0) ratio, incorporating metallicity and star formation activity, to improve molecular gas estimates across different galaxy types and redshifts.

Contribution

It introduces a new multivariate conversion function for CO-to-H2 that accounts for metallicity and star formation activity, applicable up to redshift 2.5.

Findings

The [CII]/CO ratio depends mainly on metallicity and star formation offset.

The new conversion function reduces uncertainties in molecular gas mass estimates.

Low mass galaxies' low CO luminosity is due to photodissociation, not star formation efficiency.

Abstract

We present Herschel PACS observations of the [CII] 158 micron emission line in a sample of 24 intermediate mass (9<logM$_\ast$/M$_\odot$<10) and low metallicity (0.4< Z/Z$_\odot$<1.0) galaxies from the xCOLD GASS survey. Combining them with IRAM CO(1-0) measurements, we establish scaling relations between integrated and molecular region [CII]/CO(1-0) luminosity ratios as a function of integrated galaxy properties. A Bayesian analysis reveals that only two parameters, metallicity and offset from the star formation main sequence, $\Delta$MS, are needed to quantify variations in the luminosity ratio; metallicity describes the total dust content available to shield CO from UV radiation, while $\Delta$MS describes the strength of this radiation field. We connect the [CII]/CO luminosity ratio to the CO-to-H$_2$ conversion factor and find a multivariate conversion function $\alpha_{CO}$, which can be used up to z~2.5. This function depends primarily on metallicity, with a second order dependence on $\Delta$MS. We apply this to the full xCOLD GASS and PHIBSS1 surveys and investigate molecular gas scaling relations. We find a flattening of the relation between gas mass fraction and stellar mass at logM$_\ast$/M$_\odot$<10. While the molecular gas depletion time varies with sSFR, it is mostly independent of mass, indicating that the low L$_{CO}$/SFR ratios long observed in low mass galaxies are entirely due to photodissociation of CO, and not to an enhanced star formation efficiency.

Figures

Tables

Table 1. Table 1: Derived observational quantities for the xCOLD-GASS galaxies.

xCOLD GASS ID	Log M∗	Redshift	SFR	Metallicity	Log [CII]	Log CO(1-0)	Log sSFR	f${}_{\text{[C}\text{II}\text{] mol}}$
	[M⊙]		[M⊙ yr-1]	[O/H]	[L⊙]	[L⊙]	[yr-1]	%
107026	9.37	0.0167	0.34	8.55	6.65	2.92	-9.84	67.21
108064	9.56	0.015	0.75	8.64	7.05	3.45	-9.68	64.4
108080	9.31	0.0153	0.05	8.7	5.81	2.99	-10.61	76.8
108093	9.94	0.0168	0.68	8.78	7.48	3.98	-10.11	71.3
108113	9.88	0.0199	5.12	8.79	6.99	3.75	-9.17	53.62
108129	9.47	0.0175	0.19	8.45	6.39	3.11	-10.19	72.37
108142	9.61	0.019	0.71	8.7	7.27	3.36	-9.76	65.84
108147	9.64	0.0197	0.62	8.67	7.14	3.43	-9.85	67.37
101037	9.31	0.0164	0.27	8.61	6.93	3.22	-9.88	67.86
108050	9.75	0.0162	0.76	8.82	7.32	3.76	-9.87	67.7
108021	9.64	0.0173	0.37	8.82	7.2	3.86	-10.07	70.75
109066	9.99	0.0149	0.14	8.71	6.79	3.38	-10.84	78.42
109101	9.01	0.0104	0.15	8.52	6.31	2.31	-9.83	67.04
109038	9.21	0.0116	0.66	8.4	7.24	3.11	-9.39	58.62
109139	9.59	0.0188	0.93	8.78	7.25	3.6	-9.62	63.28
109092	9.7	0.018	0.58	8.68	7.12	3.58	-9.94	68.82
109072	9.68	0.0189	0.73	8.79	7.22	3.6	-9.82	66.87
109010	9.75	0.0102	0.69	8.66	6.87	3.54	-9.91	68.34
109102	9.38	0.0117	0.44	8.33	6.82	2.93	-9.74	65.49
110038	10.0	0.0166	0.51	8.71	7.03	3.72	-10.29	73.6
109109	9.08	0.0196	0.27	8.61	6.63	—	-9.65	63.85
108054	9.93	0.0125	0.65	8.66	6.85	3.65	-10.12	71.44
108045	10.08	0.015	1.76	8.69	7.62	4.22	-9.83	67.04
109028	10.07	0.0178	1.36	8.73	7.56	4.17	-9.94	68.82

Table 2. Table 2: Derived observational quantities for the seven DGS galaxies.

DGS Name	Log M∗	Redshift	SFR	Metallicity	Log [CII]	Log CO(1-0)	Log sSFR	f${}_{\text{[C}\text{II}\text{] mol}}$
	[M⊙]		[M⊙ yr-1]	[O/H]	[L⊙]	[L⊙]	[yr-1]	%
Haro 11	10.24	0.0206	37.53	8.32	8.08	3.7	-8.67	40.3
Mrk 930	9.44	0.0183	2.78	8.08	7.44	2.81	-9.0	49.39
Haro 3	9.49	0.0031	0.55	8.12	6.49	2.42	-9.75	65.66
Mrk 1089	10.28	0.0134	5.72	8.3	7.59	3.8	-9.52	61.32
UM 448	10.62	0.0186	9.53	8.1	8.14	3.93	-9.64	63.66
Haro 2	9.56	0.0048	1.24	8.39	6.89	3.29	-9.47	60.3
II Zw 40	8.6	0.0026	0.57	7.92	6.41	1.92	-8.84	45.13

Equations

$$\text{centering}{M}_{H_2}={\alpha }_{CO}{L}_{CO(1-0)}^{{}^{\prime }}.\text{@add@centering}$$

$$\text{centering}{L}_{CO(1-0)}^{{}^{\prime }}=3.25\times 10^7{I}_{CO}\Delta v{\nu }_{obs}^{-2}{D}_L^2(1+z{)}^{-3}\text{@add@centering}$$

$$12+\log (\text{mbox}O/\text{mbox}H)=8.73-0.32\log \left(\frac{\text{mbox}[OIII]}{\text{mbox}{H}_{\beta }}/\frac{\text{mbox}[NII]}{\text{mbox}{H}_{\alpha }}\right)$$

$$\centering\mu_{*}=\frac{\mbox{M${}_{*}$}}{2\pi\mbox{R}^{2}_{50,z}}\@add@centering$$

$$\text{centering}{A}_{\text{mbox}IRX}={\left(\frac{\text{mbox}{SFR}_{\text{mbox}IR}}{\text{mbox}{SFR}_{\text{mbox}UV}}+1\right)}^{2.5}\text{@add@centering}$$

$$\text{centering}\log (\text{mbox}{sSFR}_{ms}(z,M_{\ast }))=-1.12+1.14z-0.19z^2-(0.3+0.13z)\times (\log M_{\ast }-10.5)[\text{mbox}{Gyr}^{-1}]\text{@add@centering}$$

$$\text{centering}\Delta (\text{mbox}MS)=\frac{\text{mbox}{sSFR}_{\text{mbox}measured}}{\text{mbox}{sSFR}_{ms}(z,M_{\ast })}\text{@add@centering}$$

$$\text{centering}{ϵ}_{obs}=\frac{{\sigma }_{rms}W50_{CO}}{\sqrt{W50_{CO}\Delta w_{ch}^{-1}}}\text{@add@centering}$$

$$\text{centering}{f}_{\text{mbox}[CII],mol}=-6.224-1.235\psi -0.0543{\psi }^2\text{@add@centering}$$

$$\text{centering}{y}_i=\mathcal{N}(y_{true,i},{\sigma }_{y_i}^2)\text{@add@centering}$$

$$\text{centering}{y}_i=10^{\alpha }(x_i^j{)}^{{\beta }_j}\text{@add@centering}$$

$$\text{centering}P(y_i|y_{true,i},{\sigma }_{y_i})=\frac{1}{\sqrt{2\pi {\sigma }_{y_i}^2}}\exp \left(-\frac{(y_i-10^{\alpha }(x_i^j{)}^{{\beta }_j}{)}^2}{2{\sigma }_{y_i}^2}\right)\text{@add@centering}$$

$$\text{centering}\mathcal{L}=-\frac{N}{2}\text{mbox}ln(2\pi )-\frac{1}{2}\sum _{i=1}^N\text{mbox}ln({\sigma }_{y_i}^2)-\frac{1}{2}\sum _{i=1}^N\frac{(y_i-10^{\alpha }(x_i^j{)}^{{\beta }_j}{)}^2}{{\sigma }_{y_i}^2},\text{@add@centering}$$

$$\text{centering}AIC=-2\mathcal{L}+2p+\frac{2p(p+1)}{N-p-1}\text{@add@centering}$$

$$\text{centering}BIC=-2\mathcal{L}+p\text{mbox}Log(N)\text{@add@centering}$$

$$\log \left(\frac{\text{mbox}[CII]}{\text{mbox}CO}\right)(\pm 0.223\text{mbox}dex)=11.204+0.230\log \Delta (\text{mbox}MS)-0.905[\text{mbox}12+\log (\text{mbox}O/H)]$$

$$\log \left(\frac{\text{mbox}[C\text{sc}ii]}{\text{mbox}CO}\right)=12.21+0.266\log \Delta (\text{mbox}MS)-1.0189[\text{mbox}12+\log (\text{mbox}O/H)],$$

$$\text{centering}\text{mbox}CO+\gamma \to \text{mbox}O+\text{mbox}C\text{@add@centering}$$

$$\text{centering}\text{mbox}C+\gamma \to \text{mbox}{C}^{+}+\text{mbox}{e}^{-}\text{@add@centering}$$

$$\text{centering}\text{mbox}He+\text{mbox}CRP\to \text{mbox}{He}^{+}+\text{mbox}{e}^{-}\text{@add@centering}$$

$$\text{centering}\text{mbox}{He}^{+}+\text{mbox}CO\to \text{mbox}O+\text{mbox}{C}^{+}+\text{mbox}He\text{@add@centering}$$

$$\text{centering}\text{mbox}{He}^{+}+\text{mbox}CO\to \text{mbox}{O}^{+}+\text{mbox}C+\text{mbox}He\text{@add@centering}$$

$$\text{centering}\text{mbox}CO+\text{mbox}CRP\to \text{mbox}O+\text{mbox}C.\text{@add@centering}$$

$$\text{centering}\log {\alpha }_{\text{mbox}CO}=0.742\log \frac{L_{\text{mbox}[CII]}}{L_{\text{mbox}CO(1-0)}}-0.944[12+\log (\text{mbox}O/\text{mbox}H)]-0.109\log \Delta (\text{mbox}MS)+6.439\text{@add@centering}$$

$$\text{centering}\log {\alpha }_{\text{mbox}CO}(\pm 0.165\text{mbox}dex)=14.752-1.623[12+\log (\text{mbox}O/\text{mbox}H)]+0.062\log \Delta (\text{mbox}MS),\text{@add@centering}$$

$$\text{centering}\log {\alpha }_{\text{mbox}CO}(\pm 0.165\text{mbox}dex)=15.623-1.732[12+\log (\text{mbox}O/\text{mbox}H)]+0.051\log \Delta (\text{mbox}MS).\text{@add@centering}$$

Full text

Accepted Date Received Date; in original form Date

1 Introduction

Observations of the cold interstellar medium (ISM) of normal star-forming galaxies at low and high redshifts have recently helped to establish an “equilibrium model” for galaxy evolution. Under this framework, galaxy growth is controlled by the total available gas reservoir through the interplay of gas inflows and outflows, and the efficiency of the star formation process (e.g. White & Frenk, 1991; Bouché et al., 2010; Davé et al., 2012; Lilly et al., 2013). Trying to understand what triggers and drives star formation, how the chemical enrichment of the ISM proceeds, and what causes the growth of galaxies is intricately linked to understanding their total gas content. Being able to measure the amount and properties of cold gas in a very wide range of galaxies is therefore of great importance to further our understanding of galaxy evolution.

It has been well established that most star formation in the Milky Way and nearby galaxies occurs in dense giant molecular clouds (GMCs), and that most of this molecular gas is in the form of cold H2 (e.g. Kennicutt & Evans, 2012; Solomon et al., 1987; Young & Scoville, 1991; Leroy et al., 2008). However, with the molecule lacking a permanent electric dipole moment, cold H2 is not directly observable, and it is common practice to instead trace cold molecular gas through the low-lying rotational transitions of the second most abundant molecule, 12CO (Scoville, 2013). The ground rotational transition of CO has a low excitation temperature, 5.53K, and a low critical density, 700 cm*-3*, making it easily excited in cold molecular clouds (Dickman et al., 1986). At a wavelength of 2.6mm the ground state falls within the Earth’s atmospheric window allowing it to be easily observed from ground based facilities. Given all this, it has become the workhorse tracer to quantify the total molecular gas reservoir in galaxies both near and distant (Bolatto et al., 2013).

Although the CO rotational transitions are optically thick in typical conditions, the information on the total molecular gas mass is contained in the width of the line under the assumption that GMCs are virialized and that the line emission is the superposition of a number of such virialized clouds. The correlation between the velocity-integrated line luminosity L${}_{CO}^{{}^{\prime}}$ in the $1$ $\rightarrow$ [math] transition and the total molecular gas mass is given by the empirical relation (Dickman et al., 1986; Obreschkow & Rawlings, 2009):

[TABLE]

Here M${}_{H_{2}}$ has units of M*⊙* and L${}_{CO(1-0)}^{{}^{\prime}}$ (K km s*-1* pc2), the integrated line luminosity, is related to the observed velocity integrated flux density, I${}_{CO}\Delta v$ (Jy km s*-1*), following Solomon et al. (1997):

[TABLE]

where the frequency is in GHz (Solomon et al., 1987) and the luminosity distance, DL, is in Mpc. Thus $\alpha_{CO}$, the CO-to-H2 conversion factor, can be considered a mass to light ratio. Across observations of the Milky Way and nearby star-forming galaxies with approximately solar metallicities, the empirical CO(1-0) conversion factors are consistent with a typically value of 4.36 M*⊙* (K km s${}^{-1})^{-1}$, which includes a 36$\%$ correction for Helium gas (Strong & Mattox, 1996; Abdo et al., 2010).

However CO has been difficult to detect in local low mass galaxies. Sensitivity has limited CO detections to galaxies with $Z/Z_{\odot}\geq$ 0.1 (Leroy et al., 2009; Schruba et al., 2012; Rémy-Ruyer et al., 2014). Does this mean that galaxies with even lower metallicities have very little molecular gas, or does CO become a poor tracer of the molecular ISM in low metallicity conditions? Indeed there is strong evidence that not all of the H2 is traced by CO emission; UV radiation from massive stars destroys CO to a cloud depth of a few $A_{V}$, which can correspond to a significant fraction of the total gas column in low metallicity clouds (van Dishoeck & Black, 1986; Wolfire et al., 2010). While H2 is self-shielded from this UV radiation, CO relies on dust shielding and therefore, in low metallicity star forming galaxies which have hard radiation fields and lower dust-to-gas mass ratios, CO is easily photodissociated into C*+* and O (Poglitsch et al., 1995; Röllig et al., 2006).

The regions where molecular hydrogen undergoes a dissociation transition into neutral hydrogen are suitably named photodissociation regions (PDRs hereafter) and it is here that CO is also photodissociated. In this case, the CO flux per fixed hydrogen column is less than in high metallicity environments and, in turn, the Galactic conversion factor globally underestimates the true molecular hydrogen content (Arimoto et al., 1996; Leroy et al., 2011). This molecular gas, not traced by CO emission, has been referred to as the “dark gas” (Wolfire et al., 2010; Grenier et al., 2005) and emits brightly in other fine-structure PDR tracers such as [OI] and [C ii]. In particular, [C ii] is a promising tracer to quantify the total dark molecular gas reservoir. It was first used in a low metallicity dwarf galaxy by Madden et al. (1997) and is the focus of this work.

The [C ii] 158$\mu$m emission line is one of the strongest coolants of the interstellar medium and can contribute up to a few percent of the total FIR emission from a galaxy (Tielens & Hollenbach, 1985). Ionised carbon has a lower ionisation potential than hydrogen (11.3 eV instead of 13.6 eV), and the [C ii] line lies 92K above the ground state with a critical density for collisions with neutral hydrogen of $3\times 10^{3}$ cm*-3* (Kaufman et al., 1999). This means [C ii] is produced not only in PDRs, but also in the ionised and atomic phases of the ISM. Measurements of the [C ii] emission originating from the PDRs combined with measurements of CO (which can only arise in PDRs), can be used a) to trace the total molecular reservoir on galaxy wide scales and b) to quantify the nature of variations of the conversion function. We here use the nomenclature of a conversion ‘function’ because this quantity does vary as a function of ISM properties.

There have been a significant effort, both observationally and theoretically, made to establish a robust prescription for the conversion function (e.g. Wilson, 1995; Israel, 1997; Boselli et al., 2002; Israel & Baas, 2003; Leroy et al., 2011; Genzel et al., 2012; Schruba et al., 2012; Bolatto et al., 2013; Sandstrom et al., 2013). However a consensus has not yet been reached, and some empirical methods rely on the assumption of a universal gas depletion timescale. To disentangle possible systematic variations in depletion time from the issues of CO in low metallicity environments, a possibility is to use multi-wavelength observations to investigate the properties of the ISM (e.g. Magnelli et al., 2012).

By measuring carbon both in its molecular and ionised form in PDRs in a statistically robust sample of galaxies, with a full suite of multi-wavelength observations, it is possible to investigate any secondary dependencies on the conversion function. Such a dataset would allow us to fully parametrise and quantify variations in the conversion function with global galaxy properties. In this paper we present CO(1-0) and [C ii] observations of 24 intermediate mass ($9<\log M_{\ast}/M_{\odot}<10$) and low metallicity ($0.4<Z/Z_{\odot}<1.0$) galaxies from the xCOLD GASS survey. We combine this with spectroscopic CO(1-0) and [C ii] data from the Dwarf Galaxy Survey (DGS hereafter) (Cormier et al., 2014; Madden et al., 2013; Rémy-Ruyer et al., 2014) to probe galaxies that are even more metal-poor.

In Section 2 we present an overview of the survey and details regarding the data reduction of the IRAM CO(1-0), Herschel [C ii] observations and auxiliary data from WISE, GALEX and SDSS. In Section 3 we provide scaling relations between integrated L${}_{\mbox{[C\scriptsize II]}}$/L${}_{\mbox{CO(1-0)}}$ and global galaxy parameters. In Section 3.1 we explore the contribution of the [C ii] emission from non-molecular regions and proceed to use a Bayesian Inference method and radiative transfer modelling to retrieve a prescription for the conversion function in Sections 4 and 5. Finally we discuss the implications of our new prescription on molecular gas scaling relations in Section 6.

Throughout this paper we use a standard flat $\Lambda\mbox{CDM}$ cosmology with $H_{0}=70$ km s*-1* Mpc*-1*.

2 Survey Description and Sample Selection

The extension to COLD GASS is a randomly selected sample of galaxies from the regions of overlap between the SDSS (Stoughton et al., 2002), GALEX (Martin et al., 2005), WISE (Wright et al., 2010) and ALFALFA HI (Giovanelli et al., 2005) surveys. It is an unbiased sub-sample of all the galaxies in the redshift range 0.01$<$ z $<$0.02 and is stellar mass selected (9$<$ log M*∗/M⊙* $<$ 10). There is no other selection criteria based on colour, star formation rate, etc. The SDSS data provide us with optical imaging and spectroscopy over the central 3” of our galaxies and with the GALEX and WISE data we have FUV, NUV, 3.4$\mu$m, 4.5$\mu$m, 12$\mu$m and 22$\mu$m photometry. Moreover, we have HI fluxes from observations at the Arecibo observatory as part of the GASS survey (Catinella et al., 2010; Catinella et al., 2013). The full xCOLD GASS sample contains 133 galaxies, which will be presented in full alongside all the CO(1-0) IRAM observations in Saintonge et al. 2017.

Our target selection for Herschel PACS observations involved removing from the total sample of 133 xCOLD GASS objects the passive elliptical galaxies, where Saintonge et al. (2011a) showed that molecular gas contributes insignificantly to the mass budget (M${}_{\mbox{H_{2}}}$/M*∗* $<$ 0.2%), and one galaxy which already had adequate PACS observations. After this a mass-selected sample of 103 galaxies remain which allows for roughly 20 galaxies in 5 stellar mass bins between 9$<$ log M*∗/M⊙* $<$ 10. However the Herschel telescope exhausted its supply of liquid helium coolant midway through our observations and so, out of those 103 galaxies which were originally proposed, we here present the 24 galaxies for which the [C ii] 158$\mu$m line was observed with PACS.

In order to provide a statistically robust prescription of the conversion function it is imperative to probe as large a parameter space as possible, hence we combine our data with literature data if possible. Unfortunately multi-wavelength data sets with accurate galaxy parameter measurements and observations of both CO and [C ii] are rare. Cormier et al. (2014) do present data for a compilation of galaxies from the DGS which have good CO(1-0) and [C ii] observations, which we therefore add to the xCOLD GASS objects. We are interested in objects with total galaxy-wide integrated detections in CO(1-0) and [CII]; we therefore have to eliminate all objects which have optical diameters greater than 47” (the PACS IFU map size), and objects which have resolved interferometric CO observations as aperture corrections are impossible due to the unresolved nature of the Herschel PACS data.

Overall we are able to add seven extremely low metallicity galaxies, $0.2Z_{\odot}<Z<0.5Z_{\odot}$, from the DGS survey; these are Haro 11, Mrk 930, Haro 3, Mrk 1089, UM 448, Haro 2 and II Zw 40. Derived galaxy parameters are found within Madden et al. (2013) while observed [CII]/CO(1-0) ratios, with aperture corrections, and star formation rates are calculated here. We fold these galaxies into our sample and derive all necessary measurements consistently as detailed below.

2.1 Optical, UV & IR data

For the xCOLD GASS galaxies parameters such as redshifts, sizes and magnitudes are retrieved from the SDSS DR7 database (Abazajian et al., 2009). We retrieve stellar masses and emission line fluxes from the MPA-JHU catalogue111The data catalogues are available from http://www.mpa-garching.mpg.de/SDSS/ where calculations were performed using the methods presented in Tremonti et al. (2004) and Kauffmann et al. (2003). With these, we then calculate the gas-phase metallicity of the galaxies using the prescription from Pettini & Pagel (2004) (PP04 hereafter):

[TABLE]

which was calibrated down to the low metallicity range probed here. On this scale a value of $12+\log(O/H)=8.69$ corresponds to solar metallicity following the oxygen abundance of Asplund et al. (2009). The stellar mass surface density is defined as:

[TABLE]

where R50,z is the z-band 50% flux intensity Petrosian radius, in kiloparsecs, taken from the SDSS database.

For the DGS galaxies stellar masses and redshifts are taken from the catalogue presented in Madden et al. (2013), with an explanation on how these quantities were derived, while metallicities are taken from the literature and converted to PP04 units using the methods from Kewley & Ellison (2008).

To account for unobscured star formation we use GALEX FUV and NUV images retrieved from the public GALEX All-sky and Medium Imaging surveys (AIS and MIS respectively, Martin et al., 2005), and for the obscured star formation we make use of WISE imaging. For both GALEX and WISE maps we perform aperture photometry using a similar data reduction technique to Wang et al. (2010). Total SFRs from the combination of UV and IR photometry are then calculated as presented in Saintonge et al. (2016) and we do this for both the xCOLD GASS and DGS samples.

From this it is possible to calculate, for both xCOLD GASS and DGS galaxies, the effective UV attenuation, AIRX, where the log quantity is defined in Saintonge et al. (2013) as:

[TABLE]

. This is another important quantity which may correlate with L${}_{\mbox{[C\scriptsize II]}}$/L${}_{\mbox{CO(1-0)}}$ - we express it in this form, as opposed to the log quantity, for mathematical convenience when performing the statistical analysis described later in Section 4. We also want to allow for the possibility of a redshift dependence in our conversion function and so we measure the distance off the main sequence for each galaxy. Using the analytical definition of the main sequence by Whitaker et al. (2012), where:

[TABLE]

with $z$ and M*∗* denoting redshift and stellar mass, we can then define the distance off the main sequence as:

[TABLE]

which is applicable up to z$\sim$2.5, as stated in Whitaker et al. (2012).

2.2 DSS observations and data reduction

To measure r-band magnitudes in a homogenous way between the xCOLD GASS and DGS samples we perform aperture photometry here; the DGS galaxies do not have SDSS photometric measurements unlike the former and so we utilise the ESO DSS (Digital Sky Survey). R-band images were downloaded from the ESO DSS Online Archive222http://archive.eso.org/dss/dss. Photometry was then carried out using the same apertures as for the GALEX and WISE images, allowing for the calculation of consistent, aperture-matched NUV$-r$ colours.

2.3 Herschel observations and data reduction

Twenty-four of the galaxies within the xCOLD GASS survey were observed with the PACS spectrometer (Poglitsch et al., 2010) onboard Herschel (Pilbratt et al., 2010) as part of the programme OT2_asainton_1, P.I. A. Saintonge. The seven DGS galaxies were observed as part of a guaranteed time key programme.

The PACS array is composed of 5x5 square spatial pixels each of side 9.4”, covering a total field of view of 47”. The observations were carried out in line-spectroscopy mode. The data were reduced using the Herschel Interactive Processing Environment (HIPE) (Ott, 2010). We used standard scripts of the PACS spectrometer pipeline to reduce the data from the raw product to its level 2 processed form. From this, line fitting was done using IDL scripts, written by the first author, employing routines within the IDL Astronomy User’s Library (http://idlastro.gsfc.nasa.gov/). The lines are well fitted with a single gaussian using the IDL fitting routine GAUSSFIT. The signal from each spaxel of the PACS array is fitted with a second order polynomial with the addition of a Gaussian for the baseline and emission line respectively. From this we calculate the S/N in each spaxel. To calculate the total flux for each individual galaxy, we add up the integrals of the fitted Gaussians from all the spaxels which have a S/N $>$ 3. For total uncertainties we assume a 30% error of the total [C ii] line flux, due to calibration errors (Poglitsch et al., 2010). With Herschel all 24 xCOLD GASS galaxies show a clear detection in [C ii] with a S/N $>$ 3, while all seven DGS galaxies also show clear detections.

We want to directly compare the [C ii] line luminosities with the CO(1-0) luminosity and to do this we convolve the PACS spaxel array with the IRAM CO(1-0) beam, detailed in Section 2.4. We approximate the 2-D IRAM beam point spread function as a Gaussian, with a FWHM of 21.4”, where the peak is normalised to one333The IRAM beam width and beam efficiencies used in these calculation can be found at http://www.iram.es/IRAMES/mainWiki/Iram30mEfficiencies.. We ensure the peak of this Gaussian is set to one to give maximum weighting to the central spaxel. We do this for all the xCOLD GASS objects. For the DGS objects we do different aperture corrections depending on the size of the beam PSF used for each individual CO(1-0) measurement Rémy-Ruyer et al. (2014).

2.4 IRAM observations and data reduction

Observations were carried out at the IRAM 30m Telescope, as part of the xCOLD GASS survey, using the Eight Mixer Receiver (EMIR) (Carter et al., 2012) to observe the CO(1$-$0) line, which has a rest frequency of 115.27 GHz. EMIR allows for observations on two sidebands with 8GHz bandwidth per sideband per polarisation. The second band was tuned to cover the redshifted CO(2-1)444The analysis of this paper will focus on the CO(1-0) data only. The CO(2-1) data will be presented in a future publication.. The observing strategy was identical to that described in Saintonge et al. (2011b) for the initial COLD GASS survey.

The xCOLD GASS data are reduced using the Continuum and Line Analysis Single-dish Software (CLASS555http://www.iram.fr/IRAMFR/GILDAS/doc/html/class-html/class.html) which is part of the GILDAS software. All scans are visually inspected and those with anomalous features, such as distorted baselines or increased noise due to poor atmospheric conditions or high water vapour levels, are discarded. The individual scans for a single galaxy are baseline-subtracted, using a first order fit, and then combined. This average spectrum is then binned to achieve a resolution of $\sim$20kms*-1* and the rms is obtained and recorded. The flux of the CO(1$-$0) line is measured by adding the signal within an appropriately defined frequency window; in the case of a detection this window is set by hand to match the observed line profile. For the null detection scenario the window is set to a width of 200kms*-1* or to the full width of the HI line. Uncertainties on the CO(1-) line flux is calculated as:

[TABLE]

where $\Delta$ wch, the width of each spectral channel, is equal to $\sim$21km s*-1* with mild variations due to differing redshifts. The rms noise per spectral channel is denoted by $\sigma$rms and W50CO is the FWHM of the CO(1$-$0) line.

As mentioned in Section 2.3, we convolved the PACS spaxel array with the IRAM CO(1-0) beam to measure the $L_{\mbox{[C\scriptsize II]}}$/$L_{\mbox{CO(1-0)}}$ ratio over the footprint of the IRAM beam (or for the different beam FWHMs for the CO(1-0) observations in the DGS sample). As the angular size of the xCOLD GASS galaxies is larger than the IRAM beam, small aperture corrections need to be applied to obtain total CO fluxes. We calculated these aperture corrections as in Saintonge et al. (2012), and apply them to obtain total molecular gas masses when investigating the molecular gas scaling relations, in Section 6.

We plot the IRAM and Herschel spectroscopic data for an example galaxy in Figure 1. Similar plots for all of our 24 galaxies can be found in Appendix A. For the DGS galaxies, we use the CO(1-0) fluxes compiled by Rémy-Ruyer et al. (2014). The measured [C ii] and CO(1-0) line luminosities, as well as key physical parameters are given for the xCOLD GASS objects in Table 1, and for the DGS galaxies in Table 2.

3 Observational Results: [C ii]/CO Scaling Relations

In Figure 2, we present the L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ scaling relations for the 23 galaxies from the xCOLD GASS sample and the seven galaxies from the Dwarf Galaxy Survey, against a range of parameters describing the global physical properties of the galaxies. Although there is evident scatter, a clear dependence of L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ on many of the parameters is observed. We use the Pearson correlation coefficient, $r$, to quantify the strength of the dependence; a more refined statistical approach is presented later in Section 4.

The L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ depends most strongly on parameters which a) describe the strength of the UV radiation field and b) describe the ability of the CO molecule to shield itself, via dust, from the UV radiation impinging on the surface of the the molecular regions deep inside the PDR. A clear dependence of L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ on the colour of the systems, parametrised by NUV-r photometry, is observed, as well as the specific star formation rate, distance from the main sequence, UV field hardness, and the gas-phase metallicity. The two former quantities are directly linked because NUV-r is a good proxy for sSFR since it relates a quantity tracing ongoing star formation activity (NUV magnitudes) and a quantity sensitive to the older stellar population (r-band magnitudes). Their differences arise as the sSFR takes into account internal dust attenuation while NUV-r is not dust corrected. Furthermore the gas-phase metallicity can be seen as a proxy for the total dust to gas mass fraction (via a metallicity dependent dust to gas ratio). It is this dust which shields the CO molecule from the UV radiation, with decreasing metal content the systems become more dust poor allowing the UV radiation to penetrate deeper into the molecular clouds.

Overall the strongest dependencies in Figure 2 are on quantities that are sensitive to the amount of UV radiation penetrating into the molecular clouds, photodissociating CO into its ionised form. Of the four parameters mentioned above the dependence on metallicity is strongest and will be further justified with the full multi-parameter Bayesian treatment in Section 4.

Conversely the L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ ratio does not depend strongly on parameters which describe the masses and structural properties of the galaxies. A low correlation is observed with the stellar mass of the systems, and with the morphology as measured by the stellar mass surface density, $\mu_{*}$. This implies that CO photodissociation is happening on the small scales of molecular clouds as the large scale global properties of galaxies have a low correlation; global properties are likely less important than local ones.

Finally, rather interestingly, the L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ ratio does not depend strongly on the extinction, calculated from Equation 5 which is a measure of dust emission versus emission from young stars. Importantly it is a tracer of gas density and metallicity; the clear dependance on metallicity, explained above, suggests that the gas density may not correlate strongly with L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$. We speculate that the reason for this is that on small scales these low mass, low metallicity galaxies are flocculent with very clumpy structures, hence an averaged density across the galaxy washes out average variations of the L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ ratio on small scales. An alternative explanation for the low correlation with extinction is because, in particular at low stellar mass, NUV and 22$\mu$m emission do not originate from the same HII regions (Galametz et al., 2010), and so we are not tracing the extinction of the same star-forming clouds. It is difficult to disentangle these two effects, assuming both play a role, or to rule either one out.

Our results qualitatively imply that metallicity, colour, sSFR, $\Delta$(MS) and/or hardness of the UV radiation field are responsible for variations in the L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ ratio as the latter three parameters are responsible for the total amount of UV radiation impinging onto the surface layer of molecular clouds. As these three parameters are correlated we may not need all three to statistically parametrise a prescription for L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ and hence the conversion function; a rigorous statistical approach, employing Bayesian methods, is shown below in Section 4.

3.1 Contamination of the [C ii] emission by non-PDR sources

The scaling relations in Figure 2 show integrated measurements across all phases of the ISM, however to accurately measure variations of the conversion function we must restrict our analysis to variations in the L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ ratio from molecular regions only. The CO(1-0) emission line only originates from molecular regions while [CII] emission originates from all phases of the ISM; therefore we must clean for the contaminant [C ii] emission arising from non-molecular regions.

Galaxy wide observations of the [C ii] fraction emerging from molecular regions are uncommon and unfortunately an empirical prescription for this quantity does not exist, even though it is becoming apparent that the fraction of [C ii] emission across galaxies can vary. Pineda et al. (2013) have studied this quantity across the Milky Way, as part of the GOT C*+* survey, and found the fraction to be 75%. Further observations have been done which accurately measure the fraction of [C ii] emerging from ionised gas regions, using the method of Oberst et al. (2006). This utilises the [CII]/[NII]205$\mu$m and [NII]122$\mu$m/[NII]205$\mu$m ratios. A study of NGC 891 by Hughes et al. (2015) used this method and found that between 15$\%$-65$\%$ of [CII] emission can arise from ionised regions. A further study of the star forming region BCLMP 302 in M33 by Mookerjea et al. (2011) found that 20$\%$-30$\%$ of [CII] emission arose from ionised regions while Kramer et al. (2013) find that the contribution from the HI gas rises to about 40% in the outskirts of M33. Therefore, to study variations in the conversion function, it is crucial to remove non-molecular [CII] emission from our analysis.

Recent work by Accurso et al. (2017) provides a prescription to predict the fraction of [CII] emission originating from molecular regions as a function of galaxy properties. A new multi-phase 3D radiative transfer interface was used through the coupling of starburst99 (Leitherer et al., 1999, 2010), a Stellar Spectrophotometric code, with mocassin (Ercolano et al., 2003; Ercolano et al., 2005), a photoionisation code, and 3d-pdr (Bisbas et al., 2012), an astrochemistry code. A detailed explanation of the work can be found in Accurso et al. (2017), but a summary follows.

Individual spherically symmetric star forming regions, with a central population of stars surrounded by gas, were simulated as follows. First, the stellar radiation density field, coming from the stellar population centrally located with the star forming regions, is created using sb99. From this output stellar spectrum, the luminosity, temperature and number of ionising photons of the source are calculated; these quantities are then used as in- put parameters for the 3D photoionisation code mocassin., which is used to simulate the ionised gas surrounding the stellar population. The important output of the mocassin code is the SED of the ionised gas, dust and stars emerging from the ionised region. This can be used to calculate the strength of the radiation field, G0, which impinges onto the PDR surface, by integrating the mocassin SED in the far-UV range between 912 and 2400Å. This value of G0 is used as an input into 3d-pdr which simulates the ionised and molecular gas. Constant hydrogen number density profiles were assumed in each ISM region and the dust to gas ratio scaled with the input metallicity.

To jump from these simulated star forming regions to the ISM of entire galaxies it was assumed that the physical conditions found in each cloud, for a given set of input parameters, can represent the average physical conditions found on galaxy-wide scales for galaxies with similar physical properties. Under this assumption a whole galaxy can be considered to be built up from a number of identical star forming regions, hence it provides a larger whole-galaxy model of the ISM. A summary of the input parameters was given in Table 2 of the aforementioned paper and they do consider the metallicity and sSFR regimes probed here.

Through a Bayesian Inference method, several analytical functions were provided which can be used for this work. Specifically, as we do not have HII region electron densities for the xCOLD GASS and DGS samples, we will use the prescription involving sSFR only, namely:

[TABLE]

. where $f_{\mbox{\scriptsize{[C\tiny II\scriptsize], mol}}}$ is the fraction of [CII] emission originating from molecular regions where at least 1% of gas is in the form of H2, and $\psi=$ log(sSFR). The error on this quantity is $\sigma_{f_{\mbox{\scriptsize{[C\tiny II\scriptsize], mol}}}}=$ 0.063, ( the rms uncertainity on the fitting). The physical motivation for equation 9 can be found in Accurso et al. (2017). We apply this to the 23 xCOLD GASS and seven DGS galaxies used in this work to estimate the relative fraction of [CII] emission arising from molecular regions only. These values are in the range of 79-48% as shown in Tables 1 and 2. In these tables we also provide the integrated CO(1-0) and uncorrected [CII] luminosities for readers who are trying to model galaxies and want to match the uncorrected results.

The simulations of Accurso et al. (2017) do not consider the diffuse neutral medium owing to its extremely low density, however they do have an atomic gas component which represents the CNM. The [C ii] emission originating from this atomic component will not be included in the $f_{\mbox{\scriptsize{[C\tiny II\scriptsize], mol}}}$ term, calculated from the simulations, and therefore Accurso et al. (2017) does consider a CNM like component when transitioning from ionised to molecular states; however this is not of interest to this paper as we only need the fraction of [C ii] emission originating from molecular regions.

3.2 Corrected [CII]/CO Scaling Relations

In Figure 3 we present the L${}_{\mbox{[C {\sc ii}]}}$/L${}_{\mbox{CO(1-0)}}$ scaling relations for molecular regions by using the estimated fractions of ionised carbon emission arising from this phase of the ISM. As can be seen the correction for contaminant ionised carbon emission does not effect the qualitative trends shown in Section 3. There are mild changes in the statistical trends seen as, in most cases, the measure of correlation decreases as expected as the L${}_{\mbox{[C {\sc ii}]}}$/L${}_{\mbox{CO(1-0)}}$ numerical values have now decreased. Our corrected scaling relation still imply that metallicity, colour, sSFR, $\Delta$(MS) and/or hardness of the UV radiation field are responsible for variations in the conversion function. Note that the Pearson correlation coefficients and best-fititng linear relations in Fig. 3 are only intended to help with the interpretation; only the data themselves enter in the Bayesian analysis described below.

4 Bayesian Inference

Following on from the observational scaling relations, shown in Figure 3, we now want a) to parametrise an analytic expression for the [C ii]/CO luminosity ratio from molecular regions as a function of galaxy parameters and b) to determine the minimum number of parameters needed to provide a statistically robust fit to our data. With eight galaxy parameters666We do not include sSFR in our fitting procedure as this is very similar to $\Delta$(MS) which itself can be used to include a redshift dependence applicable up to z$\sim$2.5, the highest redshift probed in the $\Delta$(MS) prescription of Whitaker et al. (2012). available we fit models with different number of free variables using all the possible combinations of parameters e.g we have ${{}^{8}C_{k}}$ number of models when we are fitting for k number of parameters. We perform a Bayesian interference method to find the best fit relations and to find the minimum number of variables needed to fit the data. Bayesian interference fitting methods have been successfully employed in several, wide-ranging, astrophysical scenarios from the derivation of the extinction law in the Perseus molecular cloud (Foster et al., 2013) and Type Ia supernova light curve analysis (Mandel et al., 2011) to the extragalactic Kennicutt-Schmidt relation (Shetty et al., 2013) and the formation and evolution of Interstellar Ice (Makrymallis & Viti, 2014). For a more in depth explanation of the Bayesian regression fitting method we refer the reader to Kelly (2007) and restrict ourselves here to the basic concepts.

The first step is to assume that the measurement uncertainty associated with each [CII]/CO observation, for each galaxy, ($y_{i}$ hereafter) is normally distributed. Therefore $y_{i}$ is a random variable distributed like:

[TABLE]

. where $\sigma_{y_{i}}$ is the measurement uncertainty associated with the observable y on the ith galaxy. As can be seen in Figure 3 all the observed scaling relations show evidence of either no correlation or a linear correlation between the log variables. We therefore use power law models in linear space as none of the above plots, in log space, show higher polynomial behaviour such that:

[TABLE]

. Where $\alpha$ and $\beta_{j}$ are free variables to be found. We can say that the probability of observing our data, given the true value of the observables and the measurement uncertainties is:

[TABLE]

.

The next assumption to make is that all our [CII]/CO observables are independent, e.g each $y_{i}$ and $x^{j}_{i}$ are independent, which is perfectly reasonable as the result from one galaxy will not affect the result from another. With this, and using the definition of independent probabilities, we can simply multiply all the individual probabilities together. Therefore this product of probabilities is our Likelihood, denoted L. By taking the log-likelihood, $\mathcal{L}$, the product returns back to a sum so:

[TABLE]

where N is the sample size. Maximising this log-likelihood for all of our models, at fixed number of free parameters, will provide us with the best fitting models for a given number of free parameters.

4.1 Model comparison and sampling methods

We aim to maximise the likelihood for our fits where the number of degrees of freedom varies. To compare likelihoods from models with different numbers of free parameters we use two different methodologies. Firstly we employ the Akaike Information Criterion (AIC) (Akaike, 1981):

[TABLE]

. where p is the number of free parameters and N is the sample size and the preferred model is that which minimises AIC. We also use a second method by employing the Bayesian Information Criterion (BIC) (Schwarz, 1978):

[TABLE]

where p and N are defined as above. We compare the results of both these two methodologies to check for agreement and to ensure our results are not biased depending on what information criterion is employed.

A direct derivation/solution for the parameters which maximise the likelihood function, Equation 13, is computationally expensive and so, to efficiently and effectively sample the full parameter space, we use the well tested Python implementation of the affine-invariant ensemble sampler for Markov Chan Monte Carlo (MCMC) emcee 777An example of the code can be found at http://dan.iel.fm/emcee/current/ (Goodman & Weare, 2010).

4.2 Statistical results

In total we run $\sum\limits_{i=1}^{8}$ ${{}^{8}C_{i}}$ models with the results presented below. By maximising the likelihood for each number of free parameters (from one to eight), and then comparing models with different sample sizes, and using the BIC and AIC, we find that only two parameters are necessary to fully explain the trends seen in Figure 3. We show in Fig 4 variations of the BIC and AIC with different number of free parameters. Both the AIC and BIC method agree that only two free parameters are necessary to provide a statistically robust fit to the data and so our choice of information criterion is extraneous.

For the two parameter case both the AIC and BIC retrieve metallicity and $\Delta$(MS) as the two galaxy parameters needed to provide a good fit to the data (the two parameter model which has the highest likelihood value in both cases). The best fitting analytical prescription is:

[TABLE]

which has a regression correlation coefficient of 0.801 meaning that these two parameters alone account for 80.1% of the correlation, with an error on the predicted $\log\left(\frac{\mbox{[C\scriptsize II\normalsize]}}{\mbox{CO}}\right)$ value of 0.223. As the number of free parameters increases (from one to eight) the log posterior and log likelihoods also decrease; hence doing a likelihood ratio analysis would have been insufficient to recover our best fitting two parameter relation as this would have returned the eight parameter case as the best relation as these have the smallest posteriors. We stress that our aim was to find a relation which can explain the correlation with as few observables as possible. Only by employing the AIC and BIC, as above, were we able to find our best fitting relation with only two parameters.

We also fit solely to the xCOLD GASS objects to demonstrate that our result is secure even when including the DGS galaxies. For the xCOLD GASS objects we find:

[TABLE]

which differs from 16 by no more than 25%, well within the observational scatter.

5 Radiative Transfer modelling - connecting [CII]/CO and $\alpha_{CO}$

The objective in obtaining the above scaling relations is to be able to derive a parametrisation for the conversion function. There are two main parameters responsible for variations in the L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ fraction, namely metallicity and $\Delta$(MS), and so it is necessary to understand how to go from a L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ parametrisation to one for \alpha$${}_{\mbox{CO}}.

To better understand the astrochemical reactions involved in the photodissociation of CO we inspect the chemical reaction database employed in the PDR code 3d-pdr (Bisbas et al., 2012). The main reactions involved in photodissociation of CO arise from the interaction of molecular species with cosmic rays and UV photons. The dominant reactions888By this we mean the reactions with the highest reaction rates which provide the main routes for creation and destruction of molecular species involving UV photons and CO are:

[TABLE]

[TABLE]

indicating that, whenever CO interacts with UV photons, neutral carbon is formed first as an intermediate species, but then quickly gets ionised as the reaction rates are of the same order999The reaction rates for these two are 2.0x10*-10* s*-1* and 3.0x10*-10* s*-1* respectively.; hence only ionised carbon forms when CO reacts with UV photons. Moreover the reaction of CO with cosmic rays (denoted CRP) happens either via the formation of a free radical, namely He*+*, as follows,

[TABLE]

[TABLE]

[TABLE]

or directly as,

[TABLE]

These show how, via cosmic rays, the photodissociation of CO can lead to the formation of neutral and ionised carbon species. Recent work by Bisbas et al. (2015) has shown that ionised carbon is the main product formed through the photodissociation of CO, but knowledge of how much neutral carbon forms is crucial to accurately constrain a prescription of $\alpha_{\mbox{CO}}$ as we only have observed CO and [CII], lacking CI observations.

To describe the formation of neutral carbon we rely on the multi-phase ISM numerical simulations performed in Accurso et al. (2017). We take their simulated low-redshift cloud results and here plot in Figure 5 the variations of the molecular region L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ against $\alpha_{CO}$ for the 2160 simulated clouds with colours indicating varying metallicities and $\Delta$(MS) of the clouds. Accurso et al. (2017) do not discuss $\alpha_{CO}$, hence we calculate this here. For each simulated cloud we calculate the total molecular hydrogen mass and total CO(1-0) luminosity and hence can accurately obtain an $\alpha_{CO}$ value. We find that these two parameters, along with L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$, account for 98.4% of the total correlation with $\alpha_{CO}$. We therefore perform a three dimensional linear fit to the above data, using the Python SciPy routine curve-fit, and find

[TABLE]

. The errors on the estimated parameters are less than 0.01 dex and hence are negligible, owing to the high level of correlation retrieved from the fitting, and therefore can be ignored. We emphasise that this fit was done to the simulated data, hence the high degree of certainty in the $\alpha_{CO}$ calculation.This has been normalised to the known Milky Way values for $\alpha_{\mbox{\scriptsize{CO}}}^{solar}$ = 4.36 M*⊙* (K km s${}^{-1})^{-1}$ (helium corrected) and (L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$)solar = 1400 from Stacey et al. (1991). We then combine this with Equation 16 to finally arrive at a prescription for the CO conversion function:

[TABLE]

where 12 +log(O/H) is the metallicity in PP04 units, with an error on the predicted $\log\alpha_{\mbox{CO}}$ estimate of $\pm 0.165$. Metallicity predominantly drives variations in $\log\alpha_{\mbox{CO}}$ with $\Delta$(MS) playing a secondary, though statistically important, role.

The result of this study is therefore a new CO-to-H2 conversion function which involves metallicity and $\Delta$(MS), and which is applicable applicable on galaxy-wide scales for star forming galaxies up to z$\sim$2.5 (the highest redshift constrained by Eq. 7). The dependence on metallicity is consistent with other conversion functions as will be shown further in Section 5.3. The secondary dependence of $\alpha_{CO}$ on $\Delta$(MS), or more generally speaking on a parameter relating to the strength of the UV radiation field, does not come as a surprise as it has been previously observed or predicted by simulations and theoretical models (e.g. Israel, 1997; Wolfire et al., 2010; Sandstrom et al., 2012; Clark & Glover, 2015). The interesting result here is that the dependence of $\alpha_{CO}$ and both metallicity and $\Delta$(MS) was not forced, it came out naturally from the data through the Bayesian analysis technique.

5.1 Caveat - Suitability across the main sequence

There are two different well-studied regimes where $\alpha_{CO}$ is know to show deviations from the Galactic value. The first is the situation of low metallicities discussed above, when CO suffers from photodissociation from the UV radiation field. The second regime is the one where the gas is affected by dynamical effects (mostly galaxy mergers), making the ISM mostly molecular and changing gas temperatures and line widths. This has been well studied in the case of local ULIRGS, with Solomon et al. (1997) reporting that an $\alpha_{CO}$ value of 1.0 M*⊙* (K km s${}^{-1})^{-1}$ is more appropriate than the Galactic value. There will be galaxies which “suffer” from both effects, and therefore the two should ideally be studied simultaneously. While this has been done in a small number of studies (see Bolatto et al., 2013, and references thereif), it is unfortunately not possibly here as we do not have a sample of galaxies which simultaneously have low metallicities and dense ISMs, and the appropriate CO and [C ii] data products.

Therefore, the conversion factor given as Eq. 25 should only be applied to galaxies within the parameter space constrained here, i.e for those with $7.9<12+\log(\mbox{O/H})<8.8$ and $-0.8<\log\Delta(\mbox{MS})<1.3$. For galaxies with $12+\log(\mbox{O/H})>8.8$ we recommend using the value obtained when using $12+\log(\mbox{O/H})=8.8$, and for galaxies with $\log\Delta(\mbox{MS})<-0.8$ we recommend using the value obtained by setting $\log\Delta(\mbox{MS})=-0.8$. Because of a lack of data for galaxies below the main sequence, this is the most conservative thing to do, as opposed to allowing the trends to continue into regimes where it was not constrained. The conversion function will level out to a constant value, roughly equal to the galactic conversion factor value, consistently with other studies of early-type passive galaxies (Davis et al., 2014).

Similarly, Eq. 25 should not be applied to galaxies well above the main sequence, where dynamical effects may be compressing the ISM, leading to higher gas densities and dust temperatures. For galaxies which show evidence of being in this regime (e.g. Saintonge et al. (2012) use the selection criteria of L${}_{\mbox{FIR}}$/L${}_{\odot}>$ 11.0 and S60$\mu$m/S100$\mu$m $>$ 0.5), we recommend the use of a lower value of $\alpha_{CO}=1$ M*⊙* (K km s${}^{-1})^{-1}$. A major consequence of this is a discontinuity of $\alpha_{CO}$; when moving up through the main sequence $\alpha_{CO}$ increases and suddenly decreases as $\log\Delta(\mbox{MS})>1.3$. This has also been recently suggested by Sargent et al. (2014) who find that $\alpha_{CO}$ increases with increasing sSFR and then suddenly decreases when entering the starburst regime. We speculate this is because ULIRGs and starburst galaxies have different physical environments to main sequence galaxies; while their sSFRs are much higher their densities are also significantly higher, hence CO is more easily shielded in these denser environments, leading to a rapid decline in $\alpha_{CO}$.

5.2 Caveat - SDSS Metallicities

Metallicity gradients of different amplitudes have been observed in galaxies similar to those studied here (Magrini et al., 2016; Tissera et al., 2016; Wuyts et al., 2016; Ho et al., 2015) meaning that the SDSS fiber spectroscopy for our galaxies, which only probes the central 3”, may not provide reliable metallicity measurements to use in the analysis of our [CII]/CO data. Star formation histories for these low mass objects are known to be bursty while ongoing star formation activity is known to be inhomogeneous with activity unevenly scattered across the galaxy (Guo et al., 2015; Domínguez et al., 2015). Hence, to accurately measure metallicity gradients we would need IFU data to get accurate spatial sampling out to the extended regions of our galaxies without a priori information of the location of the HII regions. IFU spectroscopy would be essential to determine the accurate light-weighted, integrated metallicity measurements over the area of the galaxies probed by the IRAM and Herschel observations. Obtaining IFU spectroscopy of our targets may reduce the scatter found in the metallicity dependent scaling relations and may further refine our conversion function.

5.3 Comparison with previous studies

Studies concerning the CO-to-H2 conversion function have been ongoing for decades with a plethora of varying prescriptions; some of the earliest work was carried out by Solomon et al. (1987), Tielens & Hollenbach (1985) and van Dishoeck & Black (1988) and so, to interpret our prescription, we compare it with a selection of other calibration derived using different methodologies. We plot our two parameter conversion function in Figure 6, where colour in the multi-coloured filled surface denotes varying $\Delta(\mbox{MS})$. As can be seen, metallicity predominantly drives variations in our conversion function with $\Delta(\mbox{MS})$ playing a minor, though statistically important, role. Some well-established $\alpha_{CO}$ prescriptions are also shown in Figure 6.

We first compare our prescription to the single variable metallicity dependent conversion functions of Schruba et al. (2012) and Genzel et al. (2012). These were derived by assuming a constant molecular gas depletion time and then using an inverse star formation law. Our prescription is approximately consistent with both around the high metallicity range, but then predicts different conversion factors for 12 + log(O/H)$\lesssim$8.25. The divergence between the two aforementioned prescriptions is because in addition to using different galaxy samples as calibrators, they assumed different values for the molecular gas depletion time (1.8 and 1.0 Gyrs respectively). While diverging at the low metallicity end, these metallicity-dependent prescriptions would be adequate in the metallicity range covered by xCOLD GASS sample.

The novelty of our approach means that we are able to attribute the scatter found in these previous metallicity-only dependent relations to the position off the main sequence. Those systems with stronger radiation fields (due to higher $\Delta(\mbox{MS})$ values), at fixed metallicity, will have higher conversion factors. Also galaxies which are metal poor, at a fixed position on the main sequence, will have higher conversion factors than their metal rich counterparts.

We next compare our prescription to another metallicity dependent conversion function based on the PDR modelling of Wolfire et al. (2010). The lowest metallicity employed in their grid of models was 12 + log(O/H) $=$ 8.38. The trend in both conversion functions are in broad agreement for metallicities above 12 + log(O/H) $=$ 8.38 to solar, however they are offset from one another with our prescription predicting higher conversion factors. Their grid of models allowed for a varying ionisation field (for which the sSFR/$\Delta$(MS) contributes towards) and varying metallicities, but the parametrisation only involved metallicity.

We now move onto a comparison with Wilson (1995) and Arimoto et al. (1996) who both measured virialised masses of molecular clouds to infer the dependence of their conversion function with metallicity. The resulting prescription is also in broad agreement with the trends seen in the Leroy et al. (2011) sample, who used a dust-to-gas ratio method. As can be seen in Figure 6, the trend in both conversion functions are in broad agreement for metallicities above 12 + log(O/H) $=$ 8.25 to solar. However, they are offset from one another with our prescription predicting conversion factors two to three times higher. We speculate that this is due to the integrated nature of our observations versus the highly resolved, cloud scale-resolution, of Wilson (1995) and Arimoto et al. (1996). On galaxy-wide scales, it is possible to see more diffuse H2 molecular gas which is traced by [CII] and CI, as opposed to CO, hence the need for higher conversion factors to account for this.

Overall our prescription agrees well with others in the literature that have used integrated observations and assume a constant depletion time. Our prescription does not agree as well with those that have used a) a dust to gas ratio method, b) numerical modelling c) virialised gas mass estimates with observational resolution down to cloud scales. This is because our prescription accounts for the diffuse H2 gas, not found in individual clouds and GMCs, which is better traced by ionised and neutral carbon versus its molecular counterpart.

6 Bridging low and high redshift molecular gas studies

We apply our new prescription to the full xCOLD GASS sample (including the new low mass galaxy sample and the higher mass objects presented in Saintonge et al. (2011a)) and the PHIBSS1 sample (Tacconi et al., 2013). Star formation rates, stellar masses, CO luminosities and redshifts for the PHIBSS1 sample are taken from Tacconi et al. (2013). Metallicities are estimated using the Fundamental Mass-Metallicity Relation (FMR) from Mannucci et al. (2010). The aim of this is to investigate the consequences of our conversion function on gas scaling relations and how these change when going from a standard constant conversion factor to our prescription presented in Equation 25. We do not include the DGS objects in this analysis as we are only interested in the statistical trends of the molecular gas scaling laws; we need a complete sample of galaxies at different redshifts and so use the xCOLD GASS and PHIBSS1 samples. We leave a more thorough and statistically robust treatment of the scaling relations for a future publication.

We first plot in Figure 7 a histogram showing the distribution of conversion factors for the high and low mass xCOLD GASS galaxies and the PHIBSS1 samples. For the original COLD GASS galaxies the distribution peaks slightly below the galactic value ($\sim$ 4.36 M*⊙* (K km s${}^{-1})^{-1}$) with other values clustered closely around this, hence we expect the low redshift molecular gas scaling relations to remain relatively unchanged for M*∗$>$1010*. However for the lower mass galaxies the conversion factors extend to high values. Furthermore, the predicted conversion factors for the high redshift galaxies are at most a factor of two larger or smaller than the galactic conversion factor value so we do not expect a major change in the gas scaling laws at high redshifts.

6.1 Molecular gas fractions up to z$\sim$2.0

We present scaling relations for the CO-traced molecular gas mass fraction as a function of stellar mass, specific star formation rate and $\Delta$(MS) (from left to right) for the full xCOLD GASS and PHIBSS1 samples in Figure 8. We aim to showcase our new prescription for $\alpha_{CO}$ and so, for now, only use the CO-detected galaxies in both samples. In a future publication a more thorough, quantitative, statistical treatment of the whole sample, including detections and non-detections, will be performed.

The top row in Figure 8 presents the results obtained when using a constant galactic conversion factor; this had been previously explored for the COLD GASS sample in Saintonge et al. (2011a), only for galaxies with stellar masses $>$ 1010 M*⊙* and so we are extending the sample to include galaxies with stellar masses $>$ 109 M*⊙*. This relation for the PHIBSS1 sample was first presented in Tacconi et al. (2013) and Genzel et al. (2015). The bottom row presents the results obtained when our prescription is used. We fit quadratic polynomials to the xCOLD GASS data to qualitatively show the differences in the trends; we also present the product moment correlation co-efficient and scatter for each fit, calculated using the using the Python SciPy routine linregress.

As can been seen, the correlation with stellar mass remains unchanged for galaxies with M*∗$>$1010* M*⊙* at low redshift as the standard galactic conversion factor is a good approximation to the conversion function predicted from this work. However, due to the emergence of lower metallicities in galaxies with stellar masses less than 1010 M*⊙* the trend changes. We start to see a flattening of the molecular gas mass fraction for stellar masses between 109 -1010 M*⊙* for the low redshift galaxies, which is in excellent agreement with the trends found in Grossi et al. (2016) and the star formation models of Krumholz et al. (2009, 2008) demonstrated in Lu et al. (2015). In the star formation models, this flattening occurs as depletion timescales are independent of stellar mass, something which we shall explore in Section 6.2. This flattening also validates part of the ideal gas regulator model (sSFR = constant; Lilly et al. (2013)) for which no or very little dependence of gas fractions on stellar mass is expected, as we observe up to M*∗$\sim$1010M⊙. This result is also of interest given the trend of increasing HI gas mass fraction with decreasing stellar mass, down to M∗$\sim$109M⊙* (Brown et al., 2015), implying that low mass galaxies are less efficient at converting their HI gas reservoirs into molecular gas. For the high redshift objects the trend of increasing molecular gas fraction with decreasing stellar mass is still evident; this is true across all redshifts for M*∗$>$1010* (Tacconi et al., 2013). The decreasing gas fractions as stellar mass increases, at all redshifts, is driven by the flattening of the SFR-M*∗* relation at high stellar masses.

The correlation with specific star formation rate strengthens when using the prescription presented here. A tight correlation between molecular gas and star formation is expected, confirming that the variable $\alpha_{CO}$ presented here outperforms a constant Galactic value. This is in agreement with the relations found in Saintonge et al. (2016) showing that star formation activity in a galaxy is controlled by the total available gas. Interestingly, the high redshift galaxies simply extend the trends found from the low redshift sample.

We also plot the correlation with $\Delta$(MS) and find similar results to the trends with specific star formation rate for local universe galaxies. This is expected as sSFR${}_{\mbox{MS}}$ is roughly constant in the local universe and therefore $\Delta$(MS)$\propto$sSFR at low redshift. The high redshift galaxies are offset from the trend seen in the low redshift sample as a direct consequence of the evolution of sSFR on the main sequence with redshift (e.g Karim et al. (2011)) and the increase in the gas supply rate.

6.2 Molecular gas depletion times up to z$\sim$2.0

We present scaling relations for the CO-traced molecular gas depletion times, where tdep = M${}_{H_{2}}$/SFR, as a function of stellar mass, specific star formation rate and $\Delta$(MS) for the full xCOLD GASS and PHIBSS1 samples in Figure 9. While a positive correlation between tdep and M*∗* is observed when applying a Galactic conversion factor, this trend becomes statistically insignificant with a correlation coefficient of 0.289 once our conversion function is used. The important consequence of this is that tdep (or equivalently, star formation efficiency) does not depend on stellar mass, as also seen in the PHIBSS1 sample. This is all in agreement with previous trends found in the low and high redshift universe (Leroy et al., 2013; Tacconi et al., 2013, 2017; Santini et al., 2014; Huang & Kauffmann, 2014; Sargent et al., 2014; Genzel et al., 2015). This important result validates part of the equilibrium model (Davé et al., 2012; Lilly et al., 2013) which states that gas depletion times are independent or have very little dependence on stellar mass.

The correlation with sSFR, albeit weaker, is in agreement with the relations found in Saintonge et al. (2011b), and confirms the conclusion from Bothwell et al. (2014) that it is robust against the choice of a specific metallicity-dependent $\alpha_{CO}$ conversion function. Our shallower trend is expected because our prescription is dependent on $\Delta$(MS), which is closely linked to sSFR. The trend observed here is in excellent agreement with that found in Hunt et al. (2015). Moreover the high redshift galaxies are offset as previously reported in Saintonge et al. (2012). The trend between tdep and sSFR is redshift independent once accounting for the redshift evolution of the main sequence (see bottom right panel of Fig. 9 and Genzel et al., 2015; Tacconi et al., 2017). This suggests that the process of star formation on the main sequence is driven by similar physical mechanisms across cosmic time, independently of stellar mass.

Finally, we plot in the left panel of Figure 10 CO(1-0) luminosity versus star formation rate for the whole xCOLD GASS sample. Low redshift galaxies with stellar mass $<$ 1010 M*⊙* have much lower CO luminosities per unit star formation than their higher mass counterparts. Is this due to a higher star formation efficiency in lower mass galaxies or due to the photodissociation of CO leading to higher conversion factors? To this end we plot M${}_{H_{2}}$ versus star formation rate, using our conversion function. The divergence away from the trend in the left plot is accounted for by the photodissociation of CO and low mass galaxies are, on average, as efficient as higher mass galaxies at forming stars.

It now becomes apparent why our trend for the depletion time versus sSFR, presented in Figure 9, is in agreement with that presented in Hunt et al. (2015); their conversion function was derived under the assumption that the linear relation between LCO and SFR observed in high mass, solar metallicity galaxies extends to the regime of low mass, low metallicity galaxies, with any deviations being attributed to CO photodissociation effects. The results of this study find this assumption to have been correct.

7 Summary & Conclusions

We present here results from xCOLD GASS, a legacy survey of CO(1-0) observations form the IRAM observatory, combined with Herschel observations of ionised carbon. We provide the first scaling relations for the integrated $L_{\mbox{[C\scriptsize II]}}$/$L_{\mbox{CO(1-0)}}$ ratio as a function of several galaxy parameters for a sample of low metallicity galaxies over 2 dex in $L_{\mbox{[C\scriptsize II]}}$/$L_{\mbox{CO(1-0)}}$. From this, we corrected for contaminant [CII] emission from non-molecular phases of the ISM and provide scaling relations for the molecular region $L_{\mbox{[C\scriptsize II]}}$/$L_{\mbox{CO(1-0)}}$ ratio.

We show that the integrated and molecular-only $L_{\mbox{[C\scriptsize II]}}$/$L_{\mbox{CO(1-0)}}$ ratio depends most strongly on parameters which describe the strength of the UV radiation field and the ability of the CO molecule to shield itself, via dust, from the UV radiation impinging on the surface of the the molecular regions deep inside the PDR. A clear dependence on the colour (NUV-r) of the galaxies, specific star formation rate, UV field hardness, gas-phase metallicity and $\Delta$(MS) is observed. On the other hand, the L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ ratio does not depend strongly on large scale parameters which describe the mass and structural properties of the galaxies.

Through a Bayesian analysis we establish that only two parameters, metallicity and $\Delta$(MS), are needed to robustly parametrise a prescription for variations in L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$ across our combined xCOLD GASS and DGS sample. We use our parametrisation of L${}_{\mbox{[C\scriptsize II\normalsize]}}$/L${}_{\mbox{CO(1-0)}}$, alongside radiative transfer modelling, to present a novel conversion function where:

[TABLE]

The prescription should only be applied to galaxies with $7.9<12+\log(\mbox{O/H})<8.8$ and $-0.8<\log\Delta(\mbox{MS})<1.3$ and which are not in the high-pressure “merger” regime where a lower $\alpha_{CO}$ value should be use. Our $\alpha_{CO}$ conversion function is consistent with previous metallicity-only prescriptions, however we are now able to attribute deviations from this relation to a second order dependence on the offset of a galaxy from the star-forming main sequence. The main interest of this new $\alpha_{CO}$ conversion function is that it is calibrated independently of any assumption on the molecular gas depletion timescale, making it perfectly suited for the study of gas and star formation scaling relations.

We apply the new $\alpha_{CO}$ to the full xCOLD GASS and PHIBSS1 samples and investigate such gas scaling relations. We find a complete flattening of the relation between molecular gas mass fraction and stellar mass as stellar mass decreases; this is expected based on the linearity of the M∗-SFR relation in this mass regime, and the close relation between molecular gas and star formation. As previously reported, there are systematic variations of the gas depletion timescale across the low and high redshift galaxy populations, although the dependence on sSFR is slightly weaker than previously reported. The depletion time, however, does not vary significantly with stellar mass. Instead, the low LCO/SFR ratios in low mass galaxies are entirely due to photodissociation of CO, meaning that on average lower mass galaxies are equally efficient at forming stars than their high mass counterparts.

Acknowledgments

GA would like to thank the UK Science and Technologies Facilities Council (STFC) for their support via a postgraduate Studentship. AS acknowledges the support of the Royal Society through the award of a University Research Fellowship and of a Research Grant. BC is the recipient of an Australian Research Council Future Fellowship (FT120100660). BC and LC acknowledge support from the Australian Research Council s Discovery Projects funding scheme (DP150101734). This work is based on observations carried out with the IRAM 30m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). We thank the staff of the telescope for their help in conducting the COLD GASS observations. PACS, aboard Herschel, has been developed by MPE (Germany); 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 BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain).

Appendix A Herschel PACS and IRAM beam sizes and reduced spectra

Below we present the available Herschel and IRAM data we have assembled for each of the 24 xCOLD GASS galaxies, similar to that in Figure 1.