<figure><img src="image_1.png" /></figure>

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

Caitlin M. Casey¹, Jorge A. Zavala¹, Manuel Aravena², Matthieu Béthermin³, Karina I. Caputi⁴,⁵, Jaclyn B. Champagne¹, David L. Clements⁶, Elisabete da Cunha⁷,⁸,⁹, Patrick Drew¹, Steven L. Finkelstein¹, Christopher C. Hayward¹⁰, Jeyhan S. Kartaltepe¹¹, Kirsten Knudsen¹², Anton M. Koekemoer¹³, Georgios E. Magdis⁵,¹⁴, Allison Man¹⁵, Sinclaire M. Manning¹, Nick Z. Scoville¹⁶, Kartik Sheth¹⁷, Justin Spilker¹, Johannes Staguhn¹⁸,¹⁹, Margherita Talia²⁰, Yoshiaki Taniguchi²¹, Sune Toft⁵,¹⁴, Ezequiel Treister²², and Min Yun²³

¹Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA; cmcasey@utexas.edu

²Núcleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile

³Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France

⁴Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700AV Groningen, The Netherlands

⁵Cosmic Dawn Center (DAWN), Denmark

⁶Imperial College London, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK

⁷International Centre for Radio Astronomy Research, University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia

⁸Research School of Astronomy and Astrophysics, The Australian National University, Canberra ACT 2611, Australia

⁹ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia

¹⁰Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA

¹¹School of Physics and Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA

¹²Department of Earth and Space Sciences, Chalmers University of Technology, Onsala Space Observatory, SE-43992 Onsala, Sweden

¹³Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA

¹⁴Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, DK-2100 Copenhagen, Denmark

¹⁵Dunlap Institute for Astronomy & Astrophysics, 50 St. George Street, Toronto, ON M5S 3H4, Canada

¹⁶California Institute of Technology, MC 249-17, 1200 East California Boulevard, Pasadena, CA 91125, USA

¹⁷NASA Headquarters, 300 E Street SW, Washington, DC 20546, USA

¹⁸The Henry A. Rowland Department of Physics and Astronomy, Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA

¹⁹Observational Cosmology Lab, Code 665, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

²⁰Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, I-40129, Bologna, Italy

²¹The Open University of Japan, 2-11 Wakaba, Mihama-ku, Chiba 261-8586, Japan

²²Instituto de Astrofísica and Centro de Astroingeniería, Facultad de Física, Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile

²³Department of Astronomy, University of Massachusetts Amherst, 710 N. Pleasant Street, Amherst, MA 01003, USA

Received 2019 September 3; revised 2019 October 25; accepted 2019 October 29; published 2019 December 11

Abstract

We present a physical characterization of MM J100026.36+021527.9 (a.k.a. “MAMBO-9”), a dusty star-forming galaxy (DSFG) at $z = 5.850 \pm 0.001$. This is the highest-redshift unlensed DSFG (and fourth most distant overall) found to date and is the first source identified in a new 2 mm blank-field map in the COSMOS field. Though identified in prior samples of DSFGs at 850 µm to 1.2 mm with unknown redshift, the detection at 2 mm prompted further follow-up as it indicated a much higher probability that the source was likely to sit at $z > 4$. Deep observations from the Atacama Large Millimeter and submillimeter Array (ALMA) presented here confirm the redshift through the secure detection of $^{12}$CO($J = 6 \to 5$) and p-H$_2$O ($2_{1,1} \to 2_{0,2}$). MAMBO-9 is composed of a pair of galaxies separated by 6 kpc with corresponding star formation rates of $590 M_{\odot} \text{yr}^{-1}$ and $220 M_{\odot} \text{yr}^{-1}$, total molecular hydrogen gas mass of $(1.7 \pm 0.4) \times 10^{11} M_{\odot}$, dust mass of $(1.3 \pm 0.3) \times 10^{9} M_{\odot}$, and stellar mass of $(3.2_{-1.5}^{+1.0}) \times 10^{9} M_{\odot}$. The total halo mass, $(3.3 \pm 0.8) \times 10^{12} M_{\odot}$, is predicted to exceed $10^{15} M_{\odot}$ by $z = 0$. The system is undergoing a merger-driven starburst that will increase the stellar mass of the system tenfold in $\tau_{depl} = 40–80$ Myr, converting its large molecular gas reservoir (gas fraction of $96_{-2}^{+1}$%) into stars. MAMBO-9 evaded firm spectroscopic identification for a decade, following a pattern that has emerged for some of the highest-redshift DSFGs found. And yet, the systematic identification of unlensed DSFGs like MAMBO-9 is key to measuring the global contribution of obscured star formation to the star formation rate density at $z \gtrsim 4$, the formation of the first massive galaxies, and the formation of interstellar dust at early times ($\lesssim 1$ Gyr).

Unified Astronomy Thesaurus concepts: Starburst galaxies (1570); Infrared galaxies (790); High-redshift galaxies (734); Blank fields (163)

1. Introduction

The most extreme star-forming galaxies in the universe pose unique challenges for galaxy formation theory (e.g., Fardal et al. 2001; Baugh et al. 2005; Lacey et al. 2008; González et al. 2011; Narayanan 2015). Because dust is a byproduct of star formation, the ubiquity of galaxies with high star formation rate at $z \sim 2$ means that dust-rich systems were common and dominated the cosmic star-forming budget for several billions of years (e.g., Casey et al. 2014; Madau & Dickinson 2014). However, the identification of these dusty star-forming galaxies

(DSFGs) out to higher redshifts ($z \gtrsim 4$), in the first 2 Gyr after the Big Bang, has proven exceedingly difficult. While extraordinary discoveries of DSFGs exist out to $z \sim 7$ (SPT0311 being the highest-$z$ DSFG found to date; Strandet et al. 2017; Marrone et al. 2018), their total contribution to the cosmic star formation budget is unconstrained during this early epoch (Casey et al. 2018b). Contradictory results have been presented in the literature, with some claiming that DSFGs play an insignificant role in $z > 4$ star formation with less than
10% of the total (e.g., Koprowski et al. 2017), while others suggest DSFGs may dominate cosmic star formation at a level exceeding 90% in the first gigayear (Rowan-Robinson et al. 2016). Several other works suggest the truth might lie between these two extremes (e.g., Bethermin et al. 2017; Zavala et al. 2018; Williams et al. 2019), though data to constrain this epoch are sparse, leaving estimates highly uncertain.

Identifying individual DSFGs at early epochs in the
universe's history is critical to our understanding of how
massive galaxies assemble and, independently, how vast dust
reservoirs are formed so early in a galaxy's history, whether it
be from asymptotic giant branch (AGB) stars, supernovae, or
efficient interstellar matter (ISM) grain growth (e.g., Matsuura
et al. 2006, 2009; Zhukovska et al. 2008; Asano et al. 2013;
Jones et al. 2013; Dwek et al. 2014; Lagos et al. 2019).

In this paper, we describe the detection and characterization
of the highest redshift, unlensed DSFG found to date,
confirmed at *z* = 5.85 (see also Jin et al. 2019 for an
independent analysis of this source). This galaxy was identified
as a submillimeter-luminous source by the Max-Planck
Millimeter BOlometer (MAMBO), AzTEC, and SCUBA-2
(Bertoldi et al. 2007; Aretxaga et al. 2011; Casey et al. 2013;
Geach et al. 2017), though it lacked a secure redshift
identification for many years. The source was identified
independently by many groups as a high-redshift candidate
and was recently spectroscopically observed with the Atacama
Large Millimeter/submillimeter Array (ALMA) as presented in
Jin et al. (2019). We corroborate their proposed redshift
solution through independent ALMA observations in this
paper. Here we present a multiwavelength characterization of
the source in order to constrain its physical drivers and
characteristics. Section 2 presents our observations, Section 3
presents calculations of critical physical quantities like
dynamical, gas, stellar, and dust mass, and Section 4 presents
our interpretation of this galaxy's physical drivers and broader
context. Throughout we assume a Planck cosmology (Planck
Collaboration et al. 2018) and assume a Chabrier initial mass
function for the purpose of calculating stellar masses and star
formation rates (Chabrier 2003).

## 2. Data and Observations

2.1. Source Selection and Prior Identification

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 \pm 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\sigma$ 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 \pm 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_{850} = 5.55 \pm 1.11$ mJy and no corresponding 450 µm

