An ALMA Spectroscopic Survey of the Brightest Submillimeter Galaxies in the SCUBA-2-COSMOS Field (AS2COSPEC): Physical Properties of z = 2–5 Ultra- and Hyperluminous Infrared Galaxies

We report the physical properties of the 18 brightest (S 870 μm = 12.4–19.2 mJy) and not strongly lensed 870 μm–selected dusty star-forming galaxies (DSFGs), also known as submillimeter galaxies (SMGs), in the COSMOS field. This sample is part of an ALMA band 3 spectroscopic survey (AS2COSPEC), and spectroscopic redshifts are measured in 17 of them at z = 2–5. We perform spectral energy distribution analyses and deduce a median total infrared luminosity of L IR = (1.3 ± 0.1) × 1013 L ⊙, infrared-based star formation rate (SFR) of SFRIR = 1390 ± 150 M ⊙ yr−1, stellar mass of M * = (1.4 ± 0.6) × 1011 M ⊙, dust mass of M dust = (3.7 ± 0.5) × 109 M ⊙, and molecular gas mass of M gas = (α CO/0.8)(1.2 ± 0.1) × 1011 M ⊙, suggesting that they are one of the most massive, ISM-enriched, and actively star-forming systems at z = 2–5. In addition, compared to less massive and less active galaxies at similar epochs, SMGs have comparable gas fractions; however, they have a much shorter depletion time, possibly caused by more active dynamical interactions. We determine a median dust emissivity index of β = 2.1 ± 0.1 for our sample, and by combining our results with those from other DSFG samples, we find no correlation of β with redshift or infrared luminosity, indicating similar dust grain compositions across cosmic time for infrared luminous galaxies. We also find that AS2COSPEC SMGs have one of the highest dust-to-stellar mass ratios, with a median of 0.02 ± 0.01, significantly higher than model predictions, possibly due to too-strong active galactic nucleus feedback implemented in the model. Finally, our complete and uniform survey enables us to put constraints on the most massive end of the dust and molecular gas mass functions.


Introduction
The population of submillimeter galaxies (SMGs; Smail et al. 1997;Barger et al. 1998;Hughes et al. 1998;Eales et al. 1999) was first discovered more than two decades ago by the Submillimeter Common User Bolometer Array (SCUBA) mounted on the single-dish James Clerk Maxwell Telescope at 850 μm.Subsequent observations using higher-resolution interferometers such as the Atacama Large Millimeter/ submillimeter Array (ALMA), the NOrthern Extended Millimeter Array, and the Submillimeter Array, have revealed the nature of these dusty star-forming galaxies (DSFGs).Studies have found that SMGs are dust-rich (Swinbank et al. 2013;da Cunha et al. 2015;Donevski et al. 2020;Dudzevičiūtė et al. 2020;Pantoni et al. 2021; M dust  10 8 M e ) and gas-rich (Bothwell et al. 2013;Birkin et al. 2021; M gas  10 10 M e ) galaxies with active star formation (Gullberg et al. 2019;Frias Castillo et al. 2023; SFR ∼ 10-1000 M e yr −1 ) and high dust attenuation (Dudzevičiūtė et al. 2020;Shim et al. 2022; A V  3) that span a wide redshift range (z ∼ 1-6) where the number density peaks at cosmic noon (z ∼ 2-3; Chapman et al. 2005;Danielson et al. 2017;Chen et al. 2022b).SMGs typically exhibit infrared luminosities, integrated from 8 to 1000 μm in the rest frame, higher than 10 12 L e , which result from the reprocessing of radiation by dust that absorbs the UV/optical light and reemits energy in the infrared, thereby usually classified as ultra-or hyperluminous infrared galaxies (ULIRGs18 or HyLIRGs 19 ).
The similarity between the redshift distribution of SMGs (Chapman et al. 2005;Wardlow et al. 2011;Simpson et al. 2014;Chen et al. 2016;Dudzevičiūtė et al. 2020) and that of the cosmic star formation density (Madau & Dickinson 2014;Zavala et al. 2021) suggests that SMGs are located at the key epochs of cosmic stellar mass growth.Indeed, several studies reveal that SMGs could account for ∼20%-60% of the cosmic star formation density at z  1 (Barger et al. 2012;Swinbank et al. 2013;Cowie et al. 2017;Dudzevičiūtė et al. 2020).However, it is still an open question as to what triggers their extreme star formation.The high SFRs may be induced by major mergers of gas-rich galaxies, similar to ULIRGs in the local Universe (Sanders et al. 1988;Narayanan et al. 2009;Chen et al. 2015;Perry et al. 2023).However, ample gas accretion from the intergalactic medium may trigger disk instability (Narayanan et al. 2015;Tacconi et al. 2020).One way to address this issue is to measure scaling relations with statistical samples of SMGs for gas depletion time and gas fractions against redshift and offset from the star formation main sequence (MS) and compare the results to those based on less active and less massive star-forming galaxies (Genzel et al. 2015;Tacconi et al. 2020), where such correlations have been explained through Toomre stability criteria (Toomre 1964) by Tacconi et al. (2020).
In addition to the star formation triggering mechanisms, understanding the role of dust is crucial for comprehending galaxy formation, particularly since the properties of dust-rich SMGs by observations are poorly reproduced by state-of-theart simulations (Popping et al. 2017;Hou et al. 2019;Li et al. 2019;McAlpine et al. 2019).A commonly used method for studying the dust properties is fitting the far-infrared (FIR) spectral energy distribution (SED) with a modified blackbody (MBB) model, via which characteristic dust temperature, dust mass, and dust emissivity index can be estimated.For studies that lack sufficient photometry, β, which is degenerate with temperature, is commonly fixed.While β is typically assumed to be between 1.5 and 2.0 (Scoville et al. 2017;Kaasinen et al. 2019;Dudzevičiūtė et al. 2020), recent studies suggest that β is consistent with or possibly greater than 2 (Magnelli et al. 2012;Casey et al. 2021;da Cunha et al. 2021;Cooper et al. 2022;Bendo et al. 2023;McKay et al. 2023), which is similar to the value in the local Universe (Smith et al. 2013) and suggests a dust assembly scenario by dust grain growth in the interstellar medium (ISM).A detailed investigation of β is useful not only in constraining dust grain growth models (Hirashita et al. 2014) but also in analyses of even higher-redshift galaxies that usually have sparse sampling of millimeter photometry.
Finally, an effective constraint to testing theoretical models is the mass functions, which represent the number density of objects within a given mass range and a redshift bin.This comparison can help us understand how different mass budgets assemble over cosmic time.Recent interferometric blind spectroscopic surveys conducted by ALMA and the Very Large Telescope (VLA), such as the ASPECS (Decarli et al. 2019) and the VLA COLDz survey (Riechers et al. 2019), have presented the evolution of the CO luminosity function, which is analogous to the gas mass function when assuming a constant α CO , from z ∼ 6 to z = 0.While effective, however, due to a small field of view for interferometers, the constraints at the massive end of the gas mass functions are limited.A semianalytical model developed by Béthermin et al. (2022) predicts that the ¢ L CO function can reach to the brightest end at ∼2 × 10 11 L e , successfully describing the observational results up to this limit.Observational measurements of dust mass functions have been presented at epochs up to z ≈ 3 (Magnelli et al. 2019;Pozzi et al. 2019).However, state-of-the-art simulations (Popping et al. 2017;Hou et al. 2019;Li et al. 2019;McAlpine et al. 2019) have struggled to reproduce the observed high number density of high dust mass sources, highlighting the need for a different treatment of dust production in the model but also more measurements at the massive end to better constrain the dust production of early massive galaxies.
Recently, a catalog of approximately 1000 sources with 850 μm flux (S 850 ) ranging from 2 to 20 mJy was presented by Simpson et al. (2019) from the SCUBA-2 Cosmology Legacy Survey in COSMOS (S2COSMOS), covering an area of 1.6 deg 2 .Follow-up observations of the 182 brightest S2COSMOS sources (S 850 > 6.2 mJy) were carried out using the ALMA band 7 continuum, and the results were presented as the AS2COSMOS sample by Simpson et al. (2020), where they deblended the single-dish sources into 260 SMGs with 870 μm flux (S 870 ) ranging from 0.7 to 19.2 mJy with precise localization.To further exploit the complete selection of the brightest SMGs in the AS2COSMOS sample, we have recently carried out an ALMA band 3 spectroscopic survey of the 18 brightest AS2COSMOS SMGs.The first results regarding line detection of CO and [C I], redshift distributions, and their lensing nature have been presented in Chen et al. (2022b).In this study, we perform detailed analyses of the CO and [C I] emission lines, which are tracers of molecular gas, and conduct SED fittings of the X-ray-to-radio photometry to further estimate other physical properties of the AS2COSPEC sample.
The paper is structured as follows.In Section 2, we provide an overview of the AS2COSPEC sample and describe the ancillary multiwavelength photometry.In Section 3, we present analyses to estimate the physical properties of AS2COSPEC SMGs, including SED fittings using Code Investigating GALaxy Emission (X-CIGALE) and MBB models.In Section 4, we discuss various aspects of these properties, including active galactic nucleus (AGN) fraction, gas depletion time, dust emissivity index, and gas and dust mass functions.Finally, in Section 5, we summarize our results.We adopt a flat ΛCDM cosmological model with H 0 = 67.7 km s −1 Mpc −1 , Ω M = 0.31, and Ω Λ = 0.69 (Planck Collaboration et al. 2020) throughout this paper.

