M S S Bothwell
email: matthew.bothwell@gmail.com

J E E Aguirre

M Aravena

M Bethermin

T G G Bisbas

S C Chapman

C De Breuck

A H H Gonzalez

T R R Greve

Y Hezaveh

J Ma

M Malkan

D P Marrone

E J Murphy

J S Spilker

M Strandet

J D Vieira

A Weiß

ALMA observations of atomic carbon in z ∼ 4 dusty star-forming galaxies

Keywords: galaxies, high-redshift -galaxies, ISM -gravitational lensing, stronggalaxies, evolution -galaxies, formation -

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Understanding the properties and behaviour of the interstellar medium (ISM) of galaxies in the early universe is a cornerstone of modern galaxy evolution studies. Galaxies at early epochs show significantly elevated gas fractions relative to their local analogues (Tacconi et al. 2010), and it is these massive gas reservoirs that drive the enhanced star formation rates characteristic of the high-z Universe (i.e., Madau et al. 1996;Hopkins & Beacom 2006;[START_REF] Madau | [END_REF].

A variety of techniques have been used to observe the gas reservoirs in distant galaxies. Traditionally, the detection of intergalactic molecular gas (particularly at high-z) has relied on observations of various molecular emission lines of carbon monoxide ( 12 CO; Solomon & Vanden Bout 2005;Greve et al. 2005;Bothwell et al. 2013;Carilli & Walter 2013). CO exists in the centres of molecular clouds, andwith the aid of a CO-to-H2 conversion factor -observations of the CO line luminosity can be converted into a mass of molecular gas. This conversion is non-trivial however, with the abundance of CO relative to H2 varying as a function of ISM metallicity (increased metal abundance results in larger quantities of dust to protect CO from photodissociation) and even galactic behaviour (in merging systems the star-forming ISM no longer consists of discrete molecular clouds, making CO a more efficient tracer of molecular gas). Moreover, recent results have suggested that CO may also be destroyed by cosmic rays (produced indirectly by star formation, via supernovae, as well as AGN activity). The ISM of intensely star-forming galaxies, which produce a large cosmic ray flux, may be a hostile environment for CO molecules making it a less effective tracer of molecular gas than previously thought (Bisbas et al. 2015;Clark & Glover 2015).

Another method for tracing molecular gas is via the long-wavelength dust emission. Observations of the dust continuum can be converted into a total dust mass, which, with the aid of an assumed gas-to-dust ratio, can be used to calculate a gas mass (see Santini et al. 2010;Scoville et al. 2014;Scoville et al. 2016). This technique has some advantages over the use of CO lines: dust continuum observations are generally less time intensive, allowing for larger samples to be assembled. The dust mass method is not without its disadvantages, however. The dust temperature must be constrained (or assumed) in order for a dust mass to be measured. In addition, the gas-to-dust ratio remains a relatively poorly studied quantity, which may vary by up to a factor of ∼ 20 in bright star-forming galaxies (Zavala et al. 2015), as well as potentially varying at high redshift (Dwek et al. 2014;Micha lowski 2015). Furthermore, observations of the dust continuum provide no kinematic information, which is available when observing CO emission lines.

In recent years, the emission line of atomic carbon ([Ci]( 3 P1 → 3 P0); [Ci](1 -0) hereafter) has been found to be an excellent alternate tracer of the cold molecular ISM, being closely associated with low-J CO emission across a wide range of environments and redshifts (Papadopoulos et al. 2004;Walter et al. 2011;Alaghband-Zadeh et al. 2013;Israel et al. 2015). This conclusion is supported by both detailed studies of nearby Galactic molecular clouds (in which CO and [Ci] are found to co-exist throughout the bulk of the cold molecular component; Papadopoulos et al. 2004), as well as hydrodynamic simulations (Tomassetti et al. 2014).

Using [Ci] as a molecular gas tracer offers a number of advantages, compared to observations of both CO and the dust continuum. Due to the simplicity of the quantum fine structure level of [Ci], many physical parameters (including excitation temperature and total carbon mass) can be calculated with minimal uncertainty. And while CO becomes an increasingly poor tracer of molecular gas as metallicity decreases (due to a lack of dust grains being available to shield CO molecules from photodissociation), [Ci] is affected far less severely -significantly reducing the uncertainty on the molecular gas mass introduced by unknown metal abundances. Furthermore, Bisbas et al. (2015) found that a high cosmic-ray ionisation rate will destroy CO molecules, dissociating them into [Ci], while leaving the underlying H2 unaffected. Additionally, observations of the [Ci] line provide the same valuable kinematic information as CO, offering a distinct advantage over dust-based methods. As a result [Ci] can be a powerful and effective tracer of the molecular ISM in distant galaxies.

In this work, we present observations of the [Ci](1 -0) emission line in a sample of strongly lensed Dusty Star Forming Galaxies (DSFGs) which were identified via the South Pole Telescope (SPT; Carlstrom et al. 2011) wide-field survey (Vieira et al. 2010). DSFGs are extremely luminous starforming galaxies (typical SFRs ∼ 1000 M /yr), which are thought to be the high-z progenitors of the most massive galaxies in the z ∼ 0 Universe. The SPT has proven to be an efficient machine for finding the brightest (stronglylensed) DSFGs in the Universe (Vieira et al. 2013;Hezaveh et al. 2013). The DSFGs in this work were all taken from the 26-galaxy sample targeted for spectroscopic redshift identification by the Atacama Large Millimeter Array (ALMA) in Cycle 0 (Weiß et al. 2013), which confirmed redshifts via the identification of a number of atomic and molecular emission lines. Thirteen DSFGs were found to lie at redshifts 3.24 < z < 4.85, shifting the [Ci](1 -0) into ALMA Band 3. In this work we present an analysis of the [Ci] properties of these 13 sources.

This paper falls broadly into two halves. In the first half, we present an analysis of the [Ci] properties of our sample. We describe the sample and the various ancillary data in §2, and in §3 we present our analysis of the [Ci](1 -0) data (including calculations of the mass and cooling contribution of atomic carbon, analysis of the kinematic properties of the sample, and discussions of the use of [Ci](1 -0) as a tracer of molecular gas and star formation mode). Moving to the second half of the paper, in §4 we combine the [Ci](1 -0) line with a variety of atomic and molecular emission lines in order to constrain the physical conditions in the ISM of our galaxies using a modern Photodissociation Region (PDR) modelling code, 3d-pdr. We present a discussion of our results in §5, and we present our conclusions in §6. Throughout this work we adopt a standard Λ-CDM cosmology with parameters taken from Planck Collaboration et al. ( 2016); h = 0.678, Ωm = 0.308, and ΩΛ = 0.692.

Carbon in the ISM of lensed dusty galaxies at z ∼ 4 3

SAMPLE, OBSERVATIONS AND REDUCTION

The 13 sources presented in this work were originally targeted as part of our ALMA blind redshift search programme (Weiß et al. 2013), in which 26 SPT DSFGs were observed across the entirety of ALMA Band 3 (= 84 -116 GHz) as part of the Cycle 0 'early compact array' setup. Each observation consisted of 5 distinct 7.5 GHz tunings, spaced to cover the band. Each source was observed for 120s in each tuning configuration. Further details of the observing programme and data reduction can be found in Weiß et al. (2013). 13 of the 26 DSFGs lie in the redshift range 3.24 < z < 4.85, causing the [Ci](1 -0) emission line (νrest = 492.161 GHz) to be redshifted into Band 3. Of these 13 sources, just one (SPT0345-47) was not detected in [Ci](1 -0). A further two sources (SPT0300-46 and SPT2103-60) are tentatively detected at the ∼ 3σ level (i.e., 3 times the rms channel noise, which we measure to be 2.1 mJy and 2.6 mJy per 50 km/s channel respectively). The remaining 10 sources all show clearly detected [Ci](1 -0) emission at > 4σ significance.

For SPT0345-47, we calculate a 3σ upper limit on the intensity of the [Ci](1 -0) emission line:

ICO < 3 RMS channel √ ∆VCO dv , (1)
where RMS channel is the RMS channel noise in the spectrum of SPT0345-47 (which we measure to be 2.4 mJy per 50 km/s channel), ∆VCO is the mean linewidth of the detected sample (=410 km/s), and dv is the bin size in km s -1 (=50 km/s). We calculate an upper limit on the [Ci] line intensity for SPT0345-47 of I [CI](1-0) < 1.03 Jy km/s. Figure 1 shows spectral cutouts at the position of the [Ci](1 -0) line for the 13 SPT DSFGs analysed in this work. For reference, we have overlaid (where available) the CO(2 -1) emission line (scaled arbitrarily in flux for ease of comparison).

Due to a lack of X-ray data for our sample we cannot rule out a possible AGN component in any of our 13 DSFGs. However it is unlikely that any of our sources contains a significant AGN. DSFGs as a class are star-formation dominated (i.e., the bolometric output of the galaxy originates predominantly from young massive stars, rather than accretion onto a central compact object); even DSFGs with some measurable AGN activity tend to be primarily starformation-driven objects (Alexander et al. 2005). Furthermore, Chandra X-ray observations of the most compact and IR-luminous SPT-DSFG, SPT0346-52, found no sign of AGN activity (Ma et al., 2016). We proceed with the assumption that our the DSFGs in our sample are starformation dominated objects.

All 13 galaxies presented in this work are strongly gravitationally lensed, and have detailed lens models based on ALMA 870µm observations, which allow their lensing magnifications to be calculated (Spilker et al. 2016). We discuss the use of these lens models to remove the effects of gravitational lensing in §3.1 below.

Ancillary data

The SPT-DSFG sample has been the target of several followup programmes designed to survey a number of ISM diagnostics. In particular, in addition to the [Ci](1 -0) emission lines presented in this work, 9/13 sources have observations covering the CO(2 -1) emission line, 10/13 have observations covering the CO(4 -3) emission line (the remaining 3 have observed CO(5 -4) emission, which can be converted to (4 -3) with the aid of an assumed CO Spectral Line Energy Distribution (Bothwell et al. 2013;Spilker et al. 2014)). In addition, 10/13 have the [Cii] emission line observed.

The mid-J lines of CO come from the ALMA spectra used in this work. The [Cii] observations were obtained from the APEX First Light APEX Submillimetre Heterodyne receiver (FLASH), and the Herschel 'SPIRE' FTS spectrometer. Observations were made at 345 GHz (for sources with redshifts 4.2 < z < 5.7) and 460GHz (3.1 < z < 3.8), with system temperatures of 230 K and 170 K respectively. Further details of the observations and data reduction can be found in Gullberg et al. (2015).

The low-J CO observations were taken by the Australia Telescope Compact Array (ATCA), as part of a targeted program aimed at obtaining CO(1 -0) or CO(2 -1) for SPT DSFGs with secure redshifts. ATCA was used in its H214 hybrid array configuration, with the Compact Array Broadband Backend in wide bandwidth mode. The mean rms noise for the CO(2 -1) observations used in this work was 0.5 mJy per 50 km/s channel. Further details of the observations and data reduction can be found in Aravena et al. (2016).

This suite of molecular emission lines enables a more detailed treatment of the conditions of the ISM, as multiple line ratios can be used to independently constrain various parameters of interest.

RESULTS AND ANALYSIS

Line luminosities and ratios

Throughout this work, we calculate emission line luminosities using the formulae below. Luminosities in solar units (L ), which represent the true energy output carried by the emission line (used, for example, for calculating cooling contributions) are calculated following Solomon & Vanden Bout (2005):

L line = 1.04 × 10 -3 S line ∆v νrest (1 + z) -1 D 2 L , (2)
where S line ∆v is the velocity-integrated line flux in Jy km s -1 , νrest is the rest frequency in GHz, and DL is the luminosity distance in Mpc. Alternatively, line luminosities in units of K km s -1 pc 2 are calculated using

L line = 3.25 × 10 7 S line ∆v ν -2 obs (1 + z) -3 D 2 L , (3)
which gives line luminosities proportional to brightness temperature. Table 1 lists line luminosities (in solar units) for the DSFGs used in this work. In addition, we also list in Table 4 Bothwell et al. 1 the bolometric far-IR luminosity, calculated using greybody SED fits to far-IR/millimetre photometry (Greve et al. 2012). Note that all values in Table 1 are observed quantities, which have not been corrected for the effects of gravitational lensing.

ID RA DEC z L [CI](1-0) L CO(2-1) L CO(4-3) L CO(5-4) L [CII] L FIR L [CI] /L FIR µ [J2000] [J2000] [

Correction for magnification due to gravitational lensing

The DSFGs discovered by the SPT are typically strongly lensed by an intervening massive galaxy (Hezaveh et al. 2013;Vieira et al. 2013;Spilker et al. 2016). In the most simple form, strong gravitational lensing both distorts the lensed source and boosts its apparent luminosity by a magnification factor (µ), which is dependent on both the mass of the intervening lens and the source/lens configuration. The effects of gravitational lensing need to be measured and accounted for in order to study the intrinsic properties of the source.

The absolute magnification needs to be corrected for if we are to discuss any innate source properties (such as the total molecular gas mass). Lens modelling, carried out based on our ∼ 0.5 870µm ALMA imaging, has been presented by Hezaveh et al. (2013) and Spilker et al. (2016). All of the sources in this work have their magnification factors measured, spanning a range 3.6 < µ < 27 (with a sample mean and standard deviation of 13.1±9.3). We note that these values have been calculated based on the mm-wavelength dust emission -we make the assumption that these lens models also apply to the cold molecular gas traced by [Ci].

In addition, the presence of any inhomogeneity in the source can potentially result in differential lensing, by which some regions of the source lying close to a caustic are magnified by a disproportionate amount (Hezaveh et al. 2012). This can distort the apparent source properties: a differentially lensed region with a higher than average temperature will cause the source as a whole to appear hotter. In general, differential lensing tends to selectively apply a magnification boost to compact regions, relative to more extended components (though this is not always true; see Hezaveh et al. 2012, who demonstrate a source-lens geometry which boosts an extended component relative to a more compact component).

In this work we wish to use [Ci] to trace molecular gas, as well as studying ratios of various emission lines emitted by our DSFGs. If the [Ci] emission is not conterminous with the underlying H2, or if some emission components are systematically more compact/extended than others, the effect of differential gravitational lensing could be (depending on the source-lens geometry) to distort our results.

Overall, it is likely that differential lensing does not significantly affect the ability of [Ci] to trace molecular gas. Studies of [Ci] in the local Universe reveal a [Ci] distribution which is throughly mixed with the molecular ISM as a whole (see [START_REF] Keene | IAU Symposium[END_REF]Ikeda et al. 2002;Ojha et al. 2001;Papadopoulos et al. 2004). As such, the [Ci] and underlying molecular gas likely have similar surface brightness profiles, and as such their ratios will not be changed by differential lensing. Serjeant (2012) discusses the effect of differential lensing on line ratios in sub-mm selected sources. They find that while some observational properties of DSFGs are affected by differential lensing, many ratios remain robust enough to allow physical interpretations to be drawn. One such robust parameter is the ratio between low-J CO and FIR luminosity. Given that both low-J CO and [Ci] are effective tracers of the cold molecular gas component (see §3.5), it is likely that the ratio between [Ci] and FIR luminosity is similarly robust.

Is is possible that the other line ratios in this work are susceptible to differential lensing effects -as calculated by Serjeant (2012), the ratio between high-J and low-J CO may have an added uncertainty of ∼ 30% purely due to differential lensing. This effect may affect the ratio L [CI](1-0) /L CO(4-3) analysed in §3.3, and we discuss this possibility in that section. Additionally, some ratios used as input to PDR models in §4 may be subject to differential lensing effects -in the case of the PDR models, the potential added uncertainty estimated by Serjeant (2012) is far smaller than the underlying model uncertainties.

Fig. 2 shows the intrinsic (i.e., de-lensed) [Ci] luminosities and far-IR luminosities for our sample. For comparison, we have also plotted a sample of DSFGs and AGN observed in [Ci] taken from Alaghband-Zadeh et al. (2013). The mean (±SD) [Ci] luminosity for the literature sample of DSFGs is (9.9 ± 3.9) × 10 9 K km s -1 pc 2 , while the same quantities for our sample of SPT-DSFGs (1.1 ± 0.8) × 10 10 K km s -1 pc 2 . After correction for gravitation magnification, it is clear that SPT-DSFGs have [Ci] luminosities comparable to other DSFGs in the literature.

The [Ci]-CO and [Ci]-L(FIR) line ratio

It is possible to use a combination of line ratios as probes of the conditions within the ISM of our sources. Here we use a combination of line ratios in order to compare the conditions within our DSFGs with both the DSFGs presented by Alaghband-Zadeh et al. (2013), and a sample of local (z < 0.05) galaxies. For our local sample, we use a sample observed with the Herschel/SPIRE Fourier Transform Spectrometer (FTS), presented by Kamenetzky et al. (2016). The sample as presented is a compilation of all extragalactic proposals listed in the Herschel Science Archive (HSA) with at least one reported FTS line measurement or upper limit. We have further selected galaxies with available 8 Bothwell et al. 2013). For comparison, we have also included a sample of local galaxies observed with Herschel/SPIRE, presented by Kamenetzky et al. (2016). The ratio L [CI](1-0) /L FIR is an approximate (but non-linear) tracer of the UV field strength, while the ratio L [CI](1-0) /L CO(4-3) is an approximate (but again non-linear) tracer of the gas density.

[Ci] and CO(4 -3) line fluxes, in order to compare to our high-z DSFGs.

We consider the ratio L [CI](1-0) /LFIR, which is a tracer of the strength of the interstellar UV radiation field, and the ratio L [CI](1-0) /L CO(4-3) , which is a tracer of the average gas density (Kaufman et al. 1999; see also Alaghband-Zadeh et al. 2013). It must be noted that both of these ratios are non-linear with the physical conditions they trace: a more rigorous treatment of line ratios tracing physical conditions is presented in §4 below. We have converted the 40 -120µm FIR luminosities given by Kamenetzky et al. (2016) into 8 -1000µm luminosities (to match our DSFG samples) by multiplying by a factor of 1.9, following Eq. 4 in Elbaz et al. (2002). We have also restricted the local sample to galaxies with redshifts (z > 0.05), to ensure that the ∼ 40 FWHM SPIRE beam covers physical scales of > 4kpc, therefore capturing flux beyond the galaxy centres (which are preferentially dense).

Figure 3 plots the ratio L [CI](1-0) /LFIR against the ratio L [CI](1-0) /L CO(4-3) for the three samples. We indicate on the plot the way the physical conditions (UV field strength, gas density) vary with these ratios. We firstly note that these parameters are correlated for our combined sample, with the sources exhibiting the highest densities also having the strongest interstellar UV radiation fields (and vice-versa). We also note that our SPT-DSFG sample is skewed towards the 'upper end' of the distribution, representing galaxies with the strongest UV fields and the densest gas. 66% of the SPT-DSFG sample -8/12 -have L 'significance threshold'. Any difference between the distribution of the samples must therefore be regarded as tentative. We discuss the ISM density of our sample further in §5.1.

Mass and cooling contribution of atomic carbon

The luminosity of the [Ci] line, in solar units (Eq. 1) gives cooling contribution of [Ci] (i.e., the amount of energy radiated away by the line). We calculate for our sample the ratio between the cooling contribution of [Ci] and the total FIR luminosity, L [CI](1-0) /LFIR. These ratios are listed in Table 1. Note that this ratio is unaffected by lensing magnification (in the absence of differential lensing, which is likely to be negligible for this particular ratio). 2013) sample (for the AZ13 sample, the FWHM of the CO(3 -2) line is shown). The solid and dashed lines respectively show a ratio of unity and a factor of 2 variation. There is a close kinematic correspondence between the [Ci] and low-J CO lines -in the majority of cases being equivalent within the uncertainties. There is also no systematic difference between the two (< FWHM CI /FWHM CO >= 1.03 ± 0.40).

