Physical Characterization of an Unlensed, Dusty Star-forming Galaxy at z = 5.85

The galaxy MM J100026.36+021527.9 first appeared in the literature in Bertoldi et al. (2007) as "ID9" detected by the MAMBO instrument at the Institut de Radioastronomie Millimétrique (IRAM) 30 m telescope at a wavelength of 1.2 mm with S_1.2 = 4.9 ± 0.9 mJy. We adopt the shorthand name Mambo-9 throughout this paper. Plateau de Bure Interferometer imaging of Mambo-9 exists at 1 mm with 4σ significance (its analysis was included in the Ph.D. thesis of Manuel Aravena, 2009²⁴ ). The redshift was not known at the time. The detection was independently corroborated by Aretxaga et al. (2011) as "AzTEC/C148" using the AzTEC instrument on the the Atacama Submillimeter Telescope Experiment (ASTE) at 1.1 mm with S_1.1 = 4.6 ± 1.2 mJy, in agreement with the earlier MAMBO measurement. The source was then later identified in the SCUBA-2 850 μm map of Casey et al. (2013) as "850.43" (with S₈₅₀ = 5.55 ± 1.11 mJy and no corresponding 450 μm counterpart) and further as "COS.0059" in Geach et al. (2017) with S₈₅₀ = 5.84 ± 0.87 mJy.