Sample and ALMA Band 3 Data
The AS2COSPEC sample presented in this work represents the 18 brightest (S 870 = 12.4-19.2mJy) SMGs drawn from the parent sample of 260 AS2COSMOS SMGs.These were originally detected by the ALMA band 7 follow-up continuum observations of the submillimeter sources uncovered by the SCUBA-2 survey at 850 μm in the COSMOS field (Simpson et al. 2019(Simpson et al. , 2020)).Data and data reduction of our ALMA band 3 spectral scan observations were presented in detail in Chen et al. (2022b).Here, we briefly describe the relevant points.
All 18 AS2COSPEC SMGs are detected in the band 3 (3 mm) continuum with high significance (signal-to-noise ratio, S/N  10) in the naturally weighted maps (Figure 1).Given the typical resolution of 4″, the 3 mm continuum is found to be not spatially resolved except for the three sources that have close companions, AS2COS0001.1, AS2COS0008.1, and AS2COS0028.1.For these three sources, their continuum images are produced with Briggs weighting, resulting in sufficient spatial resolutions for separating the main source and its companion.The 3 mm photometry is simply deduced from the IMFIT task in CASA, which yields consistent results to estimates based on fittings in the visibility space.The resulting flux densities are provided in Table 1.
Spectra are first extracted using an aperture equal to the synthesized beam, and the line intensities are measured using Gaussian profile modeling.Corrections to the total line intensities are then applied based on the curve-of-growth analyses.Line emission is detected in all but one AS2COSPEC SMG, and their line intensities, line widths, and corresponding luminosities are measured via line moments (Chen et al. 2022b).The results are provided in Table 1.
Spectroscopic redshifts are determined based on the measured emission lines.The majority (14 out of 18) of the sample SMGs have multiple emission lines from the original band 3 data, as well as other follow-up observations (Mitsuhashi et al. 2021;Frias Castillo et al. 2023;D. J. Taylor et al. 2023, in preparation).Three SMGs with single-line detection (AS2COS0044.1,AS2COS0066.1, and AS2COS0090.1) have their redshifts assigned to one that is closest to their photometric redshifts.Estimates based on sources with multiple line detections suggest that such an assignment of redshifts for single-line sources should be correct in ∼80% of the cases.For the remaining ∼20% of the cases, their photometric redshifts do not accurately capture the spectroscopic values (|z phot − z spec | > 1).The reason for this mismatch is currently unclear, possibly due to the fact that the modeling of the photometric redshifts does not fully capture the properties of this heavily obscured population.The only source that does not have any line detection (AS2COS0037.1) is assumed to be at the redshift gap (z = 1.74-2.05)due to the frequency coverage of the band 3 scan.Although our sample of SMGs is considerably brighter than the typical ones, strong gravitational lensing with magnification of greater than 2 can only be found in at most one source (AS2COS0002.1;Table 1).Since the evidence for lensing is weak, throughout this paper, we simply ignore lensing effects.We confirm, however, that including lensing corrections or not would not change our results significantly.

Multiwavelength Data
The COSMOS field (Scoville et al. 2007) has been a prime target field covering a wide wavelength range from X-ray to radio.We exploit these exquisite data sets and use them to measure and model the SED of our sample SMGs, which allows us to obtain estimates of their physical properties, such as stellar masses and SFRs.In Figure 1, we show postage stamps of all 18 sample SMGs using infrared imaging and overlay the 3 mm continuum from ALMA observations in contours.For our work, we mainly select and adopt the multiwavelength photometry collected by Simpson et al. (2020) for the parent AS2COSMOS sample.In the following, we briefly summarize their analyses and refer to Simpson et al. (2020) for more details.
The photometry measurements from Simpson et al. (2020) can be summarized as follows.First, the 260 AS2COSMOS SMGs were cross-matched to the COSMOS2015 catalog (Laigle et al. 2016) with a chosen matching radius of 0 85, which was shown via simulations to yield a false match rate of ∼7%.Second, Simpson et al. (2020) also included a few broad bands that had deeper images compared to those included in the COSMOS2015 catalog, so the photometry of those bands was updated.In particular, the YJHK s photometry was extracted based on the fourth data release of the UltraVISTA survey (McCracken et al. 2012), instead of the second release that was used in the COSMOS2015 catalog.In addition, the photometry of the grizY bands was also replaced by the catalog provided by the Hyper Suprime-Cam (HSC) Subaru Strategic Program (SSP; Aihara et al. 2018).Necessary aperture-to-total corrections were applied in order to obtain total flux densities for these bands.Duplicated observations were taken by different instruments for a few common filters, such as the r and i bands from HSC-SSP and Suprime-Cam, as well as the H and K s bands from VIRCAM and WIRCAM.In general, we adopt the deeper measurements if multiple choices are available.The 260 AS2COSMOS SMGs were also cross-matched to the Chandra COSMOS-Legacy Survey (Civano et al. 2016) within the 3σ positional uncertainty of the X-ray location.Finally, we only use the broadband measurements, since most of our SMGs are not detected in the narrow and medium bands, and including this photometry or not does not affect the results significantly.
In total, four of our 18 SMGs (AS2COS0014.1,AS2COS0 028.1, AS2COS0031.1, and AS2COS0066.1) can be matched to the catalog of the Chandra COSMOS-Legacy Survey (Civano et al. 2016), in which they employed the maximumlikelihood algorithm for source extraction.Two of our 18 SMGs (AS2COS0037.1 and AS2COS0065.1) have robust detections (S/N 3) in the u-band data from the Canada-France-Hawaii Telescope (CFHT).Approximately 60% of the sources are detected in the optical bands (B, V, g, r, i, z + ) from Subaru Suprime-Cam and HSC, and about 70% of the sources are detected in Y, J, H, and K s from the DR4 UltraVISTA survey.
For mid-infrared and FIR photometry, various deblending techniques were employed to address the issue of source confusion, which is caused by the modest spatial resolution of the corresponding images.For images taken by the Infrared Array Camera (IRAC; Fazio et al. 2004) on the Spitzer Space Telescope (hereafter Spitzer), IRACLEAN (Hsieh et al. 2012) was used for deblending, with an updated prior catalog by including the AS2COSMOS SMGs in the original prior made from a stacked zYHK s image.For the 24 μm image taken by the Multiband Imaging Photometer (MIPS; Rieke et al. 2004) on Spitzer, as well as the 100 and 160 μm images taken by the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) on the Herschel Space Telescope (hereafter Herschel), the "super-deblended" photometry was adopted, which was made by deblending those images with a K s + radio prior (Jin et al. 2018).Finally, for the images at 250, 350, and 500 μm taken by the SPIRE instrument on Herschel, deblending was performed following the procedures described by Swinbank et al. (2013), where they adopted prior combining sources detected at 24 μm, VLA 3 GHz, and ALMA 870 μm.
A good fraction of the sample was targeted by previous ALMA and/or Plateau de Bure Interferometer observations at 1.2-1.3mm (Brisbin et al. 2017;Smolcić et al. 2017), and the photometry is included in Simpson et al. (2020).Finally, for radio, the 3 GHz photometry was obtained by again crossmatching the AS2COSMOS SMGs to the catalog of the VLA-COSMOS 3 GHz Large Project (Smolcić et al. 2017).For the 3 GHz faint sources, the flux densities are deduced by taking the pixel values at the positions of the SMGs, which were then corrected to total flux densities based on morphological analyses of the stacked images.In addition to the 3 GHz photometry obtained by Simpson et al. (2020), we further cross-match AS2COSMOS SMGs to the VLA-COSMOS 1.4 GHz catalog presented by Schinnerer et al. (2010) using a matching radius of 1″.We identify 95 1.4 GHz counterparts among the 260 AS2COSMOS SMGs, and 10 of them are our AS2COSPEC SMGs.For the remaining eight SMGs in our sample that are not detected, we determine their flux density upper limit as 48 μJy, which is the detection limit of the VLA-COSMOS Deep Project.
Overall, the number of AS2COSPEC SMGs that is detected in each of these bands described above ranges from 1 to 18.In Table 2, we provide a summary of the number of detections for each band.

X-CIGALE SED Fitting
We employ the 2022.1 version of X-CIGALE20 (Boquien et al. 2019;Yang et al. 2020Yang et al. , 2022) ) to model the SED of our SMGs.Briefly, X-CIGALE is a package written in Python that aims at fitting multiwavelength photometry, from X-ray to radio, with the goal of estimating the physical properties of galaxies by the method of Bayesian analysis.Several emitting sources are taken into account in the fitting, including stellar emission, ionized nebular emission, AGN emission dust emission, and synchrotron emission.In Table 3, we present the selected models and the chosen range of parameters.We use the photometry described in Section 2 and summarized in Table 2 in the fitting.Any of the photometry that has an S/N lower than 3 is treated as an upper limit in X-CIGALE.
The basic principle that we follow is simply to assemble a set of parameter space enabled by X-CIGALE that produces the best fits with the lowest reduced χ 2 across the whole sample.As a check, we compare our results to those obtained from fittings using MAGPHYS, following the procedures presented in Dudzevičiūtė et al. (2020), where they established their methods based on testing against simulated galaxies.In the following, we provide justifications and highlight a few findings on a number of key selections.We present the testing results at the end of this section.

