Comprehensive Gas Characterization of a z = 2.5 Protocluster: A Cluster Core Caught in the Beginning of Virialization?

The data were taken with the VLA in D configuration in the Ka band (33 GHz). The observations as part of VLA program 15B-210 (PI: C. Casey) were taken in eight epochs in 2015 November, with seven 6 hr integrations and one 4 hr integration, for a total of 46 hr. Each of the two pointings was observed for 12.5 hr on source. We used the four 8-bit samplers (each covering a bandwidth of 1024 MHz) with 15 MHz binning in 128 MHz sidebands. The sampler pairs, A0/C0 and B0/D0, were tuned to 33.25 and 33.20 GHz, respectively. Thus, the effective frequency coverage is 32.35−33.29 GHz, corresponding to 2.46 < z < 2.56 for the ground-state transition of carbon monoxide, ¹²CO (J = 1–0).

Pointing 1 is centered at α = 10:00:57.30 δ = +02:20:13.3, selected to include a Scuba-2-selected galaxy named 450.28 (Casey et al. 2013). Prior to high-resolution ALMA follow-up, 450.28 was identified at a redshift of z = 2.47 via Hα emission (Casey et al. 2015); we now know that this spectroscopic redshift is only one of multiple sources at similar redshifts. 450.28 is also associated with an extremely bright (S_350μm ∼ 90 mJy), unresolved Herschel/SPIRE counterpart. We combine the total VLA time on source for this pointing (12.5 hr) with four sets of observations from VLA program 15B-290 (PI: T. Wang) for an effective integration time of 26 hr. Pointing 2 is centered at α = 10:00:34.6 δ = +02:22:01.20 to include another Scuba-2-selected galaxy from Casey et al. (2013) named 450.58. Both pointings were selected for their relative spatial overdensities of DSFGs and LBGs with reliable redshifts at 2.46 < z < 2.55. We use the Common Astronomy Software Applications (CASA; McMullin et al. 2007) ⁸ and the 2016 VLA pipeline to reduce the full data set. The data cubes have a channel width of 70 km s⁻¹ and a pixel scale of 042 pixel^–1 ( ≈ 3.5 kpc at z = 2.5). The images are produced via the task clean. We use natural weighting to maximize the signal-to-noise ratio (S/N) of expected point sources, with a final synthesized beam size of 287 × 256. We trim the images to 40% of the primary beam response (FWHP size $=1\buildrel{\,\prime}\over{.} 35$). We later find that the CO-identified sources are very close together with complex morphology, so we repeat the analysis in both cleaned and dirty maps to account for potential cleaning artifacts. The typical rms per channel is 10 μJy beam⁻¹.

3.2.1. Band 6 Data

3.2.2. Band 3 Data

Since this protocluster lies in the COSMOS field, the CO and dust continuum observations presented here are complemented by a wealth of multiwavelength data taken as part of the COSMOS survey (Capak et al. 2007; Scoville et al. 2007). Where the sources can be resolved individually, we match CO-identified sources to the COSMOS 2015 catalog (Laigle et al. 2016) with a match radius of 15 for more complete optical/near-infrared (OIR) constraints.

In Section 4, we model the individual OIR spectral energy distributions (SEDs) using magphys (da Cunha et al. 2008, 2015) using the following photometry: Subaru Suprime-Cam broadband and narrowband filters; UltraVISTA YJHK_s ; Spitzer/IRAC 3.5, 4.5, 5.8, and 8.0 μm; and Spitzer/MIPS 24 μm.

We also search for individual counterparts in the two major radio continuum surveys included in COSMOS: the VLA-COSMOS 1.4 GHz Deep project (with a spatial resolution of 15 and depth 7 μJy beam⁻¹; Schinnerer et al. 2007), and the 3 GHz project (angular resolution 07 and depth 2.3 μJy beam⁻¹; Smolčić et al. 2017). These radio counterparts are meant to account for the possible contribution to the infrared flux from radio AGNs, as well as star formation-generated radio emission.

Additionally, there are Herschel PACS (100 μm and 160 μm) and SPIRE (250 μm, 350 μm, and 500 μm) data available that we use to model the FIR SED of the protocluster core. However, the individual protocluster sources are confused by SPIRE's large beam (36'' at 500 μm), so we fit the FIR SED treating the core as a single source and calculate associated physical properties in aggregate.

Finally, we check the Chandra point-source catalog from Civano et al. (2016) for potential X-ray AGNs. We also check for extended X-ray emission using the merged COSMOS XMM+Chandra maps for consistency with the analysis of Wang et al. (2016).

We search for CO (1–0) sources using the positions of previously identified sources in COSMOS as a prior. Where available, we use the Hα-detected sources from Casey et al. (2015) and Band 6 continuum sources (Zavala et al. 2019) as positional priors, for a total of 12 priors. We searched for additional sources (i.e., not based on these priors) in the dirty cube by iteratively creating moment 0 maps integrated along five consecutive channels at a time, corresponding to an effective line width of 350 km s⁻¹ (typical for CO lines at high redshift; Bothwell et al. 2013). In some cases, we find that the CO detections overlap spatially when integrated, so we use their offsets in frequency to distinguish between sources, and we count the number of detections by the number of discrete bright peak pixels.

In total, we report seven CO (1–0) detections in Pointing 1 and one in Pointing 2 (Table 1). Note that all CO sources reported in Pointing 1 have also been confirmed at shallower depth in other publications (Wang et al. 2016, 2018; Gómez-Guijarro et al. 2019), though each study uses different S/N and size thresholds, so the number of independent sources may vary. At this spatial resolution, the sources are mostly unresolved, but some flux is missed in a beam-sized aperture. We therefore extract spectra using a manually placed aperture that encompasses flux from adjacent pixels that all have S/N > 3 (typically the size of 1–2 VLA beams), taking care to shape the apertures such that we do not double-count flux from nearby sources. In Appendix A.1 we describe in more detail our source extraction method and test whether the sources are spatially resolved.

Table 1. Details of CO Data and Derived Gas Properties

Name a	R.A.	Decl.	z	ICO (1−0)	FWHM	S/N b	log(L ${}_{\mathrm{CO}(1-0)}^{{\prime} }$)	M ${}_{\mathrm{gas}}$ c	ICO (3−2)	FWHM	S/N	log(L ${}_{\mathrm{CO}(3-2)}^{{\prime} }$)	r ${}_{31}$ d
	(J2000)	(J2000)		(mJy km s−1)	(km s−1)		(K km s−1 pc2)	(1010 M⊙)	(mJy km s−1)	(km s−1)		(K km s−1 pc2)
