UBIQUITOUS GIANT Lyα NEBULAE AROUND THE BRIGHTEST QUASARS AT z ∼ 3.5 REVEALED WITH MUSE^∗

Our total sample of 17 radio-quiet quasars, complemented by two radio-loud systems (see Table 1), is composed of two subsamples reflecting two different MUSE GTO programs (094.A-0396, 095.A-0708, 096.A-0345 PI: S. Lilly; 094.A-0131, 095.A-0200, 096.A-0222 PI: J. Schaye). The main subsample of 12 radio-quiet quasars has been constructed specifically for this study (094.A-0396, 095.A-0708, 096.A-0345) using the catalog of VeronCetty & Veron (2010) and selecting the brightest radio-quiet quasars known in the redshift range of to maximize MUSE throughput and minimize redshift dimming of the surface brightness (SB). We have removed quasars in proximity of bright stars and fields with high galactic extinction. Among the remaining quasars, we have then selected objects with available UVES spectroscopy to maximize "transverse" science projects within the MUSE GTO program.

Table 1.  Quasar Sample and Observation Log

Number	Quasar	R.A.	decl.	zcata	zsysb	c	d	e	Classf	Seeingg	Skyh
		(J2000)	(J2000)			(erg s−1)	(mag)	(mag)		(arcsec)	Conditions
1	CTS G18.01	00:41:31.4	−49:36:11.9	3.240	3.207		16.621	−30.22	RQ	1.08	CL
2	Q0041−2638	00:43:42.7	−26:22:10.9	3.053	3.036		18.328	−28.42	RQ	1.14	WI-CL/TN
3	Q0042−2627	00:44:33.5	−26:11:25.9	3.289	3.280		18.663	−28.23	RQ	1.18	WI-PH
4	Q0055−269	00:57:58.1	−26:43:15.8	3.662	3.634		17.475	−29.52	RQ	1.02	PH/CL
5	UM669	01:05:16.7	−18:46:41.9	3.037	3.021		17.811	−28.92	RQ	1.31	CL/PH
6	J0124+0044	01:24:04.0	00:44:33.5	3.810	3.783		18.099	−28.98	RQ	0.82	PH
7	UM678	02:51:40.4	−22:00:28.3	3.205	3.188		18.014	−28.81	RQ	0.72	PH
8	CTS B27.07	04:45:33.1	−40:48:42.8	3.270	3.132		17.978	−28.82	RQ	0.59	PH
9	CTS A31.05	05:17:42.1	−37:54:45.9	3.020	3.020		17.824	−28.91	RQ	0.72	CL
10	CT 656	06:00:08.7	−50:40:30.1	3.130	3.125		17.549	−29.24	RQ	0.70	CL
11	AWL 11	06:43:26.9	−50:41:12.9	3.090	3.079		18.078	−28.69	RQ	0.63	PH
12	HE0940−1050	09:42:53.6	−11:04:26.0	3.093	3.050		16.630	−30.12	RQ	0.74	CL
13	BRI1108−07	11:11:13.7	−08:04:03.0	3.910	3.907		18.312	−28.86	RQ	0.98	TK/TN
14	CTS R07.04	11:13:50.1	−15:33:40.2	3.370	3.351		17.601	−29.36	RQ	0.94	TN
15	Q1317−0507	13:20:29.8	−05:23:34.2	3.700	3.701		17.525	−29.51	RQ	0.94	CL
16	Q1621−0042	16:21:16.7	−00:42:48.2	3.700	3.689		17.079	−29.95	RQ	0.85	TK
17	CTS A11.09	22:53:10.7	−36:58:15.9	3.200	3.121		17.815	−28.98	RQ	0.76	CL
R1	PKS1937−101	19:39:57.4	−10:02:39.9	3.787	3.769		16.727	−30.35	RL	0.75	CL
R2	QB2000−330	20:03:24.1	−32:51:45.9	3.783	3.759		17.302	−29.77	RL	0.96	CL

Notes.

^aTaken from the catalog VeronCetty & Veron (2010). ^bMeasured from MUSE spectra from the peak of the quasar C iv emission and correcting for luminosity-dependent velocity shifts using Shen et al. (2016). The level of the intrinsic uncertainty of the C iv correction relative to the systemic redshift is ∼415 km s⁻¹ (). ^cSpecific monochromatic continuum luminosity in the observed frame used to compute the correction for the estimated C iv redshift in Shen et al. (2016). ^dComputed from MUSE datacubes with circular aperture photometry using a radius of 3 arcsec and assuming an SDSS i-band filter. No correction for galactic absorption has been applied. ^eAbsolute i-band magnitude normalized at z = 2 using Ross et al. (2013). ^fRQ—radio-quiet quasar; RL—radio-loud quasar. Classification is based on radio flux measurements from Carilli et al. (2001) and Condon et al. (1998) for RQ and RL respectively. ^gSeeing measured in the combined 1 hr datacubes as the FWHM of a Gaussian profile. ^hSky conditions during the night of observations. The meaning of the labels is the following: PH-photometric night; CL-Clear night; WI-Strong winds; TN-thin clouds; TK-thick clouds.

Download table as:  ASCIITypeset image

The second subsample of five radio-quiet and two radio-loud quasars is part of a MUSE GTO program to study the connection between absorption line systems in quasar spectra and emission-line galaxies (094.A-0131, 095.A-0200, 096.A-0222). The selection criteria of this program required quasars at higher redshift, i.e., , in order to maximize the available path-length of the forest within the MUSE wavelength range. Moreover, these quasars have very deep UVES spectroscopy. Apart from the different redshift, the five radio-quiet quasars of this subsample are very similar to the lower redshift objects in terms of luminosity. Although not originally part of this study, we have also included two radio-loud quasars from this subsample out of the three observed so far. In particular, we have excluded one radio-loud quasar (B1422+2309) from our sample because it is gravitationally lensed (Patnaik et al. 1992) into multiple components making PSF subtraction challenging and because the field is very crowded with foreground galaxies.