Star Formation History
We adopt a delayed exponential star formation history (SFH) model with an exponential burst.This assumption of two independent components of SFH is now commonly assumed, which is also adopted in MAGPHYS, and it has been argued that it better reproduces the stellar masses of SMGs (Michałowski et al. 2014).By experimenting with different ranges of the input parameters, we find that including a starburst (SB) as recent as 10 Myr is needed for more than half of the sample SMGs, which enables significantly better fitting results on FIR photometry.On average, compared to the fitting without this 10 Myr burst included, that having 10 Myr in the Note.The extraction and the estimation of the photometry in each band is described in Section 2.2.
SB model decreases the reduced χ 2 by 20% and boosts the L IR by 90%.

Simple Stellar Population
We employ the stellar population model from Bruzual & Charlot (2003) with a Chabrier initial mass function (IMF).The stellar metallicity is set to the solar value.We note that some recent studies have found evidence in favor of a non-Chabrier IMF (Zhang et al. 2018;Cai et al. 2020), which, however, has not been confirmed by other studies (Lagos et al. 2020;Lovell et al. 2021).Since this is still a topic of debate, we choose to adopt the Chabrier IMF, which has also been adopted by recent studies of DSFGs (Donevski et al. 2020;Dudzevičiūtė et al. 2020;Cardona-Torres et al. 2023).

Nebular Emission
We adopt nebular emission in X-CIGALE, which is based on the template from Inoue (2011).This includes hydrogen continuum emission and line emission from He II at 30.38 nm to [N II] at 205.4 μm.

Dust Attenuation
We adopt the two-component dust model (Charlot & Fall 2000, CF00) as the attenuation model in our SED fitting.This is motivated by the recent discoveries that the optical/ near-infrared and FIR emissions of SMGs are not necessarily colocated (Chen et al. 2015;Hodge et al. 2016;Smail et al. 2023).We set a stronger attenuation strength for the birth clouds than ISM, with slopes of −1.3 and −0.7 for the birth clouds and ISM, respectively.They are linked by the μ parameter, ranging from 0.001 to 0.44 of A V,ISM = 0.3-3.8.We also repeat the SED fitting by adopting the SB attenuation law from Calzetti et al. (2000) but do not find overall better results based on reduced χ 2 .

Dust Emission
We adopt the dust emission model from Draine et al. (2014).Briefly, this model separates the dust emission into two components: the first one is the dust emission from the stellar light absorption of the general stellar population with a single radiation field U ; min the second one is the dust emission from the star-forming region with the radiation field ranging from U min to U max following the power-law index α.There are four free parameters in this model: U min , α, q PAH , and γ. q PAH is the mass fraction of polycyclic aromatic hydrocarbon (PAH) of the total dust mass, which we set to q PAH = 0.47, 3.90, and 7.32, ranging from the minimum to maximum acceptable values in X-CIGALE.We set a wide tolerance of q PAH due to the fact that the observed 24 μm is located at the redshifted PAH emission region so that the 24 μm is sensitive to the PAH emission strength.Recent results from simulations also suggest a wide range of q PAH (Narayanan et al. 2023).

AGN
AGN are thought to exist in some SMGs (Alexander et al. 2005) and could play a significant role in galaxy evolution.We adopt the AGN template from Stalevski et al. (2012Stalevski et al. ( , 2016)), which models the central source emission in the torus with most of the dust in high-density clumps by considering the 3D radiative transfer.It also models the AGN emission from different lines of sight and derives the luminosity contribution from AGN to the total dust luminosity (fracAGN).We set a wide range of input fracAGN, 0.01, 0.05, 0.10, 0.30, and 0.50, to let the model select the best value.We find that all sources prefer the fracAGN solution between 0.01 and 0.10 (see Table 5), which means that AGN play a minor role of contributing to the infrared emission from SMGs.

Radio
The model for radio emission is comprised of a simple power-law synchrotron emission with four parameters.Two of them are the slopes of the power-law synchrotron emission for star formation (alpha_sf) and AGN (alpha_agn); one is the value of the FIR/radio correlation coefficient for star formation (qir_sf), which is and the last one is the radio-loudness parameter for AGN (R_agn).We set a range of qir_sf and R_agn for X-CIGALE to find the best solution, and alpha_sf and alpha_agn are assumed to be the same.Ten of the SMGs in our sample have both 3 and 1.4 GHz photometry measurements; thus, their spectral slopes are measured individually.For the remaining eight SMGs that only have 3 GHz photometry, we assume a slope of 0.7, following the original study of the 3 GHz catalog (Smolcić et al. 2017).Table 4 summarizes the adopted values for alpha_sf and alpha_agn of each source, and we adopt these values as the input parameters in the radio model.

Fitting Results and Comparisons with MAGPHYS
The left panel of Figure 2 shows an example SED fit of one of our SMGs.The SED fitting results of all the sources are presented in Appendix B. Each physical model is plotted with a different colored line, and the total SED fitting model is the black curve, which is the summation of all the solid lines.The relative residual ((Observation − Model)/Observation) of each photometry is shown in the lower plot.Table 5 presents the physical properties of each source, and the statistical values (mean and median) are also listed.The physical properties listed in Table 5 are the stellar mass (M * ), V-band attenuation (A V ), total infrared luminosity (L IR ), fraction of AGN luminosity to total dust luminosity (fracAGN), SFR, and dust mass (M dust ).
To confirm that our results are not sensitive to the chosen SED fitting codes, we also run the analyses using MAGPHYS, following the same procedures as those described in Dudzevičiūtė et al. (2020).The results of MAGPHYS are compared with those of X-CIGALE.To make a more general assessment, we also include fainter SMGs presented in Birkin et al. (2021), where spectroscopic redshifts and multiwavelength data are also available.We find that overall, the physical parameters are comparable between the two codes, except that the X-CIGALE estimated M dust is, on average, higher by a factor of 1.5 and 1.3 than that in MAGPHYS for our AS2COSPEC sample and the ALESS/AS2UDS sample from Birkin et al. (2021).The 100 Myr averaged SFRs from X-CIGALE are better matched to those from MAGPHYS.However, the 10 Myr averaged SFRs from X-CIGALE are better aligned to the L IR -based SFRs.
We take proper consideration of these differences when making comparisons with other works that employed MAG-PHYS (see Appendix A for more details regarding the comparisons).For discussion, we mainly focus on the results based on X-CIGALE, where the inclusion of the AGN component that allows modeling the X-ray photometry is useful in determining the physical properties of our more luminous SMGs.

MBB SED Fitting
In addition to fitting the X-ray-to-radio SED with multiple degree-of-freedom templates, we also employed a singletemperature MBB model to fit the FIR photometry.The MBB fitting method offers the advantage of simplicity and facilitates direct comparisons with results from other samples, in contrast to SED fitting with complex templates.
For the MBB fitting, we utilize photometry from PACS at 100 and 160 μm; SPIRE at 250, 350, and 500 μm; and ALMA at 870, 1250, and 3000 μm.It is worth noting that because of our focus on bright SMGs, more than 80% (15 out of a total of 18) of our sources have at least six FIR photometric data points.Together with spectroscopic redshift measurements, this allows us to obtain constraints that are much tighter than literature samples with only photometric redshifts.
The MBB model is described in the following.The flux density emitted by the dust with a characteristic temperature T dust and taking the radiative transfer mechanism into account is described by For the optically thin case, the dust emission can be approximated as where the factor (1 + z) is the k-correction term, ΔΩ is the solid angle of the source, B ν (T dust ) is the Planck function from dust, B ν (T CMB ) is the Planck function from the cosmic microwave background (CMB), and τ ν is the optical depth, which can be written as dust where Σ dust is the dust column density.In our fitting, we assume that the dust geometry is distributed in a homogeneous sphere (Inoue et al. 2020); therefore, the dust column density is  Σ dust = M dust /(4/3πR 2 ), where M dust is the dust mass, which is for the optically thin case.Different assumptions of the dust distribution geometry, such as a foreground circular obscurer (πR 2 ), a spherical shell (4πR 2 ), or a homogeneous sphere (3/4 × πR 2 ), would lead to a factor of about 4 difference in dust mass estimates (Inoue et al. 2020).We adopt a homogeneous sphere dust distribution because it yields the closest M dust values to those obtained from X-CIGALE.κ ν is the frequency-dependent dust opacity, which can be written as where β is the dust emissivity index and κ 0 is the emissivity of dust grains per unit mass at a reference frequency ν 0 .In this study, we adopt κ 0 = 0.4699 m 2 kg −1 at ν 0 = 353 GHz, which is consistent with the value reported by Draine & Li (2007), to ensure a fair comparison with the M dust estimated from X-CIGALE.
In our analyses, we considered both the general and optically thin cases for fitting, with M dust , β, and T dust as the free parameters.For the general case, we also estimated an additional parameter, λ thick , which is given by This quantity represents the wavelength at which the opacity t = l 1 thick (da Cunha et al. 2021) and provides a measure of the optical depth in galaxies, with higher values of t l thick indicating a more optically thick scenario.
To quantify the uncertainties and correlations among the parameters, we employ the Markov Chain Monte Carlo (MCMC) method to obtain probability distributions and parameter-space correlations through sampling.We present the best-fit results and corresponding uncertainties from both the optically thin and general model fitting in Table 6.We show an example result in the left panel of Figure 2, where the red and blue curves represent the best-fit MBB models for AS2COS0065.1 for the optically thin and general cases, respectively.Although the shape of the dust emission models appears similar, the best-fit parameters differ; the optically thin model exhibits higher values for M dust , while the general model has a higher characteristic dust temperature (T dust ).The fitting results for all sources are given in Appendix B.
For the optically thin case, the best-fit β ranges from approximately 1.5 to 3.2, with a median value and bootstrap uncertainty of 2.1 ± 0.1.The dust temperature spans a range from 21 to 43 K, with a median dust temperature of 29 ± 2 K.We also estimate the effective dust temperature (T eff ), obtained from Wien's displacement law for the best-fit MBB model, and find the median value that is the same as the characteristic temperature.The dust masses of our sample SMGs are all on the order of 10 9 M e , with a median value of (3.3 ± 0.6) × 10 9 M e .
On the other hand, the dust in the general model is, on average, about 8 K hotter compared to the optically thin model, as the optically thin model assumes that continuum emission escapes from the core, while the general model assumes that only the optically thin surface contributes wavelengths shorter than λ thick .However, the effective dust temperature is almost the same as that of the optically thin model, which is reasonable considering the similarity in the shape of the best-fit optically thin and general MBB models.The median value of the best-fit Notes.The uncertainty of the mean is determined by employing the error propagation methodology, which considers the uncertainty associated with each source, and the uncertainty of the median is estimated by the bootstrap method.a The uncertainty of each source is the standard deviation of total of 3500 models from the MCMC method.b Rest-frame wavelength at which the opacity t = In our analysis, we conducted an exercise to gradually exclude millimeter photometry from the fitting process.We discovered that the best-fit parameters are highly sensitive to these longer-wavelength data, which reside on the Rayleigh-Jeans side of the blackbody radiation spectrum.For instance, when the millimeter photometry is omitted, the model constraints are weakened, resulting in lower values of β and higher values of T dust .In contrast, including the millimeter photometry leads to more tightly constrained best-fit parameters, with a median β thin = 2.1 ± 0.1, in close agreement with recent studies (da Cunha et al. 2021).We also applied this exercise to the ALESS and AS2UDS SMGs of Birkin et al. (2021) and found consistent results with our AS2COSMOS SMGs.As illustrated in Figure 3, the fitting without millimeter photometry is unable to accurately constrain the model on the Rayleigh-Jeans tail.Therefore, we caution that comparisons of dust parameters obtained from MBB fitting with other literature should be approached with care, as the fitting outcomes are highly sensitive to the degrees of freedom, especially the longer-wavelength data.We note that similar findings have been reported in da Cunha et al. (2021) for their "well-sampled subset" and the"rest of the sample."

