An ALMA Spectroscopic Survey of the Brightest Submillimeter Galaxies in the SCUBA-2-COSMOS Field (AS2COSPEC): Survey Description and First Results

Our sample is drawn from AS2COSMOS, a parent sample of 260 SMGs identified by ALMA band 7 continuum follow-up observations of a highly complete (99.5%), flux-limited (S₈₅₀ > 6.2 mJy) sample of the brightest submillimeter sources uncovered by the single-dish S2COSMOS over an area of 1.6 deg² (Simpson et al. 2019, 2020). The ALMA band 3 observations presented in this paper targeted the brightest 18 AS2COSMOS SMGs that are at the flux range of S₈₇₀ = 12.4–19.3 mJy based on the ALMA measurements (Simpson et al. 2020). The main goals are to measure the redshift distributions and gravitational lensing magnifications and study the ISM properties via molecular lines such as CO and [C i]. Comparing to the previous large samples of ALMA-identified SMGs from follow-up observations of single-dish detected submillimeter sources, in Figure 1, we show that there are only two SMGs that are similarly bright as our sample in the AS2UDS sample (Stach et al. 2018), and there are none in the ALESS sample (Hodge et al. 2013). Note that about one-third of the sample was observed in the millimeter by NOEMA and other ALMA data sets, and their measured redshifts were reported by Jiménez-Andrade et al. (2020), Simpson et al. (2020), and Birkin et al. (2021), which are consistent with our results.

The data were taken in 2019 November and December under project number 2019.1.01600.S. The observations were split into 11 execution blocks, each of which amounts about 1.5 hr worth of data for both science and calibrations. The default spectral scan setup was adopted in the ALMA Observing Tool, meaning that each execution block contains data on all 18 SMGs with five spectral tunings continuously covering from 84.1 to 113.3 GHz (Figure 2), where each tuning with four spectral windows was carried out for exactly 26 s. That is, across the whole frequency range observed, each spectral channel received a total of about 4.8 minutes (26 s times 11 executions) of science data on each SMG, except in the range of 97.7–101.3 GHz, as well as overlapping edge channels, where the on-target exposure time roughly doubles. Each spectral window is tuned to a standard time division mode with 128 channels and a frequency width of 15.625 MHz, corresponding to a velocity width of about 40–55 km s⁻¹, depending on the observed frequency. At least 43 and up to 49 antennas were used for each execution block, with a baseline length ranging from 15 to 313 m, hence a nominal C43-2 configuration. Since the target SMGs are all located in the COSMOS field, the phase calibrator was always J1008+0029, and J1037–2934 was always used for the amplitude and bandpass calibrations. The weather conditions during these observations were standard band 3 weather, with precipitable water vapor ranging from 1.7 to 5.8 mm.

Zoom In Zoom Out Reset image size

We adopt the calibrations performed in the second level of quality assurance (QA2). Flagging and calibrations were done using casa pipeline version 5.6.1-8, which is also the version used for imaging. Additional manual flags were put in by the QA2 team to remove antennas or scans with poor performance, in particular, a line at 91.66 GHz leaking from the local oscillator water vapor radiometer, which is marked in Figure 2. We reviewed the weblog record and confirmed the quality of the calibration results.

To image the visibilities, the calibrated visibilities were first continuum-subtracted. To do that, the channels that contain line emissions need to be excluded for a proper first-order continuum fitting. The line channels to be masked were determined from the spectra extracted from the delivered pipeline-produced imaging cubes through typical single Gaussian fitting spanning ±3σ centered at the Gaussian peaks based on the fits. The spectra were extracted using circular apertures centered at the 870 μm continuum positions (the phase centers), and their radius were determined by a curve-of-growth method.

We then used tclean to inverse Fourier transform and clean the continuum-subtracted visibilities. The visibilities were transformed to image cubes of 270 × 270 pixels in the x- (R.A.) and y- (decl.) axes with a pixel size of 042 (∼6–8 pixels per synthesized beam) and 1870 channels in the z- (frequency) axis. We adopt natural baseline weighting across all spectral windows in order to maximize the signal-to-noise ratio (S/N) for line detection. To increase the efficiency of clean, we adopted the automasking approach (Kepley et al. 2020) in tclean, which allows casa to generate masks for each channel in an iterative process based on a few parameters. To do that, we set the usemask parameter to "auto-multithresh." For each channel, the first round of mask creation only includes pixels with an S/N of 4 and above (noisethreshold = 4), and then the mask is allowed to expand in order to include lower-S/N regions (lownoisethreshold = 2.5). The masked regions determined by automasking are then cleaned down to the 2σ level (nsigma = 2 in tclean). Examples of the cleaned channels and their associated masks are shown in Figure 3. Finally, given the same baseline weighting across the whole frequency range, the spatial resolution of each frequency channel differs slightly, up to ∼35% end-to-end. To allow more straightforward analyses and understanding of the data, we applied smoothing on the reduced cubes by using the CASA routine imsmooth, setting the kernel to the common resolution, which is normally the largest beam size of the presmoothed cube. After this step, the spatial resolution of all cubes is uniform, with a synthesized beam FWHM of 43 × 35 and a small P.A. range of 69°–72°. The final spectral sensitivity achieved is uniform across all sources to within a few percent, thanks to the observing strategy; however, it differs among channels mainly due to the tuning coverage. We show the pixel-to-pixel rms in Figure 2, and the overall median is 0.33 mJy beam⁻¹ per 15.6 MHz (∼40–55 km s⁻¹). The continuum of each source was also imaged with natural weighting, resulting in images with a synthesized beam FWHM of 36 × 30 and a small P.A. range of 77°–80°. The sensitivity achieved for the continuum is also uniform across different sources, 12–13 μJy beam⁻¹ at the representative frequency of 98.7 GHz.

Zoom In Zoom Out Reset image size

To systematically search for line emission and quantify its significance, we ran the publicly available code LineSeeker, which was first written and used to search for line emission for the ALMA frontier field survey (González-López et al. 2017) and later applied to the data taken for the ALMA Spectroscopic Survey in the Hubble Ultra Deep Field (Walter et al. 2016; González-López et al. 2019). LineSeeker utilizes a matched filtering technique that combines spectral channels based on Gaussian kernels with a range of widths. The widths are chosen to match lines detected in real observations. The line candidates are determined to be significantly detected if their S/Ns are larger than the most significantly detected negative signals in the cube, and we find that a significance threshold of S/N > 10 satisfies the aforementioned requirement in all cases. In this paper, we focus on the spectra of the known SMGs in the AS2COSMOS catalog (Simpson et al. 2020), and the results from the search of serendipitous detection will be presented in a future work. For the 18 primary SMGs, 17 have yielded significant line detection based on LineSeeker, and AS2COS0037.1 is the only source that has none.

