Mapping the Lyα Emission around a z ∼ 6.6 QSO with MUSE: Extended Emission and a Companion at a Close Separation

The study of QSOs at z > 5.6 plays a central role in our understanding of how supermassive black holes (SMBHs) and galaxies form in the early universe. Currently, there are more than 170 known QSOs at z > 5.6 (e.g., Fan et al. 2006; Jiang et al. 2009; Willott et al. 2010; Bañados et al. 2014, 2015a, 2016; Carnall et al. 2015; Reed et al. 2015; Venemans et al. 2015a; Jiang et al. 2016; Matsuoka et al. 2016), only 12 of which are located at z > 6.5 (Mortlock et al. 2011; Venemans et al. 2013, 2015b; Matsuoka et al. 2016, 2017; Mazzucchelli et al. 2017a). The host galaxies of these very first QSOs are actively forming stars, with prodigious star-formation rates SFR > 100 ${M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Venemans et al. 2012, 2016), and were able to grow SMBHs with masses exceeding ${M}_{\mathrm{BH}}$ = 10⁹ ${M}_{\odot }$ in less than 800 Myr (De Rosa et al. 2014; Venemans et al. 2015b). The assembly of such massive SMBHs at early cosmic time requires that they accrete at the Eddington (or even super-Eddington) limit throughout a large fraction of their lifetimes (e.g., Yoo & Miralda-Escudé 2004; Volonteri & Rees 2005; Volonteri 2010, 2012; Madau et al. 2014; Volonteri et al. 2015). To sustain such vigorous accretion and intense star formation, the first QSOs need the presence of copious amounts of gas in their surroundings (e.g., Di Matteo et al. 2012; Dubois et al. 2012). Possibly, this gas is aggregated in dense flows able to penetrate into the virial radius of the halo and to funnel gas onto the central SMBH (e.g., Di Matteo et al. 2012; Feng et al. 2014). If the gas in the host galaxy and in the circumgalactic medium of a QSO is illuminated by the SMBH ionizing radiation and/or the intense starburst, then it may be observable as an extended "fuzz" of fluorescent Lyα emission (Rees 1988; Haiman & Rees 2001; Alam & Miralda-Escudé 2002).

Several of these Lyα nebulae have been reported in the literature at intermediate redshifts (z ∼ 2–4), leading to the general consensus that QSOs are frequently (50%–70%) embedded in nebulae with sizes of 10–100 kpc (e.g., Heckman et al. 1991a, 1991b; Christensen et al. 2006; Hennawi & Prochaska 2013; Roche et al. 2014; Herenz et al. 2015; Arrigoni Battaia et al. 2016). In the past years, gigantic nebulae, with projected sizes ≳300 kpc, have also been detected (Cantalupo et al. 2014; Martin et al. 2014; Hennawi et al. 2015; Borisova et al. 2016a), suggesting the presence of large amounts of cold gas around intermediate-redshift QSOs. Complementary, the analysis of absorption features in close projected QSO pairs confirms that QSO host galaxies are surrounded (with a covering fraction ∼60% within the virial radius) by cold (T ∼ 10⁴ K), metal-enriched (Z ≳ 0.1 ${Z}_{\odot }$) gas (Bowen et al. 2006; Hennawi et al. 2006a; Hennawi & Prochaska 2007, 2013; Decarli et al. 2009; Prochaska & Hennawi 2009; Farina et al. 2013, 2014; Prochaska et al. 2013a, 2013b, 2014; Johnson et al. 2015; Lau et al. 2016, 2017).

Despite the aforementioned achievements, the detection of these structures at z ∼ 6 is challenging due to the rapid decrease of the surface brightness (SB) with redshift $[\mathrm{SB}\propto {(1+z)}^{-4}]$. In recent years, huge efforts have been made to probe the extended Lyα emission around z ∼ 6 QSOs with contrasting results. Decarli et al. (2012), using the Wide Field Camera 3 (WFC3) narrowband filters on the Hubble Space Telescope (HST), put strong limits on the Lyα extended emission in the proximity of the highly star-forming host galaxies of the QSOs SDSS J1148+5251 (z = 6.42, $L(\mathrm{Ly}\alpha )$ < 2.5 × 10⁴⁴ erg s⁻¹) and SDSS J1030+0524 (z = 6.31, $L(\mathrm{Ly}\alpha )$ < 3.2 × 10⁴⁴ erg s⁻¹). Conversely, the presence of a Lyα nebula has been reported in narrowband Suprime-Cam/Subaru images of the QSO CFHQS J2329−0301 at z = 6.42 (Goto et al. 2009) and subsequently confirmed with long-slit spectroscopy with the Echelle Spectrograph and Imager and with the DEep Imaging Multi-object Spectrograph (DEIMOS) spectrographs at the Keck II telescope (Willott et al. 2011; Goto et al. 2012). Whereas all of these observations consistently report the presence of a bright Lyα halo extending on scales of 15 kpc in the proximity of CFHQS J2329−0301, its luminosity is not well-constrained, with values that range from $L(\mathrm{Ly}\alpha )$ ≳ 1.7–7.5 × 10⁴³ erg s⁻¹ (Willott et al. 2011; Goto et al. 2012) up to $L(\mathrm{Ly}\alpha )$ = 3.6 × 10⁴⁴ erg s⁻¹ (Goto et al. 2009; see footnote 10 in Decarli et al. 2012). Recently, Roche et al. (2014) presented long-slit spectroscopic observations with the Optical System for Imaging and low Resolution Integrated Spectroscopy (OSIRIS) mounted on the Gran Telescopio Canarias of a sample of QSOs at z > 2, including one at z = 5.95 the radio-loud QSO SDSS J2228+0110. Although bright, extended Lyα emission appears ubiquitous in the Roche et al. sample at z = 2–3, only a tenuous detection is reported for SDSS J2228+0110, with a luminosity of $L(\mathrm{Ly}\alpha )$ ≳ 7.8 × 10⁴² erg s⁻¹ extending up to a scale of ≳10 kpc. These values should be conservatively considered as lower limits. Indeed, they were computed by extracting the spectrum of the nebula over a small stripe close to the QSO emission. Additionally, a proper subtraction of the point spread function (PSF) was hindered by strong sky lines and by the low signal-to-noise ratio (S/N) of the spectrum.

The Integral-field Spectrograph Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) on the Very Large Telescope (VLT) is the obvious game changer in this kind of study. It produces a spatially resolved spectrum with a spatial sampling of 02 × 02 and a nominal spectral resolution (R = λ/${\rm{\Delta }}\lambda$) ranging from R = 1750 at 465 nm to R = 3750 at 930 nm, allowing the technical limitations of previous spectroscopic and narrowband investigations, such as uncertainties in the systemic redshift of the QSOs, filter and slit losses, and difficulties in performing a proper PSF subtraction, to be overcome.

In this paper, we present a deep MUSE integration aimed at detecting the Lyα nebular emission around the high-redshift QSO J0305−3150 (z = 6.61; with an absolute magnitude at 1450 Å of ${M}_{1450}=-25.96\pm 0.06$) discovered by Venemans et al. (2013) using the VISTA Kilo-Degree Infrared Galaxy (VIKING) survey. Sensitive near-infrared spectroscopy observations obtained with the Folded-port InfraRed Echellette spectrograph (FIRE) mounted on the Magellan Telescope revealed the presence of an SMBH with ${M}_{\mathrm{BH}}$= (0.95–1.20) × 10⁹ ${M}_{\odot }$ accreting with an Eddington ratio ${\lambda }_{\mathrm{Edd}}={L}_{\mathrm{Bol}}$/${L}_{\mathrm{Edd}}=0.68\mbox{--}0.74$ (De Rosa et al. 2014). The [${\rm{C}}\,{\rm{II}}$] 158 μm emission line $({L}_{[{\rm{C}}{\rm{II}}]}=(3.9\pm 0.2)\times {10}^{9}\,\,{L}_{\odot }$, FWHM = 255 ± 12 km s⁻¹) and the underlying far-infrared continuum were detected by Venemans et al. (2016) using the Atacama Large Millimeter/submillimeter Array (ALMA). From these measurements, the precise systemic redshift and SFR of the QSO host galaxy can be inferred. Although QSO redshifts derived from broad emission lines are subject to systematic shifts and large uncertainties (e.g., Richards et al. 2002; Bonning et al. 2007; Hewett & Wild 2010), especially at z > 6 (Venemans et al. 2016; Mazzucchelli et al. 2017a), the narrow [${\rm{C}}\,{\rm{II}}$] line accurately traces the systemic redshift (${z}_{\mathrm{sys}}=6.6145\pm 0.0001$). In addition, by fitting the FIR continuum spectrum with a modified blackbody with a spectral index β = 1.6 (after correcting for the impact of the cosmic microwave background; da Cunha et al. 2013), Venemans et al. (2016) obtained a dust temperature of 30 K and a total far-infrared luminosity of ${L}_{\mathrm{TIR}}$ = 2.6 × 10¹² $\,{L}_{\odot }$. The investigation of the spectral energy distribution (SED) of high-z QSOs suggests that the dust emission is predominantly powered by star formation especially at λ $\gtrsim$ 100 μm rest frame (e.g., Leipski et al. 2014; Barnett et al. 2015; see also Valiante et al. 2011). This implies⁹ SFR_TIR = 545 ${M}_{\odot }\,{\mathrm{yr}}^{-1}$. Aside from PSO J036+03 (Bañados et al. 2015b; Venemans et al. 2015b), J0305−3150 shows the highest SFR and Eddington ratio among the z > 6.5 QSOs known to date. It is thus an excellent target to constrain the properties of the gas reservoir that is expected to surround the first QSOs, together with its close environment.

