Low star-formation activity and low gas content of quiescent galaxies at $z=$ 3.5-4.0 constrained with ALMA

The discovery in deep near-infrared surveys of a population of massive quiescent galaxies at $z>3$ has given rise to the question of how they came to be quenched so early in the history of the Universe. Measuring their molecular gas properties can distinguish between physical processes where they stop forming stars due to a lack of fuel versus those where star-formation efficiency is reduced and the gas is retained. We conducted Atacama Large Millimeter/sub-millimeter Array (ALMA) observations of four quiescent galaxies at $z=$ 3.5-4.0 found by the Fourstar Galaxy Evolution Survey (ZFOURGE) and a serendipitous optically dark galaxy at $z=3.71$. We aim to investigate the presence of dust-obscured star-formation and their gas content by observing the dust continuum emission at Band-7 and the atomic carbon [C I]($^3P_1$-$^3P_0$) line at 492.16 GHz. Among the four quiescent galaxies, only one source is detected in the dust continuum at $\lambda_{\rm obs} = 870 {\rm \mu m}$. The sub-mm observations confirm their passive nature, and all of them are located more than four times below the main sequence of star-forming galaxies at $z=3.7$. None of the targets are detected in [C I], constraining their gas mass fractions to be $<$ 20%. These gas mass fractions are more than three times lower than the scaling relation for star-forming galaxies at $z=3.7$. These results support scenarios where massive galaxies at $z=$ 3.5-4.0 quench by consuming/expelling all the gas rather than by reducing the efficiency of the conversion of their gas into stars.