Molecular Gas Mass
The fuel of star formation is the cold molecular gas.Our data allow us to estimate molecular gas masses (M gas ) using two different tracers: CO luminosity and [C I] luminosity.The summary of our results is provided in Table 7.In the following subsections, we describe how these estimates are obtained.

CO-H 2 Conversion
It is common to estimate molecular gas masses via CO(1-0) luminosities, which can be obtained following gas,CO CO CO 1 0 where α CO is the CO-to-H 2 conversion factor in units of ( )  --M K km s pc 1 2 1 and the factor of 1.36 accounts for the helium mass.
Note.M gas,CO and M gas,[C I] represent the gas mass estimated from the CO and [C I] line luminosity, respectively.The first uncertainty represents the measurement uncertainty from line luminosity, while the second uncertainty is the systematic uncertainty from α CO and r J1 .The uncertainty of the M gas,[C I] median is left blank because there are few data to be calculated.For AS2COS0028.1, which has double line detections in our ALMA band 3 survey, we adopt the CO(3-2) measurement for M gas estimation, since the uncertainty of the CO(3-2) measurement is better than that of CO(4-3).For those having CO(1-0) observations, we adopt the ( ) ¢ - L CO 1 0 measurements from Frias Castillo et al. (2023) to estimate their M gas .
For α CO , we adopt α CO = 0.8 ± 0.6, which is obtained from the dynamical measurements with a 15% dark matter fraction assumption (Rivera et al. 2018).Overall, 17 of our 18 sample SMGs can have their molecular gas masses estimated via CO luminosity.The M gas,CO measurements are given in Table 7, and we find a median value of M gas,CO = (1.2 ± 0.1) × 10 11 M e .

[C I]-H 2 Conversion
In addition to the CO lines, we also detect two [C I](1-0) lines in two SMGs, AS2COS0011.1 and AS2COS0031.1, and this could also be used as a gas mass tracer.We estimate the gas mass using , and the factor of 1.36 accounts for the helium contribution.We adopt the α [C I] from Birkin et al. (2021), where they find α [CI] /α CO = 4.4 ± 0.6 from the linear fitting in L CO space.The [C I]-derived gas masses are (2.5 ± 0.2) × 10 11 and (1.5 ± 1.2) × 10 11 M e for AS2COS0011.1 and AS2COS0031.1, respectively, which are consistent with the CO-derived gas masses.

Discussion
Before discussing the physical properties of our 870 μm bright SMGs, we first select the physical parameters we adopt in the following subsections.For dust masses, we select those estimated by X-CIGALE.For gas masses, we adopt those estimated by CO emission lines.We do not include AS2COS0037.1 for gas mass analyses, since it does not have any CO line detection.For SFRs, we adopt the infrared luminosity inferred values (Wuyts et al. 2011) in order to make a fair comparison to a sample of MS galaxies presented by Tacconi et al. (2020), where the IR-based SFRs are also adopted.In addition, IR-based SFRs are much less modeldependent compared to the SED-inferred SFRs.

CO Line Properties
We first discuss the CO line properties of the sample.In Figure 4, we plot line luminosities converted to CO(1-0) using ratios described in Section 3.3.1 versus redshifts, line widths, and infrared luminosities.As a comparison, we also plot similar measurements of fainter SMGs reported in the literature (Bothwell et al. 2013;Birkin et al. 2021), as well as those of a local ULIRG sample (Chung et al. 2009).
First of all, since our AS2COSPEC SMGs are selected to be the brightest in their 870 μm flux densities, it is not surprising that they are among the most luminous in CO and infrared luminosites.Second, the median FWHM line width is 610 ± 50 km s −1 , slightly higher than the measurements of the fainter SMGs, which are reported to be 540 ± 40 and 500 ± 60 km s −1 by Birkin et al. (2021) and Bothwell et al. (2013), respectively.Finally, we observe that the CO luminosity of our sample does not strongly correlate with any of the plotted properties, where we fit the three properties versus CO line luminosity and find that the best-fit slopes are consistent with zero.This is understandable, since both theory and observations have shown that 870 μm flux density is a good tracer of dust mass (Dudzevičiūtė et al. 2020;Cochrane et al. 2023), so our sample selection of a pure 870 μm flux cut effectively means a selection of dust mass and therefore likely a selection of molecular gas mass and thus CO luminosity.The narrow dynamical range in flux density of our sample could also contribute to the flat trend.Indeed, by considering our SMGs and those in the literature together, we find positive correlations between the CO luminosity and the properties plotted in Figure 4.