Throughout this paper, we assume a concordance cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, ${{\rm{\Omega }}}_{{\rm{M}}}=0.3$, and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=1-{{\rm{\Omega }}}_{{\rm{M}}}\,=0.7$. In this cosmology, at z = 6.6145, the universe is 0.808 Gyr old, and an angular scale θ = 1'' corresponds a proper transverse separation of 5.4 kpc.

J0305−3150 was observed with MUSE on 2014 December 15 and 2015 January 15.¹⁰ The total time on source of 2.4 hr was divided into three observing blocks (OBs) of 48 m each, pointing to three different locations within 4'' from the QSO. For each OB, three 960 s exposures were taken, dithering using random offsets within a 75 box. The Differential Image Motion Monitor (DIMM) seeing during the observations was mostly sub-arcsecond, ranging from 07 to 11 (with a median of 08).

We employed the MUSE Data Reduction Software (version 1.0.1; Weilbacher et al. 2012, 2014) to perform, on each of the individual exposures, bias subtraction and flat-field, twilight, and illumination corrections, as well as wavelength and flux calibrations (the latter using the standard stars HD 49798 and GD 71 observed at the beginning of each observing night). Flat-field correction and sky subtraction of each exposure were improved using the CubExtractor package (S. Cantalupo 2017, in preparation). Exposures of each OB were sampled to a common grid (02 × 02 × 1.25 Å) and then average-combined. The absolute flux calibration for the resulting datacubes was obtained by rescaling the flux of the bright sources present in the MUSE field of view to our i-band images collected with the ESO Faint Object Spectrograph and Camera 2 (EFOSC2) on the New Technology Telescope (NTT; see Venemans et al. 2013 for details). Uncertainties on the absolute flux calibration resulting from this procedure are on the order of 5%. As a final step, the three datacubes were average-combined. The astrometry solution was refined by matching sources in the datacube with the first data release of the VIKING catalog (Edge et al. 2013).

The MUSE pipeline also provides a datacube containing errors formally propagated throughout the reduction process. As observed by Bacon et al. (2015), however, this process does not take into account correlations between neighboring voxels, ending up underestimating the real noise properties of the datacube. To have more realistic uncertainties, for each wavelength slice, the average of the variance delivered by the pipeline was rescaled to match the variance of the background (i.e., after removing the contribution from bright sources; see Borisova et al. 2016a for a similar approach).

In the final datacube, at 9256 Å (i.e., at the wavelength slice where the Lyα emission of the QSO is redshifted to), the FWHM of the PSF is 058, corresponding to 3.1 kpc at the QSO's redshift. The 5σ SB limit estimated after collapsing five wavelength slices centered at 9256 Å (i.e., from 9253.5 Å to 9258.5 Å) is ${\mathrm{SB}}_{5\sigma ,\lambda }^{1}=1.9\times {10}^{-18}$erg s⁻¹ cm⁻² arcsec⁻² over a 1 arcsec² aperture (see Section 3.2). The MUSE false color RGB image and the spectrum of J0305−3150 extracted over a radius of 3.7 spaxel (074) are shown in Figure 1. Flux errors on the spectrum (${\sigma }_{A,\lambda }$) are calculated from the final datacube as

$$\begin{eqnarray}&&{\sigma }_{A,\lambda }=\sqrt{\displaystyle \sum _{i\in A}{\sigma }_{\lambda ,i}^{2}},\end{eqnarray} \tag{ 1 }$$

where A is the area over which the spectrum is extracted in each wavelength slice (λ) and ${\sigma }_{\lambda ,i}^{2}$ is the variance in the corresponding spaxels.

Zoom In Zoom Out Reset image size

In the following sections, we describe the procedure adopted to investigate the presence of an extended Lyα emission (Section 3.1) and to estimate the sensitivity reached in the reduced datacube (Section 3.2). Finally, in Section 3.3, we present the results of this analysis.

3.1. PSF Subtraction

An accurate PSF subtraction is necessary to recover the faint signal of the diffuse Lyα halo emerging from the PSF wings of the bright unresolved nuclear component. To perform this task, we created two empirical PSF models: PSF_QSO, constructed directly from the QSO emission by collapsing regions of the spectrum virtually free from any extended emission¹¹ (i.e., away from the Lyα emission; see Figure 1), and PSF_STAR, obtained from the bright star located ∼25'' northwest from the QSO by summing up its emission over the wavelength range where the extended Lyα emission is expected to fall (see Figure 1). These two PSF models are subject to different systematics, allowing us to check for the reliability of a possible detection of extended emission. The first model allows us to directly subtract the PSF contribution from the QSO without any spatial shift. However, PSF_QSO has a relatively low S/N due to the faintness of the source and the small range in wavelength used (starting 325 km s⁻¹ away from the QSO systemic redshift, i.e., at 9266 Å, up to the wavelength where the presence of strong sky emission lines drastically increase the variance, i.e., 9301 Å; see Figure 1). In addition, it may be contaminated by the wings of the possible extended emission if they are particularly broad and/or redshifted. PSF_STAR instead benefits from a higher S/N, and it minimizes PSF changes with wavelength. On the other hand, it is subject to resampling to centroid the PSF model on the QSO, and to the the spatial variation of the PSF.

To subtract the unresolved QSO emission and to recover the Lyα nebula, we adapted the technique used by Hennawi & Prochaska (2013) and Arrigoni Battaia et al. (2015a) to the three-dimensional structure of the MUSE data. First, at each wavelength slice, the PSF model is rescaled to the QSO's flux estimated in a circle with radius 2 spaxel (04). The underlying assumption is that the QSO is dominating the emission in this central region. A bright, centralized nebular component may, however, lead to an overestimate of the QSO's emission and thus to an underestimate of the total flux of the possible extended emission. Then, we defined the ${\chi }_{\lambda ,i}$ datacube,

$$\begin{eqnarray}&&{\chi }_{\lambda ,i}=\displaystyle \frac{{\mathrm{DATA}}_{\lambda ,i}-{\mathrm{MODEL}}_{\lambda ,i}}{{\sigma }_{\lambda ,i}},\end{eqnarray} \tag{ 2 }$$

where the indices λ and i indicate the wavelength slice and the 2D spaxel position, respectively; ${\mathrm{DATA}}_{\lambda ,i}$ is the datacube; ${\mathrm{MODEL}}_{\lambda ,i}$ is the rescaled PSF model; and ${\sigma }_{\lambda ,i}$ is the square root of the variance datacube. If our model accurately describes the PSF (and in the absence of systematics), at each wavelength slice the distribution of ${\chi }_{\lambda ,i}$ values should follow a Gaussian centered on zero with unit variance. Under this condition, this datacube thus permits the statistical significance of any putative detection to be assessed. Note that both PSF_QSO and PSF_STAR have, by construction, a much higher S/N than the QSO in a single wavelength slice; therefore, the contribution of the PSF model to the variance budget is negligible. We also constructed a smoothed datacube, $\mathrm{SMOOTH}[{\chi }_{\lambda ,i}]$, that is helpful for identifying the possible presence of extended emission,

