Neutral carbon and highly excited CO in a massive star-forming main sequence galaxy at z=2.2

TL;DR

This study uses ALMA observations of multiple CO and neutral carbon lines in a z=2.2 star-forming galaxy to analyze its gas properties, star formation activity, and excitation conditions, revealing complex gas components and challenging the use of CO(7-6) as a star formation tracer.

Contribution

First detailed multi-line CO and [CI] study of a main sequence galaxy at z=2.2, providing new insights into its gas mass, excitation, and star formation conditions.

Findings

Gas mass estimated at ~2x10^11 M_sun

Presence of a second, highly excited gas component

CO(7-6) not a reliable star formation tracer

Abstract

We used the Plateau De Bure Interferometer to observe multiple CO and neutral carbon transitions in a z=2.2 main sequence disk galaxy, BX610. Our observation of CO(7-6), CO(4-3), and both far-infrared(FIR) [CI] lines complements previous observations of H$\alpha$ and low-J CO, and reveals a galaxy that is vigorously forming stars with UV fields (Log($G$ G$_0^{-1}) \lesssim3.25);$ although less vigorously than local ultra-luminous infrared galaxies or most starbursting submillimeter galaxies in the early universe. Our observations allow new independent estimates of the cold gas mass which indicate $M_\textrm{gas}\sim2\times10^{11}$M$_\odot$, and suggest a modestly larger $\alpha_{\textrm{CO}}$ value of $\sim$8.2. The corresponding gas depletion timescale is $\sim$1.5 Gyr. In addition to gas of modest density (Log($n$ cm$^3)\lesssim3$ ) heated by star formation, BX610 shows evidence for a significant second gas component responsible for the strong high-J CO emission. This second component might either be a high-density molecular gas component heated by star formation in a typical photodissociation region, or could be molecular gas excited by low-velocity C shocks. The CO(7-6)-to-FIR luminosity ratio we observe is significantly higher than typical star-forming galaxies and suggests that CO(7-6) is not a reliable star-formation tracer in this galaxy.

Tables

Table 1. Table 1: BX610 observed quantities

Parameter	Unit	Value	Reference
RA	hh:mm:ss.ss	23:46:09.43	(Erb et al. 2006b)
Dec.	dd:mm:ss.ss	12:49:19.21	(Erb et al. 2006b)
redshift	z	2.2103	(Förster Schreiber et al. 2009)
CO(1-0)	Jy km s-1	0.105 $\pm$ 0.015	(Bolatto et al. 2015)
FWHM	km s-1	294 $\pm$ 49	(Bolatto et al. 2015)
CO(3-2)∗	Jy km s-1	0.87 $\pm$ 0.11	(Bolatto et al. 2015)
FWHM	km s-1	266 $\pm$ 35	(Bolatto et al. 2015)
CO(4-3)	Jy km s-1	1.82 $\pm$ 0.26	This work
FWHM	km s-1	301 $\pm$ 33	This work
CO(7-6)	Jy km s-1	1.88 $\pm$ 0.28	This work
FWHM	km s-1	258 $\pm$ 44	This work
$\left[\text{C}\text{i}\right]$(1-0)	Jy km s-1	0.80 $\pm$ 0.15	This work
FWHM	km s-1	349 $\pm$ 52	This work
$\left[\text{C}\text{i}\right]$(2-1)	Jy km s-1	1.34 $\pm$ 0.24	This work
FWHM	km s-1	247 $\pm$ 54	This work
$S_{150\text{GHz}}$	mJy	0.41 $\pm$ 0.09	This work
$S_{250\text{GHz}}$	mJy	2.39 $\pm$ 0.31	This work

Table 2. Table 2: BX610 derived parameters

Parameter	Unit	Value	Reference
Log(L${}_{\text{IR}}$ L${}_{}{}^{-1}{}_{\odot }{}^{}$)	-	12.63$\pm$0.10	This work
Log(L${}_{\text{FIR}}$ L${}_{}{}^{-1}{}_{\odot }{}^{}$)	-	12.27$\pm$0.10	This work
Log(SFR (M⊙/yr)-1)	-	2.15$\pm$0.10	This work
$M_{\text{gas(CO)}}$	1011M⊙	1.1$\pm$0.1∗	(Bolatto et al. 2015)
$M_{\text{gas([C}\text{i}\text{])}}$	1011M⊙	2.0$\pm$0.4∗	This work
$M_{\text{gas(dust)}}$	1011M⊙	2.4$\pm$0.5∗	This work
$M_{\ast }$	1011M⊙	1.0$\pm$0.3	(Tacconi et al. 2013)
${\tau }_{\text{depletion}}$	Gyr	1.5$\pm$0.3	This work
$T_{\text{ex([C}\text{i}\text{])}}$	K	31.8$\pm$6.9	This work
LVG model parameters
$T$	K	30$\pm$4	This work
Log($n$ cm3)	-	4.3$\pm$0.1	This work
Log([CO]/[C])	-	0.47$\pm$0.17	This work
Two-component PDR model parameters
Log($G_1$ G${}_{}{}^{-1}{}_0{}^{}$)	-	2.25$\pm$1	This work
Log($n_1$ cm3)	-	5.0$\pm$0.5	This work
Log($G_2$ G${}_{}{}^{-1}{}_0{}^{}$)	-	3.125$\pm$0.375	This work
Log($n_2$ cm3)	-	3.125$\pm$0.625	This work
PDR and shock model parameters
Log($G$ G${}_{}{}^{-1}{}_0{}^{}$)	-	3.25$\pm$0.25	This work
Log($n_{\text{PDR}}$ cm3)	-	3.875$\pm$0.375	This work
$v_{\text{shock}}$	km s-1	10†	This work
Log($n_{\text{shock}}$ cm3)	-	4.3†	This work