INTRODUCTION
Galaxies show a clear bimodality in physical properties, with distinct star-forming and quiescent populations. This is seen in observations such as the colormagnitude diagram (e.g., Baldry et al. 2004), the restframe UVJ diagram (e.g., Williams et al. 2009), and the stellar mass-star formation rate (SFR) diagram (e.g., Renzini & Peng 2015). Galaxies are considered to transform, via a process known as "quenching", from young blue star-forming galaxies to older and redder galaxies with little or no star-formation. Revealing how such quenching occurs and changes across cosmic time is crucial to understand the physical mechanisms behind galaxy formation and evolution in the Universe.
Although various physical mechanisms have been proposed, it remains unclear what is the main driver of galaxy quenching (see Man & Belli 2018, and references therein). One of the key observables to distinguish quenching processes is molecular gas, because gas mass fraction (f gas = M gas /(M * + M gas ))s and gas depletion timescale (t dep = M gas /SFR) are expected to change depending on the quenching mechanisms. In massive galaxies it is thought that a principle process is feedback by an active galactic nucleus (AGN) expelling gas and thus stopping star-formation (e.g., Di Matteo et al. 2005). Galaxies could also reduce the star-formation activity by virial shock heating (e.g., Kereš et al. 2005;Dekel et al. 2009) and halo gas stripping (i.e., starvation mechanisms; Larson et al. 1980;Peng et al. 2015). In both cases, the cold gas cannot fall into galaxies anymore and galaxies quench when they consume all their remaining gas. Gas-rich major mergers or violent disk instabilities can induce starbursts with short gas depletion timescales which quench galaxies by quickly consuming all the gas (e.g., Dekel & Burkert 2014). Morphological transformation, such as the formation of a central bulge, can dynamically stabilise gas against star formation (Martig et al. 2009) which would result in quiescent galaxies with high gas fractions and long gas depletion times.
Thus observations of molecular gas and star formation in quiescent galaxies can narrow the range of possibilities of quenching mechanisms. Studies of the gas properties of quiescent galaxies have been conducted up to z ∼ 3 directly using CO lines (Sargent et al. 2015;Spilker et al. 2018;Belli et al. 2021;Caliendo et al. 2021;Williams et al. 2021) or indirectly using dust continuum emission (Gobat et al. 2018;Magdis et al. 2021;Whitaker et al. 2021). Some of these studies revealed low gas mass fractions (f gas 0.1) and relatively short gas depletion timescales (t dep 1 Gyr) of quiescent galaxies at z ∼ 0.5-3.0, suggesting that these galaxies quench by strong feedback or rapid gas consumption due to active starformation (e.g., Sargent et al. 2015;Spilker et al. 2018;Whitaker et al. 2021). Gobat et al. (2018) obtained a gas mass fraction of ∼ 0.1 and a longer gas depletion timescale of 2-3 Gyr for quiescent galaxies at z ∼ 1.8 from a stacking analysis, suggesting that the quenching is caused by the reduced star-formation efficiency. The reported gas properties of quiescent galaxies at z 1 appear to show a large variety (Spilker et al. 2018;Belli et al. 2021;Williams et al. 2021), and there is no consensus on the typical gas properties of quiescent galaxies at high redshift yet.
In the last few years massive quiescent galaxies at z > 3 have been confirmed using near-infrared (NIR) spectroscopic observations (Glazebrook et al. 2017;Schreiber et al. 2018b;Tanaka et al. 2019;Forrest et al. 2020;Valentino et al. 2020a;Kubo et al. 2021). Investigating the gas properties of quiescent galaxies at z > 3, where there is less cosmic time for quenching to happen, allows insight into the range of physical quenching mechanisms as the quenching mechanisms operate on different timescales. Observations at sub-millimeter (mm) wavelengths are necessary to investigate the molecular gas contents of quiescent galaxies at z > 3. Furthermore, because the current level of star-formation in these galaxies is primarily constrained by rest-frame optical emission lines and spectral energy distribution (SED) fitting with the optical-NIR data, observations at longer wavelengths are needed to directly confirm whether they have dust obscured star-formation activity (Simpson et al. 2017;Schreiber et al. 2018c;Santini et al. 2019).
CO molecular lines have been used to measure the gas masses of quiescent galaxies at high redshifts (Sargent et al. 2015;Spilker et al. 2018;Belli et al. 2021;Caliendo et al. 2021;Williams et al. 2021). At z > 3, mid/high-J CO lines (J upper ≥ 4) are only observable with Atacama Large Millimeter/sub-millimeter Array (ALMA). One of the major uncertainties when calculating molecular gas masses from mid/high-J CO lines is the CO gas excitation state, which has been investigated only for dust-obscured galaxies such as sub-mm galaxies (SMGs) especially at z > 3 (e.g., Riechers et al. 2020;Liu et al. 2021). Atomic carbon [C i]( 3 P 1 -3 P 0 ) line (ν rest = 492.16 GHz) has been proposed as a good tracer of molecular gas in galaxies with advantages over CO molecular lines at higher redshift (e.g., Bisbas et al. 2017). Previous observational studies have reported [C i] detections from normal starforming galaxies as well as SMGs at z > 1 (e.g., Walter et al. 2011;Bothwell et al. 2017;Popping et al. 2017;Valentino et al. 2020b). Given that no assumption on the CO gas excitation state is required, the [C i] line can be a good tracer of molecular gas for quiescent galaxies at high redshift and has the additional advantage that it can reduce the uncertainty on the molecular gas mass measurement.
In this paper, we report results obtained from ALMA observations of quiescent galaxies at z = 3.5-4.0, and one optically dark galaxy at z = 3.71 which is a companion of one of the quiescent galaxies and is in a transition phase to quiescence (Schreiber et al. 2018c(Schreiber et al. , 2021, hereafter, S18c and S21). We investigate their dust-obscured star-formation activity and molecular gas content by observing the dust continuum emission at λ obs = 870µm and the [C i] line, respectively. We aim to confirm the passive nature of the quiescent galaxies at sub-mm wavelengths and to understand how massive galaxies at z > 3.5 stop star-formation by deriving a constraint on their molecular gas properties. This paper is organized as follows. In Section 2, we briefly explain the target galaxies. In Section 3, we describe the ALMA observations and explain the data reduction and stacking analyses. We show our results on the star-formation activity and gas properties of the quiescent galaxies at z = 3.5-4.0 in Section 4. We summarize the main findings in this study in Section 5. Throughout this paper, we use the cosmological parameters: H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3 and Ω Λ = 0.7. We assume the Chabrier (2003) initial mass function (IMF).
Jekyll is a massive quiescent galaxy at z spec = 3.715 with strong Balmer absorption lines, which was initially spectroscopically confirmed in G17. ZF-UDS-8197 is spectroscopically confirmed as z spec = 3.547 with the [Oiii]λλ5007,4959 doublet (S18b). The [Oiii] emission line in the MOSFIRE spectrum of this source has a broad line width of σ v = 530 ± 53 km s −1 . S18b argued that such a broad [Oiii] emission line is likely to be emitted from shock-excited gas rather than from narrow-line regions of the AGN. ZF-COS-18842 and ZF-COS-19589 are classified as sources with uncertain z spec in S18b. ZF-COS-18842 is identified with a single emission line with 5.6σ, which is attributed to the [Oiii]λ5007 line. ZF-COS-19589 is assigned z spec = 3.715 with a low probability of p = 32% based on weak Balmer absorption features. Given that this redshift lies within ∆z < 0.01 of Jekyll and this source is located only 23 arcsec away from Jekyll on the sky, there is a strong possibility that this galaxy is physically associated with Jekyll. Note that the spectroscopic redshifts of ZF-COS-18842 and ZF-COS-19589 are consistent with their photometric redshifts (S18b). None of these four galaxies is detected with X-ray observations either by Chandra Marchesi et al. 2016) or XMM-Newton (XMM-COSMOS: Cappelluti et al. 2007;Hasinger et al. 2007;Brusa et al. 2010, Subaru-XMM Newton Deep Survey: Akiyama et al. 2015).
G17 and S18b conducted SED fitting analyses for the quiescent galaxies with multi-wavelength photometry and MOSFIRE spectra to obtain their global physical quantities and star-formation histories. Their basic global physical quantities obtained by S18b are summarized in Table 1. The star-formation histories inferred from the SED fitting indicate that the four quiescent galaxies have formed at z ∼ 5.1-6.7, experienced a main star-formation period for 0.1-0.7 Gyr and then quenched 0.3-0.5 Gyr prior to the epoch of observation (S18b).