$$\begin{eqnarray}&&\mathrm{SMOOTH}[{\chi }_{\lambda ,i}]=\displaystyle \frac{\mathrm{CONVOL}[{\mathrm{DATA}}_{\lambda ,i}-{\mathrm{MODEL}}_{\lambda ,i}]}{\sqrt{{\mathrm{CONVOL}}^{2}[{\sigma }_{\lambda ,i}^{2}]}},\end{eqnarray} \tag{ 3 }$$

where $\mathrm{CONVOL}$ indicates a convolution in the spatial axis with a 2D Gaussian kernel with σ = 1 spaxel, while ${\mathrm{CONVOL}}^{2}[{\sigma }_{\lambda ,i}^{2}]$ is the convolution of the variance datacube (${\sigma }_{\lambda ,i}^{2}$) with the square of the same kernel. In the absence of systematics, spaxel values in this smoothed images should still follow a Gaussian distribution, but with smaller variance due to the increased correlation among pixels. The extended Lyα nebulae around radio-quiet QSOs typically show quiescent kinematics with FWHM ≲ 600 km s⁻¹ (Arrigoni Battaia et al. 2015b; Borisova et al. 2016a). We therefore repeat the same procedure, binning over 3, 5, 10, and 15 wavelength slices in order to maximize the S/N of the emission with FWHMs of ∼120, 200, 400, and 600 km s⁻¹. The resulting ${\chi }_{\lambda ,i}$ and $\mathrm{SMOOTH}[{\chi }_{\lambda ,i}]$ datacubes were then visually inspected to search for the presence of any extended emission.

3.2. Detection Limits