Equations

$$T_{\text{ex}}=38.8\times \text{ln}(2.11\times \frac{L_{\text{[CI]}(1-0)}^{\prime }}{L_{\text{[CI]}(2-1)}^{\prime }}{)}^{-1},$$

$$M_{\text{CI}}=5.706\times 10^{-4}Q(T_{\text{ex}})\frac{1}{3}{\text{e}}^{23.6/T_{\text{ex}}}{L}_{\text{[CI]}(1-0)}^{\prime }[{\text{M}}_{\odot }],$$

Full text

1 Introduction

Investigating the interstellar medium (ISM) has long been a crucial part of understanding the processes at work in galaxies and characterizing them. The multitude of easily accessible rotational CO lines have been particularly useful in revealing the density and temperature of molecular gas in galaxies in the local universe.

Recent advances have made it possible to observe these lines at increasingly early epochs, and there are several samples of galaxies at significant redshift with detections of CO (e.g., Carilli & Walter 2013). However, currently only a handful of these high-redshift sources have more than a few CO line detections, and those are mostly limited to extreme, high-luminosity systems often undergoing unsustainable starbursts (e.g., Alaghband-Zadeh et al. 2013; Carilli & Walter 2013; Casey et al. 2014; Bothwell et al. 2017; Yang et al. 2017; Wardlow et al. 2018; Dannerbauer et al. 2018). The existence of a galaxy main sequence has emerged as a useful way of contextualizing galaxy star formation activity. Most star-forming galaxies cluster along a positive trend line in star formation rate (SFR) and stellar mass parameter space with a secondary population of starburst galaxies with SFRs enhanced by a factor of a few to tens of times larger than their main sequence counterparts with similar stellar masses. The galaxy main sequence persists throughout the local universe and out to high redshifts with an evolving relationship between SFR and stellar mass (e.g., Brinchmann et al. 2004; Salim et al. 2007; Elbaz et al. 2007; Daddi et al. 2007). Given the abundance of galaxies on the main sequence and the dominance of main sequence galaxies in cosmic luminosity density and star formation rate density during the epoch of peak star formation, many works have examined the relationship between star formation and molecular gas content in main sequence galaxies (e.g., Schmidt 1959; Kennicutt 1998; Gao & Solomon 2004; Krumholz et al. 2009).

In this paper we present observations of four emission lines, CO(7-6), CO(4-3), [Ci] 3P2-3P1, and [Ci] 3P1-3P0 (the latter two hereafter denoted Ci and Ci, respectively), in the massive star-forming galaxy [ESS2003] Q2343-BX610 at z=2.2103, hereafter referred to as BX610. The lines complement previous detections of CO(1-0), and CO(3-2), allowing us to investigate the ISM conditions in this system (Tacconi et al. 2010; Aravena et al. 2014; Bolatto et al. 2015). Unlike most previously detected high-redshift galaxies with an established CO Spectral Line Energy Distribution (SLED), BX610 is thought to be undergoing vigorous star formation (SFR$\sim$60-200 M⊙ yr*-1*) while lying squarely on the z$\sim$2 galaxy main sequence with a stellar mass $M_{*}\sim 10^{11}$M⊙. Few z$\gtrsim$2 main sequence galaxies have been observed in similar gas tracers so far, and even among those few there are likely diverse populations (e.g., Daddi et al. 2015; Popping et al. 2017; Talia et al. 2018). Therefore, understanding the ISM in BX610 is a crucial step in understanding this abundant yet under-investigated galaxy population in the early universe.

In Sect. 2  we present our new observations. In Sect. 3  we present the initial results of our observations and briefly review supporting data from the literature. In Sect. 4  we analyze and discuss these observations and their implications for star formation within the galaxy, and in Sect. 5  we summarize our findings. Throughout the paper we adopt a flat $\Lambda$CDM cosmology based on Planck 2018 results with $\Omega_{\Lambda}$=0.6847, $\Omega_{\textrm{M}}$=0.3153, and $H_{0}$=67.36 km s*-1* Mpc*-1* (Planck Collaboration et al. 2018).

2 Observations

Observations of the CO(4-3), Ci, CO(7-6), and Ci emission lines in BX610 were conducted with the Plateau de Bure Interferometer (PdBI). Observations of the CO(4-3) and Ci lines were performed on July 11, 15, 21, 26, and 27, 2013, in the 5Dq array configuration (maximum baseline of 97m). Observations of the CO(7-6) and Ci lines were performed on July 23, 2014, with six available antennas in the D array configuration (maximum baseline of 144m). All observations used the WideX correlator and were pointed toward the position (RA, DECL.)=(23:46:09.43, 12:49:19.21).

Observations of the CO(4-3) and Ci lines were performed with the 2mm band receivers in two independent setups, tuned to 143.586 and 153.278 GHz, respectively. Observations of the CO(7-6) and [Ci] 2-1 lines were done with the 1mm receivers, in a single tuning at 251.689 GHz, as both emission lines fall within the 3.6 GHz bandwidth provided by the WideX correlator.