MCI = 5.706 × 10 -4 Q(Tex) 1 3 e (23.6/Tex) L CI(1-0) , (4)
Q(Tex) = 1 + 3e -T 1 /Tex + 5e -T 2 /Tex , (5)
and T1 = 23.6K, T2 = 62.5K are the excitation energy levels of atomic carbon.

As our observations only cover the [Ci](1 -0) emission line, we cannot directly calculate the excitation temperature Tex (which would also require the [Ci](2 -1) line). Instead, we adopt a 'typical' value of 30K (Weiß et al. 2013) We note that for a wide range of Tex the derived carbon masses depends very weakly on the assumed value of temperature (as shown in Fig. 4).

We calculate a mean observed [Ci] mass (corrected for lensing magnification as discussed in §3.2 above) of (1.2 ± 0.3) × 10 7 M ; individual atomic carbon masses for the DSFGs in our sample are given in Table 2.

[Ci] as a tracer of the total gas mass

Many authors have pointed out that [Ci](1 -0) emission is a good tracer of the bulk of the cold ISM, and therefore makes an excellent proxy for the (unobservable) H2 mass. If this is indeed the case, then the [Ci] line emission should be emitted primarily by the cool, extended gas component traditionally traced by low-J CO emission. In order to test whether this is likely the case, we compare the linewidths of the The results are shown in Fig. 5. The majority of our combined sample have low-J CO and [Ci] linewidths consistent with each other, given the observational and fitting uncertainties on each. Importantly, there is no systematic difference between the FWHMs of the low-J CO and [Ci] lines -the mean ratio between the two, for the combined sample, is < FWHMCI/FWHMCO >= 1.03 (with a standard deviation of 0.40). A Kolmogorov-Smirnov test comparing the CO and [Ci] linewidths returns P = 0.88, suggesting that the two are consistent with each other. This kinematic correspondence between the two lines suggests that the [Ci] emission is tracing the same gas component as the low-J CO emission -a line which is thought to be emitted entirely from the cold reservoir of molecular gas in the ISM. We proceed with the assumption that the [Ci] line is an effective tracer of the molecular gas reservoir in our DSFGs -though this assumption may break down in certain situations (i.e., in very dense environments, the [Ci] line can become optically thick and therefore a less effective tracer of gas). Papadopoulos & Greve (2004) give an expression for calculating the total H2 mass from the luminosity of the