The detection limit for a Lyα nebular emission in the MUSE datacube depends on its FWHM and physical size: lower SB levels can be reached by averaging in space and/or in velocity. Under the (erroneous) assumption that spaxels are independent, the theoretical detection limits for an extended source in a single wavelength slice is given by ${\mathrm{SB}}_{1\sigma ,\lambda }^{A}={\mathrm{SB}}_{1\sigma ,\lambda }$/$\sqrt{{\#}_{A}}$, where ${\mathrm{SB}}_{1\sigma ,\lambda }$ is the 1σ SB detection limit per 02 × 02 spaxel at the wavelength slice λ and ${\#}_{A}$ is the number of spaxels in the isophotal area of the source A. The corresponding limit on the total flux (and thus on the luminosity) can be written as ${F}_{1\sigma ,\lambda }^{A}={\mathrm{SB}}_{1\sigma ,\lambda }\sqrt{{\#}_{A}}{\mathrm{PS}}^{2}$, where $\mathrm{PS}$ is the pixel scale of MUSE: PS = 02 spaxel⁻¹. Binning the MUSE datacube over 3, 5, 10, and 15 wavelength slices centered at 9256 Å (i.e., at the expected position of the Lyα emission), the formal 5σ SB detection limits calculated over an aperture of 1 arcsec² are ${\mathrm{SB}}_{5\sigma ,\lambda }^{1}=[1.3$, 1.9, 3.4, 5.0] × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻², respectively.

The noise properties of the MUSE datacube are, however, not uniform. In addition, cross-talk between voxels and systematics introduced during data reduction and PSF subtraction will alter these theoretical detection limits. We tested the reliability of the calculated SB limits by introducing a set of synthetic sources in the PSF-subtracted datacubes and visually estimated the level of a convincing detection. For this purpose, we focused on the wavelength region where a Lyα line redshifted to ${z}_{\mathrm{sys}}$ would fall (i.e., at 9256 Å), and we binned the datacube over the same wavelength slices used for the sources detection (see Section 3.1). We randomly placed mock circular sources (including Poisson noise) with a top-hat SB in different locations in each pseudo-narrowband images. These mock sources have total integrated luminosities L_mock = 10^{[41.5, 42.2, 42.9, 43.6, 44.3]} erg s⁻¹, diameters d_mock = [3, 5, 25, 50, 100] kpc, and, in the wavelength space, a Gaussian distribution with FWHM_mock = [80, 120, 200, 400, 600] km s⁻¹ (where 80 km s⁻¹ is, roughly, the nominal resolution limit and the other values match the binning considered in the PSF subtraction process; see Section 3.1). A diameter of ${d}_{\mathrm{mock}}=3$ kpc corresponds to the seeing measured in the datacube and is hence unresolved.

As an illustrative example, we show in Figure 2 how synthetic sources with FWHM_mock = 120 km s⁻¹ and different ${d}_{\mathrm{mock}}$ and ${L}_{\mathrm{mock}}$ would appear in a pseudo-narrowband image obtained by collapsing three wavelength slices around the expected position of the Lyα line (i.e., from 9254.75 Å to 9257.25 Å). The corresponding ${\chi }_{i}$ and $\mathrm{SMOOTH}[{\chi }_{i}]$ images (Figure 3) show that we should be able to visually detect these sources down to $L(\mathrm{Ly}\alpha )$ ∼ 10^42.0 erg s⁻¹ if unresolved and down to $L(\mathrm{Ly}\alpha )$ ∼ 10^44.0 erg s⁻¹ if the emission is more extended.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If sources fall on the top of the QSO emission, the PSF subtraction process will hinder the achievement of the theoretical detection limits. However, the QSO is relatively faint in the datacube, and PSF wings quickly drop below the 2σ SB limit. For instance, in the pseudo-narrowband obtained by collapsing the datacube between 9238.5 Å and 9242.0 Å, the QSO radial profile becomes consistent with zero (at 2σ) at a separation of ∼12 (i.e., ∼6.5 kpc at ${z}_{\mathrm{sys}}$; see Figure 4). Contamination due to imperfect PSF subtraction are therefore expected to impact only sub-arcsecond separations.

Zoom In Zoom Out Reset image size

3.3. Tentative Detection of an Extended Emission

Figure 5 shows the result of the PSF subtraction procedure described above. A Lyα extended emission with a size of ∼16 (∼9 kpc) toward the southwest of J0305−3150 is tentatively detected in the ${\mathrm{DATA}}_{i}-{\mathrm{MODEL}}_{i}$ pseudo-narrowband image obtained by collapsing the three wavelength slices from 9259.75 Å to 9262.25 Å. After removing a circle of 04 radius centered on the QSO's position (to avoid possible contamination due to imperfect PSF subtraction,) we measured the flux integrated over an elliptical aperture, with semiminor and semimajor axes of 09 and 14. We obtained $F(\mathrm{Ly}\alpha )$ = (6.1 ± 0.8) × 10⁻¹⁸ erg s⁻¹ cm⁻² ($F(\mathrm{Ly}\alpha )$ = (5.6 ± 0.8) × 10⁻¹⁸ erg s⁻¹ cm⁻²⁾ using PSF_QSO (PSF_STAR) as our model, yielding a 7.6σ (7.0σ) detection. The inferred luminosity of the extended emission is $L(\mathrm{Ly}\alpha )$ = (3.0 ± 0.4) × 10⁴² erg s⁻¹ (using PSF_QSO), more than one order of magnitude fainter with respect to the characteristic luminosity of LAEs at z ∼ 6.6 (e.g., Hu et al. 2010; Matthee et al. 2015).

Zoom In Zoom Out Reset image size

In order to assess the reliability of this detection, we empirically estimate the effects of systematics in the pseudo-narrowband image obtained by collapsing the ${\mathrm{DATA}}_{\lambda ,i}-{\mathrm{MODEL}}_{\lambda ,i}$ datacube over the wavelength range 9259.75 Å to 9262.25 Å (see Arrigoni Battaia et al. 2015a for a similar test performed on narrowband images). From this image, we extracted the fluxes and variances for a set of elliptical apertures, randomly picking background locations (i.e., avoiding bright sources) in an annulus with internal and external radii of 5'' and 25'' from the QSO, respectively. In the absence of systematics, the flux over noise ratio of these apertures should follow a Gaussian distribution centered on zero and with σ = 1. Figure 6 shows that the actual distribution is nearly Gaussian with an offset of −0.28 and a sigma 1.42× broader. These deviations could be due to poor sky subtraction and/or to 3D correlations present in the MUSE datacube. Taking into account systematics, the significance of our detection is reduced to 5.3σ (4.9σ) using PSF_QSO (PSF_STAR) as the PSF model.

The spectrum of this possible nebula extracted from the ${\mathrm{DATA}}_{\lambda ,i}-{\mathrm{MODEL}}_{\lambda ,i}$ datacube over the same elliptical aperture considered above is shown in Figure 7. The possible Lyα halo appears as the strongest feature present within ±1000 km s⁻¹ from the QSO's systemic redshift. This narrow emission-line peaks at λ = 9261.5 Å and has an FWHM = 65 km s⁻¹. The fitted width is slightly smaller than the nominal resolution limit of MUSE at these wavelengths, but consistent given the low S/N per spectral bin of the line. Under the assumption that we are probing the Lyα emission, the line appears redshifted by 155 km s⁻¹ with respect to the systemic redshift traced by the [${\rm{C}}\,{\rm{II}}$] line. This is smaller than the 445 km s⁻¹ shift observed between the Hα emission (originating from H ii regions in the galaxy) and the Lyα line in low-redshift UV-selected galaxies (e.g., Steidel et al. 2010), but in agreement with the 175 km s⁻¹ offset seen in strong Lyα emitters (LAE) at z ∼ 2 (e.g., Hashimoto et al. 2013). A similar shift is also observed at higher redshift. For instance, Pentericci et al. (2016) observed shifts of 100–200 km s⁻¹ between the Lyα and [${\rm{C}}\,{\rm{II}}$] emission line in four galaxies at redshifts between 6.6 and 7.2. Likewise, the Lyα emission of the z ∼ 6 star-forming galaxy A383-5.2 appears to be shifted by 120 km s⁻¹ with respect to the systemic redshift probed via the ${\rm{C}}\,{\rm{III}}$] λ1909 metal line (Stark et al. 2015).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The rapid formation of SMBHs in the early universe may imply that the first QSOs are tracers of galaxy overdensities. In fact, the black hole growth may be fostered by rich environments, where interactions and mergers are more likely to occur (see Volonteri 2012 for a review). To test this hypothesis, we performed a search of LAEs in the proximity of J0305−3150.

To identify LAEs associated with the QSO, we focused our attention on the wavelength range between ∼9226 Å and ∼9287 Å, corresponding to ±1000 km s⁻¹ from the QSO systemic redshift. This region was recursively sliced in 10 Å wide (i.e., ∼325 km s⁻¹) pseudo-narrowband images. The nominal 5σ SB limit reached in the image centered at 9256 Å is ${\mathrm{SB}}_{5\sigma ,\lambda }^{1}=2.3\times {10}^{-18}$ erg s⁻¹ cm⁻² arcsec⁻² over a 1 arcsec² aperture. Assuming a spatially unresolved line emission with FWHM of 270 km s⁻¹ (i.e., the FWHM of the z ∼ 6.6 LAE spectral template presented in Ouchi et al. 2010), the corresponding 5σ luminosity limit is ${L}_{5\sigma }=6.7\times {10}^{41}$ erg s⁻¹. Each image was then processed using SExtractor (Bertin & Arnouts 1996) requiring a minimum detection area of 7 spaxel and a detection threshold of 1.5σ. Identified sources were considered only if located within a 50'' × 50'' box centered on the QSO. This search box, smaller than the full MUSE field of view, was chosen to avoid issues related to the shorter exposure times experienced by the peripheral regions due to the dithering. The neutral hydrogen in the intergalactic medium is expected to suppress virtually all of the flux blueward of the Lyα line in high-redshift galaxies. We thus excluded from our analysis sources also present in the pseudo-narrowband image obtained by collapsing ${\mathrm{DATA}}_{\lambda ,i}$ over the wavelength slices 4750 Å–8500 Å. Finally, single-exposure datacubes were visually inspected to identify false detections associated with highly deviating pixels present in only one OB.

This procedure allowed us to reveal the presence of an LAE in the immediate proximity of J0305−3150 (R.A. = 03:05:16.80, decl. = −31:50:57.3, epoch = J2000; see Figure 8). This appears as an 8.3σ detection in the pseudo-narrowband images created by summing up slices in the ${\mathrm{DATA}}_{\lambda ,i}$ datacube from 9269.75 Å to 9277.25 Å, and as a 4.2σ, 3.5σ, and 5.1σ detection in the single OB datacubes collapsed in the same wavelength range (see Figure 8). Figure 9 shows the spectrum extracted over a circular aperture with a radius of 3 spaxel. By fitting a Gaussian function over the most prominent emission line, we derive a redshift of ${z}_{\mathrm{LAE}}=6.629$ (i.e., redshifted by ∼560 km s⁻¹ with respect to the QSO's systemic redshift), an ${\mathrm{FWHM}}_{\mathrm{LAE}}=240$ km s⁻¹, and a luminosity $L(\mathrm{Ly}\alpha )$ = (2.1 ± 0.2) × 10⁴² erg s⁻¹. Consistently, the luminosity directly estimated from the pseudo-narrowband image is $L(\mathrm{Ly}\alpha )$ = (1.8 ± 0.2) × 10⁴² erg s⁻¹. No significant continuum emission is detected redward of the Lyα line.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To quantify the effect of systematics, we performed the same empirical test used in Section 3.3. From the pseudo-narrowband image, we extracted an ensemble of circular apertures with 3 spaxel radius, avoiding bright sources and image edges. The distribution of the flux over noise ratio of these apertures is well-reproduced by a Gaussian with average μ = 0.31 and sigma σ = 1.21 (see Figure 10). The significance of the detection of the LAE emitter, once systematics are taken into account, is thus 6.8σ. In Figure 9, we compare the spectra extracted from the circular apertures considered above with the spectrum of the LAE. Despite the relatively low S/N per voxel of the line, seven consecutive slices have S/N > 2 and the whole emission is the brightest line detected within ±1000 km s⁻¹ from the QSO's systemic redshift. These results corroborate the reliability of the LAE's detection. Given the high redshift of this source, MUSE is not able to cover additional rest-frame UV line diagnostics (e.g., $\mathrm{He}\,{\rm{II}}$, ${\rm{C}}\,{\rm{IV}}$), which would allow us to better determine the strength (and nature) of this faint companion.

Zoom In Zoom Out Reset image size

In the following, we discuss the implications of the tentative detection of the extended Lyα emission associated with the high-redshift QSO J0305−3150 (Sections 5.1 and 5.2) and of its companion galaxy (Section 5.3). It is worth mentioning that the real luminosity of the Lyα emission could be slightly underestimated due to the presence of neutral hydrogen in the proximity of the QSO. This effect is, however, negligible in the context of the forthcoming discussion.

5.1. The Host Galaxy of J0305−3150

In Section 3.3, we present a possible detection of extended Lyα emission around the QSO J0305−3150. The significance of this emission is strengthened by its spatial and redshift positions (close to what is expected from the [${\rm{C}}\,{\rm{II}}$] emission). In this section, we link this discovery with the properties of the QSO's host galaxy.

Figure 11 maps the distribution of [${\rm{C}}\,{\rm{II}}$] (as observed with ALMA) overplotted on the Lyα halo detected with MUSE. The two emissions do not appear to be co-spatial. A displacement of the Lyα and [${\rm{C}}\,{\rm{II}}$]/dust emission has been commonly observed in other dusty sources at high redshift. For instance, Hodge et al. (2015) reported a ∼4 kpc offset between the rest-frame UV and both the FIR and the CO emission in the z = 4.05 submillimeter galaxy GN 20. Decarli et al. (2016) and Aravena et al. (2016) observed a comparable shift in the compact star-forming galaxy ASPECS C.1 (z = 2.54). In a similar fashion, the Lyα extended emission not associated with the radio jet, the rest-frame UV continuum, and the dust emission (albeit the coarse spatial resolution) appear to be offset in the z = 4.11 radio galaxy TN J1338−1942 (e.g., De Breuck et al. 2004; Zirm et al. 2005; Venemans et al. 2007; Swinbank et al. 2015).

Zoom In Zoom Out Reset image size

For J0305−3150, we can estimate the obscuration due to the copious amounts of dust detected with ALMA (M_dust = (4.5–24) × 10⁸ ${M}_{\odot }$; Venemans et al. 2016) by comparing the luminosity of the observed extended Lyα emission with the theoretical (unobscured) Lyα emission expected due to the UV photons coming from the intense starburst detected at millimeter wavelengths. Assuming a case B recombination, the relation between $L({\rm{H}}\alpha )$ and SFR (Kennicutt & Evans 2012) becomes

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

The star-formation rate of the host galaxy (SFR_TIR = 545 ${M}_{\odot }\,{\mathrm{yr}}^{-1}$) yields a Lyα luminosity of $L(\mathrm{Ly}\alpha )$ = 8.8 × 10⁴⁴ erg s⁻¹, i.e., a factor ∼300× brighter than the observed luminosity. This implies high extinction with ${A}_{\mathrm{UV}}\sim 6.2$ mag.

Resonance scattering may trap Lyα photons and dim the expected luminosity of the extended emission. The bouncing of Lyα photons between optically thick clouds increases the total path length traveled in the dusty medium, incrementing its extinction (see the discussion in Decarli et al. 2012 and the Appendix in Hennawi & Prochaska 2013). The relative importance of this effect strongly depends on various QSO host galaxy properties, such as neutral hydrogen column density, neutral fraction, dust-to-gas ratio, and geometry, among others. However, there is no clear evidence for a broad double-peaked kinematics as expected from resonantly trapped Lyα photons (e.g., Cantalupo et al. 2005), possibly due to absorption from neutral hydrogen in the intergalactic medium. In addition, a high level of obscuration is in contrast with the low extinction observed toward the QSO: Figure 12 plots the broadband photometry for J0305−3150 taken from Venemans et al. (2013). The SED matches well the QSO composite spectrum of Selsing et al. (2016), leaving room for only little extinction. Assuming an SMC extinction curve (Gordon et al. 2003), which is appropriate for non-BAL QSOs at high redshift (e.g., Gallerani et al. 2010), we can infer a stringent limit on the obscuration of ${A}_{\mathrm{UV}}\lt 0.1$ mag.

Zoom In Zoom Out Reset image size

These results suggest that the copious amounts of dust detected with ALMA (M_dust = (4.5–24) × 10⁸ ${M}_{\odot }$; Venemans et al. 2016) can effectively obscure the UV emission and prevent the ionizing photons from escaping the QSO's host galaxy. However, UV radiation is able to leak through small openings in the dust cocoon. These leakages may be associated with short and extreme bursts of star formation that have been proposed to explain the intrinsically blue rest-frame UV slope observed in dusty, star-forming galaxies (SFR > 50 M_⊙ yr⁻¹) up to z ∼ 5 (e.g., Casey et al. 2014). A similar scenario was also proposed by Decarli et al. (2012) to explain the lack of extended Lyα emission from the host galaxy of the two highly star-forming z > 6 QSOs: SDSS J1030+0524 and SDSS J1148+5251 (see also Mechtley et al. 2012).

5.2. Large-scale Lyα Emission

Is the SB limit reached with MUSE sufficient to probe the presence of Lyα nebular emission on scales of 5–100 kpc? To test this possibility, we compare our detection limits with the circularly averaged SB profiles of bright QSOs at intermediate redshift. As a comparison sample, we consider the recent work of Borisova et al. (2016a), who investigated 17 radio-quiet QSOs at 3.0 ≲ z ≲ 3.9. We stress that, with a median absolute magnitude at 1450 Å of ${M}_{1450}=-29.2$, these are among the brightest QSOs at z ∼ 3.5, hence they are ∼3 mag brighter than J0305−3150. The nebulae surrounding these QSOs extend on scales of 50–100 kpc and their circularly averaged SB profiles are well-represented by a power-law decline in the majority of the cases. In Figure 13, we plot the power-law best fits of these extended Lyα emissions as they would appear if moved at $z={z}_{\mathrm{sys}}$, which means we corrected for redshift dimming and for different angular diameter distances. We also show the profile of the extended Lyα emission associated with J0305−3150 averaged over annuli evenly spaced in logarithmic space together with our nominal 5σ limit on the SB. It is apparent that if J0305−3150 were surrounded by an extended emission similar to bright z ∼ 3.5 QSOs, our MUSE data would have been deep enough to detect it at angular scales θ ≳ 1''.

Zoom In Zoom Out Reset image size

The nebular emission around J0305−3150 appears to have a significantly lower SB and lower total luminosity than commonly observed at intermediate redshifts. In Figure 14, we present a compilation of Lyα nebulae detected around QSOs from the literature.¹² The 5σ upper limits on the luminosity set by Decarli et al. (2012) on SDSS J1030+0524 (z = 6.31) and SDSS J1148+5251 (z = 6.42) are also plotted. These are derived by rescaling their 5σ SB limit (${\mathrm{SB}}_{5\sigma }^{1}=1.0\times {10}^{-17}$ erg s⁻¹ cm⁻² arcsec⁻² over a 1 arcsec² aperture) to a circular aperture with a diameter of 10 kpc. On the right side of Figure 14, we mark the 5σ upper limits on the luminosity for top-hat sources with diameters d = [10, 25, 50, 100] kpc and FWHM = 200 km s⁻¹. These are calculated considering the SB limit obtained by collapsing the MUSE datacube over five wavelength slices (${\mathrm{SB}}_{5\sigma ,\lambda }^{1}=1.9\,\times {10}^{-18}$ erg s⁻¹ cm⁻² arcsec⁻² over a 1 arcsec² aperture; see Section 3.2). These data suggest a decline in the total luminosity of the Lyα nebulae as a function of redshift. This may indicate a change in the gas properties and/or in the powering mechanisms at different epochs. A more quantitative interpretation, however, is hampered by the differing methodologies, ambiguities in detection criteria, and lack of a statistical sample of QSOs investigated at z > 4.

Zoom In Zoom Out Reset image size

It is of interest to compare our observational results with predictions on the luminosity of the Lyα for different emission mechanisms. In Section 5.2.1, we consider the case of recombination from optically thick clouds. The case in which the QSO radiation highly ionizes the surrounding gas which thus becomes optically thin is addressed in Section 5.2.2. Finally, in Section 5.2.3, we comment on the possibility that the gas is in a multiphase status. These calculations closely follow the formalism described in Hennawi & Prochaska (2013).

5.2.1. Optically Thick Gas

Under the assumption that the surrounding of the QSO is filled with cool, optically thick clouds, self-shielding generates a thin, highly ionized envelope around individual clouds that acts as a mirror converting a fraction of the ionizing radiation into Lyα photons (Gould & Weinberg 1996). In this scenario, the powering mechanism is the QSO ionizing radiation (${L}_{{\nu }_{\mathrm{LL}}}$, where ${\nu }_{\mathrm{LL}}$ is the frequency at the Lyman edge). $L(\mathrm{Ly}\alpha )$ is thus proportional to ${L}_{{\nu }_{\mathrm{LL}}}$ as

$$\begin{eqnarray}&&\displaystyle \frac{L(\mathrm{Ly}\alpha )}{{10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}=7.8\,{f}_{c}^{\mathrm{thick}}\displaystyle \frac{{L}_{{\nu }_{\mathrm{LL}}}}{{10}^{30}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}},\end{eqnarray} \tag{ 5 }$$

where ${f}_{c}^{\mathrm{thick}}$ is the optically thick cloud-covering factor. If we assume ${f}_{c}^{\mathrm{thick}}=0.1$ as estimated for the small-scale Lyα emission observed in z = 2–3 QSOs by Hennawi & Prochaska (2013) and ${L}_{{\nu }_{\mathrm{LL}}}=4.9\times {10}^{30}$ erg s⁻¹ Hz⁻¹ (obtained by rescaling the composite spectrum from Lusso et al. 2015 to the QSO luminosity at 1350 Å), we obtain $L(\mathrm{Ly}\alpha )$ = 3.8 × 10⁴⁴ erg s⁻¹. This discrepancy of a factor ∼130× with respect to the observed luminosity may be due to the geometry of the emission: if the UV photons break only through a small solid angle ${{\rm{\Omega }}}_{{\rm{e}}}$, the expected Lyα luminosity would be reduced by a factor ${f}_{{{\rm{\Omega }}}_{{\rm{e}}}}=\left(\tfrac{{{\rm{\Omega }}}_{{\rm{e}}}}{4\pi }\right)$. Assuming, for the sake of simplicity, that the anisotropic emission occurs in a cone, an opening angle of 30° corresponds to ${f}_{{{\rm{\Omega }}}_{{\rm{e}}}}\sim 0.07$. However, this would imply an unrealistic fraction of obscured AGNs of ${f}_{\mathrm{Obs}}=\left(1-\tfrac{{{\rm{\Omega }}}_{{\rm{e}}}}{4\pi }\right)=0.85$ (e.g., Treister et al. 2008; Lusso et al. 2013; Merloni et al. 2014). Alternatively, a factor of 100× lower optically thick cloud-covering fraction may explain the faintness of the emission. Such a low covering fraction, however, is in contrast with results from z ∼ 2 to 3 QSOs (Prochaska et al. 2013a; see also Section 5.2.3). This result pushes for different emission mechanisms to explain the observed emission.

5.2.2. Optically Thin Gas

If the gas surrounding the QSO is optically thin, the QSO radiation would be sufficiently intense to keep the gas highly ionized (i.e., the hydrogen neutral fraction is ${x}_{{\rm{H}}{\rm{I}}}=\tfrac{{n}_{{\rm{H}}{\rm{I}}}}{{n}_{{\rm{H}}}}\ll 1$). As shown in Hennawi & Prochaska (2013), in the optically thin regime, $L(\mathrm{Ly}\alpha )$ can be expressed in terms of the area-averaged neutral column density ($\langle {N}_{{\rm{H}}{\rm{I}}}\rangle$) and of the ionizing luminosity,

$$\begin{eqnarray}&&\displaystyle \frac{L(\mathrm{Ly}\alpha )}{{10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}=0.9\,\displaystyle \frac{\langle {N}_{{\rm{H}}{\rm{I}}}\rangle }{{10}^{17.2}\,{\mathrm{cm}}^{-2}}\displaystyle \frac{{L}_{{\nu }_{\mathrm{LL}}}}{{10}^{30}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}},\end{eqnarray} \tag{ 6 }$$

where the normalization of $\langle {N}_{{\rm{H}}{\rm{I}}}\rangle$ is set by the requirement that, to be optically thin, clouds must have ${N}_{{\rm{H}}{\rm{I}}}\ll {10}^{17.2}$ cm⁻². Plugging into this equation the observed luminosity of the nebula and the ${L}_{{\nu }_{\mathrm{LL}}}$ estimated above, we obtain $\langle {N}_{{\rm{H}}{\rm{I}}}\rangle \,\sim {10}^{15.0}$ cm⁻². It is thus plausible that the extended emission arises from optically thin clouds illuminated by the QSO. However, it is worth a reminder that while $\langle {N}_{{\rm{H}}{\rm{I}}}\rangle \gt {10}^{17.2}$ cm⁻² implies that the gas is in the optically thick regime, clouds with $\langle {N}_{{\rm{H}}{\rm{I}}}\rangle \lt {10}^{17.2}$ cm⁻² could be either optically thin or optically thick (see Hennawi & Prochaska 2013).

If we assume an optically thin scenario, we can relate the observed Lyα SB to the hydrogen total column density (${N}_{{\rm{H}}}$) and volume density (${n}_{{\rm{H}}};$ e.g., Hennawi & Prochaska 2013),

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{SB}(\mathrm{Ly}\alpha )}{{10}^{-19}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{arcsec}}^{-2}}=5.9{\left(\displaystyle \frac{1+z}{7.6145}\right)}^{-4}\\ &&\quad \times \,{f}_{c}^{\mathrm{thin}}\left(\displaystyle \frac{{n}_{{\rm{H}}}}{1.0\,{\mathrm{cm}}^{-3}}\right)\left(\displaystyle \frac{{N}_{{\rm{H}}}}{{10}^{20.5}\,{\mathrm{cm}}^{-2}}\right),\end{eqnarray} \tag{ 7 }$$

where ${f}_{c}^{\mathrm{thin}}$ is the covering fraction of the optically thin gas. From this equation, we can derive an estimate of the volume density of the gas giving rise to the extended emission. In fact, the total column density of the hydrogen in the proximity of z ∼ 2 QSOs has been constrained from the study of absorption systems. Photoionization models of these absorbers suggest that ${N}_{{\rm{H}}}$ is almost constant within an impact parameter of 200 kpc with a median value ${N}_{{\rm{H}}}={10}^{20.5}$ cm⁻² (e.g., Prochaska & Hennawi 2009; Hennawi et al. 2015; Lau et al. 2016). Assuming ${f}_{c}^{\mathrm{thin}}=0.5$, to explain the observed SB SB(Lyα) = (1.8 ± 0.2) × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² (calculated over the elliptical aperture considered in Section 3.3), a high gas volume density of ${n}_{{\rm{H}}}=6.1$ cm⁻³ is required. Remarkably, similarly high ${n}_{{\rm{H}}}$ were proposed to explain the emission of the giant Lyα nebulae associated with the QSOs UM 287 (Cantalupo et al. 2014; Arrigoni Battaia et al. 2015b) and SDSS J0841+3921 (Hennawi et al. 2015).

5.2.3. A Multiphase Scenario

In the previous sections, we estimated the expected emission from the gas surrounding J0305−3150 considering that it is either optically thick or optically thin to QSO radiation. However, we can also consider a multiphase scenario where low-density clouds with a high covering fraction are responsible for the absorption systems observed at z ∼ 2, whereas the observed emission rises from optically thick gas with a low covering fraction (and hence rarely intercepted in absorption). We also stress that disentangling the different emission mechanism is challenging. At a given separation from the QSO, optically thick and optically thin gas clouds could result in a similar emission with opportune combinations of ${n}_{{\rm{H}}}$, ${N}_{{\rm{H}}}$, and covering fraction. From Equations (6) and (7), it is clear that, for a given ${N}_{{\rm{H}}}$, one can consider increasing ${n}_{{\rm{H}}}$ to increase the emission in the optically thin regime. However, when the area-averaged neutral column density reaches $\langle {N}_{{\rm{H}}{\rm{I}}}\rangle \,\sim {10}^{17.2}$ cm⁻², the cloud becomes optically thick and recombination would occur only in the self-shielding layer. We can thus roughly calculate the density required for a cloud to become optically thick as a function of the distance from the QSO by matching the emission estimates in Sections 5.2.1 and 5.2.2 (for the same covering factor). We obtain that, in order to shelf-shield at a separation of 25 kpc from J0305−3150, a cloud should have a high density of ${n}_{{\rm{H}}}\sim 200$ cm⁻³.

5.3. Overdensity of LAEs around the J0305−3150

The first search for LAEs around z > 5.5 QSOs was performed by Decarli et al. (2012) using a combination of narrowband filters of the WFC3 on HST. No companion galaxies were found around the two QSOs SDSS J1030+0524 and SDSS J1148+5251. However, the field of view of WFC3 allowed the authors to probe only a small cosmological volume and a relatively bright point source detection limit was reached (∼23.4 mag). Bañados et al. (2013) and Mazzucchelli et al. (2017b) used sensitive narrow- and broadband images from the FOcal Reducer/low dispersion Spectrograph 2 (FORS2) on the VLT to investigate the environment of the two QSOs ULAS J0203+0012 (z = 5.72) and PSO J215.1512−16.0417 (z = 5.73) over a much larger area (∼37 arcmin²). Both studies report a number of LAEs consistent with (or even lower than) expectations from a blank field. Goto et al. (2017), using the Subaru Prime Focus Camera (Suprime-cam) on the Subaru telescope, similarly reported an underdensity of LAEs around the QSO CFHQS J2329−0301 (z = 6.4). It is worth mentioning that at slightly lower redshift, McGreer et al. (2014) discovered an LAE in the immediate proximity (with a separation of only 11.4 kpc) of the z ∼ 4.9 QSO SDSS J0256+0019.

The detection of one LAE in the MUSE datacube at only ∼12.5 kpc and 560 km s⁻¹ from J0305−3150 may represent the first spectroscopically confirmed evidence of the connection between high-density environments (and thus high merger rate) and z > 6 QSOs, which was postulated to explain the rapid assembly of the first SMBHs (e.g., Volonteri 2012). To confirm this scenario, we have to calculate first the probability of finding such LAEs in the proximity of the QSO.

5.3.1. Comparison to Blank Field

In order to estimate how many LAEs are expected from the blank field, we integrated the z ∼ 6.6 LAE luminosity function from Matthee et al. (2015)¹³ down to our 5σ luminosity limit for an unresolved source (i.e., ${L}_{5\sigma }=6.7\times {10}^{42}$ erg s⁻¹; see Section 4). The derived number density of LAEs is $\phi (L\,\gt {10}^{41.8}$ erg s⁻¹) = (${5.8}_{-1.2}^{+1.9}$) × 10⁻⁴ cMpc⁻³. This means that ≲0.1 LAEs are expected within the total volume explored in our analysis (i.e., 50'' × 50'' × 2000 km s⁻¹ or ∼80 cMpc³; see Section 4) and a mere ∼10⁻³ within a separation of 12.5 kpc (see Figure 15). Such a low incidence of LAEs due to the blank field supports the idea that the detected LAE has to be physically linked to the presence of the QSO.

Zoom In Zoom Out Reset image size

5.3.2. Comparison to QSO–LAE Clustering

The study of the QSO–galaxy and QSO–QSO clustering showed that bright QSOs at low and intermediate redshifts are biased tracers of massive dark matter halos with ${M}_{\mathrm{DM}}\gtrsim {10}^{12.5}\,{M}_{\odot }$ (e.g., Myers et al. 2007; Shen et al. 2007; Padmanabhan et al. 2009; Ross et al. 2009; Farina et al. 2011; Richardson et al. 2012; White et al. 2012; Shen et al. 2013; Zhang et al. 2013; Karhunen et al. 2014; Sandrinelli et al. 2014; Eftekharzadeh et al. 2015, 2017; Garcia-Vergara et al. 2017). At high redshift, constraints are given by the discovery of a QSO pair with a projected separation of only 130 kpc at z = 5.02 (McGreer et al. 2016). The inferred correlation length is r₀ > 29 cMpc, which is consistent with the r₀ ∼ 30 cMpc estimated by Shen et al. (2007, 2010) from a sample of QSO pairs at 3.5 < z < 4.5. To predict the expected number of LAEs in the presence of clustering, we follow the formalism proposed by Hennawi et al. (2006b). In summary, the real space two-point QSO–LAE cross-correlation function ${\xi }^{\mathrm{QG}}$ is integrated along the line of sight to eliminate distance distortions in redshift space, and the expected number of companions is calculated in cylindrical shells (V_shell) centered on the QSO. In practice, ${\xi }_{\mathrm{QG}}$ can be expressed as

$$\begin{eqnarray}&&{\xi }^{\mathrm{QG}}(r,R)={\left[\displaystyle \frac{{({r}^{2}+{R}^{2})}^{\tfrac{1}{2}}}{{r}_{0}^{\mathrm{QG}}}\right]}^{-\gamma },\end{eqnarray} \tag{ 8 }$$

where r and R are comoving coordinates, perpendicular to and along the line of sight, respectively; ${r}_{0}^{\mathrm{QG}}$ is the QSO–LAE cross-correlation length; and γ is the slope of the correlation function (e.g., Peebles 1980). In the presence of clustering, the number of LAEs expected within a cylindrical shell of volume ${V}_{\mathrm{shell}}$(r) = π(r² − ${r}_{\min }^{2}$) ${\rm{\Delta }}R$, where ${r}_{\min }$ is the minimum radius (in comoving Mpc) that can be probed in the MUSE datacube (i.e., the seeing radius), r is the comoving transverse distance from the QSO, and ${\rm{\Delta }}R$ is the comoving line-of-sight distance corresponding to the velocity range ±1000 km s⁻¹ from the QSO's systemic redshift, is

$$\begin{eqnarray}&&{{\rm{N}}}_{{\rm{C}}}({L}_{\mathrm{Lim}},r)=\phi ({L}_{\mathrm{Lim}}){V}_{\mathrm{shell}}(r)[1+{\bar{W}}_{p}^{\mathrm{QG}}({r}_{\min },r)],\end{eqnarray} \tag{ 9 }$$

where the volume-averaged projected cross-correlation function (${\bar{W}}_{p}^{\mathrm{QG}}$) is given by

$$\begin{eqnarray}&&{\bar{W}}_{p}^{\mathrm{QG}}({r}_{\min },r)=\displaystyle \frac{{\int }_{-\tfrac{{\rm{\Delta }}R}{2}}^{+\tfrac{{\rm{\Delta }}R}{2}}{{dR}}^{{\prime} }{\int }_{{r}_{\min }}^{r}2\pi {r}^{{\prime} }{{dr}}^{{\prime} }{\xi }^{\mathrm{QG}}({r}^{{\prime} },{R}^{{\prime} })}{{V}_{\mathrm{shell}}(r)}.\end{eqnarray} \tag{ 10 }$$

To infer ${r}_{0}^{\mathrm{QG}}$, we consider that, by definition, ${\xi }_{\mathrm{QG}}$ could be expressed as

$$\begin{eqnarray}&&{\xi }^{\mathrm{QG}}(r,R)=\langle {\delta }_{{\rm{Q}}}(r,R){\delta }_{{\rm{G}}}(r,R)\rangle ,\end{eqnarray} \tag{ 11 }$$

where ${\delta }_{{\rm{Q}}}$ and ${\delta }_{{\rm{G}}}$ are the relative density contrasts of QSOs and LAEs, respectively, and the angular brackets denote averaging over a distribution. Under the assumption that QSOs and LAEs trace the same underlying dark matter distribution and considering a linear bias (${\delta }_{{\rm{Q}}}={b}_{{\rm{Q}}}{\delta }_{\mathrm{DM}}$ and ${\delta }_{{\rm{G}}}={b}_{{\rm{G}}}{\delta }_{\mathrm{DM}}$), Equation (11) can be rewritten as

$$\begin{eqnarray}{\xi }^{\mathrm{QG}}(s,R) & = & {b}_{{\rm{Q}}}{b}_{{\rm{G}}}{\xi }^{\mathrm{DMDM}}(r,R)\\ & = & {[{\xi }^{\mathrm{QQ}}(r,R){\xi }^{\mathrm{GG}}(r,R)]}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 12 }$$

