ALMA OBSERVATIONS OF Lyα BLOB 1: HALO SUBSTRUCTURE ILLUMINATED FROM WITHIN

SSA22-LAB01 was observed in two projects (2013.1.00922.S [PI Geach] and 2013.1.00704S [PI Matsuda]) with a similar configuration in Band 7: the full 7.5 GHz of bandwidth centered at 347.59 GHz was used to measure the continuum emission at approximately the same frequency as the SCUBA-2 detection of the target. A total of 36–40 12 m antennas were used in the observations, with baselines spanning 21–918 m, which were capable of recovering emission on scales of up to 5'', with a maximum resolution of 04. Antenna T_sys temperatures were approximately 100–150 K and the mean precipitable water vapour column was 0.249 mm. Observations of the phase calibrator source J2206-0031 confirm consistent flux scaling between the two projects in the calibration phase. All calibrations were performed using the casa software and the visibilities from each project concatenated into a single measurement set for which the total integration time is 2860 s.

The dirty image reveals the two main components of submillimeter emission, and we use this information to supply circular masks (each 2'' in radius) at 22^h17^m260, +00°12'363 and 22^h17^m261, +00°12'324 (J2000) to the casa clean task. We use a natural weighting of the visibilities and clean down to a threshold of 70 μJy within 1000 iterations. We adopt multi-scale cleaning, allowing model Gaussian components of width 0'' (delta function), 05, 1'', 2'', and 3''. The clean is run twice: at full resolution and with an outer taper of 1''. To both maps, we "feather" the single dish SCUBA-2 (regridded) map to improve the recovery of the total flux of the source and potentially improve sensitivity to extended emission. The depth of the map is 40 μJy beam⁻¹, increasing to 90 μJy beam⁻¹ after tapering. Final maps are corrected for primary beam attenuation for analysis. The synthesized beams are 040 × 038 (PA 23°) and 101 × 098 (PA 112°) in the feathered full resolution and 1'' tapered maps, respectively.

To measure the integrated flux densities and sizes, we threshold the maps at 3σ, where σ is the root-mean-square noise measured around the sources and sum the flux enclosed by each 3σ contour. In the full resolution map, there are three distinct components: a, b, and c. Integrated flux density uncertainties are estimated from the root-mean-squared value, scaled to the number of beams subtended by the 3σ contour. The sizes of each source in the tapered map are measured as the maximum width of the 3σ contour, with an uncertainty $\delta \theta =0.6\times {\theta }_{\mathrm{beam}}/{\rm{S}}/{\rm{N}}$, where θ_beam is the full width half maximum of the beam and the S/N = 3.

Note that there is a discrepancy between the single dish flux and the total flux measured in the ALMA map: ${\rm{\Delta }}{S}_{850}=2.9\pm 1.1$ mJy. Although this is not significant at the 3σ level, we should discuss this potential "missing" flux. First, it is important to note that single dish flux densities in the submillimeter can be statistically boosted; Geach et al. (2016) measured the flux boosting as a function of the S/N in the SCUBA-2 Cosmology Legacy Survey (where the SCUBA-2 detection of this target originates), finding an average boosting of approximately 30% for SSA22-LAB01. The array configuration allows us to recover emission on scales of up to 5''. It is unlikely that there is an 850 μm emission component on scales larger than this (40 pkpc at this redshift), but given the sensitivity of our observations, it implies that, if real, the missing submillimeter emission is extended on scales of at least 3'' (25 pkpc). One possibility is that the emission is spread over a large number of faint clumps around the central sources—an idea we explore later in the paper. Deeper observations with a slightly more compact array configuration would be beneficial to address this.