AGN Properties
Since SMGs are proposed to be the progenitors of local massive quiescent galaxies and coevolve with AGN or quasars at cosmic noon, it is useful to quantify the level at which the SMG and AGN populations overlap.Here, we look at the relationship between AGN and the AS2COSPEC SMGs from the perspectives of SED fitting, X-ray detection, and radio excess.
First, while the fraction of AGN contribution to the total IR luminosity (fracAGN) is allowed up to 50% in our X-CIGALE SED fitting, all of the sources in our sample have a best-fit fracAGN of ∼0%-10%.Only two sources have Bayesian values higher than 5%, while all others are consistent with zero.Overall, the median Bayesian value of fracAGN = (1.0 ± 0.2)%, suggesting that the infrared luminosity is dominated by the dust emission heated by the stellar emission.
Second, as mentioned in Section 2.2, there are four counterparts found in the Chandra COSMOS-Legacy survey (Civano et al. 2016), which indicates that the X-ray-detected fraction is - + 22 18 13 %.To investigate whether the X-ray is mainly powered by AGN or star formation in these sources, we deduce the absorption-corrected rest-frame 0.5-8 keV luminosity (L 0.5-8 keV,corr ), following the method described in Marchesi et al. (2016), using our CO-based spectroscopic redshift and assuming a power-law energy distribution.The L 0.5-8 keV,corr of these four sources spans a range of (0.3-1.4) × 10 44 erg s −1 .As shown in Figure 5(a), the median X-ray-to-IR luminosity ratio L 0.5-8 keV,corr /L IR = (1.7 ± 0.8) × 10 −3 is consistent with the AGN-classified SMG in Alexander et al. (2005).These suggest that the X-ray emission of these sources is mainly powered by AGN rather than star formation.
By taking into account the difference in the X-ray sensitivity, the X-ray-detected fraction appears marginally higher in our sample compared to the less luminous SMG samples such as ALESS and AS2UDS (Wang et al. 2013;Stach et al. 2019).This correlation between IR luminosity and X-ray detection fraction is also reported in the literature using other FIRselected samples (Kartaltepe et al. 2010;Juneau et al. 2013).
Specifically, in Figure 5(b), we plot the fraction of X-ray counterpart identification rate against infrared luminosity, including our sample, fainter ALESS and AS2UDS SMG samples (Wang et al. 2013;Stach et al. 2019), and the 70 μmselected galaxies (Kartaltepe et al. 2010).We also plot the correlation curve presented by Kartaltepe et al. (2010), where they cross-match their sources to both the XMM and Chandra surveys.To make a fairer comparison with the SMG samples, we reconstruct the correlation curve of Kartaltepe et al. (2010) by first cutting the sample to only include z > 1 sources and then matching the X-ray depth to the Chandra COSMOS-Legacy survey (Civano et al. 2016).The matching of the X-ray depth is also applied to the ALESS and AS2UDS samples to assess the X-ray detection fractions.The reconstructed correlation curve of the 70 μm-selected sources is broadly consistent with the results based on SMGs.
Third, we identify the AGN by checking the excess of radio emission.Specifically, we compute the 1.4 GHz luminosity to the IR-based SFR ratio following the method described in Smolcić et al. (2017).By adopting the redshift-dependent threshold, which is the 3σ deviation to the VLA-COSMOS 3 GHz sources (Delvecchio et al. 2017), only two sources in our sample (AS2COS0001.1 and AS2COS0028.1) are identified as radio-excess sources, as shown in Figure 5(c).
Overall, from the points of view of multiwavelength SED fitting, X-ray counterpart, and excess of radio emission, we find that AGN play a minor role in each perspective.We summarize and present the best fit with a gray line.For all panels, the blue, purple, and red dots represent our AS2COSPEC SMGs with different CO transitions, the gray hexagons are the SMGs from Birkin et al. (2021) and Bothwell et al. (2013), the orange stars are the local ULIRGs from Chung et al. (2009), the gray lines are the best-fit functions including all SMG samples, and the gray bands denote the 68% confidence interval of the fitting.
Figure 5. (a) Absorption-corrected rest-frame X-ray luminosity vs. infrared luminosity.Our four sources that have X-ray counterparts are black dots, the ASAGAO millimeter galaxies (Ueda et al. 2018) are purple upward triangles, the UDF sources (Dunlop et al. 2016) are green downward triangles, the ALESS SMGs from Wang et al. (2013) are blue diamonds, the SCUBA SMGs (Alexander et al. 2005) are red stars (SB-classified sources are open and AGN-classified sources are filled), and local galaxies are yellow squares classified as AGN-dominated (labeled with "A") or star formation-dominated (labeled with "S") by Rigopoulou et al. (1999) and Tran et al. (2001).The dotted line shows the mean relation of the SB-classified SMGs from Alexander et al. (2005), the solid line represents the median relation of the AGN-classified SMGs from Alexander et al. (2005), and the dashed line demonstrates the median relation of quasars studied by Elvis et al. (1994).(b) X-ray identification rate as a function of infrared luminosity.Our AS2COSPEC sample is the black dot, the ALESS sample from Wang et al. (2013) is the blue diamond, and the AS2UDS sample from Stach et al. (2019) is the red hexagon.The error bars are calculated by Poisson statistics (Gehrels 1986), and the asymmetric error propagation is implemented (Gobat 2022).The purple and green trends are the X-ray-detected fraction of the 70 μm-selected galaxies obtained from Kartaltepe et al. (2010) with different redshift cuts.(c) Rest-frame 1.4 GHz luminosity to IR-based SFR as a function of redshift.Our AS2COSPEC SMGs are the black dots (sources with X-ray counterparts are boxed), the VLA-COSMOS 3 GHz sources (Smolcić et al. 2017) are the gray dots, and the 3σ deviation threshold of the 3 GHz sources is shown as the red curve (Delvecchio et al. 2017).
the identification of AGN in Table 8, where there is no overlap for all sources except AS2COS0028.1 and AS2COS0031.1.However, the low detection rate in the X-ray and the low excess rate in the radio do not preclude the existence of heavily obscured AGN, which may have been missed significantly by X-ray surveys as suggested by recent studies (e.g., Carroll et al. 2023).More mid-infrared data such as those from the James Webb Space Telescope (JWST) would help improve our understanding with regard to the relationship between AGN and SMGs.

Gas Fraction and Depletion Time
To quantify gas properties, we derive two parameters: the molecular gas fraction, defined as gas gas * which describes the size of the molecular gas reservoir normalized to the stellar mass, and the gas depletion timescale, defined as which is the timescale for a galaxy to convert its molecular gas to stars through star formation, assuming no gas replenishment.Overall, for the AS2COSPEC SMGs, we find a range of gas fraction from 0.1 to 3.3, with a median of μ gas = 0.9 ± 0.3, indicating that the molecular gas mass is, on average, comparable to the stellar mass.In contrast, the fainter SMGs in the same redshift range as ours from Birkin et al. (2021) have a median μ gas = 0.5 ± 0.1, which is about 2σ lower.This is a combined effect of fainter SMGs being, on average, ∼20% more massive in stellar mass and ∼20% less massive in molecular gas mass.On the other hand, the range for the depletion time for the AS2COSPEC SMGs spans from 20 to 270 Myr, with a median t depl of 90 ± 10 Myr, which is about 25% lower than the median value of 120 ± 30 Myr found for the fainter SMGs from Birkin et al. (2021) at the same redshift range.Our finding that both the increase of gas supply and the increase of star formation efficiency play a role in driving galaxies to move above the MS is consistent with recent population studies of star-forming galaxies at similar epochs (Tacconi et al. 2020;Scoville et al. 2023).
In Figure 6, we plot the evolution of the gas fraction and depletion time deduced by Tacconi et al. (2020), which is based on a sample of 2052 star-forming galaxies with d = MS ( ) log SFR SFR MS ranging from −2.6 to +2.2 at z = 0-5.5. Figure 6.Scaling relations of the gas fraction (μ gas = M gas /M * ) with redshift (a) and the gas depletion timescale (t depl = M gas /SFR) with redshift (b).In both plots, black dots are our AS2COSPEC SMGs, red triangles are the ALESS/AS2UDS SMGs (Birkin et al. 2021), and purple lines represent the predictions from the best-fit results of the MS galaxies (Tacconi et al. 2020) by adopting the median δMS and M * of our sample.The purple downward arrows demonstrate the systematic shift of the MS relations when adopting α CO = 0.8.
Note.AGN identification from multiwavelength SED fitting, X-ray counterpart, and excess of radio emission.
They claim that μ gas and t depl can be effectively described as a function of three parameters, z, δMS, and M * , and the correlations can be generally explained by the Toomre stability criteria (Toomre 1964).By comparing the functional curves with SMGs, we find that SMGs are broadly consistent with the relationships of the MS galaxies.Upon closer inspection, some data points appear to deviate from the scaling relationship of the MS galaxies.However, it is important to note that since the correlation functions are dependent on three parameters (z, δMS, and M * ), the correlation curves in Figure 6 are plotted by assuming median values of the SMG samples for δMS and M * .In addition, these scaling relations of the MS are obtained by using α CO = 4.36 for estimating the gas mass, which is systematically higher than what we calculated for the SMGs.It is thus nontrivial in this case to compare data points against the relationships using simple statistical tests.
To facilitate a more straightforward comparison, for each source, we use the relationships reported by Tacconi et al. (2020) and compute the expected μ gas and t depl given its z, δMS, and M * .We plot the results against the measured μ gas and t depl in Figure 7.In spite of the normalization caused by the α CO , we find a good agreement between the expected and measured μ gas , where the best-fit slopes for our AS2COSPEC sample, the ALESS/AS2UDS sample (Birkin et al. 2021), and the combination of both samples are consistent with unity.This infers that the correlation functions for μ gas proposed by Tacconi et al. (2020) are applicable to the SMGs.As for t depl , the best-fit slopes to the three subsets are all sublinear, which suggests that the gas depletion timescale for the SMG population is shorter compared to the one deduced based on the correlation proposed by Tacconi et al. (2020).These results suggest that while on average, SMGs share a common behavior with fainter and less massive star-forming galaxies in terms of their gas fractions, SMGs are more efficient in star formation, possibly due to being preferably located in dense environments where dynamical interactions are more common (Chen et al. 2022a;Rujopakarn et al. 2023;Smail et al. 2023).
We note that if the 100 Myr averaged SFR deduced from X-CIGALE is adopted, the median t depl of our AS2COSPEC SMGs and the fainter SMGs in the same redshift range as ours from Birkin et al. (2021) are 200 ± 60 and 270 ± 80 Myr, systematically higher than that derived using IR-based SFR.However, our AS2COSPEC sample is still 25% lower than the SMGs from Birkin et al. (2021).As for the t depl comparison to the MS star-forming galaxies from Tacconi et al. (2020), the best-fit slopes of our sample, the sample from Birkin et al. (2021), and the total SMG samples are 0.8 ± 0.5, 0.9 ± 0.2, and 0.9 ± 0.2.While this could suggest that SMGs may have a similar star formation efficiency as the less massive and less active star-forming galaxies, we stress that the adoption of different methods of deducing SFRs could make this conclusion misleading.We therefore make our conclusion based on the IR-based SFRs.

