A Long Stream of Metal-Poor Cool Gas around a Massive Starburst Galaxy at z = 2.67

We present the first detailed dissection of the circumgalactic medium (CGM) of massive starburst galaxies at z>2. Our target is a submillimeter galaxy (SMG) at z = 2.674 that has a star formation rate of 1200 $M_\odot$/yr and a molecular gas reservoir of $1.3\times10^{11} M_\odot$. We characterize its CGM with two background QSOs at impact parameters of 93 kpc and 176 kpc. We detect strong HI and metal-line absorption near the redshift of the SMG towards both QSOs, each consisting of three main subsystems spanning over 1500 km/s. The absorbers show remarkable kinematic and metallicity coherence across a separation of 86 kpc. In particular, the cool gas in the CGM of the SMG exhibits high HI column densities ($\log N_{\rm HI}/{\rm cm}^{-2} = 20.2, 18.6$), low metallicities ([M/H] $\approx$ -2.0), and similar radial velocities ($\approx$ -300 km/s). While the HI column densities match previous results on the CGM around QSOs at z>2, the metallicities are lower by more than an order of magnitude, making it an outlier in the line width$-$metallicity relation of damped Ly$\alpha$ absorbers. The large physical extent, the velocity coherence, the high surface density, and the low metallicity are all consistent with the cool, inflowing, and near-pristine gas streams predicted to penetrate hot massive halos at z>1.5. We estimate a total gas accretion rate of ~100 $M_\odot$/yr from three such streams, which falls short of the star formation rate but is consistent with simulations. At this rate, it takes about a gigayear to acquire the molecular gas reservoir of the central starburst.


