Molecular gas in the Herschel-selected strongly lensed submillimeter galaxies at z ~ 2–4 as probed by multi-J CO lines

C. Yang; A. Omont; A. Beelen; Y. Gao; P. van der Werf; R. Gavazzi; Z.-Y. Zhang; R. Ivison; M. Lehnert; D. Liu; I. Oteo; E. González-Alfonso; H. Dannerbauer; P. Cox; M. Krips; R. Neri; D. Riechers; A. J. Baker; M. J. Michałowski; A. Cooray; I. Smail

doi:10.1051/0004-6361/201731391

Home

All issues

Volume 608 (December 2017)

A&A, 608 (2017) A144

Full HTML

Free Access

Issue		A&A Volume 608, December 2017


Article Number		A144
Number of page(s)		41
Section		Extragalactic astronomy
DOI		https://doi.org/10.1051/0004-6361/201731391
Published online		15 December 2017

A&A 608, A144 (2017)

Molecular gas in the Herschel-selected strongly lensed submillimeter galaxies at z ~ 2–4 as probed by multi-J CO lines^⋆,^⋆⋆,^⋆⋆⋆

C. Yang (杨辰涛)¹^,2^,3^,4^,5, A. Omont⁴^,5, A. Beelen², Y. Gao (高煜)¹, P. van der Werf⁶, R. Gavazzi⁴^,5, Z.-Y. Zhang (张智昱)⁷^,8, R. Ivison⁷^,8, M. Lehnert⁴^,5, D. Liu (刘岱钟)⁹, I. Oteo⁷^,8, E. González-Alfonso¹⁰, H. Dannerbauer¹¹^,12, P. Cox¹³, M. Krips¹⁴, R. Neri¹⁴, D. Riechers¹⁵, A. J. Baker¹⁶, M. J. Michałowski¹⁷^{, 7}, A. Cooray¹⁸ and I. Smail¹⁹

¹ Purple Mountain Observatory/Key Lab of Radio Astronomy, Chinese Academy of Sciences, 210008 Nanjing, PR China
e-mail: yangcht@pmo.ac.cn
² Institut d’Astrophysique Spatiale, CNRS, Univ. Paris-Sud, Université Paris-Saclay, Bât. 121, 91405 Orsay Cedex, France
³ Graduate University of the Chinese Academy of Sciences, 19A Yuquan Road, Shijingshan District, 10049 Beijing, PR China
⁴ CNRS, UMR 7095, Institut d’Astrophysique de Paris, 75014 Paris, France
⁵ UPMC Univ. Paris 06, UMR 7095, Institut d’Astrophysique de Paris, 75014 Paris, France
⁶ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
⁷ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
⁸ European Southern Observatory, Karl Schwarzschild Straße 2, 85748 Garching, Germany
⁹ Max Planck Institute for Astronomy, Konigstuhl 17, 69117 Heidel- berg, Germany
¹⁰ Universidad de Alcalá, Departamento de Física y Matemáticas, Campus Universitario, 28871 Alcalá de Henares, Madrid, Spain
¹¹ Instituto de Astrofísica de Canarias (IAC), 38205 La Laguna, Tenerife, Spain
¹² Universidad de La Laguna, Dpto. Astrofísica, 38206 La Laguna, Tenerife, Spain
¹³ Joint ALMA Observatory, 3107 Alonso de Córdova, Vitacura, Santiago, Chile
¹⁴ Institut de Radioastronomie Millimétrique (IRAM), 300 rue de la Piscine, 38406 Saint-Martin-d’ Hères, France
¹⁵ Astronomy Department, Cornell University, 220 Space Sciences Building, Ithaca, NY 14853, USA
¹⁶ Department of Physics and Astronomy, Rutgers, The State University of New Jersey, 136 Frelinghuysen Road, Piscataway, NJ 08854-8019, USA
¹⁷ Astronomical Observatory Institute, Faculty of Physics, Adam Mickiewicz University, ul. Słoneczna 36, 60-286 Poznań, Poland
¹⁸ Department of Physics and Astronomy, University of California, Irvine, CA 92697, USA
¹⁹ Centre for Extragalactic Astronomy, Durham University, Department of Physics, South Road, Durham DH1 3LE, UK

Received: 18 June 2017
Accepted: 13 September 2017

Abstract

We present the IRAM-30 m observations of multiple-J CO (J_up mostly from 3 up to 8) and [C I](³P₂ → ³P₁) ([C I](2–1) hereafter) line emission in a sample of redshift ~2–4 submillimeter galaxies (SMGs). These SMGs are selected among the brightest-lensed galaxies discovered in the Herschel-Astrophysical Terahertz Large Area Survey (H-ATLAS). Forty-seven CO lines and 7 [C I](2–1) lines have been detected in 15 lensed SMGs. A non-negligible effect of differential lensing is found for the CO emission lines, which could have caused significant underestimations of the linewidths, and hence of the dynamical masses. The CO spectral line energy distributions (SLEDs), peaking around J_up ~ 5–7, are found to be similar to those of the local starburst-dominated ultra-luminous infrared galaxies and of the previously studied SMGs. After correcting for lensing amplification, we derived the global properties of the bulk of molecular gas in the SMGs using non-LTE radiative transfer modelling, such as the molecular gas density n_H₂ ~ 10^2.5–10^4.1 cm^-3 and the kinetic temperature T_k ~ 20–750 K. The gas thermal pressure P_th ranging from~10⁵ K cm^-3 to 10⁶ K cm^-3 is found to be correlated with star formation efficiency. Further decomposing the CO SLEDs into two excitation components, we find a low-excitation component with n_H₂ ~ 10^2.8–10^4.6 cm^-3 and T_k ~ 20–30 K, which is less correlated with star formation, and a high-excitation one (n_H₂ ~ 10^2.7–10^4.2 cm^-3, T_k ~ 60–400 K) which is tightly related to the on-going star-forming activity. Additionally, tight linear correlations between the far-infrared and CO line luminosities have been confirmed for the J_up ≥ 5 CO lines of these SMGs, implying that these CO lines are good tracers of star formation. The [C I](2–1) lines follow the tight linear correlation between the luminosities of the [C I](2–1) and the CO(1–0) line found in local starbursts, indicating that [C I] lines could serve as good total molecular gas mass tracers for high-redshift SMGs as well. The total mass of the molecular gas reservoir, (1–30) × 10¹⁰M_⊙, derived based on the CO(3–2) fluxes and α_CO(1–0) = 0.8 M_⊙ ( K km s^-1 pc²)^-1, suggests a typical molecular gas depletion time t_dep ~ 20–100 Myr and a gas to dust mass ratio δ_GDR ~ 30–100 with ~20%–60% uncertainty for the SMGs. The ratio between CO line luminosity and the dust mass L′_CO/M_dust appears to be slowly increasing with redshift for high-redshift SMGs, which need to be further confirmed by a more complete SMG sample at various redshifts. Finally, through comparing the linewidth of CO and H₂O lines, we find that they agree well in almost all our SMGs, confirming that the emitting regions of the CO and H₂O lines are co-spatially located.

Key words: galaxies: high-redshift / galaxies: ISM / infrared: galaxies / submillimeter: galaxies / radio lines: ISM / ISM: molecules

^⋆

Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.

^⋆⋆

Based on observations carried out under project number 076-16, 196-15 and 079-15 (PI: C. Yang); 252-11 and 124-11 (PI: P. van de Werf) with the IRAM-30 m Telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain).

^⋆⋆⋆

The reduced spectra (FITS files) are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/608/A144

© ESO, 2017

1. Introduction

The strongest starbursts throughout the star formation history of our Universe are the high-redshift hyper- and ultra-luminous infrared galaxies (HyLIRGs and ULIRGs). With infrared luminosities integrated over 8–1000 μm L_IR ≥ 10¹³ L_⊙ and 10¹³L_⊙ > L_IR ≥ 10¹² L_⊙, respectively, and star formation rate (SFR) around 1000 M_⊙ yr^-1, they approach the limit of maximum starbursts (Barger et al. 2014). Despite having comparable or slightly higher luminosities than the local ULIRGs (Tacconi et al. 2010), these submillimeter (submm) bright galaxies (SMGs, see reviews of Blain et al. 2002; Casey et al. 2014) are different, being more extended and unlike nuclear starbursts of local ULIRGS. This population of dusty starburst galaxies was first discovered in the submm band using Submillimeter Common-User Bolometer Array (SCUBA, Holland et al. 1999) on the James Clerk Maxwell Telescope (Barger et al. 1998; Hughes et al. 1998; Smail et al. 1997), and later the spectroscopy observations revealed a median redshift of ~2.5 (Chapman et al. 2005; Danielson et al. 2017). Their extremely intense star formation activity indicates that these “vigorous monsters” generating enormous energy at far-infrared (FIR) are in the critical phase of rapid stellar mass assembly. They are believed to be the progenitors of the most massive galaxies today (e.g. Simpson et al. 2014). Nevertheless, theoretical models of galaxy evolution have been challenged by the observed large number of high-redshift SMGs (e.g. Casey et al. 2014).