Dust Fraction
We begin the discussion about dust by deriving the dust fraction (μ dust = M dust /M * ), meaning the dust-to-stellar mass ratio.The dust production mechanisms include ejection from supernovae (SNe) and asymptotic giant branch stars, accretion from the ISM, and infalling from the intergalactic medium, while the mechanisms of dust destruction include SN shock, stellar feedback, and AGN feedback.The observed dust-tostellar mass ratio is therefore related to the balance between dust production and destruction and could potentially be used for constraining feedback physics in theoretical models.
In our sample, we find a median dust fraction of μ dust = (2.1 ± 1.0) × 10 −2 .This is in contrast to the median dust Figure 7.The μ gas and t depl comparisons between the observational measurements (our AS2COSPEC SMGs in black and the ALESS/AS2UDS SMGs from Birkin et al. 2021 in red) and predictions from the MS star-forming galaxies (Tacconi et al. 2020) are plotted in panels (a) and (b).For both plots, individual measurements are plotted as small transparent symbols, while the binned median values are shown as the large symbols.There are two and three bins for our sample and the literature sample, respectively.Uncertainties for the binned data are derived from the bootstrap method.The linear best fits to our sample (black), the literature sample (red), and the combination of both samples (green) are also plotted, with the upward arrows pointing to the systematic offset if adopting α CO = 4.36.fraction of μ dust = (6.7 ± 1.4) × 10 −3 obtained from the ALESS and AS2UDS SMGs of Birkin et al. (2021), where we apply the same analysis techniques to derive the dust properties.Similarly, Donevski et al. (2020) analyzed a sample of 300 Herschel-selected DSFGs up to z ≈ 5 and found a median dust fraction of μ dust = (6.6 ± 0.5) × 10 −3 using the same SED package, X-CIGALE, as our analysis.Furthermore, Dudzevičiūtė et al. (2020) determined a median dust fraction of μ dust = (8.2± 0.4) × 10 −3 for 707 AS2UDS SMGs, where the systematic difference between X-CIGALE and MAGPHYS in M dust is corrected.Comparing our sample to the SMG samples in the literature, our dust fraction is approximately four times higher.This difference is likely due to sample selection, as sources with brighter submillimeter flux tend to have higher dust masses.
Figure 8(a) displays μ dust plotted against redshift for our sample, the AS2UDS SMGs (Dudzevičiūtė et al. 2020), and the Herschel-selected DSFGs (Donevski et al. 2020).Following the method in Donevski et al. (2020) and adopting the MS correlation from Speagle et al. (2014), each sample is split into MS and SB subsamples, where the boundary between two subsets is defined as ΔMS(= SFR/SFR MS ) = 4. Overall, we find a good agreement in μ dust between SMG samples, but it appears higher than the Herschel-selected DSFGs in the SB subset.This could be understood as that the 850 μm selection tends to preferentially select higher dust mass sources.
To compare with model predictions, we also include two simulated data sets, which are also split into SB and MS subsamples in Figure 8(a).They are the cosmological galaxy formation and evolution model SIMBA (Davé et al. 2019) and the phenomenological model SIDES (Béthermin et al. 2022).While SIMBA provides dust mass, we calculate dust masses for SIDES sources using the model of Draine & Li (2007) and the provided infrared photometry.To ensure a fair comparison between observation and simulation, we only plot the simulated galaxies that have similar physical properties to our sources in Figure 8(a), namely, unlensed sources with log (M * /M e ) 10.44 and log (sSFR/yr −1 ) −9.32.The binned averages of the MS and SB samples from both simulations are shown as the solid and dashed lines in Figure 8.
We find that, for both MS and SB galaxies, the μ dust predictions from SIMBA are systematically lower than the observed values.The lack of dust-rich galaxies in SIMBA was also presented and discussed in Li et al. (2019).By comparing the observed gas-to-dust ratio and gas fraction with SIMBA predictions, Donevski et al. (2020) suggest that this discrepancy could be due to too long of a timescale for dust growth.The dust growth timescale is inversely proportional to gasphase metallicity, and most galaxies in SIMBA, especially at z  2, have subsolar metallicity (Figure 3 in Li et al. 2019), which suggests inefficiency in dust growth.The reason for SIMBA having too low of a gas-phase metallicity is unclear.On the other hand, the dust growth timescale is also inversely proportional to the cold gas surface density.According to a recent investigation conducted by Cochrane et al. (2023), it is proposed that the intensity of AGN feedback has a notable influence on the FIR sizes of central star-bursting regions, consequently affecting the cold gas surface density.Specifically, stronger AGN feedback is linked to larger sizes and lower cold gas densities.This could suggest that AGN feedback in SIMBA is too strong, which drives too low of a cold gas density, leading to too long of a dust growth timescale.Comparing the observed submillimeter sizes with those predicted by SIMBA would help confirm this hypothesis.
On the other hand, we find that SIDES predictions are in good agreement with Herschel-selected DSFGs, confirming the results reported by Donevski et al. (2020).However, they do not fully agree with SMGs.First, the predicted values are, on average, at the lower end of the measurements for the SB sources.Second, the trend in redshift for the MS sources goes in the opposite direction between the predicted and observed values; as a result, the deviation between the SIDES predictions and the observed values becomes more significant at z  4. To test if these discrepancies are simply caused by the 850 μm selection, in Figure 8(b), we plot SIDES model curves based on various 850 μm flux density cuts.We find that indeed, by imposing the 850 μm flux cuts, the predicted values are in better agreement with the measurements for the SB sources.However, for MS sources, SIDES predicts higher μ dust , up to ∼0.5 dex for the brightest SMGs.By comparing the distributions of dust and stellar mass, we conclude that the overestimating μ dust is mainly driven by about a factor of 2 lower stellar mass in SIDES SMGs, which could be caused by the large uncertainties of the most massive end of the stellar mass functions.Lastly, the decreasing trend in redshift predicted by SIDES can still be observed in Figure 8(b) despite adopting 850 μm flux density cuts.This could be due to the fact that SIDES has a limited predicting power at z  4, since it relies on observed correlations and stellar mass functions, which are mostly not well constrained at z  4.

Gas-to-dust Mass Ratio
Our data also allow us to estimate the gas-to-dust ratio, δ gdr , which has been observed to correlate with redshift (Saintonge et al. 2013;Péroux & Howk 2020) and metallicity (Leroy et al. 2011;De Vis et al. 2019).From the Spitzer Infrared Nearby Galaxy Survey, the δ gdr is roughly 110 (Draine et al. 2007).By adopting α CO = 1 and using the MBB FIR fitting method, Swinbank et al. (2013) obtain an average δ gdr = 90 ± 25 for the ALESS SMGs.If adopting the same dust distribution geometry of our analyses, δ gdr would become 68 ± 19 for the ALESS SMGs in Swinbank et al. (2013).
In our analyses, we estimated M dust using three different methods: X-CIGALE X-ray-to-radio SED fitting (Section 3.1), optically thin fitting, and general MBB fitting to the FIR photometry (Section 3.2).Based on these three dust mass estimations, we calculated the corresponding δ gdr , and the median values are 32 ± 3, 44 ± 10, and 56 ± 10, respectively.After rescaling the α CO to 0.8 and correcting for the ∼1.5 times difference in dust mass between X-CIGALE and MAGPHYS (Appendix A), the MAGPHYS-based median, δ gdr = 34 ± 4, of the ALESS and AS2UDS SMGs (Birkin et al. 2021) is consistent with our X-CIGALE-based median (δ gdr = 32 ± 3).Our MBB-based median values are also consistent with the estimates of the ALESS and AS2UDS SMGs when the restframe 870 μm luminosity is used as a tracer of M dust (Birkin et al. 2021).The relatively low gas-to-dust ratio could suggest supersolar gas-phase metallicity, which is in line with SMGs being heavily dust-enriched and consistent with some recent results (Birkin et al. 2023;Eales et al. 2023;Peng et al. 2023).

Dust Temperature
In Figure 9(a), we present the redshift dependence of L IR .We find that, compared to fainter SMGs, our AS2COSPEC SMGs are more luminous in L IR at given redshifts, as expected, and we do not find significant redshift dependency in L IR in our sample, which might be attributed to selection effects.
In Figure 9(b), we plot dust temperature (T dust ) against L IR , and we also include results from fainter SMGs from AS2UDS (Dudzevičiūtė et al. 2020) and Herschel-selected infrared luminous galaxies (L IR > 10 10 L e ) at 0.1 < z < 2 (Symeonidis et al. 2013).To ensure proper comparisons to the AS2UDS SMGs, we follow the methodology presented by Dudzevičiūtė et al. (2020) for the MBB fittings, meaning that they are not corrected for the CMB, and a fixed β = 1.8 and optically thin case are assumed.
Within a given L IR bin, the T dust of our bright SMGs is lower than that of the AS2UDS SMGs.This difference in temperatures likely arises from a selection bias in our sample, as our selection favors sources with higher dust masses due to the positive correlation between M dust and S 870 .Consequently, we preferentially select sources with colder T dust at a fixed L IR .
To demonstrate this selection effect more clearly, in Figure 9(b), we plot the selection function as the gray hatched area.This selection function is constructed using MBB models that have the same S 870 range as our sample while fixing the redshift at the median value (z = 3.3) and setting β to 1.8.
To discuss the correlation between dust temperature and redshift, we focus our analysis on a specific region of the parameter space shown in Figure 9(a), selecting sources within an L IR slice (L IR = (4-12) × 10 12 L e ) and a redshift bin (z = 2.0-4.0), which represents the region where most of our sources are located and exhibits less dependence on both L IR and redshift.This allows us to investigate the evolution of T dust more effectively.In our sample, we find no significant evolution of T dust , as illustrated in Figure 9(c).Note that the correlation observed for the AS2UDS sample shown in Figure 9 is likely influenced by the L IR -z dependence; therefore, it is not solely driven by the T dust -z correlation.