where ${\xi }^{\mathrm{DMDM}}$, ${\xi }^{\mathrm{QQ}}$, and ${\xi }^{\mathrm{GG}}$ are the real space two-point auto-correlation functions for dark matter halos, QSOs, and LAEs, respectively (e.g., Schneider 2006). Assuming a power-law form for both ${\xi }^{\mathrm{QQ}}$ and ${\xi }^{\mathrm{GG}}$ (with the same γ), the QSO–LAE cross-correlation length could be estimated as

$$\begin{eqnarray}&&{r}_{0}^{\mathrm{QG}}={({r}_{0}^{\mathrm{QQ}}{r}_{0}^{\mathrm{GG}})}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 13 }$$

where ${r}_{0}^{\mathrm{QQ}}$ and ${r}_{0}^{\mathrm{GG}}$ are the QSO–QSO and LAE–LAE auto-correlation lengths. As a proxy for ${r}_{0}^{\mathrm{QG}}$ at z ∼ 6.6, we considered ${r}_{0}^{\mathrm{GG}}={10.3}_{-8.6}^{+4.7}$ cMpc, which is the upper limit on the auto-correlation length derived by Ouchi et al. (2010) from the clustering analysis of bright LAEs at z = 6.6 and ${r}_{0}^{\mathrm{QQ}}={17.4}_{-2.8}^{+2.5}$ cMpc obtained by imposing γ = 1.8 on the study of the clustering properties of z > 2.9 QSOs of Shen et al. (2007). It is worth noticing that, at z ∼ 1–2, the QSO–QSO auto-correlation function appears to get steeper at sub-Mpc separations (e.g., Djorgovski 1991; Hennawi et al. 2006; Djorgovski et al. 2007; Myers et al. 2007, 2008; Hennawi et al. 2010; Kayo & Oguri 2012; Farina et al. 2013). However, this enhancement becomes less prominent at z > 2.9 where small- and large-scale clustering amplitudes are comparable (Shen et al. 2010).

