Probing cold gas in a massive, compact star-forming galaxy at z=6

Observations of low order CO transitions represent the most direct way to study galaxies' cold molecular gas, the fuel of star formation. Here we present the first detection of CO(2-1) in a galaxy lying on the main-sequence of star-forming galaxies at z>6. Our target, G09-83808 at z=6.03, has a short depletion time-scale of T_dep~50Myr and a relatively low gas fraction of M_gas/M_star=0.30 that contrasts with those measured for lower redshift main-sequence galaxies. We conclude that this galaxy is undergoing a starburst episode with a high star formation efficiency that might be the result of gas compression within its compact rotating disk. Its starburst-like nature is further supported by its high star formation rate surface density, thus favoring the use of the Kennicutt-Schmidt relation as a more precise diagnostic diagram for starbursts. Without further significant gas accretion, this galaxy would become a compact, massive quiescent galaxy at z~5.5. In addition, we find that the calibration for estimating ISM masses from dust continuum emission satisfactorily reproduces the gas mass derived from the CO(2-1) transition (within a factor of ~2). This is in line with previous studies claiming a small redshift evolution in the gas-to-dust ratio of massive, metal-rich galaxies. In the absence of gravitational amplification, this detection would have required of order ~1000h of observing time. The detection of cold molecular gas in unlensed star-forming galaxies at high redshifts is thus prohibitive with current facilities and requires a ten-fold improvement in sensitivity, such as that envisaged for the ngVLA.