The spectra used for later analyses are extracted based on the coordinates presented in Simpson et al. (2020), which are based on ALMA 870 μm continuum observations. All of our targets have their 3 mm continuum detected, and we confirm that via 2D elliptical Gaussian fitting using imfit. The coordinates of the 3 mm continuum best-fit peaks are consistent with those of the 870 μm continuum. We therefore adopt the 870 μm coordinates, given the significantly lower positional uncertainties. The results of imfit also suggest that the 3 mm continuum is not spatially resolved, which is expected, since no strongly lensed SMGs like our targets are typically found to have subarcsecond sizes in the rest-frame submillimeter (e.g., Ikarashi et al. 2015; Simpson et al. 2015b; Hodge et al. 2016, 2019; Chen et al. 2017; Fujimoto et al. 2018; Gullberg et al. 2019).

To determine the optimal aperture size for extracting the spectra and to derive the total line flux density, we employ curve-of-growth analyses at the frequency of the emission lines and using the continuum-subtracted cubes. For each source, spectra are extracted based on a set of apertures that are centered at the 870 μm continuum position, shaped as a synthesized beam, and scaled by a range of factors in beam size from 0.5 to 4.5 with steps of 0.25. To estimate the line flux density, we conduct χ² model fitting for each spectrum extracted, and two models, single and double Gaussian, are considered. Given the channel-to-channel variations in sensitivity, the fittings are weighted by the noise of each spectral channel, which is estimated by taking the standard deviation of a collection of 1000 flux density measurements at random positions at same frequency but with the same aperture as that used to extract the target spectra. The line flux density of each extracted spectrum is then computed by integrating the best-fit model curves.

Because of the two adopted Gaussian models, the above approach results in two curve-of-growth curves for each target SMG. The selected examples are shown in Figure 4, in which we can first observe that the integrated line flux densities are not sensitive to the different underlying models. Second, while some curves converge to a certain flux density level, as expected from a curve-of-growth analysis, some appear not to be converging. We suspect that this is due to larger-scale noise structures, which has also been shown to be the case in other ALMA data sets (e.g., Novak et al. 2020). Indeed, by taking the median of all of the curves (one from each target SMG) normalized at a given aperture size, one beam in this case, we find that this median curve nicely converges beyond around three times the beam size. In addition, this converging median curve is stable to within a few percent against the choice of aperture size to which all of the curves are normalized, suggesting that this median curve can be used as a common profile for correcting the line-integrated flux densities to total. Note that there is a hint of extended emission by comparing the median curve to the curve of growth of the synthesized beam, although the difference is not significant (Figure 4). The detailed size analyses will be presented in a forthcoming paper. However, to account for this possible extended emission, we adopt the median curve instead of the curve of the synthesized beam for the correction to total (the difference is ∼20% ± 10%). We adopt the aperture equivalent to one synthesized beam partly because all of the curves appear to be stable until this point, beyond which they start to diverge (Figure 4). It also yields higher-S/N spectra compared to those extracted with larger aperture sizes, which helps to determine a more appropriate model for the spectral profile (Section 3.2).

Zoom In Zoom Out Reset image size

Given the above justifications, the spectra that we use for analyses from this point on are therefore extracted using the aperture equivalent to the synthesized beam, and their estimated line flux densities are corrected to total based on the median of the curve-of-growth curves. For AS2COS0001.1, AS2COS0008.1, and AS2COS0028.1, where the nearby SMGs are also detected in a few cases and close enough to potentially affect the spectra, we reimage the cubes using the Briggs baseline weighting with a robust parameter of −1.5, ¹⁹ resulting into an ∼35% reduction of the beam size, allowing clear separations to the companions with at least one synthesized beam distance (∼3''). The spectra of these three SMGs are then extracted based on the smaller beam and again corrected to total using the corresponding median curve. ²⁰ Finally, for AS2COS0037.1, the only target SMG for which we do not obtain significant line detection, the spectrum is extracted using the synthesized beam, but no correction to total is applied. An overview of all 18 spectra is shown in Figure 5, where we also show the error-weighted stacked spectrum in S/N. No additional emission lines are detected in the stacked spectrum. We have also experimented with luminosity-weighted and luminosity-normalized stacking; however, the results remain unchanged.

Zoom In Zoom Out Reset image size

The first analysis we perform on the extracted spectra is to determine which model, single or double Gaussian, better describes the data. To make a quantitative assessment, we employ both the Bayesian information criterion (BIC) and a version of the Akaike information criterion (AIC) that is modified for the limited sample size (AICc). The analyses are documented in Appendix A, and the best-fit models are shown in Figure 6. In general, we find that the results based on BIC and AICc are consistent with each other in 80% of the cases, and they all point toward the fact that about 90% of the primary SMGs have lines that are better described by double Gaussian profiles. The implications of this result are addressed in the discussion section.

Zoom In Zoom Out Reset image size

Moment calculations are adopted to compute line properties using best-fit single and double Gaussian models, meaning

$\begin{eqnarray*}&&\begin{array}{l}\bar{I}=\displaystyle \int {I}_{v}{dv}\\ \bar{v}=\displaystyle \frac{\displaystyle \int {{vI}}_{v}{dv}}{\displaystyle \int {I}_{v}{dv}}\\ \bar{\sigma }=\displaystyle \frac{\displaystyle \int {\left(v-\bar{v}\right)}^{2}{I}_{v}{dv}}{\displaystyle \int {I}_{v}{dv}}\end{array}\end{eqnarray*}$

and the first moment ($\bar{v}$), i.e., intensity-weighted central frequencies, are used to determine redshifts. The zeroth and second moments will be used in a forthcoming paper to discuss the line properties in detail.

To determine redshifts, the target SMGs are grouped into three categories: SMGs with multiple line detection, those with only one line detection, and those that do not have emission line detection. Other works in the literature are referenced to aid the redshift determinations, in particular those from Jiménez-Andrade et al. (2020), Mitsuhashi et al. (2021), and Daddi et al. (2022). The category to which each SMG belongs is given in Table 1. In total, 10 (56%) SMGs have their redshifts determined by multiple emission lines, seven (39%) by a single emission line, and one (6%) by its photometric redshift.