ID M([Ci]) M(H 2 ) [CI] α CO [×10 7 M ] [×10 10 M ] [K km s -1 pc 2 ) -1 ] SPT0113-
[Ci](1 -0) line: M(H2) [CI] = 1375.8 D 2 L (1+z) -1 X [CI] 10 -5 -1 A10 10 -7 s -1 -1 × Q -1
10 S [CI] ∆v, (6) where XCI is the [Ci]/H2 abundance ratio. This does not include a contribution from helium. Here, following Papadopoulos & Greve (2004), we adopt a literature-standard [Ci]/H2 abundance ratio of 3 × 10 -5 , and the Einstein A coefficient A10 = 7.93 × 10 -8 s -1 . Q10 is the excitation factor, which we take to be 0.6. The value of Q10 is dependent on the specific conditions within the gas. Local ULIRGs have typical measured Q10 ∼ 0.5 -we have chosen Q10 = 0.6 to ensure consistency between Eq. 5, and Eq. 4. Using Eq. 5, we measure a mean (corrected for gravitational magnification) M(H2) CI = (6.6 ± 2.1) × 10 10 M . We note that this value is dependent on an assumption of a [Ci]/H2 abundance ratio (analogous to the αCO, the CO-to-H2 conversion factor used to derive gas masses from 12 CO luminosities). 2013) with both [Ci] and CO observations. To derive CO-based H 2 masses, we have assumed α CO = 0.8. The dashed diagonal line indicates equal masses, while the dotted lines show a factor of two variation each way. DSFGs for which a higher-J CO line has been used to infer a gas mass have been plotted using darker colours. Right Panel: Ratio between [Ci]-derived molecular gas mass and CO(1 -0) luminosity, for the same galaxies. This ratio can be interpreted as giving the CO-to-H 2 conversion factor α CO . All (but one) SMGs in both samples have implied values of α CO far above the canonical 'ULIRG' value, and several SPT-DSFGs have an α CO comparable to that of the Milky Way.

In recent years, several authors have developed models aiming to examine the behaviour of the [Ci] abundance. Both Offner et al. (2014) and Glover & Clark (2016) present post-processed hydrodynamical simulations of star-forming clouds, finding that the [Ci]/H2 ratio varies as a function of a range of galactic parameters, including H2 column density, the strength of the interstellar radiation field (ISRF), and metallicity. Glover & Clark (2016) find that increasing the ISRF by factors of 10 2 -10 3 raises the [Ci] abundance at low AVs by 30-50% (and probably has the same effect on the CO-to-H2 conversion factor). At high AVs, the [Ci] abundance is raised by cosmic rays. Likewise, Glover & Clark (2016) present evidence that the [Ci] abundance increases as metallicity decreases, with a scaling ∝ Z -1 . There are also indications that dense, star-forming environments will show elevated vales of XCI -Papadopoulos & Greve 2004 report a 'typical' value of XCI = 3 × 10 -5 , but an elevated value of XCI = 5 × 10 -5 in the centre of the local starburst M82. To derive gas masses here, we have taken the 'standard' value of 3 × 10 -5 , but in §5.3 below we discuss the possibility of variation in the [Ci] abundance.

Comparing to CO-based gas masses

We now compare our [Ci]-derived molecular gas masses to measurements using a more common tracer, the luminosity of 12 CO (which is converted to a molecular gas mass via the CO-to-H2 conversion factor αCO). The advantage of performing this comparison lies in the fact that the 'conversion factor' required to convert a [Ci] flux into a molecular gas mass is potentially less uncertain than the CO-to-H2 conversion factor (Papadopoulos & Greve 2004;Papadopoulos et al. 2004), being only linearly dependent on the metallicity of the gas. This is opposed to the more commonly used COto-H2 conversion factor, which depends roughly quadratically on metallicity -see Bolatto et al. ( 2013) for a recent discussion of issues surrounding the CO-to-H2 conversion factor. As such, our [Ci]-based H2 masses may be used to estimate the value of the CO-to-H2 conversion factor. This, in turn, gives insight into the conditions in the ISM of the galaxies in question -low, 'ULIRG'-like values of the conversion factor imply a dense (possibly merger-compressed) ISM, while higher values of αCO imply a more extended molecular phase.

We calculate our CO-based gas masses using the standard equation

M(H 2 ) CO = αCO L CO(1-0) , (7)
where αCO is the CO-to-H2 conversion factor in units of M (K km s -1 pc 2 ) -1 . 1