2.2.
Hyde: an optically dark galaxy at z = 3.71 Jekyll is accompanied by a source with sub-mm emission at a distance of ∼ 0 .4 (Simpson et al. 2017).
Hyde is considered to be a galaxy transitioning to quiescence because of its low SFR as compared to typical main sequence star-forming galaxies with similar stellar masses at z ∼ 3.7 (S21 and Section 4.2). Because Hyde could be representative of an important population relevant to understanding the quenching of massive galaxies at z > 3.5, we investigate its gas properties together with the four quiescent galaxies in this study.

Band-7 observation to trace dust continuum
The Band-7 observations of three quiescent galaxies, namely, ZF-COS-19589, ZF-COS-18842, and ZF-UDS-8197, were conducted in October 2019 as a part of the Cycle 7 observing program. The integration time is 16-17 min for the three targets. The central frequency of four spectral windows was set to be 343.5 GHz (∼ 870µm). We calibrated the raw data using the Common Astronomy Software Application package (casa; McMullin et al. 2007). We ran the clean algorithm with Briggs weighting (robust parameter=0.5). When there are sources with ≥ 5σ significance in the constructed maps, we ran clean again by masking the sources. The final sensitivity of the continuum maps at 870µm is 0.04 mJy/beam for ZF-COS-19589 and ZF-COS-18842 and 0.03 mJy/beam for ZF-UDS-8197. The average beam size is 0 .40 × 0 .35.
The Band-7 data for Jekyll and Hyde was taken as a part of the Cycle 6 program. Four spectral windows were set so that the [Nii] emission line at ν rest = 1461.13 GHz from Hyde can be covered (S21). The integration time is 48 min. We created the Band-7 continuum map for Jekyll and Hyde using the three spectral windows free from the emission line. The observed frequency is 302.34 GHz (991.6µm). The achieved sensitivity is 0.02 mJy/beam and the beam size is 0 .24 × 0 .20.
The Band-7 continuum maps of all the targets are shown in Figure 1 together with the Hubble Space Telescope WFC3/F160W images from CANDELS (Grogin et al. 2011). We conducted the source detection and flux measurement for the continuum maps using the casa task imfit. We fixed the central positions to the coordinates from the ZFOURGE catalog, which were determined with optical-NIR images (Straatman et al. 2016). Our detection criterion is that the peak flux from imfit has ≥ 3σ significance as done in Suzuki et al. (2021). ZF-COS-19589 and Hyde are detected with > 5σ. The peak flux obtained with imfit is regarded as a total dust continuum flux in the following analyses. The other three quiescent galaxies are not detected with > 3σ, and we assigned 3σ upper limits. The obtained continuum fluxes and upper limits are summarized in Table 1.
The offset between the coordinate from the ZFOURGE catalog and that of the 870µm continuum flux peak for ZF-COS-19589 is 0 .08. This is similar as a typical positional offset between ALMA sources and their Ks-band counterparts found in the COSMOS field reported in S18c. The observed positional offsets are considered to come from the positional accuracy of ALMA and VISTA/Ks-band images determined by the S/N of a source and the size of beam or point spread function (S18c). We can clearly infer that the detected continuum emission is associated to ZF-COS-19589, unlike the case of Jekyll and Hyde ( Figure 1) as discussed in S18c.

Band-3 observation to detect the [CI] line
In the Cycle 7 project, we observed the four quiescent galaxies with Band-3 to detect the [C i]( 3 P 1 − 3 P 0 ) line (ν rest = 492.16 GHz). The Band-3 observations of Jekyll and ZF-COS-19589 were conducted in October 2019. The observations of ZF-COS-18842 were conducted in March 2020 and March 2021. ZF-UDS-8197 was observed in March 2020, April, and June 2021. We set the frequencies of the spectral windows so that we cover the [C i] line from the targets. The integration time is 29, 80, 187, and 234 min for Jekyll, ZF-COS-19589, ZF-COS-18842, and ZF-UDS-8197, respectively. The minimum spectral resolution is set to be 7.81 MHz. We calibrated the raw data using casa and created the data cubes. The RMS level over 400 km s −1 veloc- Table 1. Summary of the observed and physical quantities of the four quiescent galaxies at z = 3.5-4.0 and Hyde (S21). All the galaxies are observed with ALMA Band-7 and Band-3.  The Band-3 observations of Hyde were conducted in September 2019 as part of the Cycle 6 program. The integration time is 198 min, and the minimum spectral resolution is set to be 7.81 MHz. The spectral data cube was created with casa/tclean after subtracting the continuum emission. The RMS level measured over 800 km s −1 width is 0.030 mJy beam −1 . We use the velocity width of 800 km s −1 for Hyde, which comes from the [Cii] line width measured by S18c. The beam size of the data cube is 0 .42 × 0 .27.
For each target, we created the data cubes with velocity widths of 100, 200, 400, and 800 km s −1 . We then ran the imfit task around the frequency corresponding to the redshifted [C i] line by fixing the central positions as done for the Band-7 data. We found no line feature with > 3σ in the four data cubes with different velocity widths for any of the sources. Figure 2 shows the onedimentional (1D) spectra extracted at the optical-NIR positions. We assign the 3σ upper limits on the [C i] line flux based on the RMS measured over 400 km s −1 velocity width for the four quiescent galaxies and 800 km s −1 velocity width for Hyde. The 3σ upper limits on the velocity-integrated [C i] line fluxes are summarized in Table 1.