The observing times were 12.5h, 14.3h, and 4.8h for the CO(4-3), Ci, and CO(7-6)/Ci observations, respectively, for a total of 31.6 hours. All observations were obtained under good weather conditions and phase stability for these frequencies, with typical seeing values between 0.74-1.50” at 2mm and 0.35” at 1mm. The bright quasars 2251+158 and 3C454.3 were used for bandpass calibration. The sources 2251+158 and J2327+096 were used for gain calibration at 2mm and 1mm, respectively. The standard calibrator MWC349 and Uranus were used for flux calibration. Extragalactic millimeter calibrator source fluxes can vary significantly over time, and in some cases may be estimated based on fluxes from recent epochs and nearby wavelengths. These methods occasionally result in larger systematic errors than expected, and so calibrator fluxes should be checked carefully. In our examination of the calibration we found a systematic discrepancy between the calculated calibrator fluxes at 150 GHz and the fluxes indicated by the ALMA calibrator archive, with the fluxes used in the PdBI calibration enhanced by 20%. To correct for this discrepancy we manually corrected our 150 GHz data values by a factor of 0.8. At the observing frequencies, the primary beam FWHMs correspond to 33” and 20” at 2mm and 1mm, respectively.

All the data were calibrated using the standard pipeline available in the GILDAS software. The visibilities were inverted and cleaned using a tight box around the initially detected source, down to a threshold of 2.5-sigma. The final cubes yield synthesized beam FWHMs of 3.5”x3.3” (P.A. 158 deg), 3.3”x3.0” (P.A. 103 deg), and 2.4”x1.2” (P.A. 197 deg) with rms noise levels of 0.9, 0.6, and 1.1 mJy beam*-1* per channel, with channel widths of 50, 50, and 45 km s*-1* for the CO(4-3), [Ci] 1-0, and CO(7-6)/[Ci] 2-1 observations, respectively.

3 Results

We analyzed our data in MAPPING, part of the GILDAS software suite. Based on UV plane source fitting, we determined all line and continuum emission to be consistent with a point source. The continuum flux density at 250 GHz is determined based on the channels below the CO(7-6) frequency and the channels above the Ci frequency ($\nu\leq$251.009 and $\nu\geq$252.256 GHz). The CO(7-6) and Ci lines cover spectral ranges 251.085$\leq\nu\leq$251.632 GHz and 251.632$<\nu\leq$252.218 GHz, respectively. Fit results for the continuum and continuum-subtracted emission lines are given in Table 1. The line spectra of the peak pixel (with continuum subtracted) are shown in Fig. 1. We also fit the peak line emission with two Gaussians to characterize their velocity distributions. Central frequencies, FWHMs, and amplitudes were allowed to vary independently for both lines. In both cases we find FWHM$\sim$250 km s*-1* and central frequencies consistent with the expected redshift of z=2.2103.

The continuum emission is significantly weaker at 150 GHz and is difficult to characterize in the CO(4-3) band, and therefore we use only the higher sensitivity data from the Ci observation to estimate continuum. The continuum estimate is based on channels outside the observed line emission ($\nu\leq$ 152.971 and $\nu\geq$ 153.636 GHz). The line emission shown in Fig. 1 was fit in the same manner as CO(7-6) and Ci. Fitting Gaussians to the peak pixel spectra of Ci and CO(4-3), we find FWHM of 349 $\pm$ 52 and 301 $\pm$ 33 km s*-1*, respectively.

In addition to the new observations presented here, we also make use of prior CO observations. BX610 was observed and detected in CO(3-2) as part of the Plateau de Bure HIgh-z Blue Sequence Survey (PHIBSS; Tacconi et al. 2010, 2013), and was reanalyzed by Bolatto et al. (2015) as part of a larger sample. The CO(1-0) line was strongly detected by Aravena et al. (2014) and then again detected by Bolatto et al. (2015) using the Karl G. Jansky Very Large Array. These line estimates are presented with our new observations in Table 1. We have included an additional 10% uncertainty on all line and continuum estimates, added in quadrature with the statistical uncertainty, to account for possible flux calibration differences between surveys.

4 Analysis

4.1 Infrared luminosity and SFR estimates

Previous estimates of SFR in BX610 vary by factors of approximately three but consistently confirm its place on the main sequence. Tacconi et al. (2013) estimated the SFR of BX610 to be 212$\pm$74 M⊙ yr*-1*, based on combined H$\alpha$ and UV observations. More recently the SFR has been estimated as 60 M⊙ yr*-1* based on evolutionary synthesis modeling of the UV-near-infrared(NIR) spectral energy distribution (SED), and 115 M⊙ yr*-1* based on newly acquired adaptive-optics-assisted H$\alpha$ observations (Förster Schreiber et al. 2014, 2018). Although these estimates have all accounted for dust extinction through correction factors or implicit dust models, they should be further compared and tested for consistency with an SED-derived estimate which includes observations of the submillimeter continuum arising from the dust itself. We do a joint UV–submillimeter photometry fit with the high-redshift Multi-wavelength Analysis of Galaxy Physical Properties (MAGPHYS) SED analysis tool (da Cunha et al. 2008) using our continuum measurements at 150 and 250 GHz, U, G, R, J, and Ks photometry from Erb et al. (2006b, a), HST NIC2 F160W photometry (Förster Schreiber et al. 2011a), and Spitzer IRAC (4.5 and 8 $\mu$m) and MIPS 24 $\mu$m data extracted from the NASA/ IPAC Infrared Science Archive. The resulting best-fit SED is shown in Fig. 2. The MAGPHYS-derived SFR is Log(SFR (M⊙/yr)-1)=2.15$\pm$0.1 and infrared luminosity (8-1000 $\mu$m) is Log(L${}_{\textrm{IR}}$ L${}_{\odot}^{-1}$)=12.63$\pm$0.10. We also make use of the 30-1000 $\mu$m luminosity, estimated to be Log(L${}_{\textrm{30-1000}}$ L${}_{\odot}^{-1}$)=12.50$\pm$0.10 and the FIR luminosity (42.5-122.5 $\mu$m) Log(L${}_{\textrm{FIR}}$ L${}_{\odot}^{-1}$)=12.27$\pm$0.10. In this case the formal errors derived from MAGPHYS are very small due to the discrete nature of the SED templates included in the fitting library. The resulting likelihood distributions of the relevant parameters have a single value dominating the total integrated likelihood power. Therefore, we have taken the sampling spacing of 0.1 in log space as the uncertainty. These results are summarized in Table 2 along with other derived parameters.