For the SMGs with multiple line detections, the measured line frequencies are cross-matched to combinations of redshifted emission lines, considering species that typically emit strongly at the submillimeter, such as carbon monoxide (CO) and atomic carbon ([C i]). We find a unique combination for each SMG, and the results are listed in Table 1. The final adopted redshifts are determined in two ways. If the SMGs have multiple line detections in this reported ALMA data set, their redshifts are determined by taking a weighted average of the redshift solutions of each line. If the second line detection originated from other work in the literature, we adopt our measurements, since they tend to have higher precision.

For the SMGs that only have one line detection, we find the solutions that are best matched to the photometric redshifts reported in S. Ikarashi et al. (in preparation), where, as in Dudzevičiūtė et al. (2020), the SED code magphys is employed. We assume the line to be one of the CO transitions, which is justified by the fact that the frequency coverage of our data would have covered a CO line if the line was instead [C i]. Since in SMGs, [C i] is typically weaker than the nearby CO lines (J_up = 4, 5; Spilker et al. 2014; Birkin et al. 2021), if the line is [C i], then the nearby CO lines would have been detected. We note that, compared to other methods, magphys has the advantage of taking into account far-infrared and submillimeter photometry, allowing redshift estimates for all AS2COSMOS SMGs, which are often faint or undetected in OIR wave bands, which the known COSMOS catalogs such as COSMOS2015 (Laigle et al. 2016) and COSMOS2020 (Weaver et al. 2022) are based on. However, given the rich information provided in the COSMOS catalogs, a good number of matches can be found on our target SMGs; we nevertheless exploit them later for detailed assessments.

Lastly, for the only SMG that has no line detection, AS2COS0037.1, the photometric redshift reported by S. Ikarashi et al. (in preparation), as well as the latest COSMOS2020 catalog (1.93^+0.04 _−0.07 and ${1.98}_{-0.15}^{+0.03}$ based on the LePhare and EAZY codes, respectively; Weaver et al. 2022), suggest that this SMG could be located within the known redshift gap of the ALMA band 3 spectral scan, which is z = 1.74–2.05. Indeed, given the depth reached by our data (∼0.3 mJy beam⁻¹), assuming a single Gaussian model of a 3σ peak and a typical line width of 500 km s⁻¹, the line intensity limit is 0.5 Jy km s⁻¹. Assuming the missing line is CO(3–2), the closest transition based on the photometric redshift, the line intensity limit corresponds to a limit on the CO(1–0) luminosity of ${L}_{\mathrm{CO}(1-0)}^{{\prime} }=(2\mbox{--}4)\times {10}^{10}$ K km s⁻¹ pc², assuming the latest ${r}_{31}\equiv {L}_{\mathrm{CO}(3-2)}^{{\prime} }/{L}_{\mathrm{CO}(1-0)}^{{\prime} }=0.63$ (Birkin et al. 2021). Given the infrared luminosity of ∼10¹³ L_⊙ (C.-L. Liao et al. 2022, in preparation), this luminosity limit is well below the measured ${L}_{\mathrm{IR}}\mbox{--}{L}_{\mathrm{CO}(1-0)}^{{\prime} }$ correlation and would have allowed detection of the line in ≳90% of the literature reported detection on SMGs (Birkin et al. 2021). We therefore adopt 1.90 ± 0.05 for AS2COS0037.1, which is simply the mean of the redshift gap with an uncertainty that makes the probability of it being outside the redshift gap equivalent to 3σ.

In Table 1, we report the redshift solutions of our sample on individual lines, as well as the final adopted values for each SMG, based on the methods stated above. We adopt the redshift results based on the selected model for each SMG listed in Table A1, although we note that these results are not sensitive to the chosen single or double Gaussian models for the line fitting.

In Figure 7, we compare the photometric and adopted spectroscopic redshifts for all target SMGs except for AS2COS0037.1. We find that based on the 10 SMGs that have multiple line detection, three sources have their spectroscopic redshifts differ from photometric redshifts by more than 1, effectively meaning that if only a single line were detected, it would have been assigned to an incorrect CO transition based on photometric redshifts. This suggests a 70% success rate for our strategy of using photometric redshifts to aid in determinations of spectroscopic redshifts for SMGs with a single line detection. However, note that this success rate is estimated based on the higher-redshift subset, which tends to have less reliable photometric redshifts. Thus, the success rate for our sample SMGs that have a single line detection could be higher.

Zoom In Zoom Out Reset image size

Our success rate is similar to that of a recent blind survey of emission lines on fainter SMGs (Birkin et al. 2021), in which they found about an 80% success rate. Excluding the three outliers, we find that the median redshift difference is 〈Δz/(1 + z_spec)〉 = 0.04 ± 0.04 with a bootstrapped error and a scatter of ${\sigma }_{{\rm{\Delta }}z/(1+{z}_{\mathrm{spec}})}$ = 0.06, suggesting that the photometric redshifts estimated via SED fitting using magphys can be accurate, on average, with good precision. However, we also find that the scatter is much larger than what can be inferred from the quoted uncertainties of the photometric redshifts, indicating that the uncertainties of the photometric redshifts are underestimated.

Based on the assessments above, we expect that for the seven SMGs that have only one line detection, about one to two sources may have their real redshifts different from, and likely higher than, the adopted values (details given in Appendix B). We take these into account when comparing our results to model predictions.

Interferometric follow-up observations in the past decade have convincingly shown that submillimeter sources detected by single-dish observations tend to break up into multiple sources with a probability strongly correlated with the flux density of single-dish measurements, from ∼10% for ∼3 mJy SMGs to >50% for SMGs at ≳10 mJy (Wang et al. 2011; Barger et al. 2012; Hodge et al. 2013; Simpson et al. 2015a, 2020; Cowie et al. 2018; Hill et al. 2018; Stach et al. 2019). While the number of constituent SMGs can range from two to four, it has also been shown that the flux densities measured from single-dish observations are dominated by the brightest constituent SMG, with a 70% contribution, on average.

