Galaxy Formation through Filamentary Accretion at z = 6.1

The formation of the first galaxies in the universe, those at $z\gt 6$, or within 1 Gyr of the Big Bang, is a question at the forefront of modern astronomy. How these galaxies accrete the gas that drives early star formation, and the relative dust attenuation, remains highly uncertain. These issues have attained new pertinence as constraints on the evolution of the neutral fraction of the intergalactic medium (IGM) have solidified. Current observational constraints based on, e.g., the cosmic microwave background (CMB), the Gunn–Peterson effect, and the effect of the neutral IGM on Lyα emission from early galaxies, suggest that the universe was highly ionized at $z\lt 6$ and significantly neutral at $z\gt 7$, with a plausible neutral fraction of 50% at $z\sim 7.5$ (Fan et al. 2006b; Robertson et al. 2015; Greig & Mesinger 2017; Ota et al. 2017). This evolution in the cosmic neutral fraction suggests substantial star formation at redshifts $z\sim 6$ to 8, if the IGM is reionized by star-forming galaxies (Fan et al. 2006a).

The Atacama Large Millimeter/submillimeter Array (ALMA) has opened a new window on the early universe, with the sensitivity and resolution to perform detailed studies of the dust and cool gas in the earliest galaxies. Of particular interest are observations of the [C ii] 158 μm line at $z\gt 6$. This line is typically the brightest emission line from star-forming galaxies, at meter through FIR wavelengths, tracing both the cool and warm interstellar medium (ISM) in galaxies (e.g., Carilli & Walter 2013; Pineda et al. 2013). Observation of dust in the most distant galaxies provides key constraints on the extinction properties, while the [C ii] line provides a means of tracing gas kinematics and obtaining information on the structure of the ISM. These quantities relate to the star formation rate of the galaxy.

Cosmological simulations of early galaxy formation have shown that the [C ii] line can be used not only to study the ISM of these sources (Vallini et al. 2013, 2015; Pallottini et al. 2017), but it might also probe the overall galaxy dynamics, including supernova-driven outflows (Gallerani et al. 2016). In the local universe, a well-established [C ii]-star formation rate (SFR) relation holds for a wide range of galaxy types, from metal-poor dwarf galaxies, to starbursts, ultra-luminous infrared galaxies, and active galactic nucleus hosting galaxies (De Looze et al. 2014; Pineda et al. 2014; Herrera-Camus et al. 2015). However, many of the high-z Lyman-Break Galaxies (LBGs) and Lyα Emitting Galaxies (LAEs) deviate from this relation for yet unclear reasons, likely related to their different metallicity, ionization state, or more efficient supernova feedback (Capak et al. 2015).

The galaxy WMH5 at z = 6.0695 is a useful test-bed to investigate detailed processes in the formation of early galaxies. This galaxy was discovered in the LBG search of Willott et al. (2013) as one of the most luminous LBGs at $z\gt 6$ (near-IR magnitude AB ∼ 24), with detected Lyα emission. Willott et al. (2015) discovered both thermal emission from dust, and [C ii] line emission, from this galaxy in Cycle 2 ALMA observations at 05 resolution. They presented a detailed analysis of the integrated spectral energy distribution (SED) of the galaxy from the optical through submm. The inferred SFR values are between 33 and 66 M_⊙ yr⁻¹, depending on the tracer, and the stellar mass is $2.3\times {10}^{10}\,{M}_{\odot }$. The overall SED has a dust attenuation of $E(B-V)=0.05$. However, the main component of the galaxy seen in the dust and [C ii] is spatially offset 04 west of the bright NIR component. The [C ii] at the NIR position is also significantly offset in velocity relative to the main [C ii] component (∼200 km s⁻¹). Willott et al. (2015) suggest a merging galaxy system; however, the spatial resolution and sensitivity of these early ALMA observations were insufficient to determine the detailed structure of this complex system.

In this paper, we present subarcsecond angular resolution and high sensitivity observations of the [C ii] and dust emission from WMH5 using ALMA in Cycle 3. These observations present a clear picture of the processes occurring in WMH5. We assume (${{\rm{\Omega }}}_{{\rm{\Lambda }}}$, ${{\rm{\Omega }}}_{m}$, h) = (0.692, 0.308, 0.678) (Planck Collaboration et al. 2016) throughout. At this distance, 1 arcsecond corresponds to 5.81 proper kpc at the redshift of WMH5 ($z=6.0695;$ Willott et al. 2015).