Mambo-9 has no clear counterpart from either Spitzer or Herschel in the range 24–500 μm (Le Floc'h et al. 2009; Lutz et al. 2011; Oliver et al. 2012). The lack of detection in these bands implies that the spectral energy distribution (SED) traces unusually cold dust or, alternatively, a very high redshift solution. This prompted a number of teams to pursue ALMA follow-up observations of the source, including the 3 mm spectral scan presented by Jin et al. (2019).

Our interest in Mambo-9 stems from the new ALMA 2 mm blank field in the COSMOS field (Cycle 6 program 2018.1.00231.S, PI Casey). The scientific objective of the 2 mm blank-field map is to constrain the volume density of DSFGs at z ≳ 4. This is made possible because 2 mm detection is an effective way to "filter out" lower-redshift DSFGs at 1 ≲ z ≲ 3, as detailed in the modeling work of Casey et al. (2018a, 2018b). Blank-field maps at shorter wavelengths (e.g., 870 μm and 1.1 mm) identify more sources than at 2 mm per given solid angle, but such work then suffers from the "needle in the haystack" problem of identifying which sources sit at z ≳ 4 (e.g., as described in Casey et al. 2019). Analysis of the 2 mm blank-field map data set will follow in a later paper.

Mambo-9 was the brightest source identified in the first 9.4 arcmin² of delivered 2 mm map data, and the ratio of 850 μm flux density to 2 mm flux density (${S}_{850\mu {\rm{m}}}/{S}_{2\mathrm{mm}}=8.3\pm 0.9$) implied a high-redshift solution, where a higher value (${S}_{850\mu {\rm{m}}}/{S}_{2\mathrm{mm}}=15\pm 3$) would be expected for DSFGs at z ∼ 1−4. An independent analysis of data from Cycle 5 program 2017.1.00373.S (PI Jin) identified a 4σ emission line consistent with the measured redshift solution as well as other possible high-z solutions (and corresponding candidate ∼4σ emission peaks), which led to the proposal for the data described herein. This line and spectroscopic redshift have since been reported in Jin et al. (2019) based on the identification of the tentative 101 GHz line as ¹²CO(J = 6→5) and a line at the edge of ALMA band 3, ∼84 GHz, as ¹²CO(J = 5→4). Our independent analysis of the same data did not lead to a significant detection of the ¹²CO(J = 5→4) line. Because the 84 GHz line sits at the very edge of ALMA band 3 (whose lower limit frequency is 84 GHz) and at low signal-to-noise ratio (S/N), additional tunings were needed to elucidate the redshift solution and characterize Mambo-9.

Observations with ALMA were obtained under program 2018.1.00037.A split into three scheduling blocks tuned to three different frequencies: two in band 3 (3 mm) and one in band 7 (870 μm). The frequencies were chosen to secure the redshift of Mambo-9, which was determined to have multiple viable solutions²⁵ at 4 ≲ z ≲ 9. The ALMA data were reduced and imaged using the Common Astronomy Software Application²⁶ (CASA) version 5.4.0 following the standard reduction pipeline scripts and using manually defined clean boxes during the cleaning process. For band 7 observations, Mambo-9 is detected at very high S/N (111σ) such that we performed self-calibration. Also in band 7, a few noisy channels in one spectral window were identified and flagged in the frequency range 331.35–333.33 GHz after visual inspection of the calibrated data. No additional flagging was required for band 7 or band 3 observations.

Band 7 observations covered frequencies 329.5–333.5 GHz and 341.5–345.5 GHz. They were taken on 2019 May 2 in the C43-4 configuration, with a synthesized beam of 036 × 030 (using natural weighting), a total integration time of 6383 s, and a mean precipitable water vapor (PWV) of 0.9 mm. The continuum rms reached over the 7.5 GHz bandwidth is 26.9 μJy/beam. We explored different visibility weights for imaging, using both natural and Briggs weighting (robust = 0.0–0.5), the latter to spatially resolve the distribution of dust in the primary components of Mambo-9.

Band 3 data were obtained in two tunings. The first covers frequencies 86.6–90.3 GHz and 98.6–102.4 GHz. These data were taken on 2019 April 30 and 2019 May 1 in C43-4 with a total integration time of 14,668 s, resolution of 119 × 109 (using natural weighting), and an average PWV ranging from 1.3 to 2.3 mm. The rms reached over a 50 km s⁻¹ channel width was 0.124 mJy/beam. The second band 3 tuning covers frequencies 94.8–98.5 GHz and 106.6–110.4 GHz. These data were taken on 2019 April 30, 2019 May 1, and 2019 May 2, with a spatial resolution of 114 × 093 (using natural weighting), a total integration time of 15,010 s, and a mean PWV ranging from 0.9 to 2.0 mm. The rms reached for a 50 km s⁻¹ channel width was 0.088 mJy/beam. All band 3 data presented in this paper also are coadded in the visibility plane with the archival spectral scan from Jin et al. (2019; 2017.1.00373.S), which contributes a total of ∼1500 s of integration time to the total (9%). Our 3 mm continuum flux density is consistent with the measurement from the Jin et al. data. In an effort to spatially resolve the components of Mambo-9, we explore different weightings with different synthesized beam sizes. The overall continuum rms achieved in the coaddition of all of the band 3 data is 5.2 μJy using natural weighting to maximize the line S/N.

Band 6 continuum data also exist for Mambo-9 from the 2016.1.00279.S program (PI Oteo), which achieved a continuum sensitivity of 0.16 mJy/beam at a representative frequency of 233 GHz; the synthesized beam size in the band 6 data is 081 × 068 (with Briggs weighting and robust = 0.0). Band 4 continuum data, from our separate 2 mm blank-field map program, 2018.1.00231.S (PI Casey), has a continuum sensitivity of 0.11 mJy/beam, a representative frequency of 147 GHz, and a synthesized beam of 183 × 143 (natural weighting). More analysis of the band 4 data will be presented in a forthcoming work. We use Briggs weighting with robust = 0.0 where the S/N is sufficiently high to provide improved spatial resolution, while we use natural weighting to measure sources' integrated flux densities, especially when the detection S/N is near the 5σ threshold.

Analysis of the band 7 data reveals two distinct point sources separated by ∼1'' oriented in a north–south direction, as shown in Figures 1 and 2. We call the northern, brighter source component A (or Mambo-9-A) and the southern, fainter source component B (or Mambo-9-B).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Mambo-9 sits in the central portion of the Cosmic Evolution Survey Field (COSMOS) covered by the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011) and thus benefits from some of the deepest ancillary data available. This source has no counterpart in the deep imaging catalog of Laigle et al. (2016). The source is detected in the deep S-CANDELS Spitzer IRAC data (Ashby et al. 2015) and appears to be marginally resolved in the north–south direction, consistent with the positions and orientation of components A and B with respect to one another. There is a marginal detection (∼3σ) of a portion of the source in the deep Hubble Space Telescope (HST) F160W imaging data near component A. However, using a 06 extraction aperture centered on the ALMA 870 μm dust map reduces this potential marginal detection to <1σ significance. There is also a detection of a faint source 1'' to the south of component B in UltraVISTA Ks band, Hubble F125W and F160W imaging, though we believe it is unassociated with Mambo-9 based on the different optical and infrared colors (e.g., the Ks-band magnitude K_s = 26.36 ± 0.35 yet there is no associated IRAC emission) and lack of an ALMA counterpart in the extraordinarily deep 870 μm image. Mambo-9 is not detected in the Spitzer 24 μm imaging (Le Floc'h et al. 2005), nor Herschel PACS 100 μm or 160 μm (Elbaz et al. 2011), or SPIRE 250 μm, 350 μm, or 500 μm (Oliver et al. 2012), which would be expected for sources of similar 850 μm flux densities at z ≲ 2–3 (Casey 2012; Casey et al. 2012; Gruppioni et al. 2013). Note that Jin et al. (2019) report photometry for this source in the Herschel SPIRE bands using the "super-deblended" extraction technique (Jin et al. 2018), although examination of the Herschel map shows no detection or contamination by nearby neighbors; we adopt upper limits for the SPIRE bands instead. Our upper limits come from the confusion-noise rms, which dominates the uncertainty of flux calibration at low S/N (Nguyen et al. 2010, upper limits in the Herschel PACS bands are limited by instrumental noise, Lutz et al. 2011). Mambo-9 is not detected in the deep 1.4 GHz radio imaging of Schinnerer et al. (2007), and though not formally detected above the 5σ threshold in the 3 GHz Very Large Array (VLA) map, there is a 3.2σ-significance peak near Mambo-9-A at 3 GHz (Smolčić et al. 2017). There is no X-ray detection. Across all measured data sets, Mambo-9 is only detected above >3σ significance in seven of 30+ different wavebands. The constraining photometric data are presented in Table 1.

We conclude that Mambo-9 is unlikely to be strongly gravitationally lensed. This is due to the lack of foreground galaxies detected at other wavelengths in the optical and near infrared. The nearest possible lensing galaxy is offset 41 to the north and has a photometric redshift of z = 2.4 and estimated stellar mass of 4 × 10⁹ M_⊙. Assuming a halo mass of 4 × 10¹¹M_⊙, this would then lead to a maximum value of lensing magnification of μ = 1.15 − 1.3 based on conservative assumptions as to the lensing Einstein radius. Strong gravitational lensing by foreground galaxies does affect several other well-known high-z DSFGs, including the three DSFGs known to sit at higher redshifts than Mambo-9: G09 83808 at z = 6.03 (Zavala et al. 2018), HFLS3 at z = 6.34 (Cooray et al. 2014), and SPT0311 at z = 6.9 (Strandet et al. 2017). All of these systems have bright optical/near-infrared sources within 2'' of the millimeter source center. This leads us to conclude that Mambo-9 is the highest-redshift unlensed DSFG identified to date.

Figure 3 shows the full band 3 data set for Mambo-9 in context; the gray background spectrum is the spectral scan from Jin et al. (2019), with reported detection of lines at 101 GHz and 84.2 GHz corresponding to ¹²CO(J = 5→4) and ¹²CO(J = 6→5). Our independent analysis of the Jin et al. data set, before publication of their reported lines, did not lead to significant detection of the 84.2 GHz line. Thus, our band 3 data (shown in yellow) were tuned to frequencies that would rule in or out other possible solutions in the range 4 < z < 9. The detection of emission features at 101 and 109.8 GHz independently corroborates the reported redshift solution in Jin et al. (2019). Note that our derived redshifts for components A (z = 5.850) and B (z = 5.852) differ slightly from the reported redshift in Jin et al. of z = 5.847. We attribute this to the difference in S/N on the ¹²CO(J = 6→5) feature.

The aggregate band 3 continuum has a flux density of 131.4 ± 5.9 μJy (20.7σ significance) in the full bandwidth and in many individual channels of our data set, so analysis of molecular line emission requires continuum-subtracted data. Note that in Table 1, the band 3 continuum flux density is split into two independent measurements by frequency, given the high S/N on the aggregate data set.

The integrated ¹²CO(J = 6→5) line flux is 0.48 ± 0.03 Jy km s⁻¹ (14.7σ significance), and the p-H₂O (2_1,1 → 2_0,2) line flux is 0.09 ± 0.02 Jy km s⁻¹ (4.3σ significance). Components A and B are only distinguishable in band 3 data when using Briggs weighting (robust = 0.0), because natural weighting results in a synthesized beam slightly larger than the 1'' separation between the sources; the disadvantage of Briggs weighting is the potential for resolving out emission. Using Briggs weighting, component A has an integrated ¹²CO(J = 6→5) line flux of 0.43 ± 0.03 Jy km s⁻¹ (13.3σ significance), while component B has a line flux of 0.07 ± 0.02 Jy km s⁻¹ (3.8σ significance). The ¹²CO(J = 6→5) and 870 μm continua are spatially aligned for component B (see Figure 1). The p-H₂O (2_1,1 → 2_0,2) detection only corresponds to component A. Within measurement uncertainties, we find that the sum of band 3 measurements of component A and component B separately (using Briggs weighting and robust = 0.0) is in agreement with the total integrated quantities as measured with natural weighting.

Figure 4 shows the line spectra for both ¹²CO(J = 6→5) and p-H₂O (2_1,1 → 2_0,2), with ¹²CO(J = 6→5) broken down into the two components. While the total line emission appears roughly Gaussian, the spectrum of component A alone appears double peaked and possibly indicative of rotation. This suggestive rotation is also seen in the position–velocity diagram shown in Figure 5, as the source is resolved across ∼2.4 beams.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The integrated line fluxes are given in Table 2. We measure the FWHM and estimate uncertainties of both the ¹²CO(J = 6→5) and p-H₂O (2_1,1 → 2_0,2) features by using Monte Carlo simulations where noise is injected and the line width remeasured. For single-profile Gaussian fits to the integrated line luminosities, we measure widths of 700 ± 70 km s⁻¹ in ¹²CO(J = 6→5) and 900 ± 200 km s⁻¹ in p-H₂O (2_1,1 → 2_0,2). When analyzing the data for components A and B separately, we measure single-profile Gaussian widths of 260 ± 40 km s⁻¹ for component A and 280 ± 130 km s⁻¹ for component B in ¹²CO(J = 6→5). However, we note that component A is best fit by a double Gaussian separated by 400 km s⁻¹ and individual line widths FWHM = 370 km s⁻¹.

Table 2.  Derived Properties of the Mambo-9 System

Derived	Units	Component	Component	Total
Property		A	B	A+B
R.A.	⋯	10:00:26.356	10:00:26.356	⋯
Decl.	⋯	+02:15:27.94	+02:15:26.63	⋯
From ALMA Spectroscopy:
z	⋯	5.850	5.852	5.850
ICO(6-5)	Jy km s−1	0.43 ± 0.03	0.07 ± 0.02	0.48 ± 0.03
σv(CO)	km s−1	260 ± 40a	280 ± 130	700 ± 70
${L}_{\mathrm{CO}(6-5)}^{{\prime} }$	K km s−1 pc2	(1.4 ± 0.9) × 1010	(2.3 ± 0.7) × 109	(1.56 ± 0.10) × 1010
${I}_{{{\rm{H}}}_{2}{\rm{O}}({2}_{\mathrm{1,1}}-{2}_{\mathrm{0,2}})}$	Jy km s−1	0.09 ± 0.02	⋯	0.09 ± 0.02
${\sigma }_{v}({{\rm{H}}}_{2}{\rm{O}})$	km s−1	900 ± 200	⋯	900 ± 200
${L}_{{{\rm{H}}}_{2}{\rm{O}}({2}_{\mathrm{1,1}}-{2}_{\mathrm{0,2}})}^{{\prime} }$	K km s−1 pc2	(2.5 ± 0.5) × 109	⋯	(2.5 ± 0.5) × 109
Mgas(CO)	M⊙	(3.3 ± 1.7) × 1011	(5 ± 3) × 1010	(3.7 ± 1.8) × 1011
Mdyn	M⊙	(2.0 ± 0.8) × 1011	(7 ± 6) × 1010	⋯
From ALMA Dust Continuum:
Mdust	M⊙	(1.3 ± 0.3) × 109	(${1.9}_{-0.8}^{+1.3}$) × 108	(${1.6}_{-0.3}^{+0.4}$) × 109
${M}_{\mathrm{gas}}(3\,\mathrm{mm})$	M⊙	(1.4 ± 0.4) × 1011	(1.2 ± 0.5) × 1010	(1.7 ± 0.4) × 1011
FWHMmaj(870 μm)b	''	015 ± 001	030 ± 005	⋯
Axis Ratio (870 μm)b b/a	⋯	0.87 ± 0.15	0.37 ± 0.24	⋯
Reff(870 μm)	pc	380 ± 30	760 ± 130	⋯
FWHMmaj(3 mm)b	''	029 ± 011	⋯	⋯
Reff(3 mm)	pc	700 ± 300	⋯	⋯
${{\rm{\Sigma }}}_{{M}_{\mathrm{dust}}}$	M⊙ pc−2	1400 ± 400	50 ± 30	⋯
ΣSFR	M⊙ yr−1 kpc−2	640 ± 170	61 ± 35	⋯
From Broadband SED Fitting:
λrest at which τν = 1	μm	$\equiv 200$	Opt. Thin	⋯
LIR	L⊙	(${4.0}_{-0.7}^{+0.9}$) × 1012	(${1.5}_{-0.5}^{+1.1}$) × 1012	(${6.3}_{-0.9}^{+1.1}$) × 1012
SFR	M⊙ yr−1	${590}_{-100}^{+140}$	${220}_{-70}^{+150}$	${930}_{-130}^{+160}$
λpeak	μm	87 ± 7	${100}_{-80}^{+30}$	⋯
Tdust	K	${56.3}_{-5.7}^{+5.9}$	${29.7}_{-6.6}^{+8.5}$	⋯
β	⋯	1.95 ± 0.11	${2.66}_{-0.34}^{+0.22}$	⋯
M⋆ (OIR-only)	M⊙	⋯	⋯	(${3.2}_{-1.5}^{+1.0}$) × 109
LUV(1600 Å)	L⊙	<7.7 × 109	(8.8 ± 3.6) × 109	<3.4 × 1010
IRX	⋯	>510	${160}_{-100}^{+280}$	⋯
AUV	⋯	>6.2	${5.0}_{-1.1}^{+1.0}$	⋯
qIR	⋯	0.4 ± 1.0	⋯	⋯
Mhalo	M⊙	⋯	⋯	(3.3 ± 0.8) × 1012

Notes. Positions measured from 870 μm dust continuum. Note that quantities derived for the total Mambo-9 system (A+B) are derived independently from the measurements of the two individual systems. In other words, Total is not simply the sum of the two, but a direct independent measurement of the integrated quantity. A brief guide to derived properties follows: I_CO(6–5) is the integrated line flux of the ¹²CO(J = 6→5) line, σ_v(CO) is the velocity dispersion of the ¹²CO(J = 6→5) feature, and ${L}_{\mathrm{CO}(6-5)}^{{\prime} }$ is the ¹²CO(J = 6→5) line luminosity. All three quantities are similarly derived for the p-H₂O (2_1,1 → 2_0,2) line. M_gas(CO) is the gas mass as derived from the ¹²CO(J = 6→5) line, while ${M}_{\mathrm{gas}}(3\,\mathrm{mm})$ is the gas mass as derived from the 3 mm dust continuum (and M_dust is the dust mass derived from the 3 mm continuum). Both assume α_CO = 6.5 M_⊙(K km s⁻¹ pc²)⁻¹. M_dyn is the dynamical mass as estimated from the ¹²CO(J = 6→5) line width and 870 μm dust continuum size. FWHM_maj is the measured FWHM size measured from dust continuum images at 870 μm or 3 mm (in the image plane), and the axis ratio indicates the relative elongation on the plane of the sky. R_eff is the circularized half-light radius in parsecs. ${{\rm{\Sigma }}}_{{M}_{\mathrm{dust}}}$ and Σ_SFR are the dust mass surface density and star formation surface density, respectively. λ_rest is the wavelength at which τ = 1 for the dust SED, L_IR is the derived IR luminosity integrated from 8 to 1000 μm, and SFR is the star formation rate converted directly from L_IR using the Kennicutt & Evans (2012) scaling. λ_peak is the rest-frame wavelength where the dust SED peaks, and T_dust is the underlying dust temperature used in the model fit to the photometry. β is the emissivity spectral index. M_⋆ is the stellar mass of the aggregate system, while L_UV is the rest-frame 1600 Å luminosity. IRX is the ratio L_IR/L_UV, A_UV is the absolute magnitude of attenuation inferred at 1600 Å, q_IR is the implied far-infrared (FIR)-to-radio ratio as in Yun et al. (2001), and M_halo is the total inferred halo mass.

^aComponent A is best characterized by double-peaked emission for which the line width of each component has width of 370 km s⁻¹. ^bFWHM of the major axis measured from a two-dimensional Gaussian fit in the image plane, and the axis ratio is from the same fit.

Download table as:  ASCIITypeset image

The p-H₂O (2_1,1 → 2_0,2) line at rest-frame 752 GHz is a medium-excitation (E_up = 136 K) transition of para-H₂O, commonly seen in emission in galaxies as a tracer of dense ($n(H)\sim {10}^{5}\mbox{--}{10}^{6}\,{\mathrm{cm}}^{-3}$), star-forming gas in the ISM (González-Alfonso et al. 2010; Liu et al. 2017; Jarugula et al. 2019; Apostolovski et al. 2019; Yang et al. 2019). Though rare in the gas phase of non-star-forming molecular clouds (Caselli et al. 2010), water is the third-most-common molecule after H₂ and CO in shock-heated regions of the ISM that trace star-forming regions (Bergin et al. 2003). Furthermore, the velocity structure of H₂O emission in nearby galaxies tends to mirror that of CO (Liu et al. 2017), suggesting that water is widespread throughout the bulk molecular gas reservoir of galaxies. This particular transition tends to be relatively bright compared to most water emission features.

The continuum-subtracted line luminosity of the p-H₂O (2_1,1 → 2_0,2) feature is (3.6 ± 0.8) × 10⁷ L_⊙, precisely on the L_IR–${L}_{{{\rm{H}}}_{2}{\rm{O}}}$ relation found in Yang et al. (2013) using our best-constrained L_IR for component A from Section 3.7.2; this corroborates earlier results that suggest H₂O might be a particularly good star formation tracer. Furthermore, the ratio of line flux between p-H₂O (2_1,1 → 2_0,2) and ¹²CO(J = 6→5) is ∼30%, consistent to within 10% of the composite DSFG millimeter spectrum from the South Pole Telescope (SPT) survey (Spilker et al. 2014).

Dust continuum detections on the Rayleigh–Jeans tail of blackbody emission can be used to directly infer Mambo-9's total dust mass (and also ISM mass by proxy, as discussed in the next subsection). Dust mass is proportional to dust temperature and flux density along the Rayleigh–Jeans tail, where dust emission is likely to be optically thin (at λ_rest ≳ 300 μm). As we discuss later in Section 3.7.2, we estimate that this is a safe assumption to make in the case of Mambo-9, where we do not think the SED is optically thick beyond λ_rest ≈ 300 μm.²⁷ Because Mambo-9 sits at a relatively high redshift, cosmic microwave background (CMB) heating of the dust is nonnegligible, and the subsequent measurement of dust mass is impacted.

For the general case of a galaxy at sufficiently high z like Mambo-9, dust mass can be calculated using the following:

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

Here, observations are acquired at ν_obs (measured in Hz) with flux density ${S}_{{\nu }_{\mathrm{obs}}}$ (measured in erg s⁻¹ cm⁻² Hz⁻¹), and ν_rest = ν_obs(1+z). The frequency at which the dust mass absorption coefficient is known is ν_ref, for example, a value of $\kappa (450\,\mu {\rm{m}})\,=1.3\pm 0.2$ cm² g⁻¹ (Li & Draine 2001; Weingartner & Draine 2001), which is observed-frame 3 mm at z = 6. Also, D_L is the luminosity distance (converted to cm), and B_ν is the Planck function evaluated at a given frequency and for a given temperature in units of erg s⁻¹ cm⁻² Hz⁻¹. For example, B(ν_rest, T_CMB(z)) is the Planck function evaluated at ν_rest for the temperature of the CMB at the measured redshift z. Here, M_dust is in units of g, which can be converted to M_⊙; β is the emissivity spectral index, which we set to β = 1.95 (we derive this value from a fit to our data in Section 3.7.2); Γ_RJ represents the Rayleigh–Jeans (RJ) correction factor, or the deviation from the RJ approximation, following the framework of Scoville et al. (2016):

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{RJ}}(\nu ,{T}_{{\rm{d}}},z)=\displaystyle \frac{h\nu (1+z)/{{kT}}_{{\rm{d}}}}{{e}^{h\nu (1+z)/{{kT}}_{{\rm{d}}}}-1}\end{eqnarray} \tag{ 2 }$$