Emission from mid-J CO lines, which usually originates in denser molecular gas, has also been linked to star formation. Several authors have noted the particularly strong correlation between CO(7-6) and infrared luminosity (Lu et al. 2014, 2015; Liu et al. 2015; Yang et al. 2017). From the relationship presented by (Lu et al. 2015; their equation 1) the CO(7-6) in BX610 would indicate a much larger infrared luminosity, L${}_{\textrm{IR}}$=1.2$\times$1013L⊙, and a corresponding SFR$\sim$1200 M⊙ yr*-1* after correcting for a Chabrier initial mass function — strikingly different from ours and previous SFR estimates. Although the MAGPHYS fit, which uses energy balance considerations to consistently interpret UV through millimeter SEDs, should provide the best estimate, without a probe of the peak of the FIR dust emission we cannot completely rule out the possibility of additional obscured star formation. As a simple test we have also investigated the mid-infrared (MIR) through submillimeter SED including Spitzer MIPS and ALMA photometry using the SED template library from Draine & Li (2007). The full template library allows a breadth of infrared luminosities, including values exceeding the CO(7-6) predicted infrared luminosity, and is not well constrained with only our three photometric points. Based on this same SED template library, Béthermin et al. (2015) determined that massive main sequence galaxies at 2$<z<2.5$ are best characterized by SEDs with a mean ionization parameter of ¡U¿=22.6. Fitting our photometry with the same template indicates a modestly larger infrared luminosity, Log(L${}_{\textrm{IR}}$ L${}_{\odot}^{-1}$)=12.89$\pm$0.03, but this still falls short of what would be expected from the CO(7-6) emission. This SED template is also shown in Fig. 2. Alternatively, the extreme discrepancy between the CO(7-6) emission and SFR indications from H$\alpha$ and SED fitting might be caused by non-SFR-related gas heating contributing to the strong CO(7-6) emission. From our MAGPHYS-derived FIR luminosity we find the luminosity ratio, Log($\frac{L_{\textrm{CO}(7-6)}}{L_{\textrm{FIR}}}$)=-4.07$\pm$0.12. This is well above all but a few outliers in the local and high-redshift samples presented in the literature. The relationship from Lu et al. (2015) is based on observations of local galaxies, including both main sequence and starburst galaxies. Although there are fewer observations constraining the relationship at high redshifts, the trend persists and, if anything, shifts to even smaller CO(7-6)-to-L${}_{\textrm{FIR}}$ ratios, making BX610 even more of an outlier (Yang et al. 2017). One of the only sources with a comparable or higher ratio and that has been analyzed with multiple molecular lines is NGC 6240, which is believed to be undergoing significant shock heating (Meijerink et al. 2013; Lu et al. 2014). Due to the extreme discrepancy between SFR as estimated by CO(7-6) and other methods, as well as the unusually high CO-to-infrared luminosity ratio and the possibility of shock excitation contributions, we suggest that CO(7-6) is not a reliable SFR diagnostic in this source.

4.2 Gas mass estimates

Cold gas mass is the most fundamental property governing star formation in galaxies (e.g., Schmidt 1959; Kennicutt 1998; Gao & Solomon 2004; Krumholz et al. 2009; Bothwell et al. 2016). Composed predominantly of H2 with a modest component of helium, cold gas is both the ultimate fuel for star formation and a motivating force by way of gravitationally compelling molecular cloud collapse. Since it is impossible to directly observe H2 from the ground, alternative tracers are used to determine cold gas masses. All of these methods require various assumptions to arrive at the final gas mass value, and these assumptions can come under question in different environments. Our access to several different gas tracers in BX610 provides us with the opportunity to estimate molecular gas masses in multiple ways, and potentially minimize uncertain assumptions.

Molecular CO, which occurs throughout molecular cores, has long been used as an easily observed tracer of molecular hydrogen (Solomon & Vanden Bout 2005). Although the abundance of CO relative to H2 varies significantly across different galactic environments (due largely to variations in metallicity and geometry allowing more or less exposure to dissociating radiation), by assuming a conversion factor, $\alpha_{\textrm{CO}}$, the line luminosity of CO(1-0) is routinely converted into star-forming gas mass (often including helium in addition to H2). This gas mass estimation is directly proportional to $\alpha_{\textrm{CO}}$, which is known to vary widely over different environments, from $\sim$1 in local ultra-luminous infrared galaxies (ULIRGs) to $\gtrsim$50 in low-metallicity systems (e.g., Bolatto et al. 2013). In treating BX610, Bolatto et al. (2015) assume a conversion factor of $\alpha_{\textrm{CO}}=4.36$M⊙(K km s*-1* pc2)-1, which yields a cold gas mass of (1.1$\pm$0.1)$\times 10^{11}\textrm{M}_{\odot}$. This conversion factor (which includes the contribution from helium) is similar to Galactic values, often used as an estimate for high-redshift main sequence galaxies (Tacconi et al. 2010; Daddi et al. 2010; Genzel et al. 2012).