The [C ii] 158 μm and thermal dust emission from WMH5 were observed with ALMA in Cycle 3 in configuration C40-3 between 2016 August 3 and 7. Between 40 and 43 antennas were used for each observation. The total time was ∼6.5 hours, out of which 4 hours were spent on source. A total bandwidth of 7.5 GHz was used, split into four ∼2 GHz sub-bands spanning 253–257 GHz and 268–272 GHz (ALMA band 6). One sub-band was centered on the [C ii] emission, redshifted to 268.836 GHz, while the other three were used to measure the rest frame of ∼160 μm continuum emission. The Cycle 3 data were calibrated using the standard calibration script prepared by the Joint ALMA Observatory staff.

Using only the Cycle 3 data results in detections of both line and continuum emission. However, to improve the quality of the data, we combined these data with the Cycle 2 observations of Willott et al. (2015). In short, they observed in 2014 June with 29–32 antennas, 7.5 GHz of bandwidth, and had 95 minutes of on-source integration time.

An issue arose while combining the two data sets, based on the relative weights of the uv data. The definition of the weight factor as a function of channel width, integration time, and system temperature changed on 2014 September 4, when CASA updated from version 4.2.1 to 4.2.2.⁹ Because the Cycle 2 observations were made before this time and the Cycle 3 observations after, the average weights differed by multiple orders of magnitude (i.e., $w\sim 0.025$ for Cycle 2 and $w\sim 20$ for Cycle 3). When using the CASA task clean to combine these two data sets, this caused the Cycle 3 data to completely dominate.

Two methods were tested to correct this. The first method set all weights equal to one, so that all data were given equal weight. The other used the CASA task statwt, which resets weights based on the intrinsic scatter of the data. We compared images created using these approaches with the image made using the original weights. Both re-weighting methods resulted in images with lower rms noise levels than the original image and similar channel maps. Because the statwt image had the lowest rms noise level and had weights based on properties of the data, we will use this weighting method for the following analysis.

While different visibility weighting schemes were applied to the data (i.e., Briggs weighting with robust parameters of 2.0, 1.5, and 1.0), the best results were found for natural (robust 2.0) weighting. The synthesized beam was 030 × 025, with major axis position angle = 74°. Using the three sub-bands that were not centered on the [C ii] line, we created a continuum image with an rms noise level of $9\,\mu \mathrm{Jy}$ beam⁻¹. The remaining sub-band was used to create a line emission cube with $15.6\,\mathrm{MHz}\sim 17$ km s⁻¹ channels and an rms noise level of 0.1 mJy beam⁻¹.

Figure 1 shows the [C ii] channel images. We see emission from ∼+260 to −270 km s⁻¹, mostly going from west to east with increasing frequency (decreasing velocity). We also note that, in addition to this higher level emission, there is substantial (i.e., $\geqslant 2\sigma$) emission present to the east in select channels from $\sim +250$ to +110 km s⁻¹. We use these images and those that follow to define two positions. The centers of two-dimensional Gaussian fits to a map using only channels where emission is concentrated in the west, named the "core" (2^h26^m270025 −4°52'384042), and of the map made using the channels where line emission extends to the east, named the "tail" (2^h26^m27.0208^s −4°52'384007). In this and subsequent analyses, we adopt zero velocity as [C ii] emission redshifted to z = 6.0695, consistent with Willott et al. (2015).

Zoom In Zoom Out Reset image size

The channel images in Figure 1 show that most of the emission is concentrated in the core, but there are two distinct knots of emission in the tail to the east at ∼−110 and ∼−250 km s⁻¹. These knots are more elongated east–west, extending E–W by about 06 in individual channels.

Figure 2 shows the integrated [C ii] spectrum of the system and the spectra taken at the locations of the core and tail, over a single spatial pixel in the cleaned image. While the core is well fit by a single Gaussian, the tail is split into three different contributions: two narrow peaks, which we call SG-1 (or Sub-Galaxy-1) and SG-2, and a spectrally broad component, which we name "Diffuse." We performed Gaussian fitting to the apparent three components of the tail. Figure 3 shows the resulting fits for the tail. Details of the Gaussian fits for the core and tail are given in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  One-dimensional Gaussian Fits to [C ii] Components