counterpart) and further as “COS.0059” in Geach et al. (2017)
with $S_{850} = 5.84 \pm 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 \gtrsim 4$. This is made possible because 2 mm detection is an effective way to “filter out” lower-redshift DSFGs at $1 \lesssim z \lesssim 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 \gtrsim 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\,µm}/S_{2\,mm} = 8.3 \pm 0.9$) implied a high-redshift solution, where a higher value ($S_{850\,µm}/S_{2\,mm} = 15 \pm 3$) would be expected for DSFGs at $z \sim 1-4$. An independent analysis of data from Cycle 5 program 2017.1.00373.S (PI Jin) identified a $4\sigma$ emission line consistent with the measured redshift solution as well as other possible high-$z$ solutions (and corresponding candidate ~$4\sigma$ 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 $^{12}$CO($J=6 \to 5$) and a line at the edge of ALMA band 3, ~84 GHz, as $^{12}$CO($J=5 \to 4$). Our independent analysis of the same data did not lead to a significant detection of the $^{12}$CO($J=5 \to 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.

## 2.2. ALMA Data

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<fnref n="25" /> at $4 \lesssim z \lesssim 9$. The ALMA data were reduced and imaged using the Common Astronomy Software Application<fnref n="26" /> (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

<footnote n="24">²⁴ Aravena (2009) is available at http://hss.ulb.uni-bonn.de/2009/1687/1687.pdf.</footnote>

<footnote n="25">²⁵ These observations were planned before the results of Jin et al. (2019) were known, though we did consider $z = 5.85$ as one of four possible redshift solutions.</footnote>

<footnote n="26">²⁶ http://casa.nrao.edu</footnote>
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 $0''36 \times 0''30$ (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 $1''19 \times 1''09$ (using natural weighting), and an average PWV ranging from 1.3 to 2.3 mm. The rms reached over a $50 \text{ km s}^{-1}$ 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 $1''14 \times 0''93$ (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 \text{ km s}^{-1}$ 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 $0''81 \times 0''68$ (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 $1''83 \times 1''43$ (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).

## 2.3. Ancillary Archival COSMOS Data Sets

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;

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Three-color rendition of the dust continuum emission in MAMBO-9: blue represents 870 µm emission, green is 1.3 mm, and red is 3 mm emission with similar beam sizes. Briggs weighting with robust = 0 is used in all bands with beam sizes of $0''36, $0''75, and $0''65$, respectively. Integer multiples of $\sigma$ above three are shown in contours for the $^{12}$CO($J=6 \to 5$) line emission (at intermediate spatial resolution using Briggs weighting) in context. The northern source is component A and southern component B of MAMBO-9. The $^{12}$CO($J=6 \to 5$) emission in component B is spatially coincident with the 870 µm dust emission. The difference in millimeter color between the two components is real; in other words, component B would have been detected at 3 mm in the dust continuum if it had an SED similar to that of component A.</figcaption></figure>

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 $0''6$ 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 \pm 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 \lesssim 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
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Ten-arcsecond image cutouts of MAMBO-9 from ALMA data sets, including 870 µm, 1.3 mm, 2 mm, and 3 mm continuum (bands 7, 6, 4, and 3), as well as continuum-subtracted moment-0 maps of the $^{12}$CO(J = 6→5) and p-H$_2$O (2$_{1,1}$ → 2$_{0,2}$) lines. At 870 µm, contours follow five times integer powers of two (from 5 to 80σ); all other maps follow odd integer multiples of σ between 3 and 21σ. The peak S/N is 111σ at 870 µm continuum, 20.7σ at 3 mm continuum, 6.7σ at 2 mm, 9.2σ at 1.3 mm, 14.7σ in the $^{12}$CO(J = 6→5) moment-0 map, and 4.3σ in the p-H$_2$O (2$_{1,1}$ → 2$_{0,2}$) moment-0 map. In reference to the 870 µm image, the northern, brighter source is component A and the southern, fainter source is component B. While component A is significantly detected in all maps, component B is only significantly detected (>3σ) at 870 µm and in $^{12}$CO(J = 6→5) (though there are marginal detections at 1.3 and 3 mm).</figcaption></figure>

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 4′′1 to the north and has a photometric redshift of $z = 2.4$ and estimated stellar mass of $4 \times 10^9 M_\odot$. Assuming a halo mass of $4 \times 10^{11} M_\odot$, this would then lead to a maximum value of lensing magnification of $\mu = 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><img src="image_4.png" /><figcaption>**Figure 3.** Aggregate band 3 spectrum of MAMBO-9 from 84 to 110 GHz extracted over both components A and B. The original spectrum of Jin et al. (2019) is shown in gray, offset by 0.7 mJy, with the identification of ¹²CO(J = 5→4) at 84.2 GHz and ¹²CO(J = 6→5) at 101 GHz. Our data are shown as a yellow histogram and have confirmed the detection of the ¹²CO(J = 6→5) line at 101 GHz and detection of the p-H₂O (2₁,₁ → 2₀,₂) line at 109.8 GHz, confirming the redshift as z = 5.850. Vertical red lines mark the expected frequencies of CO and dense gas tracers in the observed frequency range. Inset plots zoom in on the line detections at 101 GHz and 109.8 GHz.</figcaption></figure>

<table><caption>**Table 1**
MAMBO-9 Photometry</caption><thead><tr><th>Band</th><th>Wavelength</th><th>Units</th><th>Component A</th><th>Component B</th><th>Total (A+B)</th><th>Data Reference</th></tr></thead><tbody><tr><td>HST-F606W</td><td>606 nm</td><td>nJy</td><td>(3.7 ± 8.8)</td><td>(-0.5 ± 8.8)</td><td>(10.2 ± 25.7)</td><td>Koekemoer et al. (2011)</td></tr><tr><td>HST-F814W</td><td>814 nm</td><td>nJy</td><td>(9.2 ± 11.5)</td><td>(-2.6 ± 11.5)</td><td>(-3.2 ± 31.7)</td><td>Koekemoer et al. (2011)</td></tr><tr><td>HST-F125W</td><td>1.25 µm</td><td>nJy</td><td>(1.8 ± 14.0)</td><td>(34.3 ± 14.0)</td><td>(36.6 ± 47.0)</td><td>Koekemoer et al. (2011)</td></tr><tr><td>HST-F160W</td><td>1.60 µm</td><td>nJy</td><td>(5.3 ± 13.8)</td><td>(15.4 ± 13.8)</td><td>(50.0 ± 41.0)</td><td>Koekemoer et al. (2011)</td></tr><tr><td>IRAC-CH1</td><td>3.6 µm</td><td>nJy</td><td>...</td><td>...</td><td>87 ± 29</td><td>Ashby et al. (2015)</td></tr><tr><td>IRAC-CH2</td><td>4.5 µm</td><td>nJy</td><td>...</td><td>...</td><td>186 ± 37</td><td>Ashby et al. (2015)</td></tr><tr><td>MIPS24</td><td>24 µm</td><td>µJy</td><td>...</td><td>...</td><td>(10 ± 18)</td><td>Le Floc'h et al. (2009)</td></tr><tr><td>PACS</td><td>100 µm</td><td>µJy</td><td>...</td><td>...</td><td>(48 ± 152)</td><td>Lutz et al. (2011)</td></tr><tr><td>PACS</td><td>160 µm</td><td>µJy</td><td>...</td><td>...</td><td>(-56 ± 276)</td><td>Lutz et al. (2011)</td></tr><tr><td>SPIRE</td><td>250 µm</td><td>mJy</td><td>...</td><td>...</td><td>(2.9 ± 5.8)</td><td>Oliver et al. (2012)</td></tr><tr><td>SPIRE</td><td>350 µm</td><td>mJy</td><td>...</td><td>...</td><td>(2.9 ± 6.3)</td><td>Oliver et al. (2012)</td></tr><tr><td>SCUBA-2</td><td>450 µm</td><td>mJy</td><td>...</td><td>...</td><td>(2.32 ± 5.82)</td><td>Casey et al. (2013)</td></tr><tr><td>SPIRE</td><td>500 µm</td><td>mJy</td><td>...</td><td>...</td><td>(4.9 ± 6.8)</td><td>Oliver et al. (2012)</td></tr><tr><td>SCUBA-2</td><td>850 µm</td><td>mJy</td><td>...</td><td>...</td><td>5.84 ± 0.87</td><td>Geach et al. (2017)</td></tr><tr><td>ALMA-B7</td><td>871 µm</td><td>mJy</td><td>4.032 ± 0.048</td><td>1.486 ± 0.280</td><td>5.908 ± 0.052</td><td>THIS WORK</td></tr><tr><td>ALMA-B7</td><td>876 µm</td><td>mJy</td><td>3.938 ± 0.042</td><td>1.410 ± 0.285</td><td>5.262 ± 0.041</td><td>THIS WORK</td></tr><tr><td>ALMA-B7</td><td>902 µm</td><td>mJy</td><td>3.851 ± 0.050</td><td>1.180 ± 0.246</td><td>5.220 ± 0.049</td><td>THIS WORK</td></tr><tr><td>ALMA-B7</td><td>908 µm</td><td>mJy</td><td>3.853 ± 0.065</td><td>1.650 ± 0.315</td><td>5.666 ± 0.069</td><td>THIS WORK</td></tr><tr><td>AzTEC</td><td>1100 µm</td><td>mJy</td><td>...</td><td>...</td><td>4.6 ± 1.2</td><td>Aretxaga et al. (2011)</td></tr><tr><td>MAMBO</td><td>1200 µm</td><td>mJy</td><td>...</td><td>...</td><td>4.9 ± 0.9</td><td>Bertoldi et al. (2007)</td></tr><tr><td>ALMA-B6</td><td>1287 µm</td><td>mJy</td><td>1.39 ± 0.09</td><td>0.31 ± 0.09</td><td>2.05 ± 0.11</td><td>THIS WORK</td></tr><tr><td>ALMA-B4</td><td>2038 µm</td><td>µJy</td><td>556 ± 83</td><td>(15 ± 83)</td><td>630 ± 74</td><td>THIS WORK</td></tr><tr><td>ALMA-B3</td><td>2880 µm</td><td>µJy</td><td>171.6 ± 7.9</td><td>24.1 ± 7.9</td><td>190.9 ± 8.5</td><td>THIS WORK</td></tr><tr><td>ALMA-B3</td><td>3287 µm</td><td>µJy</td><td>79.8 ± 5.6</td><td>9.0 ± 6.4</td><td>103.5 ± 7.5</td><td>THIS WORK</td></tr><tr><td>VLA-3 GHz</td><td>10 cm</td><td>µJy</td><td>7.34 ± 2.29</td><td>(3.69 ± 2.26)</td><td>10.6 ± 4.1</td><td>Smolčić et al. (2017)</td></tr></tbody></table>

**Note.** Measurements with <3σ significance are enclosed in parentheses, denoting a formal nondetection; this includes measurements that have negative flux density consistent with no detection. Note that optical/near-infrared constraints for each component A and B are measured using a 0.′′6 aperture centered on the ALMA 870 µm resolved components.

### 3. Results

Joint analysis of our ALMA data leads to a spectroscopic confirmation of MAMBO-9 at z = 5.850 through the detection of ¹²CO(J = 6→5) at 14.7σ significance and p-H₂O (2₁,₁ → 2₀,₂) at 4.3σ significance. The CO line is detected in both components A and B, while the p-H₂O (2₁,₁ → 2₀,₂) is only detected in component A. This implies component B only has a single spectral line detection. However, we determine that it is extremely unlikely for component B to sit at a redshift that

is physically unassociated with component A because of their proximity on the sky, similarity in optical/near-infrared photometry, and detection of ¹²CO(J = 6→5) emission. The aggregate photometry for both components is given in Table 1.

Figures 1 and 2 show the ALMA continuum and line moment-0 maps of MAMBO-9 overlaid together and individually. The 870 µm and 3 mm maps use Briggs weighting (robust = 0) to maximize spatial resolution; 1.3 mm, ¹²CO(J = 6→5) moment-0, and 2 mm continuum maps are
shown with a weighting between Briggs and natural (robust = 0.5), and the p-H₂O (2₁,₁ → 2₀,₂) map is shown using natural weighting (robust = 2) to maximize line S/N. Figure 3 shows the aggregate ALMA band 3 spectrum of MAMBO-9 in the range 84–110 GHz. What follows is analysis of the detection of the spectral features ¹²CO(J = 6→5) and p-H₂O (2₁,₁ → 2₀,₂), a discussion of continuum-derived properties, and then the physical characteristics derived from SED fitting.

### 3.1. Millimeter Spectral Line Measurements

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₁,₁ → 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₁,₁ → 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₁,₁ → 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.

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₁,₁ → 2₀,₂) features by using

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Continuum-subtracted source spectra of MAMBO-9 in the ¹²CO (J = 6→5) and p-H₂O (2₁,₁ → 2₀,₂) lines with velocity relative to a central redshift of *z* = 5.850. The integrated spectrum (yellow histogram) is analogous to that shown in Figure 3, that is, the spectrum from naturally weighted band 3 data, but with continuum emission subtracted (the orange line indicates the level of the continuum). The rms per channel is indicated by the gray horizontal stripe. The ¹²CO(J = 6→5) line is then further separated into two components A and B using the highest-spatial-resolution reduction (as shown in Figure 2 using Briggs-weighted data with robust = 0.0). The coaddition of the spectra of components A and B is, within uncertainty, in agreement with the total integrated spectrum from the improved-sensitivity weighting, suggesting very little ¹²CO(J = 6→5) emission is resolved out. The integrated line S/N for ¹²CO(J = 6→5) is 13.3σ in component A and 3.8σ in component B. The p-H₂O (2₁,₁ → 2₀,₂) line is detected at 4.3σ in component A only (shown here are data with natural weighting where the two components are not spatially resolved).</figcaption></figure>

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₁,₁ → 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⁻¹.

### 3.2. p-H₂O (2₁,₁ → 2₀,₂) Emission

The p-H₂O (2₁,₁ → 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) ~ 10⁵–10⁶ cm⁻³), 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₂
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Position–velocity diagram of the $^{12}$CO(J = 6→5) line in the highest-resolution Briggs-weighted (robust = 0.0) data cube overlaid with the 3σ-significance solid black contours of a slightly lower resolution processing (robust = 0.5) of the same data. The image color scale is the same as in the on-sky projection in Figure 2. Both are extracted using a 0.′′7-wide “slit” with orientation position angle of 0°, as shown in the lower left inset plot. Component A spans a spatial extent 0.′′85 ± 0.′′20 (=5.0 ± 1.2 kpc total extent) and velocity $V_{max}$ = 350 ± 50 km s$^{-1}$. The position–velocity kinematics are suggestive of rotation (white dashed line) with $V_{max}$ = 300 km s$^{-1}$, though they do not rule out more complex interaction dynamics at the given spatial resolution. Component B (spatially offset 1′′ to the south of component A) is barely detected in this high-resolution data cube but detected at higher significance at lower resolution and in the moment-0 line map.</figcaption></figure>

and CO in shock-heated regions of the ISM that trace star-forming regions (Bergin et al. 2003). Furthermore, the velocity structure of H$_2$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$_2$O (2$_{1,1}$ → 2$_{0,2}$) feature is (3.6 ± 0.8) × 10$^7$ L$_{⊙}$, precisely on the L$_{IR}$–L$_{H_2O}$ 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$_2$O might be a particularly good star formation tracer. Furthermore, the ratio of line flux between p-H$_2$O (2$_{1,1}$ → 2$_{0,2}$) and $^{12}$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).

### 3.3. Dust Mass

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 $\lambda_{rest} \gtrsim 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 $\lambda_{rest} \approx 300$ µm.<fnref n="27" /> 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:

$$ M_{\text{dust}} = \frac{S_{\nu_{\text{obs}}} D_{L}^{2}(1+z)^{-(3+\beta)}}{\kappa(\nu_{\text{ref}}) B_{\nu}(\nu_{\text{ref}}, T_{\text{dust}})} \left(\frac{\nu_{\text{ref}}}{\nu_{\text{obs}}}\right)^{2+\beta} \left(\frac{\Gamma_{\text{RJ(ref,0)}}}{\Gamma_{\text{RJ}}}\right) \times \left(1 - \frac{B(\nu_{\text{rest}}, T_{\text{CMB}(z)})}{B(\nu_{\text{rest}}, T_{\text{dust}})}\right)^{-1} \quad (1) $$

Here, observations are acquired at $\nu_{\text{obs}}$ (measured in Hz) with flux density $S_{\nu_{\text{obs}}}$ (measured in erg s$^{-1}$ cm$^{-2}$ Hz$^{-1}$), and $\nu_{\text{rest}} = \nu_{\text{obs}}(1+z)$. The frequency at which the dust mass absorption coefficient is known is $\nu_{\text{ref}}$, for example, a value of $\kappa(450~\mu\text{m}) = 1.3 \pm 0.2$ cm$^2$ g$^{-1}$ (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_\nu$ is the Planck function evaluated at a given frequency and for a given temperature in units of erg s$^{-1}$ cm$^{-2}$ Hz$^{-1}$. For example, $B(\nu_{\text{rest}}, T_{\text{CMB}(z)})$ is the Planck function evaluated at $\nu_{\text{rest}}$ for the temperature of the CMB at the measured redshift z. Here, $M_{\text{dust}}$ is in units of g, which can be converted to $M_{\odot}$; $\beta$ is the emissivity spectral index, which we set to $\beta = 1.95$ (we derive this value from a fit to our data in Section 3.7.2); $\Gamma_{\text{RJ}}$ represents the Rayleigh–Jeans (RJ) correction factor, or the deviation from the RJ approximation, following the framework of Scoville et al. (2016):

$$ \Gamma_{\text{RJ}}(\nu, T_{\text{d}}, z) = \frac{h\nu(1+z)/kT_{\text{d}}}{e^{h\nu(1+z)/kT_{\text{d}}} - 1} \quad (2) $$

and $\Gamma_{\text{RJ(ref,0)}} = \Gamma_{\text{RJ}}(\nu = \nu_{\text{ref}}, T_{\text{d}} = T_{\text{dust}}, z = 0)$. Here, h is the Planck constant (6.63×10$^{17}$ erg s); k is the Boltzmann constant (1.38 × 10$^{-16}$ erg K$^{-1}$); $T_{\text{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 $\nu_{\text{rest}} = \nu_{\text{obs}}(1+z)$, the CMB temperature at the given redshift, $T_{\text{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$^9$ M$_{\odot}$ for component A and (1.9$^{+1.3}_{-0.8}$) × 10$^8$ M$_{\odot}$ 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

<footnote n="27">²⁷ If the SEDs were optically thick at these wavelengths, the dust mass would be underestimated using this technique.</footnote>
<table><caption>**Table 2**

Derived Properties of the MAMBO-9 System</caption><thead><tr><th>Derived Property</th><th>Units</th><th>Component A</th><th>Component B</th><th>Total A+B</th></tr></thead><tbody><tr><td>R.A.</td><td>...</td><td>10:00:26.356</td><td>10:00:26.356</td><td>...</td></tr><tr><td>Decl.</td><td>...</td><td>+02:15:27.94</td><td>+02:15:26.63</td><td>...</td></tr><tr><td colspan="5" style="text-align: center; text-decoration: underline;">From ALMA Spectroscopy:</td></tr><tr><td>z</td><td>...</td><td>5.850</td><td>5.852</td><td>5.850</td></tr><tr><td>I<sub>CO(6-5)</sub></td><td>Jy km s<sup>-1</sup></td><td>0.43 ± 0.03</td><td>0.07 ± 0.02</td><td>0.48 ± 0.03</td></tr><tr><td>σ<sub>v</sub>(CO)</td><td>km s<sup>-1</sup></td><td>260 ± 40<sup>a</sup></td><td>280 ± 130</td><td>700 ± 70</td></tr><tr><td>L'<sub>CO(6-5)</sub></td><td>K km s<sup>-1</sup> pc<sup>2</sup></td><td>(1.4 ± 0.9) × 10<sup>10</sup></td><td>(2.3 ± 0.7) × 10<sup>9</sup></td><td>(1.56 ± 0.10) × 10<sup>10</sup></td></tr><tr><td>I<sub>H2O(2,1-20,2)</sub></td><td>Jy km s<sup>-1</sup></td><td>0.09 ± 0.02</td><td>...</td><td>0.09 ± 0.02</td></tr><tr><td>σ<sub>v</sub>(H<sub>2</sub>O)</td><td>km s<sup>-1</sup></td><td>900 ± 200</td><td>...</td><td>900 ± 200</td></tr><tr><td>L'<sub>H2O(2,1-20,2)</sub></td><td>K km s<sup>-1</sup> pc<sup>2</sup></td><td>(2.5 ± 0.5) × 10<sup>9</sup></td><td>...</td><td>(2.5 ± 0.5) × 10<sup>9</sup></td></tr><tr><td>M<sub>gas</sub>(CO)</td><td>M<sub>☉</sub></td><td>(3.3 ± 1.7) × 10<sup>11</sup></td><td>(5 ± 3) × 10<sup>10</sup></td><td>(3.7 ± 1.8) × 10<sup>11</sup></td></tr><tr><td>M<sub>dyn</sub></td><td>M<sub>☉</sub></td><td>(2.0 ± 0.8) × 10<sup>11</sup></td><td>(7 ± 6) × 10<sup>10</sup></td><td>...</td></tr><tr><td colspan="5" style="text-align: center; text-decoration: underline;">From ALMA Dust Continuum:</td></tr><tr><td>M<sub>dust</sub></td><td>M<sub>☉</sub></td><td>(1.3 ± 0.3) × 10<sup>9</sup></td><td>(1.9<sup>+1.3</sup><sub>-0.8</sub>) × 10<sup>8</sup></td><td>(1.6<sup>+0.4</sup><sub>-0.3</sub>) × 10<sup>9</sup></td></tr><tr><td>M<sub>gas</sub> (3 mm)</td><td>M<sub>☉</sub></td><td>(1.4 ± 0.4) × 10<sup>11</sup></td><td>(1.2 ± 0.5) × 10<sup>10</sup></td><td>(1.7 ± 0.4) × 10<sup>11</sup></td></tr><tr><td>FWHM<sub>maj</sub>(870 µm)<sup>b</sup></td><td>"</td><td>0<sup>"</sup>15 ± 0<sup>"</sup>01</td><td>0<sup>"</sup>30 ± 0<sup>"</sup>05</td><td>...</td></tr><tr><td>Axis Ratio (870 µm)<sup>b</sup> b/a</td><td>...</td><td>0.87 ± 0.15</td><td>0.37 ± 0.24</td><td>...</td></tr><tr><td>R<sub>eff</sub>(870 µm)</td><td>pc</td><td>380 ± 30</td><td>760 ± 130</td><td>...</td></tr><tr><td>FWHM<sub>maj</sub>(3 mm)<sup>b</sup></td><td>"</td><td>0<sup>"</sup>29 ± 0<sup>"</sup>11</td><td>...</td><td>...</td></tr><tr><td>R<sub>eff</sub>(3 mm)</td><td>pc</td><td>700 ± 300</td><td>...</td><td>...</td></tr><tr><td>Σ<sub>Mdust</sub></td><td>M<sub>☉</sub> pc<sup>-2</sup></td><td>1400 ± 400</td><td>50 ± 30</td><td>...</td></tr><tr><td>Σ<sub>SFR</sub></td><td>M<sub>☉</sub> yr<sup>-1</sup> kpc<sup>-2</sup></td><td>640 ± 170</td><td>61 ± 35</td><td>...</td></tr><tr><td colspan="5" style="text-align: center; text-decoration: underline;">From Broadband SED Fitting:</td></tr><tr><td>λ<sub>rest</sub> at which τ<sub>ν</sub> = 1</td><td>µm</td><td>≡200</td><td>Opt. Thin</td><td>...</td></tr><tr><td>L<sub>IR</sub></td><td>L<sub>☉</sub></td><td>(4.0<sup>+0.9</sup><sub>-0.7</sub>) × 10<sup>12</sup></td><td>(1.5<sup>+1.1</sup><sub>-0.5</sub>) × 10<sup>12</sup></td><td>(6.3<sup>+1.1</sup><sub>-0.9</sub>) × 10<sup>12</sup></td></tr><tr><td>SFR</td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>590<sup>+140</sup><sub>-100</sub></td><td>220<sup>+150</sup><sub>-70</sub></td><td>930<sup>+160</sup><sub>-130</sub></td></tr><tr><td>λ<sub>peak</sub></td><td>µm</td><td>87 ± 7</td><td>100<sup>+30</sup><sub>-80</sub></td><td>...</td></tr><tr><td>T<sub>dust</sub></td><td>K</td><td>56.3<sup>+5.9</sup><sub>-5.7</sub></td><td>29.7<sup>+8.5</sup><sub>-6.6</sub></td><td>...</td></tr><tr><td>β</td><td>...</td><td>1.95 ± 0.11</td><td>2.66<sup>+0.22</sup><sub>-0.34</sub></td><td>...</td></tr><tr><td>M<sub>*</sub> (OIR-ONLY)</td><td>M<sub>☉</sub></td><td>...</td><td>...</td><td>(3.2<sup>+1.0</sup><sub>-1.5</sub>) × 10<sup>9</sup></td></tr><tr><td>L<sub>UV</sub>(1600 Å)</td><td>L<sub>☉</sub></td><td>&lt;7.7 × 10<sup>9</sup></td><td>(8.8 ± 3.6) × 10<sup>9</sup></td><td>&lt;3.4 × 10<sup>10</sup></td></tr><tr><td>IRX</td><td>...</td><td>&gt;510</td><td>160<sup>+280</sup><sub>-100</sub></td><td>...</td></tr><tr><td>A<sub>UV</sub></td><td>...</td><td>&gt;6.2</td><td>5.0<sup>+1.0</sup><sub>-1.1</sub></td><td>...</td></tr><tr><td>q<sub>IR</sub></td><td>...</td><td>0.4 ± 1.0</td><td>...</td><td>...</td></tr><tr><td>M<sub>halo</sub></td><td>M<sub>☉</sub></td><td>...</td><td>...</td><td>(3.3 ± 0.8) × 10<sup>12</sup></td></tr></tbody></table>

**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 $^{12}$CO($J=6 \to 5$) line, $\sigma_v(\text{CO})$ is the velocity dispersion of the $^{12}$CO($J=6 \to 5$) feature, and $L'_{CO(6-5)}$ is the $^{12}$CO($J=6 \to 5$) line luminosity. All three quantities are similarly derived for the p-H$_2$O ($2,1 \to 2,0$) line. $M_{\text{gas}}(\text{CO})$ is the gas mass as derived from the $^{12}$CO($J=6 \to 5$) line, while $M_{\text{gas}}(3 \text{ mm})$ is the gas mass as derived from the 3 mm dust continuum (and $M_{\text{dust}}$ is the dust mass derived from the 3 mm continuum). Both assume $\alpha_{\text{CO}} = 6.5 M_{\odot}$ (K km s$^{-1}$ pc$^{-2}$)$^{-1}$. $M_{\text{dyn}}$ is the dynamical mass as estimated from the $^{12}$CO($J=6 \to 5$) line width and 870 µm dust continuum size. FWHM$_{\text{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_{\text{eff}}$ is the circularized half-light radius in parsecs. $\Sigma_{M_{\text{dust}}}$ and $\Sigma_{SFR}$ are the dust mass surface density and star formation surface density, respectively. $\lambda_{\text{rest}}$ is the wavelength at which $\tau = 1$ for the dust SED, $L_{\text{IR}}$ is the derived IR luminosity integrated from 8 to 1000 µm, and SFR is the star formation rate converted directly from $L_{\text{IR}}$ using the Kennicutt & Evans (2012) scaling. $\lambda_{\text{peak}}$ is the rest-frame wavelength where the dust SED peaks, and $T_{\text{dust}}$ is the underlying dust temperature used in the model fit to the photometry. $\beta$ is the emissivity spectral index. $M_*$ is the stellar mass of the aggregate system, while $L_{\text{UV}}$ is the rest-frame 1600 Å luminosity. IRX is the ratio $L_{\text{IR}}/L_{\text{UV}}$, $A_{\text{UV}}$ is the absolute magnitude of attenuation inferred at 1600 Å, $q_{\text{IR}}$ is the implied far-infrared (FIR)-to-radio ratio as in Yun et al. (2001), and $M_{\text{halo}}$ is the total inferred halo mass.

attributed to the differences in SED fitting including the best-fit $\beta$ (our dust mass uncertainties do not account for the measurement uncertainties in $\beta$).

### 3.4. Gas Mass

We derive the mass of molecular gas in each galaxy from the dust continuum as well as from $^{12}$CO($J=6 \to 5$) line luminosity. In both cases, we adopt a value for the CO-to-H$_2$

conversion factor of $\alpha_{\text{CO}} = 6.5 M_{\odot}$ (K km s$^{-1}$ pc$^{-2}$)$^{-1}$ as in Scoville et al. (2016). This is in line with the Galactic value and accounts for the mass of both H$_2$ and He gas.$^{28}$

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

<footnote><sup>a</sup> Component A is best characterized by double-peaked emission for which the line width of each component has width of 370 km s$^{-1}$.</footnote>

<footnote><sup>b</sup> FWHM of the major axis measured from a two-dimensional Gaussian fit in the image plane, and the axis ratio is from the same fit.</footnote>

<footnote>$^{28}$ Without including He, the Galactic value is $\sim 4.5 M_{\odot}$ (K km s$^{-1}$ pc$^{-2}$)$^{-1}$.</footnote>
effect on the dust mass calculation:

$$M_{\text{gas}} = \frac{4\pi D_{\text{L}}^2 S_{\nu,\text{obs}}}{\alpha(\nu_{850})(1+z)^{3+\beta}} \left(\frac{\nu_{850}}{\nu_{\text{obs}}}\right)^{2+\beta} \left(\frac{\Gamma_{\text{RJ(ref,0)}}}{\Gamma_{\text{RJ}}}\right) \times \left(1 - \frac{B(\nu_{\text{rest}}, T_{\text{CMB}}(z))}{B(\nu_{\text{rest}}, T_{\text{dust}})}\right)^{-1} \quad (3)$$

Here, $\alpha_{850} = (6.7 \pm 1.7) \times 10^{19}$ erg s$^{-1}$ Hz$^{-1}$ M$_{\odot}^{-1}$ 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$_2$ 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_{\text{gas}}$ as given in Equation (3) is in units of M$_{\odot}$. With this approach, we constrain masses of molecular gas using the dust continuum to be $M_{\text{gas}} = (1.4 \pm 0.4) \times 10^{11} M_{\odot}$ for component A and $M_{\text{gas}} = (1.2 \pm 0.5) \times 10^{10} M_{\odot}$ 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 $^{12}\text{CO}(J=6 \to 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 $^{12}\text{CO}(J=6 \to 5)$ as an independent check against what we have calculated using the dust continuum.

We use the $^{12}\text{CO}(J=6 \to 5)$ line luminosity$^{29}$ to derive a molecular gas mass from $^{12}\text{CO}(J=6 \to 5)$. This requires an assumption as to the value of the gas excitation spectral line energy distribution (SLED) to convert from $^{12}\text{CO}(J=6 \to 5)$ to the ground-state $^{12}\text{CO}(J=1 \to 0)$ and then the value of the CO-to-H$_2$ 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\sigma$ uncertainty on $I_{\text{CO}(6-5)}/I_{\text{CO}(1-0)} = 10^{+30}_{-5}$). We calculate gas masses of $M_{\text{gas}}(A) = (3.3 \pm 1.7) \times 10^{11} M_{\odot}$ and $M_{\text{gas}}(B) = (5 \pm 3) \times 10^{10} M_{\odot}$ 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_{\text{CO}} = 6.5 M_{\odot}$ (K km s$^{-1}$ pc$^{-2}$)$^{-1}$. Both calculations of the gas masses account for measurement uncertainties in flux density, $\alpha_{850}$ or $I_{\text{CO}(6-5)}/I_{\text{CO}(1-0)}$, but not uncertainty in $\alpha_{\text{CO}}$. If we were to instead take $\alpha_{\text{CO}} = 1 M_{\odot}$ (K km s$^{-1}$ pc$^{-2}$)$^{-1}$, 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_{\text{gas}} = (2.2 \pm 0.6) \times 10^{10} M_{\odot}$ for component A and $M_{\text{gas}} = (1.8 \pm 0.8) \times 10^{9} M_{\odot}$ for component B

(both scaled down from the dust-continuum-estimated gas masses).

### 3.5. Size Measurements

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_c(A) = 0''068 \pm 0''002 = 408 \pm 12$ pc and $R_c(B) = 0''110 \pm 0''010 = 660 \pm 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_c(A) = 0''064 \pm 0''004 = 380 \pm 30$ pc and $R_c(B) = 0''128 \pm 0''021 = 760 \pm 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_c(A) = 0''123 \pm 0''047 = 700 \pm 300$ pc, only $1\sigma$ 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_{\text{eff}}$ sizes due to the high S/N near the peak of dust emission and the measured dust masses to calculate $\Sigma_{M_{\text{dust}}}(A) = 1400 \pm 400 M_{\odot}$ pc$^{-2}$ and $\Sigma_{M_{\text{dust}}}(B) = 50 \pm 30 M_{\odot}$ pc$^{-2}$ 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 $\Sigma_{\text{SFR}}(A) = 640 \pm 170 M_{\odot}$ yr$^{-1}$ kpc$^{-2}$ and $60 \pm 35 M_{\odot}$ yr$^{-1}$ kpc$^{-2}$. Neither is near the hypothetical Eddington limit for starbursts that are factors of

<footnote>$^{29}$ This uses the standard L'$_{\text{CO}}$ definition as in Equation (3) of Solomon & Vanden Bout (2005), the first equation of Carilli & Walter (2013), and Equation (19) of Casey et al. (2014).</footnote>
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.

### 3.6. Dynamical Mass

We derive the dynamical masses of MAMBO-9-A and MAMBO-9-B using the $^{12}$CO($J=6 \to 5$) kinematic profile, comparing a few different methods. First, for a galaxy with an unresolved velocity field, the dynamical mass is best estimated by

$$M_{\text{dyn}} = \frac{5R_e\sigma^2}{G}, \qquad (4)$$

where G is the gravitational constant, σ is the measured velocity dispersion of the kinematic feature measured (in our case, $^{12}$CO($J=6 \to 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^2i$).

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 ($^{12}$CO($J=6 \to 5$) in this case). Second, using a high-J tracer like $^{12}$CO($J=6 \to 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 \pm 30$ pc and $R_e(B) = 760 \pm 130$ pc. The velocity dispersions as measured from $^{12}$CO($J=6 \to 5$) (and as shown in Figure 4) are $\sigma_V(A) = 260 \pm 40$ km s$^{-1}$ and $\sigma_V(B) = 280 \pm 130$ km s$^{-1}$. This gives dynamical mass estimates of $M_{\text{dyn}}(A) = (5 \pm 2) \times 10^{10} M_{\odot}$ and $M_{\text{dyn}}(B) = (7 \pm 6) \times 10^{10} M_{\odot}$, 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 $^{12}$CO($J=6 \to 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$_2$O ($2_{1,1} \to 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$_2$O ($2_{1,1} \to 2_{0,2}$) line is

much more broad than the $^{12}$CO($J=6 \to 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 $^{12}$CO($J=6 \to 5$) kinematics using

$$M_{\text{dyn}} = \frac{V_{\text{max}}^2 R_{\text{max}}}{G}, \qquad (5)$$

which then similarly needs to be corrected for unknown inclination. Note that $R_{\text{max}}$ here denotes the maximum radius at which $V_{\text{max}}$ is measured and differs from $R_e$, the circularized half-light radius, used above. Using both size $R_{\text{max}} = 0''85 \pm 0''20 = 5.0 \pm 1.2$ kpc and $V_{\text{max}} = 350 \pm 50$ km s$^{-1}$ measurements from the position-velocity diagram in Figure 5, we derive an alternate dynamical mass of component A of $(2.0 \pm 0.8) \times 10^{11} M_{\odot}$. 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 \pm 0.8) \times 10^{11} M_{\odot}$, 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_{\text{dyn}}$ of $(7 \pm 6) \times 10^{10} M_{\odot}$.

### 3.7. SED Fitting

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.

#### 3.7.1. Optical/Near-infrared SED Fit

We explore what constraints can be set using the optical and mid-infrared photometry alone. As the stellar component is detected only in the deepest near-infrared imaging, we use only the HST imaging from COSMOS (Scoville et al. 2007) and CANDELS (Grogin et al. 2011; Koekemoer et al. 2011) data, in addition to S-CANDELS Spitzer IRAC measurements (Ashby et al. 2015) for the optical/near-infrared fit; all other existing data are not deep enough to provide meaningful constraints.

We measure photometry using Source Extractor (Bertin & Arnouts 1996), using a combined [3.6]+[4.5] image as the detection image. We optimize the detection parameters such that the isophotal region corresponds to an ellipse that includes the majority of the bright IRAC emission and is tuned to enclose both ALMA 870 µm continuum peaks (see Figure 2). We up-sample the IRAC images to the same 0''06 pixel scale as the HST photometry (altering the zero-point appropriately), and we run Source Extractor with the HST F606W, F814W, F125W, F140W, and F160W images as the measurement
<figure><img src="image_7.png" /><figcaption>**Figure 6.** At top, 6" × 6" cutouts of MAMBO-9 in two HST bands, the IRAC bands, ALMA 870 µm, and VLA 3 GHz. Contours in each frame denote the 5σ significance contours on the 870 µm image (also shown in Figure 2); the white dotted line shows the aperture (based on IRAC emission) used to measure photometry in the optical and near-infrared. The only significant (>3σ) detections come from *Spitzer* IRAC, ALMA, and VLA at 3 GHz. Below, the aggregate composite SED for both components A+B is shown in black, made of three primary components: the stellar and nebular line emission (dark blue), the thermal dust emission (orange), and synchrotron radio emission (purple). All three components are independently fit to data in their respective regimes. The light blue curve shows the modeled unattenuated stellar and nebular emission (described in the text); the dotted orange line shows the dust SED fit to component A, the dashed orange line is for component B, and the gray line shows the best-fit MAGPHYS SED. Our SED does not include emission from polycyclic aromatic hydrocarbons in the mid-infrared due to the existing dearth of data in that regime. Upper limits are shown as 2σ. The inset plot shows the probability distributions of stellar mass derived for MAMBO-9 from the Optical InfraRed (OIR)-only fit (blue) and MAGPHYS fit (gray).</figcaption></figure>

images. Though this area is covered by shallow F140W from the 3D-HST Survey, MAMBO-9 falls in a coverage gap (Momcheva et al. 2016). As expected, we find a significant detection in both IRAC bands, with no significant flux measured in any of the HST bands.

We use this isophotal photometry to estimate the stellar population modeling parameters following the methodology of Finkelstein et al. (2015). In brief, we use the EAZY software (Brammer et al. 2008) to fit the photometry using the updated templates derived from Flexible Stellar Population Synthesis (FSPS) models (Conroy & Wechsler 2009; Conroy et al. 2010; see a forthcoming paper by S. Finkelstein for more details on the templates). In the absence of a spectroscopic identification, such little photometric information would result in a photo-z probability distribution that is very broad, consistent with z > 2. We then performed SED fitting using the spectroscopic redshift, performing $\chi^2$ minimization of a set of Bruzual & Charlot (2003) stellar population models, with added nebular emission and dust attenuation. The results are illustrated in the optical portion of the full SED shown in Figure 6 (blue lines). The inset plot shows the inferred distribution of stellar mass for the best-fit “OIR” SEDs (in blue), with a median of $(3.2_{-1.5}^{+1.0}) \times 10^9 M_\odot$. From the limited optical/near-infrared data alone, the absolute magnitudes of attenuation estimated in the rest-frame V band are $A_V = 3.1 \pm 0.2$ and SFR = $63.4_{-6.8}^{+6.5} M_\odot \text{yr}^{-1}$ (corrected for “dust” that is estimated from the OIR fit). Both

are underestimated relative to the measured characteristics of the far-infrared SED.

The inferred UV luminosity from the best-fit SED is extrapolated to be $L_{1600} \approx 7 \times 10^8 L_\odot$, though observationally it is only strictly limited to $L_{1800} \lesssim 7 \times 10^{10} L_\odot$ (at 2σ, based on the F125W nondetection). To set more stringent constraints for the individual components A and B, we extract 0.6 circular aperture photometry on the HST data. While component A is not detected, there is a 2σ marginal detection in component B. We use these measured UV constraints to also constrain IRX, defined as the ratio $L_{IR}/L_{UV}$, and the absolute magnitudes of attenuation in the UV, $A_{UV}$, in the samples. These measurements use the $L_{IR}$ calculated in the next section, Section 3.7.2, and given in Table 2. For component A, we measure IRX > 510 and $A_{UV} > 6.2$, while for component B we measure IRX = $160_{-100}^{+280}$ and $A_{UV} = 5.0_{-1.1}^{+1.0}$. Even though the difference between the two components may be substantial, both constitute extremely obscured systems.

### 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 $\alpha$, as described in Casey 2012). Here we fix the value of $\alpha$ to $\alpha_{MIR} = 4$. Physically, a lower value of $\alpha$ 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 $\alpha$ asymptote to the pure modified blackbody solution. While no direct constraints can be made for $\alpha_{MIR}$, it should be noted that values less than $\alpha_{MIR} = 4$ violate the upper limits in the mid-infrared as measured by Spitzer and Herschel. If $\alpha_{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 $\lambda_{peak}$, the rest-frame peak wavelength, $L_{IR}$, the integrated 8–1000 µm IR luminosity, and $\beta$, 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 $\lambda_{peak}$, $L_{IR}$, and $\beta$ (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 $\tau = 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 $\beta$, the emissivity spectral index. For component B, we set an upper limit to $\beta = 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 $\lambda_{peak}$, $L_{IR}$, and $\beta$. 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 $\lambda_{peak}$ values in context against (a) the $z \sim 1–2$ $L_{IR}-\lambda_{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 \sim 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 \gtrsim 5$ should be warmer than those at $z \sim 2$, we do not see compelling evidence that this holds for either component of MAMBO-9.

The parameter $\lambda_{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 $\lambda_{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 $\lambda_{peak}$ in Figure 7 between opacity models is purely due to the different levels of effects of the CMB.

While the CMB does affect $\lambda_{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 $\Sigma_{M_{dust}} = 1400 \pm 400 M_{\odot} \text{ pc}^{-2}$ in component A suggests an optically thick SED to ~200–300 µm rest frame, while the lower dust column density in component B of $\Sigma_{M_{dust}} = 50 \pm 30 M_{\odot} \text{ pc}^{-2}$ 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 $\tau = 1$ at $\lambda_{rest} = 200 \text{ µ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_{\odot} \text{ yr}^{-1}$ for component A and $220_{-70}^{+150} M_{\odot} \text{ yr}^{-1}$ for component B. Note that the total SFR fit to the aggregate SED ($930_{-130}^{+160} M_{\odot} \text{ yr}^{-1}$) is lower than, though not fully inconsistent with, the derived SFR of the system in Jin et al. (2019) of $1283 \pm 173 M_{\odot} \text{ yr}^{-1}$. 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 $\lambda_{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 $\lambda_{peak}$) than the expected median $L_{IR}-\lambda_{peak}$ relation from simulations at $z \sim 5.9$ (Ma et al. 2019).

There is a slight discrepancy between the measured emissivity spectral index, $\beta$, of components A and B of $\beta(A) = 1.95 \pm 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 $\beta$ and $\lambda_{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 $\beta$ 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
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Far-infrared SED fitting details for both components of MAMBO-9, adopting two different model assumptions: optically thin dust (purple), and a more general opacity model (orange) that asserts $\tau = 1$ at a rest frame of 1.5 THz (200 µm). Corner plots are shown for the converged MCMC chains at left constraining $L_{IR}$, $\lambda_{peak}$, and $\beta$ for each component. The upper middle panels show $L_{IR}-\lambda_{peak}$ for each source in context against expectation from cosmological simulations (green dashed line; Ma et al. 2019), as well as against the average $L_{IR}-\lambda_{peak}$ derived for lower-redshift DSFGs (gray band; Casey et al. 2018b). DSFGs at $z > 4$ from the SPT survey are shown as navy squares (Strandet et al. 2016). The lack of detections shortward of the peak imply that the mid-infrared power law cannot be directly constrained, and here we fix $\alpha_{MIR} = 4$, which minimally affects the measured $L_{IR}$. Detections above $5\sigma$ significance are shown in black, $2.5\sigma < S/N < 5\sigma$ detections in gray, and $2\sigma$ and $3\sigma$ upper limits as light and dark gray arrows. The “super-deblended” Herschel photometry shown in Jin et al. (2019) is shown in light open squares but are not used for these fits. The best-fit MAGPHYS fit is shown as a thick blue line, while random draws from the accepted MCMC trials are shown in either light orange (general opacity) or purple (optically thin), with the median-value SEDs shown in dark orange or purple. At $z = 5.85$, the effect of the CMB on the SED is model dependent: the intrinsic emitted SEDs we would expect in the absence of the CMB are shown as dashed purple and orange lines, measurably different only for the optically thin model.</figcaption></figure>

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 $\beta = 2.89 \pm 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 $\beta$ consistent with the often-used assumption in the literature of $\beta = 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}) \times 10^{12} L_{\odot}$, SFR = $590_{-100}^{+140} M_{\odot} yr^{-1}$, rest-frame peak wavelength of $\lambda_{peak} = 87 \pm 7$ µm, dust temperature $T_{dust} = 56.3_{-5.7}^{+5.9}$ K, and emissivity spectral index $\beta = 1.95 \pm 0.11$. We infer that component B is optically thin at the peak, has a measured $L_{IR} = (1.5_{-0.5}^{+1.1}) \times 10^{12} L_{\odot}$, SFR = $220_{-70}^{+150} M_{\odot} yr^{-1}$, rest-frame peak wavelength $\lambda_{peak} = 100_{-80}^{+30}$ µm, dust temperature $T_{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

We employ the updated library of star formation histories described in da Cunha et al. (2015) and stellar population synthesis models from Bruzual & Charlot (2003) to fit the full SED (both components together) using MAGPHYS. The advantage of MAGPHYS comes through the use of energy balance, whereby energy absorbed in the rest-frame UV and optical is reradiated by dust at long wavelengths. However, this technique can break down for sources whose stellar emission is fully decoupled from long-wavelength emission, as is often the case in DSFGs (e.g., Casey et al. 2017).

We fit both components of MAMBO-9 (A+B) jointly using MAGPHYS because of the difficulty of differentiating IRAC fluxes. The best-fit MAGPHYS solution is shown in Figure 6 as well as Figure 7 for comparison against our adopted best-fit OIR-only and obscured-only SED. The MAGPHYS fit does well fitting the total IR photometry and upper limits in the mid-infrared and, likewise, provides a sensible solution to the limited OIR photometry. The probability distribution of stellar mass as derived by MAGPHYS is shown as the inset plot on Figure 6 (in gray), with $M_* = (2.1 \pm 0.7) \times 10^{11} M_\odot$. The predicted $A_V = 5.64_{-0.20}^{+0.02}$ is well aligned with the measured constraint on $A_{UV}$.

There is a marked difference between this stellar mass estimate and that provided by the analysis of the fit to the OIR component only, roughly a factor of ~60× different, with estimated uncertainties far smaller than the relative offset. This can largely be attributed to the lack of constraint near the rest-frame 1.6 µm stellar bump (redshifted to ~11 µm). The SEDs in Figure 6 show a stark difference near 10 µm, which results in these disparate predictions. Whether or not nebular emission lines are included in the stellar population model also affects the results but to a lesser degree. Several works have shown that high equivalent-width emission lines can enhance broad-band flux in the IRAC channels by a factor of ~2–3× (Shim et al. 2011; Salmon et al. 2015). The OIR-only SED fitting from Section 3.7.1 included them, which would potentially drive the stellar mass estimate lower for a given set of photometry, and the stellar population models used by MAGPHYS do not include them. At this redshift, $H\alpha$ 6563 Å emission falls into the [4.5] IRAC band, while O III 5007 Å falls into the [3.6] IRAC band, and sources with significant star formation and ionization radiation will be proportionally brighter in these bands as a result (e.g., Chary et al. 2005; Smit et al. 2014), if those lines are not substantially obscured by dust themselves (which is often the case for DSFGs; Swinbank et al. 2004; Casey et al. 2017).

A detailed analysis of the stellar masses of submillimeter galaxies was performed by Michałowski et al. (2014), who used synthetic photometry of simulated DSFGs to test different literature tools used to infer their stellar mass (among other characteristics). They found that tools like MAGPHYS that model two independent star-formation histories (continuous and starburst) best reproduce the simulated DSFG stellar masses, while single-component star-formation history models tend to underestimate the stellar mass substantially. A comparison of the two stellar mass estimates to other measured quantities in the system, like gas mass and dynamical mass, also suggests that the MAGPHYS-predicted stellar mass may be the more likely of the two (i.e., roughly equal stellar masses

and gas mass, rather than a gas mass that exceeds the stellar mass by ~60×).

After an analysis of the galaxy's total mass budget and predicted halo mass, given in Section 4, we find that the stellar mass is more likely consistent with the OIR-only value, though formally unconstrained, with the possibility of being within the very large range of a few $\times 10^9 M_\odot$ to a few $\times 10^{11} M_\odot$.

### 3.8. Radio Continuum Emission

In the radio regime, MAMBO-9 component A is marginally detected at $3.2\sigma$ 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 $\alpha_{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 \pm 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 $\alpha = -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.

## 4. Discussion

### 4.1. Mass Budget and Halo Mass Rarity

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 $\alpha_{CO} = 6.5 M_\odot (\text{K km s}^{-1} \text{pc}^2)^{-1}$, and the most massive stellar mass uses that derived using MAGPHYS. The least massive set of assumptions uses the more modest value of $\alpha_{CO} = 1 M_\odot (\text{K km s}^{-1} \text{pc}^2)^{-1}$ 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 $\kappa_{\nu}$, the dust mass
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Mass budget of the MAMBO-9 system, considering both the most massive assumptions (left), the least massive assumptions (middle), and the adopted masses (right). This violin plot shows the probability density distribution on its side for each measured variable. The most massive assumptions are derived using $\alpha_{CO} = 6.5 M_{\odot} (\text{K km s}^{-1} \text{pc}^2)^{-1}$ for the gas mass, the MAGPHYS-derived stellar mass, and the kinematically derived dynamical mass from Figure 5. The least massive assumptions are derived using $\alpha_{CO} = 1 M_{\odot} (\text{K km s}^{-1} \text{pc}^2)^{-1}$ for the gas mass, the OIR-only derived stellar mass, and the unresolved dynamical mass estimate. The dust mass is the same in all cases. The halo mass is extrapolated from the total baryonic mass assuming a conservative 1:20 baryonic-to-halo mass ratio; the lower limit on halo mass is the dotted red line, set by the cosmic baryon fraction. The gray shaded regions represent the $1\sigma$, $2\sigma$, and $3\sigma$ exclusion curves for finding a galaxy above a given mass at $z = 5.85$ in our 2 mm blank-field survey. The most massive assumptions predict a halo mass that is very unlikely to be detected by our survey ($<0.5\%$). The least massive assumptions produce a mass that is well within the bounds of our survey but has a highly unusual gas-to-dust ratio. The adopted mass—dominated by molecular gas and dark matter—predicts that MAMBO-9 is 87% likely to be the most massive galaxy in our 2 mm survey.</figcaption></figure>

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.

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_{\text{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_{\text{halo}} = (6.0 \pm 1.4) \times 10^{11} M_{\odot}$ (least massive assumptions) to $M_{\text{halo}} = (7.2 \pm 1.8) \times 10^{12} M_{\odot}$ (most massive assumptions).

The gray regions in Figure 8 denote the $1\sigma$, $2\sigma$, and $3\sigma$ 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

$\sim 155 \text{ arcmin}^2$, with $135 \text{ arcmin}^2$ at the necessary depth required to have identified MAMBO-9 at $5\sigma$ significance or greater.<fnref n="30" /> The exclusion curves represent the rarest, most massive halos that should exist in the $135 \text{ arcmin}^2$ 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 \times 10^{12} M_{\odot}$, $4.3 \times 10^{12} M_{\odot}$, and $8.0 \times 10^{12} M_{\odot}$ for $1\sigma$, $2\sigma$, and $3\sigma$, respectively.

The most massive assumptions provide an estimated halo mass that exceeds the maximum detectable halo mass by $2.8\sigma$; 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 \pm 1.4) \times 10^{12} M_{\odot}$ 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 \pm 0.8) \times 10^{12} M_{\odot}$, 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 $\sim 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 $\alpha_{\text{CO}} = 6.5 M_{\odot} (\text{K km s}^{-1} \text{pc}^2)^{-1}$, the low stellar mass $(3.2_{-1.5}^{+1.0}) \times 10^9 M_{\odot}$, and halo mass of $(3.3 \pm 0.8) \times 10^{12} M_{\odot}$, as reflected in Table 2.

The implied gas-to-dust ratio (GDR $= M_{\text{gas}}/M_{\text{dust}}$) for the system is 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_{\text{gas}} \equiv M_{\text{gas}}/(M_{\text{gas}} + M_{\star} + M_{\text{dust}})$) is extraordinarily high at $f_{\text{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 \sim 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$_2$ conversion factor holds for MAMBO-9?

<footnote n="30">³⁰ Note that the full volume of the 2 mm map to the depth where MAMBO-9 is detectable is $135 \text{ arcmin}^2$, but MAMBO-9 was originally identified in a much smaller survey area, $9.4 \text{ arcmin}^2$, constituting the data delivered in one Scheduling Block. The rest of the 2 mm map was released to us slowly over the course of the preparation of this manuscript.</footnote>
## 4.2. Implications of MAMBO-9 on the Prevalence of High-z DSFGs

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\sigma_{rms} \lesssim 0.12$ mJy/beam at 2 mm,<fnref n="31" /> 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 \lesssim 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 \gtrsim 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 \sim 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 \sim 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 \sim 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.

## 4.3. Physical Drivers of the MAMBO-9 System

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_{depl} = 38_{-12}^{+16}$ Myr, while component B will deplete its gas in $\tau_{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_{halo} = 1.6 \times 10^{15} M_{\odot}$ 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.

## 5. Conclusions

The ALMA data presented herein provide a uniquely detailed snapshot of MAMBO-9, the fourth-highest-redshift DSFG and highest-redshift unlensed DSFG to date, 1 Gyr after the Big Bang. This system is composed of (at least) two galaxies that are separated by 6 kpc likely undergoing a merger or interaction.

The northern component, MAMBO-9-A, is forming stars at nearly ≈600 $M_{\odot}$ yr⁻¹ with an incredibly dense ISM, optically thick to the peak of dust emission near ~200 µm. The southern source, MAMBO-9-B, is less extreme (in both star formation surface density and dust column density) than its northern cousin but may be of similar mass and size. Both components have very high attenuation in the rest-frame UV, with $A_{UV} > 6.2$ and $A_{UV} = 5.0_{-1.1}^{+1.0}$, respectively. MAMBO-9-A also shows some hint at rotational motion in CO kinematics, though still potentially dominated by dispersion (higher resolution observations are needed for more firm kinematic diagnostics). The measured gas depletion time for the system is less than 100 Myr, signaling the short-lived period of rapid stellar growth as is often found with high-z DSFGs.

The system's total halo mass is estimated to be $M_{halo} = (3.3 \pm 0.8) \times 10^{12} M_{\odot}$, which is 87% likely to be the most massive galaxy detectable in our 2 mm blank-field survey. This presumes that the baryonic content of MAMBO-9 is dominated by gas (96$_{-2}^{+1}$%) and that the stellar mass is low, but likely to grow by an order of magnitude in a short period of 40–80 Myr. Despite the multiple data sets providing keen insight into the physical nature of MAMBO-9, underlying systematic uncertainties exist: the unconstrained CO-to-H₂ conversion factor, and the unconstrained stellar mass chief among them, which could significantly shift our inferred halo mass, gas mass, and gas depletion time. Further analysis of the larger 2 mm-selected sample, in addition to future JWST constraints on the stellar emission (i.e., stellar mass and metallicity), will be valuable for understanding the nature of the ISM in MAMBO-9 and its relative rarity.

We have presented a detailed fit to the long-wavelength photometry in order to derive characteristics of dust emission at an epoch where the CMB temperature is nonnegligible; the objective of this analysis was to illustrate the effect of different choices on the derived results, primarily the assumed opacity of the dust. This is most important for systems that lack photometric constraints shortward of the dust emission peak (at rest-frame $\lambda_{peak} \approx 100$ µm), which will be most if not all galaxies identified at such epochs, with the exception of already-identified gravitationally lensed DSFGs detected by Herschel. The only future telescope capable of constraining this regime to sufficient sensitivity, providing direct insight into the

<footnote n="31">³¹ Upon later delivery of more 2 mm data, this rms was pushed deeper to the value reported in Table 1.</footnote>
dust opacity in high-z galaxies, would be NASA's future *Origins Space Telescope*.

While these ALMA data have provided the basis for a physical understanding of MAMBO-9, much has yet to be constrained. Further ALMA observations at higher spatial resolution could be used to investigate the internal dynamics on less than kiloparsec scales in both CO and FIR fine-structure lines. The VLA could constrain the full molecular gas reservoir in low-J CO and radio continuum to provide more direct constraints on gas mass and a possible AGN, respectively. It is important to emphasize the role of near-future facilities like *JWST* in unlocking the stellar emission in high-z sources like MAMBO-9. As Figure 6 shows, the spatial resolution of existing *Spitzer* data is insufficient to spatially resolve the pair at z = 5.85. Such systems are too faint to have been detected by *Herschel*, and they might sit just at the edge of detectability for powerful facilities like *HST*, which is only capable of probing rest-frame UV emission as it is heavily obscured. *JWST* can not only provide a much-needed stellar mass estimate, it can do it at much higher spatial resolution than has been possible in the past, while it will also probe rest-frame optical nebular line diagnostics, shedding light on its metal content and strength of the ionizing radiation field, crucial factors to measure in order to understand the unusually dust-rich environment.

MAMBO-9 remained hidden in plain sight for years as one of hundreds of 850 µm to 1 mm detected DSFGs without a secure redshift, and even after a sensitive ALMA spectral scan and suspicions of its high redshift, only tentative lines of marginal significance were used to identify its redshift as z = 5.85, first in Jin et al. (2019). Its identification is similar to other unlensed high-z DSFGs in the literature, in particular GN 20 at z = 4.05 (Daddi et al. 2009), AzTEC-1 at z = 4.34 (Yun et al. 2015), and HDF 850.1 at z = 5.18 (Walter et al. 2012), all of which took many years of effort before being spectroscopically confirmed by the Plateau de Bure Interferometer (PdBI) or the Redshift Search Receiver (RSR) at the Large Millimeter Telescope. Though ALMA is much more sensitive than PdBI (now the Northern Extended Millimeter Array, NOEMA) and RSR, it still requires a time investment of a few hours on-source to provide an unequivocal identification; such time investments per source have been rare in the first eight years of ALMA operations. This source, like those before it, highlights the severe need for systematic precision source follow-up of promising high-z DSFG candidates to secure accurate constraints on the early, obscured universe.

The authors thank Shuowen Jin for helpful discussions in the preparation of this manuscript, as well as the anonymous reviewer who provided valuable comments and suggestions. This paper makes use of the following ALMA data: ADS/JAO.ALMA #2018.1.00037.A, ADS/JAO.ALMA #2018.1.00231.S, ADS/JAO.ALMA #2017.1.00373.S, and ADS/JAO.ALMA #2016.1.00279.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. C.M.C. thanks the National Science Foundation for support through grants AST-1714528

and AST-1814034, and additionally C.M.C. and J.A.Z. thank the University of Texas at Austin College of Natural Sciences for support. In addition, C.M.C. acknowledges support from the Research Corporation for Science Advancement from a 2019 Cottrell Scholar Award sponsored by IF/THEN, an initiative of Lyda Hill Philanthropies. K.I.C. acknowledges funding from the European Research Council through the award of the Consolidator Grant ID 681627-BUILDUP. S.L.F. thanks the National Science Foundation for support through grant AST-1518183 and NASA through grant 80NSSC18K0954. K.K. acknowledges support from the Knut and Alice Wallenberg Foundation. M.T. acknowledges the support from grant PRIN MIUR 2017. S.T. acknowledges support from the ERC Consolidator Grant funding scheme (project ConTExT, grant No. 648179). The Cosmic DAWN Center is funded by the Danish National Research Foundation under grant No. 140. E.T. acknowledges support from CONICYT-Chile grants Basal-CATA AFB-170002, FONDECYT Regular 1160999 and 1190818, and Anillo de Ciencia y Tecnología ACT1720033.

## ORCID iDs

Caitlin M. Casey https://orcid.org/0000-0002-0930-6466

Jorge A. Zavala https://orcid.org/0000-0002-7051-1100

Manuel Aravena https://orcid.org/0000-0002-6290-3198

Karina I. Caputi https://orcid.org/0000-0001-8183-1460

Jaclyn B. Champagne https://orcid.org/0000-0002-6184-9097

David L. Clements https://orcid.org/0000-0002-9548-5033

Elisabete da Cunha https://orcid.org/0000-0001-9759-4797

Patrick Drew https://orcid.org/0000-0003-3627-7485

Steven L. Finkelstein https://orcid.org/0000-0001-8519-1130

Christopher C. Hayward https://orcid.org/0000-0003-4073-3236

Jeyhan S. Kartaltepe https://orcid.org/0000-0001-9187-3605

Kirsten Knudsen https://orcid.org/0000-0002-7821-8873

Anton M. Koekemoer https://orcid.org/0000-0002-6610-2048

Nick Z. Scoville https://orcid.org/0000-0002-0438-3323

Kartik Sheth https://orcid.org/0000-0002-5496-4118

Justin Spilker https://orcid.org/0000-0003-3256-5615

Johannes Staguhn https://orcid.org/0000-0002-8437-0433

Margherita Talia https://orcid.org/0000-0003-4352-2063

Yoshiaki Taniguchi https://orcid.org/0000-0003-2247-3741

Ezequiel Treister https://orcid.org/0000-0001-7568-6412

Min Yun https://orcid.org/0000-0001-7095-7543

## References

* Apostolovski, Y., Aravena, M., Anguita, T., et al. 2019, A&A, 628, A23

* Aretxaga, I., Wilson, G. W., Aguilar, E., et al. 2011, MNRAS, 415, 3831

* Asano, R. S., Takeuchi, T. T., Hirashita, H., & Inoue, A. K. 2013, EP&S, 65, 213

* Ashby, M. L. N., Willner, S. P., Fazio, G. G., et al. 2015, ApJS, 218, 33

* Baugh, C. M., Lacey, C. G., Frenk, C. S., et al. 2005, MNRAS, 356, 1191

* Behrens, C., Pallottini, A., Ferrara, A., Gallerani, S., & Vallini, L. 2018, MNRAS, 477, 552

* Behroozi, P., Wechsler, R., Hearin, A., & Conroy, C. 2018, MNRAS, 488, 3143

* Bergin, E. A., Kaufman, M. J., Melnick, G. J., Snell, R. L., & Howe, J. E. 2003, ApJ, 582, 830

* Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

* Bertoldi, F., Carilli, C., Aravena, M., et al. 2007, ApJS, 172, 132
Bethermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Bothwell, M. S., Aguirre, J. E., Chapman, S. C., et al. 2013, ApJ, 779, 67

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Cappellari, M., Bacon, R., Bureau, M., et al. 2006, MNRAS, 366, 1126

Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105

Caselli, P., Keto, E., Pagani, L., et al. 2010, A&A, 521, L29

Casey, C., Capak, P., Staguhn, J., et al. 2019, BAAS, 51, 212

Casey, C. M. 2012, MNRAS, 425, 3094

Casey, C. M., Berta, S., Béthermin, M., et al. 2012, ApJ, 761, 140

Casey, C. M., Chen, C.-C., Cowie, L. L., et al. 2013, MNRAS, 436, 1919

Casey, C. M., Cooray, A., Killi, M., et al. 2017, ApJ, 840, 101

Casey, C. M., Hodge, J., Zavala, J. A., et al. 2018a, ApJ, 862, 78

Casey, C. M., Narayanan, D., & Cooray, A. 2014, PhR, 541, 45

Casey, C. M., Zavala, J. A., Spilker, J., et al. 2018b, ApJ, 862, 77

Chabrier, G. 2003, PASP, 115, 763

Chary, R.-R., Stern, D., & Eisenhardt, P. 2005, ApJL, 635, L5

Chiang, Y.-K., Overzier, R. A., Gebhardt, K., & Henriques, B. 2017, ApJL, 844, L23

Condon, J. J. 1992, ARA&A, 30, 575

Conley, A., Cooray, A., Vieira, J. D., et al. 2011, ApJL, 732, L35

Conroy, C., Schiminovich, D., & Blanton, M. R. 2010, ApJ, 718, 184

Conroy, C., & Wechsler, R. H. 2009, ApJ, 696, 620

Cooray, A., Calanog, J., Wardlow, J. L., et al. 2014, ApJ, 790, 40

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Daddi, E., Dannerbauer, H., Stern, D., et al. 2009, ApJ, 694, 1517

Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4

de Rossi, M. E., Rieke, G. H., Shivaei, I., Bromm, V., & Lyu, J. 2018, ApJ, 869, 4

Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615

Dwek, E., Staguhn, J., Arendt, R. G., et al. 2014, ApJL, 788, L30

Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119

Fardal, M. E., Katz, N., Weinberg, D. H., Davé, R., & Hernquist, L. 2001, arXiv:astro-ph/0107290

Finkelstein, S. L., Ryan, R. E., Jr., Papovich, C., et al. 2015, ApJ, 810, 71

Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789

González, J. E., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2011, MNRAS, 413, 749

González-Alfonso, E., Fischer, J., Isaak, K., et al. 2010, A&A, 518, L43

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23

Harrison, I., & Hotchkiss, S. 2013, JCAP, 7, 022

Hodge, J. A., Swinbank, A. M., Simpson, J. M., et al. 2016, ApJ, 833, 103

Jarugula, S., Vieira, J. D., Spilker, J. S., et al. 2019, ApJ, 880, 92

Jin, S., Daddi, E., Liu, D., et al. 2018, ApJ, 864, 56

Jin, S., Daddi, E., Magdis, G. E., et al. 2019, arXiv:1906.00040

Jones, A. P., Fanciullo, L., Köhler, M., et al. 2013, A&A, 558, A62

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Koprowski, M. P., Dunlop, J. S., Michałowski, M. J., et al. 2017, MNRAS, 471, 4155

Kroupa, P., & Weidner, C. 2003, ApJ, 598, 1076

Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2008, MNRAS, 385, 1155

Lagos, C. d. P., Robotham, A. S. G., Trayford, J. W., et al. 2019, MNRAS, 489, 4196

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Le Floc'h, E., Aussel, H., Ilbert, O., et al. 2009, ApJ, 703, 222

Le Floc'h, E., Papovich, C., Dole, H., et al. 2005, ApJ, 632, 169

Li, A., & Draine, B. T. 2001, ApJ, 554, 778

Liang, L., Feldmann, R., Kereš, D., et al. 2019, MNRAS, 489, 1397

Liu, L., Weiß, A., Perez-Beaupuits, J. P., et al. 2017, ApJ, 846, 5

Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, A&A, 532, A90

Ma, X., Hayward, C. C., Casey, C. M., et al. 2019, MNRAS, 487, 1844

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Marrone, D. P., Spilker, J. S., Hayward, C. C., et al. 2018, Natur, 553, 51

Matsuura, M., Barlow, M. J., Zijlstra, A. A., et al. 2009, MNRAS, 396, 918

Matsuura, M., Wood, P. R., Sloan, G. C., et al. 2006, MNRAS, 371, 415

Michałowski, M. J., Hayward, C. C., Dunlop, J. S., et al. 2014, A&A, 571, A75

Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27

Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569

Narayanan, D. 2015, in IAU Symp. 309, Galaxies in 3D across the Universe, ed. B. L. Ziegler (Cambridge: Cambridge Univ. Press), 214

Nguyen, H. T., Schulz, B., Levenson, L., et al. 2010, A&A, 518, L5

Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 596, A100

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, arXiv:1807.06209

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31

Rowan-Robinson, M., Oliver, S., Wang, L., et al. 2016, MNRAS, 461, 1100

Salmon, B., Papovich, C., Finkelstein, S. L., et al. 2015, ApJ, 799, 183

Schinnerer, E., Smolčić, V., Carilli, C. L., et al. 2007, ApJS, 172, 46

Scoville, N. 2003, JKAS, 36, 167

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83

Scoville, N. Z., Polletta, M., Ewald, S., et al. 2001, AJ, 122, 3017

Shim, H., Chary, R.-R., Dickinson, M., et al. 2011, ApJ, 738, 69

Simpson, J. M., Smail, I., Swinbank, A. M., et al. 2015, ApJ, 807, 128

Smit, R., Bouwens, R. J., Labbé, I., et al. 2014, ApJ, 784, 58

Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1

Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677

Spilker, J. S., Marrone, D. P., Aguirre, J. E., et al. 2014, ApJ, 785, 149

Spilker, J. S., Marrone, D. P., Aravena, M., et al. 2016, ApJ, 826, 112

Strandet, M. L., Weiss, A., De Breuck, C., et al. 2017, ApJL, 842, L15

Strandet, M. L., Weiss, A., Vieira, J. D., et al. 2016, ApJ, 822, 80

Swinbank, A. M., Smail, I., Chapman, S. C., et al. 2004, ApJ, 617, 64

Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246

Thompson, T. A., Quataert, E., & Murray, N. 2005, ApJ, 630, 167

Toft, S., Zabl, J., Richard, J., et al. 2017, Natur, 546, 510

Walter, F., Decarli, R., Carilli, C., et al. 2012, Natur, 486, 233

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Williams, C. C., Labbe, I., Spilker, J., et al. 2019, ApJ, 884, 154

Yang, C., Gao, Y., Omont, A., et al. 2013, ApJL, 771, L24

Yang, C., Gavazzi, R., Beelen, A., et al. 2019, A&A, 624, A138

Yun, M. S., Aretxaga, I., Gurwell, M. A., et al. 2015, MNRAS, 454, 3485

Yun, M. S., Reddy, N. A., & Condon, J. J. 2001, ApJ, 554, 803

Zavala, J. A., Montaña, A., Hughes, D. H., et al. 2018, NatAs, 2, 56

Zhukovska, S., Gail, H.-P., & Trieloff, M. 2008, A&A, 479, 453