Rayleigh-Jeans dust continuum has also been shown to be a useful indicator of total gas mass (Santini et al. 2010; Scoville et al. 2014, 2016a). Dust emission in the optically thin Rayleigh-Jeans tail is proportional to the dust mass, and assuming typical values for dust temperature and opacity allows direct calculation of $M_{\textrm{dust}}$ and therefore $M_{\textrm{gas}}$ with an assumed dust-to-gas ratio. We use the formula presented in Scoville et al. (2016b) with our 150 GHz continuum to estimate the gas mass in BX610. We take the recommended values from Scoville et al. (2016a) for the mass-weighted dust temperature, $T$=25 K, and the dust conversion factor $\alpha_{850\mu\textrm{m}}\equiv L_{\nu 850\mu\textrm{m}}/M_{\textrm{mol}}=(6.7\pm 1.7)\times 10^{19}\textrm{erg s}^{-1}\textrm{Hz}^{-1}\textrm{M}_{\odot}^{-1}$. We find $M_{\textrm{mol}}=(2.4\pm 0.5)\times 10^{11}\textrm{M}_{\odot}$, which is significantly higher than the estimated cold gas mass from CO(1-0).

Neutral carbon occurs throughout cold molecular gas, and the optically thin [Ci] transition has been recognized as a useful indicator of molecular gas in galaxies both locally and at high redshift (e.g., Papadopoulos et al. 2004; Bothwell et al. 2017; Valentino et al. 2018). With our detections of both [Ci] lines we are able to directly determine the [Ci] excitation temperature. From Eq. 3 of Walter et al. (2011):

[TABLE]

where $L^{\prime}$ indicates the line luminosity in units of K km s*-1* pc2. We find $T_{\textrm{ex}}=31.8\pm 6.9$K. From this we can determine the neutral carbon mass. From Eq. 4 of Walter et al. (2011):

[TABLE]

where $Q(T_{\textrm{ex}}$) is the temperature-dependent partition function, $Q(T_{\textrm{ex}}$)=1+3$\textrm{e}^{-T_{1}/T_{\textrm{ex}}}$+5$\textrm{e}^{-T_{2}/T_{\textrm{ex}}}$, and $T_{1}$ and $T_{2}$ are the the [CI] transition energy levels above ground state, 23.6 and 62.5 K, respectively. Based on these equations we find $M_{\textrm{CI}}$=(1.4$\pm$0.3)$\times$107 M⊙. To determine a cold gas mass we need to use a conversion factor, $X[CI]/X[H_{2}]\equiv M_{\textrm{CI}}/(6M_{\textrm{H}_{2}})$. This conversion factor may vary with environment, and in particular with the level of metal enrichment. Although there have been fewer studies constraining this ratio across various environments, Valentino et al. (2018) found a ratio of $X[CI]/X[H_{2}]$ in the range 1-13$\times 10^{-5}$ in their sample of z$\sim$1.2 main sequence galaxies with a weighted mean of 1.55$\times 10^{-5}$. This value results in an H2 mass of $(1.5\pm 0.3)\times 10^{11}$M⊙, and a cold gas mass of $(2.0\pm 0.4)\times 10^{11}$M⊙ (introducing an additional factor of 1.36 to account for helium to be consistent with previous estimates).

Our neutral-carbon- and dust-derived gas masses are consistent with each other, although they are both significantly higher than the CO-derived cold gas mass, suggesting that either $\alpha_{\textrm{CO}}$ is higher in BX610 than the value used here, or $\alpha_{850}$ is lower and $X[CI]/X[H_{2}]$ is higher than the values used here. Given the wide range over which $\alpha_{\textrm{CO}}$ is known to vary, and the agreement between dust- and [Ci]-derived masses, we suggest the former case is most likely. If the true cold gas mass of BX610 is given by the weighted mean of the dust- and [Ci]-derived cold gas masses, $(2.1\pm 0.3)\times 10^{11}$M⊙, then the corresponding $\alpha_{\textrm{CO}}$ value is 8.2. This is still within the range of demonstrated values in star-forming galaxies. Taking our weighted mean gas mass and our MAGPHYS-based SFR, we find a gas depletion time of 1.5$\pm$0.3 Gyr consistent with local main sequence star-forming galaxies (e.g., Saintonge et al. 2011; Sargent et al. 2014; Madau & Dickinson 2014).

4.3 Emission line modeling