The relationship among these constituent SMGs is, however, unclear, in particular whether they are physically associated and undergoing dynamical interactions or unrelated SMGs that are simply a result of line-of-sight chance projections. The determination of this relationship has profound implications in our understanding of the triggering mechanisms of star formation for SMGs. Observationally, statistical approaches using photometric redshifts and number density analyses have suggested that ∼30% of these pairs are physically associated (e.g., Stach et al. 2018; Simpson et al. 2020), which is in agreement with spectroscopic measurements of small samples of paired SMGs (Hayward et al. 2018; Wardlow et al. 2018). Theoretically, models have instead predicted that the majority of these pairs are physically unassociated, especially for the bright submillimeter sources with a single-dish measured flux density of S₈₅₀ ≳ 10 mJy (Hayward et al. 2013a; Cowley et al. 2015; Muñoz Arancibia et al. 2015). To definitely settle this issue, spectroscopic follow-up observations are needed, and our sensitive spectral measurements on a complete sample of bright SMGs are suitable for providing further constraints.

Our sample include five pairs of SMGs, and their 870 μm continuum images and 3 mm spectra are shown in Figure 8. As stated in Section 3.1 and seen in Figure 8, the spatial resolutions of the naturally weighted cubes are too coarse to separate the SMG pairs of AS2COS0001.1/1.2, AS2COS0008.1/8.2, and AS2COS0028.1/28.2. We therefore reimage those cubes using Briggs weighting with a robust parameter of −1.5, which stands as a good balance between sensitivity and resolution after experimenting with different baseline weightings. To increase the detectability of line emission, we also binned the channels to ∼150 km s⁻¹ channel⁻¹, appropriate for the typical line width of SMGs, which is about 500 km s⁻¹ in FWHM (e.g., Bothwell et al. 2013; Birkin et al. 2021). The spectra are extracted using apertures centered on their reported 870 μm locations (Simpson et al. 2020) with the sizes and shapes of their corresponding synthesized beams, shown in the corners of Figure 8.

Zoom In Zoom Out Reset image size

We obtain a significant detection from AS2COS0008.2 and a marginal detection from AS2COS0001.2. The best-fit single Gaussian model suggests that the redshift of AS2COS0008.2 is 3.5739 ± 0.0009, assuming that it is the same line as AS2COS0008.1, and the significance of the detection is at >4σ based on the velocity-integrated intensity over the FWHM of the best-fit model (Figure 8). Although marginal (∼2σ over the FWHM) in our data, the same line from AS2COS0001.2 (and AS2COS0001.1) was detected at high significance (>5σ) by Jiménez-Andrade et al. (2020) with much deeper data from NOEMA, and Simpson et al. (2020) and Mitsuhashi et al. (2021) both reported detection of the [C ii] line from AS2COS0001.2. Their reported redshifts are consistent with our measurements. For the following analyses, we therefore adopt the redshift measured by Jiménez-Andrade et al. (2020) for AS2COS0001.2, which is 4.633 ± 0.001.

The velocity differences between these two pairs are within 500 km s⁻¹. At the projected distances of 20–30 kpc, their velocity differences would suggest that these two pairs are gravitationally bound systems, assuming that SMGs as suggested by clustering analyses reside within halos with masses of ∼10¹³ M_⊙ (Hickox et al. 2012; Chen et al. 2016; Wilkinson et al. 2017; An et al. 2019; Lim et al. 2020; Stach et al. 2021). The redshift differences of the two pairs are both Δz < 0.02, a definition adopted by some models for physically associated pairs (Hayward et al. 2013a; Muñoz Arancibia et al. 2015). This suggests that 40% of the paired SMGs in our sample are physically associated. In addition, given that our observations are designed to detect lines from bright SMGs, deeper observations may reveal associated line emissions from the secondary SMGs of other pairs, and the true fraction of physically associated pairs could be higher. Indeed, as seen in Figure 8, among these four secondary SMGs, those not detected in line by our data are the faintest in 870 μm flux densities (<2 mJy), and based on recent measurements (Dudzevičiūtė et al. 2020; Birkin et al. 2021), we do not expect to detect lines on these faint SMGs given the sensitivity of our data. Our results are also consistent with Simpson et al. (2020), where, based on the number counts analyses, they estimated that 62% ± 7% of the >3 mJy secondaries are physically associated with the primary SMGs.

While our results suggest a relatively high fraction (≥40%) of physically associated SMG pairs, the current sample size is admittedly too small to have much constraining power. Furthermore, it is also not straightforward to compare our results against model predictions. First of all, the adopted criteria in defining pairs are different in each model. For example, Cowley et al. (2015) made predictions based on the brightest two constituent SMGs. On the other hand, Hayward et al. (2013a) considered all SMGs with 850 μm flux densities of >1 mJy, which is the flux level that is highly incomplete in the parent AS2COSMOS catalog. Second, our spectral measurements are incomplete in terms of including all pairs or multiples of submillimeter sources above a certain flux density limit measured by SCUBA-2, ²¹ which is what model predictions are typically based on. To make fair comparisons, experiments that are tailored for these model predictions are needed.

In spite of these caveats, however, our results do provide constraints on this issue from another angle. Since our target SMGs are complete to their flux cuts, it can be stated that our results suggest that for SMGs brighter than 12.4 mJy at 870 μm, their brightest companion SMGs are physically associated with the primary SMGs ≥40% of the time.

Taking the adopted redshifts in Table 1, the overall median of the 18 primary target SMGs is 3.3 ± 0.3 with a bootstrapping uncertainty that is slightly larger than that of the photometric redshifts (3.1 ± 0.2) simply due to the three outliers whose real redshifts are all higher than their photometric redshifts. Three target SMGs, AS2COS0001.1, AS2COS0002.1, and AS2COS0006.1, have been shown to be part of a z = 4.6 coherent megaparsec-scale structure (Mitsuhashi et al. 2021).

The uncertainty for the median caused by the measurement errors is on the order of 3 × 10⁻⁴, which is the standard deviation of a set of 1000 medians, each of which is a result of a set of randomly perturbed redshifts generated according to the measurements. The uncertainty for the median caused by the wrongly assigned line transitions for the SMGs with a single line detection is estimated to be 0.1, which is again the standard deviation of a set of 1000 medians, each of which is a result of a set of redshifts that are generated based on the adopted redshifts but having randomly chosen two SMGs with a single line that are randomly assigned to one transition up or down from the adopted one, providing that the new transition would not allow the coverage of another CO transition within the bandwidth. In summary, the dominant contributor to the uncertainty of the median is the bootstrapping uncertainty, essentially reflecting the sample size and the underlying distribution.