INTRODUCTION
The evolution of the star formation activity and, consequently, the cosmic mass assembly history and the metal enrichment of the universe are tightly linked to the molecular gas content of galaxies throughout cosmic time (e.g. Walter et al. 2020). How much gas was present in these 'cosmic ecosystems' at different epochs is a fundamental question in galaxy formation and evolution studies.
H 2 , the most abundant molecule in the Universe, lacks dipole moments and therefore is a very inefficient radiator; thus the cool molecular gas' best observable tracer is thought to be carbon monoxide ( 12 C 16 O; hereafter, CO), the second most abundant molecule in the universe. Since a few years after its first astronomical detection -half a century ago (Wilson et al. 1970;Penzias et al. 1971), it was clear that this rotational emission line was a useful tool to constrain the molecular gas content in extragalactic sources and to map its distribution and kinematics (e.g. Rickard et al. 1975;Solomon & de Zafra 1975) since it is easily excited even at very low temperatures.
The estimation of the total molecular gas from the CO line luminosity relies on the CO-to-H 2 conversion factor, α CO , usually calibrated for the CO(1 → 0) transition 1 (see review by Bolatto et al. 2013). When only higher order transitions are available, an assumption on the excitation ladder (or a modeling if several lines are available) has to be done as an additional step to first infer the CO ground transition line luminosity. The shape of the spectral line energy distribution (SLED), however, depends on the excitation of the CO molecules (and therefore on the temperature and density of the gas), which might vary significantly between galaxies and even between regions of the same galaxy. This introduces significant uncertainty when inferring the low-J line luminosities from high-J transitions, at the level of making gas mass estimations unreliable, particularly when only one high-order line is available 2 . This is particularly true at z > 1 given the sensitivity and frequency coverage of current facilities and the typical SLED of star-forming galaxies, which bias the observations towards the brighter high-J CO emission lines.
Despite these difficulties, hundreds of CO studies in galaxies from z = 0 to z ≈ 3 − 4 have been published to-date, allowing us to constrain the physical properties and redshift evolution of galaxies' interstellar medium (see recent review by Tacconi et al. 2020 and references therein). Moreover, these observations (in addition to other dynamical studies) have been used to calibrate alternative tracers of the molecular gas component such as far-infrared (FIR) atomic fine structure lines (e.g. Valentino et al. 2018;Dessauges-Zavadsky et al. 2020) and dust continuum emission (e.g. Magdis et al. 2012;Scoville et al. 2016).
At higher redshifts, however, our understanding of the physical conditions of the cold star-forming interstellar medium is much more limited. And, while recent observational programs, such as the ALPINE survey Le Fèvre et al. 2020) have made a significant leap forward in our understanding of these quantities up to z ∼ 5.5 (e.g. Dessauges-Zavadsky et al. 2020), primarily through observations of [CII], the number of galaxies with low-J CO detections at these high redshifts is limited to a handful (e.g. Strandet et al. 2017;Pavesi et al. 2018). At z > 6, the situation is even more critical. The only detections of CO(1 → 0) and CO(2 → 1) come from a very extreme starburst galaxy, whose star formation rate is thought to be among the highest observed at any epoch (Riechers et al. 2013). The gas properties of normal star-forming galaxies at z > 6 are thus completely unexplored.
The reason behind this relies mainly on the limitations of existing instrumentation since the current facilities covering the radio wavebands at which these low-J CO transitions are redshifted (ν obs ≈ 15 − 35 GHz) lack the required sensitivity to detect the expected CO emission from typical star-forming galaxies at z 5. Future large radio arrays are currently being designed to specifically address this issue. The Next-Generation Very Large Array (ngVLA), for example, will provide an order-of-magnitude improvement in depth and area (Carilli et al. 2015;Murphy et al. 2018;Selina et al. 2018) allowing us to conduct surveys of cold gas from these critical low-J transitions in normal galaxies within the Epoch of Reionzation (Casey et al. 2015Decarli et al. 2018).
Nevertheless, such facilities are not expected to be operational until the end of this decade, at the earliest. Furthermore, the capability of such a facility to detect cold gas at high redshifts (z > 5) has been called into question due to the potentially strong dimming effects caused by the Cosmic Microwave Background (e.g. da Cunha et al. 2013).
Here, we exploit the gravitational amplification effect, to demonstrate the feasibility of detecting cold gas in normal galaxies within the Epoch of Reionization via the first detection of CO(2 → 1) in a main-sequence star-forming galaxy at z = 6. The detection of this line, from which the ground CO(1 → 0) line luminosity can be estimated with little extrapolation, enable us to characterize the cold ISM in a galaxy formed < 1 Gyr after the Big Bang. This allows us to test the applicability of scaling relations derived for lower redshift galaxies and test whether or not other widely used methods to infer gas masses are also valid for this z = 6 system.
This work provides us the first insights of what is expected to come with future facilities such as the ngVLA (Carilli et al. 2015) and the new Band-1 ALMA detectors (Huang et al. 2016).
The target, G09-83808, was first detected by the Herschel space telescope and identified as a high-redshift galaxy candidate in Ivison et al. (2016) Figure 1. G09-83808 in the context of the main-sequence of star-forming galaxies. The best-fit SFR and stellar mass of G09-83808, and their associated uncertainties, are represented by the red solid circle. The gray and gold regions represents the main-sequence relationships of Schreiber et al. (2015) and Pearson et al. (2018), respectively, extrapolated to z = 6. For the sake of comparison we also include z ≈ 6 star-forming galaxies from the IllustrisTNG-300 simulation (Pillepich et al. 2018). Our target, G09-83808, lies on the high-mass end of the main-sequence over a parameter space that overlaps with the the most massive star-forming galaxies from IllustrisTNG. red far-infrared colors (i.e. S 250µm < S 350µm < S 500µm ). The source was then followed-up with several telescopes including the Atacama Large submillimeter/Millimeter Array (ALMA), the James Clerk Maxwell Telescope (JCMT), the Large Millimeter Telescope (LMT), NOEMA, Spitzer, and the Submillimeter Array (SMA). The galaxy's redshift was determined to be z = 6.0269±0.0006 through the detection of multiple emission lines (including CO(5 → 4) , CO(6 → 5) , and H 2 O(2 11 − 2 02 )) in its LMT 3 mm spectrum and a subsequent detection of the [CII]( 2 P 3/2 − 2 P 1/2 ) transition with the SMA (Zavala et al. 2018; see also Fudamoto et al. 2017). Zavala et al. (2018) also presents highangular resolution ALMA observations of the dust continuum emission at λ obs ∼ 890 µm. The high resolution dust continuum observations were used for modelling the gravitational lensing effect in the uv plane using the visilens code (Spilker et al. 2016b). The best-fit gravitational magnification of the source was found to be µ 890µm = 9.3 ± 1.0.
The 3.6 and 4.5 µm Spitzer/IRAC observations were also used along with the FIR photometry to constrain the stellar mass of the galaxy through a spectral energy distribution (SED) fitting technique using the energy balance code magphys (da Cunha et al. 2008). First, since the emission of G09-83808 is blended with that from the foreground z = 0.776 lensing galaxy in the IRAC bands (see Zavala et al. 2018), the light distribution of the foreground galaxy was modeled using galfit (Peng et al. 2002) and a Sérsic profile (Sérsic 1963). Then, the emission from the foreground galaxy was subtracted from each image, and the photometry of the background source was measured from the galfit generated residual images using SExtractor (Bertin & Arnouts 1996). Finally, combining the deblended Spitzer photometry, which probes the rest-frame optical stellar emission, with the FIR data, the best-fit SED for G09-83808 was derived, from which a stellar mass of M = 7.8 +8.4 −4.2 × 10 10 M was inferred (after correcting for gravitational magnification).
With these measurements in hand, we can place our target in the context of the main-sequence of star forming galaxies.
As it can be seen from Figure 1, G09-83808 lies on the high-mass end of the main-sequence of star forming galaxies when compared to the extrapolated relation of Schreiber et al. (2015), or slightly above (with SFR/SFR MS ≈ 2 − 3) if compared to the relation of Pearson et al. (2018) or Khusanova et al. (2021) -note that the cuts at high SFRs in the figure are imposed to represent the typical turnover at high masses (e.g. Tomczak et al. 2016).
Since this M −SFR relationship is still not well determined at z ∼ 6, we make use of results from simulations to further explore the place of our galaxy in the context of typical star-forming galaxies. To do this, we plot in Figure 1 the z ∼ 6 star-forming galaxies from the Il-lustrisTNG project (Pillepich et al. 2018). These galaxies are in relatively good agreement with the adopted main-sequence relationships, ruling out any significant bias in our comparison. As shown in this figure, the properties of our target are similar to those of the most massive galaxies in IllustrisTNG. Therefore, we can conclude that G09-83808 probes the high-end of the main sequence of star-forming galaxies.
This contrasts with the only other two z > 6 galaxies discovered in blind (sub-)millimeter surveys, SPT0311-58 and HFLS3, which show SFRs of ∼ 1, 000 s of M yr −1 and are considered to be extreme starburst galaxies (Riechers et al. 2013;Marrone et al. 2018).
Finally, we highlight that the metallicity of this galaxy has recently been constrained to be Z ≈ 0.5 − 0.7 Z Right ascension  Figure 2. The CO(2 → 1) transition from G09-83808. Left: Velocity-integrated intensity (moment-zero) map of the CO(2 → 1) line over ≈ 800 km s −1 centered at 32.81 GHz. A double-arc structure is visible due to the gravitational lensing. The aperture used to extract the 1D spectrum is indicated with the dashed black line, while the beam-size of the observations is shown with the yellow ellipse in the bottom left. Right: Extracted spectrum around the expected frequency of the line at z = 6.029 (gray dashed line). The dotted red line represents the best-fit Gaussian function (with ν obs = 32.810 ± 0.002 GHz and FWHM = 0.040 ± 0.005 GHz). The solid line represents the best-fit with two Gaussian components adopted to better describe the double-peak profile of the line, which is also noticeable in other transitions (see Appendix A).