The CO observations of the sample of DSFGs included in this work do not include the ground-state (1 -0) line required to 'directly' derive a gas mass. Instead, we convert our higher-J line luminosities down into an equivalent (1-0) luminosity by assuming a typical DSFG SLED (Bothwell et al. 2013;Spilker et al. 2014). The majority (9/13) of our sample have the CO(2-1) emission line observed, which can be easily converted into CO(1-0) with minimal uncertainty (as CO(2 -1) is also an excellent tracer of the total cold molecular gas). We assume a CO(2-1)/CO(1-0) brightness temperature ratio of r 21/10 = 0.8 (Aravena et al., 2016). The remaining 4 sources have only higher-J lines -either CO(4 -3) or CO(5 -4) -which can still be extrapolated down to CO(1 -0), albeit with greater uncertainty.

Fig. 6 (left) shows a comparison of H2 masses derived using [Ci] (and XCI = 3 × 10 -5 ), with those derived using CO (and αCO = 0.8). The SPT DSFGs with only higher-J CO lines observed (and therefore with more uncertain CO-1 We hereafter omit the units of α CO in the interest of brevity.

Carbon in the ISM of lensed dusty galaxies at z ∼ 4 11 based H2 masses) are plotted with darker colours. We have also plotted DSFGs and AGN from Alaghband-Zadeh et al. ( 2013) with both CO and [Ci] observations. As Fig. 6 shows, molecular gas masses derived via CO emission (assuming αCO = 0.8) are systematically lower than those derived using [Ci](assuming XCI = 3×10 -5 ). For SPT DSFGs, we find <M(H2) CI,X=3e-5 >= (6.6±2.1)×10 10 M , and <M(H2) CO,α=0.8 >= (2.8 ± 1.7) × 10 10 M . Only a single DSFG, taken from Alaghband-Zadeh et al. ( 2013), has a CO-based molecular gas mass in excess of the value measured using [Ci]. If [Ci] does indeed have higher accuracy as a molecular gas tracer (as discussed above, it has a reduced metallicity dependence relative to CO), it seems that calculating gas masses from CO, and adopting αCO = 0.8, is underestimating the molecular masses of our DSFGs.

It is possible to invert this problem: we can compare M(H2) CI to L CO(1-0) directly, in order to estimate the value of αCO needed to bring the two molecular gas mass measurements into agreement. Fig. 6 (right) shows the ratio of M(H2) CI to L CO(1-0) , plotted against the observed FIR luminosity. The units of L CO result in the ratio M(H2) CI /L CO(1-0) being equal to the value of αCO required to force agreement.

It can be seen that in all cases but one, SPT-DSFGs and literature DSFGs have implied αCO values far in excess of the value generally adopted for both local ULIRGs and high-z DSFGs (typically, αCO = 0.8. For SPT-DSFGs, we find a mean CO-to-H2 conversion factor of αCO = 2.0 ± 1.0. Calculating this value for just those DSFGs with low-J CO mass measurements (i.e., excluding the three DSFGs only observed in Jup 4, and therefore with more uncertain gas masses), we find a slightly higher value, αCO = 2.4 ± 0.8. Such a high value of αCO would imply that the ISM is likely not tidally-compressed as a result of a merger, but exists in a more evenly distributed form. This model has some theoretical support - Narayanan et al. (2015) use hydrodynamical simulations to model a luminous DSFG powered entirely by gas inflow, with no merger needed.

The fact that our [Ci] observations suggest values of αCO a factor of ∼ 2 -3 times higher than are generally assumed for starburst-mode galaxies will have the effect of increasing the derived gas fractions of our objects. DSFGs are known to number amongst the most gas-rich systems in the Universe: Bothwell et al. (2013) found a mean baryonic gas fraction -defined as fgas = Mgas/(Mgas + M * ) -for DSFGs of fgas = 0.43 ± 0.05, using a CO-to-H2 conversion factor of αCO = 1. Unfortunately, only two of the galaxies in our sample have the measured stellar masses required to calculate a gas fraction (both taken from Ma et al. 2015); SPT2146-55 (M * = 0.8 +1.9 -0.6 × 10 11 M ) and SPT2147-50 (M * = 2.0 +1.8 -0.9 × 10 10 M ). Using standard αCO = 0.8 gas masses, we calculate gas fractions for these two galaxies of ∼ 20% and ∼ 70%, respectively. However, their [Ci]-based gas masses suggest higher gas fractions, of ∼ 60% and ∼ 80% (respectively).

We can also calculate a 'typical' gas fraction for SPT DSFGs, by comparing the mean gas mass with the mean stellar mass (keeping in mind that these values were calculated mostly for different individual galaxies, and this applies to the sample as a whole only in an average sense). Using αCO = 0.8 gas masses, we find a 'sample average' gas fraction of ∼ 40% (in agreement with the larger un-lensed sample presented by Bothwell et al. 2013). Adopting [Ci]based gas masses, we find a typical gas fraction of ∼ 60%, significantly higher than some previous CO-based estimates. These elevated gas fractions are not a unique feature of [Ci] observations -using dust-based gas masses, Scoville et al. (2016) find gas fractions of 50 -80% for the most massive, high-sSFR galaxies at z > 2.

These high gas fractions are all reliant on our assumption of a [Ci] abundance of XCI = 3 × 10 -5 . In §5.3 below we discuss the effect of challenging this assumption.

Dynamical mass and potential lensing bias

In their sample of 40 DSFGs, Bothwell et al. (2013) find a close correlation between the luminosity of 12 CO, L CO , and the full width half max (FWHM) of the CO line -a trend which Bothwell et al. (2013) ascribe to the baryondominated dynamics in the central few kpc of their DSFGs. Here we perform a similar analysis, comparing the total gas mass with the FWHM of the observed line for our sample of strongly-lensed DSFGs.

Figure 7 shows M(H2) (derived using [Ci]), plotted against the FWHM of the [Ci] line, for our SPT-DSFGs. We also plot the same quantities for 6 DSFGs taken from Alaghband-Zadeh et al. (2013). In order to compare to the wider DSFG population, we also plot CO-derived gas masses and CO linewidths for un-lensed DSFGs taken from Bothwell et al. (2013). For this latter sample, we have only used sources observed in low-J CO lines (either 2 -1 or 3 -2), as at Jup 4 CO lines become increasingly poor dynamical tracers of the total gas reservoir.

We have over-plotted a simple dynamical mass model:

M(dyn) = v 2 rot r G , (8)
where M(dyn) is the total dynamical mass, vrot is the rotation velocity of the galaxy, r is the radius, and G is the gravitational constant. To convert from our observed line FWHMs into intrinsic (dynamical mass-tracing) velocities is challenging, and requires knowledge of the geometry, kinematics, and velocity anisotropy of the galaxy. To account for these unknowns, Erb et al. ( 2006) describe a dimensionless constant, C, which ranges from C < 1 for flat rotating disks, to C > 5 for virialised spherical systems. For illustrative purposes we adopt C = 4, a value between a pure disk and a purely virialised spherical merger (see §5.2).

Taking a typical gas fraction fgas = 0.5, and assuming that the dynamics in the regions of the galaxies we observe are baryon-dominated (i.e., we assume M(H2) = M(dyn)/2), we can write an expression for M(H2) in terms of line FWHM;

M(H2) = C (FWHM/2.35) 2 r 2G (9)
The remaining unknown is the radius of the gas disk, r. Following the source size analysis of the SPT-DSFGs sample presented by Spilker et al. (2016), we take a typical radius of 2kpc. The relation in Eq. 8, taking r = 2kpc and C = 4, is overplotted on Figure 7.

Figure 7 shows that the combined sample of DSFGs (both the lensed objects observed in [Ci] and the unlensed objects observed in CO) are reasonably well described by the simple dynamical mass model in Eq. 8, following the trend M(H2) ∝ FWHM 2 . The scatter around this relation can be primarily attributed to two potential physical causes: size variation, and differing kinematics and inclinations (that is, the extent to which the observed FWHM is effectively tracing the true rotational velocity of the galaxy).

As pointed out by Hezaveh et al. (2012), lensing-selected samples of galaxies can display a size bias relative to the general un-lensed population. Members of a flux-selected lensed sample will be biased towards being systematically compact (as a compact source lying close to a lensing caustic will be magnified more than a similarly-positioned extended source). Despite this, however, based on sizes derived from lens model fits to ALMA 870µm data Spilker et al. (2016) find that the observed angular 870µm size distribution of SPT-DSFGs is statistically consistent with the distribution of 870µm angular sizes displayed by members of un-lensed comparison samples. With the higher mean redshift of the SPT-DSFG sample, and the redshift evolution of the angular scale, this implies that DSFGs in the higher-redshift SPT sample are slightly more physically (i.e., 1 at the mean redshift of SPT-DSFGs, z = 3.5, corresponds to 7.47 kpc, while 1 at the mean redshift of unlicensed DSFG samples, z = 2.2, corresponds to 8.42 kpc). This potential size bias is smaller than the uncertainties, suggesting that the dispersion in 7 is due to kinematic and inclination effects.

As Fig. 7 shows, though the dispersion around the dynamical mass relation is large, our sample of SPT-DSFGs are not systematically offset in the FWHM vs. M(H2) plane relative to the general population of un-lensed DSFGs. If any systematic offset were present, it would suggest that SPT-DSFGs have distinctly different kinematics relative to unlensed DSFGs (e.g., having predominantly virialised rather than predominantly rotational kinematics). The fact that our sample of lensed DSFGs shows no such offset strongly suggests that this bias has been introduced as a result of our lensing selection. That is, the molecular gas reservoirs in our sample are kinematically consistent with those in the underlying DSFGs population as a whole.

PDR MODELLING

The large number of lines observed in our sample of DSFGs (including multiple transitions of 12 CO, as well as atomic and ionised carbon species) allow for the use of models to help constrain the conditions in the ISM regions emitting the lines. One such class of models are 'PDR' (Photodissociation Region) models, which model galactic regions where photons are the dominant driver of the interstellar heating and/or chemistry.

PDR models are invoked in order to use line intensities (and ratios of intensities) to constrain conditions within the ISM -specifically, gas density (n cm -3 , the volume density of hydrogen gas) and UV-field strength (normally expressed in terms of G, the Milky Way UV field strength in Habing units, 1.6 × 10 3 erg s -1 cm -2 ). Certain ratios are good tracers of either the density or the UV field strength, and by combining a number of lines and ratios, the density and UV field strength can be constrained simultaneously, giving a window into the typical ISM conditions emitting the lines in question. PDR models (i.e. Bisbas et al. 2014) show that fine structure lines (e.g., the [Cii] 158µm and [Ci] lines we analyse in this work) are primarily emitted at low AVs (their peak of local emissivity is always AV <7, with the [Cii]-158µm line emitted from regions with AV < 1 -2 mag). Therefore the [Cii]/[Ci](1 -0) ratio is a good tracer of the PDR conditions at low AVs, such as the UV radiation strength (assuming that both lines are emitted from the same region and they suffer from the same beam dilution effects). Low-J transitions of CO are emitted primarily from regions with high AV, providing information about the cold/molecular gas of each object. As the J transition increases, the peak of local emissivity occurs at lower and lower AV, which can then provide information about the state of PDRs in these conditions.

Previous work presenting PDR analyses of molecular and atomic emission lines in DSFGs (e.g., Danielson et al. 2011;Alaghband-Zadeh et al. 2013) have generally used the well-known Kaufman et al. (1999) PDR models. There are, however, a number of critical parameters that affect the PDR modelling results, which the Kaufman et al. (1999) models do not account for. Cosmic rays, for example, can also contribute to the chemistry and heating, but in normal star-forming galaxies will only have a non-negligible effect at high optical depths, in the centres of molecular starforming cores. Cosmic rays in galaxies are produced by both AGN and star formation (cosmic rays are produced in supernova remnants, so the density of ionising cosmic rays will depend on the SFR, averaged over a ∼ 20 Myr timescale; see Papadopoulos 2010; Papadopoulos et al. 2011). In DSFGs, therefore, where the SFR can be many hundreds or even thousands of solar masses per year, the effect of cosmic rays must be included in the modelling. This requires the use of modern PDR codes which take this effect into account.

Carbon in the ISM of lensed dusty galaxies at z ∼ 4 13 14 Bothwell et al.

SPT0300-46 SPT0441-47 SPT0418-47 SPT0345-47 SPT0113-46 SPT0459-59 SPT0529-54 SPT0532-50 SPT2146-55 SPT2132-58 SPT2103-60 SPT2147-50 log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] log G [Habing]