4.2.1. Comparisons with Previous Measurements

Previously, redshift measurements of SMGs have been mostly focused on those with S_850/870 ∼ 2−10 mJy (e.g., Chapman et al. 2005; Ivison et al. 2007; Wardlow et al. 2011; Simpson et al. 2014; Cowie et al. 2017; Stach et al. 2019; Birkin et al. 2021). This is simply a result of the typical survey depth and area, coupled with the shape of the number counts, that maximizes the number of detections in this flux range. The median redshifts deduced differ slightly in different studies, but they generally lie in the range of 2.4–2.6, with an uncertainty of about 0.1–0.2. Together with our results, this hints that brighter SMGs are, on average, located at higher redshifts.

Since these previous studies used various methods to estimate or measure redshifts, they often covered different flux ranges. To better assess the correlation between flux density and redshift in Figure 9, we plot the most recent results from self-similar and well-defined samples, namely ALESS (da Cunha et al. 2015), AS2UDS (Stach et al. 2019), and AS2COSMOS (Simpson et al. 2020; S. Ikarashi et al. in preparation). We investigate the photometric and spectroscopic redshifts separately. In Figure 9, the photometric redshifts of all three samples are combined and shown as density contours. In both cases, the correlation between flux density and redshift can be observed. Indeed, the Spearman correlation coefficients are 0.3 and 0.4 for photometric and spectroscopic redshifts, respectively, and the p-values are ≲0.001. The maximum-likelihood linear fitting yields z_median = (2.2 ± 0.3) + (0.09 ± 0.02) × S₈₇₀ for the sample with spectroscopic redshift, consistent with the best-fit based on photometric redshifts. Our results are also consistent with previous assessments regarding either the general trends or the slope (Ivison et al. 2002; Cowie et al. 2017; Stach et al. 2019; Simpson et al. 2020; Birkin et al. 2021).

Zoom In Zoom Out Reset image size

Finally, we compare our results to the South Pole Telescope (SPT) SMGs presented by Reuter et al. (2020), where complete measurements of spectroscopic redshifts are reported on a well-defined millimeter-selected SMG sample. The reported SPT SMGs have apparent flux densities of S₈₇₀ > 25 mJy, and they are mostly strongly lensed. In Figure 9, we show their redshift distributions with respect to the intrinsic flux densities where gravitational magnifications are accounted for. Interestingly, the distribution of the SPT SMGs appears to follow our unlensed SMGs at S₈₇₀ ≳ 10 mJy, below which the known selection effect of lensing biasing sources toward higher redshifts can be clearly seen.

4.2.2. Comparisons with Models

We now compare our measurements to model predictions. Currently, there are various schools of models that can reproduce basic properties of SMGs, such as number counts, and have made predictions for other properties, such as redshifts. They include empirical models (e.g., Béthermin et al. 2012, 2017; Casey et al. 2018; Zavala et al. 2021); semiempirical models (e.g., Hayward et al. 2013b; Popping et al. 2020); semianalytical models (e.g., Lacey et al. 2016; Lagos et al. 2020); hybrid models, including both analytical and empirical recipes (e.g., Negrello et al. 2007, 2017; Cai et al. 2013); and hydrodynamical simulations with various treatments of dust and its emissions (e.g., McAlpine et al. 2019; Hayward et al. 2021; Lovell et al. 2021).

In general, while good progress has been made over the years, it is still difficult to find in hydrodynamical simulations as many bright SMGs as are detected in observations. For example, imposing the flux cut of our sample, there are only two such sources in the model presented by Lovell et al. (2021), and similar or fewer numbers are seen in other hydrodynamical models (McAlpine et al. 2019; Hayward et al. 2021). The lack of SMGs brighter than our flux cut could be due to the fact that current simulations do not have sufficient volumes, which are typically about 100 cMpc on a side. As a result, for a meaningful comparison, we do not consider predictions based on hydrodynamical simulations in this section.

We list the model catalogs used for comparisons and their details in Table 2, including the area from which the catalogs are drawn and the methods used to generate these catalogs. Briefly, all of these models employ certain methods to predict SFRs, and thus total infrared luminosity (L_IR), and adopt various templates or relations to model the infrared and submillimeter SEDs. For example, analytical solutions of SFRs are calculated in physically motivated models like SHARK and Cai–Negrello, where the predicted submillimeter flux densities are modeled using selected SED templates. For models like SIDES and Popping, SFRs are estimated by fitting their dark matter base empirical models to stellar mass functions with considerations of the star formation main sequence. SIDES then adopts an SED template to model submillimeter fluxes, where Popping applies a relation between the 850 μm flux density and SFRs and dust masses based on dust radiative transfer codes. Casey–Zavala obtains L_IR directly by modeling the infrared luminosity functions, and a modified blackbody plus mid-infrared power-law SED model was assumed, together with a relation between dust temperature and L_IR. The full descriptions of the methodologies used to generate these model catalogs can be found in the latest reference for each model listed in Table 2. However, note that for Cai–Negrello, the model has been slightly updated with a new version of the SED template, resulting in a better agreement with the recently reported number counts in the far-infrared and submilllimeter (M. Negrello et al. 2022, in preparation).

Table 2. Model Catalogs Used for Comparisons

Model	Method	Area [deg2]	Lensing	References
Cai–Negrello	Analytical + empirical	500	Yes	Cai et al. (2013), Negrello et al. (2017)
SIDES	Empirical	2	Yes	Béthermin et al. (2017)
Casey–Zavala	Empirical	10	No	Casey et al. (2018), Zavala et al. (2021)
Popping	Semiempirical	22.2	No	Popping et al. (2020)
SHARK	Semianalytical	107.9	No	Lagos et al. (2019, 2020)

Download table as:  ASCIITypeset image

We assess their redshift distributions as follows. For SIDES, since it is drawn from an area that has a size similar to our initial SCUBA-2 survey, we simply impose the flux cut of our sample and plot the cumulative distribution in Figure 10. For model catalogs that are drawn from areas that are much larger than 2 deg², we randomly draw the model catalogs 100 times based on the size of our surveyed area. For each draw, we assess the cumulative distribution, and in Figure 10, we plot the median of all 100 draws, as well as their 10th–90th percentile range, to roughly indicate the fluctuations of these distributions, which are mainly due to their field-to-field variations and shot noise.

Zoom In Zoom Out Reset image size