Integral field spectroscopy was performed with the Very Large Telescope MUSE instrument. The observations were carried out under clear or photometric conditions, between 2014 November and 2015 September. We used the extended blue setting, resulting in a minimum wavelength of 4650 Å. Each integration was 1500 s, and the field was rotated by 90° between each exposure in order to reduce fixed pattern noise from the integral field units and residual flat-field errors. Data were reduced with the MUSE pipeline (version 1.0.5), following standard procedures for bias subtraction, dark current removal, flat-fielding, and basic calibration of the individual integrations. Each individual spectrum was fully post-processed to output an image of the field of view, data cubes, and pixel tables, which we examined individually. We computed shifts between the individual exposures using standard centering tools in iraf. Finally, we used the muse_exp_combine task to drizzle and stack all the individually reduced pixtables.

The HST Imaging Spectrograph (STIS) observations (project 9174, Chapman et al. 2004) used the 50CCD clear filter, which provides a wide response in the optical, with an effective central wavelength of 5850 Å and a full width half maximum throughput of 4410 Å over a 52'' × 52'' field of view. Six exposures of SSA22-LAB01 were taken in "LOW-SKY" time for a total of 7020 s of integration. HST Legacy Archive calibrated images were combined into a single deep co-add with 01 pixels using the DrizzlePac software (version 2) from the Space Telescope Science Institute. The 5σ sensitivity limit is 27.6 mag (Chapman et al. 2004).

K-band (1.92–2.30 μm) spectroscopy of SSA22-LAB01 was obtained on a single slit of a multislit mask observed during science verification of the MOSFIRE (McLean et al. 2012) on the W. M. Keck Observatory Keck I 10 m telescope. The slit width was 07, resulting in a spectral resolving power of R ∼ 3600. The total integration time was 5040 s, in a sequence of 28 180 s exposures using an ABAB nod pattern with a nod amplitude of 30. The seeing during the observation was estimated to be 035 full width half maximum. The data were reduced by the MOSFIRE Data Reduction Pipeline (Steidel et al. 2014).

At full (04) resolution, we resolve three main components that we label "a," "b," and "c" within the SSA22-LAB01 (Figure 1). Component "a" is the brightest submillimeter source, located close to the geometric center of the LAB and is 1'' from component "b"; in the tapered map these effectively blend into a single source, and we consider a+b jointly hereafter. Component "c" is located 45 to the southeast and is associated with strong stellar continuum emission, with a precise redshift measured from [O iii] and Hβ lines, z = 3.1000 ± 0.0003, placing it within the Lyα nebula (Kubo et al. 2015). Component a+b is also coincident with stellar continuum emission of complex, clumpy morphology—it is not clear whether this is a single bound system and we are simply detecting bright clumps of submillimeter and optical emission, or if we are witnessing the coalescence of several independent low-mass systems.

Zoom In Zoom Out Reset image size

The total flux density of the signal above the 3σ level in the tapered map is S₈₅₀ = 1.68 ± 0.06 mJy; Table 1 gives the positions and flux density of the two main components a+b and c. Assuming that cool dust emission traces the dense ISM (Scoville et al. 2015), the corresponding total molecular gas mass is ${M}_{\mathrm{gas}}\approx 4\times {10}^{10}$ M_⊙ and the average molecular gas surface density is Σ_gas ≈ 140 ${M}_{\odot }$ pc⁻². If Σ_SFR scales linearly with Σ_gas (Wong & Blitz 2002), we estimate that the submillimeter sources are forming stars at rates of the order of 100 ${M}_{\odot }$ yr⁻¹. We have poor constraints on the shape of the spectral energy distribution (SED) of these sources, but we can estimate their integrated total infrared luminosities by assuming an appropriate range of templates describing the thermal dust emission of active galaxies. We adopt the Dale & Helou (2002) set of templates, characterized by a parameter α, describing the power-law index for the distribution of dust mass over heating by the interstellar radation field (U): ${dM}(U)\propto {U}^{-\alpha }{dU}$, with most galaxies falling within the range of α = 1–2.5; we conservatively consider this full range when estimating the total infrared luminosities of the sources. Redshifting to z = 3.1 and normalizing to the observed 850 μm flux density, we estimate total infrared (8–1000 μm) luminosities in the range of ${L}_{\mathrm{IR}}\approx (0.2\mbox{--}1.9)\times {10}^{12}\,{L}_{\odot }$ and ${L}_{\mathrm{IR}}\approx (0.2\mbox{--}1.5)\times {10}^{12}\,{L}_{\odot }$ for sources a+b and c respectively. The median infrared luminosity considering the full range of templates is ${L}_{\mathrm{IR}}\approx 6.0\pm 0.3\times {10}^{11}\,{L}_{\odot }$ and ${L}_{\mathrm{IR}}\approx 4.6\pm 0.3\times {10}^{11}\,{L}_{\odot }$ for the two components, respectively, where the error bar reflects the measurement uncertainty on the integrated ALMA flux. Assuming a star-formation rate calibration following Kennicutt & Evans (2012) these average luminosities correspond to SFR ≈ 90 M_⊙ yr⁻¹ and $\mathrm{SFR}\approx 70\,{M}_{\odot }$ yr⁻¹, respectively, consistent with the ISM mass scaling value derived above, though the SFRs could be a factor of three higher if the galaxies are better represented by the more "active" templates of the Dale & Helou library.