and ${{\rm{\Gamma }}}_{\mathrm{RJ}(\mathrm{ref},0)}={{\rm{\Gamma }}}_{\mathrm{RJ}}(\nu ={\nu }_{\mathrm{ref}},{T}_{{\rm{d}}}={T}_{\mathrm{dust}},z=0)$. Here, h is the Planck constant (6.63×10¹⁷ erg s); k is the Boltzmann constant (1.38 × 10⁻¹⁶ erg K⁻¹); T_d is the galaxy's mass-weighted dust temperature, not the same quantity as fit in Section 3.7.2, which is the luminosity-weighted dust temperature. We adopt a mass-weighted dust temperature of 25 K throughout to be consistent with Scoville et al. (2016). The last multiplicative factor in Equation (1) represents the correction for suppressed flux density against the CMB background, as described in da Cunha et al. (2013). This factor is a function of ν_rest = ν_obs(1+z), the CMB temperature at the given redshift, T_CMB = 2.73 K(1+z), and the CMB-corrected mass-weighted dust temperature, as given in Equation (12) of da Cunha et al. (2013). An assumption of this formulation is that the dust (at temperatures similar to the CMB temperature) is optically thin, which holds in almost all environments, with the exception of the densest cores of local ultraluminous infrared galaxies.

We derive dust masses from our 3 mm continuum photometry centered on a wavelength of 3085 μm (rest-frame 450 μm). We infer a dust mass of (1.3 ± 0.3) × 10⁹M_⊙ for component A and (${1.9}_{-0.8}^{+1.3}$) × 10⁸M_⊙ for component B. Note that the sum of these values is ∼5× higher than the dust mass derived for this system in Jin et al. (2019), with the difference attributed to the differences in SED fitting including the best-fit β (our dust mass uncertainties do not account for the measurement uncertainties in β).