To compare our observational results to the models, we generate a set of 100 randomly perturbed redshift distributions based on our measurements and their associated uncertainties, as well as the probability of wrongly assigned redshifts for SMGs with only a single line detection, following the procedure described in the previous paragraph. The results are plotted in Figure 10, again with a 10th–90th percentile range. We then perform a two-sample Kolmogorov–Smirnov (K-S) test between each of the 100 perturbed observational distributions and each of those 100 drawn from each of the model catalogs, for a total of 10,000 comparisons, except for SIDES, which only has 100. For each model catalog, we record the p-values of each comparison and compute the percentage of them having p < 0.05, a standard value that suggests that the difference between the two samples is highly significant. Changing this p-value cut does not affect the conclusion. We show the results of these comparisons in Figure 10. Strikingly, all models predict lower redshifts compared to our measurements, although from the statistical point of view, the significance of their differences differs. Our measurements agree best with SIDES; however, SIDES only provides one sight line. With a much larger volume, perhaps the comparison against SIDES would yield a result that is similar to that from the comparison against Cai–Negrello. The most significant difference comes from the comparison against SHARK, where all sight lines that we draw do not agree with our measurements.

We can also compare them in median redshifts, which are 1.7 ± 0.3 for SHARK, 2.9 ± 0.7 for SIDES, 2.6 ± 0.3 for Casey–Zavala, 2.1 ± 0.4 for Popping, and 2.9 ± 0.4 for Cai–Negrello. The error budget again includes both statistical and systematic uncertainties caused by field-to-field variations and the shot noise of the simulations. Our measured median redshift of 3.3 ± 0.3 is larger compared to all models, with the least agreement with SHARK, which is in 3.8σ tension. These are the same conclusions we make based on more detailed two-sample K-S tests. Measurements based on a much larger sample size (i.e., a factor of a few more) would allow meaningful tests on the surviving models.

4.2.3. Implications

Strong gravitational lensing (lensing magnification μ > 2) is found to be responsible for the exceptional apparent brightness (S_500μm ≳ 100, S_1.4mm ≳ 10, or S_850μm ≳ 30 mJy) of the majority of the submillimeter and millimeter sources that are uncovered by shallow but wide surveys such as those conducted by the SPT and Herschel (Negrello et al. 2010, 2017; Bussmann et al. 2013; Hezaveh et al. 2013; Wardlow et al. 2013; Spilker et al. 2016). Most models predict that the fraction of sources experiencing strong lensing drops significantly, from about ≳70% to almost zero, within a narrow range of flux density, which is about 80–100 mJy at 500 μm and 10–30 mJy at 850 μm (Wardlow et al. 2013; Negrello et al. 2017). Our statistical sample SMGs have their flux densities within this range and therefore represent an opportunity to test model predictions.

To estimate the lensing magnifications of our sample SMGs, we employ lenstool (Kneib et al. 1996; Jullo & Kneib 2009). For each SMG, we utilize the COSMOS15 catalog (Laigle et al. 2016) and construct a list of neighboring galaxies that are located within 30'' of the SMG and have had reported estimates of photometric redshifts and stellar masses. This angular distance is chosen because beyond it, the results do not change significantly. A gallery of the selected massive subsets of these neighboring galaxies is plotted for each SMG in Figure 11. To construct a model of the foreground gravitational potential for each SMG, we assume a regular Navarro–Frenk–White (NFW) density profile (α = 1; Navarro et al. 1997) for each neighboring galaxy in the list. Besides sky position, other key input parameters for the model of each neighboring galaxy include ellipticity, position angle, scale radius r_s , and characteristic velocity σ_s (Limousin et al. 2005). We describe the procedures for obtaining these parameters in the following.

Zoom In Zoom Out Reset image size

We first run SExtractor (Bertin & Arnouts 1996) to obtain the structural parameters ellipticity and position angle, and the primary image considered is the K-band image from UltraVista DR4, except for AS2COS0031.1, where it is outside the UltraVista footprint, so we adopt the 3.6 μm image from Spitzer. The scale radius and characteristic velocity are related to the assumed NFW profile, which is expressed as

$\begin{eqnarray*}&&\rho (r)=\displaystyle \frac{{\rho }_{s}}{(r/{r}_{s}){\left(1+r/{r}_{s}\right)}^{2}},\end{eqnarray*}$

where r_s is the scale radius and ρ_s is the characteristic density so that the characteristic velocity is σ_s = 4/3Gr_s ² ρ_s . Here ρ_s can be related to critical density ρ_c as ρ_s = δ_c ρ_c , where the density contrast ${\delta }_{c}=200{c}^{3}/(3\mathrm{ln}(1+c)-3(c/1+c))$, and c is the concentration factor defined as the ratio between the halo radius r₂₀₀ (the mean density within which is 200ρ_c ) and r_s . The halo mass can therefore be deduced as ${M}_{200}=800/3\pi {r}_{200}^{3}{\rho }_{c}$. That is, r_s and σ_s can be derived for a given combination of concentration factor and halo mass.

For each neighboring galaxy, we estimate its halo mass by adopting the stellar mass and redshift provided by the COSMOS15 catalog and applying the latest stellar-to-halo mass relationships given by Legrand et al. (2019). For concentration, we refer to the classical work of Wechsler et al. (2002) based on dark matter N-body simulations of ΛCDM cosmology, where the relationships between concentration, halo mass, and redshift are shown.

Each of the adopted input parameters has its own reported uncertainty or scatter, and they should be taken into account in order to properly assess the probability distributions of lensing magnifications. To do so, for each SMG, we create 100 potential models, where each constituent potential is created by randomly perturbing each of their parameters based on their quoted uncertainties. We run lenstool on these 100 potential models for each SMG and calculate the median and 1σ percentile. The results are given in Table 3.

Table 3. Lensing Magnifications

ID	μ
AS2COS0001.1	${1.4}_{-0.1}^{+0.1}$
AS2COS0002.1	${3.0}_{-0.7}^{+1.4}$
AS2COS0006.1	1.1
AS2COS0008.1	1.1
AS2COS0009.1	1.0
AS2COS0011.1	1.1
AS2COS0013.1	1.1
AS2COS0014.1	${1.7}_{-0.2}^{+0.5}$
AS2COS0023.1	${1.7}_{-0.1}^{+0.1}$
AS2COS0028.1	1.1
AS2COS0031.1	1.1
AS2COS0037.1	1.1
AS2COS0044.1	1.1
AS2COS0054.1	1.1
AS2COS0065.1	1.0
AS2COS0066.1	1.1
AS2COS0090.1	1.1
AS2COS0139.1	1.1