Stacking of the dust continuum maps
We conducted a stacking analysis for the Band-7 continuum maps of the two non-detected quiescent galaxies, namely, ZF-COS-18842 and ZF-UDS-8197, with similar stellar masses and RMS levels. We stacked the two continuum maps by weighting each map according to its RMS level. The central position of each map is fixed at the coordinate from the ZFOURGE catalog. The RMS level of the stacked continuum map is 0.026 mJy. We conducted the source detection with imfit at the center of the stacked map as done for the individual galaxies (Section 3.1), and found no emission with ≥ 3σ. The dust continuum flux upper limit for the stacked sample is < 0.078 mJy (3σ).

Stacking of the [CI] data cubes
We conducted a stacking analysis for the [C i] line of the five galaxies including Hyde. First of all, we created the uv-tapered maps of the sources except for ZF-COS-18842 to create the data cubes with a beam size of ∼ 1 .0 × 0 .9. We consider that all the five sources are at the z spec determined by S18b. We then stacked the data cubes of the five sources. We here used the data cubes with the velocity width of 100, 200, and 400 km s −1 (Spilker et al. 2018). When stacking, the individual spectra were shifted to the rest-frame frequency and weighted according to the RMS level of the (tapered) data cubes at the frequency of the [C i] line.
The source detection in the stacked data cubes was done as described in Section 3.2. There is no detection with ≥ 3σ in the three stacked data cubes with different velocity widths. We determined the upper limit on the [C i] line luminosity using the data cube with 400 km s −1 width. The [C i] line luminosity upper limit for the stacked sample is log(L [CI](1−0) /L ) < 6.41 (3σ).
Note that only Jekyll has the stellar mass of log(M * /M ) > 11 among our targets (Table 1). Given the stellar mass dependence of gas mass fraction (e.g., Tacconi et al. 2018 obtained above, and thus, the gas mass estimated with the stacked [C i] line luminosity (Section 4.3) may be inappropriate for the less massive galaxies. In order to check this, we stacked the Band-3 data cubes for the four galaxies except for Jekyll. We confirm that the [C i] line luminosity upper limit, and thus, the gas mass upper limit for the stacked sample (Section 4.3) do not largely change depending on whether Jekyll is included or not. This means that the stacking result for the five galaxies is not heavily weighted toward Jekyll and appropriate for the less massive galaxies.

870µm flux density
In the left panel of Figure 3, we compare the ∼870 µm continuum fluxes of our targets at z ∼ 3.7 with galaxies at z = 3-4 individually detected with Band-7 from the literature (Liu et al. 2019a;Dudzevičiūtė et al. 2020) and quiescent galaxies at z = 3-5 with Band-7 data from Santini et al. (2019). Santini et al. (2019) investigated the dust-obscured star-formation activity of their photometrically identified quiescent galaxies with ALMA archival data (Band-6 and 7). In Figure 3, we show those quiescent galaxies that were confirmed with high confidence (Santini et al. 2019).
We find that our four quiescent galaxies, irrespective of being detected with Band-7 or not, have ∼ 1 dex fainter ∼ 870µm continuum fluxes than the SMGs at z = 3-4. Hyde also has a fainter sub-mm flux than most of the SMGs with similar stellar masses. This indicates low dust-obscured star-formation activities in our targets as discussed in more details in the next section. By conducting a deeper observation for a spectroscopically confirmed sample than the previous work (Santini et al. 2019), we can give strong upper limits on the dust continuum emission and thus SFR IR of quiescent galaxies at z > 3.5.

Dust-obscured star-formation activity
We converted the ∼ 870µm continuum flux to a total IR luminosity (L IR ) using a dust SED library by Schreiber et al. (2018a, hereafter, S18a). The required parameters are the dust temperature (T dust ) and the mid-to-total IR color, IR8 ≡ L IR /L 8µm . The SED templates have a dust emissivity index of β 1.5 (S18a). We assume IR8 = 7.37, a typical value of star-forming galaxies at z > 2 (S18a). When we change the IR8 value by −0.5 (+0.5) dex, the change of L IR is only +0.06 (−0.02) dex. This means that this parameter has a much smaller impact on the L IR estimates than T dust as discussed below.
Dust temperature is one of the major uncertainties when estimating L IR . Because it is difficult to estimate T dust of our targets with the available data, we here assume two dust temperatures, namely, T dust = 40 K and 20 K. T dust = 40 K is a typical value of star-forming galaxies at z = 3.5-4.0 from the relation between redshift and T dust given by S18a. T dust = 20 K is motivated by recent observations of quiescent galaxies up to z ∼ 2 (Gobat et al. 2018;Magdis et al. 2021). These studies reported that quiescent galaxies tend to have lower dust temperatures than star-forming galaxies at the same epoch. ) and quiescent galaxies at z = 3-5 with Band-7 data (Santini et al. 2019). The four quiescent galaxies in this work are fainter in sub-mm by 1 dex than SMGs at z = 3-4. (Right) SFRIR as a function of stellar mass of our targets together with quiescent galaxies at z = 3-5 with Band-7 and/or 6 data from Santini et al. (2019). We show the SFRs estimated with both T dust = 40 K and 20 K for each source except for Hyde. The data points with the lower SFRs represent the results assuming T dust = 20 K. The black solid line shows the main sequence of star-forming galaxies at z = 3.7 from Speagle et al. (2014). The dash-dotted line shows the star-forming main sequence at z = 3.7 from Tomczak et al. (2016), which is established with the ZFOURGE galaxy catalog. The four quiescent galaxies including ZF-COS-19589 with the continuum detection have more than four times lower SFRs than the main sequence galaxies even when assuming the higher T dust . All the upper limits are 3σ upper limits.
By comparing the observed flux densities at Band-7 and the model SED normalized with L IR from S18a, we calculated L IR of the four quiescent galaxies. Using the two dust temperatures, we obtained two different estimates of L IR (both are included in Table 1). The 3σ upper limits on L IR obtained from the stacking analysis (Section 3.3) are log(L IR /L ) < 11.09 and 10.31 with T dust = 40 K and 20 K, respectively.
As a test, we also calculated L IR with a IR SED template (T dust ∼ 20 K) obtained by stacking analyses for massive quiescent galaxies at z < 2 in Magdis et al. (2021). The obtained L IR are broadly consistent with those obtained with the IR SED template from S18a assuming T dust = 20 K. The difference between the two measurements is less than 0.1 dex. This means that the applied IR SED templates do not significantly affect our results on L IR , and thus SFR IR , under a similar T dust .
We converted L IR to SFR using the equation from Kennicutt (1998). We divided the SFRs by a factor of 1.7 to take into account the difference between the Salpeter (1955) IMF and Chabrier (2003) IMF (Pozzetti et al. 2007). The right panel of Figure 3 shows the rela-tion between stellar mass (M * ) and SFR IR for the four quiescent galaxies and Hyde. The SFR IR of Hyde was estimated in S21 using the same IR SED library. The 3σ upper limit on SFR IR of the stacked sample becomes < 12.3 M yr −1 and < 2.1 M yr −1 with T dust = 40 K and 20 K, respectively.
In the right panel of Figure 3, we find that all the galaxies and the stacking result are located below the M * -SFR relation of star-forming galaxies (the so-called star-forming main sequence) at z = 3.7 (Speagle et al. 2014;Tomczak et al. 2016) by a factor of > 4. This result is valid irrespective of the assumed T dust . With T dust = 20 K, the SFRs of the four quiescent galaxies are more than 20 times lower than those of star-forming galaxies on the main sequence at z ∼ 3.7. Interestingly, ZF-COS-19589, which is detected at 870 µm, still appears to have a significantly lower SFR than the main sequence galaxies at the same redshift. Given that ZF-COS-19589 is associated to dust emission but has a weak star-formation activity like Hyde, ZF-COS-19589 may be a similar transitioning galaxy (S21). S18b measured the SFRs of the four quiescent galaxies with several methods, namely, the optical-NIR SED fitting, Hβ and [Oii] emission line flux estimated from the spectral fitting. The SFR IR (upper limits) obtained in this study are consistent with the SFRs estimated with the other methods within the uncertainties. Our ALMA observations confirm the low star-formation activity of the four quiescent galaxies at the sub-mm wavelengths in addition to at the optical-NIR wavelengths (S18b).
We note that the old stellar populations in quiescent galaxies can be an additional heating source of dust in the general interstellar medium (ISM). The old stellar populations in ISM can produce IR emission even in the absence of young and massive stars heating dust in birth clouds. Including the IR emission heated by old stars can lead to an overestimate of SFR IR (e.g., da Cunha et al. 2008). Indeed, S21 estimated that in Hyde, 60% of the total L IR is contributed to by old stars. However, the contributions from old stars, if any, do not significantly change our conclusions about the passive nature of the four quiescent galaxies.