VLA observations
Observations were taken using the Karl G. Jansky Very Large Array (VLA) in the C array configuration as part of project 20A-386 (PI: J. Zavala). Three different executions were performed on 2020 June 10, 14, and 19 for a total integration time of 15 h.
The WIDAR correlator setup was designed for simultaneous continuum and spectral line observations in the Ka-band, with mixed 3-bit and 8-bit samplers, resulting in a total bandwidth of 5.12 GHz subdivided into forty 128 MHz dual-polarization sub-bands with 1 MHz channels. During each execution, the source J1331+305 (aka 3C-286) served as band-pass and absolute flux calibrator, while J0909 + 0121 was used as pointing source and complex gain calibrator. A few scans from some antennas were flagged before data calibration due to phase or amplitude issues or because they were affected by radio frequency interference (although they do not represent more than a few percent of all data). Data reduction and calibration were done using the VLA pipeline following the standard procedures. Then, the data from the three executions were combined during the imaging procedure, which was done using natural weighting of the visibilities in order to maximize the sensitivity at the expense of angular resolution (producing a synthesized beam size of 0.97 × 0.68 , PA= −30 deg). This results in an r.m.s noise level of ≈ 60 µJy/beam for a ∼ 45 km s −1 channel width for the sub-band centered at the expected position of the CO(2 → 1) line.