The three giant nebulae discovered hitherto with NB imaging (Cantalupo et al. 2014; Martin et al. 2014; Hennawi et al. 2015) are all characterized by bright extended emission with surface brightness values larger than 10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻², a value that is easily reachable within 1 hr of integration time with MUSE. Therefore, for this exploratory survey (094.A-0396, 095.A-0708, 096.A-0345), we use a total exposure time of 1 hr for each quasar split into 4 × 900 s exposures. Between each individual exposure, we rotate the FoV by 90° and apply a small random dithering pattern of less than 1 arcsec. Before starting each observation, we offset the quasar position by a few arcsec away from the center of the FoV to avoid regions with higher systematics. A similar strategy has been applied to every quasar field observed as part of the MUSE GTO. Some of the quasars observed as a part of programs 094.A-0131, 095.A-0200, and 096.A-0222 have a total integration time longer than 1 hr. In order to keep our sample homogeneous in depth, we have selected only the first 4 × 900 s exposures for these quasars.

Data was collected with the MUSE/VLT instrument (Bacon et al. 2010) between 19 September 2014 and 9 November 2015. All but two observation blocks (OBs) were executed continuously during the same nights. The seeing varied in the range 0.59–1.31 arcsec (FWHM of the Gaussian at 7000 Å , measured in the combined 1 hr datacubes). The information about the quasar fields is summarized in Table 1. Together with name and coordinates, we provide (1) the original redshift from the catalog, (2) our systemic redshift estimate based on the C iv emission corrected for blueshift according to the quasar luminosity as in Shen et al. (2016), (3) the luminosity at the 1700 Å rest-frame, (4) the i-magnitude and corresponding absolute magnitude normalized to z = 2 as in Ross et al. (2013), and (5) the radio class, seeing, and sky conditions during each night.

We used the standard MUSE pipeline v1.0 (Weilbacher et al. 2012, 2014; P. M. Weilbacher 2016, in preparation) for the basic steps of the data reduction with the default (recommended) parameters. For each of the individual exposures, we performed bias subtraction, flat-fielding, twilight, and illumination correction, and wavelength calibration. Sky subtraction was not performed with the pipeline but was done at a later stage with custom developed software as discussed below. The response curve and telluric correction were obtained from one of the spectrophotometric standards observed during the same night except for field #17, where we had to use a standard star from the previous night. Finally, flux calibrated data was drizzled onto a 3D grid using the information from geometry and astrometry tables in order to produce the final datacubes.

It is known that MUSE cubes reduced with the standard pipeline might have small astrometric offsets in the coordinate system because of a small "derotator wobble" (Bacon et al. 2015). To correct for this effect, we registered the datacubes using the position of point sources in MUSE white-light images (obtained by collapsing the datacubes along the wavelength direction) in different exposures.

The final steps of the data reduction were performed with custom tools for flat-fielding correction and sky-subtraction that are part of the CubExtractor package (S. Cantalupo 2016, in preparation) and that have been specifically developed to improve data quality for the detection of faint and diffuse emission in MUSE datacubes. In particular, the flat-fielding correction is performed as a self-calibration on each individual datacube using the sky-continuum and the sky-lines as a spatially uniform source to re-calibrate each individual slice (part of an Integral Field Unit, IFU) and IFU as a function of wavelength. Sources are masked with an iterative procedure to ensure that self-calibration errors are minimized. With this procedure, the typical striped pattern of broadband and NB images of MUSE cubes is totally removed and any visible residual is typically at a level much less than 0.1% of the sky (a more detailed description of this procedure, called CubeFix, will be presented in S. Cantalupo 2016, in preparation). In some rare cases, there are not enough spatial and spectral elements in a slice to find a suitable correction or there is clear variation in a single slice that cannot be corrected with a simple rescaling factor. In these cases, we mask the volume pixels (voxels) in the datacube corresponding to the slices or region without a proper correction factor. Because of dithering and FoV rotation, these masked regions do not typically affect the final combined datacube.

Sky-subtraction is then performed on each individual, flat-field-corrected cube using CubeSharp (S. Cantalupo 2016, in preparation). CubeSharp uses a local and flux-conserving procedure to empirically correct the sky line spread function (LSF) and therefore remove sky lines minimizing the residuals due to LSF shifts and variation across the MUSE FoV—the major sources of systematic errors in MUSE cubes. Because the algorithm is flux-conserving by design, no residuals are introduced when sky-subtraction is performed with CubeSharp.

Finally, the corrected and sky-subtracted cubes are combined using an average -clipping algorithm. After this first iteration, a white-light image is created and continuum sources are identified using CubExtractor (S. Cantalupo 2016, in preparation; see Section 3.3). Using the positions and spectra of continuum sources from the combined cube, another iteration of CubeFix and CubeSharp is performed on individual exposures to improve the removal of self-calibration effects. Typically one iteration is sufficient to substantially improve the data reduction process before the individual cubes are combined again.

To reveal the presence of extended in the vicinity of a quasar, we need to remove the contribution from the quasar PSF. The advantage of integral-field spectroscopy (IFS) is that we have images of the quasar and its surroundings at different wavelengths, including regions in the spectra for which we are sure that no extended line emission should be present.

We have explored two different approaches of quasar PSF estimation and removal. One way is to fit a known function, for example, a Gaussian or Moffat profile. Unfortunately, the MUSE PSF in reconstructed and combined cubes is complex due to the nature of the instrument and the wings of the PSF profiles cannot be easily modeled by these simple analytic profiles. Indeed, both Gaussian and Moffat fitting produced clear residuals of the quasar PSF in subtracted images at a level that would preclude the detection of faint and extended emission unrelated to the quasar PSF (see also Christensen et al. 2006 and references therein).