Source	Line Peak	Velocity	FWHM	$S{\rm{\Delta }}{v}_{[{\rm{C}}{\rm{II}}]}$
	[mJy beam−1]	[km s−1]	[km s−1]	[Jy beam−1 km s−1]
Core	0.76 ± 0.03	−1 ± 5	270 ± 13	0.22 ± 0.01
SG-1	0.50 ± 0.05	−249 ± 5	94 ± 12	0.05 ± 0.01
SG-2	0.4 ± 0.1	−111 ± 5	72 ± 18	0.03 ± 0.02
Diffuse	0.24 ± 0.03	28 ± 34	238 ± 74	0.06 ± 0.03

Download table as:  ASCIITypeset image

The spectrally broad "Diffuse" component exhibits a similar centroid velocity as the core emission, suggesting that the diffuse emission is the spatial edge of the core emission. Alternatively, the Diffuse emission may originate from numerous discrete, spectrally narrow components that we cannot separate.

Figure 4 shows the velocity integrated (moment zero) [C ii] emission contours overlaid on the ground-based near-IR image. A 2D Gaussian fit of the [C ii] emission returns a deconvolved size of $(0.65\pm 0.08)^{\prime\prime} \times (0.30\pm 0.04)^{\prime\prime}$ at $(81\pm 6)^\circ$, an integrated flux density of 0.70 ± 0.08 Jy km s⁻¹, and a peak flux density of 0.19 ± 0.02 Jy beam⁻¹ km s⁻¹. However, the [C ii] emission is non-Gaussian in shape due to the multiple components already discussed. The center of the near-IR emission is located close to the tail position, not at the peak of the core, as previously discussed by Willott et al. (2015). The spatial resolution of the near-IR image is insufficient to determine the fraction of flux coming from each region.

Zoom In Zoom Out Reset image size

Figure 5 shows the $\sim 160\,\mu {\rm{m}}$ continuum emission. A Gaussian fit yields a deconvolved size of $(0.5\pm 0.1)^{\prime\prime} \times (0.3\pm 0.1)^{\prime\prime}$ at $96^\circ \pm 67^\circ$, a peak flux density of $53\pm 9\,\mu \mathrm{Jy}$ beam⁻¹, an integrated flux density of $0.15\pm 0.03\,$ mJy, and a central position of 2^h26^m270099 −4°52'383770. The morphology is obviously different from the beam shape, with an extension to the east and a $3\sigma$ extension to the north. The peak of the continuum is located between the [C ii] core and tail positions, but because this emission is a blend of the core and tail components, this is to be expected.

Zoom In Zoom Out Reset image size

4.1. [C ii] Luminosities and SFRs

We first analyze the line emission, integrated over the velocity ranges corresponding predominantly to the core and sub-galaxies. The left side of Figure 6 shows the [C ii] emission integrated over velocity ranges corresponding to the core galaxy (169 to −77 km s⁻¹), while the right-hand figure shows the [C ii] emission integrated over the velocity range of the sub-galaxies (−94 to −305 km s⁻¹). We then fit Gaussian models to the surface brightness distributions. We note that the spatial/spectral separation of the core and sub-galaxies is not completely clean, so the flux densities of the sub-galaxy complex may be overestimated.

Zoom In Zoom Out Reset image size

The core moment zero image shows a higher peak, but is spatially more compact, with a deconvolved size of $(0.32\,\pm 0.03)^{\prime\prime} \times (0.14\pm 0.01)^{\prime\prime}$, and a velocity integrated flux density of 0.32 ± 0.03 Jy km s⁻¹. The moment zero image of the sub-galaxies shows a larger spatial size, extending toward the core galaxy. The fit size is $(0.6\pm 0.1)^{\prime\prime} \times (0.38\pm 0.07)^{\prime\prime}$ and the velocity integrated flux density is 0.30 ± 0.04 Jy km s⁻¹. These [C ii] line fluxes imply [C ii] luminosities of: ${L}_{[{\rm{C}}{\rm{II}}],\mathrm{CORE}}=(3.2\pm 0.3)\times {10}^{8}\,{L}_{\odot }$ and ${L}_{[{\rm{C}}{\rm{II}}],\mathrm{SG}}=(3.0\,\pm 0.4)\times {10}^{8}\,{L}_{\odot }$.