Constraint on the gas mass fractions
We estimate the upper limits on the molecular gas masses of the four quiescent galaxies and Hyde using the upper limits on the [C i] line fluxes and the following equation (Papadopoulos & Greve 2004;Bothwell et al. 2017): (1) where D L is the luminosity distance in Mpc, the [C i]/H 2 abundance ratio X [Ci] = 3 × 10 −5 , the Einstein A coefficient A 10 = 7.93 × 10 −8 s −1 , and the excitation factor Q 10 = 0.6 (Bothwell et al. 2017). We then calculated the upper limits on the molecular gas mass fractions (f gas = M gas /(M gas + M * )) with the molecular gas mass upper limits estimated above and the stellar masses from S18b. The obtained 3σ upper limits on the molecular gas masses and gas mass fractions of the individual targets are summarized in Table 1. All of our targets have gas mass fractions of < 0.2. The molecular gas mass and gas mass fraction upper limit (3σ) of the stacked sample (Section 3.4) are log(M gas /M ) < 9.69 and f gas < 0.09.
We compare the gas mass upper limits from [C i] of ZF-COS-19589 and Hyde with the gas masses estimated from L IR . In the IR SED library of S18a, each SED template has a ratio of L IR and dust mass (M dust ). With this ratio, we can convert the obtained L IR to M dust , and . Relation between stellar mass and molecular gas mass fraction upper limit of our targets. Here we show the stacking result for all five galaxies. We also show quiescent galaxies at z = 1-3 from the literature (Sargent et al. 2015;Belli et al. 2021;Caliendo et al. 2021;Magdis et al. 2021;Williams et al. 2021) for comparison. The solid line shows the scaling relation between M * and fgas for galaxies on the starforming main sequence at z = 3.7 from Tacconi et al. (2018). The two dashed lines represent the scaling relations when galaxies are located 5× and 10× below the main sequence. Our individual galaxies and the stacking result at z = 3.5-4.0 have more than three times lower molecular gas mass fractions as compared to the scaling relation for the main sequence galaxies at z = 3.7.
then convert M dust to M gas using the gas-to-dust mass ratio as done by Magdis et al. (2021). Here we use the gas-to-dust mass ratio at solar metallicity (Magdis et al. 2021). Note that M gas obtained with this method is the atomic+molecular hydrogen gas mass. The gas mass of ZF-COS-19589 inferred from its L IR is log(M gas /M ) = 9.23 and 10.12 with T dust = 40 K and 20 K, respectively. These values are consistent with the gas mass upper limit from [C i] (Table 1). In the case of Hyde, the gas mass inferred from L IR (T dust = 31 ± 3 K) is log(M gas /M ) = 10.45 ± 0.3, which is 0.25 dex larger than the gas mass upper limit from [C i] (Table 1). This could be partly due to a non-negligible contribution from the dust emission heated by old stars on L IR of Hyde (S21). Given the uncertainty, however, the gas mass inferred from L IR is consistent with the gas mass upper limit from [C i] for Hyde. Figure 4, we show the scaling relation between M * and f gas for galaxies on the star-forming main sequence at z = 3.7 from Tacconi et al. (2018). Comparing with this scaling relation, we find that the gas mass fractions of the five massive galaxies in this work are more than three times lower than main sequence galaxies at the same epoch with similar stellar masses (at 3σ). The stacking result is located 6× below the scaling relation for main sequence galaxies. These results suggest that massive galaxies at z = 3.5-4.0 consume or expel most of the gas before or during the quenching phase.