Dust Emissivity Index
The dust emissivity index (β) is reflected in the slope of the MBB model in the Rayleigh-Jeans tail.In the Rayleigh-Jeans regime, the slope is determined by the contribution from the Planck function (α pl ), given by n µ n a S pl , and the dust emissivity (β), given by S ν ∝ ν β .A steeper slope corresponds to a higher value of β.We derive a median β of 2.1 ± 0.1 for our sample, which is consistent with the value of 2.0 ± 0.1 obtained in Birkin et al. (2021).The slightly higher β in our sample may be influenced by the selection bias toward lower T dust sources.Despite this, β is comparable between our brightest SMGs and the literature SMGs within the uncertainties (Magnelli et al. 2012;Casey et al. 2021;da Cunha et al. 2021).
Figures 9(e) and (f) present the variations of β as a function of L IR and redshift for our sample, the literature SMGs (Birkin et al. 2021;Bendo et al. 2023;McKay et al. 2023), and the local Herschel-selected galaxies from Smith et al. (2013).We perform linear fits to the parameters and find no correlations with either L IR or redshift, suggesting that the dust grain properties do not significantly correlate with L IR or redshift in galaxies.We again emphasize the importance of millimeter photometry in constraining the MBB model, as presented in Figure 3.
Several dust grain properties, such as dust grain size and composition, can influence the value of β.Theoretical studies have demonstrated an inverse relationship between β and dust grain size for different types of dust compositions (Ysard et al. 2019).For instance, β is approximately 1 for millimeter-to centimeter-sized amorphous silicate (a-Sil) grains, while β is around 2 for 0.1-10 μm sized a-Sil grains.This trend aligns with observational findings, where protostellar disks exhibit β ≈ 1 for millimeter-sized dust (Guilloteau et al. 2011), and molecular clouds show β ≈ 2 for micron-sized dust (Sadavoy et al. 2013).In the context of galactic-scale observations, the emission in the Rayleigh-Jeans region mainly originates from submicron-to micron-sized dust grains, as illustrated in Jones et al. (2013).This suggests that the β value we found in our sample is less dependent on dust grain size.However, the composition of the dust grains plays a more significant role in determining β.Simulations conducted by Hirashita et al. (2014) indicate that a β value of 2 can result from emissions by either graphite or silicate grains.Experimental studies, such as Inoue et al. (2020), have demonstrated significant variations in β across different materials.The issue of β and composition does not have clear-cut conclusions, and there is a degeneracy between β and the exact dust grain composition.Therefore, in our findings of β ≈ 2 for SMGs at cosmic noon, we cannot definitively identify the dust grain composition.However, because there is no correlation between β and redshift shown in Figure 9(f), it is possible that the dust grain composition may not undergo strong evolution.
To gain a comprehensive understanding of dust grain properties, relying solely on β is inherently limited.While β provides some insights, it primarily informs us about dust emission in the Rayleigh-Jeans region and allows us to rule out certain dust materials based on their associated β values.To conduct a more comprehensive study of dust properties, it is crucial to incorporate the mid-infrared data.In this regard, the JWST emerges as a powerful tool that can significantly contribute to our understanding of this issue.The mid-infrared observations enabled by JWST can provide valuable information about smaller dust grains and PAHs (Spilker et al. 2023), facilitating further investigations into the intricacies of dust grain properties.

Mass Function
Since our sample is a result of a uniform S 870 -selected survey, we can derive mass functions based on our measurements and compare to the model predictions.

Dust Mass Function
We derived the dust mass function at a given redshift bin by where i represents the index of sources in a given redshift bin.We divided our sample into two redshift bins: a lower redshift bin with eight sources and a higher redshift bin with nine sources.The redshift bin boundaries are given as z = [2.26, 3.29, 4.78], with the centers of each bin located at z = 2.78 and 4.04.N i represents the number of galaxies within a certain redshift bin, ΔV corresponds to the comoving volume of the corresponding redshift bin, ΔM denotes the dust mass bin width for sources within the redshift bin, and C stands for the completeness.Table 9 presents the estimated dust mass function values for the lower and higher redshift bins.Since our sample is flux-limited, it is important to account for the incompleteness of our selection in mass in order to include all sources within the given dust mass bins.As illustrated in the M dust -S 870 plot provided in Dudzevičiūtė et al. (2020), while the correlation is strong, there is still an ∼0.4 dex dispersion in M dust for a given flux.Consequently, our fluxlimited selection would miss a portion of likely warm and high M dust sources.To estimate the completeness in mass, for each redshift bin, we first compute the number of sources that are supposed to be above the measured minimum dust mass of our AS2COSPEC SMGs, based on the M dust -S 870 correlation (Dudzevičiūtė et al. 2020), scaled by the systematic offset between X-CIGALE and MAGPHYS, and the AS2COSMOS number counts (Simpson et al. 2020).The completeness is then the ratio between the total number of our sample SMGs and the total number of sources above this dust mass limit.As a result, we determine the completeness and Poisson uncertainties for the lower and higher redshift bins to be - + 15 % 6 8 and - + 9 % 3 4 .We then apply this correction factor to adjust the gas mass function.
In Figure 10(a), we compare our estimates to literature data that cover similar redshift ranges.Our measurements align well with the IRAM GISMO 2 mm survey sources from Magnelli et al. (2019), who modeled the dust mass of galaxies observed by the IRAM GISMO instrument by fitting mid-infrared to millimeter photometry using the Draine & Li (2007) model.Furthermore, we include measurements of the Herschelselected galaxies (Pozzi et al. 2019) and the rest-frame ∼180 μm-selected SMGs from STUDIES and AS2UDS (Dudzevičiūtė et al. 2021), where the dust mass estimations are derived through MBB fitting to the FIR and millimeter photometry.For the former, the dust temperature T dust was  (Magnelli et al. 2019), the Herschel-selected galaxies (Pozzi et al. 2019), and the rest-frame ∼180 μm-selected SMGs (Dudzevičiūtė et al. 2021).The x-axis error bars demonstrate the M dust bin width, while the yaxis error bars are the combination of Poisson errors and the uncertainties from completeness.Predictions from simulations (Popping et al. 2017;Hou et al. 2019;Li et al. 2019;Béthermin et al. 2022)   This suggests that SMGs have a cold gas reservoir that can be scaled from typical star-forming galaxies but are much more efficient in forming stellar mass, possibly caused by dynamical interactions.5.For the AS2COSPEC SMGs, the median dust mass derived from X-CIGALE modeling is (3.7 ± 0.5) × 10 9 M e .Given their stellar masses, AS2COSPEC SMGs have one of the highest dust-to-stellar mass ratios, with a median of (2.1 ± 1.0) × 10 −2 , four times more than any other DSFG samples, suggesting that bright SMGs are undergoing a phase of rapid dust mass production.Physically motivated models underpredict the observed values, possibly due to too low of a gas-phase metalicity or too-strong AGN feedback.6. Due to 850 μm selection, the dust temperature of our sample is biased toward lower values at a given infrared luminosity, with a median optically thin T dust of 29 ± 2 K.The median dust emissivity index β = 2.1 ± 0.1 agrees with measurements from some previous studies of less luminous sources (Magnelli et al. 2012;Casey et al. 2021;da Cunha et al. 2021).By combining recent β measurements of DSFGs in the literature, we find a lack of correlation of β with redshift and infrared luminosity.This may suggest common dust grain compositions for infrared luminous galaxies across a large fraction of cosmic time.7. Finally, we divide our sources into two redshift bins at z = 2-5 and observe no significant evolution in the dust and molecular gas mass functions at the high-mass end.These findings align with measurements from other surveys targeting sources at a similar redshift range.However, physically motivated models tend to underpredict the number density of these massive sources, again possibly suggesting too-strong feedback implemented in these models.
Our analyses provide valuable insights into the fundamental relations and ISM properties of the 18 brightest AS2COSPEC SMGs, contributing to the understanding of the evolutionary studies related to gas, dust, and star formation.These findings promote our knowledge of the unique characteristics of SMGs and their role in galaxy evolution.Future analyses of additional AS2COSPEC sources will be beneficial in validating the scaling relations by expanding the sample size.

Appendix B Best-fit Results of X-CIGALE and MBB
We include all the best-fit X-CIGALE and MBB fitting results of our AS2COSPEC SMGs in Figures 12 and 13.

Figure 1 .
Figure 1.Thumbnail images of our sample SMGs sized 20″ × 20″, where the background color map uses IRAC 4.5 μm, IRAC 3.6 μm, and UVISTA K s bands as red, green, and blue, respectively.The contours show the ALMA band 3 continuum at the 5σ, 7σ, 9σ, and 11σ levels.The spectroscopic redshift (Chen et al. 2022b) and 870 μm flux density(Simpson et al. 2020) of each source are also labeled on the image.We note that some of the sources have very weak detection or even no detection in the near-infrared due to dust extinction, but they are robustly detected in the FIR.