Since the initial discovery of SMGs at 850 μm with SCUBA at the end of the last century, Chapman et al. (2005) carefully studied the properties of this 850 μm-selected SMG population and concluded that those with S_{850 μm} > 1 mJy contribute a significant fraction to the cosmic star formation around z = 2–3, that is ≳10%. Several other works have also confirmed that SMGs play a key role in the cosmic star formation at high-redshift (e.g. Murphy et al. 2011; Magnelli et al. 2013; Swinbank et al. 2014; Michałowski et al. 2017). For the ULIRGs studied with a median redshift of 2.2, it can be >65% according to Le Floc’h et al. (2005) and Dunlop et al. (2017, see ALMA counts by Karim et al. 2013; Oteo et al. 2016; Aravena et al. 2016a; Dunlop et al. 2017, and references therein for updated SMG counts; and Casey et al. 2014, for redshift distributions of SMGs selected at 850–870 μ.

It is important to understand the extreme star-forming activity within SMGs through studying their molecular gas content which serves as the basic ingredient for star formation, especially those at the peak of the star formation history (i.e. z ~ 2–3, Madau & Dickinson 2014). Nevertheless, due to their great distances, the number of well-studied high-redshift SMGs with several CO transitions at different energy levels is limited (see reviews of Solomon & Vanden Bout 2005; Carilli & Walter 2013) and this is mostly achieved through strong gravitation lensing and/or in quasi-stellar objects (QSOs), including IRAS F10214+4724 (Ao et al. 2008), APM 08279+5455 (Weiß et al. 2007), Cloverleaf (Bradford et al. 2009), SMM J2135-0102 (Danielson et al. 2011), G15v2.779 (Cox et al. 2011) and in the weakly lensed SMG, HFLS3 at z = 6.34 (Riechers et al. 2013; Cooray et al. 2014). Our knowledge of the detailed physical and chemical properties and processes related to star formation within these high-redshift Hy/ULIRGs is still limited.

Tacconi et al. (2008) found that high-redshift SMGs have large reservoirs of molecular gas about 10^10–11 M_⊙ (see also Ivison et al. 2011; Riechers et al. 2011c). CO rotational lines are contributing a significant amount of cooling of the molecular gas. By measuring the multiple-J CO lines, we can constrain the kinetic temperature and the gas density of the emitting regions (e.g. Rangwala et al. 2011) using non-local thermodynamic equilibrium (non-LTE) models. From the observations of the aforementioned individual high-redshift galaxies, the variety of CO spectral line energy distribution (SLED) shows that multiple molecular gas components in terms of their different gas densities and kinetic temperatures are required to explain the entire CO SLEDs. The mid/high-J CO emission can be explained by a warm component with molecular gas volume density of 10³–10⁴ cm^-3 which is more closely related to the ongoing star formation, while there is also an extended cool component dominating the low-J CO (e.g. Ivison et al. 2010; Danielson et al. 2011). Recent works with Herschel SPIRE/FTS spectra of 167 local galaxies by Liu et al. (2015) and 121 local LIRGs by Lu et al. (2017) also favour the presence of multiple CO excitation components. Daddi et al. (2015) reached similar conclusions for z ~ 1.5 normal star-forming galaxies. The differences in the J_up > 6 part of the CO SLEDs reveal different excitation processes (e.g. Lu et al. 2017): in most cases, the CO emission is insignificant for J_up> 7 CO lines; in the few cases (≲10%) where L_IR is dominated by an active galactic nucleus (AGN) there is a substantial excess of CO emission in the J_up> 10 CO lines (van der Werf et al. 2010), likely associated with AGN heating of molecular gas; there could also be a small number of exceptional cases, like NGC 6240, where shock excitation dominates (Meijerink et al. 2013).

Thanks to the extra-galactic surveys at FIR and submm bands like the Herschel-Astrophysical Terahertz Large Area Survey (H-ATLAS, Eales et al. 2010), the Herschel Multi-tiered Extragalactic Survey (HerMES, Oliver et al. 2012) and South Pole Telescope (SPT) survey (Vieira et al. 2013), large and statistically significant samples of SMGs have been built. It was found that with a criterium of source flux at 500 μm, namely S_{500 μm} > 100 mJy (galaxy–galaxy) strongly lensed high-redshift SMGs can be efficiently selected (e.g. Negrello et al. 2007, 2010, 2017; Vieira et al. 2010; Wardlow et al. 2013; Nayyeri et al. 2016). The strong lensing effect not only boosts the sensitivity of observations but also improves the spatial resolution so that we can study the high-redshift galaxies in unprecedented detail (see e.g. Swinbank et al. 2010).

The spectroscopy redshifts (mostly determined from CO lines) have now been determined in more than 24 Herschel-selected, lensed H-ATLAS SMGs thanks to the combined use of various telescopes; for example, Herschel itself, using the SPIRE/FTS (George et al. 2013; Zhang et al., in prep.), CSO with Z-Spec (Scott et al. 2011; Lupu et al. 2012), APEX (Ivison et al., in prep.), IRAM/PdBI (Cox et al. 2011; Krips et al., in prep.), LMT, ALMA (Asboth et al. 2016) and especially the Zpectrometer on the GBT (Frayer et al. 2011; Harris et al. 2012, in prep.) and CARMA (Riechers et al. 2011b, and in prep.).

In their parallel work on strongly lensed SMGs (Vieira et al. 2010, 2013; Hezaveh et al. 2013; Spilker et al. 2016, and the references therein), the SPT group used a selection based on the 1.4 mm continuum flux density. The ALMA blind redshift survey of these 1.4 mm-selected SMGs shows a flat redshift distribution in the range z = 2–4, with a mean value of ⟨ z ⟩ = 3.5, being in contrast to the 850–870 μm SCUBA/LABOCA-selected sample (Weiß et al. 2013; Spilker et al. 2016). This can be explained by the different flux limits of the two samples, namely, the SPT-selected sources are intrinsically brighter than the classic 850–870 μm SCUBA/LABOCA-selected SMGs (Koprowski et al. 2014).

Efficient CO detection in lensed SMGs has significantly enlarged the sample size of multi-J CO detections, with the aim of allowing statistical studies. Thus, we present here our observations of multi-J CO emission lines in 16 H-ATLAS lensed SMGs at z ~ 2−4, for a better understanding of the physical conditions of the ISM in high-redshift SMGs on a statistical basis.

Although there is a large number of CO observations in high-redshift sources, only a few high-density tracers with high dipole, for example, HCN, have so far been detected, most of which in QSOs (e.g. Gao et al. 2007; Riechers et al. 2010), and even fewer detections in SMGs (e.g. Oteo et al. 2017). Submm H₂O lines, another dense gas tracer, have been reported in 12 H-ATLAS lensed SMGs (Omont et al. 2011; Omont et al. 2013, O13 hereafter Yang et al. 2016, Y16 hereafter) using IRAM NOrthern Extended Millimeter Array (NOEMA), and also in other galaxies (see the review by van Dishoeck et al. 2013). An open question is whether or not the submm H₂O emission lines trace similar regions as traced by mid/high-J CO and HCN. The difficulty of the comparison is coming from the currently limited high-resolution mapping of the submm H₂O lines. However, by comparing line profiles of unresolved observations of lensed SMGs, Y16 argue that the mid-J CO lines originate in similar conditions to the submm H₂O lines. This can be further tested by a larger sample from this work, and more directly, the high angular-resolution mapping of the emissions: see, for example, the cases of SDP 81 as probed by ALMA (ALMA Partnership 2015), NCv1.143 observed by NOEMA and of G09v1.97 through ALMA observations (Yang et al., in prep.).

In this paper, we study the physical properties of the molecular gas in a sample of 16 lensed SMGs at z ~ 2–4 by analysing their multiple-J CO emission lines. This paper is organised as follows: we describe our sample, the observations and data reduction in Sect. 2. The observed properties of the multi-J CO emission lines are presented in Sect. 3. The global properties of the SMGs together with the differential lensing effect is discussed in Sect. 4. A detailed discussion of the CO excitation is given in Sect. 5. Sect. 5.3 describes the discussion of molecular gas mass and star formation. We compare the emission lines of CO and submm H₂O in Sect. 5.4. Finally, we summarise our results in Sect. 6. A spatially-flat ΛCDM cosmology with H₀ = 67.8 ± 0.9 km s^-1 Mpc^-1, Ω_M = 0.308 ± 0.012 (Planck Collaboration XIII 2016) and Salpeter’s (1955) initial mass function (IMF) has been adopted throughout this paper.

2. Sample, observations, and data reduction

2.1. Selection of the lensed SMGs

Unlike the previously studied SMGs, our sample is drawn from shorter wavelengths using Herschel SPIRE photometric data at 250, 350, and 500 μm. In order to find the strongly lensed SMGs, all of our targets were selected from the H-ATLAS catalogue (Valiante et al. 2016) with a criterion of S_{100 μm} > 100 mJy based on the theoretical models of the submm source number counts (e.g. Negrello et al. 2010, 2017). Then, a Submillimeter Array (SMA) subsample was constructed based on the availability of previously spectroscopically confirmed redshifts obtained by CO observations (Bussmann et al. 2013, hereafter Bu13); it includes all high-redshift H-ATLAS sources with F_{500 μm} > 200 mJy in the GAMA and NGP fields (300 deg²). From SMA 880 μm images and the identification of the lens deflectors and their redshifts, Bu13 built lensing models for most of them.

Our sample was thus extracted from Bu13’s H-ATLAS-SMA sources with the initial goal of studying their H₂O emission lines (see Table 6 of Y16). It consists of 17 lensed SMGs with redshift from 1.6 to 4.2. We have detected submm H₂O emission lines in 16 sources observed with only one non-detection from the AGN-dominated source, G09v1.124 (O13; Y16, Table 2). However, for this CO follow-up observation, we dropped three sources among the H₂O-detected 16: SDP 11 due to its low redshift z < 2, NCv1.268 because of its broad linewidth that brings difficulties for line detection in a reasonable observing time, and G15v2.779 because it has already been well observed by Cox et al. (2011). Nevertheless, we included G15v2.779 in discussing the main results to have a better view of CO properties for the whole sample. Our CO sample of 14 sources (13 observed with the IRAM’s Eight Mixer Receiver, for example, EMIR, in this work plus G15v2.779 studied by Cox et al. 2011) is thus a good representative for the brightest high-redshift H-ATLAS lensed sources with F_{500 μm} > 200 mJy and at z> 2 (except SDP 81 with F_{500 μm} ~ 174 mJy). Besides these 14 sources, we also include two slightly less bright sources, G12v2.890 and G12v2.257, down to F_{500 μm} > 100 mJy. In the end, as listed in Table 3, the entire sample includes 16 lensed SMGs from redshift 2.2 to 4.2.

Table 1

Basic information on the CO rotational lines and [C I] ³P fine structure lines used in this paper.

Table 2

Observation log.

The lensing models for twelve of the SMGs are provided by Bu13 through SMA 880 μm continuum observations. Table 3 lists the magnification factors (μ₈₈₀) and inferred intrinsic properties of these galaxies together with their CO redshifts from previous blind CO redshift observations. After correcting for the magnification, their intrinsic infrared luminosities are ~4–20 × 10¹² L_⊙. Since the lensed nature of these SMGs and their submm selection may bias the sample, we will compare their properties with other SMG samples later from Sect. 3 to Sect. 5.3.

In this work, in order to explore the physical properties of the bulk of the molecular gas, we targeted the rotational emission lines of CO, mostly from J_up = 3 to 8 and up to 11 in a few cases. [C I](2–1) line is also observed “for free” together with CO(7–6) . Basic information such as the frequencies, upper-level energies, Einstein A coefficients and critical densities of the CO and [C I] lines are listed in Table 1. The targeted CO lines are selected based on their redshifted frequencies so that they could be observed in a reasonably good atmospheric window in EMIR bands. In total, we observed 55 CO lines, with 8 [C I](2–1) lines acquired simultaneously with CO(7–6) in 15 sources (Table 2).

Table 3

Previously observed properties of the entire sample.

Fig. 1

Distribution of the observed velocity-integrated CO line flux density versus the rotational quantum number J_up for each transition, i.e. CO SLEDs. Black dots with error bars are the velocity integrated flux densities from this work. Red dots are the data from other works: all CO(1–0) data are from Harris et al. (2012); CO(4–3) in G09v1.97 is from Riechers et al. (in prep.); CO(6–5) , CO(7–6) and CO(8–7) in SDP 17b are from Lupu et al. (2012); CO(8–7) and CO(10−9) in SDP 81 are from ALMA Partnership (2015); CO(3–2) in G12v2.30, CO(4–3)in NCv1.143 and CO(3–2) in NBv1.78 are from O13; CO(4–3) in NAv1.56 is from Oteo et al. (in prep.). For a comparison, we also plot the CO SLED of G15v2.779 (Cox et al. 2011). We mark an index number for each source in turquoise following Table 3 for the convenience of discussion.

2.2. Observation and data reduction

The observations were carried out from 2011 June 30th to 2012 March 13th, and from 2015 May 26th to 2016 February 22nd using the multi-band heterodyne receiver EMIR (Carter et al. 2012) on the IRAM-30 m telescope. Bands at 3 mm, 2 mm, 1.3 mm and 0.8 mm (corresponding to E090, E150, E230 and E330 receivers, respectively) were used for detecting multiple CO transitions. Each bandwidth covers a frequency range of 8 GHz. We selected the wide-band line multiple auto-correlator (WILMA) with a 2 MHz spectral resolution and the fast Fourier Transform Spectrometer with a 200 kHz resolution (FTS200) as back ends simultaneously during the observations. Given that the angular sizes of our sources are all less than 8′′, observations were performed in wobbler switching mode with a throw of 30′′. Bright planet/quasar calibrators including Mars, 0316+413, 0851+202, 1226+023, 1253-055, 1308+326 and 1354+195 were used for pointing and focusing. The pointing model was checked every two hours for each source using the pointing calibrators, while the focus was checked after sunrise and sunset. The data were calibrated using the standard dual method. The observations were performed in average weather conditions with τ_{225 GHz} ≲ 0.5 during 80% of the observing time.

Data reduction was performed using the GILDAS¹ packages CLASS and GREG. Each scan of the spectrum was inspected by eye and the bad data (up to 10%) were discarded. The baseline-removed spectra were co-added according to the weights derived from the noise level of each. We also note that due to the upgrade of the optical system of the IRAM-30 m telescope in November 2015, the telescope efficiency has been changed by small factors for lower band receivers (see the EMIR commissioning report by Marka & Kramer 2015², for details). All our sources are a factor of 3–7 smaller compared with the beamsize of IRAM-30 m at the observing frequencies, so that they can be treated as point sources. Accordingly, we apply the different point source conversion factors (in the range of 5.4–9.7 Jy/K depending on the optics and the frequency) that convert $T_{A}^{*}$ $\hbox{$T^*_\mathrm{A}$}$ in units of K into flux density in units of Jy for the spectra. A typical absolute flux calibration uncertainty of ~10% is also taken into account. We then fit the co-added spectra with Gaussian profiles using the Levenberg-Marquardt least-square minimisation code MPFIT (Markwardt 2009) for obtaining the velocity integrated line fluxes, linewidths (FHWM), and the line centroid positions.

Fig. 2

Observed CO(3–2) -normalised CO SLED (without lensing correction) of the H-ATLAS SMGs, in which both J_up = 1 and J_up = 3 CO data are available. The inset shows a zoom-in plot of the flux ratio of CO(1–0) /CO(3–2). The grey histogram shows the ratio distribution, while the grey line shows the probability density plot of the line ratio (considering the error). A mean ratio of I_CO(1–0) /I_CO(3–2) = 0.17 ± 0.05 has been found for our lensed SMGs. This is 1.3 ± 0.4 times smaller than that of the unlensed SMGs of Bo13. For comparison, we also plot the SLED of the Milky Way and the Antennae Galaxy.

3. Observation results

3.1. Observed CO line properties

We have detected 47 out of 55 J ≥ 2 CO and 7 out of 8 [C I](2–1) observed emission lines in 15 H-ATLAS lensed SMGs (signal to noise ratio S/N ≳ 3, see Table B.1). The observed spectra are

Fig. 3

Upper panel: linewidths with errors from three different samples, with probability distributions obtained by adaptive kernel density estimate (Silverman 1986): black symbols and line are from this work, orange symbols and dashed-dotted line are the J_up ≥ 2 CO linewidth distribution in unlensed SMGs (Bo13) and the green symbols and dashed line represent the linewidth from the J_up ≤ 2 CO lines of the lensed SPT sources (Aravena et al. 2016b). Our lensed sources with μ> 5 are indicated with open circles while the other sources are shown in filled circles. We note that although there is no lensing model for G12v2.43 and NAv1.144, it is suggested that their μ are likely to be ~10 (see Sect. 4.2 and Fig. 6). Thus, they are also marked with open circles. Lower panel: cumulative distribution of ⟨ ΔV_CO ⟩ for the three samples with the same colour code.

Table 4

Dynamical masses of the sample.

displayed in Fig. A.1 and the fluxes are also shown in the form of CO SLEDs in Fig. 1, indicated by black data points. Detected multi-J CO lines are bright with velocity-integrated flux densities ranging from 2 to 22 Jy km s^-1. To further compare the CO SLEDs, the CO(3–2) normalised CO SLEDs are plotted in Fig. 2 for all the H-ATLAS sources with CO(3–2) detections, overlaid with those of the Milky Way (Fixsen et al. 1999) and the Antennae Galaxy (Zhu et al. 2003). The CO SLEDs are mostly peaking from J_up = 5 to J_up = 8. The histogram of the flux ratio between CO(1–0) and CO(3–2) shows that the average I_CO(1–0) /I_CO(3–2) ratio is 0.17 ± 0.05, which is 1.3 ± 0.4 times smaller than that of the unlensed SMGs (Bothwell et al. 2013, hereafter Bo13). This is likely to be related to differential lensing, in that the magnification factor of CO(3–2) is larger than that of CO(1–0) due to the differences in their emitting sizes. The resulting ratio of I_CO(3–2) /I_CO(1–0) is thus larger in our lensed sources compared to the unlensed SMGs. We further discuss this in Sect. 4.3. Here we define the ratio between the lensing magnification factor of CO(3–2) (assumed to be equal to the magnification factor μ₈₈₀ derived from SMA 880 μm images) and CO(1–0) to be $\frac{μ_{CO (3 - 2)}}{μ_{CO (1 - 0)}} = \frac{μ_{880}}{μ_{CO (1 - 0)}} = 1.3 \pm 0.4.$ $\begin{equation} \label{eq:diff_lens} \frac{\mu_\mathrm{CO(3\text{--}2)}}{\mu_\mathrm{CO(1\text{--}0)}}=\frac{\mu_{880}}{\mu_\mathrm{CO(1\text{--}0)}}=1.3\pm0.4. \end{equation}$ (1)We correct for differential lensing for CO(1–0) data using this factor as described in Sect. 3.3.

One of the most important characteristics of the CO lines is its linewidth. The CO linewidth (FHWM) distribution of our lensed SMGs is displayed by the black solid line in the upper panel of Fig. 3 with the corresponding cumulative fraction shown in the lower panel. This curve shows that the linewidths are distributed between 208 and 830 km s^-1 (see Table 4 for the weighted average values of the linewidth). Around 50% of the sources have linewidths close to or smaller than 300 km s^-1. The median of the whole distribution is 333 km s^-1 and its average value 418 ± 216 km s^-1. Figure 3 also displays the linewidth distributions and the cumulative curves of two other samples of unlensed SMGs (orange dash-dotted lines) and lensed SPT-selected SMGs (green dashed lines) for comparison as discussed in Sect. 3.2.

Among our 16 sources, 12 of them show a single Gaussian CO line profile. SDP 81, NBv1.78 and G15v2.235 have double Gaussian CO line profiles. Although G09v1.97 might show a single Gaussian line profile, it is likely that there is a weak component in the blue wing, that we have confirmed by a higher sensitivity ALMA observation (Yang et al., in prep.). The high S/N PdBI spectrum of CO(4–3) line in NAv1.56 (Oteo, in prep.) also shows a line profile consisting of a narrow blue velocity component and a broad red component. However, due to the limited S/N, we can only identify the CO(5–4) line observed by EMIR with a single Gaussian profile.

The CO line profiles between different J_up levels within each source may vary, since their critical density and excitation temperature are different. However, by checking our CO spectral data as displayed in Fig. A.1, we find the differences between the line profiles (mostly by checking the linewidth) are insignificant given the current S/N. Their linewidths generally agree with each other within their uncertainties.

3.2. Comparing our sample to the general SMG population

If we wish to use our sample of lensed sources and the increased sensitivity allowed by magnification to infer general properties of the SMG population, it is important to investigate whether or not it is representative of this population and to recognise the possible biases introduced by lensing selection. For this purpose, we may compare it, especially for CO emission, with the sample of unlensed SMGs of the comprehensive CO study by Bo13. Thanks to early redshift determination, this sample of 32 SMGs initially detected at 850 μm was the object of a large program at IRAM/NOEMA detecting multiple low/mid-J CO lines. As discussed by the authors, although not completely free from possible biases, the sample appears to be a good representative of the whole SMG population. Compared to ours, its redshift distribution is similarly concentrated in the redshift range 2 to 3, with a similar extension up to ~3.5, but it also extends below 2 down to z ~ 1 in contrast to our sample. Both samples have very comparable distributions of their FIR luminosity L_FIR (a typical ratio between L_IR and L_FIR is 1.9; e.g. Dale et al. 2001), from a few 10¹² L_⊙ to just above 10¹³ L_⊙, with a mean value of 6.0 × 10¹² L_⊙ for the Bo13 sample and 8.3 × 10¹² L_⊙ for ours. As expected from the Herschel selection of our sample, its dust temperature T_d (Bu13) is slightly higher (⟨ T_d ⟩ = 37 K) than for typical samples of 850 μm-selected SMGs such as that of Bo13; but there is no obvious evidence of any bias in our lensed sample with respect to the whole Herschel SPIRE SMG population.

An important parameter is the extension radii of the dust emission at submm, which is believed to be comparable to that of high-J CO emission as discussed by Bo13 (note that the CO(1–0) line is expected to be more extended, see below). Values of this radius for our sources are reported in Table 3 as computed in Bu13 lens models. All values remain <~3 kpc, with a mean value of ~1.5 kpc. A similar distribution was found by Spilker et al. (2016) for a larger sample of similar strongly lensed sources found in the SPT survey. These authors have compared the intrinsic size distribution of the strongly lensed sources (including Bu13 ones) to a similar number of unlensed SMGs and found no significant differences.

Fig. 4

L_IR vs. $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ from local star-forming galaxies to high-redshift SMGs. The low-z data shown in grey including galaxies with 10⁹ ≤ L_IR ≤ 10¹² L_⊙ (only <20% of the local sources are ULIRGs) are from Liu et al. (2015) and Kamenetzky et al. (2016), with the typical error shown by the grey error-bars. The high-redshift SMG data in blue are from Carilli & Walter (2013, including HFLS3 from Riechers et al. 2013. The red data points represent the H-ATLAS SMGs from this work. Solid light blue lines are linear fits to the local galaxies, showing the average ratios of $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ /L_IR, with the ±2σ limits indicated by the dashed green lines. The insets show the histograms of the distribution of the ratio between $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ and L_IR for the three samples. It is clear from the correlation plots and the histograms that the high-redshift SMGs are above the low-redshift correlation for J_up = 3 and J_up = 4, with a significant smaller ratio of $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ /L_IR. Our H-ATLAS SMGs are located in the same region as other SMGs.

In contrast with these similarities of lensed and unlensed SMG samples, the CO linewidths of our lensed flux-limited sample appear anomalously low on average as quoted above. This is obvious from the comparison with the Bo13 sample: see Fig. 3 and the comparison of the distribution of the linewidth, the mean values (±1σ) are 418 ± 216 km s^-1 for our H-ATLAS flux-limited sample, 502 ± 249 km s^-1 for Bo13 sources with z ≥ 2, and 430 ± 140 km s^-1 for the SPT lensed SMG based on CO(1–0) and CO(2–1) observations by Aravena et al. (2016b). The median values of linewidth for the three samples are 333 km s^-1, 445 km s^-1 and 420 km s^-1, respectively, while the mode values are 264 km s^-1, 346 km s^-1 and 328 km s^-1, respectively. The range of the CO linewidths of our lensed SMGs are similar to those of the unlensed Bo13 sample, although the former has a concentration towards a narrower linewidth; more precisely, 50% of them have linewidths ≲333 km s^-1. In order to further compare these three samples, KS-tests were performed. The value of KS probability P_KS will be small if the two comparing data sets are significantly different. For the linewidth of our sample and the unlensed SMG sample, P_KS = 0.23 with a maximum deviation of 0.3; while for comparing our sample with the SPT lensed SMG sample, P_KS = 0.30 and the maximum deviation equals 0.3. These values of P_KS show that the differences among the samples are not statistically significant, indicating that they could arise from similar distributions. Nevertheless, the shapes of the probability distributions and the accumulative distributions of the linewidth for the three samples show some differences as displayed in Fig. 3. The difference between our lensed sample and the SPT one might be expected since the lensed SPT linewidths come from CO(1–0) and CO(2–1) observations which likely trace a larger velocity range of the gas, and thus tend to have larger linewidths compared with mid/high-J CO lines. However, linewidths of the Bo13 SMG sample are also from low/mid-J CO observations. The difference between this unlensed sample and our H-ATLAS flux-limited sample is rather likely coming from differential lensing, as discussed in the subsequent subsection. We note, nevertheless, that the percentage of double-peak CO profiles appears consistent (~25%) for our sources and those of Bo13.

3.3. Intrinsic CO emission properties

We derive the apparent line luminosities, for example, μL_line (in units of L_⊙ ) and $μ L_{line}^{'}$ $\hbox{$\mu L^{\prime}_\mathrm{line}$}$ (in units of K km s^-1 pc²), from the observed line flux densities using the classical formulae as given by Solomon et al. (1992): $L_{line} = 1.04 \times 10^{-3} I_{line} ν_{rest} (1 + z)^{-1} D_{L}^{2}$ $\hbox{$L_\mathrm{line}=1.04 \times 10 ^{-3} I_\mathrm{line} \nu_\mathrm{rest} (1+z)^{-1}D_\mathrm{L}^2$}$ and $L_{line}^{'} = 3.25 \times 10^{7} I_{line} ν_{obs}^{-2} (1 + z)^{-3} D_{L}^{2}$ $\hbox{$L^{\prime}_\mathrm{line}=3.25 \times 10 ^{7} I_\mathrm{line} \nu_\mathrm{obs}^{-2} (1+z)^{-3}D_\mathrm{L}^2$}$ . The resulting line luminosities are listed in Table B.1. The range of the apparent line luminosities is $μ L_{line}^{'}$ $\hbox{$\mu L^{\prime}_\mathrm{line}$}$ ~ 2–48 × 10¹⁰ K km s^-1 pc². After correcting the lensing magnification, the range of the intrinsic CO line luminosities is ~(1–60) × 10⁷ L_⊙ or ~(2–170) × 10⁹ K km s^-1 pc². As usual, the value of $L_{CO}^{'}$ $\hbox{$L^{\prime}_\mathrm{CO}$}$ decreases with increasing J_up of the CO lines. Besides CO, we have also derived the intrinsic luminosities of the [C I](2–1) line, observed together with CO(7–6) , to be ~(3–23) × 10⁷ L_⊙ or ~(2–13) × 10⁹ K km s^-1 pc².

In the following analysis, we have included multi-J CO data found in the literature for our sources, especially CO(1–0) from Harris et al. (2012), compensating for the absence of this line in our observations (see caption of Fig. 1). However, due to the differential lensing effect on the CO(1–0) data as discussed in Sect. 4.3, we only use these CO(1–0) fluxes for the CO line excitation modelling, after applying a factor of 1.3 ± 0.4 to correct the differences between the magnification factors of mid/high-J CO and that of CO(1–0) following Eq. (1) (as argued in Sect. 4.3, we assumed the magnification of mid/high-J CO lines is equal to μ₈₈₀, and we use μ as μ₈₈₀ hereafter if not specified).

After correcting for the lensing magnification, Fig. 4 shows the correlation between the intrinsic values of L_IR and $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ lines from J_up = 3 to J_up = 11, over-plotted on the local correlations (Liu et al. 2015, see also Greve et al. 2014; Kamenetzky et al. 2016; Lu et al. 2017). One should note that >80% of the local sources in Liu et al. (2015) are galaxies with L_IR ≤ 10¹² L_⊙, that is, luminous infrared galaxies (LIRGs) and normal star-forming galaxies. As found previously, most of these local sources can be found well within a tight linear correlation between L_IR and $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ for the mid-J and high-J CO lines, although for the low-J CO lines, the local ULIRGs seem to be lying above the correlation at a ≳2σ level, having larger L_IR / $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ ratios (e.g. Arp 220). As shown by the histograms of the L_IR / $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ ratios in Fig. 4, comparing with local galaxies (mostly populated by galaxies with L_IR = 10⁹–10¹² L_⊙), both our H-ATLAS SMGs and the previously studied SMGs are slightly above the correlation with larger L_IR / $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ ratios for J_up = 3 to J_up = 5 CO lines. In contrast, for the J_up ≥ 6 CO lines, both the local galaxies and the high-redshift SMGs with L_IR from 10⁹ L_⊙ to a few 10¹³ L_⊙ can be found within tight linear correlations. The H-ATLAS SMGs show no difference with other previously studied SMGs. Among the CO transitions, CO(7–6) has the tightest correlation across different galaxy populations (~0.17 dex), which agrees well with Lu et al. (2015). This again indicates that the dense warm gas traced by the J_up ≥ 6 CO lines is more tightly correlated with on-going active star formation (without considering AGN contamination to the excitation of CO), and CO(7–6) may be the most reliable star formation tracer among the CO lines.

We have also compared the CO line ratios in local ULIRGs with those in our lensed SMGs, by taking CO(5–4) and CO(6–5) for example. The ratios of $L_{CO (5 - 4)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(5\text{--}4)}$}$ / $L_{CO (6 - 5)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(6\text{--}5)}$}$ from the two sub-samples turn out to be similar within the uncertainties. Their mean values are 1.6 and 1.4 with the standard deviations of 0.35 and 0.37 for local ULIRGs and high-redshift lensed SMGs, respectively. This suggests that the differential lensing is unlikely introducing a large bias of choosing molecular gas with very different gas conditions.