The integrated CO SLED, comprising our $J_{\rm up}$=4 and 7 observations as well as $J_{\rm up}$=1 and 3 from Bolatto et al. (2015), is shown in Fig. 3. Qualitatively we see that the SLED rises strongly through $J_{\rm up}$=4 and appears to plateau by $J_{\rm up}$=7. Carbon monoxide SLEDs of widely differing shapes have been observed in different sources (e.g., Carilli & Walter 2013; Daddi et al. 2015). Although it is not clear that galaxies of a similar class necessarily have similar CO SLEDs, the presence of strong CO emission at J_{\rm up}$$\geq7 is consistent with intense star formation. The sample of BzK galaxies at z$\sim$1.5 from Daddi et al. (2015) shown in Fig. 3 for example tend to plateau at lower J levels, as do the z$\gtrsim$1 main sequence galaxies in the Hubble Ultra Deep Field detected by Decarli et al. (2016). This contrast suggests BX610 has more hot, dense, and highly excited molecular gas, possibly due to enhanced star formation or other mechanisms.

4.3.1 Large velocity gradient modeling

As a first attempt to characterize the molecular and atomic gas, we fit the CO SLED and the neutral carbon emission lines with large velocity gradient (LVG) models produced by RADEX (van der Tak et al. 2007). Large velocity gradient models are useful tools for determining bulk gas properties such as density, temperature, and chemical abundances, without specifying any particular power source (such as star formation, AGNs, or shocks for example). While most commonly used to model CO SLEDs, LVG models can also include neutral carbon. Under the assumption that the neutral carbon is tracing the same gas and predominantly excited by molecular hydrogen, as is CO, a single LVG model can simultaneously characterize line emission from both CO and [Ci] (see e.g, Israel et al. (2015)). We use a CO abundance per velocity gradient of 10*-5* pc (km s*-1*)-1, a line width of 300 km s*-1*, and a background temperature of 8.75 K(=2.73 K$\times(1+\textrm{z})$) to produce a grid of CO and [Ci] emission line models. Our grid varies the hydrogen density from Log($n$ cm3)=1 to 7 in steps of 0.05, temperature from $T$=10 to 130 in steps of 1 K, and Log([CO]/[C]) abundance ratio from 0.3 to 1.3 in steps of 0.067 (corresponding to [CO]/[C] abundances of 2-20). We compare each emission line model with our observed SLED to find $\chi^{2}$ values. For the full set of models we convert these $\chi^{2}$ values to normalized likelihood values ($\mathcal{L}\equiv\frac{\textrm{e}^{-\chi^{2}/2}}{\sum\textrm{e}^{-\chi^{2}/2}}$). In Fig. 5 we show a corner plot illustrating the likelihood distribution of the LVG parameter space. As seen in the plot, the temperature and density are well constrained with 68% likelihood ranges of $T$=30$\pm$4 K and Log($n$ cm3)=4.3$\pm$0.1, respectively. A [CO]/[C] abundance ratio near the low end is preferred in our fit, with a 68% likelihood range of Log([CO]/[C])=0.47$\pm$0.17. Our best-fit model, with $n$=104.2 cm*-3*, $T$=34 K, and Log([CO]/[C])=0.43 is overlaid on our observations in Fig. 3. This model provides an adequate fit to the full set of emission lines ($\chi^{2}_{\rm reduced}=1.04$).

4.3.2 Photodissociation regions and shocks

Photodissociation regions (PDRs) powered by the UV emission from young, newly formed stars have long been recognized as an excitation source of infrared line emission and PDR models have proven quite useful in describing individual star formation regions locally, as well as integrated emission from star-formation-dominant galaxies. Furthermore, these models have proven effective at differentiating between normal and starburst galaxies based on their density and UV radiation field intensity (e.g., Malhotra et al. 2001; Luhman et al. 2003; Stacey et al. 2010; Brisbin et al. 2015). In the case of BX610, PDR models allow us to simultaneously consider both the CO and [Ci] emission, and also characterize infrared luminosity since much of the UV luminosity from young, newly formed stars is reprocessed by dust into infrared radiation.

Although multiple different PDR models exist with differing levels of complexity, for our purposes the straightforward semi-infinite slab models of the “PDR Toolbox” model suffice (Pound & Wolfire 2008; Kaufman et al. 2006). The parameter space for these models varies gas density from Log($n$ cm3)=1 to 7 in steps of 0.25, and UV radiation (6 eV$<h\nu<$13.6 eV) field intensity (measured in units of the Habing field, G${}_{0}=1.6\times 10^{-3}$ erg cm*-2* s*-1*) from Log($G$ G${}_{0}^{-1}$)=-0.5 to 6.5 in steps of 0.25. At each density and UV intensity, a line flux is given for each of our observed lines (as well as an infrared continuum flux which varies by UV field but not density), which we then scale up or down uniformly by an effective beam-filling factor to most closely match the absolute values of all our observed emission features.

A single component PDR model is not able to adequately fit all four CO lines, both [Ci] lines, and infrared luminosity. In general the observed CO lines require relatively higher densities ($n>$104 cm*-3*) while the [Ci] lines require a lower-density PDR component ($n<$104 cm*-3*). BX610 may be better represented with two PDR components. This might be the case if BX610 has a phase of higher-density star-forming molecular gas emitting high-J CO emission, and a less-dense but widespread component responsible for much of the [Ci] and infrared continuum emission. The best two-component PDR model includes a modest density and UV component, Log($n$ cm3)=3.125$\pm$0.625, Log($G$ G${}_{0}^{-1}$)=3.125$\pm$0.375, which dominates the [Ci] and infrared luminosity, and a higher-density component, Log($n$ cm3)=5.0$\pm$0.5, Log($G$ G${}_{0}^{-1}$)=2.25$\pm$1.0, which dominates the CO SLED. This two-component model suggests a mixture of gas conditions, but with generally low or modest UV fields, typical of local main sequence galaxies rather than starbursts (Malhotra et al. 2001; Luhman et al. 2003). This fit is significantly better than a single-component PDR fit, though it is sufficiently far off to suggest that it does not perfectly explain the emission in BX610 ($\chi^{2}_{\rm reduced}=2.73$). The two-component PDR SLED is plotted in Fig. 3, and the corner plot of its parameter space is shown in Fig. 5. The predicted infrared luminosity is $L_{30-1000}=3.1\times$1012, matching our estimate well.