The emergent luminosity of Lyα emission associated with obscured star formation is ${L}_{\mathrm{Ly}\alpha }=0.05{f}_{\mathrm{esc},\mathrm{Ly}\alpha }{L}_{\mathrm{IR}}$ (Dijkstra & Westra 2010), under the assumption of Case B conditions and where ${f}_{\mathrm{esc},\mathrm{Ly}\alpha }$ is the escape fraction (of Lyα photons) from the galaxy. This corresponds to ${L}_{\mathrm{Ly}\alpha }\sim {f}_{\mathrm{esc},\mathrm{Ly}\alpha }{10}^{44}$ erg s⁻¹ per ALMA source, similar to the total Lyα luminosity of SSA22-LAB01. From the models of Leitherer et al. (1999), the ionizing flux (ν = 200–912 Å) scales like ${L}_{200-912\mathring{\rm{A}} }\approx \mathrm{SFR}\times {f}_{\mathrm{esc},\mathrm{UV}}\times {10}^{43}$ erg s⁻¹ (for 100 Myr of continuous star formation, solar metallicity and Salpeter IMF). Thus, for star-formation rates of the order  of 100 M_⊙ yr⁻¹ and modest high escape fractions, the ALMA sources appear to be viable "power sources" for extended Lyα emission through photoionization. The next question follows naturally: what is the nature of the environment around these central star-forming sources?

Like other giant LABs (e.g., Prescott et al. 2012), SSA22-LAB01 is thought to trace a massive dark matter halo of the order of ${M}_{{\rm{h}}}\approx {10}^{13}\,{M}_{\odot };$ this is close to the mass where cold flows are not expected to survive inside the virial radius at this epoch (Dekel et al. 2009). Recently, moving mesh hydrodynamic simulations suggest cold streams are assimilated into the hot halo gas inside the virial radius even down to halo masses of ${M}_{\mathrm{halo}}\approx {10}^{11.5}\,{M}_{\odot }$ (Nelson et al. 2013), suggesting that cold flows are not sufficient to explain the extended Lyα luminosity of giant LABs. Following arguments put forward by Hayes et al. (2011), Geach et al. (2014), and Beck et al. (2016), we postulate that extended line emission in SSA22-LAB01 could arise from Lyα photons generated by central star formation in the submillimeter sources, escaping and resonantly scattering in neutral gas associated with a large ensemble of lower-mass galaxies in the same potential well (e.g., Cen & Zheng 2013, Francis et al. 2013). This is consistent with recent spectropolarimetry observations of SSA22-LAB01 that favor a scenario in which Lyα photons undergo long flights from central sources before scattering in H i in the CGM (Hayes et al. 2011; Beck et al. 2016, though see Trebitsch et al. 2016). In Figure 1, we show the polarization vectors of Hayes et al. (2011) that describe a circular pattern roughly centered on the submillimeter sources.