The other method is purely empirical and uses the data itself to construct PSF images that are rescaled and subtracted at each wavelength layer¹² (see also Husemann et al. 2013 and Herenz et al. 2015 for an other iterative empirical approach). In particular, for each wavelength layer, we produce a pseudo-NB image with a spectral width of 150 spectral pixels () at the position of the quasar. This spectral width is a compromise between minimizing the PSF wavelength variations and maximizing the signal-to-noise ratio (S/N) in the empirical PSF images. The flux in each empirical PSF image is then rescaled assuming that the quasar is dominating the flux in the majority of the central 5 × 5 pixels, i.e., 1 × 1 arcsec². In particular, because the flux in individual pixels may be affected by cosmic rays or other artifacts, we compute the rescaling factor between the flux in each layer and the empirical PSF images using an averaged-sigma-clip algorithm. Once the empirical PSF image has been rescaled, we cut from it a central circular region with a radius of about five times the seeing (masking any pixel with negative flux) and subtract this circular cut-out from the corresponding wavelength layer in the datacube. In some cases, the nebulae are so bright that their central parts will be visible at a low level even in a -wide band. In this case, we iterate the PSF removal procedure, masking the wavelength region associated with the nebulae to avoid over-subtraction.

This empirical PSF subtraction produces excellent results, especially on large scales around the quasar. However, we note that this method cannot provide any meaningful information on the central 1 arcsec region used for the PSF rescaling and that it is not able to treat blending of the quasar PSF with other nearby continuum sources, if present. For the purposes of this study, we will see that the method is adequate since the detected nebulae extend far beyond the central 1 arcsec region around the quasars. Moreover, we mask and remove continuum sources around the quasar as described below, so their residuals do not affect the results presented here.

After quasar PSF subtraction has been performed, we remove any other possible continuum sources for each spaxel in the cube using a fast median-filtering approach based on the following method. (1) We first resample the spectrum of each spaxel into spectral regions (bins) with a size of 150 spectral pixels () using a median filter. (2) The resampled spectrum is then smoothed with a median filter with radius equivalent to two bins to minimize the effect of line features in individual spectral regions. (3) We subtract from each voxel in the cube the estimated continuum from the corresponding spaxel and spectral region. This procedure allows a fast and efficient removal of continuum sources (a full median-filter approach for each spaxel and each spectral pixel would be more computationally expensive). In some cases, stars or background galaxies with emission lines in their spectra are not properly removed by this procedure due to the large window size and negative or positive residuals are visible in the cube. However, these residuals do not affect our results because we mask any bright star or galaxy in the cube before extraction, as described below.

The final step in our data analysis procedure is the detection and extraction of extended line emission from the reduced and processed cubes. For this task, we use newly developed three-dimensional automatic extraction software called CubExtractor (S. Cantalupo 2016, in preparation) based on a three-dimensional extension of the connected-labelling-component algorithm with union finding of classical binary image analysis (see, e.g., Shapiro & Stockman 2001). In particular, after datacubes and associated variances are being filtered (smoothed), voxels above a given user-defined S/N threshold are connected and objects are detected if they contain a certain number of connected voxels above a user-defined threshold. Three-dimensional segmentation masks produced by CubExtractor are then used for photometry and in order to obtain several lower-dimensional projections of the extracted objects, such as optimally extracted images, two-dimensional spectra, and velocity maps, as we will show in the next sections.

Because detection and extraction is based on an S/N threshold, the characterization of the noise in MUSE datacubes is a crucial aspect of the process. A propagated variance, estimated for each individual step of the reduction process, is provided by the MUSE pipeline and propagated during flat-fielding and continuum subtraction with CubeFix and CubeSharp. The variance is also propagated during exposure combination. The variance associated with MUSE resampled datacubes by the pipeline is an underestimate of the true variance because resampling introduces correlated noise that cannot be easily captured by current detection and extraction algorithms. In order to include this extra source of noise in an approximate fashion and to obtain the right "normalization" for the variance, we proceed in the following way. (1) We estimate from the cube itself the "spatial" variance of the pixel flux for each wavelength layer (this estimate is essentially a single value for each layer and therefore does not contain information on the spatial variation of the noise). (2) We compute for each wavelength layer the spatial average of the propagated variance obtained by the pipeline. (3) We rescale the propagated variance by a constant factor for each wavelength layer in order to match the average "spatial" variance estimated from the cube itself. The value of the rescaling factor is typically around 1.4.

Using an initial guess for the redshift of the quasar (from the original catalog or using the location of C iv/1550 and He ii/1640 lines), we extract various subcubes from the processed datacube with wavelength ranges capturing the expected line, C iv/1550 and He ii/1640 lines (the sizes of the subcubes correspond to ∼15,000 km s⁻¹). Because we are only interested in extended emission, we use a large number of minimum connected voxels in CubExtractor for detection (10,000). Moreover, we apply before detection a spatial Gaussian filtering with the σ value of 0.5 arcsec (without smoothing in wavelength) to bring out extended but narrow features, and we use a minimum S/N of two with respect to the rescaled variance cubes as described above. This value typically corresponds to a surface brightness of about 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² for a 1 arcsec² aperture in a single wavelength layer (i.e., ) for our cubes.

With this set of parameters and masking regions associated with bright continuum sources or sky residuals, we always find one single detection in the expected Lyα wavelength range for every quasar field. As we will show in the next section, these detections are giant Lyα nebulae extending up to several tens of arcsec around the quasars.

As discussed in Section 3, we have detected an extended source in each individual quasar field around the expected wavelength for emission. Each of these sources is characterized by more than 10,000 connected voxels (or individual spatial and spectral elements in MUSE cubes) with an S/N greater than two—after smoothing—and the confidence level of the detections, as we will show in this section, is very high.