Figure 15 shows, as a function of r, the estimated number of LAEs located within ±1000 km s⁻¹ from the QSO's systemic redshift. Considering the cross-correlation length estimated above, ≲0.1 LAEs are expected within 100 ckpc from the QSO. The presence of one LAE with such a small separation suggests a physical interaction with the QSO's host galaxy. Intriguingly, the detected LAE and the possible Lyα halo point to each other both in velocity and in position with respect to the QSO, opening up the possibility that we may be witnessing an ongoing merger and that the Lyα halo may be associated with the interaction between the LAE and the QSO's host galaxy. We are cautious, however, that these results are based on the detection of a single LAE in the proximity of one z ∼ 6.6 QSO. Previous studies (Decarli et al. 2012; Bañados et al. 2013; Mazzucchelli et al. 2017b; Goto et al. 2017), even if subjected to larger redshift uncertainties, did not report any evidence for an excess of LAEs in the proximity of z ∼ 6 QSOs (but see Decarli et al. 2017 for the detection of [${\rm{C}}\,{\rm{II}}$] bright galaxies in the proximity of z ∼ 6 QSOs). Measurements of the high-redshift QSO–galaxy clustering in a statistical fashion is thus fundamental to discern among different scenarios proposed to explain the rapid formation of SMBHs at the end of cosmic reionization.