The STIS imaging reveals many faint (f_λ ≈ 10⁻¹⁹ erg s⁻¹ cm⁻² Å⁻¹) clumps of UV emission surrounding the submillimeter sources (Figure 1, Chapman et al. 2004; Matsuda et al. 2007, Uchimoto et al. 2008, see also Prescott et al. 2012). There is no spectroscopically confirmed redshift for source a+b, and, unlike source c, it is not unambiguously associated with an STIS counterpart. Instead, there are a number of STIS clumps at (and around) the position of a+b, and the optical colors of these are broadly consistent with a z = 3.1 redshift (Kubo et al. 2015). In addition, there is a possible H i absorption feature in the Lyα spectrum at this location (Weijmans et al. 2010; Geach et al. 2014); thus, source a+b is unlikely to be a chance projection of a submillimeter source at a different redshift. In the event that this ALMA source is not associated with the LAB, then obviously the main impact on our analysis and interpretation is that there will be approximately a factor of two fewer Lyα and UV photons generated by central star formation.

If some (a fraction are likely to be chance projections) of the other faint STIS sources represent low-mass satellites of the ALMA sources, then the cold gas associated with them could be responsible for scattering Lyα photons emerging from the central active galaxies. A key test is actually confirming the "membership" of the STIS sources, and to this end we have obtained near-infrared spectroscopic confirmation of one of these faint companion sources "S1" (Figure 2), detecting the [O iii] line at z = 3.0968. S1 is 21 (17 pkpc) from ALMA source "c" with a line-of-sight velocity offset of just 250 km s⁻¹. Although not a complete survey by any means, this observation does support the view that the central ALMA sources are accompanied by a population of low-mass satellites. With this empirical evidence in hand, we now consider whether the same observations are predicted by galaxy formation and evolution models. In particular, we run a zoomed cosmological hydrodynamical simulation coupled with radiative transfer schemes to produce mock observations that can be directly compared to the data. The main objective is to test whether systems similar to SSA22-LAB01 arise in current galaxy formation schemes.

Zoom In Zoom Out Reset image size

We simulate the formation of a galaxy in a massive (${M}_{{\rm{h}}}\approx {10}^{13}\,{M}_{\odot }$) halo (by z = 2) utilizing a cosmological zoom technique with the hydrodynamics code gizmo (Hopkins 2014). We employ a pressure-entropy formulation of smoothed particle hydrodynamics, which conserves momentum, energy, angular momentum, and entropy. We first run a coarse large (144 Mpc³) dark matter-only simulation to identify a halo of interest. This simulation is run at a relatively low-mass resolution ${M}_{\mathrm{dm}}=8.4\times {10}^{8}\,{M}_{\odot }$. At z = 2, we select a massive ($2.03\times {10}^{13}\,{M}_{\odot }$) halo, and re-simulate the evolution of this halo at significantly higher resolution (Feldmann et al. 2016), including baryons with a particle mass of m_bar = 2.7 × 10⁵ M_⊙. The initial conditions are generated with the music code (Hahn & Abel 2011) and the minimum baryonic, star, and dark matter force softening lengths are, respectively, 9, 21, and 142 proper parsecs at z = 2. The halo simulated here is the exact same as the main halo in Narayanan et al. (2016).

Gas cools using a cooling curve that includes both atomic and molecular line emission, with the neutral ISM broken into atomic and molecular components following an analytic algorithm describing the molecular fraction in a cloud based on the gas surface density and metallicity (Krumholz et al. 2009). Star formation occurs in molecular gas and follows a volumetric relation $\dot{{\rho }_{\star }}={\rho }_{\mathrm{mol}}/{t}_{\mathrm{ff}}$, where ρ_mol is the volume density of molecular gas and t_ff is the freefall time. Star formation is restricted to gas that is locally self-gravitating following the algorithms of Hopkins et al. (2013). Once stars have formed, they impact the ISM via a number of feedback channels. In particular, we include models for momentum deposition by stellar radiation and supernovae, photoheating of H ii regions, and stellar winds (see Hopkins 2014 for details). gizmo tracks 11 metal species (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) with yields from SN 1a and SN II following Iwamoto et al. (1999) and Woosley & Weaver (1995).