ANALYSIS AND RESULTS
3.1. The CO(2 → 1) line and the lensing model Figure 2 shows the CO(2 → 1) line detection (moment-0 map and 1D extracted spectrum) from G09-83808. To measure the total flux density of the line we fit either one or two Gaussian profiles to the extracted spectrum (see Figure 2), and we measure it directly from the clean moment-0 map. These methods give us a statistically consistent integrated line flux of µS CO(2→1) = 0.6 ± 0.1 Jy km s −1 , which implies a line luminosity of µL CO(2→1) = 1.8 ± 0.3 × 10 11 K km s −1 pc 2 (following Solomon et al. 1992).
To measure the gravitational amplification factor on the line, we use the lens modelling code visilens (Spilker et al. 2016b), which directly models the visibilities in the uv plane. The modelling was done in a similar way as for the 890 µm dust continuum emission previously reported in Zavala et al. (2018), parameterizing the lens mass profile as a singular isothermal ellipsoid and the background source as a n = 1 Sérsic profile. While this modelling provides a good fit to the data (see Figure 3), we acknowledge that the CO(2 → 1) data alone cannot constrain the source profile. The as- sumed n = 1 Sérsic profile is however supported by the modelling of high-SNR, high-angular resolution observations of the dust continuum emission (Zavala et al. 2018;Tadaki et al. 2022. After marginalizing over the different parameter space explored by the code, we derive a best-fit magnification factor of µ = 12 ± 3 and a source-plane (intrinsic) size of R circ = 0.6 ± 0.2 kpc (0.11 ± 0.03 ). These values are in good agreement with the results presented in Zavala et al. (2018), who derived a magnification of 9.3 ± 1.0 and a size of R circ = 0.6 ± 0.1 kpc for the 890 µm dust continuum emission, and with Tadaki et al. (2022) who found µ = 8.4 +0.7 −0.3 and R eff ≈ 0.65 kpc for the 1.5 mm continuum, despite using a different lensing code. This suggests that there is no significant differential magnification between the CO(2 → 1) and the dust thermal emission and implies similar and co-spatial source plane areas.

The effect of the CMB
The systematic increase of the temperature of the Cosmic Microwave Background (CMB) with increasing redshift (T CMB ∝ T z=0 CMB (1 + z)) affects the physical conditions and detectability of galaxies' dust and molecular gas emission. As detailed in da Cunha et al. (2013), these effects are more important when the CMB temperature (≈ 19 K at z = 6) is close to the dust or gas temperature (T dust or T kin , respectively). Harrington et al. (2021) and Jarugula et al. (2021) have recently reported high kinetic temperatures (T kin /T d ≈ 2 − 3) for dusty star-forming galaxies. If this holds at higher redshifts, it would imply that the CMB effects on the measured CO line fluxes are relatively minor (note that the CMB also affects the kinetic and dust temperatures; da Cunha et al. 2013). G09-83808, the source studied in this work, indeed shows a bright CO(12 → 11) detection as expected for a high kinetic temperature. Nevertheless, a detailed modelling of its SLED (including the CO(2 → 1), CO(5 → 4), CO(6 → 5), and CO(12 → 11) transitions) suggests that the CO(2 → 1) line is almost totally dominated by a relatively cold component with T kin = 55 ± 16 K (Tsujita et al. submitted).
Assuming this kinetic temperature and following da Cunha et al. (2013), we infer a CMB correction factor of 1.14 +0.04 −0.06 (note that we used the dust temperature as the minimum kinetic temperature in this calculation). The intrinsic CO(2 → 1) flux density is then estimated to be 0.06 ± 0.02 Jy km s −1 , which translates into an intrinsic line luminosity of L CO(2→1) = (1.7±0.5)×10 10 K km s −1 pc 2 . These values, which have been corrected by the gravitational amplification, will be adopted for the rest of the analysis.

Total molecular gas content
The CO(2 → 1)-to-CO(1 →) line ratios reported in the literature show a small scatter of less than a factor of 2 (e.g. Harrington et al. 2021). This contrasts with the high dispersion in the line ratios between the CO(1 → 0) ground transition and higher-J lines, which could differ by up to a factor of ∼ 6 (see the CO SLED of DSFGs reported in Casey et al. 2014; see also Carilli & Walter 2013). This implies that the CO ground state luminosity, and thus the molecular gas mass, can be inferred with relatively little extrapolation from the CO(2 → 1) detection.
As mentioned in Section 1, the total molecular gas mass can be directly inferred from the CO(1 → 0) line luminosity via the CO-to-H 2 conversion factor, α CO . However, this conversion factor suffers from large uncertainties and might depend on several physical parameters of a galaxy's ISM (e.g. Bolatto et al. 2013). A value of α CO = 1.0 M K km s −1 pc 2 is typically adopted for nearby nuclear starburst galaxies 3 (i.e. ULIRGs; Downes & Solomon 1998), while a higher value of α CO = 4.6 M K km s −1 pc 2 is thought to be more appropriate for more moderate galaxies (including the Milky Way; see review by Bolatto et al. 2013). These two extreme values would bracket the molecular gas of G09-83808 between ∼ 2 × 10 10 M and ∼ 9 × 10 10 M .
Very high α CO could, however, lead to a molecular gas mass estimate larger than the dynamical mass (e.g. Bryant & Scoville 1999). Such dynamical information can thus be used to constrain the CO-to-H 2 conversion factor in individual systems (e.g. Bothwell et al. 2010;Engel et al. 2010). In this work we exploit the dynamical information encoded in the CO(2 → 1) emission line that traces the cold gas reservoir to restrict the CO-to-H 2 conversion factor of our target.
The dynamical mass is estimated by: where R e is the effective circularized radius, σ is the velocity dispersion, and G the gravitational constant (e.g. Casey et al. 2019). To correct for the unknown inclination, i, we multiply by a factor of 3/2 (the reciprocal of the expectation value of sin 2 (i)). Based on this equation, we infer a dynamical mass of M dyn = (2.6 +2.4 −1.5 )×10 10 M (note that a similar value is obtained if we use instead the isotropic viral estimator equation, e.g. Engel et al. 2010). This places an upper limit on the conversion factor of α CO < 2.5 M K km s −1 pc 2 (taking into account 1σ uncertainties), which justifies the adoption of α CO = 1.0 M K km s −1 pc 2 , as typically measured for ULIRGs (Bolatto et al. 2013). The total de-lensed molecular gas mass of G09-83808 is then estimated to be M H2 = (2.0 ± 0.7) × 10 10 M .

A Test on the dust continuum method
3 Historically, the ULIRG value was adopted to be α CO = 0.8 M K km s −1 pc 2 , nevertheless, a value of α CO = 1 M K km s −1 pc 2 aligns better with the more recent review by Bolatto et al. 2013. Scoville et al. 2016 presented a relationship between the specific luminosity at rest-frame 850 µm, L ν (850µm), and the CO line luminosity, L CO(1→0) . Such a relation serves as the base for an empirical calibration aimed at inferring galaxies' molecular gas content via dust continuum observations that probe the Rayleigh-Jeans tail of dust black-body emission. While this time-efficient approach has been revolutionary for the measurement of ISM masses -because dust continuum is observationally cheaper than spectral observations of CO transitions, its calibration and validity have not been tested beyond z ∼ 4 (even when studies using this method extends up to z ∼ 6; e.g. Liu et al. 2019).
Thanks to the CO(2 → 1) detection in G09-83808, here we test whether or not the L ν (850µm)-L CO(1→0) calibration holds at z = 6 for this system. Note that we decide to test this relation rather than the final molecular gas recipe of Scoville et al. 2016 to avoid possible differences in the α CO assumptions.
To mimic the typical procedure adopted in the literature, we calculate the 850 µm specific luminosity from one single data-point close to the Rayleigh-Jeans tail (in this case, the gravitationally-corrected 1.1 mm flux density of S 1.1mm = 2.2 ± 0.3 reported in Zavala et al. 2018). This calculation assumes a dust emissivity index of β = 1.8 and a mass-weighted dust temperature of 25 K (see details in Scoville et al. 2016).
As mentioned above, the alternative approach typically adopted in the literature to estimate the CO(1 → 0) line luminosity from high-J transitions might result in larger discrepancies than the one observed here. Therefore, although larger samples are imperative to derive any firm conclusions, this result supports the empirical approach of using single band continuum observations as a proxy for molecular gas content up to z ∼ 6 (modulo an appropriate α CO ) 4 . We highlight, though, that the CMB effects could be more pronounced in galaxies at higher redshift and/or with colder temperatures.
An important corollary that can be inferred from this, is that the gas-to-dust abundance ratio at z ∼ 6 is similar to that found in lower redshift galaxies, at least for this massive (and relatively metal-rich) galaxy. This is in line with some recent simulations which predict a modest redshift evolution on the gas-to-dust ratio (e.g. Li et al. 2019;Popping & Péroux 2022).

Gas depletion timescale, gas-to-dust ratio, and gas fraction
With a gas mass measurement in hand, we can now explore other parameters such as the gas-to-dust ratio and gas fraction, and the gas depletion timescale, and compare them to those measured in other galaxies. In addition, we can test whether or not the scaling relations found for lower redshift systems still hold at z ∼ 6. This is the main goal of this section.
For this comparison we use the recent compilation by Tacconi et al. (2020) of ∼ 2, 000 galaxies with gas measurements and their derived relationships, which are built upon previous significant efforts (e.g. Genzel et al. 2015;Scoville et al. 2017;Tacconi et al. 2018). Although this collection includes galaxies up to z ∼ 5, the low number of galaxies at high redshifts limits the analysis and the derived relationships to z 4. Hence, our observations allow us to investigate the evolution of these parameters within an unexplored redshift range and test the validity of the extrapolated scaling relations at z ∼ 6.
It is important to highlight that Tacconi et al. (2020) used a CO conversion factor of α CO = 4.36 M K km s −1 pc 2 and a gas-to-dust ratio of δ GDR = 67 to re-calibrate all the gas mass measurements from the literature homogeneously. This is justified since their analysis is focused on galaxies around the mainsequence (∆MS ≡ log(SFR/SFR MS ) : [−1, 1]), for which little to no systematic variation of these parameters have been found (see detailed discussion in Genzel et al. 2015;Scoville et al. 2017;Tacconi et al. 2020). Nevertheless, this has been shown not to be the case for extreme outliers with log(SFR/SFR MS ) > 1 (which represent only ∼ 10 % of their sample). Therefore, for those extreme galaxies with gas mass estimates based on α CO = 4.36 M K km s −1 pc 2 , we have scaled the gas masses down by a factor of 4.36, and propagated this to other measurements, to reflect our assumption of α CO = 1 M K km s −1 pc 2 . The appropriateness and the implications of this choice are discussed below.

Gas-to-dust ratio
To calculate the gas-to-dust ratio, we adopt the amplification-corrected gas mass reported above in §3.2 and the demagnified dust mass of M d = (1.9 ± 0.4) × 10 8 M reported by Zavala et al. (2018). Using these values, we estimate a gas-to-dust ratio of δ GDR = 105 ± 40, which is in excellent agreement with the widely adopted ratio of 100:1 for massive galaxies with close to solar metallicities (e.g. Leroy et al. 2011;Rémy-Ruyer et al. 2014) and consistent within 1σ with the value of ∼ 67 adopted in the recent review by Tacconi et al. (2020).
Given that G09-83808 was shown to have a near-solarmetallicity ISM (Tadaki et al. 2022), this result supports the conclusion by Rémy-Ruyer et al. (2014) who claimed that the metallicity is the main physical property driving the gas-to-dust ratio (see also Li et al. 2019;Popping & Péroux 2022), and it suggests that this relation does not significantly evolve with redshift (c.f. Saintonge et al. 2013).
In addition, this supports our assumption of α CO = 1 M K km s −1 pc 2 for the gas estimation. Adopting a higher value (as the typical Milky Way value) would not only result in a gas mass in excess of the dynamical mass, but also in an relatively high gas-to-dust ratio of around δ GDR ∼ 500.

Depletion time
The depletion timescale, which expresses the time in which the molecular gas reservoir would be consumed at the current SFR (in the absence of inflows or outflows of molecular gas), is calculated as τ depl = M mol /SFR (note that the inverse of this quantity is usually referred to as the star formation efficiency, SFE = 1/τ depl ). For this calculation we assume the molecular gas mass derived above ( §3.2) and the SFR of 380 ± 50 M yr −1 reported by Zavala et al. (2018). Combining these measurements we infer a gas depletion time of τ depl = 50 ± 20 Myr for G09-83808.
In Figure 4 we compare the gas depletion time of G09-83808 with the large compilation of Tacconi et al. (2020) and their scaling relation. Tacconi et al. (2020) found that the integrated depletion time scale depends mainly on redshift and offset from the main-sequence (τ depl ∝ (1 + z) −1 × ∆MS −0.5 ). This trend has been confirmed to extend up to z ∼ 5.5 thanks to the ALPINE survey (Dessauges-Zavadsky et al. 2020). The inferred depletion time of our target is indeed consistent with this trend of decreasing depletion time with redshift, although significantly shorter than the expected relation for main-sequence galaxies. Surprisingly, its depletion of τ depl = 50 ± 20 Myr is in better agreement with the extrapolated relation for starburst galaxies (see Figure  4), despite being located on the main-sequence of starforming galaxies (as discussed in §2.1).
To extend our comparison to other high-redshift dusty star-forming galaxies, we also include in Figure 4 the seventeen sources reported in Aravena et al. (2016), which cover a redshift range of z ≈ 2.5 − 5.5 with a median SFR of ∼ 1, 000 M yr −1 . All these galaxies, including G09-83808 and the most extreme galax- The derived depletion timescale of G09-83808, represented by the red solid circle, is shown against the field scaling relation derived by Tacconi et al. (2020) and their compiled sample (see §3.3 for more details). Blue circles represent main-sequence star-forming galaxies (log(SFR/SFRMS) < 1) in such compilation, where the darkness and size represent a proxy for the number density. Open red circles represent their starburst galaxies (log(SFR/SFRMS) ≥ 1) with gas masses scaled down to adjust for αCO (see §3.3). The solid and dashed lines are the best-fit relations for main-sequence (log(SFR/SFRMS) = 0) and starburst galaxies (log(SFR/SFRMS) ≈ 1.2) , respectively, extrapolated to z ∼ 7. We also include the sample of dusty star-forming galaxies of Aravena et al. (2016) selected with the SPT (yellow triangles) and highlight the locus occupied by the ALPINE main-sequence galaxies (Dessauges-Zavadsky et al. 2020; gray regions). While the scaling relation successfully captures the behavior towards shorter depletion timescales for starbursts, a significant fraction of them scatter below the relationship, suggesting even shorter depletion times (or higher star formation efficiencies). Right: Analog to the left panel (with same symbols and colors) but for the gas-to-stellar mass ratio. For the SPT sample, we only plot those sources with stellar mass measurements in Ma et al. 2015. While some z ≈ 1 − 2 starbursts have high gas fractions consistent with the best-fit scaling relation, G09-83808, the SPT galaxies, and most of the local (z ∼ 0) ULIRGs do not show enhanced gas fractions with respect to main-sequence galaxies.  Figure 4). By looking carefully at Figure 4, one can realize, though, that a significant fraction of these extreme galaxies lie even below this starburst relationship. This is true for almost all of the z ∼ 0 ULIRGs and the SPT galaxies, which suggests that all these extreme galaxies might have even higher star formation efficiencies. Aravena et al. (2016) concluded that most of the galaxies in their sample are indeed experiencing a starburst phase likely triggered by major mergers. This triggering mechanism could thus explain the shorter depletion times (i.e. higher star formation efficiencies) of these galaxies compared to galaxies around the main-sequence. Similarly, local ULIRGs and other outlier objects are also known to be almost invariably mergers with very compact nuclei (Sanders et al. 1991), which supports the existence of a starburst mode of star formation (e.g. Sargent et al. 2012;Silverman et al. 2015).
The double-peak profile in the lines of G09-83808 (see Figure 2 and Appendix A) might indeed be indicative of a merger activity, as in the case of the core of Arp 220 (Downes & Solomon 1998) and other DSFGs (e.g. Bothwell et al. 2013). Nevertheless, Tsujita et al. (submitted) showed, based on a source-plane reconstruction of the [OIII]88µm and [NII]205µm emission lines across different velocity channels, that this galaxy is well-modelled by a rotating disk with some compact star-forming clumps. It is thus likely that the high star formation efficiency of this galaxy and its starburst-like properties are the consequence of gas compression induced by its compact size (e.g. Elbaz et al. 2011;Gómez-Guijarro et al. 2022).