As noted previously, BX610 has a very high CO(7-6)-to-LFIR ratio, similar to NGC 6240 in which shocks are suspected of powering strong CO emission. Shocks should therefore also be considered as a possible power source for CO emission in BX610. To investigate this possibility we use C shock models from Flower & Pineau Des Forêts (2010). These models give CO line emission over a small grid of shock velocities ranging from 10 to 40 km s*-1* in steps of 10 km s*-1* and at pre-shock densities of 2$\times$104 and 2$\times$105 cm*-3*. We use the shock models in conjunction with PDR toolbox models to simultaneously reproduce the CO SLED. Dust is not efficiently heated through shocks, so we additionally constrain the PDR model with LFIR. Although line emission from neutral carbon could also be enhanced in shocks, the available Flower & Pineau Des Forêts (2010) models do not include [Ci] emission, so we use the observed [Ci] lines to constrain the upper limit of PDR-predicted [Ci] emission. Specifically, a PDR-predicted [Ci] line at or below the observed flux contributes nothing to the overall $\chi^{2}$, and a predicted line flux above the observed line contributes in the normal way, weighted by its residual divided by the line flux uncertainty. We emphasize that this is only a suggestive exercise to check the plausibility of shocks contributing significantly to CO emission. Due to the nonlinear nature of this model, as well as the number of free parameters and the dearth of CO lines which constrain shocks, it has no statistical authority. The best fit shown in Fig. 3 matches the observed CO SLED and infrared luminosity well ($L_{\textrm{30-1000}}=3.1\times 10^{12}$L⊙), with shocks dominating CO emission at J_{\rm up}$$\geq 3. The shock–PDR parameter likelihoods shown as a corner plot in Fig. 6 show well-constrained UV fields, Log($G$ G${}_{0}^{-1}$)=3.25$\pm$0.25, and PDR density, Log($n$ cm3)=3.875$\pm$0.375, which is consistent with our prior LVG analysis. Due to the sparseness of the C shock model grid we cannot robustly estimate shock property uncertainties. Our best-fit model corresponds to a pre-shock density of 2$\times$104 cm*-3* and a shock velocity of 10 km s*-1*, though higher velocities and densities cannot be ruled out.

4.3.3 Physical interpretation

As BX610 is one of very few examples of z$>2$ main sequence galaxies with detailed gas tracer detections, it is worth investigating the nature of its star formation, and comparing this with the other few examples of this population. Recent work has revealed a population of compact star-forming galaxies (cSFGs) in the galaxy main sequence of the early universe characterized by dense star-forming regions within the cores of massive galaxies (e.g., Nelson et al. 2014; Barro et al. 2014; Williams et al. 2014; van Dokkum et al. 2015; Barro et al. 2017; Tadaki et al. 2017; Talia et al. 2018; Puglisi et al. 2019). Although they lie on the main sequence, their star formation densities are more akin to local starbursts, and they are thought to represent a population transitioning from star-forming galaxies to compact quiescent galaxies. Compact star-forming galaxies are generally characterized by broad emission lines (FWHM$\sim$600 km s*-1*), compact emission regions ($r_{\textrm{e}}\sim$ a few kpc), and short depletion times, $\tau_{\textrm{dep}}\sim 100$ Myr. The source 3D-HST GS30274 at z=2.225 has been well studied through PDR-diagnostic lines including several CO and [CI] lines (Popping et al. 2017). Analysis by Popping et al. (2017) established the PDR parameters of the source as $n$=(6$\pm$0.5)$\times 10^{4}$ cm*-3* and $G$=(2$\pm$2)$\times 10^{4}G_{0}$, indicating starburst-like conditions in the compact core. The density and UV field suggested by our PDR+shock model are not particularly extreme; they are consistent with normal star-forming galaxies, though the upper limits of their error ranges do not positively rule out starbursts either (Malhotra et al. 2001; Luhman et al. 2003). Typical local ULIRGs usually have UV fields Log($G$ G${}_{0}^{-1}$)$\geq$3 and densities Log($n$ cm3)$\geq$4, making the PDR component in our PDR+shock model consistent with both the low-density, low-UV field range of ULIRGs as well as the modestly high-density and UV field range of main sequence galaxies. This leaves open the possibility that BX610 may be forming stars in a scaled up “normal” mode of star formation, similar to local non-starburst galaxies, providing an example of an increasingly important population of galaxies in the universe (see also Coppin et al. 2012; Brisbin et al. 2015). Our two-component PDR model does show that a component with high density, comparable to densities observed in cSFGs, might be present in BX610, though even in this case the second component does not perfectly fit the template of cSFGs since this PDR component also has a relatively low intensity UV field. The emission line widths in BX610, $\sim$300 km s*-1*, are smaller than expected in rapidly rotating cSFGs, and the depletion time is longer as well. Furthermore, the effective half light radius for BX610, $r_{\textrm{e}}$, in the H band was found to be 4.4 kpc (Förster Schreiber et al. 2011b), significantly exceeding the selection criteria for cSFGs put forward by van Dokkum et al. (2015), which corresponds to a limiting $r_{\textrm{e}}$ of 2 kpc for BX610. The combined evidence of the larger spatial extent, lower line velocity widths, longer depletion time, and less extreme PDR parameters in BX610 relative to cSFGs leads us to believe that BX610 is not a cSFG.