To model the ALMA and STIS observations we employ the dust radiative transfer code Powderday which computes the stellar SEDs of the star particles (constrained by their ages and metallicities) and propagates this radiation through the dusty ISM. The gas particles are projected onto an adaptive grid with an octree-like memory structure (using yt, Turk et al. 2011), and the dust mass is assumed to be a constant fraction (40%) of the total metal mass in a given cell (Dwek 1998; Watson 2011). Photons are emitted in a Monte Carlo fashion and absorbed, re-radiated, and scattered by dust; this process is iterated upon until the dust temperature has converged. The stellar SEDs are calculated at the run time for each stellar particle utilizing fsps (Conroy & Gunn 2010). We employ the Padova iscochrones and assume a Kroupa IMF. The dust radiative transfer employs hyperion (Robitaille 2011). We assume a Weingartner & Draine (2001) dust size distribution, with R_V = 3.15. To create mock observations, we consider the observed frame 850 μm and optical emission predicted by Powderday and hyperion projected onto a celestial grid at the position of SSA22-LAB01.

To simulate the ALMA observations, we use the casa task simalma, setting the simulation parameters to match the real observations (Section 2.1), including the same array configuration. The resulting maps have a nearly identical synthesized beam and 1σ noise as the data. For the STIS imaging, we convolve the predicted observed frame UV-through-optical emission with the 50CCD throughput curve to generate single images projected along the same view angles as the 850 μm images. These images are then convolved with the STIS point-spread function (PSF) and gridded to an identical pixel scale as the observations. Finally, we create a noise image by removing >1.5σ features from the real STIS image and replace those pixels with values randomly drawn from the remainder of the pixels in the image. The noise-free, PSF-convolved flux image is then added to the noise image to create the mock observation.

Finally, we consider Lyα radiative transfer of centrally generated photons through the simulation volume. We use the code tlac (Gronke et al. 2014, 2015) with the default core skipping parameter x = 3. The central 400 pkpc of simulation volume is interpolated onto a regular grid, and we treat two central, star-forming galaxies (SFRs 189 M_⊙ yr⁻¹ and 115 M_⊙ yr⁻¹) as sources of Lyα photons. To compute the Lyα luminosity emitted by the cells, we use the relation between star-formation rate and Lyα luminosity (applicable to solar metallicity and a Kroupa IMF)

$$\begin{eqnarray}&&{L}_{\mathrm{Ly}\alpha }/\mathrm{erg}\,{{\rm{s}}}^{-1}=1.87\times {10}^{42}\times (\mathrm{SFR}/{M}_{\odot }\,{\mathrm{yr}}^{-1}).\end{eqnarray} \tag{ 1 }$$

The photons are emitted from a random point inside the corresponding cell with an emission frequency drawn from a distribution depending on the thermal velocity of the gas in the cell (Gronke et al. 2014, 2015) using the "peeling" algorithm (Yusef-Zadeh et al. 1984; Zheng & Miralda-Escudé 2002; Laursen 2010). We follow the prescriptions of Laursen (2010) for the Small Magellanic Cloud (Pei 1992; Weingartner & Draine 2001; Gnedin et al. 2008) when calculating the dust content in the grid, and a full description of this method can be found in that work. We compute an effective dust density n_d at every cell as

$$\begin{eqnarray}&&{n}_{{\rm{d}}}=({n}_{{\rm{H}}{\rm{I}}}+{f}_{\mathrm{ion}}{n}_{{\rm{H}}{\rm{II}}})(Z/\langle {Z}_{\mathrm{SMC}}\rangle ),\end{eqnarray} \tag{ 2 }$$