Model description

We model our line intensities and line ratios using the 3dpdr code (Bisbas et al. 2012) which has been fully benchmarked against the tests discussed in Röllig et al. (2007).

We use 3d-pdr to generate a grid of different uniformdensity, one-dimensional simulations in which we vary the density (specifically the H-nucleus number density, nH) and the UV radiation field (G). The parameter space covers 2 log(nH/cm -3 ) 7 in density and -0.2 log(G/G•) 5.8 in UV field, where G• corresponds to the radiation field strength in Habing units. In general we adopt the heating and cooling functions as described in Bisbas et al. (2012) with the following updates. We use the Bakes & Tielens (1994) PAH photoelectric heating with the modifications suggested by Wolfire et al. (2003) to account for the revised PAH abundance estimate from Spitzer data. We also include the PAH scaling factor given by Wolfire et al. (2008).

We calculate the formation rate of H2 on grains using the treatment of Cazaux & Tielens (2002a,b) and Cazaux & Tielens (2004). We use a subset of the UMIST 2012 network (McElroy et al. 2013) of 33 species (including e -and 330 reactions) and we adopt solar undepleted elemental abundances (Mg= 3.981 × 10 -5 , C= 2.692 × 10 -4 , He= 0.85, O= 4.898 × 10 -4 , assuming that H = 1; see Asplund et al. 2009). The cosmic-ray ionisation rate is taken to be ζCR = 10 -15 s -1 which is approximately 100 times higher than the average ζCR of Milky Way (see Papadopoulos et al. 2010, who demonstrated that ζCR scales with SFR; we discuss the effect of varying this value below). When calculating atomic and molecular emissivities, we integrate clouds to a depth of AV = 7 (see Pelupessy et al. 2006). AV may be converted into an equivalent column density NH via the constant AV0 = AV/NH = 6.289 × 10 -22 mag cm 2 . Our model clouds have uniform density profiles, and a microturbulent velocity of v turb = 1.5 km s -1 to account for microturbulent heating.

PDR modelling results

Results for line ratios for each of our 13 galaxies are shown in Fig. 8. We have considered the ratios