Constraint on the gas depletion timescales
We estimate the upper limits on the gas depletion timescale (M gas /SFR) of ZF-COS-19589 and Hyde with confirmed SFRs (Section 4.2). The gas depletion timescale of Hyde is estimated to be < 0.3 Gyr. This is shorter than the typical depletion timescale of main sequence galaxies (∼0.4-0.6 Gyr) expected from the scaling relation of Tacconi et al. (2018). Combining with the upper limit on the gas mass fraction, Hyde appears to be losing almost all its gas content on a short timescale. This galaxy might be nearing the end of a starburst phase after which it will become quiescent. As for ZF-COS-19589, the gas depletion timescale upper limit (3σ) becomes < 0.4 Gyr and < 2.4 Gyr when assuming T dust = 40 K and 20 K, respectively. In order to give a further constraint on the depletion timescale of ZF-COS-19589, we would need to constrain its T dust . S18b characterized the star-formation histories of the quiescent galaxies at z = 3-4 with the parameter z quench , the redshift when SFR drops down to 10% of the SFR in the main formation phase (< SFR > main ), and z form , the redshift when half of the total stellar mass was formed. Here, the main formation phase is determined as the contiguous time period surrounding the time of peak SFR where 68 % of the integrated SFR took place. < SFR > main is the mean SFR during this formation phase (S18b; S18c). The interval between z form and z quench can be used as a proxy of the quenching timescale of galaxies. Our four quiescent galaxies dropped down to 10 % of < SFR > main in 0.15-0.51 Gyr since their formation. The star-formation histories from S18b and the low gas mass fractions obtained in this study, therefore, suggest that the molecular gas exhaustion in our galaxies happened on a short timescale of the order of 100 Myr. This is consistent with the gas depletion timescale of Hyde. These results imply that massive galaxies at z > 3.5 would have a gas depletion timescale with the order of 100 Myr when they start quenching and keep such a short depletion timescale at least during the quenching phase.
It is not clear how long galaxies can keep a short gas depletion timescale after quenching. Williams et al. (2021) suggested the possibility that galaxy quenching happens with a drop in gas mass fraction due to high star-formation efficiency or feedback, which is followed by a period of low gas mass fraction and long gas depletion timescale. Indeed, Martig et al. (2013) showed that the star-formation efficiency of elliptical galaxies with f gas ∼ 1.3 % and 4.3% is ∼ 5 and 1.2 times lower than that of spiral galaxies with the same gas mass fraction, respectively, in their simulations. It is suggested that the dynamical stabilization of a gas disc by a central bulge becomes effective once the gas mass fraction of a galaxy becomes sufficiently low (Martig et al. 2009(Martig et al. , 2013. If our quiescent galaxies at z = 3.5-4.0 are already gaspoor such as f gas ∼ 1 %, they may have a gas depletion timescale with the order of 1 Gyr (Gobat et al. 2018;Magdis et al. 2021) even though the depletion timescale was shorter when their quenching happened. 4.5. Comparison with quiescent galaxies at z < 3.5 We show quiescent galaxies at z ∼ 1-3 from previous studies (Sargent et al. 2015;Belli et al. 2021;Caliendo et al. 2021;Magdis et al. 2021;Williams et al. 2021) in Figure 4 for comparison. Our results on the gas properties of quiescent galaxies at z > 3.5 are consistent with Sargent et al. (2015); Caliendo et al. (2021); Williams et al. (2021) and also Whitaker et al. (2021), showing that the low star-formation activity of massive quiescent galaxies is due to lack of fuel rather than reduced star-formation efficiency. Figure 5 shows the redshift evolution of the gas mass fraction of quiescent galaxies from z = 0 to z = 4 by combining this work with previous studies (Sargent et al. 2015;Spilker et al. 2018;Davis et al. 2019;Belli et al. 2021;Caliendo et al. 2021;Magdis et al. 2021;Williams et al. 2021). In addition to the quiescent galaxies at z = 1-3 already shown in Figure 4, we show quiescent galaxies at z ∼ 0.8 (Spilker et al. 2018) and z ∼ 0 (Davis et al. 2019). We also show the evolution of the gas mass fraction of star-forming galaxies on the main sequence with log(M * /M ) = 10.8 from Tacconi et al. (2018).
In order to investigate a possible evolutionary path of quiescent galaxies at z ∼ 3.7, we calculate the evolution of the molecular gas mass with a closed-box model, i.e., no inflow and outflow, assuming a constant gas depletion timescale (Gobat et al. 2018;Spilker et al. 2018;Williams et al. 2021). We consider a galaxy with log(M * /M ) = 10.79 and log(M gas /M ) = 10.15 at z = 3.715 as ZF-COS-19589. In the closed-box model, a galaxy is considered to consume the remaining gas only by star-formation according to a constant gas depletion timescale without further gas accretion and removal. We assume that ∼40% of the mass of stars formed is returned to ISM (for the Chabrier IMF; Madau & Dick- inson 2014). In Figure 5, we show the two model tracks assuming two different gas depletion timescales, namely, t dep = 0.4 Gyr and 2.4 Gyr. These values come from the 3σ upper limits of ZF-COS-19589 with T dust = 40 K and 20 K (Section 4.4). Note that because the gas mass fractions of our galaxies at z = 3.5-4.0 are all upper limits, the model track in Figure 5 should be regarded as an upper limit. The evolutionary track with t dep = 0.4 Gyr in Figure 5 appears to be consistent with the gas mass fraction upper limits of quiescent galaxies at z ∼ 1.5-2 from Sargent et al. (2015); Caliendo et al. (2021); Williams et al. (2021). When we just focus on the gas mass fraction values, the quiescent galaxies at z ∼ 3.7 can evolve into a population of quiescent galaxies with little gas content at z ∼ 1.5-2 if they keep a short gas depletion timescale and experience no further gas accretion. The model track with t dep = 2.4 Gyr shows a flatter evolution and appears to be consistent with the results from Magdis et al. (2021). However, quiescent galaxies at z ∼ 0.5-1.0 with f gas ∼ 0.1 are more gas-rich than the model track with t dep = 2.4 Gyr. The observed gas mass fractions of quiescent galaxies at z < 3 show a large scatter at a given redshift, and we cannot explain all the data points at z < 3 with this simple model when assuming a single gas depletion timescale since z ∼ 3.7. This may suggest that quiescent galaxies at z = 3.5-4.0 have various gas depletion timescales although there would be the contributions from galaxies quenched at later epoch, i.e., z < 3.5 (Magdis et al. 2021).
Note that here we just focus on the gas mass fraction values across cosmic time and do not consider the stellar mass evolution. Indeed, some of the quiescent galaxies at z = 1-2 in the literature are systematically more massive (∼ 0.5 dex) than the quiescent galaxies at z ∼ 3.7 (Figure 4). It is unlikely that the stellar mass of a quiescent galaxy is increased by more than 0.5 dex only with residual star-formation. A fair comparison between quiescent galaxies between at z ∼ 3.7 and z ∼ 1.5 might be difficult with the current samples. It would be necessary to increase the sample size of quiescent galaxies at z > 3 (and at lower redshifts as well), and to give a stronger constraint on their gas properties in order to further discuss the evolution of the gas contents in quiescent galaxies across cosmic time.

Comparison with Jekyll analogs in a semi-analytical model
There are some attempts to identify massive quiescent galaxies at z > 3.5 in semi-analytical models (e.g., Qin et al. 2017;Rong et al. 2017). Here we focus on Qin et al. (2017), using the meraxes semi-analytical model, because this model predicts a similar number density of quiescent galaxies at z = 3-4 as the observed one within a factor of two (S18b). Qin et al. (2017) found three analogues of Jekyll in meraxes and investigated the time evolution of the analogs. They showed that the Jekyll analogs experienced an intense star-formation event and black hole growth via galaxy mergers at z ∼ 5-6. After the mergers, gas cooling in ISM is significantly suppressed by energy radiated from the central black hole. Eventually, the Jekyll analogs consume the remaining cold gas with a timescale of 100-300 Myr and then quench. The stellar masses of the Jekyll analogs at z ∼ 3.71 are log(M * /M ) = 10.95-11.01. The gas mass fractions of the Jekyll analogs in the simulation are f gas ∼ 0.05 at z ∼ 3.71, which is consistent with the observed upper limits on the gas mass fractions of the four quiescent galaxies and Hyde. The observed gas properties of our quiescent galaxies and Hyde are likely to be qualitatively consistent with the formation history of the Jekyll analogs shown in Qin et al. (2017). This result may suggest the idea that AGN feedback plays an important role in the quenching of massive galaxies at high redshift.
We note that meraxes cannot fully reproduce the observed properties of massive quiescent galaxies at z = 3.5-4.0 (Qin et al. 2017;S18b). It is shown that the Jekyll analogs in meraxes have a longer duration of the star-formation phase and quench at later times than the observed quiescent galaxies. This indicates that further improvements on the models would be required. On the observational side, a stronger constraint on the gas mass and gas depletion timescale for individual quiescent galaxies would be necessary for more quantitative comparison with theoretical models.

CONCLUSION
We showed the results obtained from sub-mm observations with ALMA of four quiescent galaxies at z=3.5-4 and one optical-dark galaxy at z = 3.7 named Hyde (S18c; S21). With Band-7 and Band-3, we investigated the presence of the dust-obscured star-formation and the molecular gas traced by the [C i] line. We find that the four quiescent galaxies, including one detected with dust continuum, are located below the main sequence of starforming galaxies at z = 3.7 by a factor of > 4 (at 3σ). We confirm that the quiescent galaxies have weak or little dust-obscured star-formation. This study demonstrates that the previous multi-wavelength photometric analyses and NIR spectroscopic follow-up observations successfully identified true quiescent galaxies among the mass-selected galaxies at z =3-4 (Spitler et al. 2014;Schreiber et al. 2018b).
None of the five targets, including Hyde, have detectable [C i] line emission. The upper limit on their gas mass fractions is estimated to be < 0.2. Comparing with the scaling relation for galaxies on the star-forming main sequence at z = 3.7 (Tacconi et al. 2018), the five massive galaxies have more than three times lower gas mass fractions than star-forming galaxies with similar stellar masses. Hyde has an upper limit on the gas depletion timescale of < 0.3 Gyr, which is shorter than a typical value of the main sequence galaxies at the same epoch. Although we cannot give a constraint on the gas depletion timescales of three out of the four quiescent galaxies, the upper limits of their gas mass fractions obtained by this study and the star-formation histories inferred from the SED fitting (G17; S18b) support a scenario where massive galaxies at z = 3.5-4.0 would quench with an abrupt exhaustion of the molecular gas rather than due to a reduction in star-formation efficiency.
Given the large variation of molecular gas properties of quiescent galaxies reported at z < 3 (e.g., Belli et al. 2021;Williams et al. 2021), deeper observations of the gas contents for larger numbers of quiescent galaxies at z > 3 will be necessary to fully understand the quenching of massive galaxies in the early Universe.