where ${n}_{{\rm{H}}{\rm{I}}}$ and ${n}_{{\rm{H}}{\rm{II}}}$ are the neutral and ionized hydrogen densities, respectively, at every cell. Z and $\langle {Z}_{\mathrm{SMC}}\rangle$ denote the metallicity at every cell and the SMC average, respectively, and f_ion accounts for the dust content in the ionized gas (Laursen 2010). We assume that f_ion is 1% and note that $\langle {Z}_{\mathrm{SMC}}\rangle$ is 0.6 dex below the solar value.

Figure 3 shows the synthetic ALMA and STIS observations. We do not expect perfect fidelity with the observations, but we can note some important similarities. As in the observation, the simulated halo contains two main star-forming galaxies with SFRs 189 and 115 M_⊙ yr⁻¹, respectively; these are well within the range of SFRs we estimate in Section 3.1. Both galaxies are submillimeter emitters, but only one has a flux density high enough to be detected in our mock ALMA map, with S₈₅₀ ≈ 0.5 mJy, within a factor of two of the observations. In the simulation, both galaxies are bright STIS sources, whereas this is only true of one of three observed ALMA components. The two central galaxies are surrounded by several STIS-detectable clumps (within ∼5'' or about 40 pkpc) that are of similar flux to those observed. Formally, we detect five sources at the >3σ level in the synthetic STIS image²⁷ , with flux densities of f_λ ≈ (0.2–4) × 10⁻¹⁹ erg s⁻¹ cm⁻² Å⁻¹. This is compared to eight sources with f_λ ≈ (0.2–2) × 10⁻¹⁹ erg s⁻¹ cm⁻² Å⁻¹ over the same area in the real image using identical detection criteria (Figure 1). In the real data, some of these sources could be chance projections, and we quantify this contamination by running the detection on the wider STIS image to measure the average surface density of sources in the same flux range. In a random 5'' radius aperture, we would expect 4 ± 1 sources by chance, so the predicted abundance of STIS sources in the simulation is in good agreement with the observations.

Zoom In Zoom Out Reset image size

What are these faint neighboring sources? In the simulation, these STIS clumps are the brightest of a population of star-forming sub-halos (that are apparent in the left panel of Figure 3) that "swarm" the central star-forming galaxy, bombarding it over a ∼1 Gyr period in a phase of hierarchical growth that accompanies the accretion of recycled gas previously heated by stellar feedback at an earlier epoch (Narayanan et al. 2015). So the simulation can broadly reproduce the ALMA and STIS observations, and we can interpret this as a demonstration that, as far as can be determined by the data, the halo population observed in SSA22-LAB01 is consistent with the prediction of this particular model of galaxy evolution for a halo of equivalent mass. What of the extended Lyα emission?

Considering the two main star-forming galaxies sources as Lyα emitters (with Lyα luminosities scaled from their SFRs, Equation (1)), we find that 58% of centrally generated photons escape the simulation box, and 16% of those scatter at least once in H i in the CGM. The majority (77%) of the H i mass in the central 400 pkpc of the simulation volume is tied to sub-halos, and thus it follows that the majority of the Lyα scattering occurs in and around satellites, which can be treated as cold clumps within the hot halo gas. Figure 4 shows the projected Lyα surface brightness maps of the simulated halo compared to the latest integral field observations from MUSE. Resonant scattering in the satellite population gives rise to extended Lyα emission with surface brightness distributions similar to observed LABs (${\mu }_{\mathrm{Ly}\alpha }\approx {10}^{-18}$ erg s⁻¹ cm⁻² arcsec⁻²). The observed Lyα surface brightness and morphology is highly orientation dependent Lyα sources—the emission need not be symmetric about the central sources that are the source of the Lyα photons. This is simiular to the situation in real LABs, for example, the next largest LAB in SSA22, SSA22-LAB02 contains an X-ray luminous AGN (Geach et al. 2009; Alexander et al. 2016) that is significantly offset from the peak and geometric center of the LAB. This has important implications for correctly identifying potential luminous galaxy counterparts to LABs, the apparent absence of which in some systems has previously been argued in favor of cold mode acceretion (e.g., Nilsson et al. 2006, but see also Prescott et al. 2015).

Zoom In Zoom Out Reset image size