While the combined flux densities of the tail components in Table 1 are considerably less than those of the core, these luminosities are comparable. This is reasonable, as these luminosities are derived from flux densities integrated over the entire emission area of each source, while the values in Table 1 represent only a single pixel. Because the tail emission is considerably more extended, its flux density at a given point will be less than that of the concentrated core, but the sum of each component may be comparable.

Using these luminosities, we can estimate the star formation rates, assuming the relationship given in Equation (8) of Vallini et al. (2015) and a metallicity. Note that this relationship is based on the results of a cosmological simulation, where discrete values of metallicity (i.e., 1.0, 0.2, and $0.05\,{Z}_{\odot }$) were assumed. This coarse sampling creates an inherent uncertainty in the relationship, which is exacerbated by additional assumptions (e.g., a Kennicutt–Schmidt power law index). Assuming a solar metallicity, the implied star formation rates for the core and sub-galaxies are $\sim 20\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. If we assume a metallicity of $0.1\,{Z}_{\odot }$, these values increase to $\sim 80\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$.

Equation (8) of Vallini et al. may also be used to solve for the metallicity of an object, given its SFR and [C ii] luminosity. To find the total SFR of WMH5, we add SFR_SED and SFR_FIR of Willott et al. (2015), yielding $76\pm 11\,{{\rm{M}}}_{\odot }\,{\mathrm{yr}}^{-1}$. This addition is due to the fact that the core galaxy (their galaxy "A") dominates the observed mm (rest-frame FIR) continuum emission, while the tail galaxies (their galaxy "B") is the source of the observed NIR (rest-frame UV) emission. Because the sources of the mm and NIR emission are spatially separated, and the constructed SED is mainly shaped by the UV data, the addition of the two SFRs is required to find the total value for the galaxy complex.

We use the total SFR and our total [C ii] luminosity of ($7.0\pm 0.8)\times {10}^{8}\,{L}_{\odot }$, yielding ∼0.2 Z_⊙ for the system as a whole. However, due to the discrete metallicities probed by Vallini et al. (see their Figure 8), we only state a metallicity limit of $0.1\mbox{--}1.0\,{Z}_{\odot }$.

4.2. [C ii]/Continuum Ratio

At this high resolution, we are able to separately analyze the continuum and line surface brightness at the position of the core and sub-galaxies. The continuum surface brightnesses at the core and tail positions are $(49\pm 9)\,\mu \mathrm{Jy}$ beam⁻¹ and $(42\pm 9)\,\mu \mathrm{Jy}$ beam⁻¹, respectively.

For the line surface brightness, we use the spectral fitting results shown in Table 1. The core [C ii] velocity integrated surface brightness is 0.21 ± 0.01 Jy beam⁻¹ km s⁻¹. For the tail position, the sum of the emission from the two sub-galaxies and the diffuse component gives: 0.15 ± 0.08 Jy beam⁻¹ km s⁻¹. The line to continuum ratio for the core is then: $(4.3\pm 0.8)\times {10}^{3}$ km s⁻¹. The value for the tail is $(4\pm 2)\,\times {10}^{3}$ km s⁻¹.

We compare these ratios to observations of other high-redshift galaxies. In the sample of Capak et al. (2015), which consisted of $z\sim 5.5$ main sequence galaxies observed in [C ii], four galaxies were detected in $\sim 158\,\mu {\rm{m}}$ (rest frame) continuum emission. These four give ratios of $(6\pm 2)\,\times {10}^{3}$ km s⁻¹ (HZ4),¹⁰ $(1.4\pm 0.4)\times {10}^{4}$ km s⁻¹ (HZ6), $(3.8\,\pm 0.3)\times {10}^{3}$ km s⁻¹ (HZ9), and $(1.25\pm 0.08)\times {10}^{3}$ km s⁻¹ (HZ10). Thus, the ratios of both segments of WMH5 fall within the scatter of the Capak et al. sample. Willott et al. (2015) show that the system as a whole has a [C ii] line luminosity to FIR continuum luminosity ratio $\sim 5\times {10}^{-3}$, which is larger than low redshift star-forming galaxies by a factor of a few. They point out that such a large ratio: "...suggests extended star formation with low metallicity and an intense radiation field."