INTRODUCTION
The global gas supply for in situ star formation is a central question in galaxy formation and evolution, because star formation and merging are the two primary channels through which galaxies grow (Oser et al. 2010). According to spherical hydrodynamical models (Birnboim & Dekel 2003) and cosmological simulations (Keres et al. 2005), stable accretion shocks are established near the virial radius when a dark matter (DM) halo grows to a mass threshold of M shock = 2 − 3 × 10 11 M ⊙ . So in massive halos, a significant fraction of the accreted gas is expected to be shock-heated to the virial temperature (T vir = 8 × 10 6 (M halo /10 13 M ⊙ ) 2/3 K) and develops an atmosphere of hot diffuse gas. The virial shock effectively cuts off the fuel supply for star formation, because of the inefficient radiative cooling of the hot gas even in the denser inner regions (Kereš et al. 2009). But at high redshift, narrow filaments of cool gas (T 10 5 K) from the cosmic web may penetrate the hot atmospheres of rare, massive halos without ever being shock-heated to the virial temperature, thanks to the lower masses of typical halos at higher redshifts that define the width of the filaments (Dekel & Birnboim 2006;Dekel et al. 2009). In fact, this cold mode accretion may dominate over the hot mode accretion (radiative cooling of shock-heated virialized gas) at all halo masses at z > 2 ( Kereš et al. 2009).
In emission, the predicted cold-mode accretion streams (or "cold streams" in short) feeding high-redshift massive galaxies may appear as giant filamentary Lyα nebulae around QSOs (Weidinger et al. 2004;Cantalupo et al. 2014;Martin et al. 2015Martin et al. , 2019 and in dense protocluster environments (Møller & Fynbo 2001;Hennawi et al. 2015;Umehata et al. 2019;Li et al. 2019;Daddi et al. 2020). However, it has been difficult to rule out outflows as the alternative interpretation, especially when QSO photoionization contributes to the Lyα emission and the chemical abundance of the nebulae cannot be easily measured. In absorption, cold streams can be detected and distinguished from other gaseous components based on neutral hydrogen (H I) column density, kinematics, and particularly chemical abundance (Fumagalli et al. 2011;Theuns 2021).
In this project, we have selected a sample of massive starbursts at high redshifts in the vicinity of background QSOs to trace the cool gas supply in these early massive halos. There, the problem of gas supply is the most acute because of the extremely short gas exhaustion timescale. We then utilize the absorption-line spectra of background QSOs to characterize the physical state of their circumgalactic medium (CGM)the gas between the inner regions of galaxies and the diffuse intergalactic medium (IGM) -and to search for large-scale cool gas reservoirs. In this section, we review our knowledge of the target galaxy sample and QSO absorption-line systems in the literature. These earlier studies have motivated this project and will provide valuable reference samples that can be compared with the system dissected in this work.

Submillimeter Galaxies
Heated dust in the interstellar medium cools by emitting a modified blackbody spectrum (MBB; S ν ∝ (1 − e −τν )B ν (T )) with temperatures in the range 10 K T 100 K, form-ing the far-infrared (IR) hump in the spectral energy distribution (SED) of a galaxy. At any given frequency along the Rayleigh-Jeans tail where the dust should be optically thin, the observed flux density (S ν,obs ) of the MBB is proportional to the dust mass (M dust ), the dust temperature (T ), and the redshift (z) - where d A (z) is the angular diameter distance at redshift z (which varies by only 22% between z = 1 and z = 4) and β ≈ 2 is the dust emissivity parameter (κ ν ∝ ν β ). Therefore, galaxies selected at long wavelengths, such as the (sub)millimeter regime, preferentially have higher dust mass, higher dust temperature, and are at higher redshift than galaxies selected at shorter wavelengths. Furthermore, holding the metallicity (Z gas ) constant, high dust mass together with high temperature implies that the galaxies are gas-rich (M gas = M dust /Z gas ) and have high star-formation efficiency (SFE = SFR/M gas , where SFR is the star formation rate), because a result from the Stefan-Boltzmann law. Indeed, follow-up observations of the brightest galaxies selected at 850 µm (S 850 3 mJy), the submillimeter galaxies (SMGs; Smail et al. 1997;Barger et al. 1998;Blain et al. 2002), have revealed a significant population of gas-rich starburst galaxies that contribute almost as much to the cosmic SFR density as UV-selected Lyman break galaxies at z = 2 − 3 (Chapman et al. 2005;Casey et al. 2014). The SMGs are mature ( M star ∼ 10 11 M ⊙ ; Hainline et al. 2011;Michałowski et al. 2012;Targett et al. 2013), metal-rich ( Z ∼ Z ⊙ ; Swinbank et al. 2004), gas-rich ( M mol ∼ 3 × 10 10 M ⊙ ; Greve et al. 2005;Tacconi et al. 2008;Ivison et al. 2011;Bothwell et al. 2013), extreme star-forming systems (SFR ∼ 500 M ⊙ yr −1 ) with a broad redshift distribution that peaks at z ∼ 2.5 (Chapman et al. 2005;Wardlow et al. 2011). The molecular gas reservoirs are turbulent, likely due to starburst-driven galactic outflows (Falgarone et al. 2017). Notably, the nearly linear relation between CO and IR luminosities implies an almost constant gas depletion timescale of τ dep ≡ M mol /SFR ∼ 0.1 Gyr (Bothwell et al. 2013), which is far shorter than that of normal star-forming galaxies on the main sequence (τ dep ∼ 0.6 Gyr at z = 2.5; Tacconi et al. 2018), justifying the usage of "starburst" in describing SMGs 9 .
On the other hand, the autocorrelation length for SMGs of ∼11 Mpc at z = 1 − 3 implies a characteristic dark matter halo mass of M halo ∼ 9 × 10 12 M ⊙ for h = 0.7 (Hickox et al. 2012). The high halo mass is consistent with the high maximum rotation velocities (V circ 500 km s −1 ) observed in several bright SMGs with spatially resolved kinematics (e.g., Hodge et al. 2012;Xue et al. 2018). For Navarro-Frenk-White (NFW) halos (Navarro et al. 1996) at z = 2.5, the halo mass is di-9 Although most of the difference in τ dep is driven by the conversion factor from CO to molecular gas (α CO ≡ M mol /L ′ CO in units of M ⊙ /(K km s −1 pc 2 )), constraints from dynamical masses and dust masses have shown that SMGs indeed have lower α CO than that of normal starforming galaxies (e.g., Hodge et al. 2012;Magnelli et al. 2012a;Xue et al. 2018) and that the value adopted from local ultraluminous IR galaxies (ULIRGs, α CO = 1.0; Downes & Solomon 1998;Papadopoulos et al. 2012) is more appropriate for SMGs than the Galactic value (α CO = 4.3; Bolatto et al. 2013). rectly related to the maximum circular velocity by a power law: M halo = 10 13 M ⊙ (V circ /500 km s −1 ) 3 (Bullock et al. 2001;Klypin et al. 2011). Such a mass is well above the threshold mass for stable virial shocks (M shock ), and atmospheres of hot gas at the virial temperature (∼ 8 × 10 6 K) are expected to fill the halo. But as previously discussed, at the early epoch of the SMGs, cool gas filaments can penetrate their halos, which could potentially deliver enough gas to build the molecular gas reservoir that supports the ongoing intense star formation.

QSO Absorption-line Systems
Ever since the discovery of multiple absorption redshifts in QSO spectra (Burbidge et al. 1968), quasar absorption-line spectroscopy has become a powerful tool to study diffuse gas at various phases in the IGM and the CGM, which account for the majority of the baryonic mass in the universe (see Péroux & Howk 2020 for a recent review). The optical depth at the Lyman limit (λ rest = 912 Å) reaches unity when the H I column density reaches log N HI = 17.2 10 Accumulating evidence suggests that these optically thick absorbers trace material in virialized structures (i.e., the CGM, Fumagalli et al. 2016;Lehner et al. 2016), while the optically thin absorbers in the Lyα forest (LYAF) likely trace the IGM (Rauch 1998). Due to their distinct physical properties, the optically thick absorbers are empirically subdivided into three categories based on their H I column densities: the Lyman limit systems (LLSs, 17.2 ≤ log N HI < 19) that are mostly ionized, the damped Lyα absorbers (DLAs; log N HI ≥ 20.3) that are mostly neutral, and lastly the super-LLSs or sub-DLAs 11 for the intermediate category of absorbers with 19 ≤ log N HI < 20.3. Unlike absorbers at lower column densities, gas in the DLAs is mostly neutral. In fact, at all epochs since z ∼ 5, the DLAs have contained most of the neutral gas that is poised to fuel star formation in galaxies (Wolfe et al. 2005).

Emission−Absorption Connection
Because the H I column density threshold of DLAs was set by the observed limit of 21 cm emission at the cutoff boundaries of nearby spiral disks (Wolfe et al. 1986), the DLAs were expected to arise from gas-rich galactic disks even at high redshifts. However, the emission counterparts (i.e., the DLA galaxies) of most DLAs have eluded detection. Among the limited detections in optical searches, it is found that the DLA galaxies are very faint (r 24) and close (∼2 ′′ ) to the QSOs (e.g., Steidel & Hamilton 1992;Fynbo et al. 2008), making it difficult to measure their redshifts. To improve efficiency, searches of the emission counterparts of DLAs have focused on DLAs that are clearly chemically enriched ([M/H] > −0.7) (e.g., Fynbo et al. 2013;Jorgenson & Wolfe 2014) or sightlines that pass through multiple DLAs (e.g., Srianand et al. 2016). But still, only 16 z > 1.9 DLA host galaxies have been identified via emission lines in the optical (see Krogager et al. 2017;Møller & Christensen 2020, for compilations) over an extensive search period of nearly three decades. The advent of (sub)millimeter interferometers such as the Atacama Large Millimeter/submillimeter Array (ALMA) and the Very Large Array (VLA) have significantly improved the success rate of identifying absorption-selected galaxies, because (1) the contrast between DLA host galaxies and the QSO is more favorable at longer wavelengths, and (2) the interferometers have an unattenuated view over a wide FoV (thus does not require lucky slit placements). In only a few years, there have been four z ∼ 4 DLA galaxies identified in [C II] 158 µm (Neeleman et al. 2017 and five z ∼ 2 DLA in CO(4-3) and CO(3-2) (Kanekar et al. 2020). Interestingly, the DLA galaxies previously identified in the optical/near-IR are not detected in CO and vice versa (Kanekar et al. 2020), indicating that observations at different wavelength are complementary to one another and that H I-absorption-selection tags gas-rich galaxies of all types. Some DLAs also have multiple emission counterparts that are consistent with the absorption redshifts, suggesting a group/cluster environment (e.g., Fynbo et al. 2018).
The opposite approach from the searches of DLA galaxies is to start from an emission-selected galaxy sample and look for corresponding absorption lines in the spectra of nearby background QSOs. This approach requires chance alignments of foreground galaxies and background QSOs, thus requires large samples of both populations. The implicit assumption is that the emission-selected galaxies have similar CGM properties, and therefore, the absorption signals obtained from different galaxy-QSO pairs can be combined to provide meaningful average properties of a typical halo in the studied galaxy population. The searches for absorbers are no longer limited to DLAs, but to all optically thick absorbers (i.e., LLSs and sub-DLAs). At z 2, the targeted galaxy populations have included Lyman Break Galaxies (LBGs) (e.g., Simcoe et al. 2006;Rudie et al. 2012;Rudie et al. 2013;Crighton et al. 2013Crighton et al. , 2015 and QSOs (e.g., Hennawi et al. 2006;Prochaska et al. 2013a;Lau et al. 2016). In addition, using a sample of projected QSO pairs where one of the QSOs intercepts a DLA, Rubin et al. (2015) have probed the CGM of the DLA galaxy without identifying the DLA galaxy in emission. These studies have mapped out the H I column density, the ion ratios, and the metallicity as a function of impact parameters (R ⊥ ). The covering fraction of optically thick H I absorbers increases from ∼30% around LBGs (Rudie et al. 2012) and DLAs (Rubin et al. 2015) to 60% around QSOs (Prochaska et al. 2013a) for sightlines out to R ⊥ = 100 − 200 kpc (comparable to the virial radius of DM halos with M halo = 10 12.5 M ⊙ at z = 2: R vir = 154 kpc). The abundance of neutral gas in the halos of QSOs is particularly puzzling. Simulations predict that such massive halos are dominated by a hot T ∼ 10 7 K virialized plasma and a significantly lower covering factor of optically thick H I absorbers (e.g., Faucher-Giguère et al. 2015).
The unexpectedly large covering factor of LLSs around z ∼ 2 QSOs and the difficulty of reproducing the SMG population in galaxy formation models, could both be symptoms of the same problem. Attempting to reduce this tension between observations and theory, more recent cosmological zoom-in simulations have implemented recipes of stronger and presumably more realistic stellar feedback, which manages to preserve cool gas reservoirs in the accreted sub-halos during earlier phases of star formation before their infall into the massive halo. The presence of these gas-rich sub-halos increases the cool gas covering factor around QSOs (Faucher-Giguère et al. 2016) and their prolonged bombardment to the central galaxy leads to a rising star formation history that eventually produces SMGs between 2 < z < 3 (Narayanan et al. 2015;Lovell et al. 2021).

Organization
QSO absorption-line spectroscopy combined with efficient emission-line mapping provides a powerful method to link star-forming galaxies with the neutral gas reservoir that may fuel future star formation. Our understanding of the formation and evolution of massive galaxies is severely limited by the lack of observational constraints of the CGM of SMGs. The advent of Herschel large-scale far-infrared surveys have provided an opportunity to use projected SMG−QSO pairs to probe the CGM of SMGs. In this paper, we focus on one particularly interesting system -GAMA J0913−0107 -where two background QSOs have revealed an unusually H I-rich CGM around a luminous SMG.
The main text of the paper is organized as follows. We first provide an orientation of the system in § 2, then proceed with a detailed study of the emission sources ( § 3) and the absorption-line systems ( § 4), before finally drawing connections between the absorbers and their emission counterparts in § 5. We conclude the paper with a summary of the main results and a discussion of the implications in § 6. To keep the main text focused on the SMG−DLA system at z ≈ 2.67, we move additional material to the Appendices. We analyze the nearby optical source to the SMG and its potential lensing effect in Appendix A, present the methodology and result of our blind search of line emitters in the ALMA band-3 data in Appendix B, give an inventory of the line-of-sight contaminating absorbers at other redshifts in Appendix C, provide tables of detailed ionic column density and metallicity measurements in Appendix D, and describe our attempt to detect CO emission from faint optical sources near QSO1 and the identification of Comp b in Appendix E.
Throughout this paper, we adopt a model optimization method that combines a heuristic χ 2 minimization algorithm with a Markov Chain Monte Carlo (MCMC) algorithm (hereafter, the "AMOEBA + MCMC" method). It begins with using the downhill simplex method AMOEBA (Press et al. 1992) with simulated annealing AMOEBA_SA to find the solution that minimizes the residual. Although computationally more expensive than other least-χ 2 solvers (e.g., the Levenberg-Marquardt technique), AMOEBA_SA has the advantage of avoiding being trapped in local minima in a multidimensional parameter space. This advantage is particularly important in more complex problems such as fitting the H I absorption profiles with many Voigt profiles ( § 4.2). Next, starting from the minimum-χ 2 solution of AMOEBA_SA, we use the Differential Evolution MCMC algorithm (Ter Braak 2006) implemented in EXOFAST_DEMC (Eastman et al. 2013(Eastman et al. , 2019 to obtain the final solution and the statistical uncertainties of the parameters. The EXOFAST_DEMC routine first determines the stepping scale of each parameter by varying it from the minimum χ 2 solution until the χ 2 increases by one. It then starts the chains from positions that are randomly offset from the minimum χ 2 solution. The routine stops when the chains are considered well-mixed and the steps in the initial "burn-in" phase are removed. The marginalized 1σ confidence interval of each parameter is determined from the values at 15.8 and 84.1 percentiles of the concatenated chains, and the median values are adopted as the formal solution. Parameters derived from the model parameters are treated likewise: their formal values and uncertainties are calculated from the 50, 15.8, and 84.1 percentiles of the array directly calculated from the chains of model parameters. We assume the ΛCDM cosmology with Ω m = 0.3, Ω Λ = 0.7, and h ≡ H 0 /(100 km s −1 Mpc −1 ) = 0.7 and quote proper/physical distances. Fig. 1 illustrates the GAMA J0913−0107 system, where we label its major components: the SMG, its CO companions, and the two background QSOs. Table 1 lists their coordinates, redshifts, CO (3−2) luminosities, along with the impact parameters of the QSO sightlines.
Our ALMA 345 GHz imaging pinpointed the position of the Herschel source (Fu et al. 2017), and Gemini near-IR and ALMA 94 GHz spectroscopy jointly determined a spectroscopic redshift of 2.674 ( § 3.1 & 3.2). We designate the SMG as ALMA J091339.55−010656.4 or SMM J0913 in short. In addition, our ALMA 94 GHz observations detected companion galaxies in CO (3−2) near the redshift of the SMG: Comp a at z = 2.6747, and Comp b at z = 2.6884, 2.6917. Both are within 23 ′′ of the SMG position. Comp b has two redshifts because it is a superposition of two galaxies (see § 3.5). Because Herschel has FWHM resolutions of 18 ′′ , 25 ′′ , and 35 ′′ at 250, 350, and 500 µm, respectively, the SMG and its companions are blended in the Herschel images. But the contribution of the companions to the Herschel fluxes should be negligible given their orders-of-magnitude lower CO line luminosities.
When comparing the CO map with the KiDS r-band image in Fig. 1, we notice an r = 21.6 optical source just 0.8 ′′ from the SMG. As we will show in Appendix A, it is a foreground galaxy with a spec-z of z = 0.055 and its gravitational lensing effect on the SMG is negligible. We also find that Comp b is just ∼3 ′′ from an elongated optical source to the NNE. In the KiDS DR4 catalog, the optical source has a designation of J091338.527−010703.60 and magnitudes of r = 23.8 ± 0.1 and H = 22.0 ± 0.3. It has a photo-z of z p = 0.79 +0.45 −0.06 . Its SED shows a clear 1-magnitude drop-off between Y -band and Z-band, corresponding to a 4000 Å-break at z ∼ 1.2 − 1.5, which is consistent with the maximum-likelihood photo-z of z p = 1.37 in the catalog. We thus conclude that the optical source is most likely a foreground galaxy, although the extraction of an ALMA spectrum near the position of the optical source led to the identification of Comp b (Appendix E).

THE SUBMILLIMETER GALAXY AND ITS COMPANIONS 3.1. ALMA Position and Near-IR Spectroscopy
Herschel/SPIRE positions have large uncertainties. A comparison between ALMA and Herschel positions showed a 1σ positional offset of ∼4.2 ′′ for sources with S/N ∼ 6 at 250 µm (Eq. 5 of Fu et al. 2017). Near-IR slit spectroscopy requires sub-arcsec positions, so we carried out 0.5 ′′ -resolution ALMA band-7 (345 GHz/870 µm) observations of GAMA J0913−0107 as part of our Cycle-3 project 2015.1.00131.S (see Fu et al. 2017, for details). GAMA J0913−0107 shared an hour-long observing session with nine other Herschel SMGs in the same H-ATLAS field. With four scans and 44 antennas, we accumulated a total on-source integration time of 189.5 s. A single high S/N source is detected in the ∼17 ′′ Full Width at Half Power (FWHP) of the primary beam. It has a 870 µm flux density of S 870 = 7.4 ± 0.5 mJy, an offset from the Herschel position by 4.1 ′′ , and a beam-deconvolved FWHM of 0.46 ′′ ± 0.04 ′′ along the major axis (which corresponds to ∼4 kpc at z = 2.5).
The accurate ALMA position enabled our follow-up near-IR slit spectroscopy. Observations of SMM J0913 were carried out with Gemini near-infrared spectrograph (GNIRS; Elias et al. 2006) on 2017 Feb 19 as part of our queue program GN-2017A-Q-31. The A0-type star HIP49125 (V = 7.19, K = 6.553 Vega mag) was observed right after the science target to provide telluric correction and flux calibration. Because our goal was to measure the redshift, we used the cross-dispersed mode with the 32 l/mm grating to achieve a continuous wavelength coverage between 0.85 and 2.5 µm (orders 3 to 8). After applying a 30 ′′ offset from an offset star to the NE, we placed the 1 ′′ -wide 7 ′′ -long slit on the ALMA 870 µm position at a position angle (PA) of 73 deg (E of N; the PA was chosen to reach a guide star). The expected spectral resolution of this configuration is R = 510 (FWHM = 590 km s −1 ), but the actual spectral resolution may be higher depending on the source size and the seeing. We took 24 exposures of 136 s with a 3 ′′ -step ABBA dithering pattern. Data reduction was performed with a modified version of Spextool (Cushing et al. 2004) for GNIRS by K. Allers. The final coadded spectrum in Fig. 2 includes all 24 frames and has a total on-source time of 54.4 min. We detected an emission line at ∼4σ-level at 2.4121 µm (heliocentric corrected, vacuum wavelength), which we identified as Hα (λ rest = 6564.63 ) at z Hα = 2.6743 ± 0.0003. The [N II] λ6585.28 line is undetected, likely due to the elevated background noise at its wavelength and its lower flux. Our best-fit Gaussian model yields a line FWHM = 230 ± 60 km s −1 (i.e., the line is unresolved) and a line flux of F Hα = (6.1 ± 1.0) × 10 −17 erg s −1 cm −2 , which are comparable to those of the SMGs identified by the VLA (Alaghband-Zadeh et al. 2012;Fu et al. 2016).

ALMA CO (3−2) Spectral Line Imaging
To detect the molecular gas reservoir that fuels the intense star formation in the SMG, we carried out deep ALMA band-3 (100 GHz/3 mm) spectral line observations of SMM J0913 on 2018 December 10 and 15 with our Cycle 6 project 2018.1.00548.S. We tuned the four 1.875 GHz-bandwidth spectral windows to center on 92.2 (BB3), 94.0 (BB4), 104.2 Middle -An ALMA map zoomed in onto the SMG showing CO (3−2) emission between 2.67 < z < 2.70. This 36.5 ′′ ×50.4 ′′ region encloses the SMG and its CO companion galaxies (red tickmarks), and the two QSOs in the background of the system (black tickmarks). This composite CO image is formed by combining the 11 channels within ±140 km s −1 of z = 2.674 (i.e., ν obs = 94.120 ± 0.039 GHz). To show the CO emission from Comp b, the 8 ′′ ×8 ′′ region centered on Comp b (dotted box) is formed by combining the two channels where CO emission is detected (ν obs = 93.7535, and 93.6676 GHz, corresponding to z = 2.6884 and 2.6917). The contours are drawn at −3 (black dotted), 3, 4, 5(black solid), 20, 40, and 60σ (white solid). The synthesized beam of 1.6 ′′ ×1.3 ′′ is shown at the lower right corner. Right -A deep r-band image of the same region from KiDS (5σ detection limit at ∼25 mag). In all images, the position of QSO1 sets the origin of the coordinates. · · · · · · · · · · · · · · · SDSS J091338.30−010708.6 QSO2 138.4096520 −1.1190520 2.7488 5.6 × 10 9 10.8 · · · 85.0 · · · NOTE. -QSO coordinates are from the KiDS DR4 catalog, the SMG coordinates are from the higher-resolution ALMA band-6 image, and the coordinates of Comp a and Comp b are from their CO emission detected in the ALMA band-3 datacube. Comp b has two rows because it is a superposition of two CO emitters likely involved in a merger. The quoted redshift of QSO2 is from the CO (3−2) line, which is slightly off from its optical redshift (zopt = 2.7498, δv = 80 km s −1 ). The column L ′ CO3−2 lists the CO (3−2) line luminosities from ALMA. The CO emission of QSO1 is outside of our spectral coverage. The columns θ 1 and θ 2 list the angular separations from QSO1 and QSO2, respectively, and the corresponding transverse proper distances at the source redshift (i.e., the impact parameters) are listed as R ⊥,1 and R ⊥,2 .
(BB1), and 106.0 GHz (BB2) in dual linear polarization mode (XX and YY). We chose a spectral averaging factor of 16 to bin the Frequency Division Mode (FDM)'s input channel spacing of 0.488 MHz to an output channel spacing of 7.8125 MHz. The spectral averaging significantly reduces the output data rate and essentially eliminates the correlation between adjacent channels introduced by the Hanning window function applied to the correlation functions (see § 5.5 in the ALMA Technical Handbook). The resulting spectral response function is basically a top-hat function with a width of 7.8125 MHz (i.e., the output channel spacing), which corresponds to a spectral resolution of 22−25 km s −1 . The two lower frequency spectral windows provide a continuous frequency range between 91.27 and 94.94 GHz, which covers the CO (3−2) line (ν rest = 345.79599 GHz and λ rest = 866.96337 µm) between 2.642 ≤ z ≤ 2.789, encompassing the SMG at z = 2.674 and QSO2 at z = 2.7498. The frequency range covers a velocity window between −2620 and +9240 km s −1 relative to the SMG redshift. The two higher frequency spectral windows provide a continuous frequency coverage between 103.27 and 106.94 GHz, which covers the CO (3−2) line between 2.234 ≤ z ≤ 2.348 and traces the continuum emission at rest-frame frequencies around ν rest = 386 GHz (λ rest = 776 µm) at z = 2.674.
The primary beam of the ALMA 12-m antennas has an FWHP of 62 ′′ at 94 GHz. We set the field center at R.A. = 09 h 13 m 38.89 s , Decl. = −01 • 07 ′ 03.6 ′′ , which is near the position of QSO1 but ∼12 ′′ offset from the SMG. Three of the four planned observing sessions were executed, accumulating a total on-source time of 143.8 min with 6.05 s integrations. Either 43 or 46 12-m antennas were operational, with baselines ranging between 15.1 m and 740.5 m. The BL Lac object J0854+2006 served as the amplitude, bandpass, and pointing  calibrator, and the flat-spectrum radio quasar J0909+0121 as the phase calibrator (Bonato et al. 2019).

ALMA Data Processing
The raw visibility data were flagged and calibrated by the ALMA pipeline (Pipeline ver. 42030M, CASA ver. 5.4.0-68). The calibrated visibilities of the three observing sessions were then combined to form the final calibrated measurement set (MS). The pipeline worked very well. After inspecting the amplitudes of the calibrated visibilities, we found that additional flagging was only necessary for a tiny fraction of data. We used the CASA task flagcmd to flag the cross-correlation data of the antenna pair DA62 and DA65 in the 94.0 GHz spectral window between channels 188 and 191. We use the CASA task tclean to image the calibrated visibilities of each spectral window into spectral data cubes. When visibilities are gridded into regularized uv-cells, we adopt natural weighting to maximize the sensitivity. The synthesized beams are on average 1.7 ′′ ×1.3 ′′ in FWHM, so we set the imaging pixel size to 0.2 ′′ . In the spectral dimension, we retain the original channel spacing of 7.8125 MHz. The data were recorded in the Topocentric (TOPO) reference frame. Due to the motion of the Earth, every observing scan has a slightly different sampling in sky frequency. We image the data to the solar System Barycenter (BARY) reference frame to be consistent with the velocities measured in the heliocentric-corrected optical and near-IR spectra.
Significant continuum emission and a strong emission line at ∼94.1 GHz (in BB4) is detected at the ALMA 870 µm position of SMM J0913. To minimize the sidelobes from this bright source, we used Clark CLEAN deconvolution algorithm with a mask consisting of a single 2 ′′ -radius circle centered on the SMG. The CLEAN depth is set to be 2× the rms listed below.
For each spectral window, we generate two datacubes: one avoids interpolation in the spectral dimension by setting interpolation = nearest and is uncorrected for the primary beam, and the other uses linear interpolation and is corrected for the primary beam. The former is better suited for blind line searches because maps in adjacent channels remain uncorrelated, while the latter is better suited for measuring line parameters such as central frequency, width, and integrated flux. The resulting spectral cubes have a dimension of 540 pixels by 540 pixels by 240 channels. At the phase center, the sensitivities of the datacubes reach rms = 0.165, 0.171, 0.158, 0.155 mJy beam −1 channel −1 for BB1, BB2, BB3, and BB4, respectively. The rms values are consistent with the visibility noise that we measured with visstat (σ ∼ 250 mJy visibility −1 channel −1 ), because rms ∼ σ/ √ n ch n pol n baseline n int , where n ch = 1 (1 channel binning), n pol = 2 (2 polarizations), n baseline = n ant (n ant − 1)/2 = 903 for n ant = 43 (903 baselines for 43 antennae), and n int = t on_source /t int ∼ 1426 (total number of integrations).
In addition, we generate one continuum image from the two higher frequency spectral windows (BB1 and BB2). We used the Multi-term Multi-Frequency Synthesis (mtmfs) CLEAN algorithm with a linear spectral model (i.e., nterms = 2). Again, we used natural weighting, 0.2 ′′ pixel size, and the same CLEAN mask centered on the SMG. The sensitivity of the continuum image reaches rms = 7.4 µJy beam −1 , consistent with the rms of the spectral cubes divided by √ n ch ≃ 22. Fig. 3 shows the ALMA band-3 spectrum of the SMG. The spectrum is extracted with an elliptical aperture matching the beam-convolved source size. A prominent emission line peaks at ν obs = 94.12 GHz, which we identify as the CO (3−2) line at z CO = 2.67399 ± 0.00003. The CO detection thus confirms the Hα redshift from the GNIRS spectrum (z Hα = 2.6743 ± 0.0003, see § 3.1). Because the CO line is detected at a higher S/N and is less affected by dust extinction, we adopt the CO redshift for the SMG throughout the paper, i.e., z SMG = z CO = 2.674 (a slightly rounded-up value for simplicity).
0.54 ± 0.17 arcsec NOTE. -We have adopted a CO (3−2)/CO (1−0) brightness temperature ratio of r 31 = 0.52 and a CO to molecular gas conversion factor of Here and in Table 3, the intrinsic source sizes are given by the beam-deconvolved major-and minoraxis FWHMs (a & b), measured by CASA task imfit.
broad (FWHM ∼ 1000 km s −1 ) emission-line component with a peak flux of ∼0.4 mJy underneath the prominent narrow component (FWHM ∼ 250 km s −1 ). We thus model the CO spectrum with two Gaussians and compare its result with a single-Gaussian model.
We find that the improvement of the double-Gaussian model over the single-Gaussian model is highly significant. The formal double-Gaussian solution achieves a χ 2 = 220.2 for a degree-of-freedom (DOF) of 209. For comparison, the formal single-Gaussian solution achieves a χ 2 = 251.8 for The errors of the Herschel photometry include confusion noises of 5.3, 6.4 and 6.7 mJy beam −1 at 250, 350, and 500 µm, respectively (Pascale et al. 2011;Rigby et al. 2011). The IR luminosity, L IR , is defined as the integrated luminosity between 8 and 1000 µm at rest-frame. The dust-obscured SFR is estimated using the commonly used conversion: SFR/M ⊙ yr −1 = 10 −10 L IR /L ⊙ (e.g., Daddi et al. 2010). DOF = 212. According to the F-test, such a difference rejects the null hypothesis, that the double-Gaussian model does not provide a significantly better fit, at a confidence level of 99.99964% or 4.6σ. The result of the double-Gaussian fit from the "AMOEBA + MCMC" method is listed in Table 2. The broad component has an FWHM of ∼900 km s −1 and ac-counts for almost a quarter of the total emission-line flux. The existence of a narrow CO line on top of a broad CO line with essentially no velocity offset indicates that the SMG's intense star-forming nucleus (the narrow component) is either embedded in a fast-rotating disk or driving a bipolar outflow (the broad component). Unfortunately, the spatial resolution of the ALMA data is inadequate to distinguish between the two scenarios. We note that such a broad CO component would not have been detected in shallower spectroscopic data that are more generally available to SMGs, so this feature may not be unique to SMM J0913.
With the redshift determined, we fit the SED between 250 µm and 3 mm from Herschel/SPIRE and ALMA with a modified blackbody curve. We adopt the general solution of the radiative transfer equation assuming local thermal equilibrium at a constant temperature T : where B ν (T ) is the Planck function at a temperature of T and a rest-frame frequency ν, πr 2 s the effective size of the dust emitting region, and d L the luminosity distance. Assuming that the dust opacity follows a power-law with a negative slope of −β at wavelengths greater than the dust size (∼10 µm), the optical depth should follow the same power-law: where ν 0 (λ 0 ) is the rest-frame frequency (wavelength) at which the dust becomes optically thick. Given the dust massabsorption coefficient of κ = 0.07 m 2 kg −1 at 850 µm for Galactic dust (Dunne et al. 2000;James et al. 2002), it can be shown that the dust mass is: M dust = 9.0 × 10 9 M ⊙ (πr 2 s /kpc 2 ) (λ 0 /850 µm) β . (5) The result of the SED fit gives us a measure of the dust temperature, the dust-obscured SFR, the dust mass, and the effective size of the dust photosphere. Fig. 4 left shows the multi-band photometry, the median MCMC model, and the 1σ spread of the models. Table 3 lists the formal parameters and their uncertainties. The dust is relatively warm (T = 44 ± 7 K) and the dust photosphere has an effective size of 9 +10 −5 kpc 2 , comparable to the intrinsic source size measured at 343.5 GHzπab/4 = 6.4 ± 1.5 kpc 2 .
Combining the results from the SED fit and the CO (3−2) spectroscopy, we found a gas depletion timescale of ∼0.1 Gyr, which is similar to other SMGs but 6× shorter than co-eval main-sequence galaxies (Fig. 4 middle). The SMG's dust emission is resolved by the ALMA band-6 data with a beamdeconvolved size of 0.46 ′′ ×0.28 ′′ . Its CO (3−2) emission is resolved by the ALMA band-3 data with a beam-deconvolved size of 0.76 ′′ ×0.54 ′′ . In both cases, we have measured the intrinsic source sizes from CLEAN'ed images using the CASA task imfit. These size measurements allow us to place the SMG on the Kennicutt-Schmit relation (Fig. 4 right). It is characterized by high surface densities of SFR and molecular gas even compared to the SMG population. But like other SMGs, it features a high star-formation efficiency that is distinct from normal star-forming galaxies at z ∼ 2 (e.g., Daddi et al. 2010;Genzel et al. 2010).

Companion Galaxies
We carried out a blind search of line emitters in the ALMA datacubes with a matched-filter algorithm and tested the fidelity of the detections with simulated noise-only interferometer data (see Appendix § B). In the spectral window centered at 94 GHz (i.e., BB4), we found two robust line emitters: the SMG at z = 2.674 (S/N = 67) and Comp a at z = 2.6747 (S/N = 6.4). The other companion, Comp b, was detected at high significance when we combine the two ALMA channels closest in velocity to the metal-line-detected absorbing clouds C1 and C2 toward QSO1 (see § 4.3). The source has a peak S/N of 4.5 when combining the two channels, while its peak S/N is only 3.6 in the two individual channels, making it confused with noise spikes. The matched-filter algorithm fails to identify Comp b because it assumes that only adjacent channels can boost the S/N above the detection limit. In other words, the detection of Comp b is possible only because (1) we have utilized the prior knowledge of the redshifts of the absorption lines (with the implicit assumption that the emission counterparts have similar redshifts), and (2) the emission counterparts of the two clouds are superimposed on the sky (increasing the S/N when they are combined).
Our blind search detected four additional high-fidelity line emitters: QSO2 at z = 2.7488 (S/N = 7.0), its companion at z = 2.7392 (δv = −770 km s −1 ; S/N = 6.1) located ∼30 ′′ to the NE of QSO2, and two additional sources at z = 2.3452 (S/N = 5.5) and z = 2.3324 (S/N = 5.2) that may correspond to the z abs = 2.345 H I and C IV absorbers that appear toward both QSOs (see Fig. C1 in Appendix C). But on the other hand, only the SMG is detected in the continuum image of the two higher-frequency spectral windows (i.e., BB1 and BB2).
While the detection of CO in the SMG and QSO2 is expected, the detection of their companion CO emitters is not. We can use the ALMA Spectroscopic Survey in the Hubble Ultra Deep Field (ASPECS; Decarli et al. 2019) to estimate a baseline level of source-detection probability in normal field environments. The ASPECS covers an area of 4.6 arcmin 2 (PB ≥ 0.5, where PB is a correction factor for the primary beam pattern) with 17 pointings in band 3 and a spectral range of 21 GHz (84−105 GHz) with five tunings. The rms sensitivity varies with frequency with a range between 0.12 and 0.4 mJy beam −1 channel −1 for a channel spacing of 7.8 MHz (the same as ours). To provide a conservative estimate, we only count the 7 sources detected at S/N > 6 between 96 and 103 GHz (González-López et al. 2019), where rms ≃ 0.135 mJy beam −1 channel −1 . Only within this spectral range is ASPECS more sensitive than our data (rms ≃ 0.16 mJy beam −1 channel −1 ). This gives a source density of 0.22 ± 0.08 arcmin −2 GHz −1 in the field. Given that our ALMA observations cover an area of 0.88 arcmin 2 where PB ≥ 0.5 and are ∼20% shallower, one would expect to identify less than 0.36 ± 0.14 sources at S/N > 6 over the 1.875 GHz bandwidth of a baseband, and only a third of these (i.e., < 0.12 ± 0.05 sources) are expected to fall within ±0.3 GHz (1000 km s −1 ) of the main galaxies to be considered as companions. In other words, one would need to increase our survey area by >8× to detect a chance "companion" galaxy of the SMG or QSO2 in the field. Yet we have detected one S/N > 6 companion within 30 ′′ of each main galaxy. Our result thus indicates that both the SMG and QSO2 inhabit overdense environments, which is consistent with their purported large halo masses (Hickox et al. 2012).
In addition to the companion galaxies detected in CO emission, both the SMG and QSO2 are also associated with absorbers of high H I column density in the spectrum of a common background QSO (QSO1), as we will show in the next section.

ABSORPTION-LINE SYSTEMS
The two QSOs in the GAMA J0913−0107 system were first identified in the Sloan Digital Sky Survey (SDSS) DR9 quasar catalog (Pâris et al. 2012). The pair has a separation of only 10.8 ′′ , and more importantly, two closely separated DLAs were immediately identified in the low-resolution SDSS spectrum of QSO1 (Noterdaeme et al. 2012a). The stronger DLA (log N HI ≃ 21.3) at z abs ≈ 2.75 is associated with QSO2 at z = 2.7488, providing an important probe of the CGM around QSOs at an impact parameter of R ⊥ = 85 kpc (see Appendix C). The other DLA (log N HI ≃ 20.5) at z abs ≈ 2.68 provides a window to probe the CGM of the SMG at z = 2.674, which is just 11.7 ′′ from the QSO.
We searched the spectral databases with specdb 12 and found that the QSO pair had accumulated an excellent set of spectroscopic data from Gemini Multiobject Spectrograph (GMOS; Prochaska et al. 2013b), VLT/X-shooter (Finley et al. 2014), and Magellan Echellette Spectrograph (MagE; Rubin et al. 2015). Finley et al. (2014) noticed the strong coincident absorption at z abs ≈ 2.68 in the SDSS spectra of both QSOs, which motivated them to obtain the higher resolution X-shooter spectra for a detailed analysis. The absorption structure toward both QSOs is resolved into three major subsystems of variable metallicities and with a total velocity span of >1700 km s −1 . The observed kinematic and metallicity coherence across sightlines is remarkable, given the 86 kpc separation between the QSOs. The authors interpreted the system as a gaseous overdensity extended by six Mpc along the line-of-sight, which is suggestive of a clumpy filamentary structure that may eventually collapse and form a proto-cluster. They attributed the two main subsystems at lower velocities (A and B) as part of the IGM because of their low metallicity ([Fe/H] < −1.9) and suspected that the third main subsystem (C) with [Fe/H] = −1.1 is likely associated with a galaxy. Now with the detection of the SMG and its companion galaxies, we will use the coincident absorptionline system to characterize the CGM of these galaxies in § 5. We will show that subsystems A and B are cool gas streams in the CGM of the SMG, and subsystem C is indeed associated with a galaxy (Comp b).
In this section, we present a re-analysis of the z abs ≈ 2.68 absorption system using a new reduction of the X-shooter spectra ( § 4.1). Finley et al. (2014) used vpfit 13 to fit Voigt profiles to the entire spectrum. Our approach is complementary to the vpfit analysis and our results show a good agree-12 https://github.com/specdb/specdb 13 https://people.ast.cam.ac.uk/~rfc/vpfit.html ment with those presented in Finley et al. (2014). The main differences between the two analyses are: 1. We fit Voigt profiles to the H I Lyman series after masking out contaminating LYAF lines and quantify the statistical and systematic uncertainties of the model parameters using an MCMC algorithm ( § 4.2).
2. We measure the ionic column densities of metals with the apparent optical depth method (AODM; § 4.3). This is a more direct and conservative technique compared with Voigt profile fitting using vpfit, because it relies only on equivalent width measurements and uses a straightforward method to detect line saturation.
3. We use ionic column density ratios to constrain the photoionization model for each cloud, which in turn provides the ionization correction factors necessary for metallicity estimates ( § 4.4). (1−2.5 µm). The total exposure times are 100 min for QSO2 and 310 min for QSO1. We downloaded the raw data from the ESO archive and reduced the data with the spectroscopy data reduction pipeline developed by George Becker 14 . The final 1D spectra were corrected for the 0.2Å (0.5 pixel) wavelength redshift of the spectra in the VIS arm noticed by Noterdaeme et al. (2012b), which is likely produced by uncompensated instrumental flexure.
We fit the QSO continuum using the Python software package linetools 15 . First, the spectrum is divided into a number of wavelength intervals, which are ∼50 Å wide shortward of the Lyα emission, narrower across strong emission lines, and wider in regions free of emission lines and longward of Lyα. Next, a spline is fit through the central wavelength and median flux of each interval (i.e., the spline "knots"). Finally, these "knots" are iteratively added, deleted, or moved until a satisfactory continuum fit is obtained. Fig. 5 shows the H I absorption profiles (Lyα through Lyδ) of the absorbers at z abs ≈ 2.68 toward the two QSOs. Although the two QSOs are separated by 10.8 ′′ (86 kpc at z = 2.674), their H I absorption profiles show strikingly similar velocity structures spanning over 1800 km s −1 , as first noted by Finley et al. (2014). The kinematic coherence indicates that the medium responsible for the absorption is extended at least 86 kpc across the sky plane. Line blending from other absorbers is evident, as indicated by the disagreement in velocity profile among the Lyman series. To measure H I column densities in such a complex situation, it is beneficial to first identify a guessed solution by iteratively varying the Voigt profiles (convolved to R = 6400) until an acceptable fit to the data is obtained. The guessed solution not only provides a good starting point for the formal minimum χ 2 approach below, but also helps to identify regions contaminated by line blending (which thus should be flagged out). During this procedure, we find that a minimum of 10 clouds are needed to adequately fit the Lyman series in each sightline. Because each cloud is described by three parameters (v, b, log N HI ), our model for each QSO spectrum has a total of 30 free parameters.

Voigt Profile Fitting of Neutral Hydrogen
To model the absorption toward QSO1, the H I Lyman lines of the DLA at z abs ≈ 2.751 (see Fig. C1a) must be included in the model because its Lyα and Lyδ blend with the Lyα and Lyγ profiles of the absorber at z abs ≈ 2.68. We find that the DLA's H I absorption is adequately modeled as two clouds separated by 290 km s −1 (z abs = 2.7502, 2.7538), each with log N HI = 21.0 and b = 40 km s −1 (see Fig. C2). These parameters for the DLA at z abs = 2.751 are fixed in the fitting process.
The χ 2 minimization is focused on the velocity range between −1500 and 2100 km s −1 for QSO1 and between −600 and 1350 km s −1 for QSO2. With the guessed solution, we also mask out the pixels that are clearly contaminated by line blending within the fitting ranges. The surviving "good" pixels are indicated by diamond symbols with error bars in Fig. 5 and the Voigt models are optimized using the AMOEBA + MCMC method described in § 1.4. The priors of central velocities and column densities are centered around the guessed solution, with bounds of ±100 km s −1 for v HI and ±0.8 dex for log N HI . On the other hand, the Doppler parameter, b HI , is allowed to vary between 5 and 70 km s −1 . For three H I components, we found it necessary to fix their velocities to those measured from low-ion metal lines, because the H I series alone do not constrain their velocities well. Specifically, these components are at −470 and −225 km s −1 toward QSO1 and at −209 km s −1 toward QSO2. The optimized models are plotted against the data as blue curves in Fig. 5 and the formal parameters and their statistical uncertainties are tabulated in Table 4.
Because of the empirical nature of our placement of the unabsorbed QSO continuum, the Voigt parameters suffer from significant systematic uncertainties. In particular, we are interested in the systematic uncertainties of log N HI , which depends on (1) the column density due to the varying gradient of the curve of growth, (2) the quality of the spectrum, and (3) the significance of line blending. To quantify this, we run the same modeling procedure as above but vary the QSO continuum model by ±10% and use the resulting offsets between the three formal solutions to estimate systematic uncertainties. All subsequent errors in log N HI and metallicities ([X/H]) include both statistical and systematic uncertainties. Fig. 5 reveals that there are three separate kinematic clumps centered around δv ≈ −300, +400, +1200 km s −1 relative to   Table 6 are highlighted in color. Vertical dotted lines strike out regions that are blended with lines from absorbers at other redshifts (see Appendix C). The error spectrum is plotted (blue) when it shows significant structures.  1526.7070 −0.8761 · · · · · · Si IV 1393.7602 −0.2899 33.49 45.14 · · · 1402.7729 −0.5952 · · · · · · Fe II 1608.4508 −1.2388 7.90 16.20 · · · 2382.7642 −0.4949 · · · · · · · · · 2600.1725 −0.6209 · · · · · · NOTE. -The columns λrest and log f list the rest-frame wavelengths and the oscillator strengths (Morton 2003). The columns IP 0 and IP 1 list the ionization potentials (IPs) to create the ion from the immediate lower state and to ionize it to the immediate higher state, respectively. z SMG = 2.674 (i.e., z abs ≈ 2.6703, 2.6789, 2.6887), which we designate as subsystems A, B, and C, respectively, following the nomenclature of Finley et al. (2014). The same clumps appear toward both QSOs, although their column densities vary between sightlines. Half of the six subsystems are optically thick (i.e., log N HI > 17.2), including two sub-DLAs (QSO1-A and QSO1-C) and one LLS (QSO2-A). Metal absorption lines from these subsystems are thus expected, as we will show in the next subsection. Fig. 6 compares the velocity profiles of H I Lyβ and a selection of metal line transitions commonly observed in LLSs and DLAs (see Table 5 for transition data). We find that five of the six H I subsystems are detected in at least one metal transition; the only exception is QSO1-B. Similar to their H I absorption, the metal-line absorptions of subsystems QSO1-A and C, the two sub-DLAs, are resolved into multiple com-ponents. Note that the X-shooter spectrum has higher resolution for metal lines in the range 1500 < λ rest < 2700 Å (R = 11000 or FWHM = 27 km s −1 ) than for the H I Lyman lines at λ rest < 1217 Å (R = 6400 or FWHM = 47 km s −1 ). For each distinct metal-line cloud, we define a velocity integration window (highlighted in Fig. 6) and name it by adding a number suffix to designate its associated subsystem. The cloud QSO1-C1 shows the strongest metal absorption with at least four blended components within ∼200 km s −1 . We treat it as a single entity here, because for the purpose of measuring the gas metallicity it is unnecessary to deblend these components with Voigt profile fitting. Because the absorption toward QSO2 spans a narrower velocity range than that toward QSO1, the former is missing the most blueshifted cloud "A1" and the most redshifted cloud "C2". For completeness, we defined QSO1-B1 based on its H I absorption because no metal lines are detected there. As a result, there are a total of eight metal-line clouds. The top section of Table 6 lists the velocity integration windows of the clouds and their H I column densities by summing log N HI of the Voigt components within the velocity windows. Each cloud contains only one H I Voigt component except QSO1-B1 and QSO2-A2, both of which contain two closely separated components.

Ionic Column Densities from the AODM Method
In Appendix D, we provide an overview of the AODM method and our measurements of ionic column densities from all of the selected transitions (Table D1). The listed uncertainties of the unsaturated and unblended detections in the Table  include both statistical and systematic errors. Column densities from the AODM method are taken as lower limits for lines with more than one saturated pixels (which we define as I obs /I 0 ≤ 0.05 or τ ≥ 3) and are taken as upper limits for lines that are blended with transitions from absorbers at other redshifts. Lastly, for undetected transitions, we quote 3σ sta upper limits on the column densities. Systematic errors are not used here because it's not meaningful to adjust the QSO continuum around an undetected transition.

Ionization Correction and Metallicities
A relative metallicity measurement of the intervening gas requires (1) H I column density, (2) ionic column density of a metal element, (3) the reference solar abundances, and (4) the ionization correction. The definition of the relative metallicity makes this explicit: where X i denotes the ionic state i of element X, f Xi ≡ N Xi /N X is the fraction of the element in the ionic state i, f HI ≡ N HI /N H is the neutral fraction of Hydrogen, [X/H] ′ ≡ log(N Xi /N HI ) − log(N X /N H ) ⊙ is the raw metallicity, and IC ≡ log f HI − log f Xi is the ionization correction. We have obtained the first two items (log N HI and log N Xi ) in § 4.2 and § 4.3 and the results are listed in Tables 4 and  D1. Combined with the elemental abundances of the presentday solar photosphere from Asplund et al. 2009, we are ready to calculate the raw metallicity [X/H] ′ . Next, we calculate the ionization correction (IC) using CLOUDY photoionization models (Ferland et al. 2017).
1 Ryd at the illuminated face. The former has been measured, but the latter needs to be constrained by comparing the column density ratios of different ionic states of the same elements and predictions from photoionization models. Therefore, for each cloud listed in Table 6, we calculate a set of photoionization models, with a termination condition set to meet the observed H I column density. For the ionizing source, we use the Haardt & Madau (2012) radiation background interpolated to z = 2.67 with contributions from both galaxies and quasars. For the cloud, we assume plane-parallel geometry, the solar relative abundance pattern, a metallicity of [M/H] = −1.5 16 , and a range of hydrogen volume densities (1.09 ≥ log n H /cm −3 ≥ −4.91) to cover ionization parameters between −6 ≤ logU ≤ 0. For each cloud, we compare the observed ionic column ratios (C IV/C II and Si IV/Si II) and the model-predicted ratios to constrain the ionization parameter. We list the constraints on logU in Table 6. It is rare to have detections of both high ions and low ions in the same cloud, leading to many upper or lower limits of logU. But we found that the most plausible ionization parameters lie between −4 logU −2, comparable to other published LLSs (e.g., Prochaska 1999;Lehner et al. 2016). Depending on the data constraint, we adopt logU values of −3.5, −3.0, or −2.0 for each cloud. Finally, once both log N HI and logU are fixed, we use the CLOUDY model of the same parameters to calculate the ICs for all of the ions (Table D2), which are then used to obtain the ionization-corrected metallicity measurements for all of the transitions (Table D3). Notice that the ionization correction only becomes important for (1) low ions in clouds with log N HI 19 and (2) intermediate or high ions (e.g., C IV and Si IV) at all column densities. For low ions in sub-DLAs, the ICs are fairly small (±0.15 dex). We also use the model to infer the total H column density (log N H ), the H neutral fraction (log f HI ), the H volume density (log n H ), and the characteristic line-of-sight depth of the cloud (log l = log N H − log n H ). The results are listed in Table 6. 16 The derived ICs are insensitive to the assumed metallicity.
The best metallicity measurement is provided by O I, because O I has the smallest ionization correction factors due to its charge-exchange reactions with hydrogen (Field & Steigman 1971) and its hydrogen-like ionization potential. Unfortunately, O I λ1302 is saturated in C1 toward QSO1 (a common issue of the transition for DLAs) and is undetected in cloud B toward QSO1 and the clouds toward QSO2. As a result, transitions from other ions need to be used. The bottom section of Table 6 lists the final adopted metallicities from our preferred α-element transitions and Fe II. Note that the ICs for QSO2-B1 and QSO2-C1 are large because only the C IV lines are detected in these clouds. We fail to obtain a reliable metallicity measurement for QSO1-B1 because of the absence of metal absorption. For the four clouds where we have both [α/H] and [Fe/H], we found a moderate level of α-enhancement ([α/Fe]) between 0.15 and 0.61 (with an inverse-variance-weighted mean of 0.2), comparable to those previously measured in z > 2 DLAs ([α/Fe] = 0.30 ± 0.16; Rafelski et al. 2012).
We opt not to correct the gas-phase metallicity for dust depletion, because (1) little depletion is expected from volatile elements such as O and C, (2) the SDSS spectra and photometry of the QSOs show no evidence of significant dust reddening, and (3) the depletion factors are largely uncertain in external galaxies. As a reference, in the Milky Way's ISM, volatile elements (e.g., C, N, O, S, Zn) show depletions of (X/H) ISM − (X/H) gas 0.3, while refractory elements (e.g., Mg, Al, Si, Fe, Ni) show 0.7 (X/H) ISM − (X/H) gas 2.0 (Savage & Sembach 1996;Groves et al. 2004). These local measurements from a metal-rich ISM provide strict upper limits on the level of depletion expected in the CGM of the SMG.

EMISSION-ABSORPTION CONNECTION
Identifying the emission counterparts of the intervening gas helps us compare the properties of the galaxies to those of their CGM. Having analyzed the ALMA CO emitters in § 3  (Table 4), (d) Al IIλ1670.8 absorption, and (e) metallicities of metal-line-defined clouds ( Table 6). The absorbers toward QSO1 and QSO2 are color-coded in red and blue, respectively. In (a), the SMG spectrum has been divided by 4× to show it together with the spectra of its companions, and the vertical dashed lines indicate the centroid velocities of the CO emission lines at 0, 58, 1171, and 1447 km s −1 . All velocities are relative to the redshift of the SMG (z SMG = 2.674) and the gray shaded regions indicate velocities beyond the escape velocity of a 10 13 M ⊙ halo (vesc ≃ 700 km s −1 ). and the QSO absorption spectra in § 4, we are now ready to draw connections between the emission and the absorption based on proximity in both spatial and redshift dimensions. Fig. 7 directly compares the absorption profiles from the QSOs to the CO emission profiles from the SMG and its companions. The H I Lyβ and Al IIλ1670.8 profiles from the two QSOs are plotted together to illustrate the striking kinematic coherence. In addition, the figure illustrates log N HI from the Voigt profile solution in Table 4 and the ionization-corrected metallicities of the metal-line-defined clouds in Table 6. In § 4, we have found that the z abs ≈ 2.68 H I absorbers have total H I column densities of log N HI = 20.53 +0.06 −0.06 and 18.60 +0.23 −0.43 toward QSO1 and QSO2, respectively. Each absorber is resolved into three main subsystems (A, B, and C) with velocity spans of ∼1500-2000 km s −1 . Although their H I column densities vary significantly between the two QSO sightlines, their radial velocities show remarkable consistency, indicating that the QSOs are intercepting three expansive sheets/filaments of gas. At the same time, the extreme velocity widths of the absorption-line systems suggest that they probe merging systems .
Results in Fig. 7 show that subsystem C is unlikely to be in the same halo as subsystems A, for several reasons. First, subsystem C is 10× more metal-enriched than subsystem A ([M/H] ≃ −1.1 vs. −2.1). Secondly, the velocity spans of ∼1950 km s −1 (QSO1) and ∼1460 km s −1 (QSO2; Table 4) and their asymmetric distributions around z SMG (absorption is centered at z abs = 2.68) are inconsistent with gravitational motions inside even a 10 13 M ⊙ halo centered on the SMG, because the escape velocity of such a halo following the NFW profile is flat at ∼700 km s −1 between 60 kpc and the virial radius of 186 kpc. Lastly, we have detected CO emission at almost exactly the redshifts of the absorbing clouds in subsystem C in Comp b, which lies much closer to the QSOs than the SMG. Therefore, we consider subsystem A part of the CGM of the SMG at z = 2.674, and subsystem C part of the CGM of Comp b at z = 2.6884 and 2.6917.
As for subsystem B (log N HI ≃ 16), although its velocity allows an association with the SMG, it is unimportant because its contribution to the CGM is negligible compared to subsystem A.
Once the emission counterparts are determined, the absorption-line measurements from the two QSO sightlines can be plotted against the impact parameter to show crude radial profiles of the CGM around the SMG and Comp b. We consolidate the results in Table 7, where we have assigned a velocity window for each subsystem that captures its major clouds. We calculate the total H I column densities from the H I Voigt components within these velocity windows. When there are multiple metal-line clouds in a subsystem, the metallicity is calculated as the N H -weighted mean [α/H], where N H is the H I + H II column density based on the adopted CLOUDY model.
The GAMA J0913−0107 system gives us the first glimpse into the CGM around an SMG. With high H I column density and low metallicity at large impact parameters, the CGM of the SMG is distinct from the CGM of QSOs and normal starforming galaxies at z ∼ 2 − 3.
The profile of H I column density is shown in Fig. 8a. We compare our measurements with literature QSO absorptionline measurements in the surroundings of z ∼ 2 − 3 QSOs (Lau et al. 2016) and Lyman Break Galaxies (LBGs; Simcoe et al. 2006;Rudie et al. 2012;Crighton et al. 2013Crighton et al. , 2015. The H I column densities of the SMG's CGM, similar to coeval QSOs, are significantly greater than those of star-forming galaxies at R ⊥ 70 kpc. The H I column density declines as we move away from the SMG, with a gradient of −2.0 ± 0.4 dex per  Rudie et al. (2012), and Crighton et al. (2013Crighton et al. ( , 2015, and SFG [M/H] data are from Simcoe et al. (2006). The horizontal dashed lines in the bottom panel indicate the range of IGM metallicities measured from the LYAF at z abs ∼ 2.5 ([M/H] = −2.85 ± 0.75; Simcoe et al. 2004). The sloped line shows the oxygen abundance profile of the giant spiral galaxy M101 measured with an electron-temperature-based method (Eq. 5 of Kennicutt et al. 2003). The solid portion is covered by H II regions, while the dotted portion is an extrapolation. 100 kpc. How do the observed column densities compare with the projected surface mass density Σ M (R) of NFW halos? To make this comparison, we first calculate Σ M (R) by integrating the NFW density profile ρ(r) up to the virial radius: We then convert it to H I column densities by assuming a baryon−dark-matter density ratio of  Fig. 8a shows that the H I column density profiles of the SMG and the QSOs are consistent with the expectation from a ∼10 13 M ⊙ halo. This is in agreement with the de-17 Comparable to the neutral fractions estimated by our photoionization models in Table 6. tection of a high column of H I (log N HI = 18.6) at an impact parameter as far as 176 kpc. The agreement also shows that the neutral gas in the absorption-line systems can account for ∼10% of the total baryonic mass in the halo if they have a filling factor close to unity. On the other hand, the log N HI profile of LBGs is more consistent with a ∼10 12 M ⊙ halo.
The metallicity profile is shown in Fig. 8b. The literature data points show that, metal-enriched gas with −1 [M/H] 0 dominates the inner part (R ⊥ 150 kpc) of the CGM around both QSOs and LBGs. By contrast, the CGM of the SMG is poor in metal, with almost a constant metallicity of [M/H]≈ −2 across two sightlines separated by 86 kpc (10.8 ′′ ). Its metallicity is near the 1σ upper bound of the metallicities measured in the LYAF (i.e., the IGM; Simcoe et al. 2004), and it is lower by ∼1.5 dex at R ⊥ = 93.1 and by ∼0.5 dex at R ⊥ = 175.5 kpc than that of the CGM of QSOs.
On the other hand, the CGM of Comp b show properties similar to that of normal star-forming galaxies at z ∼ 2 − 3. Both line emitters in Comp b have orders of magnitude lower molecular gas mass than the SMG (see Table 1). Assuming a CO-to-molecular-mass correction factor of α CO = 4.3 and a CO excitation correction of r 31 = 0.52, the emitters at z = 2.6884 and 2.6917 have molecular gas masses of M mol = 5.6 × 10 9 and 4.4 × 10 9 M ⊙ , respectively. The gas masses are comparable to those measured in the lensed normal starforming galaxies that have stellar masses of ∼ 10 10 M ⊙ (Saintonge et al. 2013). One would thus expect a CGM similar to that of LBGs. The corresponding subsystem C provides absorption-line measurements at R ⊥ = 32.3 and 58.9 kpc. It shows significantly metal-enriched gas (compared to the IGM level) with a large variation in H I column density over just a difference of 27 kpc in impact parameter. Comp b thus is surrounded by a metal-enriched clumpy medium extended to at least ∼60 kpc.
In the above, we have shown that the CGM of SMM J0913 is distinctly different from the CGM of QSOs and normal starforming galaxies. But how do our absorbers compare with other DLAs in terms of absorption-line properties only? The H I column densities of subsystems A and C toward QSO1 miss the DLA threshold of log N HI = 20.3 by merely ∼0.1 dex. Because such small differences are comparable to the measurement uncertainty, more liberal thresholds have been used to select DLAs (e.g., log N HI > 20.1 in Rubin et al. 2015). Combined with the result that the two sub-DLAs have neutral fraction of ∼50% from CLOUDY models, we believe that it is appropriate to compare their properties with those of the general DLA population.
With high-resolution spectra of 100 DLAs at z abs ∼ 2 − 4, Rafelski et al. (2012) found that their metallicity distribution is well fit by a Gaussian with a mean at [M/H] = −1.57 and a dispersion of 0.57 dex. The DLA metallicity only mildly decreases with redshift, but shows a strong correlation with the width of low-ion metal lines (e.g., ∆V 90 or W 1526 ) (Neeleman et al. 2013). This correlation between kinematics and metallicity is generally interpreted as a manifestation of the mass-metallicity relation (e.g., Tremonti et al. 2004;Erb et al. 2006), because the line width may reflect the halo mass (like in the Tully-Fisher relation), which in turn is proportional to the stellar mass (Møller et al. 2013).
Previous works have used higher-resolution spectra (R 40000) to measure ∆V 90 , the velocity interval including 90% of the optical depth of an unsaturated line. And the alternative kinematic parameter W 1526 , the rest-frame equivalent 10 100 9.-Metallicity−line-width relation of (sub-)DLAs at z ∼ 2 − 3. Gray squares are 44 DLAs between 2 < z abs < 3 compiled from the literature (Neeleman et al. 2013). The red and blue circles show respectively subsystems A and C toward QSO1. For both sub-DLAs, we have estimated the velocity widths using the separation of resolved components in the R ∼ 11000 X-shooter/VIS spectrum.
width of Si IIλ1526.7, is unsuitable for our sub-DLAs because the line is unsaturated (see Fig. 6). The equivalent width only becomes a good kinematics indicator when the line is saturated, weakening its dependency on the Si + column density (and consequently, on metallicity). So to place our sub-DLAs on the relation, we adopt the velocity separation between kinematically resolved substructures as a surrogate of ∆V 90 . QSO1-A shows two velocity components of similar strength separated by ∼250 km s −1 , and QSO1-C is dominated by the stronger C1 cloud, which appears to be a blend of four components with a velocity span of ∼200 km s −1 . These estimates of the velocity width should be considered as lower limits on ∆V 90 , because they are the separations between peak optical depths. Fig. 9 compares the two sub-DLAs against literature DLAs. To control the redshift evolution of the mass-metallicity relation, only the DLAs between 2 < z abs < 3 are plotted. While QSO1-C blends in the general trend established by DLAs, QSO1-A is a clear outlier. Firstly, very few DLAs have metallicities as low as QSO1-A: only 2 out of the 44 DLAs (4.5 +5.5 −1.5 %) have [M/H] ≤ −2.1. Secondly, its high velocity width places it significantly below the metallicity-linewidth relation of DLAs, which would have predicted a velocity width of only ∼20 km s −1 based on its metallicity at [M/H] = −2.1.
This finding suggests that most of the DLAs are likely closer to their hosts than our sub-DLA QSO1-A (R ⊥ = 93 kpc), consistent with previous observations of DLA host galaxies (see § 1.3). More importantly, the unusual combination of high velocity-width and low metallicity provides a method -based on absorption-line properties alone -to potentially separate (sub-)DLAs associated with normal starforming galaxies at low impact parameters with those associated with more massive galaxies like SMGs at large impact parameters.

SUMMARY & DISCUSSION
We have carried out a comprehensive study of the emissionabsorption system GAMA J0913−0107 at z ∼ 2.67 with a multi-wavelength data set obtained primarily from Herschel, ALMA, and VLT/X-shooter. The system consists of a bright SMG, its CO companion galaxies, and a number of optically thick H I absorbers toward two background QSOs within 22 ′′ of the SMG. Our main results are: 1. The Herschel-selected SMG at z = 2.674, with an 870 µm flux of 7.4 mJy and an IR luminosity of ∼ 10 13 L ⊙ , is one of the most luminous dusty star-forming galaxies. Its properties are similar to the general SMG population at z ∼ 2, featuring a short gas depletion timescale of ∼0.1 Gyr and compact (sub-arcsec) sizes in both dust emission and CO (3−2) emission. The high S/N spectrum reveals two CO (3−2) components at almost the same redshift: ∼1/4 of the line flux is in a broad component with a FWHM ≃ 900 km s −1 , while ∼3/4 of the flux in a narrow component with FWHM ≃ 250 km s −1 .
2. Three companion CO emitters are identified within 30 ′′ and 1500 km s −1 of the SMG. A comparison with the source counts from the ASPECS field survey indicates that the SMG lives in an over-dense environment.
3. Two nearby QSOs provide background beacons to probe the CGM of the SMG. A DLA with a total H I column density of log N HI = 20.5 is identified at z abs ∼ 2.68 in the closer QSO sightline at θ = 11.7 ′′ . The DLA is quite unusual, in terms of both the large impact parameter (R ⊥ = 93.1 kpc to the SMG) and the enormous velocity span (∼2000 km s −1 ). X-shooter resolved the DLA into three major subsystems, including two sub-DLAs with distinctly different metallicities separated by ∼1600km s −1 . Remarkably, the same subsystems are also found in the farther QSO sightline at θ = 22.1 ′′ : they have nearly the same velocities and metallicities as their counterparts at θ = 11.7 ′′ , despite lower H I column densities (total log N HI = 18.6).
4. We use the absorption-line systems to characterize the CGM of the SMG and its companion Comp b, and we compare their properties with the CGM of QSOs and normal star-forming galaxies. The CGM of the SMG forms its own category: while its high column densities at large impact parameters are similar to the massive halos inhabited by z = 2 − 3 QSOs, its metal content (∼1% solar) is an order of magnitude lower than the circum-QSO medium. On the other hand, the CGM of the much less luminous Comp b is more consistent with that of normal star-forming galaxies at z = 2 − 3: showing significant metal enrichment (∼10% solar) within R ⊥ 60 kpc.
The detection of high-column density, mostly neutral, metal-poor gas in the CGM of a massive dusty starburst galaxy at z = 2.674 has powerful implications to theories of galaxy formation and evolution. The remarkable consistency of the H I absorbers in both kinematics and metallicity across two sightlines separated by 86 kpc is at odds with CGM models that assume randomly floating H I clouds in pressure equilibrium with hot X-ray gas. Instead, it is logical to assume that the background QSOs have intercepted a single filament of cool gas permeating in the halo. Narrow filaments of cool gas and satellite galaxies can penetrate the hot CGM of massive halos without ever being shock-heated to the virial temperature, because (1) massive halos are rare and tend to form at the intersections of the cosmic web and (2) the cooling time is shorter in the filaments than in the halo (Dekel & Birnboim 2006;Dekel et al. 2009;Kereš et al. 2009). At z 2.5, such cold streams can survive even in halos more massive than ∼10 13 M ⊙ (although note that this mass limit is highly uncertain). In particular, stream-feeding is likely important in the bright SMGs with S 850 > 5 mJy because their comoving space density exceeds the expectation from minor and major mergers (Dekel et al. 2009).
The observed properties of the absorption-line systems match the simulation-predicted properties of cold streams. First, the simulations of Dekel et al. (2009) show that the inflow velocity is comparable to the virial velocity and is roughly constant along the filament. This is consistent with the velocity shift (δv = −300 km s −1 ) and the kinematic coherence we saw between the clouds in both QSO sightlines. Secondly, by post-processing gas in seven simulated halos with M halo = 10 10 − 10 12 M ⊙ between 1.5 < z < 4.0, Fumagalli et al. (2011) found that the absorption-line systems associated with the smooth stream component have systematically lower metallicity (∼1% solar). This is exactly the level of the gas metallicity we measured in the absorbers associated with the SMG.
Radially inflowing on nearly a free-fall timescale, the cold streams may account for the bulk of the baryonic accretion rate and become the dominant mechanism to feed the growth of galaxies (Kereš et al. 2009). We can crudely estimate the gas accretion rate from the filament that the QSOs intercepted. The filament has a length >176 kpc and a depth on the order of 10 kpc at R ⊥ = 93 kpc. The former is estimated from the impact parameter of QSO2, and the latter is estimated from the ratio between the column density and the volume density of Hydrogen, l = N H /n H , inferred from the photoionization model of the sub-DLA QSO1-A2. The depth-to-distance ratio is 0.11 radian or 6 • , comparable to the opening angles of 20 − 30 • seen in simulations (Dekel et al. 2009). Our photoionization models also indicate similar H I+H II column densities for QSO1-A2 (log N H = 20.4) and QSO2-A2 (log N H = 20.1), the two main clouds associated with the SMG at R ⊥ = 93, 176 kpc. By assuming an opening angle of β = 25 • , an average hydrogen column density of log N H = 20.2, and a 10 13 M ⊙ NFW halo at z = 2.674, we estimate that the mass of the filament is M fil = f He m p N H (β R 2 vir /2) ∼ 1.3 × 10 10 M ⊙ . Given an accretion timescale of τ acc = R vir /V vir = 4 × 10 8 yr, the gas accretion rate from this single filament is ∼33 M ⊙ yr −1 . Typically three main filaments are seen in the simulations, so our estimate shows that cold-mode accretion can supply gas at a rate of ∼100 M ⊙ yr −1 . Although accounting for only ∼10% of the current SFR, our estimated cool gas accretion rate is in fact comparable to the rate of total gas supply to the central galaxies in 10 13 M ⊙ halos at z = 2 from cosmological simulations (see Fig. 9 of Kereš et al. 2009), and at this rate the molecular gas reservoir of 10 11 M ⊙ can be acquired in just ∼1 Gyr. On the other hand, star formation at the current intensity seems unsustainable despite the efficient gas supply from cold streams.
In this work, we have presented the first observational evidence that supports the existence of cold streams in the CGM of a massive starburst galaxy. The GAMA J0913−0107 system has an excellent data set and the results are highly informative, but larger samples are clearly desired to draw conclusions on the general properties of the CGM. We hope that this will serve as a springboard for upcoming statistical studies of the CGM in similar galaxies. As an attempt to inform these future studies, it is worth discussing the major technical challenges we had faced in this program: 1. The large beam of Herschel (17.8 ′′ at 250 µm) makes it inefficient to identify SMG-QSO pairs with small angular separations (θ 10 ′′ ). As a result, QSO1 in GAMA J0913−0107 is the only one that probes below 100 kpc; and despite intercepting a sub-DLA it has yet to expose the chemically enriched area of the CGM of the SMG.
2. High S/N spectra with moderately high spectral resolution are needed to unambiguously detect optically thick absorbers with 17.2 < log N HI 19. For example, with the R ∼ 2000 SDSS spectrum of QSO2, we couldn't have identified the LLS associated with the SMG (i.e., QSO2-A2 with log N HI = 18.6). But the LLS is unambiguously detected in the X-shooter spectrum because of its resolved H I Lyman profiles and the clear detection of low-ion metal lines. Similarly low column densities are expected in most of our sample, because all of the other spectroscopically confirmed SMG-QSO pairs we have so far have impact parameters between 100 < R ⊥ < 300 kpc (Fu et al. 2016, and unpublished data) and none of them show obvious (sub-)DLA features (which may be expected given the large impact parameters).
3. It has been difficult to obtain spectroscopic redshifts of the Herschel sources because (1) they require subarcsec positions from interferometers to place spectroscopic slits, (2) the near-IR spectral range suffers from heavy telluric absorption and OH emission, and (3) SMGs tend to be weak in rest-frame optical lines. The latter two points are the main reasons why our redshift success rate is only ∼60%.
Possible solutions to these issues may be already on the horizon. To address the first challenge, we need to design an efficient interferometer imaging survey, because a sample of Herschel sources with less than 10 ′′ apparent offsets from optical QSOs is likely overwhelmed by contaminating sources. Two third of the Herschel-SDSS-selected pairs with apparent separations between 5 ′′ and 10 ′′ turned out to be IR-luminous QSOs instead of projected SMG-QSO pairs. A good strategy is to observe multiple sources located within a ∼10 • diameter circle in a single ALMA scheduling block (SB). In Fu et al. (2017), we managed to observe ∼10 targets in a single ∼50 min SB, achieving an on-source integration time of 200 s per source and an rms of 0.12 mJy beam −1 . Do we have enough such pairs to populate a 50-min SB? The surface density of Herschel-QSO pairs with offsets less than 10 ′′ is 0.16 deg 2 (26 pairs in the 161.6 deg 2 H-ATLAS GAMA fields), which gives ∼13 targets for a ∼10 • diameter circle.
To address the second challenge, we need QSO spectra with quality similar to the X-shooter spectra presented here to better sample the spatial profile of the CGM. The QSOs in our sample are selected to have g < 22, the majority of which requires ∼1.5 hours' exposure time with an Echellette spectrograph on a 10-m-class telescope to reach a sufficient S/N at R ∼ 8000 (e.g., see the Keck/ESI survey of DLAs by Rafelski et al. 2012).
To address the last challenge, we need a more efficient method to obtain SMG redshifts. One potential approach is to exploit the frequency scan mode offered by modern millimeter interferometers. This method has the additional advantage of bypassing the interferometer imaging step because the primary beam is usually larger than the Herschel positional uncertainty and the line detection also provides positional information. For instance, with NOEMA scans of the 2 mm and 3 mm bands (only two spectral configurations per band), Neri et al. (2020) obtained 12 secure redshifts for 13 sources. The average telescope time spent is ∼105 min per source, including ∼40 min overhead. Although the targets in this Pilot study have 500 µm flux densities greater than 80 mJy (many are strongly lensed), this technique could be applied to fainter sources (like ours with S 500 > 20 mJy) as the instrument sen-sitivity and overheads continue to improve.
We thank D. Kereš  The source is detected in all of the nine photometric bands included in the KiDS+VISTA photometry catalog, with r = 21.58 ± 0.02 and K s = 20.84 ± 0.18. Our best-fit stellar population synthesis model of the photometry reveals a stellar mass of ∼ 3 × 10 7 M ⊙ and an SFR of ∼ 0.006 M ⊙ yr −1 . We have used the Bruzual & Charlot (2003) models assuming exponentially declining star-forming histories and the Chabrier (2003) initial mass function.
Could the SMG be gravitationally magnified by this foreground dwarf galaxy? Using the halo-mass−stellar-mass relation from abundance matching (Fig. 6 of Bullock & Boylan-Kolchin 2017), we estimate a halo mass of ∼ 3 × 10 10 M ⊙ . The halo mass corresponds to a line-of-sight velocity dispersion (σ) of just ∼29 km s −1 , assuming a singular isothermal sphere (SIS; ρ(r) = σ 2 /2πGr 2 ) and the fitting function of the halo overdensity ∆ c (z) from Bryan & Norman (1998). The Einstein radius of the SIS can be calculated as (e.g., Kneib & Natarajan 2011): where c is the speed of light, D ds is the angular diameter distance between the deflector and the background source, and D s is the angular diameter distance between the observer and the background source. For our system consisting of a foreground lens at z d = 0.055 with σ = 29 km s −1 and a background source at z s = 2.674, the Einstein radius is r E ∼ 0.024 ′′ . Because r E is 33× smaller than the 0.8 ′′ offset between the SMG and the foreground galaxy, we conclude that the foreground galaxy is unlikely to have any measurable lensing effect on the SMG.
B. BLIND LINE DETECTION IN THE ALMA BAND-3 DATA To search for faint emission-line sources with unknown line widths in 3D data cubes, a matched-filtering algorithm is commonly used: e.g., in SoFiA (Serra et al. 2015), FindClump (Walter et al. 2016), LineSeeker (González-López et al. 2019, and MF3D (Pavesi et al. 2018). We wrote an IDL program to implement this simple algorithm. The convolution kernel we chose to filter the data in the spectral dimension is a top-hat function with a variable half-width between n = 1 and n = 9. For each channel, the data in the neighboring ±n channels are stacked with equal weighting. Given the average channel spacing of ∼25 km s −1 , the convolved channel widths range between ∼75 km s −1 and ∼475 km s −1 . As shown in González-López et al. (2019), the simple top-hat function is as effective as the more sophisticated Gaussian kernels in detecting low S/N line emission.
In each convoluted channel map, we measure the rms noise level with a robust sigma routine and detect unresolved sources near the highest S/N pixel. An elliptical Gaussian fixed to the shape of the core of the dirty beam is fit to the 8 ′′ ×8 ′′ subregion centered on the pixel and subtracted from the image. The parameters of the best-fit Gaussians are saved at each iteration to form the raw source catalog. The iterative process continues until the image contains no pixels above the S/N threshold of S/N pix,th = 4.5. It is worth noting that this source-detection algorithm is similar to the minor cycles of the CLEAN deconvolution algorithm, but here we subtract only the core of the dirty beam to save computing time. This simplified approach is justified by the low S/N of the sources other than the SMG.
Given that a single source can be detected in multiple channels and with multiple convolution kernels, we remove duplicated detections by iteratively looping through the raw source list from the highest to the lowest S/N and discard all detections within 2 ′′ and 0.2 GHz from the highest S/N source remaining in the list.
Because our search is restricted to point sources, the source S/N simply scales with the ratio between the peak of the best-fit Gaussian (S peak ) and the rms noise of the convolved channel map: where a scaling factor is used to account for fitting errors (Eq. 9 of Rengelink et al. 1997, see also Condon 1997. To estimate the fidelity of the detected sources, we search for sources in simulated noise-only interferometer data instead of using negative "sources" in the actual data (e.g., González-López et al. 2019). First, we introduce random thermal noise by replacing the calibrated MS's visibilities with a normally distributed random array of complex numbers generated with numpy.random. This is equivalent to the CASA simulator function setnoise in the "simplenoise" mode, but our approach is faster because it only writes the DATA column once and does not add MODEL_DATA and CORRECTED_DATA columns to the MS. Unlike the fixed random number seed (11111) adopted in the CASA simulator, each noise realization uses a different seed in our code. We set the widths of the normal distributions for the real and imaginary parts to the standard deviations of the original visibilities measured with visstat (σ ∼ 250 mJy visibility −1 with a slight dependence on the spectral window). Then, we use tclean to image the noise-only visibilities into spectral data cubes with natural weighting. The same tclean parameters are adopted except that we turn off de-convolution by setting niter = 0, because the simulated data contain no sources. For each spectral window, we generate a set of 10 simulated data cubes to provide enough source counts for the source fidelity calculation.
We run the same line search code on the simulated data with the same detection parameters. We compare the normalized cumulative distribution functions (CDF) of the detected sources in the actual data and in the simulated data to estimate the source fidelity. We use the following equation, similar to LineSeeker used in the ASPECS-LP survey (González-López et al. 2019): where F(≥ S/N | n) is the fraction of sources detected at or above the source S/N, with its detection kernel width n, and at any frequency channel in the spectral window. The subscript indicates whether it is from the actual data or the simulated noise-only data. F data is measured directly from the CDF, while F sim is obtained from the best-fit error function of the CDF to mitigate noise at high S/N due to low source counts. The fidelity is thus the probability that the detected source is not due to random noise. A source has a high fidelity near unity when F sim ≪ F data , and a low fidelity near zero when F sim ≈ F data . We identified a total of six emission-line sources with fidelity > 0.9 within the primary beam from the non-interpolated datacubes of all four spectral windows. We extracted the spectrum for each detection from the linear-interpolated and primary-beamcorrected data cubes with an elliptical aperture matched to the sizes of the restoring beams; i.e., we had assumed that the sources are unresolved. We measured the central frequency (ν 0 ), FWHM, and line flux with a single-Gaussian model and list the results in Table B1. Fig. B1 illustrates the distribution of the detected sources in the field, and Fig. B2 shows the zoomed-in version of the integrated intensity maps and their ALMA spectra.    (Chen et al. 2015;Chen et al. 2016). The X-shooter spectra confirm all of the previously known absorbers and reveal several additional absorbers toward both QSOs: z abs = 1.0855, 2.283, 2.9147 toward QSO1 and z abs = 0.886, 1.0865, 2.2525, 2.2845, 2.3065, 2.3445, 2.3595 toward QSO2. Fig. C1 shows portions of the X-shooter spectra to illustrate all of the major absorbers we have identified (8 toward QSO1 and 12 toward QSO2), omitting only the Mg II absorber at z = 1.501 toward QSO2.
The z abs ≈ 2.75 DLA toward QSO1 is apparently associated to QSO2 at z = 2.7488. The DLA has been previously analyzed as part of the Quasar Probing Quasar (QPQ) project using a lower-resolution GMOS spectrum, from which they measured log N HI = 21.3 ± 0.15 (Prochaska et al. 2013b) and rest-frame equivalent widths of 2.60 ± 0.05 Å for C IIλ1334.5 and 0.51 ± 0.05Å for C IVλ1548.2 . We obtained similar results using the X-shooter spectrum (Fig. C2). With Voigt profile fitting, we find that the H I Lyman series is adequately fit with two components separated by 290 km s −1 (z abs = 2.7502, 2.7538), each with log N HI = 21.0 and b = 40 km s −1 . We thus obtain a total column density of log N HI = 21.3 with an estimated systematic uncertainty of ∼0.1 dex. With the AODM method and ICs from a CLOUDY model with log N HI = 21.3, logU = −2.5, and the HM12 radiation background, we measure an ionization-corrected α metallicity of [C/H] = −1.2 ± 0.1 from C IVλ1550.8 and an iron metallicity of [Fe/H] = −1.6 ± 0.2 from Fe IIλ1608.5 (see Table C1). Given the impact parameter of R ⊥ = 85 kpc, this DLA fits nicely with the CGM profiles of z = 2 − 3 QSOs in Fig. 8 (Lau et al. 2016  NOTE. -All measurements were made using a velocity integration window between 0 and 550 km s −1 relative to zsys = 2.7488. We use "<" for upper limits due to non-detections, ">" for lower limits due to line saturation, and " " for upper limits due to line blending.  NOTE. -This table lists the measured column densities from each transition, where we use "<" for upper limits due to non-detections, ">" for lower limits due to line saturation, " " for upper limits due to line blending, and "..." for no measurement due to poor data quality. E. THE IDENTIFICATION OF THE EMISSION COUNTERPART OF SUBSYSTEM C In the deep r-band image from KiDS (5σ limit at ∼25 mag) shown in Fig. E1a, we labeled four faint optical sources within 7 ′′ of QSO1. Here we explore whether any of these sources is connected to the DLA at z abs ≈ 2.68 by examining their photometric redshifts and their ALMA spectra.
Three −0.06 . where the r− and H-band magnitudes are from the homogenized "Gaussian Aperture and PSF (GAaP)" photometry and z p are the 9-band photometric redshift estimates from the Bayesian photometric redshift code BPZ (Benítez 2000). The 68% confidence intervals of the photometric redshifts suggest that both sources are in the far foreground of the SMG SMM J0913 (z SMG = 2.674). The fourth object (D) is not in the catalog likely because its proximity to the bright QSO1. We measured its position directly from the image: R.A. = 09 h 13 m 39.05 s , Decl. = −01 • 07 ′ 06.5 ′′ .
We then extracted spectra at their optical positions from the ALMA band-3 datacube. For objects A, B, and D, we adopted elliptical apertures matching the synthesized beam size (1.7 ′′ ×1.3 ′′ at PA = 49 • ). We show these spectra in Fig. E1c. Even at the depth of our ALMA data (rms= 0.155 mJy beam −1 channel −1 in BB4), none of the sources show emission lines at a detectable level.
For Object C, we initially used an aperture centered on the optical position and detected hints of emission lines at the expected frequencies of absorption-line clouds C1 and C2. We then made a CO map by combining the two channels that show the most significant emission. The CO image in Fig. E1b led to the discovery of Comp b as it reveals a highly significant source ∼3 ′′ to the SSW of the optical position, which we have designated as Comp b. Guided by the CO image, we re-extracted a spectrum from a 3.4 ′′ ×1.8 ′′ elliptical aperture that matches the geometry of Comp b to optimize the line detection. This spectrum is shown in the Object C panel of Fig. E1c. Line emission is clearly detected at the expected frequencies of absorption-line clouds C1 and C2 toward QSO1. Through this exercise, we have identified Comp b as the most likely emission counterpart of absorption subsystem C toward both QSOs. .6676 GHz, where CO emission from Object C is detected. The contours are drawn at −3, −2 (dashed), 2, 3, and 4σ (solid). In both images, QSO1 sets the origin of the coordinates. In the four panels in c, we show the ALMA spectra of the objects. The vertical dashed lines indicate the CO (3−2) frequencies that correspond to the redshifts of the major absorption-line clouds toward QSO1.