[Cii]/[Ci](1 -0), CO(4 -3)/CO(2 -1), and [Ci](1 -0)/CO(4 -3). The ratio [Cii]/[Ci](1-0)
is an effective tracer of the UV field strength (being essentially insensitive to gas density), while the ratios CO(4-3)/CO(2-1) and [Ci](1-0)/CO(4-3) are highly sensitive to gas density, while being comparatively unaffected by the UV field. Using these observed line ratios, we find probability density distributions for the mean density and the mean UV field strength for each galaxy, based on our grid of one-dimensional PDR models. Using these ratios together provides a simultaneous constraint on both the gas density and the UV field strength in the emitting gas (with the related uncertainties being dictated by the uncertainties on the line ratios). Values for the best-fitting density and UV field strength are given in Table 3. From our sample, we are unable to constrain the ISM conditions in the five DSFGs SPT0125-50, SPT0459-59, and SPT0532-50 (which lack the [Cii] line observations required to measure the UV field strength), SPT0345-47 (which lacks a [Ci] detection), and SPT0529-54 (which lacks a low-J CO detection, and is poorly fit by our models). Fig. 9 shows the distribution of best fitting UV field strength and ISM density values (for the remaining 8 DSFGs for which we could constrain these parameters). For comparison, we have included the approximate ranges of log n and log G values exhibited by 'normal' (i.e., 'Main Sequence') galaxies (Malhotra et al. 2001), local starbursts (Stacey et al. Carbon in the ISM of lensed dusty galaxies at z ∼ 4 15

ID log(n) [cm -3 ] log(G 0 ) Gas Temperature [K] SPT0113-46 4.6 ± 0.4 3.3 ± 0.6 22 +18 -8 SPT0125-50 - - - SPT0300-46 5.1 ± 0.5 1.3 ± 1.2 14 +1 -1 SPT0345-47 - - - SPT0418-47 5.4 ± 0.3 4.5 ± 0.7 42 +20 -24 SPT0441-46 4.1 ± 0.2 3.4 ± 0.8 26 +18 -9 SPT0459-59 - - - SPT0529-54 - - - SPT0532-50 - - - SPT2103-60 4.9 ± 0.4 4.0 ± 0.7 38 +15 -24 SPT2132-58 6.0 ± 0.4 3.7 ± 0.8 21 +8 -5 SPT2146-55 5.8 ± 0.2 2.6 ± 0.2 15 +1 -1 SPT2147-50 5.4 ± 0.2 3.8 ± 0.8 19 +26 -5
Table 3. The inferred values of the density (n), radiation field (G 0 ), and gas temperature for the DSFGs analysed in this work. See Figs. 7 and 9 for the relevant plots. As described in the text, these values were derived using the code 3d-pdr, by integrating line emissivities up to A V = 7, and by setting the cosmic ray flux to be 100 times that of the Milky Way. Gas temperatures are defined at A V = 0.1. Some uncertainties (SPT2146-55, SPT0300-46) are unrealistically small; this is due to the flatness of the temperature map at the position of the n,G values of these galaxies. In reality the error bars on these values are certainly larger.

1991), and local ULIRGs (Davies et al. 2003). Unsurprisingly, the DSFGs in our sample display higher densities and UV field strengths than the sample of 'normal' local galaxies. They also represent a wide range in both log n and log G, and are not obviously comparable to either the class of local ULIRGs, or the starburst galaxies. Our DSFGs have a mean derived density of < log n >= 5.2 ± 0.6 cm -3 , and < log G >= 2.9 ± 1.5. This mean density is greater than the upper range of densities derived for local ULIRGs by Davies et al. (2003). As outlined above, we integrate our clouds up to a depth of AV = 7. Varying this limit typically affects the strength of the derived UV radiation field (lower AV limits resulting in lower values of log G, while derived densities are unaffected.

Gas temperatures

As 3d-pdr predicts a unique gas temperature for each value of log n and log G, it is also possible to use the best-fitting values of log n and log G as a temperature diagnostic. The temperature in the model varies as a function of AV, with high temperatures on the outside of clouds where the gas interacts directly with the incoming UV field, to lower temperatures in the shielded interiors. Here, we analyse the temperature of the gas that is the source of [Ci] and CO emission -that is, the AV at which the [Ci] and CO emissivity peaks. This is generally at Av ∼ 3. Figure 10 shows the distribution of gas temperatures, defined at Av ∼ 3, as a function of log n and log G. The n, G pairs for each DSFG are overplotted on this gas temperature map. We have also listed gas temperatures, with uncertainties, in Table 3. We find a mean gas temperature for our sources of 25 K. In Fig. 11 (2016). These temperatures are slightly lower than typical dust temperatures in SPT-DSFGs (which are 30K -50K), though the relatively large uncertainties on the PDR-derived gas temperatures, coupled with the AV-dependence of the temperatures as a whole, mean that the gas and dust temperatures are consistent within the uncertainties.

Readers will note that our specific results depend on an assumed ionising cosmic ray flux. As cosmic ray flux scales with SFR, it is clear that ζCR will be higher in our DS-FGs than for typical low-z galaxies; as discussed above, we have taken to be ζCR = 10 -15 s -1 , approximately 100 times that of the Milky Way. The choice of ζCR = 100 × ζCR,MW is a conservative lower limit, however (given the SFRs of our DSFGs are at least several hundred times the SFR of the Milky Way). It is therefore worth investigating how our derived results would change if we assumed a higher value of ζCR. Re-calculating our values of n, G with CR fluxes ζCR = 5 × 10 -15 s -1 (500 × ζCR,MW), we find that the derived values of n, G typically increase by ∼ 0.5 dex (with temperatures increasing by a factor of ∼ 2 -5). We also note, however, that at these high CR fluxes the model begins to fail to reproduce the line ratios for some DSFGs, with no possible solutions appearing in our n, G parameter space. Despite this, the qualitative direction of increasing the CR flux is clear: further increasing CR flux results in higher values of n, G, and hotter resulting gas temperatures. This is mainly driven by the [Cii]/[Ci](1 -0) ratio (increasing ζCR means that a given [Cii]/[Ci](1 -0) ratio implies higher G). Our primary PDR modelling result (that our DSFGs are characterised by dense gas and strong interstellar UV radiation fields) remains robust under a range of assumed cosmic ray fluxes.

Comparing to other PDR models

Finally, we turn to a comparison of our results with those produced by alternate PDR models. As discussed above, previous studies of PDR regions in DSFG-like objects (Daniel-son et al. 2011;Alaghband-Zadeh et al. 2013) have used PDR models presented by Kaufman et al. (1999). It is therefore illuminating to compare the results produced by 3d-pdr with those produced by the older Kaufman et al. (1999) models. Our method for deriving densities and UV field strengths from the Kaufman et al. (1999) models is identical to the method we use for 3d-pdr; we produce contour plots based on the line ratios used above ([Cii]/ [Ci](1 -0), CO(4 -3)/CO(2 -1), and [Ci](1 -0)/CO(4 -3)), and take log n and log G to be the peak of the probability distribution function defined by these line ratio constraints. Fig. 12 shows the comparison between physical parameters derived using these two models. It can be seen that the UV field strengths derived by the two models are broadly comparable, with no systematic differences (albeit with relatively large scatter). The densities, however, show a marked difference; the densities derived by 3d-pdr are systematically higher than those derived by Kaufman et al. (1999). The mean derived density produced by the Kaufman et al. (1999) models is < log n >= 4.4 ± 0.4 cm 3 , a factor of ∼ 6 lower than that produced by 3d-pdr (< log n >= 5.2 ± 0.6 cm 3 ). This is also the reason for the higher derived densities in this work, compared to densities derived for SPT-DSFGs by Gullberg et al. (2015) (who quote densities ranging from 2 < log(n)/cm -3 < 5). Gullberg et al. (2015) measured densities using the older Kaufman et al. (1999) models. Because all previous works applying PDR modelling techniques to DSFGs have tended to use Kaufman et al. (1999) models, it seems that the gas density in these extreme high-z objects may have been significantly underestimated.

DISCUSSION

Gas density

Throughout this work, results from both simple line ratio analysis ( §3.3) and more complex PDR modelling ( §4.2) have shown that our sample of strongly-lensed DSFGs exhibits very high ISM densities -higher than both local ULIRGs and unlensed DSFGs at z ∼ 2.

As shown in §3.3 (Fig. 3 (2013) sample to be dominated by non-detections. This is not the case.

The other potential selection effect that could bias the SPT sample is the effect of differential gravitational lensing, as discussed in §3.2 above. The SPT-DSFGs are offset from the Alaghband-Zadeh et al. ( 2013) DSFGs in both L [CI](1-0) /LFIR and L [CI](1-0) /L CO(4-3) . Firstly, while the ratio L [CI](1-0) /L CO(4-3) is potentially affected by differential lensing (due to observational size differences between cold and warm gas tracers; Ivison et al. 2011), the offset is also seen in L [CI](1-0) /LFIR. As discussed in §3.2 the ratio L [CI](1-0) /LFIR is unlikely to be biased due to differential lensing. Secondly, as the specific source-lens geometry varies from source to source, it is unlikely that a combination of lensing geometry and source composition could conspire to produce an offset for our entire sample (from Fig. 3, the 7 densest and strongest UV field galaxies are all SPT sources).

Of course it is possible, given the fairly low significance of the differences between Alaghband-Zadeh et al. ( 2013) DSFGs and SPT-DSFGs (p = 0.061) that the apparent difference in the distributions is due to low number statistics (just 10 and 12 points in each dataset, respectively). However, our PDR model results point to the same conclusion. If SPT-DSFGs are denser than unlensed DSFGs at z ∼ 2, this observation requires some explaining.

A higher average density for SPT-DSFGs would suggest that SPT-DSFGs have some combination of higher average ISM masses, or smaller radii, compared to unlensed DSFGs at z ∼ 2. As discussed earlier in §3.5 (and shown in Fig. 6), there is no difference in gas masses between SPT-DSFGs and the comparison sample. Furthermore, Spilker et al. (2016) found that the observed angular sizes of SPT-DSFGs (at <z>∼ 4) are consistent with the angular sizes of lower redshift unlensed DSFGs (at <z>∼ 2). Allowing for the redshift evolution of the angular scale, this implies that DSFGs in the higher-redshift SPT sample are more physically compact (1 at z = 3.5 corresponds to 7.47 kpc, while 1 at z = 2.2 corresponds to 8.42 kpc). Given the lack of difference in gas mass between our sample and the lower-z unlensed DSFGs, the smaller physical sizes of SPT-DSFGs would result in gas densities roughly (8.42/7.47) 3 ∼ 40% higher than unlensed DSFGs at lower redshift. Both our PDR model results and our line ratio analysis point to the ISMs in our sample of lensed DSFGs being the densest, most extreme star forming environments in the early Universe. et al. (2012) find that the ratio M(H2) dense /M(H2) total displays a bi-modality in the local Universe, with ULIRG/merger systems having elevated dense/total gas ratios relative to secular/disk star forming galaxies. As a result, they claim that this ratio can be used to characterise the star formation 'mode' in galaxies. Any number of dense/total gas tracers may be used for this purpose, but here (following Papadopoulos et al. 2012) we use the ratio r CO43/CI = L CO(4-3) /L CI(1-0) .

The star formation mode

Papadopoulos

[Ci](1 -0) traces total gas, with a critical density of nc ∼ 5 × 102 cm -3 , while CO(4 -3) traces dense gas with a critical density of nc > 5 × 10 5 cm -3 . Papadopoulos et al.

(2012) find that galaxies forming stars in a 'ULIRG/merger' Carbon in the ISM of lensed dusty galaxies at z ∼ 4 17 As in the text, we have taken A V = 7, and a cosmic ray flux rate 100 times that of the Milky Way. The x-axis shows values derived using the PDR models presented by Kaufman et al. (1999), which have been used by previous authors to investigate the ISM of DSFGs (e.g., Alaghband-Zadeh et al. 2013). Points have been colour-coded according to their CO(4 -3)/[Ci](1 -0) ratio, which (as discussed in §3.5) traces the dense/total gas ratio and provides an indication of the 'star formation mode', from disk-like to merger-like. We find that while the 3d-pdr code produces UV field strengths comparable to those estimated by the older K99 models, the gas densities produced by 3d-pdr are approximately 1 dex higher than those produced by the K99 models.

mode have < r CO43/CI > = 4.55 ± 1.5, while secular disk galaxies show typical values of < r CO43/CI > = 0.89 ± 0.44.

Our DSFGs have a continuous distribution in r CO43/CI , with no suggestion of bi-modality. Our DSFGs span a range 1.4 < r CO43/CI < 4.0, with a mean value of <r CO43/CI > = 2.6 ± 0.93 . Our DSFGs do not lie exclusively in either the 'ULIRG/merger' regime or the 'secular disk' regime. While our galaxies are undoubtably not forming stars in the same 'mode' as local disk galaxies (which generally exist in a quiescent equilibrium between inflows, star formation, and outflows), it seems that neither are they directly comparable to local ULIRGs/mergers, which have a shock-compressed ISM forming stars in a central compact burst. A similar conclusion was reached by Bothwell et al. (2013), who found that the dynamical properties of their sample of ∼40 DS-FGs could not be explained using a single 'disk' or 'merger' model applied to the entire population.

Dust based gas masses, and the effect of varying [Ci] abundance

Above, we compared gas masses derived using the [Ci] emission line with those derived using a more traditional method -a low-J emission line of CO. We found that, given an assumed [Ci] abundance of XCI = 3 × 10 -5 , [Ci]-based gas masses were several times higher than those deriving using CO and a CO-to-H2 conversion αCO = 0.8, implying that a higher CO-to-H2 conversion factor (αCO = 2 -2.4) may be appropriate. This is far higher than generally assumed for U/LIRGs and DSFGs, though we note that the original estimates of αCO in ULIRGs (Solomon et al. 1997) were fairly uncertain, displaying a large dispersion.

As noted in the introduction, gas masses can also be estimated using dust masses (combined with an assumption about the gas-to-dust ratio). Gas masses for a sample of SPT-DSFGs have been calculated using this method by Aravena et al. (2016). It is therefore possible to compare gas masses derived using these three independent methods: [Ci] CO, and dust.

By assuming a gas-to-dust ratio of δGDR = 100, Aravena et al. (2016) combine their observed values of L CO with derived dust masses to estimate a CO-to-H2 conversion factor of αCO ∼ 1 for that sample of SPT-DSFGs. This is in tension with our [Ci] results above, which found an average αCO ∼ 2.0 -2.4 for an almost identical sample. The discrepancy with our [Ci]-derived values must therefore be attributable to uncertainty in either the gas-to-dust ratio or the [Ci] abundance; neither of these quantities can be measured in our high-z DSFGs, and so both have to be assumed (taking cues from studies of local galaxies).

Both the gas-to-dust ratio and the [Ci] abundance are functions of metallicity, though the gas-to-dust ratio scales more steeply with metallicity (roughly quadratically; Rémy-Ruyer et al., 2014). The existence of low metallicities in our DSFGs could therefore raise the gas-to-dust ratio, thus raising the dust-based value of αCO into closer agreement with the [Ci]-based αCO presented above. The other alternative is that our assumed value of [Ci] abundance (XCI = 3 × 10 -5 ) could be too low, leading us to overestimate gas masses derived using [Ci].

We now consider this problem from both sides: firstly, what gas-to-dust ratio would be required to bring the dustbased gas masses into agreement with our [Ci] gas masses? While Aravena et al. (2016) calculated gas masses based on a δGDR = 100, raising our dust-based gas masses into agreement with our [Ci] gas masses would require δGDR = 200-240. This is not physically implausible; Rémy-Ruyer et al., (2014) report that the gas-to-dust ratio follows a power law with metallicity, log δGDR = 2.21+2.02(x -x), where x = 12+log(O/H), and x is solar metallicity (8.69). A gas-todust ratio δGDR = 200 would require a metallicity ∼ 0.1 dex below solar -not unfeasible at the redshifts of our DSFGs.

Secondly, we consider the value of the [Ci] abundance required to reduce our [Ci]-based gas masses into line with those derived with dust masses. To reduce our [Ci]-based gas masses into agreement with those derived using a gasto-dust ratio δGDR = 100 (and αCO = 1) would require increasing our assumed [Ci] abundance to XCI = 7 × 10 -5 . This is again not physically implausible; the dense nuclei of nearby starbursts show elevated carbon abundances (White et al. 1994 report a [Ci] abundance of XCI = 5 × 10 -5 in the nucleus of the starburst M82). Furthermore, Bisbas et al. (2015) describe a physical model in which CO molecules are dissociated by ionising cosmic rays -CO dissociation can increase the abundance of [Ci], suggesting that strongly star-forming galaxies (like our DSFGs) with high cosmic ray fluxes are likely to show elevated [Ci] abundances. Given this, a [Ci] abundance of XCI = 7 × 10 -5 is certainly possible.

SUMMARY AND CONCLUSIONS

There has been much discussion in recent years as to the relation between the population of mm/submm-selected DS-FGs and more familiar classes of local galaxy (such as local ULIRGs). This work adds to a growing body of evidence that DSFGs in the early Universe represent a heterogeneous population with properties that cannot be understood as being simply 'scaled up' versions of IR-luminous starbursts in the z ∼ 0 Universe. Assuming a 'standard' carbon abundance, XCI ∼ 3 × 10 -5 suggests extreme gas fractions (fgas ∼ 0.6 -0.8) in tension with those calculated using both low-J CO and dust masses. Reconciling these methods may require high carbon abundances in our DS-FGs (XCI ∼ 7 × 10 -5 ), which may result from the cosmic ray dissociation of CO into [Ci]. The objects in this work have ISMs characterised by dense, carbon rich gas unlike that in ULIRGs in the local Universe.

In this work, we have presented an analysis of a sample of 13 strongly-lensed DSFGs at z ∼ 4, selected at 1.4mm by the South Pole Telescope survey. We have used ALMA Band 3 observations of the emission line of atomic carbon to characterise the properties of, and conditions within, the ISM of these extreme galaxies. Our main conclusions are as follows:

• Using the luminosity of [Ci] as a tracer of the total gas mass, and assuming a [Ci]/H2 abundance ratio of 3 × 10 -5 , we find a mean H2 mass for our DSFGs of (6.6 ± 2.1) × 10 10 M , and a typical (sample-averaged) gas fraction of fgas ∼ 0.6.

• This gas mass is higher than the value derived from observations of CO emission lines (assuming a 'standard' CO-to-H2 conversion factor of αCO = 0.8). It is also higher than the gas mass derived using observations of the dust continuum (assuming a gas-to-dust ratio δGDR = 100). We find that a CO-to-H2 conversion factor of αCO = 2 -2.4, and a gas-to-dust ratio δGDR ∼ 200 would be needed to bring the two gas mass estimates into agreement. These values are higher than that generally adopted for extreme DSFGs. Alternatively, the [Ci] abundance may be very high: a [Ci]/H2 ratio of 7 × 10 -5 would lower the [Ci]-based gas masses into agreement with the conventional CO and dust measurements. In the latter case, the high [Ci] abundance could be driven by the dissociation of CO into [Ci] by ionising cosmic rays.

• We use a range of ancillary line observations for our galaxies, and the PDR modelling code 3d-pdr, to estimate the conditions within the ISM of our galaxies. Using a number of line ratios, we find that our DSFGs exhibit strong UV field strengths and dense gas emitting regions, comparable to (or even denser than) local ULIRGs and unlicensed DSFGs at z ∼ 2. Furthermore, we find gas densities significantly denser than those derived using older PDR models commonly used to model ISM conditions in high-z DSFGS.

• We also use this PDR code to estimate the gas temperature within the ISM of our DSFGs. We find evidence for gas with typical temperatures at AV = 3 of ∼ 25K. Within the uncertainties, these temperatures are consistent with derived dust temperatures for our DSFGs.

Figure 1 .Figure 1 .

 11 Figure1. Ci spectra for galaxies in this work (grey), binned to 50 km s -1 resolution. Where ATCA CO(2 -1) spectra are available, they have been over plotted (at matched velocity resolution) in blue. CO(2 -1) spectra have their fluxes arbitrarily normalised, and are intended for comparison of line profiles only.

Figure 3 .

 3 Figure 3. Plot showing the ratio L [CI](1-0) /L FIR plotted against the ratio L [CI](1-0) /L CO(4-3) , for the SPT-DSFGs in this paper, and the DSFGs from Alaghband-Zadeh et al. (2013). For comparison, we have also included a sample of local galaxies observed with Herschel/SPIRE, presented byKamenetzky et al. (2016). The ratio L[CI](1-0) /L FIR is an approximate (but non-linear) tracer of the UV field strength, while the ratio L[CI](1-0) /L CO(4-3) is an approximate (but again non-linear) tracer of the gas density.

Figure 4 .

 4 Figure 4. Example figure showing the dependency of the derived mass of atomic carbon on the assumed excitation temperature Tex.The track shown is for the mean intrinsic (i.e., corrected for lensing magnification) carbon luminosity of our sample, and the dashed lines show the range implied by the uncertainty on this value. It can be seen that for a wide range of excitation temperatures (> 20K), the derived mass of atomic carbon is only very weakly dependent on the specific temperature assumed.

(Figure 5 .

 5 Figure 5. The linewidths of [Ci] plotted against the linewidth of CO(2 -1) for the DSFGs in this work. Also shown are the linewidths for the Alaghband-Zadeh et al. (2013) sample (for the AZ13 sample, the FWHM of the CO(3 -2) line is shown). The solid and dashed lines respectively show a ratio of unity and a factor of 2 variation. There is a close kinematic correspondence between the [Ci] and low-J CO lines -in the majority of cases being equivalent within the uncertainties. There is also no systematic difference between the two (< FWHM CI /FWHM CO >= 1.03 ± 0.40).

  [Ci] lines to those of the CO(2-1) lines (only possible of course for the DSFGs in our sample that have both lines detected). We also include four DSFGs from Alaghband-Zadeh et al. (2013) which have CO(3 -2) line measurements. (All other DSFGs in the Alaghband-Zadeh et al. (2013) sample only have CO emission at Jup 4 observed, which are increasingly poor tracers of the cool component of the ISM.)

Figure 6 .

 6 Figure 6. Left Panel: Comparison of molecular hydrogen masses, derived via [Ci] and CO emission, for the SPT DSFGs in this paper and the DSFGs in Alaghband-Zadeh et al. (2013) with both [Ci] and CO observations. To derive CO-based H 2 masses, we have assumed α CO = 0.8. The dashed diagonal line indicates equal masses, while the dotted lines show a factor of two variation each way. DSFGs for which a higher-J CO line has been used to infer a gas mass have been plotted using darker colours. Right Panel: Ratio between [Ci]-derived molecular gas mass and CO(1 -0) luminosity, for the same galaxies. This ratio can be interpreted as giving the CO-to-H 2 conversion factor α CO . All (but one) SMGs in both samples have implied values of α CO far above the canonical 'ULIRG' value, and several SPT-DSFGs have an α CO comparable to that of the Milky Way.

Figure 7 .

 7 Figure 7. The observed linewidth (full width half maximum, FWHM), plotted against total H 2 mass for SPT-DSFGs and literature DSFGs. For SPT-DSFGs (red), and DSFGs taken from Alaghband-Zadeh et al. (2013) (blue), linewidths and H 2 masses are derived using [Ci]. For the literature DSFGs, taken from Bothwell et al. (2013), linewidths and H 2 masses are derived using either CO(2 -1) or CO(3 -2). The purple line represents a simple dynamical mass model (given in Eq. 8).

Figure 8 .

 8 Figure8. Probability distribution functions of log(n) and log(G), with tracks showing the constraints due to measured line ratios, for each of the SPT DSFGs in this work. These results were produced using the code 3d-pdr, by integrating line emissivities up to A V = 7, and by setting the cosmic ray flux to be 100 times that of the Milky Way. Contours of the resultant 'consensus value' of log(n) and log(G) are overlaid in grey. DSFGs with only a single available line ratio (such as SPT0459-59) have unconstrained densities and UV field strengths.

Figure 9 .

 9 Figure 9. Plot showing the 'best fitting' values of density (n) and UV-field strength (G 0 ) for our DSFGs. Results were calculated by integrating to A V = 7, and with a cosmic ray flux rate 100 times that of the Milky Way. For reference, the approximate values for Main Sequence galaxies (Malhotra et al. 2001), local starbursts (Stacey et al. 1991), and local ULIRGs (Davies et al. 2003) are shown. The average density for the four DSFGs without a constraint on G 0 is shown at the bottom.

Figure 10 .

 10 Figure 10. Plot showing the gas temperature implied by the density (n) and UV-field strength (G 0 ) outputted by 3d-pdr. The location of our sources is shown. Our sample has a mean gas temperature of 25K. Temperatures are calculated at a depth of A V = 3, corresponding to the peak emissivity of [Ci] and CO.

Figure 11 .

 11 Figure11. Gas temperature implied by the density (n) and UVfield strength (G 0 ) outputted by 3d-pdr (defined at A V = 3), plotted against dust temperature derived using SED fitting to far-IR photometry.

  ), SPT-DSFGs are offset from unlensed DSFGs at z ∼ 2 towards lower L [CI](1-0) /LFIR and lower L [CI](1-0) /L CO(4-3) ; this offset physically implies denser ISMs with stronger UV radiation fields. It is unlikely that the SPT sample is skewed towards the highest density, highest UV field strength objects because of any simple selection effects. Firstly, we note that the comparison sample of unlensed z ∼ 2 DSFGs presented in Alaghband-Zadeh et al. 2013 were not selected for [Ci] observation based on their CO brightness 2 . Additionally, while the [Ci] observations for the SPT sample come from ALMA (and are thus more sensitive than the IRAM-PdBI [Ci] observations presented by Alaghband-Zadeh et al. 2013), the Alaghband-Zadeh et al. (2013) sample contains only two non-detections. If the SPT sample was offset from the Alaghband-Zadeh et al. (2013) sample solely because of ALMA's increased [Ci] sensitivity, we should expect the Alaghband-Zadeh et al.

Figure 12 .

 12 Figure 12. A comparison of log n (left panel) and log G (right panel), as derived by two different PDR model codes. The y-axes in each panel show values derived by the PDR code used in this work, 3d-pdr.As in the text, we have taken A V = 7, and a cosmic ray flux rate 100 times that of the Milky Way. The x-axis shows values derived using the PDR models presented byKaufman et al. (1999), which have been used by previous authors to investigate the ISM of DSFGs (e.g.,Alaghband-Zadeh et al. 2013). Points have been colour-coded according to their CO(4 -3)/[Ci](1 -0) ratio, which (as discussed in §3.5) traces the dense/total gas ratio and provides an indication of the 'star formation mode', from disk-like to merger-like. We find that while the 3d-pdr code produces UV field strengths comparable to those estimated by the older K99 models, the gas densities produced by 3d-pdr are approximately 1 dex higher than those produced by the K99 models.

Table 1 .

 1 Line luminosities, far-IR luminosity, [Ci]/FIR ratios, and lensing magnifications (µ) for the DSFGs studied in this work. All luminosities are given in solar units, and have not been corrected for gravitational lensing. Where 12 CO(4 -3) transition line luminosities are not directly measured, they are inferred from the 12 CO(5 -4) line luminosity using the conversions derived byBothwell et al. (2013), and are shown in parentheses. The far-infrared luminosities are derived by integrating under a modified blackbody curve from 8-1000 µm. Lensing magnifications are derived from visibility-based lens models fit to ALMA 870µm observations(Spilker et al. 2016).[Cii] data are taken fromGullberg et al. (2015). CO(2 -1) data are taken fromAravena et al. (2016) -other CO lines are taken from the program described inWeiss et al. (2013).

	10 8 L ]	[10 8 L ]	[10 8 L ]	[10 8 L ]	[10 10 L ] [10 13 L ]	×10 -6

Table 2 .

 2 Masses of atomic carbon and molecular hydrogen (derived using the [Ci] flux) for our sample, and the implied CO-to-H 2 conversion factor.

	46	0.66 ± 0.13	3.81 ± 0.92	2.4 ± 0.6
	SPT0125-50	0.71 ± 0.15	4.08 ± 1.09	1.1 ± 0.3
	SPT0300-46	1.12 ± 0.50	6.48 ± 3.45	1.4 ± 0.7
	SPT0345-47	< 0.62	< 3.56	< 0.9
	SPT0418-47	0.35 ± 0.08	2.03 ± 0.60	2.3 ± 0.7
	SPT0441-46	0.73 ± 0.29	4.24 ± 2.06	2.6 ± 1.3
	SPT0459-59	3.82 ± 1.10	21.9 ± 7.60	2.7 ± 0.9
	SPT0529-54	0.70 ± 0.13	4.05 ± 0.90	1.6 ± 0.3
	SPT0532-50	1.05 ± 0.24	6.05 ± 1.71	1.1 ± 0.3
	SPT2103-60	0.55 ± 0.13	3.21 ± 0.95	2.1 ± 0.7
	SPT2132-58	0.79 ± 0.28	4.54 ± 1.97	1.2 ± 0.5
	SPT2146-55	2.17 ± 0.56	12.5 ± 3.90	3.8 ± 1.4
	SPT2147-50	1.19 ± 0.35	6.89 ± 2.46	1.5 ± 0.6

  , we compare our derived gas temperatures with dust temperatures, taken fromWeiß et al. 2013, Gullberg et al. 2015, and Spilker et al.

The criteria for observation was a 1.4 GHz continuum detection (providing an accurate position), K-band brightness (allowing a Hα line redshift to be obtained), as well as a southerly declination to allow potential ALMA followup.

This value is in good agreement with the mean r CO43/CI found by Alaghband-Zadeh et al. (2013), who found a sample mean value of 2.5 ± 0.2.

ACKNOWLEDGEMENTS

This paper makes use of the following ALMA data: ADS/JAO.ALMA #2011.0.00957.S, #2011.0.00958.S, #2012.1.00844.S, and #2012.1.00994.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This work has made use of the NASA ADS. The SPT is supported by the National Science Foundation through grant PLR-1248097, with partial support through PHY-1125897, the Kavli Foundation and the Gordon and Betty Moore Foundation grant GBMF 947. MSB is supported by STFC grants ST/M001172/1 and ST/K003119/1. We acknowledge support from the U.S. National Science Foundation under grant No. AST-1312950. M.A. acknowledges partial support from FONDECYT through grant 1140099.

APPENDIX A: LINE FLUXES