4. Galactic properties and differential lensing

4.1. Molecular gas mass

One of the most commonly used methods to derive the mass of molecular gas in galaxies is to assume that it is proportional to the luminosity $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ through a conversion factor α_CO such as $M_{H_{2}} = α_{CO} L_{CO (1 - 0)}^{'}$ $\hbox{$M_\mathrm{H_2}=\alpha_\mathrm{CO}\lpcol10$}$ , where M_H₂ is the mass of molecular hydrogen and α_CO is the conversion factor to convert observed CO line luminosity to the molecular gas mass without helium correction (see Bolatto et al. 2013, for a review). Here we adopt a typical value of α_CO = 0.8 M_⊙ (K km s^-1 pc²)^-1 which is usually found in starbursts as observed in local ULIRGs (Downes & Solomon 1998). The total mass of molecular gas M_gas is then inferred by multiplying M_H₂ by the factor 1.36 to include helium. One should also note that at z = 2.1–4.2, the cosmic microwave background (CMB) temperature reaches ~8.5–14.2 K, which is non-negligible to the low-J CO lines. A typical underestimation of the CO(1–0) luminosity could be around 10%–25% if T_k = 50 K, and for the bulk of the molecular gas, which is normally colder than 50 K and is only bright in the low-J transitions, the CMB effect may be even more severe as pointed out by Zhang et al. (2016; see also da Cunha et al. 2013). Although far from being settled, recent observations of high-redshift SMGs favour α_CO being close to the value of local ULIRGs with large uncertainties (Ivison et al. 2011; Magdis et al. 2011; Messias et al. 2014; Spilker et al. 2015; Aravena et al. 2016b).

Table 5

Observationally derived physical properties of the H-ATLAS SMGs.

Half of our sources were observed in their CO(1–0) line with the Green Bank Telescope (GBT) by Harris et al. (2012). The corresponding apparent luminosities μ $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ (not corrected for lensing) are reported in Table 5. However, it is impossible to infer the total mass of molecular gas in the absence of a detailed lensing model including the extended part of CO(1–0) emission. We may nevertheless directly compare the CO(3–2) and CO(1–0) apparent luminosities μ $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ for the seven Harris’ sources for which we observed the CO(3–2) line (Table 5). The error-weighted mean ratio of the luminosity of CO(3–2) to CO(1–0) is 0.65 ± 0.19. This is marginally larger at about the 1σ level by a factor 1.3 ± 0.4 than the median brightness temperature ratios r₃₂/r₁₀ of 0.52 ± 0.09 reported for unlensed SMGs by Bo13, and 0.55 ± 0.05 reported by Ivison et al. (2011, as described in Eq. (1. This difference seems to suggest an effect of differential lensing, the more compact CO(3–2) emission being more magnified than the extended CO(1–0) emission (see Sect. 4.3 for a detail discussion).

However, the mass of molecular gas M_gas can be directly inferred from higher J_up CO lines, mostly CO(3–2) , as for cases of other high-redshift SMGs where CO(1–0) observations are lacking. Moreover, comparing with the CO(1–0)line, the CO(3–2) line tends to be less affected by differential lensing because its spatial distribution is closer to that of the submm dust emission upon which the lensing models are built. Therefore, by assuming that our lensed SMGs are similar to the unlensed high-redshift SMGs, the brightness temperature ratio r₃₂/r₁₀ = 0.52 ± 0.09 from Bo13 yields β_CO32 = 1.36 × 0.8/0.52 = 2.09 M_⊙ (K km s^-1 pc²)^-1 for the conversion factor defined as $M_{gas} = β_{CO 32} L_{CO (3 - 2)}^{'} .$ $\begin{equation} \label{eq:beta_co} M_\mathrm{gas}=\beta_\mathrm{CO32}\lpcol32. \end{equation}$ (2)The masses of molecular gas (including He) are thus derived and reported in Table 5³ These values for M_gas are in the same range, 10¹⁰–10¹¹M_⊙, as those derived for unlensed SMGs by Bo13. This is confirmed by the direct comparison of the distributions of $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ after lensing correction (Fig. 6). But one should keep in mind the accumulation of uncertainties about our M_gas estimates: to the usual uncertainty on α_CO or β_CO32, one should add that of the lensing model, especially in the absence of high-resolution CO imaging. The derived gas mass appears exceptionally high for G15v2.235, about three times larger than for any other source and twice more massive than for any unlensed SMG of Bo13. Either the magnification factor is larger than the low value, 1.8 ± 0.3, derived by Bu13, or this source is an exceptional galaxy.

These masses of gas may be compared with the mass of dust derived, for example, through the gas to dust mass ratio $δ_{GDR} = M_{gas} / M_{dust} .$ $\begin{equation} \label{eq:gdr} \delta_\mathrm{GDR}=M_\mathrm{gas}/M_\mathrm{dust}. \end{equation}$ (3)The dust masses were taken from Bu13. We recall that they are derived by performing a single component modified black body model with the Herschel SPIRE and SMA photometric fluxes, with mass absorption coefficient κ_dust interpolated from Draine (2003). The values of δ_GDR for our sample are given in Table 5. They range from 31 ± 14 to 100 ± 33 with a mean of 56 ± 28. Our value is generally in agreement with the mean value of δ_GDR = 75 ± 10 for ALESS high-redshift SMGs (Simpson et al. 2014; Swinbank et al. 2014) within 1σ level. This range is also similar to that of the local ULIRGs (Solomon et al. 1997).

4.2. Dynamical mass

In high-redshift SMGs, an important fraction of the baryonic mass is in the form of molecular gas, and the CO linewidth can serve as a good dynamical mass indicator with an assumption about the dynamical structure and extent of the system (e.g. Tacconi et al. 2006; Bouché et al. 2007). From the measured linewidth of the CO lines, we can in principle derive the dynamical mass within the half-light radius (r_half) by assuming that the lensed SMG can be treated as either a virialised system or a rotating disk with an inclination angle i.

If the system is virialised, the dynamical mass can be calculated following the approach of Bo13 as $M_{dyn, vir} = 1.56 \times 10^{6} {(\frac{σ}{km s -1})}^{2} (\frac{r}{kpc}) M_{⊙},$ $\begin{equation} \label{eq:Mdyn_vir} M_\mathrm{dyn,vir}=1.56\times10^{6} \left({\sigma \over \mathrm{\,\kms}}\right)^2 \left({r \over \mathrm{kpc}}\right) \,\msun, \end{equation}$ (4)in which the velocity dispersion $σ = Δ V_{CO} / (2 \sqrt{2 \ln 2})$ $\hbox{$\sigma = \Delta V_\mathrm{CO}/(2\sqrt{2\ln{2}})$}$ , r is the radius of the enclosed region for calculating the dynamical mass and ΔV_CO is the CO linewidth. If the system is a rotating disk, the dynamical mass can be derived from $M_{dyn, rot} = 2.32 \times 10^{5} {(\frac{v_{cir}}{km s -1})}^{2} (\frac{r}{kpc}) M_{⊙},$ $\begin{equation} \label{eq:Mdyn_rot} M_\mathrm{dyn,rot}=2.32\times10^{5} \left({v_\mathrm{cir} \over \mathrm{\,\kms}}\right)^2 \left({r \over \mathrm{kpc}}\right)\,\msun, \end{equation}$ (5)following Wang et al. (2013) and Venemans et al. (2016), in which v_cir is the circular velocity equal to 0.75ΔV_CO/sini, r is the disk radius and i is the inclination angle of the disk on the sky in a range from 0° to 90°. By assuming an inclined thin-disk geometry, we can derive the inclination angle from the minor to major axis ratio b/a of the rotating disk as i = cos^-1(b/a). When possible, we use the minor to major axis ratio of the source image which was derived in the lensing models of (Bu13). This yields the values of i reported in Table 4. Otherwise, we assume an average inclination angle of 55° as suggested by Wang et al. (2013). However, we also note that 1/sini can take large values for galactic disks seen close to face on. We can thus calculate the two estimates of the dynamical mass enclosed in the half-light radius r_half, for example, M_dyn,vir and M_dyn,rot from the CO linewidth and/or the minor to major axis ratios following Eqs. (4) and (5).

However, it is important to note that there is certainly a significant fraction of the SMG mass distributed outside r_half. Accordingly, the value of M_dyn,vir or M_dyn,rot should serve as a lower limit of the dynamical mass of the entire region where molecular gas resides. This is a fortiori true for the total mass of the galaxy including the extended diffuse, cool component beyond ~3 kpc. As suggested by the CO(1–0) observations of a sample of SMGs at z ~ 2.4, Ivison et al. (2011) find a typical size of ~7 kpc for the CO(1–0) , with a linewidth of ~563 km s^-1. Using CO(2–1) data, Hodge et al. (2012) also suggest a rotating disk of molecular gas with a radius of ~7 kpc in GN20 and a CO linewidth equal to 575 ± 100 km s^-1. JVLA and ATCA observations of CO(1–0) emission at high-redshift reported by several other works (e.g. Greve et al. 2003; Riechers et al. 2011c; Deane et al. 2013; Emonts et al. 2016; Dannerbauer et al. 2017) also support the existence of such an extended cold gas component. Both their linewidth and the size of the emitting region are larger than those of most of our sources, suggesting an underestimation for the dynamical mass for our sample. But even the mass of the starburst, FIR-emitting core is likely to be underestimated by these formulas for M_dyn. It is challenging to estimate the ratio between the total dynamical mass within the entire CO-emitting region and the one we calculated from r_half, although Tacconi et al. (2006) suggest a value about 5 for such a ratio.

As seen in Table 4, the range of M_dyn,vir, given by Eqs. (4) and (5), varies from ~10¹⁰ M_⊙ to 3 × 10¹¹ M_⊙ while values of M_dyn,rot are comparable but slightly smaller by a factor of ~1.4 on average. This is again ~3–10 times smaller than the value of M_dyn,rot given in Ivison et al. (2011) and Hodge et al. (2012). We compare these values with the derived masses of gas and discuss them in the following subsection.

4.3. Possible lensing biases

It is well known that differential lensing may be a serious problem for galaxy-galaxy strong-lensing studies of extended objects, especially for multi-line and continuum comparison (e.g. Serjeant 2012; Hezaveh et al. 2012). Although the problem should be dimmer for the compact cores (r ≲ 1–3 kpc) emitting the continuum and high-J CO lines in our sources, it needs consideration in case of complex caustics at this scale. It may become worse for low-J CO studies since they may involve more extended SMG components (Ivison et al. 2011). In addition, the flux-limited selection of our lensed sources may bias our sample towards the most compact objects.

Because different excitation levels of CO trace predominantly regions with different gas density and different temperature (see the critical densities and the energy levels of the CO lines in Table 1), the sizes of the emitting regions of each J CO line are expected to be somewhat different. This variation of the emitting region of each CO transition will certainly bring differences in the resulting parameters, such as the total magnification factor, and derived quantities such as the molecular gas mass and line ratios with respect to the intrinsic ratios of the unlensed galaxy. The differential lensing effect could arise in a complex way from the specific spatial configuration of the caustic line with respect to the background emission. However, a detailed modelling of complex effects of differential lensing is beyond the scope of the present study for two main reasons: the low resolution of our single-dish CO data and their limited range of J_up values, mostly from 3 to 8. It is expected that the regions emitting such lines will not differ very much with J_up for most sources and remain close to that of the observed 880 μm dust continuum. This is also consistent with the similarity that we find for the linewidths of the different CO transitions within each source. We will therefore neglect the effects of differential lensing in estimating the ratios of these different mid/high-J emissions. Of course, the validity of this assumption should be verified, taking into account the particularities of each source and its lensing, when high-angular-resolution images are available. However, we can perform a first verification in the only case of our sources, SDP 81, for which ALMA high-resolution CO images have been published, noting that it is one of our most extended sources. These images from the ALMA long baseline campaign observation (ALMA Partnership 2015) show the resolved structure of the dust and CO. From lens modelling, the studies of Dye et al. (2015) and Rybak et al. (2015a,b) show that the differences between the magnification factors of the CO(5–4) and CO(8–7) lines are within 1σ and close to that of dust emission.

On the other hand, the effects of differential lensing might be much more severe and non-negligible when comparing CO(1–0) with high-J CO lines. Previous high-angular-resolution imaging studies of SMGs show evidence that the cooler, low-density emitting regions of low-J CO lines are more extended than the dust continuum emission and that of the high-J CO lines (e.g. Tacconi et al. 2008; Bo13; Spilker et al. 2015; Casey et al. 2014). Especially for CO(1–0) , JVLA images of high-redshift SMGs reveal a significant extension, usually several times larger on average, compared to high-J CO (Ivison et al. 2011, see also Engel et al. 2010; Riechers et al. 2011a). This is probably also true for the [C I](1–0) line because its spatial distribution agrees well with that of CO(1–0) (e.g. Ikeda et al. 2002; Glover et al. 2015).

With such an extension, typically ~7 kpc, substantial differential lensing seems unavoidable for most strong lensing configurations, yielding a lower magnification for CO(1–0) compared with more compact mid/high-J CO emission. Indeed, such an effect has already been directly observed in at least two strongly lensed SMGs: SDP 81 (Rybak et al. 2015b) and SPT0538-50 (Spilker et al. 2015). As suggested in Hezaveh et al. (2012) such a difference will be moderate in the low-magnified system with small μ but can be non-negligible for the highly magnified systems where μ> 10, which is likely the case of G09v1.40, SDP 81, NCv1.143, NAv1.144 and NAv1.56 in our sample (Table 3). By assuming that μ for CO(3–2) is similar to the dust continuum, the ratio of the velocity integrated flux density between CO(1–0) (from Harris et al. 2012) and CO(3–2) for our sources has a mean value of 0.17 ± 0.05 (Fig. 2), that is, a factor 1.3 ± 0.4 lower than in other high-redshift unlensed SMGs (Eq. (1)). Such a difference in the ratio of CO(1–0) over CO(3–2) can be explained by differential lensing. Nevertheless, one should note that the difference is not at a very significant level. Because we have no high-resolution maps of the CO lines, it is beyond our ability to reconstruct the exact magnification factor for each emission line. Here we assumed that the magnification factors are the same for all the J_up ≥ 3 CO lines as that of dust emission, and we applied the magnification factor μ₈₈₀ derived from the SMA 880 μm images (Table 3). For the CO(1–0) line, we thus applied the factor μ₈₈₀/1.3, as derived above for the multi-J CO line excitation modelling as described in Sect. 5.

Fig. 5

A comparison between molecular gas mass and dynamical mass of our H-ATLAS lensed sample. The grey dashed lines indicate the ratio of M_gas /M_dyn,vir. Colours are coded according to the average CO linewidth. The size of the symbol represents the value of r_half of each source. There is a clear trend that sources with smaller linewidths have large ratios of M_gas /M_dyn,vir. The source index can be found in Fig. 1 and Table B.1.

Another important aspect of differential lensing is the possible distortion of the line profile. As shown in the case of the high-resolution and high-sensitivity CO spectrum of SDP 81, the line profiles show asymmetry features with a prominent red component accompanied by a weaker blue component (ALMA Partnership 2015). By reconstructing the source in the image plane, Swinbank et al. (2015) show that SDP 81 is a clumpy rotating disk and the red part of the disk is more magnified than the blue part, which causes the line-profile asymmetry. This might also happen to our other sources, especially for the case in which the caustic lines cross only part of the galaxy. It is not impossible that such effects might lead to underestimate the wings of some lines and thus explain at least part of the excess of narrow linewidths that we observed (see the case of G09v1.97 in Sect. 5.4 and Fig. A.1). Another cause of this excess could be a possible bias between the magnification and intrinsic source size (Spilker et al. 2016) which could perhaps bias against composite broad profiles of slightly extended sources in an early merger state. However, it seems that further observation and modelling of high-resolution CO images is needed to progress in completely explaining if this excess is real.

Fig. 6

Upper panel: μ $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ plotted against CO linewidth. Orange data points represent the unlensed SMGs from Bo13, while grey data points (indexed as in Table 2) are from this work, for which the filled circles show the sources with existing lensing model and the black open squares are the ones without lensing models. The solid turquoise line shows the fit proposed by Bo13 for the relation of $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ vs. linewidth as $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ ∝ $Δ V_{CO}^{2}$ $\hbox{$\Delta V_\mathrm{CO}^2$}$ . The dashed turquoise line indicates the positions for μ = 10 assuming this relation. As shown in the plot, the magnification factors of G12v2.43 (#5) and NAv1.177 (#9) are likely to be large. Lower panel: intrinsic $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ plotted against the CO linewidth. The light orange region shows the ±2σ range of the scatter derived from the SMGs of Bo13. Our sources generally agree with the correlation.

An effect of underestimating the CO linewidth would be to underestimate dynamical masses. In order to check this, Fig. 5 shows the plot of the relation of M_gas and M_dyn changes with CO linewidth and source size. For most broad-line sources, the values of the ratio M_gas /M_dyn are from 0.2 to 0.5, which appear possible, although the high value for G15v2.235, 2.4 (source #13 as shown in Figs. 5 and 6) seems to point out a problem with its lensing model. On the other hand, for all narrow-line sources, values of M_gas /M_dyn greater than 1, and even than 3 for most of them (which is equivalent to a 1.7–2 fold underestimation of the linewidth) point out a serious problem. The Spearman’s correlation coefficient between CO linewidth and M_gas /M_dyn is −0.83 with a p-value of 0.0016. The r_half value of each source is also indicated by the symbol size. Four of the smallest sources (#1, #2, #7 and #11) have high M_gas /M_dyn values, since the differential lensing could also potentially affect the estimation of the source size. However, we find a much weaker correlation between r_half and M_gas /M_dyn, suggesting that the impact of differential lensing on the source size is much weaker compared to that of the linewidth. It is possible that the sources with narrow linewidth having higher values of M_gas /M_dyn is partly due to differential lensing since the dynamical mass is proportional to $Δ V_{CO}^{2}$ $\hbox{$\Delta{V^2_\mathrm{CO}}$}$ , as identified in SDP 81 (Rybak et al. 2015b; Dye et al. 2015; Swinbank et al. 2015) and in GO9v1.97 (Yang et al., in prep.) and perhaps in SDP 17b whose CO(4–3) and H₂O profiles are clearly asymmetric (O13). Also, it seems however that at least part of this problem reflects the fact that Eqs. (4) and (5) might underestimate the dynamical mass by a large factor, likely to be up to 5, as quoted, for example, by Tacconi et al. (2006) for unlensed SMGs. It is however obvious that further observation and modelling of high-spatial-resolution CO images is needed to progress in completely explaining such problems.

Acknowledging the possible bias from the narrow-linewidth sources, after excluding the sources with ΔV_CO < 400 km s^-1, and also G15v2.235 as mentioned before, we derived an average value of M_gas /M_dyn,vir = 0.34 ± 0.10, in line with the SPT sources (Aravena et al. 2016b), other unlensed SMGs (Bo13) and empirical model predictions (Béthermin et al. 2015). By assuming that the ISM is dominated by molecular content, and a small dark matter contribution within r_half, the ratio can serve as a proxy of molecular gas mass fraction. Then the molecular gas mass fraction of the H-ATLAS SMGs is thus ~34% with a significant uncertainty, yet it is consistent with Bo13’s average value computed from M_gas /(M_gas+M_∗), in which M_∗ is the stellar mass.

It has been proposed that there exists a simple linear correlation between $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ and ΔV of the CO(1–0) line (e.g. Harris et al. 2012; Bo13; Goto & Toft 2015; Dannerbauer et al. 2017), $L_{CO}^{'} \propto (σ^{2} R) / (α_{CO} G)$ $\hbox{$\lpco \propto (\sigma^2 R)/(\alpha_\mathrm{CO} G)$}$ , where σ is the velocity dispersion of the CO line, R is the CO emitting radius, α_CO is the CO luminosity to gas mass conversion factor and G is the gravitational constant. We recall from Eqs. (4) and (5) that the dynamical mass is proportional to σ²R and M_gas = $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ α_CO; thus this correlation simply reflects the variation of the ratio between gas mass and dynamical mass. In Fig. 6, we overlay both the apparent CO line luminosity μ $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ and the intrinsic $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ plotted against the CO linewidth on those of the Bo13’s unlensed sources. The flat distribution of the lensed sources in the upper panel of Fig. 6 shows clearly the lensed feature. After correcting for the magnification, our sources are generally within the 2 σ regions from Bo13’s fit. However, it is clear that all the sources with ΔV_CO < 400 km s^-1 are above the correlation and very close to the + 2σ limit. Again, this supports our previous argument that these linewidths are likely being underestimated.

5. Physical properties of molecular gas

5.1. Multi-J CO line excitation and LVG modelling

As indicated by the histograms of $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ /L_IR in Fig. 4, the shape of the average CO SLED of the H-ATLAS SMGs follows the trend of other high-redshift SMGs and both of them depart from the average CO SLED of local galaxies with L_IR < 10¹² L_⊙ for the low-J (J_up = 3, 4, 5) part at the ~1σ levels. Figure 7 shows the L_IR -normalised CO SLED of the high-redshift SMGs (previous detections in the literature together with H-ATLAS ones) comparing with those of the local galaxies from Liu et al. (2015).

Fig. 7

The CO SLEDs of local galaxies and high-redshift SMGs normalised by L_IR⁴. Grey symbols indicate the Gaussian mean and deviation of the ratio $L_{CO}^{'} / L_{IR} \times J_{up}^{2}$ $\hbox{$L'_\mathrm{CO}/L_\mathrm{IR} \times J_\mathrm{up}^2$}$ for local galaxies mostly with 10⁹ L_⊙ ≤ L_IR ≤ 10¹² L_⊙ (Liu et al. 2015). We also include a typical local ULIRG, Arp 220 (purple dashed line and open triangles, Rangwala et al. 2011), an AGN-dominated source, Mrk 231 (red dashed line and open diamond, van der Werf et al. 2010), and a LIRG, NGC 7771 (green dashed line and open square, Liu et al. 2015). Red symbols show high-redshift SMGs from both Carilli & Walter (2013) and from this work. The dots without error bars (J_up = 10, 11) indicate that the volume of the $L_{CO}^{'} / L_{IR} \times J_{up}^{2}$ $\hbox{$L'_\mathrm{CO}/L_\mathrm{IR}\times J_\mathrm{up}^2$}$ values is insufficient for a normal-distribution fitting. We thus only indicate their mean value here. It is clear that the ratio of $L_{CO}^{'} / L_{IR} \times J_{up}^{2}$ $\hbox{$L'_\mathrm{CO}/L_\mathrm{IR}\times J_\mathrm{up}^2$}$ decreases with increasing J_up for local galaxies while it remains flat for high-redshift SMGs and a typical starburst-dominated ULIRG, Arp 220.

Previous studies of global CO excitation in both local and high-redshift galaxies (e.g. Weiß et al. 2007; Rangwala et al. 2011; Deane et al. 2013; Papadopoulos et al. 2014; Zhang et al. 2014; Spilker et al. 2014; Liu et al. 2015; Daddi et al. 2015) show that there are most likely two excitation components dominating the CO emission from ground level up to J_up = 11, a low-excitation component peaking around J_up = 3 to J_up = 4 and a high-excitation component peaking at $\hbox{$J_\mathrm{up} \gtrapprox 6$}$ . Rosenberg et al. (2015) further quantitively classify the local galaxies into three groups based on the shape of their CO SLEDs, which provides clues towards the dominant excitation conditions within.Comparing the average L_IR -normalised CO SLED of local galaxies (dominated by normal star-forming galaxies and LIRGs with L_IR = 10⁹–10¹² L_⊙) with that of the SMGs in Fig. 7, it is clear that the low-excitation component is more prominent in local galaxies, resulting in the average CO SLED peaking at J_up = 3 or J_up = 4, and decreasing with increasing energy levels (as in the case of a local LIRG, NGC 7771, shown in Fig. 7). For the SMGs, the low-excitation component is rather weak while the high-excitation component is comparable to the local normal star-forming galaxies and LIRGs, resulting in a rather flat SLED. To compare the L_IR -normalised CO SLED of high-redshift SMGs with that of the local ULIRGs, we also overplot a typical local non-AGN-dominated ULIRG (classified as a class II galaxy which is dominated by starburst as in Rosenberg et al. 2015), that is, Arp 220. It is found that the average L_IR -normalised CO SLED of high-redshift SMGs agrees well with that of Arp 220. Since the average values of J_up = 10 and J_up = 11 of high-redshift SMGs are calculated based on only a few sources, the deviations between Arp 220 and high-redshift SMGs for these two lines are not significant from a statistical point of view. Nevertheless, for the AGN-dominated ULIRG, for example, Mrk 231 as shown in Fig. 7 (classified as a class-III galaxy which is dominated by AGN powering, Rosenberg et al. 2015), the L_IR -normalised CO SLED is below that of the high-redshift SMGs. This shows that AGN are contributing much less of the L_IR luminosity of our high-redshift SMGs compared to the AGN-dominated ULIRGs. Thus, in the high-redshift SMGs, the average CO gas excitation conditions are likely to be similar to those of local non-AGN-dominated ULIRGs.

To further investigate the CO line excitation and extract the information of physical conditions of the molecular gas, we apply a large velocity gradient (LVG) statistical equilibrium method (e.g. Sobolev 1960; Goldreich & Kwan 1974; Scoville & Solomon 1974) for modelling the fluxes of multiple CO lines. We adopt a one-dimensional (1D) non-LTE radiative transfer code developed by van der Tak et al. (2007), that is, RADEX, with an escape probability of β = (1−e^τ) /τ derived from an expanding sphere geometry. The CO collisional data are taken from the LAMDA database (Schöier et al. 2005).

As a first step, we use one excitation component in the LVG modelling. Similar to Weiß et al. (2007), the inputs of the code are the molecular gas kinetic temperature (T_k), the volume density of the molecular hydrogen (n_H₂), the column density of the CO molecule (N_CO), and the solid angle (Ω_app, note that this solid angle includes the lensing magnification factor) of the source which scales with the resulting fluxes from each CO transition equally, so that the shape of the CO SLED only depends on T_k , n_H₂ and N_CO . We fix the velocity gradient to 1 km s^-1, so that the actual input of N_CO is column density per unit velocity gradient N_CO /dv instead.

A Bayesian approach is used to fit our observed flux to the fluxes generated from RADEX models given the parameters p (model parameter p includes T_k , n_H₂ , N_CO /dv and Ω_app). We use the code emcee (Foreman-Mackey et al. 2013⁵) to perform the Markov chain Monte Carlo (MCMC) calculation with the affine-invariant ensemble sampler (Goodman & Weare 2010). The Bayesian posterior probability of the model parameters given our data I^data can thus be written as following (e.g. the notation in Wall & Jenkins 2012): $\Pr (p | I^{data}) = \frac{\Pr (p) \Pr (I^{data} | p)}{\Pr (I^{data})},$ $\begin{equation} {\rm Pr}(p|I^\mathrm{data}) = \frac{{\rm Pr}(p){\rm Pr}(I^\mathrm{data}|p)}{{\rm Pr}(I^\mathrm{data})}, \label{eq:bayesian} \end{equation}$ (6)in which Pr(p | I^data) is the posterior probability of the parameter p given: the prior probability of p as Pr(p); the likelihood of the resulting CO flux I^data given the parameter inputs p as Pr(I^data | p); and the probability of the data Pr(I^data), also called evidence, which is commonly treated as a normalising factor. By assuming the noise is independent Gaussian centred, we can write the likelihood as the product of Gaussian probability distributions, $\Pr (I^{data} | p) = \underset{i}{􏽙} \frac{1}{\sqrt{2 π σ_{i}^{2}}} \exp (- \frac{{(I_{i}^{data} - I_{i}^{model} (p)^{)}}^{2}}{2 σ_{i}^{2}}),$ $\begin{equation} {{\rm Pr}}(I^\mathrm{data}|p) = {\prod_{i}} \frac{1}{\sqrt{2 \pi \sigma_i^2}}\exp \left({-\frac{\left(I^\mathrm{data}_i-I^\mathrm{model}_i(p)\right)^2}{2\sigma_{i}^2} }\right), \label{eq:likelihood} \end{equation}$ (7)where σ_i is the error associated with each set of measured fluxes $I_{i}^{data}$ $\hbox{$I_i^\mathrm{data}$}$ , and the RADEX-generated results given a set of input parameters p is given by $I_{i}^{model} (p)$ $\hbox{$I^\mathrm{model}_i(p)$}$ . We note that we use the logarithmic form of Eq. (7) in our practical calculations for convenience, so that the resulting parameters are all in logarithmic form.

Rather than generating a grid of line fluxes for a range of input parameters from the LVG models (e.g. Kamenetzky et al. 2011; Krips et al. 2011; Spinoglio et al. 2012), we directly use the Python package emcee to call pyradex (a Python wrapper of RADEX written by A. Ginsburg⁶) in each iteration for computing the RADEX results and passing them to Python, and sample the posterior probability distribution function. This can avoid calculations in the unfavourable part of the parameter space, thus saving the total running time of the codes. Following previous works (e.g. Spilker et al. 2014), we adopted flat log-prior within physically reasonable ranges, which are boundaries of the parameter space that we explored. The prior possibilities outside the boundary are set to 0. The parameter-space boundaries are as follows: n_H₂ = 10²–10⁷ cm^-3, T_k = T_CMB–10³ K, N_CO /dv = 10^15.5–10^19.5 cm^-2 km^-1 s, in which T_CMB is the CMB temperature at the redshift of the source, which can be derived from T_CMB = 2.7315(1 + z). We also adopt the range of dv/dr to be 0.1–1000 km s^-1 pc^-1 (e.g. Tunnard & Greve 2016), which limits the range of the ratio between N_CO and n_H₂ . This prior also puts limits on the ratio between the LVG-solved dv/dr and the dv/dr derived from the virialised state, that is, K_vir = (dv/dr)_LVG/(dv/dr)_vir, in which (dv/dr)_vir = 0.65α^0.5(n_H₂/ (10³ cm^-3))^-0.5 km s^-1 pc^-1, where α = 0.5–3 depending on the density profile (Papadopoulos & Seaquist 1999). Additionally, we set the priors to limit the column length to be smaller than the diameter of the entire SMG, which is about 7 kpc. This yields a constraint of the ratio N_CO /n_H₂ that is well outside the range given by the prior from dv/dr. Lastly, the molecular gas mass traced by the CO lines should not exceed the dynamical mass, for example, ~10¹² M_⊙ (see Sect. 4.2). This yields an upper limit of N_CO /dv to be smaller than ~10²⁰ cm^-2 km^-1 s (e.g. Rangwala et al. 2011), which is well outside the parameter space as well.

Table 6

Single-component MCMC-resulting molecular gas properties of the H-ATLAS SMGs.

A total of 400 walkers have been deployed to explore the parameter space initiated from the point of solution acquired by the quasi-Newton solver. We ensured proper convergence of the MCMC chains by a burn-in period of 100 iterations and 1000 subsequent iterations. The resulting posterior probability distributions and the marginal distribution of the parameters (generated by corner.py, Foreman-Mackey 2016) are shown by the blue density-contour plots and the blue histograms in Fig. C.1. We also indicate the 39% and 68% quantiles of the marginalised probability distribution of the parameter with dashed lines. The solutions with the maximum posterior probability within the 39% and 68% quantiles (the ±1σ range) are marked with orange lines and points. The corresponding fit to the CO SLED is also shown in the figure with an orange line overlaid on the black data points. All the results, the median value, the ±1σ range and the maximum posterior probability, are summarised in Table 6, except for NAv1.195 because only one CO line of it has been observed, leading to an unreliable fitting.

From the single excitation component fitting, the range of n_H₂ is found to be ≈10^2.5–10^4.1 cm^-3, T_k is from 22 K to 750 K and N_CO /dv ≈ 10^17.13–10^18.22 cm^-2 km^-1 s for the H-ATLAS SMGs. In most cases, the values are close to those found by single-component LVG modelling of local ULIRGs (e.g. Ao et al. 2008) and high-redshift SMGs (e.g. Lestrade et al. 2010; Combes et al. 2012; Riechers et al. 2013). Nevertheless, one should note that the observed CO SLEDs are dominated by the excitation from dense and warm molecular gas as suggested by the CO SLEDs peaking around J_up = 4–7. Our observed CO SLEDs are biased towards the mid- and high-J CO lines, and thus a single component fit is biased towards the high-excitation component seen in local ULIRGs. Indeed, most of our values of T_k from the single component analysis are higher than the low-excitation component seen in local ULIRGs but close to the values of the high-excitation component, for example, the warm molecular gas as found in Arp 220 (Rangwala et al. 2011). The lack of J_up ≤ 2 data will likely lead to overestimations of the values of T_k for our SMGs. The under-presentation of the low-excitation component is also shown in the fitted CO SLED to the observed flux: the modelled fluxes of CO(3–2) and CO(1–0) are often underestimated, especially in the cases of G09v1.40, SDP 17b, G12v2.43, NAv1.144 and G12v2.890 shown in Fig. C.1.

Therefore, to fully consider both the low-excitation component with a cooler temperature and the warmer, dense high-excitation component, we perform a two-excitation-component LVG modelling with the CO SLEDs: a low-excitation component with a lower value of T_k and a high-excitation component with a higher T_k . An MCMC method similar to the aforementioned single-component LVG modelling is adopted. We assign two sets of n_H₂ , T_k , N_CO /dv and Ω_app to the two excitation components, so that the two components can have different physical conditions. Similar priors are also applied to help constrain the posterior distribution. For the two-excitation-components fit, the number of free parameters can be reduced significantly in some sources, therefore we carefully add some more informative priors. Besides the similar priors such as those used in the single-component analysis, we put additional prior constraints on the two-component LVG modelling as follows: (i) The size of the cooler low-excitation component should be larger than that of the high-excitation component (this has been suggested by the observations of the sizes of the emitting regions of different transitions of the CO lines, e.g. Ivison et al. 2011); and (ii) most importantly, we assume that the temperature of the low-excitation component should be close to the cold dust temperature. At high densities (n_H₂ ≥ 10^4.5 cm^-3), the temperatures of the dust and the gas can be well coupled (Goldsmith 2001). However, the different heating mechanisms do not necessarily produce thermally balanced dust and gas temperature, especially at lower densities (n_H₂ ≤ 10^3.5 cm^-3). So we use rather loose priors for the T_k of the cool low-excitation component, that is, a normal distribution with the value of mean and standard deviation equal to the cold dust temperature Pr(T_k)_cool ~ $\hbox{$\mathcal{N}$}$ (T_d, $T_{d}^{2}$ $\hbox{$T_{\rm d}^2$}$ ). This prior offers a reasonable guess of the T_k for the cooler component in the range of ~0 K–90 K (which will be further reduced into the range of T_CMB–90 K by the aforementioned priors).

We have applied similar MCMC walkers in the parameter space and generated the posterior probability distributions (Fig. C.1). The dark green contours and histograms are the posterior probability distribution and marginal probability distribution for n_H₂ , T_k and N_CO /dv of the low-excitation component, while those shown in light green are for the high-excitation component. The solutions with the maximum posterior probability are shown by the dotted-dash pink line and red dashed line for the low- and high-excitation components, respectively. The corresponding CO SLEDs are over-plotted in Fig. C.1 with dotted-dash pink line and the dashed red lines. The values are summarised in Table 7. We show the histograms of the derived parameters for low- and high-excitation components in Fig. 8, together with those of the single component LVG modelling. For the low-excitation component, the density ranges from ~10^2.8 cm^-3 to 10^4.6 cm^-3 with large uncertainties, the gas temperature ranges from ~20 K to 30 K, and the CO column density per unit velocity gradient ranges from 10^15.7 cm^-2 km^-1 s to 10^17.9 km^-1 s with significant uncertainties. For the high-excitation component, the density ranges from ~10^2.8 cm^-3 to 10^4.2 cm^-3, the gas temperature ranges from ~60 K to 400 K, and the CO column density per unit velocity gradient ranges from 10^17.1 cm^-2 km^-1 s to 10^18.1 cm^-2 km^-1 s. As a comparison, the gas densities of the two components are close, and the differences are within uncertainties. The gas temperatures of the high-excitation component are higher (peaking around ~200 K) than those of the low-excitation ones (peaking around ~25 K). The T_k of the cooler component is ~10–15 K lower than the dust temperature T_d as shown in Table 3.

Table 7

Two-component MCMC-resulted molecular gas properties of the H-ATLAS SMGs.

As shown in the two-component-model-produced CO SLEDs in Fig. C.1, after including the low-excitation component, the flux of the low-J CO lines can be better reproduced, especially in G09v1.40, SDP 17b, G12v2.43, and NAv1.144. This indicates that to fully explain the CO SLEDs, at least two excitation components are needed. The gas temperatures also decreased comparing the T_k from single component LVG modelling with the T_k of the high-excitation component from two-component analysis, as in the case of G09v1.97, G09v1.40, G12v2.43, G12v2.30, NCv1.143 and NAv1.56. This again suggests that bias could be introduced using only a single excitation component to explain the full CO SLED. The physical properties of both components agree well with other studies of high-redshift SMGs (e.g. Scott et al. 2011; Danielson et al. 2013; Spilker et al. 2014). One has a low gas temperature, while another is smaller in size but with a warmer gas temperature, supporting the idea that there is likely more than one molecular gas excitation component. The latter is thought to be more closely related to the on-going star-forming activity compared with the cooler component. We also note that in the cases of NBv1.78, NAv1.56, G12v2.890 and G15v2.779, due to the very limited number of data points, the two-component fitting is highly reliant upon the priors and does not produce better fitting results than the single-component fittings. Thus, those individual two-component fitting results should be used with caution. Nevertheless, here, for the purpose of a statistical analysis of the physical properties of the gas excitation, we include these results in the discussion below.

We have investigated whether or not the CO linewidth correlates with the derived gas excitation condition. By comparing the linewidth with the parameters derived from LVG modelling, no significant correlation is found, with the absolute values of the correlation coefficient ≲ 0.3.

Fig. 8

Histograms of n_H₂ , T_k and N_CO /dv derived from single-component LVG modelling (black), and those of the low- and high-excitation components from two-component LVG modelling as shown in blue and red, respectively.

To further investigate the physical properties from a statistical point of view, we plot the gas thermal pressure P_th (defined by P_th ≡ n_H₂ × T_k) sampled from the MCMC posterior probability distributions versus the star formation efficiency (SFE) proxy defined as the ratio between L_IR and molecular gas mass in Fig. 9. Because L_IR is proportional to the star formation rate (e.g. Kennicutt & Evans 2012), the ratio between L_IR and M_gas is thus a good representative of SFE as defined by $SFE \equiv SFR / M_{gas},$ $\begin{equation} \textit{SFE} \equiv \textit{SFR}/\mgas, \label{eq:sfe} \end{equation}$ (8)representing the SFR per unit molecular gas mass. As displayed in the left panel of Fig. 9, for comparison, we also include the values of the Milky Way (Draine 2011), the Tuffy galaxies (Zhu et al. 2007), the Antennae galaxies (Gao et al. 2001; Zhu et al. 2003), Arp 220 (Rangwala et al. 2011) and the z ~ 6 SMG HFLS3 (Riechers et al. 2013) for comparison. For the gas thermal pressure P_th versus SFE of the H-ATLAS SMGs only (as indicated by the grey dashed square in the left panel of Fig. 9), we find the Pearson’s correlation coefficient R_P = 0.68 (with the p-value p = 0.003), indicating the existence of a strong correlation. The values of the local ULIRG Arp 220 and the high-redshift SMG HFLS3 also follow the similar relation we find in the plot, with similar values of P_th and SFE as for the H-ATLAS ULIRGs. However, the dynamical range of the P_th is small for the SMGs and local ULIRGs; thus the correlation may be biased by a few sources with either very large or very small values of P_th. Therefore, we include some nearby galaxies with lower SFE: the Milky Way, the Tuffy galaxies and the Antennae galaxies. After including these sources, we find the Pearson’s correlation coefficient increases to R_P = 0.89 (p = 0.001). With a fit to all the data points from the single component LVG modelling, we get a slope of 1.1 ± 0.5. We have also tested whether this correlation could arise from any relation between n_H₂ and SFE or between T_k and SFE. For these two pairs of quantities, the correlation coefficient is much weaker, R_P< 0.33 (p> 0.14), compared with the one between the gas thermal pressure P_th and SFE. This rules out the possibility that either n_H₂ or T_k is dominating the correlation. All these pieces of evidence point out that there is likely a strong close-to-linear correlation between P_th of the molecular gas and SFE, suggesting that the thermal pressure of the bulk of molecular gas is playing an important role in regulating the star formation at galactic scale across a range of redshifts. This is inline with the theoretical works discussing the relation between gas pressure and SFE (e.g. Elmegreen & Efremov 1997; Wong & Blitz 2002).

Fig. 9

Left: thermal gas pressure plotted against star formation efficiency as indicated by L_IR /M_gas. Filled circles are the H-ATLAS SMGs. We also plot the values of the Milky Way, the Tuffy galaxies, the Antennae galaxies, Arp 220 and HFLS3 in filled triangles for comparison. The colour is coded based on their values of n_H₂ . The red line shows the fit to the correlation, which yields a slope of 1.1 ± 0.5 for all the sources. The grey dashed square shows the region where H-ATLAS SMGs reside. Middle: similar correlation plot as in the left panel but only for the cool component from the two-component LVG modelling. There is no correlation found (R_P = 0.01 and p = 0.97). Right: similar correlation plot as in the left panel but only for the warm component in the two-component LVG modelling. The dashed red line is an overlay of exactly the same red line plotted in the left panel. The data points follow the same correlation found in the single component fit with an R_P = 0.37 and p = 0.16. In all the three panels, the legends show the Pearson’s correlation coefficient R_P with corresponding p-value. The values of n_H₂ for each point are indicated by the colour bar.

After decomposing the CO excitation into the low-excitation component and the high-excitation one, we also plot their P_th versus SFE separately. The results of the (cool) low-excitation component and the (warm) high-excitation component are shown in the middle and right panel of Fig. 9, respectively. There is likely no correlation between P_th and SFE for the low-excitation component, as suggested by the low coefficient R_P = 0.01 (p = 0.97). Nevertheless, for the high-excitation component there is still evidence of a correlation between P_th and SFE with R_P = 0.37 (p = 0.16). Although, a large dynamical range and smaller uncertainties of P_th and SFE are needed to further confirm this correlation. We also overlay the red line from the correlation fit to P_th versus SFE of the single component modelling. It is clearly shown in the plot that, the data points of the high-excitation component follow this correlation well. This again shows that the high-excitation component is more closely related to the on-going star formation activity, while the low-excitation gas is much less tied to star formation.

5.2. The [C I](2–1) emission line in the high-redshift SMGs

The ³P fine structure [C I] lines of atomic carbon, [C I](1–0) at 492.2 GHz and [C I](2–1) at 809.3 GHz, are in a simple three-level energy system. The critical densities n_crit for the [C I](1–0) and [C I](2–1) lines are both ~0.5–1 ×10³ cm^-3, which is comparable to those of the CO(1–0) and CO(2–1) lines (Table 1). The energy levels for [C I]³P₂ and [C I]³P₁ are 23.6 K and 62.6 K, respectively. Atomic carbon is found to be well-mixed with the bulk of the H₂ gas, making it a promising molecular gas tracer together with low-J CO lines (e.g. Papadopoulos & Greve 2004). From these two optically thin [C I] lines, one can derive the excitation temperature and gas density without relying on other complementary information (e.g. Weiß et al. 2003). In the high-redshift Universe, [C I] has only been detected in a small number of systems, mainly gravitationally lensed SMGs and quasars (see the references in Walter et al. 2011; see Bothwell et al. 2017, for recent detections of [C I](1–0) lines in a sample of SPT-selected lensed SMGs; see also Wilson et al. 2017, for detections of [C I](2–1) by stacking Herschel SPIRE/FTS spectral data of galaxies at moderate redshifts). In this work, we have detected seven [C I](2–1) lines out of eight observations in our lensed high-redshift SMGs.

Papadopoulos & Greve (2004) find a good agreement between the total molecular gas mass derived from [C I] and CO lines and dust continuum in local ULIRGs. [C I] likely traces H₂ even more robustly than the low-J CO lines (plus standard conversion factor) in extreme conditions on galactic scales (e.g. Zhang et al. 2014). High-redshift observations (e.g. Weiß et al. 2003; Alaghband-Zadeh et al. 2013) also support such an agreement, although Bothwell et al. (2017) recently found that either a larger α_CO or a high [C I] abundance is needed to balance the gas mass derived from CO and from [C I] lines. Since the observations of CO(1–0) lines at high-redshift becomes observationally difficult, acquiring the intensities of the [C I] lines, which are usually brighter than the CO(1–0) lines and residing in favourable bands for observation, could help us to better determine the total molecular mass in high-redshift SMGs. Zhang et al. (2016) show that at high-redshift the CO(1–0) line will also suffer observing against the CMB, making the line more difficult to observed. Moreover, they also found that CO can be destroyed by the cosmic rays coming from the intense star-forming activities, suggesting the [C I] lines to be a better tracer of the total molecular gas mass in such environments (Bisbas et al. 2015).

Fig. 10

$L_{[C I] (2 - 1)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{[C{\scriptsize{\,I}}]}(2\text{--}1)}$}$ versus $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ including (U)LIRGs from Jiao et al. (2017; dark purple) H-ATLAS SMGs (red) and other SMGs from Walter et al. (2011; green). Two of our sources, NAv1.177 (#9) and G12v2.257 (#15) are lacking lens modelling, thus we plot the limit, and the arrows show the direction for lensing correction. A fit to the correlation from the local to high-redshift galaxies is indicated by the blue line, with a slope of 0.94 ± 0.05. We also show a fit with a fixed slope of 1, that is, a linear fit, by the orange dashed line. The source index ID is indicated in Fig. 1 and Table 3.

Using Herschel spectral data of a sample of local (U)LIRGs, Jiao et al. (2017) find a tight correlation between the CO(1–0) line luminosity and the [C I] line luminosities. The tight correlations suggest that the [C I] lines trace the total molecular gas mass as the CO(1–0) line does. Using the empirical correlations, Jiao et al. (2017) derive a conversion factor α_{[C I] (2–1)} that converts $L_{[C I] (2 - 1)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{[C{\scriptsize{\,I}}]}(2\text{--}1)}$}$ into M_gas to be α_{[C I] (2–1)} = 27.5 ± 1.3 ( K km s^-1 pc²)^-1. To test whether this correlation can be extended to high-redshift SMGs, we plot the $L_{[C I] (2 - 1)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{[C{\scriptsize{\,I}}]}(2\text{--}1)}$}$ versus $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ including our lensed SMGs in Fig. 10. The fit to all the galaxies across different redshifts indicates a tight correlation close to linear, that is, the slope equals 0.94 ± 0.05. When fixing the slope to one, we find the fitted linear correlation to be valid both for local sources and high-redshift SMGs. Nevertheless, one would need to know the abundance of [C I] and the excitation temperature for the [C I] lines to properly derive the molecular gas mass from the [C I] lines as in Weiß et al. (2003). This would require the detection of both the [C I](1–0) and [C I](2–1) lines for the H-ATLAS SMGs, while we have only detection of the 2–1 line.

5.3. Star formation and the molecular gas content in high-redshift GMGs

One of the key interests of studying star formation at galactic scale across all redshifts from the observation point of view is to quantify the “star formation law”, which is the correlation between M_gas and SFR (or the surface SFR density as originally defined by Kennicutt 1998b). This law is not only an essential input for the theoretical models of galaxy formation and evolution but also an important tool to study the relation between star formation and the molecular gas content. As we mentioned in Sect. 5.1, the SFR per unit molecular gas mass (i.e. SFE) shows how efficiently each unit of molecular gas mass is converted into stars (Eq. (8)). The inverse of SFE can be defined as the molecular gas depletion time, $t_{dep} \equiv M_{gas} / SFR,$ $\begin{equation} {t_\mathrm{dep}} \equiv \mgas /\textit{SFR}, \label{eq:t_dep} \end{equation}$ (9)which indicates the exhausting time-scale of the molecular gas mass with the current SFR. The gas depletion time is a good way of representing the variations of the star formation properties. The value of t_dep varies from ~1.3–1.5 Gyr for local star-forming galaxies (e.g. Kennicutt 1998b; Daddi et al. 2010; Genzel et al. 2010; Saintonge et al. 2011) to smaller values for high-redshift ones (“main sequence”): ~0.7 Gyr (e.g. Tacconi et al. 2013; Sargent et al. 2014) or even smaller values (Saintonge et al. 2013). Indeed, there seems to be a cosmic evolution trend for star-forming galaxies as found by Saintonge et al. (2013). In Fig. 11 (extended from Fig. 7 of Aravena et al. 2016b), we plot the gas depletion time of our sources (see values in Table 5) and compare it with different kinds of galaxies across a range of redshifts, including; the aforementioned evolution track with z of star-forming galaxies (Saintonge et al. 2013); nearby ULIRGs (Solomon et al. 1997) and z = 0.6–1 ULIRGs (Combes et al. 2013); other lensed/unlensed SMGs/ULIRGs studied (see Aravena et al. 2016b, and the references within). As shown in Fig. 11, unlike star-forming galaxies, for all types of ULIRGs, the depletion time is found to be much smaller, from ~10 to 100 Myr (roughly 10 times smaller than that of the star-forming main sequence galaxies), and there seems to be no evidence of a cosmic evolution of t_dep. The H-ATLAS SMGs show no difference in the t_dep compared with SPT lensed ones and other SMGs/ULIRGs studied. They are well below the values for star-forming main sequence galaxies and show no evidence of variation across redshifts.

Fig. 11

Extended from Fig. 7 of Aravena et al. (2016b): molecular gas depletion time of our lensed SMGs (black circles), SPT lensed SMGs (green diamonds), and other unlensed SMGs (orange stars). The ranges of local ULIRGs (Solomon et al. 1997) and z = 0.6–1 ULIRGs (Combes et al. 2013) are also included as the red and purple regions, respectively. The cyan region indicates the main sequence galaxies as described by Saintonge et al. (2013) through the formula t_dep = 1.5(1 + z)^α, in which α is from α = −1.5 (Davé et al. 2012) to α = −1.0 (Magnelli et al. 2013).

Fig. 12

Left panel: the ±1σ range of δ_GDR and α_CO derived from observed M_dust and $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ using Eq. (10). Red region is from this work while the green region is from Aravena et al. (2016b). Our results generally agree with the SPT sources’ albeit having slightly larger α_CO or smaller δ_GDR. Right panel: the ratio between $L_{CO}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}}$}$ and M_dust of SMGs across the different redshifts. The SPT lensed SMGs in green points are from Aravena et al. (2016b), while our sources are in red. We also plot other SMGs from the literature in yellow filled squares (see references in text). The black line shows a best-fit to all the data points and dashed black lines show the ±1σ ranges (for the slope only) from the fit. The resulting relation is shown in the legend.

However, one should note that the values of M_gas derived above are based on the assumption of α_CO = 0.8 M_⊙ (K km s^-1 pc²)^-1. This conversion factor may be different in various types of galaxies. To derive reliable value of α_CO, one needs to estimate the molecular gas mass from other methods rather than converting from CO fluxes. One of the most convenient ways is to compute the molecular gas mass from the dust mass by assuming a value for the gas-to-dust mass ratio δ_GDR, since the dust continuum is much easier to measure in practical observations (e.g. Scoville et al. 2014, 2017). Nevertheless, δ_GDR has no constant value for different galaxies. In fact, both δ_GDR and α_CO are functions of the metallicity of different galaxies based on both observational results (e.g. Wilson 1995; Leroy et al. 2011; Magdis et al. 2011; Genzel et al. 2012) and theoretical works (e.g. Narayanan et al. 2012). To avoid the uncertainties caused by assuming single values for δ_GDR and α_CO, we explore a combination of both, which is derived directly from observational quantities. From Eqs. (2) and (3), we can derive $δ_{GDR} / α_{CO} = L_{CO (1 - 0)}^{'} / M_{dust},$ $\begin{equation} \label{eq:alpha_co_gdr} \delta_\mathrm{GDR}/\alpha_\mathrm{CO} = \lpcol10/ \mdust , \end{equation}$ (10)in which $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ is from CO line observations (converted mostly from CO(3–2) in our cases as mentioned in Sect. 4.1), and M_dust is calculated from the observed submm/mm dust continuum flux densities based on a modified black body model (both for our work and the SPT lensed SMGs).

The left panel of Fig. 12 shows the possible ±1σ range of δ_GDR versus α_CO for H-ATLAS SMGs computed using Eq. (10), together with the one from SPT-lensed SMGs for comparison. By taking a common δ_GDR of 100, the corresponding α_CO for our H-ATLAS SMGs is ~0.7–2 M_⊙ (K km s^-1 pc²)^-1, while for the SPT sources, α_CO is slightly smaller, ~0.4–1 M_⊙ (K km s^-1 pc²)^-1. Nevertheless, one need to accurately measure either δ_GDR or α_CO in order to break the degeneracies between these two quantities. Normally, one could use the calibration of the relation between the metallicity Z and δ_GDR (e.g. Magdis et al. 2011), to derive the latter from the former by observing the optical emission lines. But for the high-redshift SMGs, due to the extreme dust obscuration, it is very challenging to measure the metallicity from optical observation, making it thus difficult to accurately pinpoint the value of δ_GDR and α_CO.

By including other published SMGs at different redshifts (Magdis et al. 2011; Combes et al. 2012; Salomé et al. 2012; Walter et al. 2012; Riechers et al. 2013), we have plotted the $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ /M_dust versus redshift for all of them in the right panel of Fig. 12 similarly to Tan et al. (2014). There seems to be an increasing trend for the SMGs across the cosmic time, from nearby to the high-redshift Universe. The best fit to this trend is $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ /M_dust ∝ (1 + z)^{0.9 ± 0.5}, in the redshift range from z ≈ 2 up to z ≈ 6. This agrees with the fact that the average SPT sources have higher values of δ_GDR/α_CO (i.e. $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ /M_dust), since the SPT sources have a higher average redshift compared with the H-ATLAS one. It has been suggested that δ_GDR is linearly anti-correlated with Z, that is, δ_GDR ∝ Z^-1 (e.g. Santini et al. 2010; Leroy et al. 2011; Magdis et al. 2011), while the dependence of α_CO seems to be steeper (e.g. α_CO ∝ Z^-1.4, calibration from Leroy et al. 2011). Combining these two calibrations, one would derive a correlation between the ratio δ_GDR/α_CO ( $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ /M_dust) and Z (e.g. Magdis et al. 2012). This offers one possible explanation for the observed correlation, if the higher-redshifted SMGs in this plot have larger values of metallicity Z. Considering the cosmic enrichment of metallicity at high-redshift (e.g. Troncoso et al. 2014), the higher-redshifted SMGs in the plot are unlikely to be evolutionally linked with the lower redshifted ones. At given redshifts, Z increases with stellar masses (see a summary in Tan et al. 2014). Therefore, this increasing value of $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ /M_dust could rather be explained by a selection bias towards the high-mass systems which correspond to those higher-redshifted SMGs in the figure.

5.4. CO gas and H₂O gas comparison

The systematic study of local galaxies (from normal star-forming galaxies to nearby ULIRGs) shows the close relation between the submm H₂O emission and the star-forming activity (Yang et al. 2013). A similar conclusion has been extended to the high-redshift by the study of a group of lensed SMGs (O13 and Y16). The submm H₂O lines are dominated by FIR pumping, which is closely related to the warm dust (T_d ~ 40–90 K, see e.g. González-Alfonso et al. 2014). One may compare the gas content traced by these submm H₂O lines with the CO lines we observed. The H₂O column density derived is around ~0.3 × 10¹⁵ cm^-2 km^-1 s to ~2 × 10¹⁶ cm^-2 km^-1 s (Y16), which is about several tens up to several thousand times smaller than the CO column density (from single component LVG modelling).

Fig. 13

Histogram showing the distribution of the ratio (indicated by black data points) between the linewidths of CO and H₂O. The H₂O linewidth are from Y16 and O13. It is clear that most of our sources have similar CO and H₂O linewidth, with an average ratio of 1.0 ± 0.1. For G09v1.97, the IRAM-30 m spectrum shows that the asymmetric CO line is ~1.5 times larger than that of the H₂O line (Y16) as shown by the dashed error bar and histogram. However, using the ALMA high S/N CO and H₂O spectral data of G09v1.97 (Yang et al., in prep.), we find a ratio of the linewidth to be very close to 1 as indicated by the red point. The source of SDP 81 is not shown here because it is a special case as discussed in the text.

Comparing the linewidths between CO and H₂O is instructive. As shown in Fig. 13, most of the H-ATLAS sources have similar CO and H₂O linewidths. This suggests that they are possibly coming from similar regions as found in O13 and Y16. However, for SDP 81 and G09v1.97, the CO linewidth obtained using the IRAM-30 m data is approximately two times larger than that of the H₂O lines. As shown in O13, the blue component of the CO line of SDP 81 is much weaker compared with the red component, and this component is not even detected by ALMA long baseline observations (ALMA Partnership 2015). Thus, here the linewidth ratio only accounts for the H₂O linewidth of the red component, which results in the CO linewidth being about twice larger than that of the H₂O line in SDP 81. However, if the S/N of the spectral data is improved, one will detect the weakly magnified velocity component as described in Sect. 4.3. From the ALMA data of CO(6–5) and J = 2 H₂O lines of G09v1.97 (Yang et al., in prep.), we find asymmetrical CO and H₂O line profiles that are very similar, showing a dominant red velocity component with an approximately six times weaker blue velocity component (see Fig. A.1 for the ALMA CO(6–5) spectrum). The linewidth ratio, acquired by the high S/N ALMA data, is found very close to 1 as indicated by the red data point in Fig. 13. The case of G09v1.97 clearly shows that the low/mid-S/N spectral data may suffer from the bias caused by differential lensing effect as mentioned in Sect. 4.3. Besides the line profiles between H₂O and CO, Oteo et al. (2017) find similar line profiles for H₂O and mid-J CO and HCN lines in a z ~ 1.6 lensed SMG, SDP 9.

The similarity of the line profile strongly suggests that the emitting regions are co-spatially located. Both the gas tracers are closely linked to the on-going star formation activities. However, further detailed studies combining the gas excitation, dust emission, and the far-infrared pumping of the H₂O lines is needed to fully incorporate the complex physical processes in these SMGs.

6. Summary

In this work, we report a survey of multiple J_up (mostly from J_up = 3 up to J_up = 8) CO lines in a sample of Herschel-selected lensed SMGs at z ~ 2–4. We have detected 47 CO lines and 7 [C I] lines in these SMGs using the IRAM-30 m telescope.

Comparing the CO linewidth with those of the SPT lensed SMGs and an unlensed sample, we find evidence of a significant bias introduced by differential lensing that distorts the line profile, resulting in underestimation of the linewidth (usually by a factor of 2). This induces underestimation of the dynamical mass if one uses the observed linewidth measured for the lensed SMGs blindly. Differential lensing also slightly affects flux estimates of the lowest-J CO lines. This is mostly due to the difference in spatial distribution, since the lowest-J lines, especially CO(1–0) , are usually several times more extended than the mid-J to high-J CO lines. By comparing the CO(1–0) to CO(3–2) flux ratio of the SMGs for our lensed SMGs and the ones of the unlensed SMGs, we find that the differential lensing could cause a ~1.3 ± 0.4 times overestimation of the magnification of the CO(1–0) .

We also calculate the dynamical mass of the H-ATLAS SMGs; if one did not correct for differential lensing effects, this would lie in the range of 0.1–7 × 10¹¹ M_⊙, with a large uncertainty due to the unknown extension and sky inclination (if assuming a rotating-disk scenario) of the CO emission. But for the sources with a narrow linewidth, due to the aforementioned differential lensing effect, their dynamical masses might be underestimated. Nevertheless, for the sources with broad linewidth, the ratio between gas mass and dynamical mass equals about 0.34 ± 0.10, that is, a gas fraction of ~34%, which is close to the values of other SMGs and empirical model predictions.

The multiple-transition CO line data allow us to study the molecular gas excitation of the SMGs. The CO SLEDs are mostly peaking around J_up = 5–7, thus dominated by the warm dense gas. Using LVG modelling by fitting the SLEDs with a single excitation component via a Bayesian approach, we derive a gas density n_H₂ ≈ 10^2.5–10^4.1 cm^-3, gas temperature T_k ≈ 20–750 K and CO column density per unit velocity gradient N_CO /dv ≈ 10^16.4–10^17.8 cm^-2 km^-1 s. But we notice that the J_up ≤ 3 CO lines are likely under-predicted by the single-component LVG model. This indicates the existence of a low-excitation component. We have thus performed a two-component LVG modelling and derived the gas density, temperature, and the CO column density per unit velocity gradient for each component. The low-excitation component has a cooler gas temperature of about 20–30 K with a gas density of about 10^2.8–10^4.6 cm^-3, while for the high-excitation component, the gas temperature is higher, ~60–400 K, with a gas density of about 10^2.7–10^4.2 cm^-3. We find a correlation between the gas thermal pressure P_th derived from single component LVG modelling and the star formation efficiency. The warm high-excitation component also follows a similar trend while the cool low-excitation component is much less correlated with the on-going star formation. This shows that the high-excitation warm dense gas is more closely related to the star formation than the cool gas in the SMGs. This agrees with previous studies of a large sample of local galaxies.

We also study the global properties of the molecular gas and their relation with star formation. By assuming a conversion from CO luminosity to gas mass, α_CO = 0.8 M_⊙ (K km s^-1 pc²)^-1, we derive gas to dust mass ratios in the approximate range from 30 to 100. The gas depletion time t_dep of the H-ATLAS SMGs shows no difference compared with other lensed/unlensed SMGs; the value is also close to that of the local ULIRG. Furthermore, no cosmic evolution trend is found for t_dep. However, in order to avoid the uncertainties from the assumptions of the values of α_CO (and δ_GDR), we study the quantity $L_{CO (1 - 0)}^{'}$ $\hbox{$L^{\prime}_{\mathrm{CO}(1\text{--}0)}$}$ /M_dust which is proportional to the ratio δ_GDR/α_CO. We find for the SMGs, this ratio is increasing with increasing redshift, which could be caused by a selection bias.

With the detections of seven [C I](2–1) lines in our H-ATLAS lensed SMGs, we extend the correlation between the luminosity of the CO(1–0) line and the [C I](2–1) line found in a sample of local (U)LIRGs. All of them can be well fitted with a single linear correlation, indicating that the [C I] lines can be good indicators of the total molecular gas mass, both for nearby (U)LIRGs and the high-redshift SMGs. However, [C I](1–0) data are needed to properly derive the gas mass from the [C I] lines.

Finally, we compare the linewidths of the CO and H₂O emission lines for a sample of 13 SMGs, using multi-J CO emission lines. We find that the linewidths of the CO lines and H₂O lines agree very well in most cases. This supports our previous argument that the emitting regions of the CO and H₂O lines are likely to be co-spatially located.

This work reports for the first time, a systematic study of the CO gas excitation in a sample of high-redshift lensed SMGs. However, deriving accurate values for the physical properties of the molecular gas is challenging. To model the CO excitation with two excitation components, eight free parameters are needed. A rather complete sample of the CO SLED is thus required to better constrain the models. Also, for strongly lensed systems, due to differential lensing, the line profile can be easily distorted (e.g. the case of SDP 81 – Dye et al. 2015), namely the magnification factor at different velocity components of the emission line can be rather different. One could overcome this disadvantage by increasing the S/N of the spectrum using a telescope with better sensitivities such as ALMA and future NOEMA.

¹

See http://www.iram.fr/IRAMFR/GILDAS for more information about the GILDAS softwares.

²

Report is available on the IRAM-30 m wiki page: https://www.iram.es/IRAMES/mainWiki/Iram30mEfficiencies

³

As stated in the caption of Table 5, when CO(3–2) observation is absent we have used another line with similar factors, β_CO43 = 2.65 M_⊙ (K km s^-1 pc²)^-1 for CO(4–3) and β_CO21 = 1.30 M_⊙ (K km s^-1 pc²)^-1 for CO(2–1) based on Bo13.

⁴

The actual value used in y-axis is $L_{CO}^{'} / L_{IR} \times J_{up}^{2}$ $\hbox{$L'_\mathrm{CO}/L_\mathrm{IR} \times J_\mathrm{up}^2$}$ . $J_{up}^{2}$ $\hbox{$J_\mathrm{up}^2$}$ is included so that the unit of $L_{CO}^{'} \times J_{up}^{2}$ $\hbox{$L'_\mathrm{CO}\times J_\mathrm{up}^2$}$ is comparable to velocity integrated flux density, which is Jy km s^-1.

⁵

https://github.com/dfm/emcee

⁶

https://github.com/keflavich/pyradex

Acknowledgments

We would like to thank the anonymous referee for constructive comments. The authors are grateful to the IRAM staff for their support. C.Y. thanks Claudia Marka and Nicolas Billot for their help of the IRAM-30 m/EMIR observation. The authors also thank Yinhe Zhao for kindly sharing the local [C I] line data. C.Y. and Y.G. acknowledge support by National Key R&D Program of China (2017YFA0402700) and the CAS Key Research Program of Frontier Sciences. C.Y., A.O. and Y.G. acknowledge support by NSFC grants 11311130491 and 11420101002. C.Y., A.O., A.B. and Y.G. acknowledge support from the Sino-French LIA-Origins joint exchange program. C.Y. is supported by the China Scholarship Council grant (CSC No. 201404910443). E.G.-A. is a Research Associate at the Harvard-Smithsonian Center for Astrophysics, and thanks the Spanish Ministerio de Economía y Competitividad for support under projects FIS2012-39162-C06-01 and ESP2015-65597-C4-1-R, and NASA grant ADAP NNX15AE56G. Z.Y.Z., R.J.I. and I.O. acknowledge support from ERC in the form of the Advanced Investigator Programme, 321302, COSMICISM. H.D. acknowledges financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the 2014 Ramón y Cajal program MINECO RYC-2014-15686. M.J.M. acknowledges the support of the National Science Centre, Poland through the POLONEZ grant 2015/19/P/ST9/04010. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 665778. I.R.S. acknowledges support from STFC (ST/P000541/1), the ERC Advanced Grant DustyGal (321334) and a Royal Society/Wolfson Merit Award. US participants in H-ATLAS acknowledge support from NASA through a contract from JPL. Italian participants in H-ATLAS acknowledge a financial contribution from the agreement ASI-INAF I/009/10/0. SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (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). Extensive use was made of the NASA IDL (Interactive Data Language) Astronomy User’s Library (https://idlastro.gsfc.nasa.gov) and the Coyote Graphics Library.

References

Alaghband-Zadeh, S., Chapman, S. C., Swinbank, A. M., et al. 2013, MNRAS, 435, 1493 [Google Scholar]
ALMA Partnership, Vlahakis, C., Hunter, T. R., Hodge, J. A., et al. 2015, ApJ, 808, L4 [NASA ADS] [CrossRef] [Google Scholar]
Ao, Y., Weiß, A., Downes, D., et al. 2008, A&A, 491, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Aravena, M., Decarli, R., Walter, F., et al. 2016a, ApJ, 833, 68 [NASA ADS] [CrossRef] [Google Scholar]
Aravena, M., Spilker, J. S., Bethermin, M., et al. 2016b, MNRAS, 457, 4406 [Google Scholar]
Asboth, V., Conley, A., Sayers, J., et al. 2016, MNRAS, 462, 1989 [NASA ADS] [CrossRef] [Google Scholar]
Barger, A. J., Cowie, L. L., Sanders, D. B., et al. 1998, Nature, 394, 248 [NASA ADS] [CrossRef] [Google Scholar]
Barger, A. J., Cowie, L. L., Chen, C.-C., et al. 2014, ApJ, 784, 9 [NASA ADS] [CrossRef] [Google Scholar]
Becker, R. H., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559 [NASA ADS] [CrossRef] [Google Scholar]
Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bezanson, J., Karpinski, S., Shah, V. B., & Edelman, A. 2012, ArXiv e-prints [arXiv:1209.5145] [Google Scholar]
Bisbas, T. G., Papadopoulos, P. P., & Viti, S. 2015, ApJ, 803, 37 [NASA ADS] [CrossRef] [Google Scholar]
Blain, A. W., Smail, I., Ivison, R. J., Kneib, J.-P., & Frayer, D. T. 2002, Phys. Rep., 369, 111 [NASA ADS] [CrossRef] [Google Scholar]
Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207 [Google Scholar]
Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047 (Bo13) [Google Scholar]
Bothwell, M. S., Aguirre, J. E., Aravena, M., et al. 2017, MNRAS, 466, 2825 [Google Scholar]
Bouché, N., Cresci, G., Davies, R., et al. 2007, ApJ, 671, 303 [NASA ADS] [CrossRef] [Google Scholar]
Bradford, C. M., Aguirre, J. E., Aikin, R., et al. 2009, ApJ, 705, 112 [NASA ADS] [CrossRef] [Google Scholar]
Bussmann, R. S., Pérez-Fournon, I., Amber, S., et al. 2013, ApJ, 779, 25 (Bu13) [NASA ADS] [CrossRef] [Google Scholar]
Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Carter, M., Lazareff, B., Maier, D., et al. 2012, A&A, 538, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45 [Google Scholar]
Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772 [NASA ADS] [CrossRef] [Google Scholar]
Combes, F., Rex, M., Rawle, T. D., et al. 2012, A&A, 538, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Combes, F., García-Burillo, S., Braine, J., et al. 2013, A&A, 550, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cooray, A., Calanog, J., Wardlow, J. L., et al. 2014, ApJ, 790, 40 [NASA ADS] [CrossRef] [Google Scholar]
Cox, P., Krips, M., Neri, R., et al. 2011, ApJ, 740, 63 [NASA ADS] [CrossRef] [Google Scholar]
da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13 [Google Scholar]
Daddi, E., Elbaz, D., Walter, F., et al. 2010, ApJ, 714, L118 [Google Scholar]
Daddi, E., Dannerbauer, H., Liu, D., et al. 2015, A&A, 577, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Dale, D. A., Helou, G., Contursi, A., Silbermann, N. A., & Kolhatkar, S. 2001, ApJ, 549, 215 [NASA ADS] [CrossRef] [Google Scholar]
Danielson, A. L. R., Swinbank, A. M., Smail, I., et al. 2011, MNRAS, 410, 1687 [NASA ADS] [Google Scholar]
Danielson, A. L. R., Swinbank, A. M., Smail, I., et al. 2013, MNRAS, 436, 2793 [NASA ADS] [CrossRef] [Google Scholar]
Danielson, A. L. R., Swinbank, A. M., Smail, I., et al. 2017, ApJ, 840, 78 [NASA ADS] [CrossRef] [Google Scholar]
Dannerbauer, H., Lehnert, M. D., Emonts, B. H. C., et al. 2017, A&A, 608, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98 [NASA ADS] [CrossRef] [Google Scholar]
Deane, R. P., Heywood, I., Rawlings, S., & Marshall, P. J. 2013, MNRAS, 434, 23 [NASA ADS] [CrossRef] [Google Scholar]
Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615 [NASA ADS] [CrossRef] [Google Scholar]
Draine, B. T. 2003, ARA&A, 41, 241 [NASA ADS] [CrossRef] [Google Scholar]
Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press) [Google Scholar]
Dunlop, J. S., McLure, R. J., Biggs, A. D., et al. 2017, MNRAS, 466, 861 [Google Scholar]
Dye, S., Furlanetto, C., Swinbank, A. M., et al. 2015, MNRAS, 452, 2258 [NASA ADS] [CrossRef] [Google Scholar]
Eales, S., Dunne, L., Clements, D., et al. 2010, PASP, 122, 499 [NASA ADS] [CrossRef] [Google Scholar]
Elmegreen, B. G., & Efremov, Y. N. 1997, ApJ, 480, 235 [NASA ADS] [CrossRef] [Google Scholar]
Emonts, B. H. C., Lehnert, M. D., Villar-Martín, M., et al. 2016, Science, 354, 1128 [NASA ADS] [CrossRef] [Google Scholar]
Engel, H., Tacconi, L. J., Davies, R. I., et al. 2010, ApJ, 724, 233 [NASA ADS] [CrossRef] [Google Scholar]
Fixsen, D. J., Bennett, C. L., & Mather, J. C. 1999, ApJ, 526, 207 [Google Scholar]
Foreman-Mackey, D. 2016, J. Open Source Softw., 1, 24 [Google Scholar]
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
Frayer, D. T., Harris, A. I., Baker, A. J., et al. 2011, ApJ, 726, L22 [NASA ADS] [CrossRef] [Google Scholar]
Fu, H., Jullo, E., Cooray, A., et al. 2012, ApJ, 753, 134 [NASA ADS] [CrossRef] [Google Scholar]
Gao, Y., Lo, K. Y., Lee, S.-W., & Lee, T.-H. 2001, ApJ, 548, 172 [NASA ADS] [CrossRef] [Google Scholar]
Gao, Y., Carilli, C. L., Solomon, P. M., & Vanden Bout, P. A. 2007, ApJ, 660, L93 [NASA ADS] [CrossRef] [Google Scholar]
Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091 [NASA ADS] [CrossRef] [Google Scholar]
Genzel, R., Tacconi, L. J., Combes, F., et al. 2012, ApJ, 746, 69 [NASA ADS] [CrossRef] [Google Scholar]
George, R. D., Ivison, R. J., Hopwood, R., et al. 2013, MNRAS, 436, L99 [NASA ADS] [CrossRef] [Google Scholar]
Giordano, M. 2016, ArXiv e-prints [arXiv:1610.08716] [Google Scholar]
Glover, S. C. O., Clark, P. C., Micic, M., & Molina, F. 2015, MNRAS, 448, 1607 [NASA ADS] [CrossRef] [Google Scholar]
Goldreich, P., & Kwan, J. 1974, ApJ, 189, 441 [NASA ADS] [CrossRef] [Google Scholar]
Goldsmith, P. F. 2001, ApJ, 557, 736 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
González-Alfonso, E., Fischer, J., Aalto, S., & Falstad, N. 2014, A&A, 567, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Goodman, J., & Weare, J. 2010, CAMCoS, 5, 65 [Google Scholar]
Goto, T., & Toft, S. 2015, A&A, 579, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Greve, T. R., Ivison, R. J., & Papadopoulos, P. P. 2003, ApJ, 599, 839 [NASA ADS] [CrossRef] [Google Scholar]
Greve, T. R., Leonidaki, I., Xilouris, E. M., et al. 2014, ApJ, 794, 142 [NASA ADS] [CrossRef] [Google Scholar]
Harris, A. I., Baker, A. J., Frayer, D. T., et al. 2012, ApJ, 752, 152 [NASA ADS] [CrossRef] [Google Scholar]
Hezaveh, Y. D., Marrone, D. P., & Holder, G. P. 2012, ApJ, 761, 20 [NASA ADS] [CrossRef] [Google Scholar]
Hezaveh, Y. D., Marrone, D. P., Fassnacht, C. D., et al. 2013, ApJ, 767, 132 [NASA ADS] [CrossRef] [Google Scholar]
Hodge, J. A., Carilli, C. L., Walter, F., et al. 2012, ApJ, 760, 11 [Google Scholar]
Holland, W. S., Robson, E. I., Gear, W. K., et al. 1999, MNRAS, 303, 659 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Hughes, D. H., Serjeant, S., Dunlop, J., et al. 1998, Nature, 394, 241 [NASA ADS] [CrossRef] [Google Scholar]
Ikeda, M., Oka, T., Tatematsu, K., Sekimoto, Y., & Yamamoto, S. 2002, ApJS, 139, 467 [NASA ADS] [CrossRef] [Google Scholar]
Ivison, R. J., Smail, I., Papadopoulos, P. P., et al. 2010, MNRAS, 404, 198 [NASA ADS] [Google Scholar]
Ivison, R. J., Papadopoulos, P. P., Smail, I., et al. 2011, MNRAS, 412, 1913 [Google Scholar]
Jiao, Q., Zhao, Y., Zhu, M., et al. 2017, ApJ, 840, L18 [Google Scholar]
Kamenetzky, J., Glenn, J., Maloney, P. R., et al. 2011, ApJ, 731, 83 [NASA ADS] [CrossRef] [Google Scholar]
Kamenetzky, J., Rangwala, N., Glenn, J., Maloney, P. R., & Conley, A. 2016, ApJ, 829, 93 [NASA ADS] [CrossRef] [Google Scholar]
Karim, A., Swinbank, A. M., Hodge, J. A., et al. 2013, MNRAS, 432, 2 [NASA ADS] [CrossRef] [Google Scholar]
Kennicutt, Jr., R. C. 1998a, ARA&A, 36, 189 [Google Scholar]
Kennicutt, Jr., R. C. 1998b, ApJ, 498, 541 [NASA ADS] [CrossRef] [Google Scholar]
Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531 [NASA ADS] [CrossRef] [Google Scholar]
Koprowski, M. P., Dunlop, J. S., Michałowski, M. J., Cirasuolo, M., & Bowler, R. A. A. 2014, MNRAS, 444, 117 [NASA ADS] [CrossRef] [Google Scholar]
Krips, M., Martín, S., Eckart, A., et al. 2011, ApJ, 736, 37 [NASA ADS] [CrossRef] [Google Scholar]
Le Floc’h, E., Papovich, C., Dole, H., et al. 2005, ApJ, 632, 169 [NASA ADS] [CrossRef] [Google Scholar]
Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12 [NASA ADS] [CrossRef] [Google Scholar]
Lestrade, J.-F., Combes, F., Salomé, P., et al. 2010, A&A, 522, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Liu, D., Gao, Y., Isaak, K., et al. 2015, ApJ, 810, L14 [NASA ADS] [CrossRef] [Google Scholar]
Lu, N., Zhao, Y., Xu, C. K., et al. 2015, ApJ, 802, L11 [Google Scholar]
Lu, N., Zhao, Y., Díaz-Santos, T., et al. 2017, ApJS, 230, 1 [Google Scholar]
Lupu, R. E., Scott, K. S., Aguirre, J. E., et al. 2012, ApJ, 757, 135 [NASA ADS] [CrossRef] [Google Scholar]
Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
Magdis, G. E., Daddi, E., Elbaz, D., et al. 2011, ApJ, 740, L15 [NASA ADS] [CrossRef] [Google Scholar]
Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6 [NASA ADS] [CrossRef] [Google Scholar]
Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Markwardt, C. B. 2009, in Astronomical Data Analysis Software and Systems XVIII, eds. D. A. Bohlender, D. Durand, & P. Dowler, ASP Conf. Ser., 411, 251 [Google Scholar]
Meijerink, R., Kristensen, L. E., Weiß, A., et al. 2013, ApJ, 762, L16 [NASA ADS] [CrossRef] [Google Scholar]
Messias, H., Dye, S., Nagar, N., et al. 2014, A&A, 568, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Michałowski, M. J., Dunlop, J. S., Koprowski, M. P., et al. 2017, MNRAS, 469, 492 [NASA ADS] [CrossRef] [Google Scholar]
Murphy, E. J., Chary, R.-R., Dickinson, M., et al. 2011, ApJ, 732, 126 [NASA ADS] [CrossRef] [Google Scholar]
Narayanan, D., Krumholz, M. R., Ostriker, E. C., & Hernquist, L. 2012, MNRAS, 421, 3127 [NASA ADS] [CrossRef] [Google Scholar]
Nayyeri, H., Keele, M., Cooray, A., et al. 2016, ApJ, 823, 17 [NASA ADS] [CrossRef] [Google Scholar]
Negrello, M., Perrotta, F., González-Nuevo, J., et al. 2007, MNRAS, 377, 1557 [NASA ADS] [CrossRef] [Google Scholar]
Negrello, M., Hopwood, R., De Zotti, G., et al. 2010, Science, 330, 800 [NASA ADS] [CrossRef] [Google Scholar]
Negrello, M., Amber, S., Amvrosiadis, A., et al. 2017, MNRAS, 465, 3558 [NASA ADS] [CrossRef] [Google Scholar]
Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Omont, A., Neri, R., Cox, P., et al. 2011, A&A, 530, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Omont, A., Yang, C., Cox, P., et al. 2013, A&A, 551, A115 (O13) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Oteo, I., Zwaan, M. A., Ivison, R. J., Smail, I., & Biggs, A. D. 2016, ApJ, 822, 36 [NASA ADS] [CrossRef] [Google Scholar]
Oteo, I., Zhang, Z., Yang, C., et al. 2017, ApJ, 850, 170 [NASA ADS] [CrossRef] [Google Scholar]
Papadopoulos, P. P., & Seaquist, E. R. 1999, ApJ, 516, 114 [NASA ADS] [CrossRef] [Google Scholar]
Papadopoulos, P. P., & Greve, T. R. 2004, ApJ, 615, L29 [Google Scholar]
Papadopoulos, P. P., Zhang, Z.-Y., Xilouris, E. M., et al. 2014, ApJ, 788, 153 [NASA ADS] [CrossRef] [Google Scholar]
Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Rangwala, N., Maloney, P. R., Glenn, J., et al. 2011, ApJ, 743, 94 [NASA ADS] [CrossRef] [Google Scholar]
Riechers, D. A., Weiß, A., Walter, F., & Wagg, J. 2010, ApJ, 725, 1032 [NASA ADS] [CrossRef] [Google Scholar]
Riechers, D. A., Carilli, L. C., Walter, F., et al. 2011a, ApJ, 733, L11 [Google Scholar]
Riechers, D. A., Cooray, A., Omont, A., et al. 2011b, ApJ, 733, L12 [NASA ADS] [CrossRef] [Google Scholar]
Riechers, D. A., Hodge, J., Walter, F., Carilli, C. L., & Bertoldi, F. 2011c, ApJ, 739, L31 [Google Scholar]
Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Nature, 496, 329 [Google Scholar]
Rosenberg, M. J. F., van der Werf, P. P., Aalto, S., et al. 2015, ApJ, 801, 72 [NASA ADS] [CrossRef] [Google Scholar]
Rybak, M., McKean, J. P., Vegetti, S., Andreani, P., & White, S. D. M. 2015a, MNRAS, 451, L40 [NASA ADS] [CrossRef] [Google Scholar]
Rybak, M., Vegetti, S., McKean, J. P., Andreani, P., & White, S. D. M. 2015b, MNRAS, 453, L26 [NASA ADS] [CrossRef] [Google Scholar]
Saintonge, A., Kauffmann, G., Wang, J., et al. 2011, MNRAS, 415, 61 [NASA ADS] [CrossRef] [Google Scholar]
Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2 [NASA ADS] [CrossRef] [Google Scholar]
Salomé, P., Guélin, M., Downes, D., et al. 2012, A&A, 545, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
Santini, P., Maiolino, R., Magnelli, B., et al. 2010, A&A, 518, L154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19 [NASA ADS] [CrossRef] [Google Scholar]
Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Scott, K. S., Lupu, R. E., Aguirre, J. E., et al. 2011, ApJ, 733, 29 [NASA ADS] [CrossRef] [Google Scholar]
Scoville, N. Z., & Solomon, P. M. 1974, ApJ, 187, L67 [NASA ADS] [CrossRef] [Google Scholar]
Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84 [NASA ADS] [CrossRef] [Google Scholar]
Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150 [NASA ADS] [CrossRef] [Google Scholar]
Serjeant, S. 2012, MNRAS, 424, 2429 [NASA ADS] [CrossRef] [Google Scholar]
Silverman, B. W. 1986, Density Estimation for Statistics and Data Analysis (Chapman and Hall) [Google Scholar]
Simpson, J. M., Swinbank, A. M., Smail, I., et al. 2014, ApJ, 788, 125 [Google Scholar]
Smail, I., Ivison, R. J., & Blain, A. W. 1997, ApJ, 490, L5 [NASA ADS] [CrossRef] [Google Scholar]
Sobolev, V. V. 1960, Moving Envelopes of Stars (Cambridge: Harvard University Press) [Google Scholar]
Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677 [NASA ADS] [CrossRef] [Google Scholar]
Solomon, P. M., Downes, D., & Radford, S. J. E. 1992, ApJ, 387, L55 [NASA ADS] [CrossRef] [Google Scholar]
Solomon, P. M., Downes, D., Radford, S. J. E., & Barrett, J. W. 1997, ApJ, 478, 144 [NASA ADS] [CrossRef] [Google Scholar]
Spilker, J. S., Marrone, D. P., Aguirre, J. E., et al. 2014, ApJ, 785, 149 [NASA ADS] [CrossRef] [Google Scholar]
Spilker, J. S., Aravena, M., Marrone, D. P., et al. 2015, ApJ, 811, 124 [NASA ADS] [CrossRef] [Google Scholar]
Spilker, J. S., Marrone, D. P., Aravena, M., et al. 2016, ApJ, 826, 112 [NASA ADS] [CrossRef] [Google Scholar]
Spinoglio, L., Pereira-Santaella, M., Busquet, G., et al. 2012, ApJ, 758, 108 [NASA ADS] [CrossRef] [Google Scholar]
Swinbank, A. M., Smail, I., Longmore, S., et al. 2010, Nature, 464, 733 [Google Scholar]
Swinbank, A. M., Simpson, J. M., Smail, I., et al. 2014, MNRAS, 438, 1267 [Google Scholar]
Swinbank, A. M., Dye, S., Nightingale, J. W., et al. 2015, ApJ, 806, L17 [Google Scholar]
Tacconi, L. J., Neri, R., Chapman, S. C., et al. 2006, ApJ, 640, 228 [NASA ADS] [CrossRef] [Google Scholar]
Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781 [Google Scholar]
Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74 [Google Scholar]
Tan, Q., Daddi, E., Magdis, G., et al. 2014, A&A, 569, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge University Press) [Google Scholar]
Troncoso, P., Maiolino, R., Sommariva, V., et al. 2014, A&A, 563, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Tunnard, R., & Greve, T. R. 2016, ApJ, 819, 161 [NASA ADS] [CrossRef] [Google Scholar]
Valiante, E., Smith, M. W. L., Eales, S., et al. 2016, MNRAS, 462, 3146 [NASA ADS] [CrossRef] [Google Scholar]
van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
van der Werf, P. P., Isaak, K. G., Meijerink, R., et al. 2010, A&A, 518, L42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
van Dishoeck, E. F., Herbst, E., & Neufeld, D. A. 2013, Chem. Rev., 113, 9043 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Venemans, B. P., Walter, F., Zschaechner, L., et al. 2016, ApJ, 816, 37 [NASA ADS] [CrossRef] [Google Scholar]
Vieira, J. D., Crawford, T. M., Switzer, E. R., et al. 2010, ApJ, 719, 763 [NASA ADS] [CrossRef] [Google Scholar]
Vieira, J. D., Marrone, D. P., Chapman, S. C., et al. 2013, Nature, 495, 344 [NASA ADS] [CrossRef] [Google Scholar]
Wall, J. V., & Jenkins, C. R. 2012, Practical Statistics for Astronomers (Cambridge University Press) [Google Scholar]
Walter, F., Weiß, A., Downes, D., Decarli, R., & Henkel, C. 2011, ApJ, 730, 18 [Google Scholar]
Walter, F., Decarli, R., Carilli, C., et al. 2012, Nature, 486, 233 [Google Scholar]
Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44 [NASA ADS] [CrossRef] [Google Scholar]
Wardlow, J. L., Cooray, A., De Bernardis, F., et al. 2013, ApJ, 762, 59 [NASA ADS] [CrossRef] [Google Scholar]
Weiß, A., Henkel, C., Downes, D., & Walter, F. 2003, A&A, 409, L41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Weiß, A., Downes, D., Neri, R., et al. 2007, A&A, 467, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Weiß, A., De Breuck, C., Marrone, D. P., et al. 2013, ApJ, 767, 88 [NASA ADS] [CrossRef] [Google Scholar]
Wilson, C. D. 1995, ApJ, 448, L97 [NASA ADS] [CrossRef] [Google Scholar]
Wilson, D., Cooray, A., Nayyeri, H., et al. 2017, 848, 30 [Google Scholar]
Wong, T., & Blitz, L. 2002, ApJ, 569, 157 [NASA ADS] [CrossRef] [Google Scholar]
Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062 [NASA ADS] [CrossRef] [Google Scholar]
Yang, C., Gao, Y., Omont, A., et al. 2013, ApJ, 771, L24 (Y13) [NASA ADS] [CrossRef] [Google Scholar]
Yang, C., Omont, A., Beelen, A., et al. 2016, A&A, 595, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Zhang, Z.-Y., Henkel, C., Gao, Y., et al. 2014, A&A, 568, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Zhang, Z.-Y., Papadopoulos, P. P., Ivison, R. J., et al. 2016, R. Soc. Open Sci., 3, 160025 [Google Scholar]
Zhu, M., Seaquist, E. R., & Kuno, N. 2003, ApJ, 588, 243 [NASA ADS] [CrossRef] [Google Scholar]
Zhu, M., Gao, Y., Seaquist, E. R., & Dunne, L. 2007, AJ, 134, 118 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Spectral data

Fig. A.1

Spatially integrated spectra of CO lines in the H-ATLAS sources. The red lines represent the Gaussian fitting to the emission lines. The detections are ≳3σ. Zero velocity is set to the CO line sky frequency according to the previously measured spectroscopy redshifts z_spec given in Table 3. For CO(6–5) in G09v1.97, we also overlaid the ALMA spectral data (Yang et al., in prep.) in purple for comparison.

Fig. A.1

continued.

Fig. A.1

continued.

Appendix B: Additional table

Table B.1

Observed parameters of the CO and [C I] emission lines.

Appendix C: Bayesian approach of the radiative transfer modelling

Fig. C.1

Upper left panel: the CO SLED (without correcting for lensing magnification) of each source is plotted in black. The flux of CO(1–0) has been corrected for differential lensing effect as discussed in the text. The solid orange curve shows the best fit from the single component model corresponding to the maximum posterior possibility, while the solid purple line shows the best fit of the two-component model. Dashed red line shows the warmer component and the dashed-dotted line shows the cooler component fit. The upper limits are shown in grey open circles with downward arrows. Upper right panel: the posterior probability distributions of molecular gas density n_H₂, gas temperature T_kin and CO column density per velocity N_CO/ dv of the source, with the maximum posterior possibility point in the parameter space shown in orange lines and points. The contours are in steps of 0.5 σ starting from the centre. Lower panels: the posterior probability distributions for n_H₂, T_kin and N_CO/ dv of the cooler (dark green) and warmer (light green) component of the two-component model.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

continued.

All Tables

Table 1

Basic information on the CO rotational lines and [C I] ³P fine structure lines used in this paper.

	Fig. 8 Histograms of n_H₂ , T_k and N_CO /dv derived from single-component LVG modelling (black), and those of the low- and high-excitation components from two-component LVG modelling as shown in blue and red, respectively.
In the text

Molecular gas in the Herschel-selected strongly lensed submillimeter galaxies at z ~ 2–4 as probed by multi-J CO lines⋆,⋆⋆,⋆⋆⋆

1. Introduction

2. Sample, observations, and data reduction

2.1. Selection of the lensed SMGs

2.2. Observation and data reduction

3. Observation results

3.1. Observed CO line properties

3.2. Comparing our sample to the general SMG population

3.3. Intrinsic CO emission properties

4. Galactic properties and differential lensing

4.1. Molecular gas mass

4.2. Dynamical mass

4.3. Possible lensing biases

5. Physical properties of molecular gas

5.1. Multi-J CO line excitation and LVG modelling

5.2. The [C I](2–1) emission line in the high-redshift SMGs

5.3. Star formation and the molecular gas content in high-redshift GMGs

5.4. CO gas and H2O gas comparison

6. Summary

Acknowledgments

References

Appendix A: Spectral data

Appendix B: Additional table

Appendix C: Bayesian approach of the radiative transfer modelling

All Tables

All Figures

Molecular gas in the Herschel-selected strongly lensed submillimeter galaxies at z ~ 2–4 as probed by multi-J CO lines^⋆,^⋆⋆,^⋆⋆⋆

5.4. CO gas and H₂O gas comparison