5.3.3. Physical Properties of the LAE

The Lyα emission of the companion galaxy could be boosted by the local enhancement of the ionizing background in the vicinity of a QSO. Star formation can power the Lyα emission up to a rest-frame equivalent width of W₀(Lyα) = 240 Å (e.g., Schaerer 2002). This value is commonly used as a limit to identify LAEs associated with a fluorescent reprocessing of QSO radiation (e.g., Cantalupo et al. 2012). Some authors consider a less stringent limit of W₀(Lyα) > 100 Å for their selection (e.g., Trainor & Steidel 2013; Borisova et al. 2016b). However, this may lead to a high level of contamination from non-fluorescent objects (e.g., Borisova et al. 2016b).

Here, we test if the properties of the identified LAEs are consistent with a Lyα fluorescence scenario. No significant continuum emission is detected redward of the Lyα emission. To constrain W₀(Lyα), we thus consider the 1σ limit on the mean continuum (${\sigma }_{{\rm{C}}}$) obtained by averaging down the errors on the extracted spectrum (${\sigma }_{\mathrm{LAE},\lambda }$) over the wavelength range between 9281 Å and 9350 Å (i.e., up to the edge of the datacube): ${\sigma }_{{\rm{C}}}^{2}={\sum }_{\lambda =9281\mathring{{\rm{A}}}}^{9350\mathring{{\rm{A}}}}1/{\sigma }_{\mathrm{LAE},\lambda }^{2}$. This leads to a 3σ lower limit on the rest-frame equivalent width of ${W}_{3\sigma }(\mathrm{Ly}\alpha )=\tfrac{F(\mathrm{Ly}\alpha )}{3\,{\sigma }_{{\rm{C}}}(1+{z}_{\mathrm{LAE}})}=7.4$ Å, which is not strict enough to rule out a non-fluorescent scenario. As a matter of fact, if powered by fluorescence, the luminosity of the LAE is proportional to ${L}_{{\nu }_{\mathrm{LL}}}$ and inversely proportional to the square of the perpendicular distance from the QSO, assuming, for the sake of simplicity, that the line-of-sight separation is negligible (see Hennawi & Prochaska 2013). We would expect $L(\mathrm{Ly}\alpha )$ ∼ 7 × 10⁴³ erg s⁻¹, where we considered the LAE to be unresolved, i.e., with a diameter matching the seeing FWHM (058). The observed luminosity is a factor ∼30× fainter. Note that this estimate depends on unknown quantities, namely, the real distance of the QSO, the luminosity of the QSO in the direction of the LAE, and the size of the surface illuminated by the ionizing radiation. These could contribute to reduce the expected fluorescence emission. On the other hand, it is possible that we are underestimating the real Lyα luminosity. Up to 90% of the flux may come from an undetected extended component (Wisotzki et al. 2016) and a fraction of the Lyα emission could be concealed by the neutral hydrogen. Within these uncertainties, we favor a scenario where star formation, rather than fluorescence, is inducing the Lyα emission. From Equation (4), we can thus derive the star-formation rate of the LAE, SFR_LAE ∼ 1.3 M_⊙ yr⁻¹ (see discussion in Rauch et al. 2013 for the possibility of QSO feedback triggering the star formation in a close-by LAE at z ∼ 3.0).