In addition to these two discrete locations, we may examine the ratio of [C ii] to FIR luminosities over the source as a whole (see Figure 7). The integrated [C ii] flux density was converted to a luminosity using Equation (1) of Solomon et al. (1992). A graybody spectrum identical to that used by Willott et al. (2015) was adopted to convert our continuum flux density to a FIR luminosity (i.e., ${T}_{\mathrm{dust}}=30$ K, $\beta =1.6$, and a FIR range of λ = 42.5–122.5 μm). Note that we assume that these parameters are uniform across the source.

Zoom In Zoom Out Reset image size

This ratio plot shows that while the core and far-east tail areas behave similarly, there is a gap of low ${L}_{[{\rm{C}}{\rm{II}}]}/{L}_{\mathrm{FIR}}$ between them. A similar area of low [C ii] emission is observed in the simulations of Pallottini et al. (2017). In their Section 3.3.3, they reason that this "middle ground" contains diffuse (low n) molecular clouds, which only weakly emit [C ii] (${L}_{[{\rm{C}}{\rm{II}}]}\propto n$). In addition, the relatively high temperature of the CMB at $z\sim 6$ (∼19 K) diminishes the [C ii] luminosity.

4.3. Velocity Behavior

Figure 8 is a position–velocity (PV) plot. Again, in addition to the core emission at zero offset, two distinct tail components are evident, plus diffuse emission throughout. The spatial offset of the two tail components from the core is 024 in projection, or about 1.4 kpc.

Zoom In Zoom Out Reset image size

This PV plot is similar to Figure 6 of Willott et al. (2015), although at higher spatial resolution. While they detected SG-1 (their component "B"), we are able to separate their component "A" into the core galaxy and SG-2. Both sub-galaxies are clearly distinct from the core galaxy and show similar positional offsets.

Our best-fit velocities (Table 1) and those of Willott et al. (2015, their Figure 1) are also in agreement, with component "A" (${v}_{[{\rm{C}}{\rm{II}}]}=0\pm 13$ km s⁻¹) and our "core" (${v}_{[{\rm{C}}{\rm{II}}]}\,=-1\pm 5$ km s⁻¹) in agreement. In addition, component "B" (${v}_{[{\rm{C}}{\rm{II}}]}=-238\pm 13$ km s⁻¹) and SG-1 (${v}_{[{\rm{C}}{\rm{II}}]}=-249\pm 5$ km s⁻¹) agree.

While the velocity of the tail components is multivalued, we may create a velocity weighted (moment 1) image of the core galaxy. The zero of velocity here is assumed to be the [C ii] emission redshifted to z = 6.0695. Figure 9 shows that the core galaxy is not highly reminiscent of either a solid-body rotator or a classic rotating disk but does show an identifiable velocity gradient.

Zoom In Zoom Out Reset image size

4.4. Dust Mass

The merger interpretation emerging from the present observations is supported by recent AMR, zoomed simulations of high-redshift galaxies with properties similar to WMH5. Pallottini et al. (2017) studied in detail Dahlia, a typical $z\approx 6$ LBG with an SFR $\approx 100\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and an effective stellar radius of ∼0.6 kpc. This model galaxy has 14 satellites distributed over a distance of 100 kpc from it. The six largest ones have a dark matter halo mass in the range ${M}_{h}=2.5-12\times {10}^{9}\,{M}_{\odot }$, three of which are located within 3 kpc from Dahlia's center (whereas our sub-galaxies are $\lt 6\,\mathrm{kpc}$ distant). The two most prominent among these satellites (at a distance of 12 and 15 kpc from the center) sustain considerable star formation rates of 30 and $10\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, and have a stellar mass of ${M}_{* }=6\times {10}^{9}\,{M}_{\odot }$ and ${M}_{* }=2\times {10}^{9}\,{M}_{\odot }$. Their star formation history is relatively bursty and increases with time. They are also the oldest among the satellites, with an age at z = 6 of 350 Myr. Compared to Dahlia (${M}_{* }=1.6\times {10}^{10}$), their stellar mass is 38% and 12% smaller, respectively. The satellites are mostly aligned along three main filaments connecting to the center (see Figure 4 of Pallottini et al. 2017); they are also much more compact than Dahlia: the largest satellite within 3 kpc has an effective radius of 0.2 kpc, i.e., about one-third that of the central galaxy.