We have used several models to investigate the observed emission in BX610. The LVG model we use yields a robust estimate of the average gas characteristics in the galaxy. It is also the simplest and the most statistically robust model with a $\chi^{2}_{\rm reduced}$ close to one. Invoking more complicated models is unmerited in that their additional complexity cannot be justified by a corresponding improvement in their fit. This is reflected in the poorer $\chi^{2}_{\rm reduced}$ of our two-component PDR model, and the shock+PDR model which eludes simple $\chi^{2}$ analysis. Nonetheless, examining these models may still yield useful physical insight. The simple LVG model for instance says nothing about what mechanism is heating the gas. Photodissociation region models and shock models impose such mechanisms as well as assumptions which further constrain the resulting fits (PDR models for instance assume a stratified temperature structure rather than a single temperature cloud as in the LVG model).

The very high CO(7-6)-to-$L_{\textrm{FIR}}$ ratio we observe, as well as the consistency of shock models with our CO SLED, suggests shocks may be present in BX610, and raises questions about the nature of these shocks. In local source NGC 6240, with a similarly high CO(7-6)-to-$L_{\textrm{FIR}}$ ratio, a recent merger event likely caused widespread shock-powered turbulence (Meijerink et al. 2013). A similar event may be the cause in BX610 as well, although there is little evidence to confirm or exclude this possibility. Our unresolved observations provide no morphological information on the dust or molecular gas distribution aside from an upper size constraint, and the emission line profiles are not especially broad, nor are they seen to contain multiple velocity components to indicate potential outflows. Förster Schreiber et al. (2014) however find evidence of a weak or obscured AGN based on MIR colors. They also find broad H$\alpha$ line emission spatially offset from the stellar bulge, which is interpreted as coming from star-formation-powered outflows. This offset could be evidence of a merger interaction in the process of relaxation, though differential dust obscuration could also yield an offset in H$\alpha$ emission. Future observations could better disentangle the nature of the shock- and star-formation-excited gas both through improved spatial resolution to look for morphological interaction signatures and by investigating higher-J CO lines and other shocked emission-line diagnostics such as H2O.

5 Summary

We used the Plateau de Bure Interferometer to observe CO(4-3), CO(7-6), Ci, and Ci line emission, as well as the underlying continuum at 150 and 250 GHz in BX610, a star-forming disk galaxy at z=2.2103. Complementing previously detected low-J CO lines, these observations provide insight into a main sequence galaxy in the early universe, which may be representative of a more common population of galaxies than typical starbursting submillimeter galaxies. Our independent measurements of gas mass based on [Ci] and dust continuum are in good agreement with each other, but are significantly elevated above previous estimates based on CO(1-0), suggesting a modestly higher $\alpha_{\textrm{CO}}$ in BX610 of $\sim$8.2. This is particularly important given the paucity of main sequence galaxies observed in more than one cold gas mass tracer at z$>$1.

The SED of BX610, newly extended with submillimeter photometry, suggests star formation rates of 140$\pm$14 M⊙yr*-1* , consistent with prior estimates and confirming the existence of this galaxy on the main sequence. Compared to BzK galaxies at z$\sim$1.5, BX610 has an enhanced mid- to high-J CO SLED. One particularly distinguishing feature of BX610 from most local main sequence galaxies is its apparent strong mid-J CO emission compared to $L_{\textrm{FIR}}$: we find Log($\frac{L_{\textrm{CO}(7-6)}}{L_{\textrm{FIR}}}$)=4.06$\pm$0.06. Such a high ratio is known to be shared with only a few other sources, most famously NGC 6240, and is suspected of being a result of shock-enhanced CO emission. Furthermore, a simple single component PDR model is unable to reproduce the observed high-J CO emission without strongly exceeding the observed neutral carbon and infrared luminosity emission. We find that either a second higher-density PDR component or an additional C shock component are needed to fit the emission lines well. Given the empirical correlation between galaxies with high $\frac{L_{\textrm{CO}(7-6)}}{L_{\textrm{FIR}}}$ ratios and shock excitation, we generally favor the PDR+shock scenario in BX610, although these scenarios should be investigated further with spatially resolved observations and other shock-diagnostic emission lines. The PDR component of the PDR+shock model is consistent with vigorous star formation occurring in BX610, though not as intense as in off-main sequence starbursts which often show significantly stronger UV fields and shorter gas depletion times. The cause for shock-induced CO emission, if present, is currently uncertain. Whether these shocks are the result of recent interactions or galaxy harassment is uncertain, but should be considered in future investigations.

6 Acknowledgements

This work is based in part on observations made with the Spitzer Space Telescope, obtained from the NASA/ IPAC Infrared Science Archive, both of which are operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with the National Aeronautics and Space Administration. D.B. acknowledges support from FONDECYT postdoctorado project 3170974. D.R. acknowledges support from the National Science Foundation under grant number AST-1614213. 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. Este trabajo contó con el apoyo de CONICYT + Programa de Astronomía+ Fondo CHINA-CONICYT.