Note. Details can be found in Section 4.3. Uncertainties less than 0.05 are omitted.

Download table as:  ASCIITypeset image

As a check, the lensing magnification of AS2COS0001.1 was also estimated by Jiménez-Andrade et al. (2020) with a different assumption of the density profile, and they found a factor of 1.5, consistent with our results. Notably, we find that if we only look at the mass distributions of much closer surroundings, the magnifications of a few SMGs, such as AS2COS0002.1 and AS2COS0014.1, would be significantly underestimated. This is because in these cases, there are a few very massive (up to 10^11.5 M_⊙) foreground galaxies, and in the case of AS2COS0002.1, a number of galaxies, including the very massive ones, are located at z ∼ 1, suggesting a group environment.

Based on Table 3, we can see by looking at their medians that only one out of the 18 SMGs in our sample can be classified as being strongly lensed with μ > 2, suggesting a relatively low strongly lensed fraction of <10% compared with the brighter sources. This is in line with what has been predicted by some observations (Chapman et al. 2002) and models. We note that our methodology in constructing foreground gravitational potentials does not account for contributions from halos larger than galactic halos (i.e., group-scale halos); thus, it is possible that the magnifications of some SMGs may have been underestimated by a moderate amount. On the other hand, our methodology also does not account for underdense regions, so the overall magnification is slightly biased high, and as a result, many unlensed SMGs have μ = 1.1. Nevertheless, these two factors are not expected to affect our results significantly.

We compare our results with predictions of magnification for the flux range of our sample from SIDES and Cai–Negrello in Figure 12 and find that our results are inconsistent with SIDES but consistent with Cai–Negrello. This can be understood by the fact that, while the Cai–Negrello model takes into account both the redshift distribution of model SMGs and the mass distributions in the foreground, the predictions from SIDES are simply drawn from some probability distributions according to the redshifts of the model SMGs, as stated in Béthermin et al. (2017).

Zoom In Zoom Out Reset image size

We now expand the comparisons by including published measurements for a wide range of 870 μm flux densities, including the Submillimeter Array (SMA) and ALMA follow-up studies of Herschel-selected sources (Bussmann et al. 2013, 2015), as well as those based on the SPT survey (Spilker et al. 2016). For each of the SPT sources, separate magnification factors were reported for each of their constituent components, so we compute their flux-weighted magnification as the nominal magnification for each SPT source (Figure 3 in Spilker et al. 2016). It is also worth noting that while the Herschel-selected sample observed by the SMA is a flux-limited sample, the one observed by ALMA was selected to have a nearby foreground source to maximize the probability of detecting strongly lensed objects. As a result, the Herschel-selected sample observed by ALMA is likely biased high in lensing magnification for the same flux range.

The results are plotted in Figure 13, where we find that despite the possibility of having a bias, the Herschel–ALMA measurements are consistent with our results at the overlapping flux density range (10–20 mJy). We can also compare to the model predictions from Cai–Negrello. Note that in the Cai–Negrello model, sources that are not strongly lensed (μ < 2) are assigned 1 for their magnification. We therefore adopt the same method for the measurements and compute the binned median with bootstrapped errors. We compare our results with the median of the model prediction with the same flux bins and find that they have good agreement overall.

Zoom In Zoom Out Reset image size

Last but not least, having lensing magnification measurements of a flux-limited sample means that we can compute the number counts of the strongly lensed SMGs. The cumulative number counts of the parent AS2COSMOS sample were presented in Simpson et al. (2020), and they are plotted in Figure 14. Since our sample is limited to 12.4 mJy at 870 μm, the cumulative counts for the strongly lensed SMGs can be simply scaled from the parent counts of the corresponding flux bin (12.1 mJy in this case) by dividing its value by our sample size. There is no source between 12.1 and 12.4 mJy, so no correction needs to be applied. To estimate the uncertainty, the Poisson error of 1 is adopted and propagated based on the uncertainty of the parent counts. The results are shown in Figure 14 as the red point. As a check, the value is very close to one divided by the size of the footprint of the primary S2COSMOS catalog, so one divided by 1.6 deg². Overall, we again find good agreement with the model prediction of Cai–Negrello.

Zoom In Zoom Out Reset image size

As shown in Section 3.2, we find that the emission lines of 15 out of the 17 (∼90% ± 30%) primary SMGs with line detections appear to be better modeled by double Gaussian profiles. This high fraction of SMGs exhibiting double Gaussian line profiles is in contrast with previous studies of other SMGs samples, which have typically reported fractions under 50% (Greve et al. 2005; Bothwell et al. 2013; Birkin et al. 2021). This could be in part due to our smaller sample size, or it could be because of our sample SMGs being brighter so that they may appear dynamically distinct from their fainter counterparts. Unless dispersion-dominated, it is expected that the lines would appear as double Gaussian in systems that exhibit velocity gradients due to either orderly rotating disks or mergers. Spatially resolved dynamical studies are needed to make a more definite conclusion on their physical origins.