Interestingly, the three nearby satellites are clearly visible in the simulated [C ii] maps, thus supporting the analogy with our observations. The satellites move at ∼80 km s⁻¹, i.e., about half the virial velocity of Dahlia. They also have very low metallicity, due to their relatively young age (≤100 Myr). Compared to Dahlia, which has a metallicity equal to $0.5\,{Z}_{\odot }$, the nearby satellites barely reach $Z={10}^{-2}\,{Z}_{\odot }$; the two distant ones instead have $Z=0.03\,{Z}_{\odot }$ and $Z=0.1\,{Z}_{\odot }$. If the dust content scales with metallicity, this should result in a substantially lower extinction. Thus, the comparison with theoretical models supports the idea that WMH5 features the build-up of a galaxy occurring via infall and merging of a number of star-forming satellite sub-galaxies.

Of course, we have to be aware that the comparison with observations is performed here only against a single, although very similar, simulated system. Thus, it is probably unsafe to push the comparison beyond the present analysis. Before drawing more robust conclusions, one would need to extend the study to a larger sample of simulated galaxies that could inform us on the frequency of multiple/merger systems expected at these early redshifts. Nevertheless, the emerging scenario is one in which an intriguing resemblance between theoretical expectations and the data analysis certainly exists.

Our geometry is also akin to that found in Katz et al. (2016, their Figure 20), in that we have a "filament" of gas falling into a galaxy, with spatially separated [C ii] and UV/Lyα emission. Katz et al. interpret this geometry as being caused by the infalling galaxy causing new star formation at the interface of the filament and galaxy, while introducing relatively pristine (low metallicity) gas. This low-Z environment is weak in [C ii] emission, but high in UV/Lyα. Deeper inside the host galaxy, the environment is more processed (higher metallicity), so while the UV/Lyα emission is obscured, [C ii] emission may escape.

New ALMA observations of the dust continuum and [C ii] 158 μm emission from the z = 6.0695 LBG WMH5 show a compact main galaxy in dust and [C ii] emission, with two, what we call, "sub-galaxies" located along a filamentary tail extending about 5 kpc (in projection) to the east. These sub-galaxies are distinct in velocity, with narrow velocity dispersions of about 80 km s⁻¹, and the tail joins smoothly into the main galaxy velocity field. The sub-galaxies themselves are extended east–west by about 06 in individual channel images.

We conjecture that the WMH5 system represents the early formation of a galaxy through the accretion of smaller satellite galaxies embedded in a smoother gas distribution, along a filamentary structure. While we cannot rule out outflow, we find that notion less likely, given the narrow velocity dispersion of the individual sub-galaxies. Such narrow velocity dispersion is more characteristic of gravitational or tidal features than strong outflows. Likewise, the significant dust masses and star formation rates suggest individual sub-galaxies, as might be expected in a coalescing galaxy system.

Our filamentary interpretation assumes a three-dimensional structure that we are unable to directly observe. However, the continuity of the velocity field from the core to each sub-galaxy, the small velocity widths of each sub-galaxy, the intrinsic E–W extension in channel images of the sub-galaxies themselves, and the similarity of these results to those found using three-dimensional cosmological simulations all provide circumstantial evidence supporting this conclusion, although the term filamentary remains presumptive.

These observations highlight the power of ALMA to image early galaxy formation. The results support the idea that early galaxies can be enriched in heavy elements and dust due to their very early history of star formation, including enrichment of the material that accretes onto main galaxy. The new Band 5 system of ALMA is currently being commissioned (163–211 GHz). This band will open the z = 8–11 range for high resolution, sensitive observations of the [C ii] 158 μm line.

This paper makes use of the following ALMA data: ADS/JAO.ALMA #2013.1.00815.S and #2015.1.00834.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in co-operation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. G.C.J. is grateful for support from NRAO through the Grote Reber Doctoral Fellowship Program. A.F. acknowledges support from the ERC Advanced Grant INTERSTELLAR H2020/740120. R.W. acknowledges supports from the Thousand Youth Talents Program of China, the National Science Foundation of China (NSFC) grants No. 11473004, 11533001, and the National Key Program for Science and Technology Research and Development (grant 2016YFA0400703). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.