In Figure 1, we present the "optimally extracted images" of the detected objects in each MUSE cube. Each image has a linear size of 44 arcsec and the original position of the quasar is marked by a black dot. These images have been obtained using the three-dimensional segmentation mask (called the 3D-mask here for simplicity) associated with the detection object. We note that the 3D-mask defines a three-dimensional iso-S/N surface in the cube (after spatial smoothing as described in Section 3) and therefore is ideal to obtain images and spectra with maximal S/N after collapsing one of the spatial or the spectral dimensions. In particular, the images presented in Figure 1 have been obtained by selecting all voxels in the PSF-subtracted and continuum-subtracted MUSE cubes, using corresponding 3D-masks of each nebula, and integrating their fluxes along the wavelength direction. These images can also be thought of as pseudo-NB images where the width of the filter is adjusted for each spaxel to maximize the S/N of the objects. The widths of this pseudo-filters vary from one layer (typically at the edges of the object) to a few tens of layers in the brightest or kinematically broader parts of the sources. The projection of the 3D-mask on the plane of the sky is indicated by the thick contours overlaid on the images in Figure 1 and typically corresponds to an SB of about 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻². We note that this deep sensitivity level—obtained in a single hour of MUSE observation—is comparable to the sensitivity of a 20h observation with a NB filter on an 8-meter telescope (e.g., Cantalupo et al. 2012). For display purposes, we have added—by means of the union operator—to the 3D-mask one wavelength layer of the cube corresponding to the central wavelength of the nebulae. The voxels associated with this layer are shown outside of the thick contours in Figure 1.

Zoom In Zoom Out Reset image size

We stress that these "optimally extracted images" are quite different from standard NB (see Appendix A) or broadband images—with which the reader may be more familiar—because the number of wavelength layers (and therefore the flux and corresponding noise) contributing to the image depends on spatial position and is thus correlated with the S/N. This gives us the unique possibility of capturing in a single image, at the best S/N threshold, the surface brightness values for both kinematically narrow and broad features that would have been either lost in the noise or underestimated in an NB image with a single width. As a consequence, the images in Figure 1 have a much larger dynamic range, despite the short integration time, with respect to a standard image. One possible drawback of this approach, however, is that the image noise cannot be simply estimated visually as the noise will be correlated with the number of layers and therefore with the flux and spatial position. To help the reader visualize the true noise in these images, we have therefore estimated the noise and the S/N for each pixel in the image by variance propagation (using the same rescaling factor derived in the CubExtractor run) taking into account the number of layers contributing to each pixel. The thin contours in Figure 1 show the image S/N contours using these propagated variances and the integrated flux. The first contour corresponds to an S/N of two, the second one to four, and further steps between the contours to ΔS/N = 6. These contours are the closest representation of the true image noise as they are, ideally, independent of the number of layers contributing to each pixel (of course, this would not be true in the case of significant correlated noise in the wavelength direction).

In Figure 2, we show different projections of the PSF and continuum-subtracted MUSE datacubes using the CubExtractor 3D-masks, i.e., the "optimally extracted" two-dimensional and one-dimensional spectra. In each case, we have used the spatial projection of the 3D-masks as an aperture, i.e., the thick contours shown in Figure 1, to extract two-dimensional spectra by integrating the associated voxels along the east direction (the spatial x-axis in Figure 1). The x-axes in Figure 2 represent the wavelength dimension and the y-axes correspond to the spatial y-axes in Figure 1 (with the same scale as in Figure 1). Below each two-dimensional spectrum, we also show the associated one-dimensional spectrum obtained by integrating all spatial pixels (black lines) and compare it to a scaled version of the original quasar spectrum (gray lines).

Zoom In Zoom Out Reset image size

Because we are using a two-dimensional projection of the 3D-mask, the number of integrated voxels, and therefore the noise level in the two-dimensional spectra, is constant across the wavelength direction at each fixed spatial position. The emergence of the extended line emission in this figure is clear, particularly in the integrated one-dimensional spectra. This figure also confirms that the detected emission for the smallest and more circular nebulae, e.g., #11 and #17, cannot be due to PSF removal artifacts, as the nebular spectrum (thick line) is always narrower than and strikingly different from the original quasar spectrum (thin lines).