ID1	150.237350	2.338078	2.494	122.0 ± 6.4	486 ± 36	16.04	10.55 ± 0.02	${22.96}_{-0.74}^{+0.74}$	1261.0 ± 143.0	532 ± 78	11.26	10.61 ± 0.05	1.14 ± 0.14
ID2	150.237333	2.336935	2.496	39.0 ± 3.9	532 ± 54	4.92	10.05 ± 0.05	${7.35}_{-0.45}^{+0.45}$	< 190	< 500	⋯	< 9.79	< 0.54
ID3	150.238529	2.336777	2.503	46.3 ± 6.9	1271 ± 143	3.96	10.13 ± 0.07	${8.76}_{-0.81}^{+0.81}$	⋯	⋯	⋯	⋯	⋯
ID4	150.236912	2.335836	2.504	92.6 ± 12.6	1095 ± 135	7.74	10.43 ± 0.06	${17.53}_{-1.47}^{+1.47}$	⋯	⋯	⋯	⋯	⋯
ID5 e	150.228896	2.329801	2.507	37.4 ± 5.1	172 ± 71	9.97	10.04 ± 0.06	${7.10}_{-0.59}^{+0.59}$	⋯	⋯	⋯	⋯	⋯
ID6	150.239017	2.336322	2.509	50.0 ± 7.0	981 ± 145	4.81	10.16 ± 0.07	${9.50}_{-0.82}^{+0.82}$	⋯	⋯	⋯	⋯	⋯
ID7	150.239917	2.336443	2.511	64.1 ± 6.4	619 ± 56	8.62	10.27 ± 0.05	${12.19}_{-0.75}^{+0.75}$	⋯	⋯	⋯	⋯	⋯
450.58	150.150250	2.364176	2.464	180.5 ± 19.1	334 ± 59	9.87	10.70 ± 0.03	32.94 ± 2.64	< 1400	< 500	< 2	< 10.65	< 0.90

Note.

^a Names beginning with "ID" are labeled in order of increasing redshift in Pointing 1; 450.58 is in Pointing 2. ^b S/N is the ratio of integrated moment 0 flux density to the rms in the moment 0 maps. ^c M_gas is calculated assuming α_CO = 6.5 M_⊙/(K km s⁻¹ pc²), as advocated for in Appendix A.3 rather than the typical value of α_CO = 1 assumed for high-redshift DSFGs (e.g., Carilli & Walter 2013). ^d Ratio of ${L}_{\mathrm{CO}(3-2)}^{{\prime} }$/L ${}_{\mathrm{CO}(1-0)}^{{\prime} }$. ^e Identified as 850.82 in Casey et al. (2013).

Download table as:  ASCIITypeset image

In Figure 2, we show zoomed-in moment 0 maps from the VLA at four velocity slices with respect to a central redshift of z = 2.503, plus ALMA continuum and HST maps. While the dust emission in each source is relatively compact, the marginally extended nature of very nearby CO (1–0) sources combined with a heterogeneous set of line profiles points to a complicated overall molecular gas morphology. The apparent extended nature of the molecular gas along the frequency/velocity axis could imply either that there is significant gas interaction between galaxies or that we are tracing individual peculiar velocities in a single halo; we discuss this further in Section 6. Figure 3 shows the full zoomed-out moment 0 maps overlaid on NIR images from UltraVISTA at the same velocity slices, along with their corresponding spectra, where one can clearly see the ambiguous and overlapping nature of the molecular gas reservoirs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We then calculate the spectral line moments to derive molecular gas masses from the spectra shown in Figure 4. None of the spectra show evidence of continuum emission, so we do not perform any continuum subtraction. Some CO lines appear to be composed of both broad and narrow components or double-peaked lines. Therefore, we directly calculate the spectral line moments, allowing us to constrain the redshift, FWHM, and total line intensity without having to assume that each CO line is composed of a single Gaussian component. In Appendix A.2 we outline these calculations in detail, roughly following the methods of Bothwell et al. (2013).

Zoom In Zoom Out Reset image size

After measuring the line intensity I_CO, the line luminosity ${L}_{\mathrm{CO}}^{{\prime} }$ is the following:

$$\begin{eqnarray}&&{L}_{\mathrm{CO}}^{\mbox{'}}=3.25\times {10}^{7}\times {S}_{\nu }{\rm{\Delta }}v\frac{{D}_{L}^{2}}{{(1+z)}^{3}{\nu }_{\mathrm{obs}}^{2}}{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2},\end{eqnarray} \tag{ 1 }$$

where S_ν Δv = I_CO is measured in Jy km s⁻¹, the luminosity distance D_L in Mpc, and ν in GHz.

Finally, we calculate the molecular gas mass, ${M}_{\mathrm{gas}}={\alpha }_{\mathrm{CO}}L{{\prime} }_{\mathrm{CO}}$, where α_CO = 6.5 M_⊙/(K km s⁻¹ pc²) is the empirical factor to convert to gas mass from line luminosity. In Appendix A.3 we motivate this choice of α_CO, which is higher than the value of α_CO ≈ 1 M_⊙/(K km s⁻¹ pc²) that has historically been assumed for high-redshift galaxies.

We repeat the search for line detections in the ALMA Band 3 data following the procedures outlined in Section 4.1. Since the observations were tuned to a central redshift of z = 2.47, only two galaxies in VLA Pointing 1 and one in Pointing 2 have the frequency coverage to detect CO (3–2). Of these, only ID1 is detected in CO (3–2), with S/N of 9. We report the detections and upper limits in Table 1. We find that the ratio of (3–2) to (1–0) line luminosities r₃₁ = 1.37 ± 0.20 for ID1 and an upper limit of r₃₁ < 0.72 for ID2, suggesting a thermalized state of the ISM in line with r₃₁ = 0.90 ± 0.40 found for combined populations of submillimeter star-forming and AGN host galaxies at 2 < z < 3 (Sharon et al. 2016; Kirkpatrick et al. 2019). The CO (3–2) emission in ID1 appears slightly resolved by ALMA and more compact than the CO (1–0) emission (see Figure 5), although since the CO (1–0) sources are unresolved, higher-resolution observations are needed to verify this. Slightly more compact emission in CO (3–2) is expected since the more highly excited gas is likely isolated to small molecular clouds throughout the galaxy, but it is more likely due to the poor resolution of the VLA, since AlMA's beam size (θ ≈ 09) corresponds to roughly 8 kpc at this redshift.

Zoom In Zoom Out Reset image size

We use the ALMA Band 6 map centered on Pointing 1 to derive millimeter continuum flux densities. Using the CO positions as priors, we search the non-primary-beam-corrected maps for sources with S/N > 3. ID5 lies outside of the Band 6 primary beam. Of the remaining six sources, five are detected at 1.2 mm. Dust continuum detections at this wavelength correspond to λ ≈ 300 μm, which is on the Rayleigh–Jeans tail of the blackbody emission, where the dust emission is assumed to be optically thin. We can then calculate dust masses directly since it is proportional to the observed millimeter flux and the dust temperature, which we do using the framework described for 850 μm sources in Dunne et al. (2000):

$$\begin{eqnarray}&&{M}_{\mathrm{dust}}=\frac{{S}_{{\nu }_{\mathrm{obs}}}{D}_{L}^{2}{\left(1+z\right)}^{-(3+\beta )}}{\kappa ({\nu }_{\mathrm{ref}}){B}_{\nu }({\nu }_{\mathrm{ref}},{T}_{\mathrm{dust}})}{\left(\frac{{\nu }_{\mathrm{ref}}}{{\nu }_{\mathrm{obs}}}\right)}^{2+\beta }\left(\frac{{{\rm{\Gamma }}}_{\mathrm{RJ}(\mathrm{ref},0)}}{{{\rm{\Gamma }}}_{\mathrm{RJ}}}\right).\end{eqnarray} \tag{ 2 }$$

The emissivity, β, is fixed at 1.8. The quantity κ_ν is the dust mass absorption coefficient evaluated at rest-frame ν_ref; here we use κ (450 μm) = 1.3 ± 0.2 cm² g⁻¹ (Li & Draine 2001), which is 1.5 mm in the observed frame. B_ν is the Planck function evaluated at mass-weighted dust temperature T_d , which we fix at 25 K as in Scoville et al. (2016). Finally, Γ_RJ is the Rayleigh–Jeans correction factor. The quantities derived from these fluxes are listed in Table 2, and the maps are displayed in the middle panels of Figure 2.

We fit OIR SEDs to those spectroscopic sources with matches within 15, in the Laigle et al. (2016) COSMOS photometric catalog. We do this using the updated form of the energy balance code magphys (da Cunha et al. 2008, 2015), ⁹ which handles obscured galaxies with a wider range of star formation histories (SFHs) appropriate for high-redshift galaxies with significant dust content. This code uses the Bruzual & Charlot (2003) stellar population synthesis code to model spectral evolution, combined with the Charlot & Fall (2000) model to compute the total IR luminosity absorbed and reradiated by dust in the ISM. The high-redshift version of this code takes into account the fact that the copious amounts of dust in early star-forming galaxies are optically thick to the light in stellar birth clouds, such that the optical and IR parts of the SED are not produced in the same parts of a galaxy (i.e., the IR SED is not simply the difference between the attenuated and unattenuated optical spectra).

We use the data only out to 24 μm because the six central CO sources in Pointing 1 are detected as a single Herschel/SPIRE source and deblending it will introduce new uncertainties in a portion of the SED that should not directly affect the stellar mass estimate. To verify this, we also test the SED fits including ALMA 1.1 mm dust continuum and the less confused 100 μm and 160 μm data, which resulted in much worse fits and derived stellar masses < 10% different from the fits that use only the OIR data. This is likely due to the issue of spatial decoupling of dust and stellar light, and the attenuation is well constrained with the IR point at 24 μm. Figure 6 shows the individual OIR SEDs for the Pointing 1 sources, which all have reduced χ² ≲ 1. The corresponding analysis (spectrum, moment 0 map, and SED) is shown for 450.58 in Figure 7.

Zoom In Zoom Out Reset image size

4.4.1. AGN Content

4.4.2. Stellar Masses

Since dust-obscured systems suffer heavy uncertainties owing to degeneracies in reddening effects, geometry of star-forming regions, and reliability of OIR photometry, we verify the stellar masses using several methods.

As a first pass, we simply assume a constant mass-to-light ratio arising from the rest-frame near-infrared luminosity, using the methods outlined in Borys et al. (2003) and Hainline et al. (2011). These studies model the M/L_H and M/L_K ratios using a sample of z ∼ 2 DSFGs fitted with the Hyper-z SED code (Bolzonella et al. 2000), averaging between the most bursty and the most continuous SFHs. K-band light has the advantage of having low sensitivity to previous SFHs and dust obscuration, the main source of degeneracy for dust-obscured systems. H band similarly traces redder stellar populations and suffers minimally from extinction from dust, and it further minimizes the influence of thermally pulsating asymptotic giant branch stars (Hainline et al. 2011).

H and K bands are both safe options here, with the main sources of uncertainty coming from the SED itself, a possible contribution from an AGN due to the power-law emission of hot dust, and a lack of constraint on the dust attenuation in the IR regime. Hainline et al. (2011) note that uncertainties in the mass-to-light ratio due to SFH alone result in, at most, a factor of 2–3 difference in the stellar mass, but they cannot be constrained further without spatially resolved observations of the unobscured stellar populations (which will become possible with the James Webb Space Telescope).

A more constraining measurement of the stellar masses is the value computed from our own SED modeling using magphys, whose inferred mass-to-light ratios have agreed with synthetic SEDs from simulations (e.g., Hayward & Smith 2015). We do find that for some galaxies the measured stellar masses seem unphysically high given the galaxies' maximum possible ages, especially when compared with the observed stellar mass functions that are typically used for abundance matching (e.g., Behroozi et al. 2010). However, this may not be surprising, as other studies have found that dusty star-forming galaxies tend to dominate the high-mass end (log (M_⋆/M_⊙) > 10.3) of the stellar mass function at z > 2 (Martis et al. 2016).

We find that the stellar masses calculated from the rest-frame H-band luminosity are on the order of unity or slightly lower than the magphys estimates, while the K-band estimates are either at unity or factors of 3–4 higher. It is possible that there is a nonstellar contribution to the rest-frame K-band luminosity, but since we have checked for an excess in the mid-infrared photometry above the SED and did not find significant differences, this could be due to the two bands tracing stellar populations of slightly different ages. We use the magphys estimates as the reported stellar masses, with the uncertainties given as a range between the lowest-mass and highest-mass results of all three methods. The stellar masses are fairly high, on the order of ∼10¹¹ M_⊙, implying that these galaxies are mature and nearing their stellar mass limits given the age of the universe at z = 2.5. The H- and K-band photometry and the final stellar mass estimates are listed in Table 3. In Section 6 we discuss the molecular gas fractions and the implication that these galaxies may be nearly quenched.

Table 3. Near-IR Photometry and Derived Stellar Masses

Name	Ks	H	M⋆
	(μJy)	(μJy)	(1011 M⊙)
ID1	4.98 ± 0.18	7.06 ± 0.20	${2.1}_{-1.5}^{+2.8}$
ID2	1.85 ± 0.18	2.10 ± 0.20	${0.1}_{-0.0}^{+0.1}$
ID3	2.59 ± 0.19	3.99 ± 0.21	${1.0}_{-0.6}^{+0.8}$
ID4	6.69 ± 0.18	10.89 ± 0.20	${2.9}_{-1.7}^{+1.7}$
ID5	2.54 ± 0.17	4.89 ± 0.17	${1.9}_{-1.1}^{+0.8}$
ID6	3.43 ± 0.19	6.44 ± 0.21	${1.6}_{-0.8}^{+1.3}$
ID7	1.97 ± 0.18	3.22 ± 0.20	${1.6}_{-0.9}^{+0.9}$
450.58	4.00 ± 0.17	7.10 ± 0.16	${1.7}_{-1.1}^{+1.8}$

Download table as:  ASCIITypeset image

Table 4. Long-wavelength Photometry Used to Calculate Aggregate Star Formation Properties

Name	S100μm	S160μm	S250μm	S350μm	S450μm	S500μm	S850μm	log LIR	SFR	τdep
	(mJy)	(mJy)	(mJy)	(mJy)	(mJy)	(mJy)	(mJy)	(M⊙)	(M⊙ yr−1)	(Myr)
Core	4.1 ± 1.7	35.0 ± 4.8	62.2 ± 13.9	86.7 ± 19.9	47.9 ± 7.7	75.9 ± 18.5	12.2 ± 1.1	13.19 ± 0.05	${2300}_{-250}^{+280}$	${340}_{-43}^{+39}$
ID5	1.1 ± 0.2	2.6 ± 0.5	14.0 ± 4.3	14.6 ± 4.4	12.8 ± 4.1	13.9 ± 1.5	4.1 ± 0.8	12.43 ± 0.08	${400}_{-70}^{+80}$	${180}_{-36}^{+30}$
450.58	<1.5	<4	11.3 ± 2.2	15.57 ± 2.7	15.3 ± 4.1	14.5 ± 3.1	5.3 ± 0.8	12.40 ± 0.11	${370}_{-80}^{+110}$	${885}_{-210}^{+264}$

Download table as:  ASCIITypeset image

To calculate the SFRs and L_FIR, we model the far-IR SED as outlined in Casey (2012), using Herschel/PACS 100 μm and 160 μm and SPIRE 250 μm, 350 μm, and 500 μm photometry. Because of the poor angular resolution of SPIRE (FWHM 36'' at 500 μm), six of the sources in Pointing 1 are spatially confused as one source, with an additional independent detection corresponding to ID5. We do not deblend the confused source owing to deblending model degeneracies. Since we are interested in the aggregate quantities for the protocluster core, we can still evaluate the total depletion timescale using the total SFR. Thus, the SFR, far-infrared luminosity, and gas depletion timescales are calculated individually for ID5 and 450.58, and as aggregate quantities for the six core galaxies, as well as any galaxies that may lie at the same redshift but are not detected in CO. After extracting the Herschel photometry, we fit a modified blackbody with a fixed slope (β = 1.8) where only λ_peak is allowed to vary. This results in a total SFR for the six core galaxies of ${2300}_{-250}^{+280}$ M_⊙ yr⁻¹, ${400}_{-70}^{+80}$ M_⊙ yr⁻¹ in ID5, and ${370}_{-80}^{+110}$ M_⊙ yr⁻¹ in 450.58. The relevant photometry and IR-derived quantities are listed in Table 4.

We dig further into the suggestion in the literature by Wang et al. (2016) that the Pointing 1 core is a virialized cluster itself, possibly embedded within the larger protocluster. To do this, we present two contrasting hypotheses. Under "Hypothesis A," the galaxies in VLA Pointing 1 make up a virialized, nascent cluster core. In this scenario, the galaxies are likely occupying a spherical volume with the different redshifts traced to line-of-sight peculiar velocities. Under "Hypothesis B," the galaxies of Pointing 1 may be extended in a filament along the line of sight and not yet virialized in a core. In this section we explore the evidence for both possibilities using the reported extended X-ray detection, our analysis of the member galaxies' dynamical states, estimates of allowable dark matter halo masses, and the relative rarity of a halo of this mass.

A 3D representation of the z = 2.503 knot is shown in the fourth panel of Figure 9. Six sources in Pointing 1 span 15 '' (140 physical kpc, or 490 comoving kpc) on the sky (with ID5 45'' away) and span 2.494 < z < 2.511 in redshift. Under Hypothesis A, this Δz corresponds to a peculiar velocity range of Δv = 1400 km s⁻¹, assuming that the galaxies are spherically distributed in one halo. In this case, the gas would not be extended along the line of sight as suggested by Figure 9, making the spatial distribution of the total molecular gas significantly less complicated. We compared the galaxies' offsets in velocity from the center with what would be expected from a point-source halo or a Navarro–Frenk–White (NFW) halo with concentration c = 5 (Prada et al. 2012). The velocity distribution is not Gaussian, but the velocities are consistent with a ∼10¹³ M_⊙ halo. It is possible that we are sparsely sampling a population of galaxies that are not detected in this study, since the core contains only six galaxies, so a single bound halo is not immediately ruled out.

Under Hypothesis B, we assume that the redshifts instead correspond to a line-of-sight filament, which has a maximum end-to-end size of 28 comoving Mpc (8 proper Mpc). Figure 9, plotted as an extended filament, shows that the 3D distribution of the sources is patchy and nonspherical. As discussed in the source extraction procedure (Section 3.1), the peak flux centroids were straightforward to identify, but the edges of the molecular gas reservoirs are more ambiguous. The VLA data have a resolution of 25 (∼20 kpc), so we cannot meaningfully resolve mergers or outflows, but the diversity of line profiles and close separations imply that these galaxies are not simple rotational disks. While we do not directly claim any of the CO detections to be formal mergers (as their optical counterparts are well resolved), it is clear from the extended nature of the overall molecular gas distribution that there is some gas interaction between galaxies. Under this hypothesis, we do not consider the galaxies to be relaxed.

Wang et al. (2016) report a 0.5–2.0 keV detection coincident with the center of Pointing 1 in the stacked COSMOS Chandra/XMM maps, which is 4σ on 16'' scales and ∼3σ at 32''. They report a flux within a 32'' aperture of (6.9 ± 2.0) × 10⁻¹⁶ erg s⁻¹ cm⁻², which we use as reported in our following calculations. We evolve a 2 keV X-ray spectrum typical of galaxy groups to estimate the K-correction (Böhringer & Werner 2010) and calculate a total luminosity of L_{0.1−2.4 keV} = (8.9 ± 2.0) × 10⁴³ erg s⁻¹. Because it is more extended than a point source arising from a single galaxy and its luminosity is of order ∼10⁴⁴ erg s⁻¹, Wang et al. (2016) conclude that this source is evidence for a hot ICM, indicating that the core is already a virialized cluster. This would make it the highest-redshift galaxy cluster known to date, consistent with our Hypothesis A.

In the event that the core is not yet virialized (Hypothesis B), the X-ray source may have plausible alternate origins, e.g., it may arise from a recently extinguished radio galaxy somewhere within the line-of-sight filament. In this scenario, the extended X-ray emission would outlive the radio detection of the galaxy. We will explore a possible alternative explanation of a so-called inverse Compton (IC) ghost in greater detail.

Based on the higher-resolution Chandra detection, we can rule out a point source, so if real it is likely extended (albeit with low S/N), though there is the possibility that the X-ray detection is not a real source at all. For the sake of arguing each hypothesis, we assume for now that the X-ray emission is real, extended, and associated with the protocluster core, as chance alignments of X-ray and radio emission in COSMOS have a probability of < 10⁻⁴ (Jelić et al. 2012). Figure 10 shows the Chandra/XMM 2σ–3σ contours reported by Wang et al. (2016), as well as the 1.4 GHz and 3 GHz radio continuum emission in the central six galaxies in Pointing 1. There are two point sources detected at 1.4 GHz that are coincident with CO detections (ID1 and ID3) totaling S_1.4 ≈ 140 μ Jy, which we will return to in our discussion of Hypothesis B.

Zoom In Zoom Out Reset image size

To further explore the X-ray emission in the context of Hypothesis A, we use the correlation between L_X and cluster halo masses (Leauthaud et al. 2010) to estimate the minimum luminosity of the ICM at z = 2.5. As described later in Section 7.3, we calculate a maximum allowable dynamical mass range of the core of log M_h /M_⊙ = 12.9–14.3. In this mass range, we can calculate the expected luminosity of an ICM-dominated X-ray signal, using the L_X − M_h relationship from Leauthaud et al. (2010), which was calibrated on groups and clusters in COSMOS. We find an expected range of log(L_X) = 44.4–45.7 erg s⁻¹, which is a factor of 3–50 times brighter than the observed X-ray luminosity here. However, it is possible that the formation of the ICM at z = 2.5 is still in its nascent stages and would not be as bright as the intracluster gas in more mature, lower-redshift clusters, even if they have similar masses.

To explore Hypothesis B, we turn to Fabian et al. (2009), which discusses an extended X-ray detection coincident with the giant elliptical galaxy HDF130 at z = 1.99 (Casey et al. 2009). This X-ray source, with luminosity L_2−10keV ≈ 5.4 × 10⁴³ erg s⁻¹ (comparable to our source), shows evidence for the IC ghost effect. Radio-loud AGNs, when active, accelerate electrons to Lorentz factors γ > 10⁴, which in turn undergo IC scattering off the CMB. Lorentz factors of γ ∼ 10⁴ are required to generate the GHz-frequency synchrotron emission in the radio band, but electrons with γ ∼ 1000 can upscatter CMB photons to the X-ray bands. The Compton cooling time is 10 times longer at X-ray wavelengths than at radio wavelengths, so the "ghost" of the extended lobes of the AGN would be detectable in X-ray long after the extended radio signal has died out.

Electron acceleration ends when the AGN shuts off, so the higher-frequency radio emission arising from synchrotron radiation will die off first. This means that with a recently extinguished radio galaxy we would not necessarily detect a luminous (≥mJy flux densities) radio AGN counterpart to the X-ray emission at 1.4 GHz (see, e.g., Simionescu et al. 2016), but rather at MHz frequencies. ¹⁰ Indeed, the 1.4 GHz radio emission in HDF130 is present as point sources only at the 400 μJy level, more consistent with a low-luminosity AGN. The most obvious candidate analog of HDF130 within our Pointing 1 is ID3, which appears to be somewhat radio luminous given its submillimeter emission; thus, its 140 μJy flux at 1.4 GHz could also be attributable to a low-luminosity AGN.

In terms of morphology, the double-lobed X-ray source in Fabian et al. (2009) is much more linearly extended (∼40'') than the single source considered here. On the most compact scales (10''), the morphology of the Wang et al. (2016) source appears mostly spherical or slightly elongated in the NW direction and does not show evidence for lobed structure. However, the detection beyond 20'' is marginal (2σ–3.3σ). In both models (e.g., Mocz et al. 2011) and observations (e.g., Jelić et al. 2012), the X-ray IC lobes arising from powerful radio galaxies are expected to span a few hundred up to 1000 kpc, or ∼25''–120'' on the sky. The source here spans only 20'' on the sky at maximum (170 proper kpc), on the smaller end of the range of IC ghost lobe sizes (see 690 kpc for HDF130 in Fabian et al. 2009).

The measured 20'' size of this source also happens to be a factor of 2–3 smaller than the typical ICM diameter of 600–700 kpc in low-redshift clusters (Piffaretti et al. 2011) of similar mass (see halo mass calculations in Section 7.3), though it is unclear whether the ICM would be physically smaller in a structure that has very recently collapsed.

Finally, the relatively low S/N of the X-ray signal at large radii causes some additional concern. For example, the Chandra/XMM flux drops off quickly—<50% of its peak value at double the extraction radius—which is not typical of the radial profiles of intracluster gas. However, the overall low S/N of the extended Chandra/XMM detection (<3σ beyond 30'' and 4σ at 16'') makes the radial profile difficult to constrain and the extended nature marginal. Though the XMM-only map has better sensitivity, the source is several arcminutes off-axis in the COSMOS maps, so it is a marginal point-source detection as well. Based on its size and luminosity, we argue that the X-ray signal, if real, is more consistent with emission from a single galaxy IC ghost than with an ICM origin. In Section 7.4, we discuss the probabilities of observing either scenario (nonthermal IC emission or ICM gas) given the source's selection in the COSMOS survey field.

Estimating the halo mass for this source necessarily depends on whether these galaxies occupy a single dark matter halo or each reside in individual halos at z = 2.5, so we use several methods to calculate halo masses using different assumptions about the available CO, X-ray, and optical data. For each of the halo mass estimators we discuss, we also calculate the equivalent z = 0 halo masses based on the halo assembly histories used by Harrison & Hotchkiss (2013) in their publicly available code to calculate relative halo rareness. These authors describe how to compute curves in mass–redshift space along which the numbers of collapsed halos at that mass and redshift are equally rare according to the selected rareness statistic. Thus, the equivalent z = 0 halo mass is that which has the same value of the rareness statistic as the z = 2.5 halo. We caution that the equivalent z = 0 halo mass is not a prediction for the descendant mass of the halo—there is significant scatter in the final halo masses of nascent protocluster halos (e.g., Chiang et al. 2013; Lovell et al. 2018), so we cannot constrain this further.

The following calculations assume that Hypothesis A holds, in which the extended morphology of the molecular gas along the line of sight would be attributable to velocity dispersion of sources within a single gravitationally bound halo. We can place a hard lower limit on the single halo mass by assuming that the baryonic concentration within the halo is not below the cosmic baryonic fraction, f_b = 0.19 (Planck Collaboration et al. 2016; Behroozi & Silk 2018), reaching a conservative minimum of log (M_h /M_⊙) ≈ 12.9. A more useful lower limit for Hypothesis A, though, is the minimum halo mass for which the CO galaxies could be bound to a single halo. As in Section 5, we compared the velocity offsets of each galaxy with those expected from an NFW halo. The minimum halo mass at which these galaxies can be considered bound is log (M_h /M_⊙) ≈ 13.2. Although it is a small sample size, the core galaxies' velocities are consistent with a bound dark matter halo above this mass. This minimum halo mass at z = 2.5 has an equivalent z = 0 halo mass of approximately 8 × 10¹³ M_⊙.

Also under Hypothesis A, we estimate the dynamical mass (M_dyn = M_halo + M_baryon) using the virial theorem:

$$\begin{eqnarray}&&{M}_{\mathrm{dyn}}=C\frac{R{\sigma }^{2}}{G},\end{eqnarray} \tag{ 3 }$$

where R is the proper size of the structure, σ is the velocity dispersion with respect to the center, and C is a constant taken to be 3 for spherical mass distributions (Ho & Kormendy 2000). We take R to be 190 kpc, which is the transverse size of the Pointing 1 core (and consistent with the X-ray size of 130 kpc reported in Wang et al. 2016). The velocity dispersion of the CO detections, σ, with respect to the average redshift of the seven sources (z = 2.503) is measured to be σ = 500 ± 120 km s⁻¹. We estimate M_h(z=2.5) = (3.3 ± 0.6) × 10¹³ M_⊙, which translates to a z = 0 halo mass of approximately 5 × 10¹⁴ M_⊙, comparable to the Coma Cluster. This is consistent with the dynamical mass estimate of ${3.2}_{-1.2}^{+1.9}\times {10}^{13}$ M_⊙ from Wang et al. (2016), which was calculated using the transverse extent of the X-ray detection.

Still assuming Hypothesis A, we also calculate the halo mass directly using the luminosity of the 0.5–2 keV X-ray detection assuming that it arises from a hot ICM in a cluster-sized halo. We use the K-correction described in Section 7.2 to derive the rest-frame luminosity at 0.1–2.4 keV and then use the L_X/M_h relation from Leauthaud et al. (2010). This results in a halo mass estimate of M_h(z=2.5) = (5.0 ± 0.2 ) × 10¹³ M_⊙ (or M_h(z=0) ∼ 9 × 10¹⁴ M_⊙), which is the highest mass estimate of any of our methods. One should note, however, that this correlation with halo mass has been measured only out to z = 1 for virialized cluster cores. In the event that there is a redshift evolution in this correlation, this halo mass would be inaccurate.

Alternatively, under Hypothesis B, the redshifts of the CO sources may instead correspond to true separations along the line of sight, and we assume they occupy completely separate dark matter halos. Under this hypothesis, we calculate the total halo mass by summing the individual halo masses of the seven galaxies. To measure the individual halo masses, we use abundance matching, using the stellar–halo mass relationships of Behroozi et al. (2010). Since our measurements of the stellar masses approach the upper limit of the predicted z ∼ 2 stellar masses according to Behroozi et al. (2010), this sometimes results in invalid extrapolations to compute halo masses. To avoid this problem, we use a Monte Carlo approach assuming a lognormal distribution of stellar masses for each galaxy, centered on the reported stellar masses in Table 3. Taken in aggregate, this results in a total halo mass of M_h(z=2.5) = (1.7 ± 0.8) × 10¹³ M_⊙. This should be taken as an upper limit for Hypothesis B because it represents the extreme scenario in which there is no gravitational interaction between any of the galaxies; here we cannot account for the possible double-counting of overlapping dark matter halos. If we assume that these galaxies—currently extended in a filament—will eventually collapse into a single halo, its equivalent halo mass at z = 0 is 2 × 10¹⁴ M_⊙. While this is the smallest mass estimate of our three methods by a factor of a few, it still implies a Coma-sized cluster at z = 0.

Finally, the absolute maximum allowed halo mass is found by calculating the exclusion curve as a function of redshift, i.e., the halo mass above which it would come into tension with our ΛCDM cosmological parameters. Using the halo mass assembly history from Harrison & Hotchkiss (2013), we find that there is a <5% chance of detecting a halo more massive than 2 × 10¹⁴ M_⊙ within the COSMOS survey at z = 2.5, so we can safely assume that this protocluster core is not more massive than this.

In Figure 11 we show a summary of each of our baryonic and halo mass estimates. The allowable halo mass ranges from log M_h /M_⊙ = 12.9 to 14.3, where the minimum represents the lowest baryonic fraction and the maximum is that allowed by our cosmological assumptions. Within the errors, all of our halo mass estimates are consistent with one another: under either hypothesis, this is a uniquely massive overdensity that will rival the largest galaxy clusters by z = 0. It is likely that the halo structure lies somewhere between the two extremes of a single coalesced halo and seven completely separate halos. Because the CO-detected galaxies are close enough (in the transverse direction, regardless of peculiar velocities) that there is probably overlap in their dark matter halos, it is likely, no matter the configuration at z = 2.5, that the seven galaxies in Pointing 1 will collapse into a single halo at lower redshift. Indeed, it is plausible that they will form a single brightest cluster galaxy within the z = 2.47 parent cluster by z = 0.

Zoom In Zoom Out Reset image size

In order to place the relative rarity of the z = 2.5 core structure in context, we compare the two formation hypotheses from a probabilistic standpoint. To do this, we calculate the probabilities that structures like the core discussed in this paper would be reasonably observable in a survey size of 2 deg², here corresponding to the COSMOS field in which the source was found. First, we estimate how many total halos we expect to exist in the COSMOS field with mass similar to or greater than the z = 2.5 core in this study that could be detectable in X-ray emission. We compare the probabilities of the X-ray emission in each halo originating from hot ICM gas, nonthermal emission like IC processes, or an AGN with radio and X-ray lobes. Independently, we then estimate the number of halos in COSMOS massive enough to host massive radio galaxies, a population that likely overlaps with the halos that could emit X-rays. Finally, we estimate simply the rareness of a spatial overdensity of massive star-forming galaxies (which may not necessarily reside in a single halo or be associated with AGN or ICM features).

7.4.1. Rarity of X-Ray-emitting Halos

We estimate the total number of halos that are as rare as the core discussed here and could be expected to be observable. Following the methods used in Marrone et al. (2018) for SPT0311-58, we use the halo rarity framework of Harrison & Hotchkiss (2013) as discussed in Section 7.3 to estimate the relative rareness of this protocluster core. By imposing a physically sensible redshift range and observational limits of a survey, one can then use this statistic to infer how many halos in that survey volume might collapse into cluster-sized halos ¹¹ of similar or greater mass by z = 0.

The necessary observational constraints on the rareness statistic are the survey area, the redshift range in which to search for objects, and the selection function or the limiting halo mass as a function of redshift.

We describe the observational selection function for such rare halos using X-ray emission; an implicit assumption of this calculation is that all similarly rare halos will be detectable in X-ray emission above a given mass threshold as a function of redshift. Whether or not this is a realistic expectation is discussed further below. The survey area is the 2 deg² covered by both XMM and Chandra in COSMOS. The volume over which we calculate the probability of observing a halo of any mass and redshift that is equally as rare as this protocluster core is enclosed by 1.5 < z < 4, motivated by the exclusion of mature clusters at z < 1.5 and the assumption that a virialized, detectable ICM would not physically exist before z = 4 (Chiang et al. 2017), nor would enough radio galaxies have turned off to produce IC ghosts in X-rays (Mocz et al. 2011). Finally, we define a limiting halo mass as a function of redshift. We assign this based on the X-ray flux limits in COSMOS such that we are selecting for X-ray detectability. In other words, we assume that the X-ray luminosity scales with halo mass, regardless of whether the halo is actually producing X-rays. Then, using the redshift-dependent K-correction function noted in Section 7.3, we estimate the expected observed luminosity of an X-ray-emitting halo as a function of redshift. This is then matched to a halo mass using the L_X–M_h relation of Leauthaud et al. (2010), resulting in a limiting halo mass function set by the Chandra flux limit in the COSMOS survey (which is 2.2 × 10⁻¹⁶ erg s⁻¹ cm⁻²).

Based on the survey limits alone, we estimate ∼10 possible halos similar to this protocluster core that are hypothetically observable via extended X-ray emission within in the COSMOS field between 1.5 < z < 4. This implies that, based on its halo mass alone, this source is not unphysically rare, and rather it might be expected given the relatively large size of the COSMOS field.

One significant caveat of this calculation, however, is that not all massive halos at high redshift will actually generate a detectable X-ray signal, since not all massive halos will host AGNs with lobes, and/or protoclusters may not have accreted sufficient gas to generate thermal bremsstrahlung emission from a hot ICM. We stress that this result of 10 similar halos is the maximum number of halos that might be massive enough to emit detectable X-rays—it is based simply on the X-ray flux limit of COSMOS and is not a predicted luminosity function of cluster ICM. Indeed, observations suggest that cluster halos do not undergo shock heating to virial temperatures until z ∼ 2 (see compilation of Overzier 2016, Figure 11); the implications of that work are that the number of detectable halos of similar rarity in COSMOS would be reduced to a total of 4 and that the z = 2.5 structure would be categorically outside of the allowable redshift range. Additionally, studies evolving the X-ray and radio luminosity functions to high redshift (Celotti & Fabian 2004; Finoguenov et al. 2010; Mocz et al. 2011) suggest that clusters at high redshift (z > 2) and luminosities greater than a few times 10⁴³ erg s⁻¹ should be exceedingly rare, while, in contrast, extended X-ray emission arising from Fanaroff–Riley Type II (FR II; Fanaroff & Riley 1974) galaxies will dominate the extended X-ray source luminosity function. The cumulative upper limit number of halos as a function of redshift is plotted in black in Figure 12.

Zoom In Zoom Out Reset image size

7.4.2. Rarity of Massive Radio Galaxy Halos

7.4.3. Probability of Detecting IC Ghosts

7.4.4. ICM or IC Ghost: The More Likely Physical Scenario?

Since the X-ray detection in question here is itself marginal, it is not this X-ray source in isolation that makes this proposed protocluster core so extraordinary. Instead, it is the X-ray signal in addition to its extreme star formation activity in the form of multiple rare star-forming galaxies with demonstrated fast duty cycles (with a total depletion timescale of ∼350 Myr; Section 5). We thus use other elements of the observational data to determine how common star-forming galaxy overdensities like in this protocluster core might be in the COSMOS field using the Laigle et al. (2016) photometric redshift catalog.

We search for the observed number of overdensities in the catalog using a near-IR selection criterion, which affords us the angular resolution that selecting by submillimeter wavelengths (e.g., SPIRE or Scuba-2) would lack and allows us to select for red star-forming galaxies. First, we imposed a photometric redshift cut of z_phot < 4, so that we are selecting for overdensities that may coincide with the X-ray-detectable halos discussed in the previous section. We also imposed a magnitude range of $21.5\lt {m}_{{K}_{s}}\lt 23.4$, which spans from the 3σ detection limit to the magnitude at which most galaxies at z > 1 do not exceed (Davidzon et al. 2017), since we also want to filter out low-redshift overdensities.

We define the surface density on the sky Σ_N by a nearest neighbor algorithm:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{N}=\frac{N}{\pi {r}_{N}^{2}},\end{eqnarray} \tag{ 4 }$$

where N is the Nth neighbor and r_N is the distance (in arcseconds) to the Nth-nearest neighbor. After determining the average 2D density of K_s -band sources in COSMOS, we decided that N = 10 is a sufficient number of neighbors to determine a spatial overdensity above the field. We use a k-d tree to construct the catalog of 10th-nearest neighbors, first in two dimensions (only sky positions) and then in three by including their photometric redshifts (spectroscopic when available). A group of galaxies is then declared a 2D overdensity if its 10th-nearest neighbor is closer than 50'', which is a distance significance of 2σ above the average r₁₀ in the catalog. A 3D overdensity is declared if σ_z ≤ 0.2 for the 10 sources in question, comparable to the uncertainties in photometric redshift overall.

Using this simple K_s -band magnitude cut, we predict 204 chance alignments along the line of sight that look like overdensities on the sky. Because there are so many galaxies in this catalog (35,000), true overdensity signals are likely washed out by outliers, because even the protocluster discussed here is not detected as a 2D overdensity. Therefore, we impose an additional color cut to select for massive red star-forming galaxies whose neighbors have similar colors, where J − K_s ≥ 1.3 (consistent with the criteria outlined in Wang et al. 2016).

The final catalog of 6000 galaxies contains 16 2D associations of 10 neighboring galaxies and 2 3D associations across the whole COSMOS field, both including the core in question. We then searched the Herschel SPIRE maps for bright (>4σ), red (S₅₀₀ > S₂₅₀) counterparts to these associations, as this has frequently been used to select for obscured star formation in galaxy overdensities (e.g., Clements et al. 2014; Greenslade et al. 2018). Only one candidate overdensity passes this criterion, which is the core discussed here.

Finally, we calculate the "probability to exceed" measure from Harrison & Hotchkiss (2013), finding that there is only a 30% chance that this is the most massive observable structure at 2 < z < 3 in COSMOS. Since no clear candidates emerge from our selection criteria, they could be found through other selection methods such as pairwise angular correlation functions (Ando et al. 2020), Lyα tomography (Li et al. 2021; Horowitz et al. 2019) and absorption surveys (Cai et al. 2017).

We summarize the probabilities of observing the X-ray detections discussed in Section 7.4 and the K_s -band galaxy overdensities in Figure 12. We expect that approximately 10 halos of similar z = 0 halo masses to this protocluster core could exist in the COSMOS field if they are detectable via extended X-ray emission. Using the COSMOS near-IR catalog, we showed that a comparable number of 2D overdensities within 50'' could exist in the same field, although many of these will be chance alignments along the line of sight. Both forms of X-ray emission may be more common than the 3D overdensities of red massive galaxies, as we find only two 3D overdensities across the full field. We find that only one of these candidates—this structure—is also associated with an unusually bright red Herschel detection associated with such prodigious star formation (2290 M_⊙ yr⁻¹), making this structure one of the most unique overdense objects in COSMOS below z = 4.

We extract spectra in two ways from each VLA 33 GHz cube: (1) a single-pixel extraction at the position of the peak flux in moment 0 maps, with the channels centered on the expected redshift of the source, and (2) extraction at the same position but with a custom aperture, where all adjacent pixels have S/N > 3 (typically the size of one to two VLA beams), taking care to place the apertures such that we do not double-count flux between nearby sources.

If the sources were resolved, we would expect there to be a significant excess in the integrated flux density measured in a large aperture versus the brightest pixel, whereas an unresolved source would show little variation in the integrated flux density regardless of the aperture size (provided that there is not a significant velocity gradient across the source).

In the top panel of Figure 13 we show an example spectrum of the brightest source (ID1) extracted with each of these two methods and show that the line profiles do not differ between an aperture and the brightest pixel, which we expect since the sources have low enough S/N that we would not be sensitive to velocity gradients across the source. In the bottom panel of Figure 13 we show that most sources are unresolved or marginally resolved by the ratio of the line intensity I_CO extracted in a custom aperture to that measured from the single brightest pixel. Most sources are roughly consistent with unity, indicating they are point-like sources, but several of the brighter sources (ID1, ID4, and 450.58) appear marginally resolved, where we measure 25%–75% more flux within a beam-sized aperture.

Zoom In Zoom Out Reset image size

Another way to test whether sources are resolved would be to fit a 2D Gaussian model to each source and compare the estimated size with the size of the VLA beam. We tried this using the CASA task imfit, but we encounter catastrophic failures or unphysically large fit sizes owing to the crowding of the sources and likely decorrelation of the signal at high observing frequencies. We additionally reimaged the cubes using uvtaper to degrade the spatial resolution and test whether we have resolved out any flux. Regardless of the spatial resolution, we find total aperture fluxes consistent within 10% and a consistent inability to fit individual 2D Gaussian models to the sources. Performing the same analysis in the u-v plane or on cleaned data does not change the resolution results.

We calculate the CO line properties by directly measuring the spectral line moments, following the procedures of Bothwell et al. (2013). The CO redshift is determined by the first moment, which is the intensity-weighted line center. The intensity-weighted second moment determines the line width, which is defined as

$$\begin{eqnarray}&&{s}_{\nu }=\frac{\int {(\nu -\bar{\nu })}^{2}{I}_{\nu }\,d\nu }{\int {I}_{\nu }\,d\nu },\end{eqnarray} \tag{ A1 }$$

such that the Gaussian FWHM is $2\sqrt{2\,\mathrm{ln}2}$ s_ν .

The calculation of spectral moments is sensitive to the noise in the spectrum, so we run a Monte Carlo simulation 2000 times, in which Gaussian noise (σ equal to the rms of the line-free spectrum) is artificially added to the spectral data and the line moments recalculated. (Artificially adding noise results in a reported error that is a factor of 2^−1/2 higher than the true noise, but this is accounted for in the final results.) In our Monte Carlo procedure, we restrict the line search to a radius of 100 km s⁻¹ from the prior guess for the redshifts. The peak, center, and FWHM are determined by the medians of the results of this procedure, with the errors as the standard deviations in these distributions. We additionally multiply the spectrum by −1 and repeat the procedure to make sure there are no negative 3σ "sources" in the spectrum, of which we find none.

Finally, the integrated line intensity I_CO is formally the velocity-integrated line flux along ± 2σ from the line center, but we simply sum the flux within the range where the signal remains positive in adjacent channels.

We use the CO (1–0) line measurements to calculate molecular gas mass, using the usual factor α_CO to convert from line luminosity to total molecular gas mass, but the appropriate value of α_CO to use can be ambiguous. The Galactic value of α_CO = 4 M_⊙/(K km s⁻¹ pc²) was calibrated using observations of individually resolved giant molecular clouds (GMCs) within the Milky Way and nearby spiral galaxies, while local luminous infrared galaxies (LIRGS) were presumed to have much lower conversion factors of α_CO ∼0.8 M_⊙/(K km s⁻¹ pc²), with CO emission stemming from a warm, diffuse ISM rather than virialized GMCs (see review of Bolatto et al. 2013). While α_CO is likely dependent on metallicity and the thermal state of molecular gas clouds, a value of α_CO = 1 M_⊙/(K km s⁻¹ pc²) is typically used for unresolved, compact submillimeter galaxies (Bolatto et al. 2013; Carilli & Walter 2013) at high redshift. If galaxies are dominated by dense, nuclear starbursts, then a low value of α_CO is appropriate, which has been a valid assumption especially in the case where higher α_CO values yield gas masses in excess of the dynamical masses of galaxies.

However, Scoville et al. (2016) broadly argue for a higher value of α_CO in DSFGs, especially in cases where the densest sites of star formation—where there is significant dust heating—are not spatially resolved. Instead, they advocate for using a conversion factor that is globally averaged, where galaxies are not dominated by these hot starbursting regions. In this case the authors advocate for the Galactic value of α_CO, which is α_CO = 4.5 M_⊙/(K km s⁻¹ pc²) for the mass of H₂, or α_CO = 6.5 M_⊙/(K km s⁻¹ pc²) if assuming a factor of 1.36 for the contribution of helium. We adopt this Galactic value of α_CO = 6.5 M_⊙/(K km s⁻¹ pc²) here to indicate the total expected molecular gas mass.