We derive the mass of molecular gas in each galaxy from the dust continuum as well as from ¹²CO(J = 6→5) line luminosity. In both cases, we adopt a value for the CO-to-H₂ conversion factor of α_CO = 6.5 M_⊙ (K km s⁻¹ pc²)⁻¹ as in Scoville et al. (2016). This is in line with the Galactic value and accounts for the mass of both H₂ and He gas.²⁸

We follow the methodology described in the appendix of Scoville et al. (2016) to derive a gas mass from the dust continuum, modified to account for CMB heating similar to the effect on the dust mass calculation:

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{gas}} & = & \displaystyle \frac{4\pi {D}_{{\rm{L}}}^{2}{S}_{\nu ,\mathrm{obs}}}{\alpha {({\nu }_{850})(1+z)}^{3+\beta }}{\left(\displaystyle \frac{{\nu }_{850}}{{\nu }_{\mathrm{obs}}}\right)}^{2+\beta }\left(\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{RJ}(\mathrm{ref},0)}}{{{\rm{\Gamma }}}_{\mathrm{RJ}}}\right)\\ & & \times \,{\left(1-\displaystyle \frac{B({\nu }_{\mathrm{rest}},{T}_{\mathrm{CMB}}(z))}{B({\nu }_{\mathrm{rest}},{T}_{\mathrm{dust}})}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 3 }$$

Here, α₈₅₀ = (6.7 ± 1.7) × 10¹⁹ erg s⁻¹ Hz⁻¹ M_⊙⁻¹ is the empirically calibrated conversion factor from 850 μm luminosity to ISM mass from Scoville et al. (2016), which bypasses use of both the uncertain dust-to-gas ratio and the dust mass absorption coefficient (used to measure dust masses above). Note that intrinsic to this calculation is the assumed value of the CO-to-H₂ conversion factor as stated above. Similar to the calculation of dust mass, we adopt a single mass-weighted dust temperature of 25 K. Note that M_gas as given in Equation (3) is in units of M_⊙. With this approach, we constrain masses of molecular gas using the dust continuum to be M_gas = (1.4 ± 0.4) × 10¹¹M_⊙ for component A and M_gas = (1.2 ± 0.5) × 10¹⁰ M_⊙ for component B.

Historically, it has been more common to use transitions of CO to infer the underlying gas mass in galaxies (Solomon & Vanden Bout 2005; Carilli & Walter 2013), yet it comes with substantial uncertainty, especially when using high-J transitions like ¹²CO(J = 6→5). High-J transitions of CO tend to trace dense gas regions of the ISM, which are a relatively poor probe of the entire molecular gas reservoir of a galaxy. Nevertheless, here we offer a calculation of the gas mass from ¹²CO(J = 6→5) as an independent check against what we have calculated using the dust continuum.

We use the ¹²CO(J = 6→5) line luminosity²⁹ to derive a molecular gas mass from ¹²CO(J = 6→5). This requires an assumption as to the value of the gas excitation spectral line energy distribution (SLED) to convert from ¹²CO(J = 6→5) to the ground-state ¹²CO(J = 1→0) and then the value of the CO-to-H₂ conversion factor (Bolatto et al. 2013). We assume that Mambo-9 has a CO SLED similar to other high-z DSFGs in the literature (summarized in Figure 45 of Casey et al. 2014, including a substantial contribution from the compilation of Bothwell et al. 2013; we adopt the blue shaded region from that figure as the 1σ uncertainty on ${I}_{\mathrm{CO}(6-5)}/{I}_{\mathrm{CO}(1-0)}={10}_{-5}^{+30}$). We calculate gas masses of M_gas(A) = (3.3 ± 1.7) × 10¹¹M_⊙ and M_gas(B) = (5 ± 3) × 10¹⁰M_⊙ for components A and B, respectively. These are broadly consistent with, yet more uncertain than, the dust-continuum-derived gas masses.

Later in Section 4.1 we discuss the implications of the rarity of this halo on the measured gas mass, in particular the assumption that ${\alpha }_{\mathrm{CO}}=6.5$ M_⊙(K km s⁻¹ pc²)⁻¹. Both calculations of the gas masses account for measurement uncertainties in flux density, α₈₅₀ or ${I}_{\mathrm{CO}(6-5)}/{I}_{\mathrm{CO}(1-0)}$, but not uncertainty in α_CO. If we were to instead take α_CO = 1 M_⊙(K km s⁻¹ pc²)⁻¹, more in line with measured constraints on low- and high-redshift dusty starbursts (e.g., Downes & Solomon 1998; Tacconi et al. 2008; Bolatto et al. 2013), the gas mass would scale down proportionally by a factor of 6.5. This would give us gas masses of M_gas = (2.2 ± 0.6) × 10¹⁰ M_⊙ for component A and M_gas = (1.8 ± 0.8) × 10⁹M_⊙ for component B (both scaled down from the dust-continuum-estimated gas masses).

We measure resolved sizes for both components A and B multiple ways to check consistency. First we fit Sérsic profiles to the two components of the highest resolution and S/N data we have: the self-calibrated 870 μm data using Briggs weighting with robust = 0.0. This probes rest-frame ∼127 μm, near the peak of the long-wavelength SED, so these sizes trace the star-forming region in Mambo-9 most closely. We perform this analysis in the uv plane following the methodology outlined in Spilker et al. (2016). Because both components are only marginally resolved, the size is measured with a fixed Sérsic index of n = 1 consistent with an exponential disk (though we found that Gaussian sizes, with n = 0.5, are consistent with those measured for n = 1). We measure circularized half-light radii of R_e(A) = 0068 ± 0002 = 408 ± 12 pc and R_e(B) = 0110±0010 = 660 ± 60 pc. Uncertainties do not account for the unconstrained Sérsic index, which was fixed to n = 1. We compare these sizes to the deconvolved two-dimensional Gaussian sizes measured in the image plane following the methodology of Simpson et al. (2015) and Hodge et al. (2016); we find that, though the image plane fits are somewhat sensitive to image weighting, they are broadly consistent with the uv-plane analysis, with R_e(A) = 0064 ± 0004 = 380 ± 30 pc and R_e(B) = 0128 ± 0021 = 760 ±130 pc (both of these quoted values are for Briggs weighting). Note that, even though the synthesized beam of these data is larger than the measured sizes, the very high S/N enables us to measure half-light radii sufficiently smaller than the beam size FWHM.

Though our 3 mm data are not at the same high S/N as the 870 μm data and also have lower spatial resolution, we fit sizes to both continuum and CO moment-0 maps of component A (using Briggs weighting with robust = 0.0) to explore possible differences using the different tracers. Unfortunately, component B is not detected at sufficiently high S/N to have measurable sizes in either 3 mm continuum or CO. We measure a 3 mm continuum circularized half-light radius in the image plane of R_e(A) = 0123 ± 0047 = 700 ± 300 pc, only 1σ discrepant from the measured size at 870 μm. The 3 mm continuum is a probe of the cold dust in the ISM. Though we infer that the CO moment-0 map is marginally resolved, we find that component A is consistent with both a point source and the measured 3 mm and 870 μm sizes; there is significant uncertainty in the CO size, due to the lower S/N and poor spatial resolution.

We also use the measured sizes to estimate the average dust column densities within the half-light radius, which further inform the conditions of the ISM in Mambo-9; we adopt the measured 870 μm R_eff sizes due to the high S/N near the peak of dust emission and the measured dust masses to calculate ${{\rm{\Sigma }}}_{{M}_{\mathrm{dust}}}(A)\,=1400\pm 400$ M_⊙ pc⁻² and ${{\rm{\Sigma }}}_{{M}_{\mathrm{dust}}}(B)=50\pm 30$ M_⊙ pc⁻² for A and B, respectively. If a Milky Way-type dust is assumed (as tabulated in Table 6 of Li & Draine 2001), these dust mass column densities imply that the dust SED should likely be opaque to rest-frame wavelengths of ∼200–400 μm in the case of component A and ∼25–70 μm in the case of component B. Similarly, we can estimate the star formation surface density (using the SFRs derived for each component later in Section 3.7.2), and we arrive at Σ_SFR(A) = 640 ± 170 M_⊙ yr⁻¹ kpc⁻² and 60 ± 35 M_⊙ yr⁻¹ kpc⁻². Neither is near the hypothetical Eddington limit for starbursts that are factors of several larger (Scoville et al. 2001; Scoville 2003; Murray et al. 2005; Thompson et al. 2005), though some have recently pointed out that the limit might be much higher yet considering starbursts are distributed over some area and are not point sources.

We derive the dynamical masses of Mambo-9-A and Mambo-9-B using the ¹²CO(J = 6→5) kinematic profile, comparing a few different methods. First, for a galaxy with an unresolved velocity field, the dynamical mass is best estimated by

$$\begin{eqnarray}&&{M}_{\mathrm{dyn}}=\displaystyle \frac{5{R}_{{\rm{e}}}{\sigma }^{2}}{G},\end{eqnarray} \tag{ 4 }$$

where G is the gravitational constant, σ is the measured velocity dispersion of the kinematic feature measured (in our case, ¹²CO(J = 6→5)), and R_e is the effective circularized radius, and the factor of five is a constant of proportionality determined to best represent galaxies in the local mass plane (Cappellari et al. 2006; Toft et al. 2017); this constant does not account for the inclination angle, i. Correcting for unknown inclination requires an additional factor of 3/2 (which is the reciprocal of the expectation value of sin²i).

Because the size measurements for component A and component B are broadly consistent and we lack data for a more detailed analysis, we use the measured 870 μm dust-emitting sizes to estimate the galaxies' dynamical masses. There are a number of potential caveats in doing this. First, the dynamical mass is best measured in the same tracer used to infer the galaxy's kinematics (¹²CO(J = 6→5) in this case). Second, using a high-J tracer like ¹²CO(J = 6→5) would likely bias the dynamical mass estimate because it only probes dense gas regions. While both of these concerns are important to keep in mind, a few facts provide reassurance that our assumptions are sufficient in this case: first, the fact that the galaxy's 3 mm size is not significantly larger than its 870 μm size, and second, the fact that the uncertainty on the measured quantities dominates the calculation of the dynamical mass. In other words, the uncertainties on R_e and σ, combined with uncertainty in unconstrained i, are significant enough to dominate over variations in tracer-dependent forms of these quantities.

Thus, we adopt circularized effective radii of R_e(A) = 380 ± 30 pc and R_e(B) = 760 ± 130 pc. The velocity dispersions as measured from ¹²CO(J = 6→5) (and as shown in Figure 4) are σ_V(A) = 260 ± 40 km s⁻¹ and σ_V(B) = 280 ± 130 km s⁻¹. This gives dynamical mass estimates of M_dyn(A) = (5 ± 2) × 10¹⁰ M_⊙ and M_dyn(B) = (7 ± 6) × 10¹⁰ M_⊙, respectively. Though the dynamical mass estimated for component B is a bit larger (due to its larger physical size), the uncertainty is quite large and consistent with being an equal or smaller mass companion. While the mass calculated for component A seems rather precise, it should be noted that the double-peaked ¹²CO(J = 6→5) spectrum of component A is poorly fit to a single Gaussian component. Note that if we instead calculate a dynamical mass from the p-H₂O (2_1,1 → 2_0,2) line width, which is much more uncertain, we get dynamical mass estimates an order of magnitude larger; as Figure 4 shows, this is not because the p-H₂O (2_1,1 → 2_0,2) line is much more broad than the ¹²CO(J = 6→5) line, but it represents the difference between a single and double component fit.

Alternatively, the dynamical mass of component A could be estimated directly from the resolved ¹²CO(J = 6→5) kinematics using

$$\begin{eqnarray}&&{M}_{\mathrm{dyn}}=\displaystyle \frac{{V}_{\max }^{2}\,{R}_{\max }}{G},\end{eqnarray} \tag{ 5 }$$

which then similarly needs to be corrected for unknown inclination. Note that R_max here denotes the maximum radius at which V_max is measured and differs from R_e, the circularized half-light radius, used above. Using both size R_max = 085 ± 020 = 5.0 ± 1.2 kpc and V_max = 350 ± 50 km s⁻¹ measurements from the position–velocity diagram in Figure 5, we derive an alternate dynamical mass of component A of (2.0 ± 0.8) × 10¹¹M_⊙. The dynamical mass for component B cannot be constrained using this method because of the low S/N of the line and no resolved rotation curve.

Given the caveats of using a different tracer for size measurements, we adopt the more conservative higher-mass dynamical constraint for component A of (2.0 ± 0.8) × 10¹¹M_⊙, but cover the implications of either dynamical mass constraint in our discussion of the total mass budget in Section 4.1. For component B, we adopt the only estimated, yet highly uncertain, value of M_dyn of (7 ± 6) × 10¹⁰ M_⊙.

We fit spectral energy distributions using three methods: the stellar component only (as in Finkelstein et al. 2015), the obscured component only (as in Casey 2012), and both together using energy balance techniques (specifically magphys; da Cunha et al. 2008, 2015). The three approaches, used to derive different physical quantities, are described below. The derived properties are given in Table 2. As Figure 6 shows, the final SED we adopt for Mambo-9 is outlined in black; the details are described throughout this section.

Zoom In Zoom Out Reset image size

3.7.1. Optical/Near-infrared SED Fit

3.7.2. Far-infrared/Millimeter SED Fit

The obscured SED (probed by rest-frame wavelengths ∼5–3000 μm) has no detections at wavelengths shortward of 850 μm. Nevertheless, due to the superb quality of the ALMA continuum detections, we can fit a somewhat well-constrained obscured SED with a single modified blackbody plus mid-infrared power law. The mid-infrared component is unconstrained, due to the dearth of detections shortward of the peak, but we adopt a model that follows ${S}_{\nu }\propto {\nu }^{-\alpha }$ (the power law joins the modified blackbody at the point where the derivative is equal to α, as described in Casey 2012). Here we fix the value of α to α_MIR = 4. Physically, a lower value of α represents a higher proportion of emission emanating from dust heated from discrete sources rather than the cooler dust heated by the ambient radiation field in the ISM. Very high values of α asymptote to the pure modified blackbody solution. While no direct constraints can be made for α_MIR, it should be noted that values less than α_MIR = 4 violate the upper limits in the mid-infrared as measured by Spitzer and Herschel. If α_MIR is constrained to values in excess of ∼4, then the total impact on the integrated L_IR is negligible. The SED is fit using a simple Metropolis Hastings Markov Chain Monte Carlo with free parameters of λ_peak, the rest-frame peak wavelength, L_IR, the integrated 8–1000 μm IR luminosity, and β, the emissivity spectral index. This is an updated fitting technique that largely follows the methodology outlined in Casey (2012) but substitutes least-squares fitting for Bayesian analysis and a contiguous function for a piecewise power law and modified blackbody (this will be described in a forthcoming paper by P. Drew). This fit embeds the impact of CMB heating on ISM dust at high redshift as prescribed in da Cunha et al. (2013); at z = 5.85, the CMB is 18.7 K. Note that λ_peak, L_IR, and β (the three free parameters of the fit) are measured from the intrinsic emitted SED, not the observed SED and observed photometry, which has been affected by the CMB.

Figure 7 shows the results of the obscured SED fit for both an optically thin and a more general opacity model for the two major components of this source. The general opacity model assumes τ = 1 at rest-frame 1.5 THz (or 200 μm; Conley et al. 2011). Due to the high S/N of the ALMA measurements, in particular the 870 μm data near the SED peak, the long-wavelength portion of the SED and the peak are more precisely constrained than for most DSFGs and also allow for an independent measurement of β, the emissivity spectral index. For component B, we set an upper limit to β = 3 based on the low S/N of the source's photometry. The lower left section of each component fit shows a corner plot of the converged MCMC chains in λ_peak, L_IR, and β. Both optically thin and general opacity models are significantly higher quality for component A than component B. The upper middle panel places the measured L_IR and λ_peak values in context against (a) the z ∼ 1−2 L_IR–λ_peak relationship derived for DSFGs in Casey et al. (2018b), (b) the measured characteristics of z > 4 DSFGs from the SPT survey (Strandet et al. 2016), and (c) in contrast to expectation for z ∼ 6 galaxies from theoretical modeling (Ma et al. 2019). Note that while several modeling papers (e.g., Behrens et al. 2018; Liang et al. 2019; Ma et al. 2019) suggest that the luminosity-weighted dust temperatures of galaxies at z ≳ 5 should be warmer than those at z ∼ 2, we do not see compelling evidence that this holds for either component of Mambo-9.

Zoom In Zoom Out Reset image size

The parameter λ_peak is preferred over a direct fit to the physical dust temperature because it is insensitive to the opacity model assumed; in other words, an SED that peaks at rest-frame 90 μm could have an intrinsic dust temperature ranging anywhere from 30 to 50 K depending on the geometry and column density of the dust. But because Mambo-9 sits at such a high redshift where CMB heating is nonnegligible, the opacity model assumptions do affect the intrinsic rest-frame peak wavelength λ_peak. Fit to the same photometry, optically thin dust SEDs will consistently have lower dust temperatures than more general opacity assumptions that allow for dust self-absorption near the peak, so they are proportionally more affected by CMB heating. Thus, the difference in measured λ_peak in Figure 7 between opacity models is purely due to the different levels of effects of the CMB.

While the CMB does affect λ_peak by way of the underlying physical dust temperature, the gap between L_IR that is fit with different opacity models is smaller than what it would be in the absence of the CMB or at lower redshifts. As shown on the right-hand panels of Figure 7, the difference between the emitted SED (dashed line) and observed SED (dark solid line) is much larger for the optically thin SED (purple) than for the general opacity model (orange), so while the CMB has little effect on the derived L_IR of the general opacity model fits, it has a small but measurable effect on L_IR of the optically thin model.

Through analysis of the dust mass surface density from the 870 μm data, we have roughly constrained the wavelength at which the SED becomes optically thick. A value of ${{\rm{\Sigma }}}_{{M}_{\mathrm{dust}}}\,=1400\pm 400$ M_⊙ pc⁻² in component A suggests an optically thick SED to ∼200–300 μm rest frame, while the lower dust column density in component B of ${{\rm{\Sigma }}}_{{M}_{\mathrm{dust}}}=50\pm 30$ M_⊙ pc⁻² suggests the SED is optically thin near the peak (these measurements assume the dust mass absorption coefficients as given in Li & Draine 2001). These different dust column densities imply that the dust SEDs of the two components should be treated differently, with component A more reminiscent of the very high column densities of dust that is ubiquitous among the brightest DSFGs at lower redshifts. Thus, we adopt the more general opacity model (with τ = 1 at λ_rest = 200 μm) for component A and the optically thin model for component B.

The implied SFRs from the dust emission, assuming the Kennicutt & Evans (2012) scaling (which uses an IMF from Kroupa & Weidner 2003), are ${590}_{-100}^{+140}$ M_⊙ yr⁻¹ for component A and ${220}_{-70}^{+150}$ M_⊙ yr⁻¹ for component B. Note that the total SFR fit to the aggregate SED (${930}_{-130}^{+160}$ M_⊙ yr⁻¹) is lower than, though not fully inconsistent with, the derived SFR of the system in Jin et al. (2019) of 1283 ± 173 M_⊙ yr⁻¹. While the dust temperature of the general opacity model may seem high compared to some DSFGs in the literature, Figure 7 shows they are fully consistent with the measured λ_peak values for both lower-redshift DSFGs (as derived in Casey et al. 2018b) and for existing measurements of z > 4 DSFGs from the SPT survey (Strandet et al. 2016). They are both colder (i.e., higher λ_peak) than the expected median L_IR–λ_peak relation from simulations at z ∼ 5.9 (Ma et al. 2019).

There is a slight discrepancy between the measured emissivity spectral index, β, of components A and B of β(A) = 1.95 ± 0.11 and $\beta (B)={2.66}_{-0.34}^{+0.22};$ while this could be a real discrepancy, the quality of the constraint for component B is weak, and there is a significant degeneracy between β and λ_peak in the absence of high S/N data on the Rayleigh–Jeans tail. Note that Jin et al. (2019) conclude that the CMB might be affecting the net SED by artificially steepening β for Mambo-9, but we do not see evidence for this in the case of component A, and we only see weak evidence that this is the case for component B. In other words, if we fit the SED of component B without accounting for CMB heating, we would measure a steeper value of β = 2.89 ± 0.40 than we do having accounted for the CMB. The difference in our conclusions regarding component A is driven by the inclusion (or not) of the single-dish photometry from Herschel (in particular the deblended photometry), AzTEC, and SCUBA-2. The band 7 data presented in this paper results in a much less steep Rayleigh–Jeans tail and derived value of β consistent with the often-used assumption in the literature of β = 1.5−2.0.

In conclusion, we infer that component A is optically thick at the peak, has a measured L_IR = (${4.0}_{-0.7}^{+0.9}$) × 10¹² L_⊙, $\mathrm{SFR}={590}_{-100}^{+140}$ M_⊙ yr⁻¹, rest-frame peak wavelength of λ_peak = 87 ± 7 μm, dust temperature ${T}_{\mathrm{dust}}={56.3}_{-5.7}^{+5.9}$ K, and emissivity spectral index β = 1.95 ± 0.11. We infer that component B is optically thin at the peak, has a measured L_IR = (${1.5}_{-0.5}^{+1.1}$) × 10¹² L_⊙, $\mathrm{SFR}={220}_{-70}^{+150}$ M_⊙ yr⁻¹, rest-frame peak wavelength ${\lambda }_{\mathrm{peak}}={100}_{-80}^{+30}$ μm, dust temperature ${T}_{\mathrm{dust}}\,={29.7}_{-6.6}^{+8.5}$ K, and emissivity spectral index $\beta ={2.66}_{-0.34}^{+0.22}$.

3.7.3. Energy-balance SED Fit

In the radio regime, Mambo-9 component A is marginally detected at 3.2σ significance in the very deep 3 GHz map (Smolčić et al. 2017) and not detected at 1.4 GHz (Schinnerer et al. 2007). In the absence of direct constraints on the radio SED, we adopt a synchrotron slope of α_rad = −0.8 (consistent with Condon 1992, often used in the literature when there is a dearth of multiple radio continuum constraints) fit to the 3 GHz photometry. The far-infrared/radio correlation for star formation would predict radio emission below the detection limit in the case of both a nonevolving relation (∼50 nJy, expected q_IR = 2.4; Yun et al. 2001) and evolving relation (∼300 nJy, expected q_IR = 2.0; Delhaize et al. 2017). The marginal 3 GHz detection implies a value of q_IR = 0.4 ± 1.0, suggestive of a possible buried active galactic nucleus (AGN).

Though suggestive of an AGN, it should be noted that the Magphys SED, which assumes a nonevolving FIR–radio correlation of q_IR = 2.34 (da Cunha et al. 2015), intersects our measured 3 GHz radio continuum measurement without the need to invoke a possible radio AGN. This is due to a higher L_IR for the Magphys fit, largely driven by excess emission in the mid-infrared regime, where we lack data. It is also due to the assumed synchrotron slope in the radio regime, closer to α = − 0.65, and incorporating emission from free–free emission. The combination of these different assumptions leads to a star-formation-dominated radio SED in line with measurements.

Whether or not there is an AGN in Mambo-9 is unclear. It is not detected in the Chandra X-ray imaging in COSMOS, though no detection would be expected at the given depth for a z = 5.85 buried AGN. Future observations of the CO SLED, other (sub)millimeter emission features, and rest-frame optical emission lines from the James Webb Space Telescope (JWST) will play a crucial role in inferring whether or not such a buried AGN exists.

Figure 8 depicts the total mass budget of the Mambo-9 system with different sets of assumptions drawn from calculations in the previous sections, from the most massive set of assumptions to the least massive set of assumptions. All analysis here accounts for both Mambo-9-A and Mambo-9-B together and primarily hinges on the adopted gas and stellar masses. The most massive gas mass derived for Mambo-9 uses α_CO = 6.5M_⊙(K km s⁻¹ pc²)⁻¹, and the most massive stellar mass uses that derived using Magphys. The least massive set of assumptions uses the more modest value of α_CO = 1M_⊙(K km s⁻¹ pc²)⁻¹ and the OIR-only derived stellar mass from Section 3.7.1. The dust masses are the same for all models, as this primarily depends on κ_ν, the dust mass absorption coefficient and dust temperature, for which there is little evidence for significant dynamic range in different environments. The baryonic masses are then derived from the sum of these three components.

Zoom In Zoom Out Reset image size

The halo masses are then extrapolated from the total baryonic masses using two methods, both rather conservative; these are the same assumptions used to estimate the halo mass of SPT0311, the most distant DSFG (Marrone et al. 2018). The first sets an absolute floor to the halo mass by assuming the cosmic baryon fraction f_b = 0.19 (Planck Collaboration et al. 2016) and converting from the median baryonic mass total calculated in all cases. This absolute lower limit is drawn by a dotted red horizontal line in Figure 8. The alternate technique assumes that baryons, at most, are 5% of the halo mass; in reality, this assumption is somewhat conservative as often baryons make up less than this, but given that this is at quite high redshift and in a regime near the tip of the mass function, 5% is an appropriate conservative assumption (e.g., Behroozi et al. 2018). The red distributions in Figure 8 show the final extrapolated halo mass distributions, ranging from M_halo = (6.0 ± 1.4) × 10¹¹M_⊙ (least massive assumptions) to M_halo = (7.2 ± 1.8) × 10¹²M_⊙ (most massive assumptions).

The gray regions in Figure 8 denote the 1σ, 2σ, and 3σ exclusion curves for our parent survey at z = 5.85 calculated using the analysis described in Harrison & Hotchkiss (2013). In other words, Mambo-9 was selected for follow-up out of our larger 2 mm blank-field survey covering a contiguous area of ∼155 arcmin², with 135 arcmin² at the necessary depth required to have identified Mambo-9 at 5σ significance or greater.³⁰ The exclusion curves represent the rarest, most massive halos that should exist in the 135 arcmin² survey area at any redshift (but whose mass is evolved backward or forward to that at z = 5.85) with 66%, 95%, and 99.7% likelihood. These upper halo mass limits are 2.5 × 10¹²M_⊙, 4.3 × 10¹²M_⊙, and 8.0 × 10¹²M_⊙ for 1σ, 2σ, and 3σ, respectively.

The most massive assumptions provide an estimated halo mass that exceeds the maximum detectable halo mass by 2.8σ; in other words, finding a halo this massive in a survey this small is only 0.5% likely. On the other hand, the least massive assumptions predict a total halo mass that is allowable in the given survey limits. Note that in an intermediate solution where the lower gas mass and higher stellar mass are adopted, the predicted halo mass is still quite high at (4.7 ± 1.4) × 10¹²M_⊙ by virtue of the high stellar mass, with only a 6% likelihood of having identified such a massive system in our 2 mm survey. If we instead adopt the higher gas mass and lower stellar mass, the total halo mass is predicted to be (3.3 ± 0.8) × 10¹² M_⊙, with a 13% likelihood of sitting in our survey volume. Although it might make the most sense, from a mass analysis standpoint, to adopt the lowest mass values, this would lead to an unusual implication: the implied gas-to-dust ratio (GDR) would be anomalously low, at GDR ∼ 16, lower than the vast majority of galaxies in the literature, even at supersolar metallicity (Rémy-Ruyer et al. 2014). For this reason, we adopt the last intermediate mass constraints (i.e., high gas mass, low stellar mass), due to the higher likelihood (13%) of occurring in the survey volume than any other permutation. In other words, we adopt the high gas mass using α_CO = 6.5M_⊙(K km s⁻¹ pc²)⁻¹, the low stellar mass (${3.2}_{-1.5}^{+1.0}$) × 10⁹ M_⊙, and halo mass of (3.3 ± 0.8) × 10¹² M_⊙, as reflected in Table 2.

The implied gas-to-dust ratio (GDR = M_gas/M_dust) for the system is $\mathrm{GDR}={122}_{-43}^{+64}$, in line with the canonical ratio of 100:1 assumed for galaxies with near-solar-metallicity environments (Rémy-Ruyer et al. 2014). Such a GDR is consistent with a solar-metallicity environment, though the metal content has not been directly constrained in Mambo-9. The gas fraction (${f}_{\mathrm{gas}}\equiv {M}_{\mathrm{gas}}/({M}_{\mathrm{gas}}+{M}_{\star }+{M}_{\mathrm{dust}})$) is extraordinarily high at ${f}_{\mathrm{gas}}={96}_{-2}^{+1} \%$ when adopting our current values; though unlike galaxies in the nearby universe, it is not outside of the realm of theory to observe such gas-rich systems at z ∼ 6 or beyond.

Confirming the high gas fraction and implied halo mass requires observations of the unobscured stellar component in the mid-infrared. These refined measurements could, in turn, lead to more physically motivated questions regarding the origins of galaxies like Mambo-9. For example, does dust at z = 5.85 have the same absorptive properties as local universe dust, even if it might have a very different origin and composition (e.g., de Rossi et al. 2018)? What CO-to-H₂ conversion factor holds for Mambo-9?

This source is the first that was found in a blank-field 2 mm map begun in ALMA Cycle 6, designed to efficiently select z > 4 dust-obscured galaxies. So far, the strategy has been effective. The initial map in which this source was found (corresponding to one ALMA scheduling block) has an effective area of ∼9.4 arcmin² with 1σ_rms ≲ 0.12 mJy/beam at 2 mm,³¹ which is the depth required to detect Mambo-9 to 5σ significance or above. Mambo-9 was the only detection above this threshold in this small map, despite several other SCUBA-2 detected DSFGs sitting in the map region, likely indicating those other DSFGs sit at lower redshifts than z ≲ 3. We use the models of Casey et al. (2018a, 2018b) to comment on the possible implications of Mambo-9 on the number density of DSFGs at z ≳ 4. The first "dust-poor" model (Model A therein) represents a universe with very few DSFGs beyond z > 4, and the "dust-rich" model (Model B) represents a universe rich with DSFGs at z > 4. These two models bracket extreme interpretations of measurements in the literature, and we refer the reader to those papers for more thorough discussion. While the dust-poor model would predict only ∼0.5 sources in a map this size and a median redshift of z ∼ 3.5 (± 1σ ranging over 3 < z < 4.5), the dust-rich model would predict four sources in this map and a median redshift of z ∼ 5.5 (± 1σ ranging over 4 < z < 7.5).

While a single source cannot rule either model in or out, or any plausible model in between because of Poisson statistics, it is interesting to note that the first source found using this mapping strategy sits at z = 5.85. This is above the median redshift predicted for both models and in the tail of the redshift distribution predicted for the dust-poor model that is consistent with claims from the rest-frame UV community on the lack of dust at early times. Though tempting to infer that there might be a substantial hidden population of DSFGs at these high redshifts based on the small survey volume, it is also important to point out that star formation in massive galaxies is thought to be more heavily clustered (Chiang et al. 2017) at higher and higher redshifts, such that at z ∼ 7, half of all star formation takes place in the progenitors of z = 0 galaxy clusters. This implies that our surveys of the distant universe will need larger areas to not be susceptible to the effects of cosmic variance. Forthcoming analysis of the full 2 mm map data set will provide a crucial next step in constraining the volume density of galaxies like Mambo-9.

Our data suggest that Mambo-9-A and Mambo-9-B make up a close pair of galaxies separated by 6 kpc with a mass ratio ranging from ∼1:1 to ∼1:10 depending on the tracer. They are likely interacting at this close proximity, and the interaction could play a substantial role in driving gas densities high enough to trigger intense star formation. At the given star formation rates, component A will deplete its star-forming gas in ${\tau }_{\mathrm{depl}}={38}_{-12}^{+16}$ Myr, while component B will deplete its gas in ${\tau }_{\mathrm{depl}}={80}_{-40}^{+160}$ Myr. This starburst episode has the potential to increase the stellar mass by an order of magnitude, though the stellar mass of this system is highly uncertain. The star formation rate surface densities differ an order of magnitude between component A and B, with A consistent with some of the densest star-forming galaxy cores known in the local universe.

Our analysis of the halo rarity of Mambo-9 suggests that it will live in a massive galaxy cluster with ${M}_{\mathrm{halo}}=1.6\times {10}^{15}$ M_⊙ by z = 0, as its halo is already sufficiently massive at z = 5.85 to constitute a node of the cosmic web. Though its fate is unknowable, it is likely that both components of Mambo-9 end up forming the very old stellar population in the core of a brightest cluster galaxy in the universe today.