In this paper, we presented a sensitive search for extended Lyα emission around the starbursting QSO J0305−3150 at z ∼ 6.61. The nominal 5σ SB limit reached with MUSE is ${\mathrm{SB}}_{5\sigma }^{1}\,=1.9\times {10}^{-18}$ erg s⁻¹ cm⁻² arcsec⁻² over a 1 arcsec² aperture (estimated by collapsing the datacube over five wavelength slices centered at λ = 9256 Å). This formally corresponds to a luminosity limit of ${L}_{5\sigma }\sim {10}^{42.0}$ erg s⁻¹ for unresolved sources with FWHM = 200 km s⁻¹, and to ${L}_{5\sigma }\sim {10}^{43.2}$ erg s⁻¹ for a circular source with a diameter of 50 kpc and the same FWHM. The primary results of this study are:

• 1.  
  After carefully subtracting the unresolved emission from the central QSO, we detect the presence of a tenuous ($L(\mathrm{Ly}\alpha )$ = (3.0 ± 0.4) × 10⁴² erg s⁻¹) Lyα halo extending over ∼10 kpc. To date, this is the first such nebula observed at z > 6.5 and one of the faintest extended emissions ever observed around a QSO at any redshift (see Figure 14). Despite the depth of our data, we do not detect the large-scale (10–100 kpc) Lyα emission frequently observed in z ∼ 2–4 bright QSOs.
• 2.  
  A comparison between the Lyα emission revealed by MUSE and the FIR properties of the host galaxy inferred from ALMA observations allows us to speculate on the geometry of the dust surrounding the QSO. In particular, a patchy geometry of the dust cocoon can explain (i) the displacement between extended Lyα and resolved [${\rm{C}}\,{\rm{II}}$] emission lines, (ii) the discrepancy between the observed Lyα luminosity and expectation from the intense star formation of the host galaxy, and (iii) the mild dust extinction present along the QSO line of sight. This configuration permits the ionizing radiation from the newly formed stars to escape the host galaxy and give rise to the observed Lyα emission.
• 3.  
  We estimate that the extended Lyα emission is too faint to arise from recombinations on the "skin" of optically thick clouds. A more plausible scenario is that the QSO radiation is sufficiently intense to maintain the surrounding gas highly ionized, hence we are observing fluorescent emission coming from optically thin clouds. Intriguingly, a consequence of this emission mechanism is a hydrogen volume density of the gas illuminated by the QSO of ${n}_{{\rm{H}}}\sim 6$ cm⁻³, similar to that estimated for the giant (≫100 kpc) Lyα nebulae recently discovered around z ∼ 2 QSOs.
• 4.  
  An LAE with $L(\mathrm{Ly}\alpha )$ ∼ 2.1 × 10⁴² erg s⁻¹ is present at 12.5 kpc and 560 km s⁻¹ from J0305−3150. Our current constraints on the rest-frame equivalent width (W₀(Lyα) > 7.4 Å) and luminosity, although not conclusive, disfavor a picture where the Lyα emission is boosted by the QSO radiation. Assuming that the Lyα line is powered solely by star formation, we derive a star-formation rate of SFR_LAE ∼ 1.3 M_⊙ yr⁻¹.
• 5.  
  We calculate the enhanced probability of finding such a close LAE due to the clustering of galaxies around QSOs. From the extrapolation of the z ∼ 6.6 QSO–LAE large-scale correlation function, we estimate this probability to be small (<10%). This supports a picture in which dissipative interaction and QSO activity in the young universe are connected.

E.P.F., B.P.V., and F.W. acknowledge funding through the ERC grant "Cosmic Dawn." Support for R.D. was provided by the DFG priority program 1573 "The physics of the interstellar medium." S.C. gratefully acknowledges support from Swiss National Science Foundation grant PP00P2_163824. E.P.F. is grateful to M. Fouesneau for providing support in the use of Python for the analysis of the MUSE datacubes, to M. Fumagalli for sharing data on the QSO Q0956+122, and to M. Rauch for providing feedback on the manuscript. We thank the members of the ENIGMA group¹⁴ at the Max Planck Institute for Astronomy (MPIA) for helpful discussions. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013); APLpy,¹⁵ an open-source plotting package for Python based on Matplotlib (Hunter 2007); and IRAF.¹⁶ Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 094.B-0893(A).