On the other hand, we find that in our sample, for lines that are better described by a double Gaussian, the median separation in velocity between the red and blue peaks is 430 ± 40 km s⁻¹, which is consistent with previous studies using their double-peaked SMGs (Bothwell et al. 2013; Birkin et al. 2021), as well as recent spatially resolved observations of CO/[C ii] lines on a few SMGs (Chen et al. 2017; Calistro Rivera et al. 2018; Litke et al. 2019). This measurement of red and blue velocity separations allow us to crudely estimate dynamical masses, provided there are some reasonable assumptions. The dynamical mass estimate can be expressed as ${M}_{\mathrm{dyn}}=C{({v}_{c}/\sin (i))}^{2}\times R/G$, where G is the gravitational constant, v_c is the circular velocity, R is the representative radius, C is the correction factor depending on the mass distributions and the adopted representative radius, and i is the inclination angle. This assumes a rotation-dominated system, which is supported by recent high-resolution and high-sensitivity dynamical measurements of SMGs (e.g., De Breuck et al. 2014; Lelli et al. 2021; Rizzo et al. 2021). Assuming an exponential disk profile with a scale length h and a half-light radius r_e , it has been shown that R = 3h = 1.79r_e ²² is a good approximation for spatially unresolved CO observations (de Blok & Walter 2014). Under this assumption of an exponential disk, C is about 2 (Binney & Tremaine 2008). We adopt r_e = 3 kpc based on recent spatially resolved CO measurements of SMGs (Chen et al. 2017; Calistro Rivera et al. 2018), and v_c is half of the median velocity separation that we measure. Since v_c in our case is an average over a sample of sources with random inclination angles, we adopt an average $\langle {\sin }^{2}(i)\rangle =2/3$ (Tacconi et al. 2008). Since the dynamical mass estimated this way only includes mass within the assumed r_e , we multiply by a factor of 2 to obtain the total dynamical mass. As a result, we estimate a total dynamical mass of M_dyn = 3.5 ± 0.6 × 10¹¹ M_⊙. Accounting for the typical dark matter fraction at high redshifts, which was recently reported as 12% by Genzel et al. (2020), we deduce a total baryon mass of M_baryon = 3.1 ±0.6 × 10¹¹ M_⊙. Note that the observational constraints on the dark matter fraction are obtained from galaxies at z ∼ 2, slightly lower than our sample SMGs. However, recent IllustrisTNG hydrodynamical simulations have shown that the dark matter fractions are expected to be lower at higher redshifts (Lovell et al. 2018). On the other hand, the simulations have predicted about a factor of 2 higher fractions than those obtained from observations. It is outside the scope of this paper to argue which is correct, although, by adopting either a factor of 2 higher dark matter fraction or no dark matter, our conclusions would not be significantly altered. Finally, in a follow-up paper, where we plan to present the properties of the ISM via CO/[C i] and dust SED measurement and modeling (C.-L. Liao et al. 2022, in preparation), we will show that by including stellar, gas, and dust mass, the average total baryon mass is consistent with the dynamical estimates presented in this work.

Following Birkin et al. (2021), the estimates of total baryon mass and measurements of velocity dispersion via molecular lines allow us to assess the connection between our sample SMGs and their proposed descendants at lower redshifts, such as compact quiescent galaxies at similar redshifts or local massive early-type galaxies (ETGs; Lilly et al. 1999; Swinbank et al. 2006; Alexander & Hickox 2012; Toft et al. 2014; Dudzevičiūtė et al. 2020). This is under the assumption that the total mass and dynamical energy are mostly preserved along the evolution.

In Figure 15, we plot our sample as a combined median, where σ is the median of all of the line moment measurements presented in Section 3.2, which is 260 ± 21 km s⁻¹. As a comparison, measurements of local ETGs that are based on the MASSIVE (Davis et al. 2019) and ATLAS^3D (Cappellari et al. 2013) surveys are shown, as well as the quiescent galaxies (D_n 4000 > 1.55; Kauffmann et al. 2003; Wu et al. 2018) at z = 0.6–1.0 based on the LEGA-C survey (van der Wel et al. 2021; Wu et al. 2021) and some compiled z ∼ 2−4 spectroscopically confirmed quiescent galaxies (Stockmann et al. 2020; Valentino et al. 2020). The results from the similar analyses of mostly fainter SMGs presented by Birkin et al. (2021) are also shown. For ETGs and high-redshift quiescent galaxies, we adopt their stellar masses for the total baryon masses, given that stellar mass dominates baryon mass in these populations, whereas for fainter SMGs in Birkin et al. (2021), stellar, gas, and dust masses are summed to deduce their baryon masses. Our sample of bright SMGs lies in the regions occupied by local ETGs and high-redshift quiescent galaxies, corroborating the proposed evolutionary link among these populations.

Zoom In Zoom Out Reset image size

In addition, given that our sample is well defined in both flux density and area, we can compute the number densities of our sample SMGs and compare the results with those based on recently confirmed quiescent galaxies at similar redshifts, especially those at z > 4. The total number of primary SMGs in the redshift bins of 2.1 < z < 3, 3 < z < 4, and 4 < z < 5.1 is 6, 6, and 5, respectively. The lower bound of the lowest redshift bin is determined by the upper bound of the redshift gap (Section 3.3), and the upper bound of the highest redshift bin is justified by the fact that, given our wavelength coverage, emission lines should have been detected from sources at z = 5.1–7.2. The sky area covered by our initial SCUBA-2 survey is 1.6 deg² (Simpson et al. 2019), so the volumes covered in the three redshift bins are all approximately 2.0 × 10⁷ cMpc³. Thus, the observed number density of SMGs with S₈₇₀ > 12.4 mJy is about 3 × 10⁻⁷ cMpc⁻³ for all three redshift bins. However, assuming a typical lifetime of 100 Myr for our sample SMGs (C.-L. Liao et al. 2022, in preparation), the duty cycle corrected number densities are 3.3^+2.8 _−1.8 × 10⁻⁶, ${1.9}_{-1.1}^{+1.6}\times {10}^{-6}$, and ${1.0}_{-0.6}^{+0.9}\times {10}^{-6}$ cMpc⁻³ for the three redshift bins, respectively. The uncertainties include Poisson errors and a 40% field-to-field variation (Simpson et al. 2020). Our number density estimate at z > 4 is similar to that presented by searches based on Herschel (Ivison et al. 2016) and the ALMA 2 mm survey (Manning et al. 2022), although the sample selections are different and thus direct comparisons are not straightforward.

There has been an active discussion regarding the progenitors of the recently confirmed z ∼ 3−4 quiescent galaxies, and extreme dusty star-forming galaxies such as SMGs at z > 4 have been proposed to be ideal candidates, given their nature (e.g., Straatman et al. 2014; Valentino et al. 2020). From the number density, our results at z > 4 appear to be in agreement with those of the z ∼ 3−4 quiescent galaxies estimated by Davidzon et al. (2017) and Girelli et al. (2019; ∼10⁻⁶ cMpc⁻³) but lower by a factor of a few to 10 compared to other estimates (∼10⁻⁵ cMpc⁻³; Straatman et al. 2014; Schreiber et al. 2018; Merlin et al. 2019). The different number densities estimated for quiescent galaxies could be due to selection effects, such as different mass limits, as well as the methods to classify quiescent galaxies. However, given that we are selecting the brightest SMGs, we expect our number density estimates to be lower limits, and fainter SMGs at z > 4 are confirmed to exist (Section 4.2.1). We can conclude that our sample of SMGs has the properties expected for the progenitors of some z ∼ 3−4 quiescent galaxies, especially the most massive ones, with stellar masses of >10¹¹ M_⊙.