Gas fraction
Here, we explore the gas fraction (in terms of the gasto-stellar mass ratio) of our target and place it in the context of previous studies.
The redshift evolution of this quantity, and its dependency with other parameters, was also studied and parameterized in the recent review by Tacconi et al. (2020). They concluded that the redshift evolution of µ mol = M mol /M is better described by a quadratic function, with a steep rise between z = 0 and z ∼ 2 − 3, the redshift at which the galaxies show the largest molecular fractions (at fixed SFR/SFR MS ), and with a subtle turnover at higher redshifts (Figure 4). This flattening (or slight turnover) at high redshift is also supported by the recent results of Dessauges-Zavadsky et al. (2020), who measured an average gas fraction of µ mol ≈ 60% over z ≈ 4.5 − 5.5 (although with a relatively large scatter driven likely by galaxies' stellar masses, as can bee seen in the figure). To lesser degree, this parameter also depends on the vertical location of a galaxy along the main-sequence plane as seen in Figure 4. Such a relation implies that the increase in galaxies' SFRs with redshift (at fixed SFR/SFR MS ) and with offset from the main sequence (at fixed redshift) is due primarily to increased gas content, and secondarily to an increased efficiency for converting gas to stars -at least for those 'normal' galaxies around the main-sequence (the focus of the Tacconi et al. 2020 study). Figure 4 reveals that some z ∼ 1 − 2 starbursts do also have high gas fractions, suggesting that their high SFRs could also be the result of having a large gas mass reservoir, as seen for modestly star-forming galaxies. Nevertheless, most of the local ULIRGs do not show enhanced gas fractions. The same is true for the few SPT-selected galaxies with stellar mass measurements (Ma et al. 2015), whose gas-to-stellar mass ratios are in agreement with 'normal' star-forming galaxies (see Figure 4). Starburst galaxies could thus have a large range of gas fractions, implying that their high SFRs are not always caused by their high gas masses but also by other mechanisms that enhance their star formation efficiencies, as discussed before in the literature (e.g. Scoville et al. 2017).
The galaxy studied in this work has, interestingly, a relatively low gas fraction of µ mol = 0.26 +0.29 −0.16 (see red solid circle in Figure 4), comparable to those of the SPT sample. This supports the scenario discussed above: G09-83808 shows a starburst mode of star formation with a high star formation efficiency driven likely by its internal properties such as gravitational instabilities or gas compression.

An alternative scenario with a Milky Way-like CO-to-H 2 conversion factor
While a low ULIRG-like CO-to-H 2 conversion factor (α CO ∼ 1 M K km s −1 pc 2 ) is favored by the current analysis thanks to the limits placed by the dynamical mass, we cannot totally discard the possibility of a higher value, particularly if the system is experiencing a merger (since the dynamical mass assumes a virialized system). The source-plane reconstruction of this galaxy is well-modelled with a disk profile though (see §3.1 and Tadaki et al. 2022 who used higher angular resolution observations with a physical spatial resolution below 1 kpc in the source plane). In addition, the [OIII]88µm and [NII]205µm emission lines show evidence of a monotonic velocity gradient consistent with a rotating disk (Tsujita et al. submitted). This backs the applicability of a dynamical estimator for the total gas of G09-838083 and the inferred upper limit on the conversion factor of α CO < 2.5 M K km s −1 pc 2 .
Assuming a larger value of α CO would not only imply a molecular gas mass exceeding the dynamical mass but also a gas-to-dust ratio of ∼ 500 (if assuming the Milky Way conversion factor reported in Bolatto et al. 2013); at odds with the typical values measured for massive and metal-rich star-forming galaxies. A galaxy with δ GDR = 500 would be expected to have a relatively low metallicity in the range of Z ≈ 0.15 − 0.35Z (according to the δ GDR -metallicity relations of Rémy-Ruyer et al. 2014;Tacconi et al. 2018;Boogaard et al. 2021). This contrasts with the metallicity of G09-838083 which was constrained to Z ≈ 0.5 − 0.7 Z by Tadaki et al. (2022). Therefore, given the current constraints, the possibility of G09-83808 having a large Milky Way-like CO-to-H 2 conversion factor seems unlikely.

DISCUSSION AND CONCLUSIONS
Taking advantage of gravitational amplification, here we present the first detection of CO(2 → 1) in a galaxy lying on the main-sequence of the M − SFR plane and the characterization of its cold ISM and star formation.
While G09-83808 lies around the scatter of the mainsequence of star-forming galaxies at z ∼ 6 (see Figure  1), it shows properties that resemble those from local ULIRGs and other starbursts, including a compact (dust and gas) size of R circ = 0.6 ± 0.2 kpc ( §3.1), a relatively high dust temperature (T d = 49 ± 3 K; Zavala et al. 2018), a low α CO of < 2.5 M K km s −1 pc 2 ( §3.2), and a short depletion time (τ dep = 50 ± 20 Myr; §3.3.2). All this suggests that galaxies with starburst-like ISM properties and star formation modes could also be hidden within the main-sequence population (see also Gómez-Guijarro et al. 2022).
The starburst-like nature of G09-83808 is further revealed by the Kennicutt-Schmidt relation ( Figure 5). This galaxy shows a high star formation rate surface density (Σ SFR = 170 ± 110M yr −1 kpc 2 ) in agreement with what is expected for merger/starburst galaxies, according to the relation of Genzel et al. 2010 (see dashed black line in the figure). Previous studies have found that galaxies with high Σ SFR show shorter depletion The molecular and star formation rate surface density of our z = 6 target, G09-83808, is shown in comparison to other galaxies (following the same color code as in Figure 4). Note that for the Tacconi et al. (2020) compilation we use the reported optical effective radius as a proxy for the size of the star forming regions and the molecular reservoirs (and our scaled gas masses for the starburst galaxies; see §3.3).
For the SPT sample we use the dust-based sizes from Spilker et al. 2016b. Despite the assumptions on the sizes, all these starburst galaxies, including G09-83808, lie on the Genzel et al. (2010) relation for merger/starburst galaxies (dashed black line), while the main-sequence galaxies lie below in better agreement with the relation for 'normal' star-forming galaxies (black solid line). This diagram reveals the starburst nature of our target, which is not captured by the mainsequence plane.
timescales, smaller gas fractions, warmer dust temperatures, and high excitation conditions (e.g. Narayanan & Krumholz 2014;Franco et al. 2020;Burnham et al. 2021;Puglisi et al. 2021) in line with the properties of our target. We thus hypothesize that the star formation rate surface density can be a better indicator of starburst-like galaxies than the main-sequence offset (or at least, alternative to). Despite its starburst-like nature, this high-redshift ULIRG analog, shows a relatively low gas fraction of M gas /M = 0.29 +0.33 −0.17 (see Figure 4 and Section 3.3.3), implying that its starburst episode is not produced by an increased gas supply, but rather by a physical triggering mechanism that enhances its star formation efficiency. Typical mechanisms cited in the literature to explain this starburst mode of star formation include major/minor mergers, compactness, and disk instabilities. While higher angular resolutions are required to precisely determine the causes of the high star formation efficiency in G09-83808, we attribute it to its internal properties, such as compactness and/or gravitational instabilities which facilitate gas compression (see further discussion in Tsujita et al. submitted, who found evidence of star-forming clumps within the rotating disk of this system). This is in line with recent results from Gómez-Guijarro et al. (2022), who found that the most compact galaxies show the shortest depletion time and a starburst-like modes of star formation.
Without further significant gas accretion, we estimate, following Casey et al. (2021), that G09-83808 will finish its star formation episode at z ≈ 5.6 − 5.7 with a final stellar mass on the order of 1 × 10 11 − 3 × 10 11 M . This is well aligned with the derived masses of massive quiescent galaxies discovered at high redshifts (e.g. Tanaka et al. 2019;Marsan et al. 2020;Valentino et al. 2020;Santini et al. 2021;Stevans et al. 2021), which also show compact sizes similar to the size of our target. This hints at a progenitor-descendant scenario between these populations of galaxies, and support other results from the literature suggesting that compact galaxies with low gas fractions and short depletion times are likely in the latest stages of a starburst episode and on the verge of quenching (e.g. van Dokkum et al. 2015;Spilker et al. 2016a;Puglisi et al. 2021;Gómez-Guijarro et al. 2022;Ikarashi et al. 2022).
The determination of the molecular gas mass via direct detection of a low-J CO line transition also allows us to test whether or not other widely used methods to infer gas masses are also valid at z = 6. We found that the relationship between the rest-frame 850 µm specific luminosity and CO line luminosity (L ν (850 µm ) − L CO(1→0) ), on which the Scoville et al. 2016 calibration for estimating ISM masses from dust continuum emission is built on, satisfactorily predicts the CO(1 → 0) line luminosity for this system within a factor of ∼ 2. Therefore, single band continuum observations seem to be a promising way to constrain galaxies' gas masses even at this early epoch (modulo the assumption of a similar α CO ). A corollary that could follow this result is that the gas-todust ratio, δ GDR , in massive (metal-rich) star-forming galaxies does not significantly evolve with redshift, as predicted by recent simulations (e.g. Li et al. 2019;Popping & Péroux 2022). This is independently confirmed for our target using the gas and dust estimates, from which we derive a δ GDR = 105 ± 40. This value is also in relatively good agreement with the expected δ GDR given its high metallicity of Z ≈ 0.5 − 0.7 Z (if assuming a linear correlation between gas-to-dust ratio and metallicity, as those reported in Tacconi et al. 2018 andBoogaard et al. 2021).
This experiment would not have been possible without the gravitational amplification effect on this galaxy, as of order ∼ 1000 h on the VLA would have been required. In the near future, the ALMA Band-1 receiver (Huang et al. 2016) might be able to detect CO(2 → 1) in similar unlensed galaxies with observations on the order of a few 10 s of hours although limited to z 5.5. The detection of cold molecular gas in 'normal', unlensed starforming galaxies at higher redshifts is thus prohibitive with current facilities and requires a ten-fold improvement in sensitivity. This stresses the necessity of new generation facilities such as the next-generation VLA (ngVLA, Carilli et al. 2015) to directly pin down the molecular gas content of galaxies within the epoch of reionization and to calibrate alternative tracers.
Given that the gravitational amplification factor in this galaxy is similar to the expected increase in sensitivity of the ngVLA, this work can be seen as a preview of the unique expected capabilities of the ngVLA to study the cold ISM in galaxies within the first billion years of the Universe. At the same time, this detection solves previous concerns about the detectability of such emission lines due to the increasing temperature of the CMB, anticipating the success of future radio facilities, although we note that the CMB impact could be higher in more extended sources with lower gas kinetic temperatures.
We thank Mark Lacy, Drew Medlin, and Aaron Lawson for their support with data reduction and the organizers of the NRAO 8th VLA Data Reduction Workshop. We also thank Kevin Harrington, Ryota Ikeda, Ikki Mitsuhashi, and Eric Murphy for fruitful discussions during the preparation of this manuscript and Gergö Popping for facilitating the catalogs of simulated galaxies. Finally, we thank the anonymous reviewer for a constructive report and valuable suggestions.
The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.2.00128.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.  Tadaki et al. (2022). The lines have been scaled to qualitatively match the line peak of the CO(2 → 1) line in order to compare the line profile and line-width. A double-peak line shape is clearly visible in both the atomic and molecular lines.