Figure 2 .
Figure 2. Left: an example source of the best-fit SED fitting from X-CIGALE (top panel) and the relative residual of each photometry (bottom panel).Photometry with S/N 3 is plotted as the black dots, and that with S/N < 3 is plotted as the green triangles.The black curve is the total best-fit SED model from X-CIGALE, and other colored curves represent the emission from different mechanisms.The dark red and dark blue curves plotted in the FIR regime are, respectively, the best fit of the optically thin and general MBB models (Section 3.2).Right: the corner plot of the MBB fitting from the MCMC of the same source shows the correlation of each parameter in the fitting.Again, red and blue denote the optically thin and general cases, respectively.Each dot represents a set of parameters that construct a model solution, and there are total of 3500 models for each source.The contours are shown at 0.5σ, 1σ, 1.5σ, and 2σ, while the vertical and horizontal lines denote the median value among the total of 3500 sets.
in the general MBB model is 2.1 ± 0.1, consistent with that deduced from the optically thin model.The median dust mass from the general MBB model is about 30% ± 20% lower but still comparable to the optically thin results when accounting for uncertainties.Again, since there is no significant difference between the best-fit optically thin and general models, the infrared luminosity integrated from the general model is comparable to the value of the optically thin model.Finally, by adopting a typical physical size of R = 2 kpc(Hodge et al. 2016), the best-fit λ thick ranges from 64 to 159 μm, with a median value and bootstrap uncertainty of 107 ± 9 μm, consistent with previous estimates for other SMGs(Riechers et al. 2013;Simpson et al. 2017).

Figure 3 .
Figure3.Distribution of the best-fit β, T dust , M dust , and L IR (from left to right) for the combination of our AS2COSPEC SMGs, as well as the ALESS and AS2UDS sources inBirkin et al. (2021).The pink and blue histograms represent the distributions from millimeter photometry-included and -excluded optically thin MBB fitting.The estimation of L IR does not significantly differ between these two sets.However, the other three dust parameters are more tightly constrained if the millimeter data are included.

Figure 8 .
Figure 8.(a) Dust fraction (μ dust = M dust /M * ) as a function of redshift.Black dots represent our AS2COSPEC MS SMGs, and red dots are the AS2COSPEC SB SMGs.Distributions of the Herschel-selected DSFGs (Donevski et al. 2020) and the AS2UDS SMGs (Dudzevičiūtė et al. 2021) are plotted as the blue and yellow shaded areas, where the MS subsamples are the nonhatched regions and the SB ones are the hatched ones.To remove the dust mass dependency due to SED fitting codes, systematic correction between MAGPHYS and X-CIGALE is applied to the AS2UDS SMGs.Simulated predictions from SIMBA (Davé et al. 2019) and SIDES (Béthermin et al. 2022) are plotted as the purple and gray lines, where the solid lines represent the MS galaxies and the dashed lines are the SB ones.(b) Same as panel (a), showing the SIDES predictions with different S 850 selection.Due to the small sample size, we do not plot the evolution of sources brighter than 8 mJy.Curves with a brighter S 850 cut are consistent with our sources, showing that, compared to literature samples, our high μ dust measurements could be due to the sample selection.

Figure 9 .
Figure 9. (a) L IR as a function of redshift of our sample colored by their T dust and the AS2UDS SMGs (Dudzevičiūtė et al. 2020; orange squares).The red box encloses sources, which are plotted in panel (c), at a given L IR bin.(b) T dust vs. L IR of our sample (black dots), the AS2UDS SMGs (Dudzevičiūtė et al. 2020; small transparent orange squares), and lower-redshift Herschel-selected galaxies (Symeonidis et al. 2013; purple diamonds).Larger yellow squares are the median values of AS2UDS SMGs from Dudzevičiūtė et al. (2020) at several L IR bins.The gray hatched area represents the selection function of our AS2COSPEC SMGs.(c) T dust vs. redshift of sources within the red box in panel (a).No significant evolution is observed in our sample.(d) Median T dust values for different subsets as depicted in panels (b) and (c).For the plots shown in panels (a)-(d) that compare to Dudzevičiūtė et al. (2020), we obtain the T dust and L IR of our sample from the non-CMB optically thin MBB fitting with β fixed to 1.8 using 100-870 μm photometry, consistent with the methodology in Dudzevičiūtė et al. (2020).(e) β vs. L IR for our sample (black dots), the ALESS/AS2UDS SMGs (Birkin et al. 2021; red triangles), the GOODS-S SMGs (McKay et al. 2023; purple hexagons), and the binned data of the local Herschelselected galaxies (Smith et al. 2013; orange stars).A linear fit is applied, revealing no correlation with a slope of 0.06 ± 0.03.(f): β vs. redshift for the same samples shown in panel (e), and measurements from the lensed Herschel BEARS program (Bendo et al. 2023) are plotted as blue pentagons.Again, a linear fit is performed, indicating no evolution with a slope of 0.03 ± 0.02.(g) Median β of different subsets in panels (e) and (f).For panels (e)-(g), the β values are derived from optically thin MBB models fitted to photometry ranging from 100 to 3000 μm.All SMG samples are concatenated together and binned into six bins, represented by green hexagons.

Figure 10 .
Figure10.(a) Dust mass function derived from our AS2COSPEC SMGs, the IRAM GISMO 2 mm survey sources(Magnelli et al. 2019), the Herschel-selected galaxies(Pozzi et al. 2019), and the rest-frame ∼180 μm-selected SMGs(Dudzevičiūtė et al. 2021).The x-axis error bars demonstrate the M dust bin width, while the yaxis error bars are the combination of Poisson errors and the uncertainties from completeness.Predictions from simulations(Popping et al. 2017;Hou et al. 2019;Li et al. 2019;Béthermin et al. 2022) are plotted in different line styles with huge variances in the low-mass region.(b) Gas mass function derived from our AS2COSPEC SMGs, AS2UDS SMGs (Dudzevičiūtė et al. 2020), the VLA COLDz survey (Riechers et al. 2019), and the ASPECS LP (Decarli et al. 2019).The x-axis error bars demonstrate the M gas bin width, while the y-axis error bars are the combination of Poisson errors and the uncertainties from completeness.Prediction of the semianalytical model is obtained from Popping et al. (2019), and those of the SIDES model (Béthermin et al. 2022) are converted from their Figure10.(a) Dust mass function derived from our AS2COSPEC SMGs, the IRAM GISMO 2 mm survey sources(Magnelli et al. 2019), the Herschel-selected galaxies(Pozzi et al. 2019), and the rest-frame ∼180 μm-selected SMGs(Dudzevičiūtė et al. 2021).The x-axis error bars demonstrate the M dust bin width, while the yaxis error bars are the combination of Poisson errors and the uncertainties from completeness.Predictions from simulations(Popping et al. 2017;Hou et al. 2019;Li et al. 2019;Béthermin et al. 2022) are plotted in different line styles with huge variances in the low-mass region.(b) Gas mass function derived from our AS2COSPEC SMGs, AS2UDS SMGs (Dudzevičiūtė et al. 2020), the VLA COLDz survey (Riechers et al. 2019), and the ASPECS LP (Decarli et al. 2019).The x-axis error bars demonstrate the M gas bin width, while the y-axis error bars are the combination of Poisson errors and the uncertainties from completeness.Prediction of the semianalytical model is obtained from Popping et al. (2019), and those of the SIDES model (Béthermin et al. 2022) are converted from their ( ) ¢ - L CO 1 0 function.

Figure 11 .
Figure 11.Comparison of the SED fitting results from CIGALE and MAGPHYS: (a) the reduced χ 2 , (b) the stellar mass, (c) the dust mass, (d) the dust luminosity, (e) the SFR averaged in 100 Myr, and (f) the SFR averaged in 10 Myr.Our SMG sample is plotted with the black dots, and the red triangles demonstrate the typical SMGs in Birkin et al. (2021).Note that for panels (e) and (f), the SFRs of MAGPHYS's estimations are both the 100 Myr averaged values.The medians of the X-CIGALE value to the MAGPHYS value (median(C/M)) for our AS2COSPEC SMGs and sources from Birkin et al. (2021) are shown in the legend in each panel.

Figure 12 .
Figure12.Best-fit models of X-CIGALE and MBB for our AS2COSPEC SMGs.Symbols are the same as in Figure2.

Figure 13 .
Figure13.Best-fit models of MBB fitting for our AS2COSPEC SMGs.The red vertical line marks the observed-frame λ thick , and other symbols are the same as in Figure2.

Table 2
Photometry Used in the SED Fitting and the Detection Situation in Each Band

Table 3
The Selected Modules and Input Parameters in CIGALE Note . Parameters that are not listed in this table are the X-CIGALE default values.The alpha_sf and alpha_agn values in the radio module vary from source to source, which depends on the radio synchrotron slope of the sources.See Section 3.1.7andTable4formore details.
(Smolcić et al. 2017)uncertainties are SMGs having both the 3 and 1.4 GHz photometry, so we derive their slopes from these two photometries.Those showing single values without uncertainties are SMGs having only the 3 GHz photometry, so we assume the slope of 0.7(Smolcić et al. 2017)in the X-CIGALE SED fitting.

Table 5
Physical Parameters from X-CIGALENotes.The uncertainty of the mean is determined by employing the error propagation methodology, which considers the uncertainty associated with each source, and the uncertainty of the median is estimated by the bootstrap method.

Table 7
Gas Mass M gas,[C I] Note.The three values of the mass function shown in each redshift bin represent the lower bound, the value, and the upper bound.Both F Mgas and F Mdust are in units of 10 −6 Mpc −3 dex −1 .