As is clear from Figure 1, the nebulae have a large diversity of shapes and sizes. The smallest nebulae tend to be more symmetric and circular (with the exception of #2) while the largest (#1, #3, extending up to 320 pkpc) show evidence of filamentary structures and multiple components. The radio-loud quasar nebulae (R1 and R2) do not look particularly different from other nebulae in these images. The brightest nebula (#6) has a peculiar morphological structure with a sudden, steep decrease of the total flux around a distance of about 50 pkpc from the quasar position. In some cases (e.g., #7) the brightest part of the nebulae is clearly offset with respect to the quasar position and in general there seems to be a small offset between the nebula flux centroid and the quasar position. There is no evidence for "bipolar ionization cone illumination" patterns, although nebulae #3 and #13 show hints of "bipolar" structure in their morphology.

We have defined the sizes of the nebulae as the maximum projected sizes of the 3D-masks (see Figure 1) and have compared the sizes to the total luminosities as is commonly done in recent literature, see Table 2 and Figure 3. For illustrative purposes, we include in the same figure the results of other studies targeting radio galaxies (Reuland et al. 2003; Villar-Martín et al. 2007), radio-loud quasars (Heckman et al. 1991b), Lyα blobs (Matsuda et al. 2004; Prescott et al. 2009), and other radio-quiet quasars (Bergeron et al. 1999; Christensen et al. 2006; North et al. 2012). We stress that this comparison is only qualitative because these studies have different sensitivity limits, redshifts, methods, and definitions of sizes.

Zoom In Zoom Out Reset image size

Table 2.  Properties of the Detected Giant Nebulae

Number	Object Name	a	b	Sizec	d	Luminosity	FWHMe	C iv/ f	He ii/ f
				(pkpc)	()	(erg s−1)	(km s−1)	()	()
1	CTS G18.01	3.207	3.248	240	47.50		570	<0.04	<0.03
2	Q0041−2638	3.036	3.078	170	23.75		320	<0.06	<0.02
3	Q0042−2627	3.280	3.306	320	47.50		510	<0.01	<0.01
4	Q0055−269	3.634	3.662	180	55.00		770	0.07	<0.03
5	UM669	3.021	3.040	160	40.00		660	<0.05	<0.03
6	J0124+0044	3.783	3.847	190	46.25		930	0.04	<0.03
7	UM678	3.188	3.208	150	31.25		580	<0.06	<0.11
8	CTS B27.07	3.132	3.155	160	35.00		590	<0.04	<0.06
9	CTS A31.05	3.020	3.050	120	41.25		780	<0.09	<0.04
10	CT 656	3.125	3.159	130	31.25		640	<0.07	<0.06
11	AWL 11	3.079	3.118	130	30.00		670	<0.07	<0.04
12	HE0940−1050	3.050	3.091	170	45.00		660	<0.06	<0.01
13	BRI1108−07	3.907	3.935	160	53.75		760	<0.04	<0.04
14	CTS R07.04	3.351	3.368	170	40.00		660	0.05	0.01
15	Q1317−0507	3.701	3.720	140	33.75		560	<0.11	<0.09
16	Q1621−0042	3.689	3.704	120	31.25		550	<0.07	<0.05
17	CTS A11.09	3.121	3.150	150	28.75		490	<0.10	<0.05
R1	PKS1937−101	3.769	3.791	110	103.75		2120	0.11	0.12
R2	QB2000−330	3.759	3.788	120	52.50		1120	0.06	<0.06

Notes.

^aSystemic redshift of the quasar as in Table 1 derived from the peak of C iv line and corrected according to Shen et al. (2016). ^bMeasured from the flux-weighted centroid of the nebular emission. ^cMaximum linear projected size measured from the spatial projection of the CubExtractor object segmentation mask ("3D-mask"). ^dMaximum spectral width of the "3D-mask" and width of the pseudo NB used for the circularly averaged SB profile. ^eSpatially averaged FWHM. ^fLimits on the line ratios correspond to 2σ.

Download table as:  ASCIITypeset image

In this section, we investigate the average radial SB profile of the nebulae as a way of constraining the origin of the emission. Despite the different morphologies and clear asymmetries in some of the nebulae, we decided to use the standard approach in the literature of measuring circularly averaged SB profiles for both simplicity and ease of comparison with previous works. For similar reasons, we decided to use standard NB images derived from MUSE cubes, i.e., using a fixed width rather than the 3D-mask discussed above. In particular, we fix the width of these pseudo-NB images to the maximum spectral width of the nebulae as defined by the 3D-mask (see Table 2). We note that the resulting profiles are noisier, especially at the edges of the nebulae, with respect to the images presented in Figure 1. However, the noise is better characterized and there are no effects due to S/N thresholding of the flux. We also computed and subtracted, if present, any residual background level from each pseudo NB.¹³ This residual background level was not higher than 1σ of the flux pixel distribution in any of the case. Therefore, its contribution is small and only affects the profiles at larger distances from the quasars.

The resulting circularly averaged SB profiles for each of the nebulae are presented as black lines in Figure 4, while the purple line is the average profile that is well fitted by a power law with slope of −1.8. The gray line is a 2σ Gaussian estimate of the noise associated with the profile. Surprisingly, despite the different morphologies, sizes and luminosities, all profiles look very similar to each other, including those of the radio-loud quasars, and are better represented by power laws in the majority of the cases rather than exponential profiles (in Appendix B we report the results of our fitting procedure including both power law and exponential profiles). The most notable exception is #6, which also has a peculiar morphology, as discussed in the previous Section.

Zoom In Zoom Out Reset image size

Figure 5 compares these profiles (black lines for individual profiles, purple line for the average profile, as in Figure 4) with other studies. In particular, we show as a orange line the profile of the Slug Nebula (Cantalupo et al. 2014). We note that this profile was not PSF-subtracted¹⁴ and therefore we use as a dashed line style to indicate in the Slug Nebula the region of the profile that is contaminated by the quasar PSF. The blue line shows the average profile of halos around emitting galaxies in the deep MUSE observation of the HDFS (Wisotzki et al. 2016) with the shaded blue area indicating the range of individual profiles. Finally, the green line is the stacked SB profile of NB images at wavelengths of 92 Lyman Break Galaxies (LBGs) at obtained by Steidel et al. (2011).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Clearly, the profiles of the MUSE nebulae, presented in the leftmost panel, are very similar to the Slug Nebula despite the different luminosity and redshift and, somewhat more surprisingly, the profiles also have similar slopes as the LBG stacked profile found by Steidel et al. (2011). On the other hand, the haloes of the average emitters are clearly characterized by a different profile. As we will discuss in Section 5, these similarities may hint to a common fluorescent origin for quasar fluorescent nebulae and the LBG haloes in the fields observed by Steidel et al. (2011). Once the different redshifts of the Slug Nebula and the nebulae found in this study are taken into account, we see in the middle and right panels that, despite its apparent larger size and luminosity, the Slug Nebula is perfectly compatible with the average detection obtained with MUSE.

Returning to the 3D-masks to define the location of the nebulae in the MUSE datacubes, we produced two-dimensional maps of the first and second moments of the flux distribution (in wavelength) to get an indication of the centroid velocity and width of the emission line for each spatial location. These resolved kinematic maps give us the opportunity to detect kinematic patterns, e.g., evidence for rotation, inflows or outflows and to identify kinematically distinct regions in the nebulae, if present (see e.g., Martin et al. 2015; Prescott et al. 2015a for some examples). Because the shape of the emission from intergalactic gas cannot be simply modeled with a single analytic function, such as a Gaussian, we decided to follow this non-parametric approach of computing the moments of the flux distribution within the 3D-mask region. Restricting the analysis to the voxels associated with the CubExtractor 3D-masks significantly reduces the effect of the noise in this non-parametric approach.

In Figure 6, we show the maps of the first moment of the flux distribution, i.e., the flux-weighted velocity centroid shift relative to the peak of the integrated emission of each nebula. While some systems, e.g., #15, show possible evidences of rotation in a disk-like structure, the majority of the nebulae do not show clear evidences from the emission of rotation or other ordered kinematic patterns. The largest nebulae show, instead, strikingly coherent kinematical structures over very large distances (e.g., #1 and #3). We note that the line is in general not the best indicator for kinematics because of radiative transfer effects, but these typically tend to disrupt coherency on large scales rather than enhance them (Cantalupo et al. 2005).

In Figure 7, we present the maps of the second moment of the flux distribution, i.e., the velocity dispersion. For consistency with previous works in the literature we show the Gaussian-equivalent FWHM derived by multiplying the second moment by 2.35. Again, we stress that the line may also be broadened by radiative transfer effects, so from this FWHM we cannot directly constrain the thermal properties of the gas. Nonetheless, the relative comparison between different objects is still informative. This figure clearly shows the main difference between radio-quiet and radio-loud systems: radio-quiet nebulae are narrower (500–700 km s⁻¹) than radio-loud systems (FWHM > 1000 km s⁻¹) in agreement with previous results (e.g.,Villar-Martín et al. 2007). There is only one radio-quiet nebula with broad-line emission, i.e., #6, and this is also the nebula that showed a clearly distinct SB profile (see Section 4.3). This is a striking confirmation that, even when ubiquitously detected, radio-quiet quasar nebula may have a different origin with respect to radio-loud systems. We will return to this point in our discussion below.

Zoom In Zoom Out Reset image size

The MUSE wavelength coverage allows us to search for other extended rest-frame UV lines, namely C iv/ and He ii/. These are typically the brightest UV lines after hydrogen for AGN-photoionized nebulae and their line ratios have traditionally been used to constrain the origin and physical properties of the emitting gas (see e.g., Dey et al. 2005; Villar-Martín et al. 2007, Humphrey et al. 2008; Prescott et al. 2009; Scarlata et al. 2009; Arrigoni Battaia et al. 2015a, 2015b). In particular, the strengths of He ii and C iv with respect to may provide information on the ionization parameter, gas density and metallicity, and eventually on the mechanism responsible for the gas emission. Giant nebulae around high-redshift radio galaxies typically show C iv and He ii as well as other lines (Villar-Martín et al. 2007), which suggests a high metallicity for the gas and therefore an origin related to outflows rather than pristine intergalactic gas accretion.

We again used the 3D-masks produced by CubExtractor to search for spatially coherent C iv and He ii lines with the following strategy. Because the emission peak may be shifted by radiative transfer effects with respect to other lines, we "scanned" around the expected location of these lines by shifting the defined 3D-mask along the wavelength direction. We then compute the total flux of the voxels associated with the 3D-mask at each spectral location. Using the same three-dimensional aperture for and the other lines guarantee that we do not have any aperture effects in the line ratios. However, because we expect the emission to always be brighter (and more extended), this approach gives us very conservative limits in the case of non-detections. In particular, we use the statistics of the "scan" to estimate the noise associated with the aperture photometry performed with the 3D-mask and we consider any line falling within ±3000 km s⁻¹ from the expected line center with a flux higher than a detection. In the case of a non-detection, we use the as a limit.

The results are summarized in Table 2 and presented in Figure 8, where we also compare them with previous studies of high-redshift radio galaxies (HzRG) from Villar-Martín et al. (2007). In all but three radio-quiet systems, we have detected neither C iv nor He ii emission. The anomalous nebula #6 (see above) is the only radio-quiet system with both (marginal) C iv and He ii detections (at and respectively), while nebulae #4 and #14 only have detectable C iv (at and respectively). For all other nebulae, we have only upper limits for C iv and He ii ( in Figure 8). Although we only have conservative limits in the majority of the cases, it is interesting to note that our nebulae seem to lie in a different part of the line-ratio diagram with respect to radio-galaxy halos, thus reinforcing the suggestion from the kinematic analysis of a different origin for these systems. We leave a more detailed analysis of the detected sources in He ii and/or in C iv and the implication of the line ratios for future work.

Zoom In Zoom Out Reset image size

In this study, we have shown that all brightest radio-quiet quasars in our sample, located at the redshift range , are surrounded by giant (pkpc) nebulae with an SB level above ∼10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻². This result seems to be in stark contrast with previous narrow-band and spectroscopic surveys for bright radio-quiet quasars at high redshift. These studies have found detection rates of giant nebula less than 10% using NB imaging at (Cantalupo et al. 2014; Martin et al. 2014; Hennawi et al. 2015; Arrigoni Battaia et al. 2016), while spectroscopic surveys at any redshift have found exclusively smaller nebulae ( pkpc) and only in about 50% of the cases (North et al. 2012; Hennawi & Prochaska 2013; Fathivavsari et al. 2015). Because all these surveys—including our MUSE observations—reached a similar depth, the origin of this discrepancy must be searched for somewhere else. In particular, we explore in this section the possibility that the different observational techniques may be playing an important role in producing such different results. As we will see, however, on the basis of current surveys we cannot completely exclude that a true evolution in the frequency of giant nebulae as a function of quasar redshift and luminosity may be present in the universe based on current surveys, and future work will be needed to further explore this possibility.

Let us assume that there is no redshift evolution in the intrinsic properties of giant nebulae, around radio-quiet quasars and compare our results to NB surveys at . As we have shown in Figure 5, this assumption does not seem unreasonable: the Slug Nebula discovered in NB imaging (Cantalupo et al. 2014) has an SB profile that is perfectly compatible with the MUSE detection once corrected for SB-redshift dimming (which should of course favor detections at lower redshift). As discussed in the Introduction, there are two main effects that will reduce the detectability of such bright nebulae in NB imaging observations: (1) filter losses, if the nebular emission line falls at the edge or outside of the NB filter, and (2) PSF losses, if the quasar PSF hides the presence of the nebula or reduces our ability to detect extended emission. Spectroscopic surveys with a single slit orientation also suffer from slit losses, but they are not affected by filter losses.

From our MUSE and quasar spectra, we can estimate how many of our nebulae would have been missed by NB filter observations because of filter losses. For this thought experiment, we will assume that quasars would have been selected for NB observations if their estimated systemic redshift from corrected broad-line emission (C iv in our case, but typically Mg ii at ) would have placed the expected emission at the peak of the NB filter transmission. In reality, the selection is made on a range of systemic redshifts (typically within 500 km s⁻¹ from the filter center), but we will assume that this range is absorbed by the intrinsic scatter in the C iv-corrected redshift estimation ( km s⁻¹ at 1σ, Shen et al. 2016). Assuming a typical NB filter width of 3000 km s⁻¹ (FWHM), we can notice from Figure 2 that two nebulae would have been totally missed by NB observations (#2, #6) while for the other five (#1, #10, #11, #12, #16) the displacement between the emission and the filter transmission peak would have resulted in a substantial flux loss that would likely have compromised the detectability of the nebulae, especially in the faintest but most extended regions. However, at least nine nebulae (#3, #4, #5, #7, #8, #9, #13, #14, #15, #17), which is at least 50% of the sample, have emission so close to the estimated systemic redshift that they would not have been missed by NB observations, including also the systemic redshift errors from C iv or Mg ii. This factor is much larger than the 10% detection rate estimated from NB imaging. Therefore, we conclude that filter losses cannot completely explain the discrepancy between MUSE detection rates and those of NB surveys, though they can certainly reduce the NB detection rate significantly. We note that radio-loud quasars typically have much broader emission and therefore filter losses would be even less important for these systems. This may partially explain the higher detection rate of radio-loud nebulae versus radio-quiet ones.

Regarding PSF losses, we notice that only about 40% of the MUSE radio-quiet nebulae have a circular shape that could have been difficult to distinguish from the quasar PSF (#5, #8, #12) or that are small and faint enough to make PSF subtraction possibly challenging (#10, #15, #16, #17) for NB imaging or spectroscopy. At least 7 out of 17 nebulae (i.e., about 40%) have clear asymmetric shapes or offsets from the quasar position that should have been easy to identify even with a very approximate, or absent, quasar PSF subtraction (#1, #2, #3, #6, #7, #13, #14). The remaining 20% of the sample (#4, #9, #11) does not fall clearly in any of the previous categories. Despite the small number statistics, we do not think that PSF losses are the main reason for the discrepancy with NB surveys, but they could contribute by reducing the observed frequency of nebulae. This suggestion is reinforced by the fact that the two radio-loud nebulae discovered with MUSE (#R1, #R2) also have circular morphologies and that their detection rate should therefore not have been different from that for radio-quiet systems in NB surveys if PSF losses are important. Taken together, filter losses and PSF losses may result in a decrease by a factor of four, at most, in the detection rate of radio-quiet nebulae in NB observations, which is still not enough to account for the discrepancy with MUSE data.

The discrepancy with the detection rates measured in spectroscopic surveys can, however, be explained by a combination of slit losses and PSF losses. Given the complex and asymmetric morphology of a large fraction of MUSE nebulae, the probability that a single slit orientation would cover a 100 pkpc patch of the extended emission is quite low. There are only a few cases in which the nebulae are symmetric enough to expect a detection over scales extending to 100 pkpc (e.g., #5); however, the difficulty of estimating the exact PSF of the quasar through slit spectroscopy may result in a substantial reduction of the detection rate. This effect may also reduce the detectability of the inner parts of the nebulae to the level found by spectroscopic surveys (North et al. 2012; Hennawi & Prochaska 2013; Fathivavsari et al. 2015). Moreover, the availability of more spatial positions and the large FoV provided by MUSE with respect to a single slit will clearly facilitate the extraction and detection of extended emission around the quasars. Finally, it is easy to explain the discrepancy with other results based on IFS that have a much smaller FoV with respect to MUSE (e.g., Christensen et al. 2006; Herenz et al. 2015): the minimum size of the MUSE detected nebulae is about 13 arcsec—which is much larger than the FoV of 8'' × 8'' of previous IFU surveys—and their SB profile is quite shallow in the central regions because separate sky observations for sky-subtraction were not obtained in these surveys, such extended nebulae could have been removed as a sky feature (see, e.g., Figure 7 in Herenz et al. 2015).

If observational techniques and their limitations are not the only sources of the different detection rates, we should consider what other intrinsic factors in the quasar sample selection may contribute to the difference. In particular, we notice that the MUSE quasar redshifts are systematically higher than those of quasars targeted by previous NB surveys, but are compatible with some of the spectroscopic observations (e.g., North et al. 2012). Similarly, the luminosities of the MUSE quasars are systematically higher than for previous surveys because of the limited availability of quasars with redshifts falling in the same NB filters (e.g., Cantalupo et al. 2014; Arrigoni Battaia et al. 2016) or quasar pairs for the largest spectroscopic survey to date (Hennawi & Prochaska 2013). In particular, the quasars selected for this study are on average two absolute magnitudes brighter (i-band rescaled to z = 2) with respect to the sample of Arrigoni Battaia et al. (2016). We have explored if any of the main nebula parameters affecting the detection rate (e.g., size, luminosity, SB profile) depend on the quasar luminosity in our sample but found no significant correlations in all cases. However, we cannot exclude that a correlation may be present at much lower quasar luminosities with respect to the objects observed in this survey. Future studies extending the range of quasar luminosities probed with MUSE or extending to lower redshifts (e.g., with the Keck Cosmic Web Imager, Morrissey et al. 2012) will help to resolve these questions.

Three mechanisms are able to produce extended and kinematically narrow emission from cosmic hydrogen around a quasar, as observed in our study: (1) recombination radiation due to quasar photo-ionization (Cantalupo et al. 2005, 2014; Kollmeier et al. 2010); (2) "continuum-pumping" or scattering from the quasar broad-line region (see, e.g., Cantalupo et al. 2014); and (3) collisional excitation (so-called "cooling radiation"; Haiman et al. 2000, Fardal et al. 2001; Dijkstra et al. 2006; Cantalupo et al. 2008, Rosdahl & Blaizot 2012). All these mechanisms require the presence of "cool" gas around the quasars (i.e., with temperatures well below 10⁵ K). In addition, recombination radiation and "continuum-pumping" also require that the gas is "illuminated" by the quasar or by some other bright source of UV photons.

Estimating the possible contribution from collisional excitation to the nebular emission is notoriously difficult, given the exponential dependence on temperature of the collisional excitation rate (see, e.g., Furlanetto et al. 2005; Cantalupo et al. 2008) and it requires radiation-hydrodinamical simulations with high spatial and temporal resolution in order to capture non-equilibrium ionization effects in addition to accurate modeling of cooling and heating mechanisms. On the other hand, if the gas is highly photoionized, collisional excitation contribution is always negligible with respect to recombination radiation that has a weak dependence on temperature (in the range of 10³–10⁵ K) and on the ionization state, and it is therefore much easier to estimate. Given the presence of hyper-luminous quasars associated with our nebulae (by construction, in our survey), we think therefore that collisional excitation is the least plausible of the three mechanisms listed above. If this was not the case, the small duty cycle of bright quasars at (e.g., 0.004−0.05, Shen et al. 2007) and our detection rate of 100% would imply that giant and bright nebulae should be about a hundred times more common than quasars at this redshift and this is certainly not observed.¹⁵

A fluorescent origin—either from recombination radiation or "photon-pumping"—for the extended emission detected with MUSE implies that we are observing cosmic gas that is both relatively "cold" ( ∼ 10⁴ ) and "illuminated" by a bright UV source, giving us the opportunity to constrain both the amount of cold gas around quasars and the geometry of their UV emission. In particular, from our 100% detection rate of giant nebulae, we can infer that "cold" gas is ubiquitous out to at least 50 projected physical kiloparsecs from any hyper-luminous quasar at , independent of the quasar emission opening angle. This distance covers about half of the virial radius for the expected dark matter halo associated with quasars at these redshifts ( , independent of luminosity, see, e.g., Shen et al. 2007), which instead are expected from numerical simulations to contain mostly hot, virialized gas at temperatures of K (see Cantalupo et al. 2014). We stress that we do not see any kinematic signatures of extended disks or outflows, and therefore we believe that the origin of this gas is intergalactic. In some cases, we have detected cold gas up to a projected distance of about 200 pkpc from the quasar (e.g., #3), extending up to at least twice the virial radius of the typical dark matter halos associated with quasars at this redshift.

Any constraints on the true total extent of the cold gas distribution are, however, degenerate with the geometry and opening angle of the quasar UV radiation. Turning this argument around, if the cold gas around quasars is extended (in all directions) up to the projected value that we measure in the images, then the quasar emission would have to be isotropic or covering at least 180° if coming from a single "cone." Assuming instead that the dense gas that we can trace in emission always extends to at least 200 pkpc from the quasars (using nebula #3 as a reference, or assuming that we trace the same optically thick gas detected in absorption by Prochaska et al. 2013), then we will have a typical opening angle of about 25% to 50% of the total solid angle for biconical emission, or half of these values for a single "cone" emission. Future studies targeting both gas in emission and absorption around individual quasars (e.g., using quasar pairs) will help to break the degeneracy.

What are the physical properties of this ubiquitous cold gas around bright quasars? From the SB profiles corrected for redshift dimming and from the morphologies, we find that the observed emission properties of the nebulae are remarkably similar to that of the Slug Nebula (Cantalupo et al. 2014), despite the different redshifts, overall sizes, and luminosities. On the one hand, our results show that these structures are common in the universe around bright quasars and, on the other hand, all the implications of the Slug Nebula apply also in the case of the nebulae detected in this study. In particular, the large values of the SB and the lack of clear He ii and C iv emission would imply large cold gas densities, unless a large clumping factor (with clump sizes of the order of few tens of parsecs, see Cantalupo et al. 2014; Arrigoni Battaia et al. 2015a) is invoked.

Remarkably, the LBG stacked SB profile obtained by Steidel et al. (2011) also has a very similar slope as to the Slug Nebula and MUSE detected nebulae. This profile is very different from the results of other surveys targeting Lyα emitting galaxies both in large areas and using deep integrations (e.g., Wisotzki et al. 2016). How can we explain this apparent discrepancy? In view of our results, and considering that the quasar phase should only be a small part of the life of a galaxy, we must conclude that "cold" gas is always present around (massive) galaxies as well, and that the only difference for the appearance or not of a large nebula is the illumination factor provided by nearby or internal quasars. At this point, it is worth noting that the LBGs observed by Steidel et al. (2011) are located in three fields of which one contains a hyper-luminous quasar (at the same redshift of the galaxies) and the other two contains two large overdensities of galaxies and AGNs that might boost the ionizing UV background by a significant factor. This observation might then also be compatible with a fluorescent scenario, and it would imply the same origin for all giant extended emission around both quasars and galaxies.