The Stacked Lyman-Alpha Emission Profile from the Circum-Galactic Medium of Quasars Based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), Ministério da Ciência, Tecnologia e Inovação (Brazil) and Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina)..
In the context of the FLASHLIGHT survey, we obtained deep narrow band images of 15 quasars with GMOS on Gemini-South in an effort to measure Ly emission from circum- and inter-galactic gas on scales of hundreds of kpc from the central quasar. We do not detect bright giant Ly nebulae (SB erg s cm arcsec at distances kpc) around any of our sources, although we routinely (%) detect smaller scale kpc Ly emission at this surface brightness level emerging from either the extended narrow emission line regions powered by the quasars or by star-formation in their host galaxies. We stack our 15 deep images to study the average extended Ly surface brightness profile around quasars, carefully PSF-subtracting the unresolved emission component and paying close attention to sources of systematic error. Our analysis, which achieves an unprecedented depth, reveals a surface brightness of SB erg s cm arcsec at kpc, with a detection of Ly emission at SB erg s cm arcsec within an annulus spanning kpc kpc from the quasars. Assuming this Ly emission is powered by fluorescence from highly ionized gas illuminated by the bright central quasar, we deduce an average volume density of cm on these large scales. Our results are in broad agreement with the densities suggested by cosmological hydrodynamical simulations of massive () quasar hosts, however they indicate that the typical quasars at these redshifts are surrounded by gas that is a factor of times less dense than the () gas responsible for the giant bright Ly nebulae around quasars recently discovered by our group.
Subject headings:galaxies: formation — galaxies: high-redshift — intergalactic medium
Over the past decade galaxy formation studies have increasingly emphasized the need to understand the gaseous phases that surround the galaxies themselves (e.g., Keres2005; Ocvirk2008; Fumagalli2011; Stinson2012; Fumagalli2014; Faucher-Giguere2015), which represent both the fuel for and product of star-formation. Specifically, the circum-galactic medium (CGM), which represents gas extending from kpc from galaxies, encodes information about the complex interplay between outflows from galaxies and accretion onto them from the intergalactic medium (IGM) (e.g., Fraternali2015; Nelson2015). The CGM has been typically studied by analyzing absorption features along background sightlines at small impact parameter from galaxies (Croft2002; Bergeron2004; Rudie2012; Crighton2015; Nielsen2015; Crighton2013) and quasars (Bowen2006; Hennawi2006; Hennawi2007; Prochaska2009; Hennawi2013; Farina2013; Farina2014; Prochaska2013; Prochaska2013b; Johnson2015).
Specifically, using this technique with projected quasar pairs, it has been shown that massive galaxies hosting quasars at are surrounded by a massive ( M), enriched ( Z), cool ( K) CGM with high covering factor ( % within the virial radius of 160 kpc; Prochaska2013; Prochaska2013b; Prochaska2014). Notwithstanding great effort, state-of-the-art simulations of galaxy formation still struggle in reproducing these observations (Rahmati2013; Fumagalli2014; Faucher-Giguere2015; Meiksin2015, but see also Rahmati2015; Faucher-Giguere2016). Nevertheless, these absorption studies are limited by the paucity of bright background sources and the inherently one-dimensional nature of the technique, and hence do not paint a complete picture of the CGM. Direct observations of the CGM in emission would be highly complementary, and this new observable would enable a more constrained comparison to the results of simulations.
In particular, it has been suggested that the Ly transition should be the easiest line to detect from the quasar CGM, and from the IGM on larger scales. Indeed, gas in the CGM and IGM could reprocess impinging UV radiation resulting in fluorescent Ly emission. Although the fluorescence signal powered by the ambient UV background (Hogan1987; Binette1993; Gould1996; Cantalupo2005), with an expected surface brightness (SB) SB erg s cm arcsec (e.g., Rauch2008), is still too faint to be detectable, this emission can be boosted up to observable levels by the additional ionizing flux of a nearby quasar (Rees1988; HaimanRees2001; Alam2002; Cantalupo2005).
Given the extreme challenge of detecting this faint quasar powered fluorescent emission, previous searches have been difficult to interpret (e.g., Fynbo1999; Francis2004; Cantalupo2007; Rauch2008). The most compelling candidates for quasar powered Ly fluorescence are the compact sources, a.k.a. “dark galaxies”, at Mpc distances from quasars discovered by Cantalupo2012. Indeed, they reported an excess of Ly emitters (LAEs) with rest-frame Ly equivalent width exceeding the maximum expected from star-forming regions, i.e. Å (e.g., Charlot1993).
At smaller distances from the quasar , in the CGM and at the CGM-IGM interface, one also expects gas to be illuminated. Motivated by this, many have searched for Ly emission close to the quasar (e.g., HuCowie1987; Heckman1991spec; Heckman1991; Christensen2006; Hennawi2009; North2012; Hennawi2013; Roche2014). While the consensus emerging from this work is a relatively high rate of detection of extended Ly emission on small-scales (radius of kpc), convincing evidence for larger scale emission remains elusive.
Recently more systematic searches with higher sensitivity have uncovered several very extended Ly nebulae around QSOs (Cantalupo2014; Jackpot; Trainor2013; Martin2014a; Martin2015). In particular, our group discovered the two most striking examples of giant Ly nebulae around the radio-quiet quasars: UM 287 (Cantalupo2014), and SDSSJ0841+3921 (Jackpot). With respective sizes of 460 kpc and 310 kpc, these constitute the largest nebulae ever observed around a QSO, signifying large reservoirs of cool ( K) gas. Their bright Ly emission has been explained as fluorescent recombination emission from the CGM/IGM powered by the intense ionizing radiation emitted by the central quasar. Under this interpretation, independent arguments based on cosmological zoom-in simulations post-processed with radiative transfer (Cantalupo2014), as well as the non-detection of He ii and C iv emission from the nebula (FAB2015), indicate the Ly emission from the UM 287 nebulae comes from a large population of compact ( pc), dense ( cm), cool gas clumps (FAB2015). Jackpot came to similar conclusions about the gas responsible for the Ly emission around SDSSJ0841+3921 based on joint modeling of the nebula in absorption and emission.
Interestingly, an increasing number of independent observations point in the same direction, suggesting high densities (comparable to the interstellar medium) and small cloud sizes in the quasar CGM (Prochaska2009; QPQ8). However, note that only a fraction of objects studied through absorption line techniques require high ISM-like densities (e.g., Prochaska2009; QPQ8), although several assumptions are required in order to constrain this parameter (e.g. photoionization models). The densities found in giant Ly nebulae are orders of magnitude higher than predictions from simulations on CGM scales ( cm; e.g., Rosdahl12), probably reflecting the inability of current simulations of galaxy formation to resolve small scale structure in the cool CGM (Crighton2015; FAB2015; Jackpot).
High density clumps have also been invoked to explain the extended emission line regions (EELR; e.g., Stockton2006) around active galactic nuclei (AGN; Stockton2002; FuStockton2007). The narrow lines ( km s) around AGN can extend out to tens of kpc, and are thought to be mainly powered by the AGNs (e.g., Husemann2014). Given their high surface brightness erg s cm arcsec, the searches for Ly emission around quasars clearly detect the Ly signal from these regions (e.g., Christensen2006).
Although the discovery of bright high surface brightness Ly nebulae on large scales around quasars provides an important new diagnostic of the CGM, revealing problems with our understanding of cool gas in massive halos, they represent a relatively rare phenomenon and it is unclear whether the extreme properties implied by these nebulae are representative of the quasar CGM as a whole. Indeed, as summarized in Jackpot, three independent surveys (FAB2013; Hennawi2013; Trainor2013) suggest that only 10 % of quasars exhibit giant Ly nebulae, i.e. with high SB emission ( erg s cm arcsec at distances kpc from the quasar). It is still unclear whether this frequency of detection results from different physical properties of the surrounding gas (e.g. higher densities in the giant Ly nebulae), and/or variations in the solid angle into which quasars emit their ionizing radiation (Cantalupo2014; Jackpot), and/or peculiar environment of these sources (Jackpot).
In this paper we conduct a stacking analysis of narrow-band (NB) images targeting the Ly line for 15 QSOs that do not show giant Ly nebulae. Our goal is to measure the average surface brightness of the quasar CGM, an important quantity for planning future integral field unit (IFU) observations with new instruments such as the Multi Unit Spectroscopic Explorer (MUSE; Bacon2010) and the Keck Cosmic Web Imager (KCWI; Morrissey2012). These data have been collected in the framework of the larger Fluorescent Lyman-Alpha Survey of cosmic Hydrogen iLlumInated by hIGH-redshifT quasars (FLASHLIGHT, see also FAB2013), which aims to better constrain the fluorescence signal on both CGM and IGM scales around QSOs. This survey has been conducted using custom-built narrow-band filters with a program using the Gemini Multi-Object Spectrograph (GMOS, Hook2004) on Gemini-South, and a second one using the Low-Resolution Imaging Spectrometer (LRIS, Oke1995) on Keck. While the on-going program on Keck resulted so far in the discovery and study of the aforementioned giant Ly nebulae around UM 287 (Cantalupo2014) and SDSSJ0841+3921 (Jackpot)111Note that this nebula was previously detected with long-slit observations by Hennawi2013., in this work we focus on the data collected at Gemini-South.
To our knowledge, similar Ly stacking analysis has never been conducted around QSOs, but have been explored in the case of LAEs and Lyman break galaxies (LBGs), which exhibit extended Ly halos on scales of tens of kpc. In particular, Rauch2008 showed that the average and median SB profiles for 27 line emitters at extend out to kpc at a level of SB erg s cm arcsec. After this initial study of individual objects based on a 92 hour ultra-deep longslit observation, subsequent workers chose to instead stack the data of several objects observed with much shorter exposure times. Specifically, Steidel2011 reported evidence for extended Ly emission out to kpc at a similar SB level around LBGs by stacking 92 objects. Matsuda2012 performed a stacking analysis for LAEs and 24 LBGs, confirming the result of Steidel2011 in the case of LBGs, and showing a probable environmental dependence of the size and luminosities of their average LAE Ly profiles, with larger and brighter Ly halos in the richest environment, though always fainter than LBGs. Momose2014 expanded the work by Matsuda2012 up to also detecting extended Ly halos. Recently Wisotzki2016, exploiting the higher sensitivity of the MUSE instrument (Bacon2010), showed that the Ly SB profiles of individual high- galaxies can be studied down to SB erg s cm arcsec (on scales of tens of kpc) if very long 27 hour exposure times are employed.
Note that several mechanisms can in principle power extended Ly halos, making it difficult to provide a unique interpretation, especially when only Ly images are available. The main mechanisms invoked in the literature, which may even act together include: (i) resonant scattering of Ly line photons (Dijkstra2008; Steidel2011; Hayes2011), (ii) photoionization by stars or by a central AGN (Geach2009; Overzier2013), (iii) cooling radiation from cold-mode accretion (e.g., Fardal2001; Yang2006; Faucher2010; Rosdahl12), and (iv) shock-heated gas by galactic superwinds (Taniguchi2001). An exemplar case for the ambiguity that can arise in interpreting Ly emission on scales is the ongoing debate about the nature of the Ly blobs (LABs), i.e. large ( kpc) luminous ( erg s) Ly nebulae at (e.g., Steidel2000; Matsuda2004; Dey2005; Prescott2009; Yang2011; FAB2015LAB) – although there is increasing evidence that LABs are frequently associated with obscured AGN (Geach2009; Overzier2013; Prescott2015). While the brightest giant Ly nebulae discovered around luminous quasars appear to be unambiguously powered by photoionization from the central source (Cantalupo2014; Jackpot), at the much lower SB erg s cm arcsec levels that we probe in this study, other emission mechanisms are also plausible. In this work we focus mainly on the photoionization scenario, however additional data and a detailed comparison with cosmological simulations of very massive systems are needed to assess the contribution from the different powering mechanisms.
This paper is organized as follows. We summarize our observations in §2. In §3, we describe the data reduction procedures, the global surface brightness limits of our data, the masking technique used before the SB profile estimate, and the characterization of the point-spread-function (PSF). In §4 we present the stacking analysis around the 15 QSOs in our sample and the resulting extended Ly SB profile. In §4.1 we test the significance of our results, while in §5 we interpret our results in the context of the simple model for emission from the CGM presented in Hennawi2013, and compare to other measurements in the literature. Finally, we summarize and conclude in §6.
Throughout this paper, we adopt the cosmological parameters km s Mpc, and . In this cosmology, 1″ corresponds to 8.2 physical kpc at . All magnitudes are in the AB system (Oke1974), and all distances are proper.
In March 2013 we successfully installed a custom narrow-band filter for the Gemini Multi Object Spectrograph (GMOS, Hook2004) on the Gemini South telescope, targeting Ly emission around quasars at ( Å, Å, peak transmission T%)222The filter is still available at the telescope for observations. Note that in June 2014, after our last run, the CCDs of GMOS-S were upgraded to the new Hamamatsu CCDs, more sensitive at longer wavelengths, e.g. Å. For this reason, the overall efficiency of the system at 4000 Å has been degraded by 25% (for the blue sensitive CCD). . In one year we observed a total of 17 quasars. Three of these have longer integrations, typically of 5 hours in a series of dithered333For all our observations executed a random dithering pattern within a 15″15″ box, in order to fill the two CCD gaps of arcsec each. 1800s exposures, achieving a depth of SB erg s cm arcsec ( in 1 square arcsec). The other 14 quasars were observed in a fast survey mode, to try to uncover new giant Ly nebulae similar to UM 287 (Cantalupo2014) and SDSSJ0841+3921 (Jackpot). The ‘fast survey’ observations were carried out using typical exposure times of hours in a series of dithered 1200s exposures, achieving an average depth of SB erg s cm arcsec ( in 1 square arcsec).
To study a representative sample of quasars we selected our targets from the SDSS/BOSS catalogue (Paris2014). In particular, we focus on the brightest quasars within our custom narrow-band filter which do not have bright stars in their vicinity, or significant Galactic extinction in their direction. As NIR spectra are not available for most of the SDSS QSOs, accurate redshifts are determined using the Mg ii line, which has a small and known offset from the systemic redshift of the quasar (Richards2002; Shen2016). Indeed, the typical uncertainty on these redshift estimates is (or equivalently km s), which is much smaller than the width of the narrow-band filter used, i.e. (or equivalently km s). Thus, to be sure that the Ly emission of our targets falls within the narrow band filter, we selected only QSOs whose redshift gives a maximum shift of Å from the filter’s center (or equivalently km s). In the Appendix we quantify the flux losses due to the uncertainty on the systemic redshift of our sample to be of the order of 3%, and thus negligible.
In addition to the NB data, we have observed the same fields in -band. For the long integrations the typical exposure time was 3 hours in a series of dithered 240s exposures, whereas in the fast survey we observed for 40 minutes in a series of dithered 300s exposures. The ‘fast-survey’ data were collected in March 2014 in a 4 night run (program ID: GS-2014A-C-2), while the longer exposures were obtained in service mode during 2013-2014 (program ID: GS-2013A-Q-36).The seeing was 0.7″ for the service mode program, while it ranged between 0.5″ and 1.9″ during the visitor run, with a median seeing of 1.2″. A 4 night run (program ID: GS-2013B-C-4) was completely lost in September 2013, due to bad weather conditions. The observations taken with GMOS/Gemini-South and used in this work are summarized in Table 1.
|Target||-mag||Exp. Time NB||Exp. Time||Depth|
3. Data Reduction
As our goal is to detect extremely low surface brightness emission for which control of systematics is crucial, we opted to write our own custom data-reduction pipeline in the Interactive Data Language (IDL), and used only the gmosaic routine in the publicly available Gemini/IRAF GMOS package 444http://www.gemini.edu/sciops/data-and-results/processing-software, to correctly mosaic the three different chips of the CCD. After mosaicing, we used our custom pipeline to bias-subtract the images, and flat-field them using twilight flats. To improve the flat-fielding, essential for detecting faint extended emission across the fields, we further corrected for the illumination patterns using night-sky flats. The night-sky flats were produced by combining the unregistered science frames with an average sigma-clipping algorithm after masking out all the objects. Satellite trails, CCD edges, bad pixels, and saturated pixels were masked. In particular, we created a bad-pixel mask using an average dark frame obtained from dark images with the same exposure time and binning as our science frames. We found that creating a custom bad-pixel mask was necessary because of the increased number of bad pixels arising from our non-standard binning of pixels (resulting in a pixel scale of 0.29″ pixel). Given the low count levels through our NB filter, this binning was necessary to minimize read-noise, resulting in nearly background limited observations for our choice of NB exposure time ( s, see §2). Each individual frame was cleaned from cosmic rays using the L.A.Cosmic algorithm (vanDokkum2001).
The final stacked images in each filter (NB and -band) were obtained by averaging the science frames with weights and masking pixels rejected via a sigma clipping algorithm. Sky-subtraction was performed on the individual frames before the final combine. To ensure that we do not mistakenly subtract any extended emission, or introduce systematics during the sky subtraction, we simply estimated an average sky value for each individual image, after masking out all sources, and subtracted this constant value from each image. Given the very small field distortions of GMOS555The field distortions of GMOS-S are regularly checked. For the E2V detector they were estimated to be 1 pixel in the x-direction and 2 pixels in the y-direction at the edge of the field at the time of observations (Pascale Hibon, private communication). Note that these values correspond to the unbinned CCD, and thus 1 pixel corresponds to 0.07″, resulting in distortions significantly smaller than the pixel scale of our binned images (0.29″ pixel). We confirmed this small level of distortion by comparing the SDSS-DR9 catalogue with stars in our individual frames., we calibrated the astrometry after the stacking666In this way, we avoided the challenge of finding a good astrometric solution for individual shallow NB frames with a limited number of sources.. The astrometry was calibrated with the SDSS-DR9 catalogue using SExtractor (Bertin1996) and SCAMP (Bertin2006). The RMS uncertainties in our astrometric calibration are 0.2″.
For the flux calibration of the NB imaging, we use the standard stars (H600, G60-54, G138-31) that were repeatedly observed during the observations, typically at the beginning and at the end of the night. To avoid systematics in the flux calibration, these spectrophotometric stars are selected to be free of any feature at the wavelength of interest, and to have a good sampling of their tabulated spectrum (at least 1 Å). All data show consistent zero points during different nights, and are stable during each single night (ZP). Regarding the -band images, the flux calibration is performed using several photometric stars in different PG-fields (PG0918+029, PG1047, PG1633+099) observed at the beginning and at the end of each night. Also in this case, all data show consistent zero points, with good agreement between different nights, and stability during each single night (ZP). Note that our -band zero point value is consistent with the GMOS tabulated values777http://www.gemini.edu/sciops/instruments/performance-monitoring/data-products/gmos-n-and-s/photometric-zero-points. The uncertainty in the derived zero-points is mag.
During the data reduction steps (e.g. bad pixel masking, satellite trails masking, etc.), our pipeline consistently propagates the errors, and produces a variance image for each stacked image, which contains contributions from the sky and background photon counting as well as read noise. An accurate final noise model is of fundamental importance for the measurement of very low surface brightness signals (Hennawi2013; FAB2015LAB; FAB2015). We compute a global surface brightness limit for detecting the Ly line using a global root-mean-square (rms) of the images. To calculate the global rms per pixel, we first mask out all the sources in the images, paying particular attention to the scattered light and halos of bright foreground stars, and then compute the standard deviation of sky regions using a sigma-clipping algorithm. We convert these rms values into the surface brightness limits per 1 sq. arcsec aperture listed in Table 1.
3.1. Image Preparation for the Profile Extraction
Figure 1 shows an example of our final stacked images for the NB (left) and -band (right) filters for SDSSJ121503.13+003450.6. These images show the full field of view (FOV) of our single field observation, i.e. 5.5′ 5.5′. From Figure 1, it is clear that our images are nearly flat in both NB and -band and lack significant large-scale residuals, indicating that our flat-fielding and background subtraction are good, even though we simply subtracted a constant background from each individual frame (see § 3). The lower corners of the images show the supports of the CCDs, and are thus noisier in the final stack. These parts of the images are not used in our analysis, but we show them here for completeness. Note that the noise at the position of the two vertical CCD gaps ( arcsec each) is slightly higher. This is more evident in the NB stack compared to the -band image. Indeed, given the typical number of only 6 dithered NB images of 1200s, the gap has been filled, but results in slightly shallower data for a small fraction of the image888Note that the targeted quasars are ″( kpc) away from the CCD gaps, i.e. they are close to the center of the GMOS-S field-of-view..
Interestingly, in the NB data presented in Figure 1 a LAB has been identified at 69″ from the quasar with an average erg s cm arcsec, with a size of 7.8″ ( kpc). This source confirms that we are able to detect extended emission at the level planned for our observations, at the characteristic SB of the giant Ly nebulae around UM 287 and SDSSJ0841+3921.
As we are interested in the Ly line emission of the gas distribution around the QSO, one might naively think to compute a continuum subtracted image, e.g. as performed in Cantalupo2014. However, continuum subtraction would inevitably increase the statistical and systematic error budget (e.g. PSF matching, diffuse light, large galaxies in the -band, unknown spectral slope of the stars/galaxies, and the use of two images instead of only one), significantly decreasing our ability to detect the expected weak diffuse Ly emission signal. Furthermore, in the halo of a QSO we do not expect to detect extended diffuse continuum emission at a level comparable to the diffuse Ly emission999Note that the two-photon continuum expected from the same gas emitting the extended Ly emission is of the order of erg s cm Å (see Figure 9 of FAB2015), and thus negligible in comparison to the Ly emission itself..
For these reasons, we decided to avoid the continuum-subtraction step, and instead mask all the continuum sources detected in the -band. To build these masks we run SExtractor on the -band images to identify all the sources down to a very low detection threshold (DETECT_THRESH1.0), and allowing the detection of very compact objects (DETECT_MINAREA5 pixels). We then use the ‘segmentation’101010The ‘segmentation’ image is already a mask of the image, in which each source is represented by its total isophotal area with flux equal to the identification number in the SExtractor catalogue. image produced by SExtractor to create the final object mask for each NB stack. Given that SExtractor unambiguously assigns an identification number to each source in the field, we can ‘switch off’ the mask for the sources we are interested in, i.e. the QSO or the stars used for the PSF comparison in our analysis.
Furthermore, to ensure that we do not mistakenly detect Ly signal from compact objects in proximity of the QSOs, i.e. the LAEs which do not have a continuum detection, we produce an analogous mask using the NB image. However, this mask targets only compact sources but does not mask diffuse emission, i.e. DETECT_MINAREA=5 pixels and DETECT_MAXAREA=15 pixels. Note that this NB masking step removes only a small fraction of pixels from our analysis, i.e. % from the whole image, or % within the 1′1′ FOV centered on the quasar, with little variation from field to field. To generate the masks, we do not convolve the images with a spatial filter in either of the SExtractor runs.
We then generate a final mask by combining the two individual masks, so that both pixels from the -band and NB masks are neglected in our analysis. After this masking procedure has been adopted, the fraction of the FOV usable around each quasar on average is about %, while within a region of 1′1′centered on the quasar it is about %, with obvious variations from field to field.
Figure 2 shows the NB image of SDSSJ121503.13+003450.6 this time after masking all the sources except the QSO, using the final mask. It is clear that our method does not mask extended emission in the NB images. Indeed, although there are compact continuum sources within the LAB detected in the NB image, the extended emission is still clearly visible after applying the final mask (compare Fig. 1 with Fig. 2). For reference, in this figure we overplot in red the annuli within which we calculate the radial SB profiles in the next sections.
3.2. Individual NB Images vs Giant Ly Nebulae and EELR
In Figure 3 we show images in the NB filter and -band for all the 15 quasars observed with GMOS-S in the FLASHLIGHT survey. We also show -images in the NB filter with the final mask for each field (black) superposed. The images are computed as , where is the final stack and is the square root of the variance image. Indeed, as previously described, our custom pipeline results in a final stack for each filter (NB and -band), and a final variance image (), which contains all the information on the error budget. In the images, emission will be manifest as residual flux, inconsistent with being Gaussian distributed noise (if the noise has been correctly propagated, and thus the image is an accurate description of the noise in the data). From these images it is clear that we have not detected any giant Ly nebulae similar to UM 287 (Cantalupo2014) or SDSSJ0841+3921 (Jackpot), i.e. SB erg s cm arcsec on scales kpc, even though we should have been able to easily detect such giant Ly nebulae given the SB limits quoted in Table 1.
Furthermore, we have also searched our data for emission on smaller scales, i.e. the emission typical of EELR (e.g., Husemann2013; Husemann2014) or so called “Ly fuzz” (e.g., Heckman1991; Christensen2006; Hennawi2013; Herenz2015), by subtracting a crude estimate of the unresolved emission of each quasar. We proceed as follows. First, as we are interested in the Ly line emission, we used the -band images to evaluate the continuum emission of the quasar and subtract it from the NB images, following Yang2009 and FAB2015LAB, as the Ly line is within the -band filter (but close to the edge, at low transmission, i.e. %):
where is the flux in the -band, is the flux in the NB filter, and represent the FWHM of the -band and NB filter, respectively, is the flux density of the continuum within the -band, and is the Ly line flux. The parameter allows a better match to the continuum, and ranges between 0.6 and 1.0 in our sample, reflecting our ignorance on the slope of the continuum for individual sources, as probed by the and NB filters. Note that this range of value is in agreement with the expected continuum slope of quasars in this wavelength range (VandenBerk2001).
After the continuum subtraction, the images contain information only on the Ly emission, but still contaminated from the unresolved line emission of the quasar. To removed the unresolved Ly emission, we have assumed that the -band image of the quasar represents a good estimate of the point-spread-function (PSF) of our observations111111As the seeing in the -band and NB are in agreement, this is true for both -band and narrow-band (see section 3.3 for more detailed analysis of the PSF). Note that in this analysis we used the masked images.. We have then rescaled the peak of this PSF model to the peak of the quasar Ly image, and subtracted it. In this way, we are left with a map of the Ly line which is cleaned from the unresolved line emission emerging from the broad-line regions (BLR; e.g, Stern2015) of the quasar, and the unresolved continuum emission, e.g. from the accretion disk or host galaxy.
In Figure 4 we show the results of this attempt. For each quasar we present a 20″20″ Ly maps in surface brightness units, and the correspective -images, with both images masked using the final mask computed using the aforementioned method. The PSF subtraction reveals Ly emission above SB erg s cm arcsec extending on scales of radius arcsec ( kpc) around 7 out of 15 quasars, i.e. % of the sample, while the other objects show more compact residuals, or residuals affected by systematics in the PSF subtraction. Indeed, the other 8 sources (indicated with a red square in the top-left corner) show asymmetrical residuals, with positive and negative fluctuations clearly visible in the images presented in Figure 4. These residuals could be due both to systematics introduced by the use of the -band of the quasar as PSF, and/or to intrinsic asymmetries in the emission on these small scales. A careful PSF subtraction to study the emission on 50 kpc scale is beyond the scope of this work, but we stress that would be more accurate in the case of integral field spectroscopy (e.g., Husemann2014). Note that, probably because of these systematics in the PSF subtraction, our detection rate for these ‘Ly fuzz’ is at the low end of the frequency range of % deduced in the literature (Christensen2006; Courbin2008; Hennawi2009; North2012), but higher, as expected, than long-slit spectroscopic searches, i.e. % (Hennawi2013).
In summary, the analysis of the individual objects do not present any giant Ly nebula similar to UM 287 (Cantalupo2014) or SDSSJ0841+3921 (Jackpot), but only Ly emission on scales kpc. We ascribe this emission to the EELR or the so called “Ly fuzz”, or in other words, to gas photoionized by the central AGN or by star-formation in the host galaxy of the quasar (Heckman1991; Husemann2014). Other mechanisms could also contribute, e.g. shock-heated gas, scattered Ly emission. In the subsequent sections we perform the stacking analysis to search for the signal from the quasar CGM, i.e. on larger scales.
3.3. Constraining the PSF of the NB Images
To quantify the extended emission in the Ly line around the quasars in our fields, i.e. not related to unresolved emission from the quasar itself, we have to carefully estimate the PSF of our NB images. We need to compute the PSF with high accuracy, in order to detect even small deviations from its shape in the average profile of the QSOs. For this purpose, we match the SDSS star catalogue with the sources in our 15 NB fields and select all the high S/N stars. This operation resulted in a sample of 115 usable stars. For each of these stars, we have created the NB masked image as explained in §3.1, and calculated the stars emission profile in radial logarithmic bins (see red circles in Figure 2) with an unweighted average of all the not-masked pixels within each annulus. We then consistently propagate the errors from the variance image121212In the remainder of the paper we use the same method for the profile extraction.. The profiles of the 115 stars are then averaged with equal weights to obtain the PSF of the NB images131313Note that we use the maximum number of available stars to minimize the errors, even if this led to a different number of stars per image. This effect could in principle bias the PSF characterization towards an image. However, we have found 5-10 usable stars per field, with only two fields without available stars for our analysis. This sampling should not affect significantly our results given that the data were taken under similar weather conditions..
Figure 5 shows the average radial profile for the 115 high S/N stars in our NB images (black points). The profile is consistent with a Moffat function, as expected for seeing limited observations, where the PSF is determined by the wiggling of the stars on the CCD (Trujillo2001). The normalized Moffat profile is defined as
where the full width at half maximum is given by FWHM, and the total flux is normalized to 1. We fit the average profile of the stars with this function and obtained best fit parameter values of FWHM=1.2″and . The FWHM value we obtained is perfectly consistent with the independent measurement of the median seeing in our observations (see §2) estimated by selecting all the stars in the images, and applying the psfmeasure task within the IRAF software package. In the absence of imperfections in telescope optics, the parameter in eqn. (3) should have a value of , as expected from turbulence theory (e.g., Saglia1993; Trujillo2001). However, it is known that PSFs typically measured in real images have larger “wings”, or equivalently, smaller values of . This is due to the fact that the real seeing also depends on the performance of the telescope optics, and not only on the atmospheric conditions (Trujillo2001). Our value is in agreement with this picture, however there is no tabulated PSF for the GMOS-S instrument (German Gimeno141414German Gimeno is the current Instrument Scientist for GMOS-S., private communication).
To test our PSF estimate, we perform the same calculation on masked -band images. We select 15 high S/N not saturated stars, i.e. one per field to include in our calculation all possible seeing variations between different stacks. We then calculate the average radial profile, as for the stars in the NB images. Figure 5 shows (in magenta) this average radial profile after normalizing it to the NB data points. The normalization factor is computed as the ratio between the peak values of the two average profiles. The shape of the Moffat profile of the NB and -band images are in remarkable agreement. Note that in the case of the -band PSF calculation, we probe only smaller radii, due to confusion arising mainly from diffuse light from bright objects. Finally, note that the Moffat profile goes smoothly to zero at large radii, and the average star profile is consistent with this within our uncertainties (see log-linear region of the plot in Fig. 5).
Having characterized the PSF of our observations, we can now ascertain whether our sample of QSOs exhibit a signal from the surrounding diffuse gas distribution. To do so, we apply the same stacking analysis to the quasar sample.
As done for the stars, we have computed the radial profile for each quasar from the NB masked images. The 15 quasar profiles are then averaged to obtain the average signal around quasars. Figure 6 shows the average radial profile of the QSOs, together with the average profile of the stars normalized to the QSOs data points. The normalization has been performed by multiplying the average star profile by the ratio between the peak values of the two average profiles, i.e. the first bin. As the first annulus for the profile extraction has a radius of 1.2″ this is also consistent with an approximate normalization in flux. In addition, for comparison purposes, in this plot we slightly shift the stars’ PSF profile to smaller radii, i.e. the profiles are estimated with the same radial bins, but are shown at different locations to avoid the overlap of the error-bars. The QSOs’ profile clearly deviates from a pure Moffat function at large radii ( kpc), and at a low surface brightness level of SB erg s cm arcsec, and, as expected, is consistent with zero at very large radii ( kpc).
This extended emission is more evident if we take the difference between the QSOs’ profile and the normalized stars’ profile, or in other words, if we remove the contribution from the unresolved emission of the QSOs. The resulting profile thus represents the average Ly emission arising from the distribution of the cool gas ( K) around a typical bright QSO. Figure 7 displays the profile of the difference between QSOs and the normalized star profile, which extends for hundreds of kpc to very low SB levels SB erg s cm arcsec. However, we stress again that we do not subtract the continuum emission from the NB images and thus, on small scales, i.e. the scale of the quasar’s host (tens of kpc), the signal we see might be due to contamination from the resolved host galaxy or is emitted within the aforementioned EELR. If this is due to the host galaxy, we estimated from our profile that the galaxy should have erg s cm Å within an aperture of 3 arcsec radius ( kpc), or equivalently a magnitude of 23.9. This value is similar to the rest-frame UV apparent magnitude of a bright Lyman break galaxy (LBG) at comparable redshift (e.g., Reddy2009), and thus brighter than a galaxy for LBGs, i.e. (; Adelberger2000).
In addition, we check if the signal we detect is dominated by a single bright object in the stack by repeating the same experiment each time removing a different QSO from our sample. Specifically, we compute 15 average profiles of 14 QSOs, each time removing a different QSO. Then, to obtain the CGM profiles, we subtract from each of them the normalized star profile. In this way, it would be clear if a bright QSO in our sample dominates the stack at any radius.
Figure 8 shows the results of this test. It is clear that the profile is not changing significantly between different calculations for kpc. For smaller radii, as explained above, differences could arise due to seeing variations or host galaxy contamination. However, this is not the focus of the current work, and we are confident that none of the QSOs in our sample dominates the average profile on the large scales kpc characteristic of the CGM. However, given that the signal is close to our detection limit we verify the significance of our results in the next section.
4.1. Testing the Significance of the Results with a Monte Carlo Analysis
Given the faintness of the emission we are looking for, we have applied an ad-hoc complicated preparation of our images before the profile extraction. Although this effort suggests a tentative detection, our results could still be dominated by systematic errors. In particular, our measurement of the CGM radial profile out to very large radii could be affected by 1) errors in the empirical PSF subtraction of the quasar profile using the PSF of stars, including errors in the re-scaling of the stellar PSF; and 2) systematics in the data reduction and in the data itself (e.g., flat-fielding, constant sky-subtraction, ccd-noise, masking). These systematics could prevent the noise from averaging down according to simple photon-counting statistics, i.e. the expected error on the average SB should be inversely proportional to the square root of the area employed in the average .
For these reasons, to firmly assess the statistical significance of our results, we perform a Monte Carlo analysis. We have built samples of 15 stars, by randomly selecting objects from the 115 stars that were used to determine the PSF of the NB images (see §3.3). For each sample of stars we have then computed the average radial profile, and subtracted from it the normalized average profile of the global sample of the stars. In this way, we are left with realizations of our same experiment, i.e. QSOs’ profile stars’ profile. These profiles should be consistent with zero at each radius, and the values, i.e. the ratio between the signal and the 1 error within each bin, should be consistent with a Gaussian distribution of unit variance. If however we find Gaussian distributions with variances greater than 1, this would give us information on the degree to which we are underestimating the noise, and/or at which radius systematics begin to dominate the error budget. Furthermore, these Monte Carlo simulations of the null experiment will allow us to robustly estimate the statistical significance of our detection in a way that also folds in both the formal statistical errors, as well as the systematic errors.
In Figure 9 we show the distribution of the values for these 1000 realizations in each radial bin, starting from the kpc radius, and compare them to our results. We focus on the histograms for the larger radii because, as discussed in the previous section, the smaller scales can be contaminated by EELR or host galaxy emission which are not the focus of the current study. It is clear that the histograms have an approximately Gaussian shape, well centered on zero (blue vertical lines). However the variances increase at large radii, indicating that the error budget is deviating from the expectation based on photon counting statistics. Indeed, for kpc, we only slightly underestimate the noise (by %), while at larger radii we incur significantly larger errors than photon counting would suggest. This is likely due to the fact that bins at larger radii, i.e. kpc, probe larger area and the implied SB limit is significantly lower. As our SB limit decreases we become more sensitive to subtle systematics in the data, e.g. correlations between neighboring pixels, which are larger than the formal photon counting error, and which may not average down as . As we noted at the beginning of this section, this deviation can be ascribed to several possible effects.
The brown vertical lines in Figure 9 indicate the measurements obtained in the previous section for comparison (Figure 7). Note that all our measurements constitute a positive fluctuation which is not expected for a spurious signal151515This statement assumes that the measurements in each radial bin are independent, which is a reasonable assumption provided that there is not a large scale systematic that correlates multiple bins.. In other words, the Monte Carlo realizations are well centered around zero, while our measurements always exhibit positive values. Given the current binning, with the exception of the bin at 362 kpc, we have a fluctuation in every bin (or in other words above the 84th percentile of the Monte-Carlo realizations).
This can be better visualized in Figure 10, where we highlight the area within the 16th and 84th percentile of the 1000 realizations in comparison to the difference profile between QSOs and stars, e.g. same profile of Figure 7. Note that in agreement with the previous histograms (Fig. 9), the 16th percentile of the 1000 realizations is negative. This can be seen in the log-lin part of Figure 10. Note however that our Monte-Carlo realizations are not well centered around zero following the uncertainties expected in some of the annuli, i.e. the error on the mean should be given by and thus fluctuations of our mean values about zero should be comparable to this. In other words, we stress that the presence of systematics also appears to give rise to larger fluctuations about zero in our Monte-Carlo realizations than simple Gaussian statistics would suggest. However, it is reassuring that these fluctuations are nevertheless both positive and negative and do not show up in all annuli.
Further, although the statistical significance in a single radial bin is not high, we can combine our data to increase the signal-to-noise ratio of our detection. We decide to consider the data in one large radial bin spanning from kpc kpc, and conduct a similar Monte Carlo analysis. Figure 11 shows the histogram of the 1000 realizations in comparison to our measurement (red vertical line) of the stacked PSF subtracted emission (i.e. QSO-star) for this large radial bin. As expected from the previous analysis, the distribution is clearly Gaussian and centered at zero161616Note that in this case the mean value of the histogram of the Monte Carlo realizations is well centered on zero and consistent with the error on the mean., and we somewhat underestimate the error. Specifically, given that we find , our error is underestimated by 30%. Taking this into account, the deviation of our measurement () from the distribution would correspond to a detection in this wide bin, or in other words, it falls in the 98.5th percentile of the Monte-Carlo realizations.
This detection corresponds to SB erg s cm arcsec within an annular region extending from kpc to kpc around the QSO. In the next section we will discuss the implications of this measurement in the context of our current understanding of the gas distribution around QSOs, both from observations and from simulations.
As discussed in detail in the literature, given the resonant nature of the Ly line one needs to be particularly cautious when interpreting observations of this emission (e.g., Gould1996; Cantalupo2005; Dijkstra2006; Verhamme2006). Further, as summarized in the introduction, multiple physical processes can potentially power extended Ly emission and disentangling them is challenging, especially given that several of them could be at work. These mechanisms are (i) resonant scattering of Ly line photons (Dijkstra2008; Steidel2011; Hayes2011), (ii) photoionization by stars or by a central AGN (Geach2009; Overzier2013), (iii) cooling radiation from cold-mode accretion (e.g., Fardal2001; Yang2006; Faucher2010; Rosdahl12), and (iv) shock-heated gas by galactic superwinds (Taniguchi2001; Dijkstra2012).
In the following discussion we assume that the dominant mechanism powering the faint Ly emission we observed around quasars is photoionization from the central AGN, resulting in a boosted Ly fluorescence signal from the surrounding gas distribution. Indeed, scattering of Ly photons from the quasar should not be a major contributor on such large scales ( kpc, Cantalupo2014). This is because the resonant scattering process results in very efficient diffusion in velocity space, such that the vast majority of resonantly scattered photons produced by the quasar itself (or by the central galaxy) escape the system at very small scales ( kpc), and hence do not propagate to larger distances (e.g., Dijkstra2006; Verhamme2006; Cantalupo2005; FAB2015). Further, in the case of Ly emission in a superwind we would expect much higher values than observed for the Ly surface brightness, i.e. at least SB erg s cm arcsec (Taniguchi2001), depending on shock velocity. However, we cannot rule out a cooling radiation scenario, and indeed, it is interesting to note that the SB that we measure seems consistent with the SB values shown for massive halos at large radii ( kpc) by Rosdahl12, even though their simulations were run to interpret the LAB phenomenon and do not include the presence of an active AGN. Despite this, they predict much higher emission (SB erg s cm arcsec) at smaller scales ( kpc) than we observe. Modeling the expected Ly SB from cooling radiation clearly requires further detailed investigation.
We thus interpret the extended emission that we measured as AGN boosted Ly fluorescence, and adopt two different approaches for our analysis. First, we consider the observed SB in the context of the simple model for cool halo gas described in Hennawi2013 (already used in FAB2015), combined with recent observational constraints on the cool gas distribution around QSOs deduced from absorption line studies. Second, we compare our emission measurement to the results of a cosmological hydrodynamical simulation post-processed with ionizing and Ly radiative transfer, which was run to interpret the UM 287 nebula (see Cantalupo2014).
By studying absorption features in the spectra of projected quasar pairs, the “Quasars Probing Quasars” series of papers (QPQ; Hennawi2006 and references thereafter) has shown that the optically thick absorbers arising from gas clouds within the halo of quasars are anisotropically distributed with respect to the quasar itself (Hennawi2007). In other words, a quasar emits radiation anisotropically or intermittently, so that the number of optically thick absorption systems is higher in the background sightlines (i.e. at a certain projected distance from the quasar) in comparison to that observed along the line-of-sight to that quasar. This result, together with the absence of fluorescent emission from these optically thick absorbers (Hennawi2013) as well as photoionization modeling of individual absorbers (Prochaska2009; QPQ8) all indicate a scenario in which the optically thick gas observed in absorption toward background sightlines is shadowed from the ionizing radiation of the quasar. Prochaska2013b have shown that these optically thick absorbers have a high covering factor on kpc in the quasar’s CGM (see also Hennawi2007; Prochaska2013; Hennawi2013).
Hennawi2013 argued that, given this high covering factor and the luminosity of the quasar, the optically thick clouds observed in absorption would result in a very bright signal if directly illuminated by the ionizing radiation of the quasar. We can estimate this SB for our sample assuming the simple model in Hennawi2013, where a quasar is surrounded by a spherical halo populated with spherical clouds of cool gas ( K) at a single uniform hydrogen volume density , which are uniformly distributed throughout the halo, with an aggregate covering factor of . By assuming this model, the Ly surface brightness for an optically thick distribution of clouds is given by (Hennawi2013)
where we used the redshift of our observations (), the virial radius of the dark matter halo hosting quasar as distance ( kpc), and the specific luminosity at the Lyman limit () corresponding to the average magnitude of our sample, i.e. -mag . This value of is obtained as in FAB2015 by integrating the Lusso2015 composite quasar spectrum against the SDSS filter curve, and choosing the normalization which gives the correct magnitude for each quasar in our sample. This signal is almost three orders of magnitude larger than the observed SB value of the Ly emission at 160 kpc from our stacking analysis, i.e. SB erg s cm arcsec, and would have been easily detected even in our individual exposures. Along the lines of the argument in Hennawi2013, we conclude that the low SB signal that we observed cannot be arising from optically thick gas clouds illuminated by the quasar.
On the other hand, if we assume that gas with physical properties similar to the gas probed by absorption studies is uniformly distributed throughout the halo, the quasar would shine on the remaining component, keeping it highly ionized, i.e. optically thin. Following Hennawi2013, we can thus show that such optically thin clouds would result in a Ly SB level comparable to our observations
where we used the redshift of our observations (), the median value for the column density found in absorption studies (QPQ8) and an expected value for the volume density from simulation on such scales (e.g. 160 kpc; Rosdahl12).
Following the foregoing discussion, it is thus plausible that the low Ly signal that we detect in our observations is due to optically thin gas illuminated by the quasars. Indeed, the only way to match our observations via emission from optically thick clouds would be to require a very low covering factor in eqn. (5), i.e. , but this would be clearly in conflict with the results in Prochaska2013b, arguing for the optically thin scenario. Note that optically thin emission is also in agreement with our previous work on giant Ly nebulae. Indeed, although the quasars in our sample are on average fainter (average -mag=18.64) than the UM 287 quasar (-mag=17.28, Cantalupo2014), Jackpot show that the quasar SDSSJ0841+3921 (-mag=19.35), which is fainter than our average quasar, is able to keep the gas highly ionized and optically thin out to scales of a few hundred kpc. Specifically, SDSSJ0841+3921 is surrounded by a large scale ( kpc) Ly nebulosity which is as bright (SB erg s cm arcsec) 171717This is the average SB within the 2 S/N isophote, after removing the contribution from the unresolved emission of the 4 AGN. as that surrounding the UM 287 quasar (SB erg s cm arcsec) 181818This is the average SB within the same S/N isophote as for SDSSJ0841+3921, after removing the contribution from the unresolved emission of the 2 AGN.. Similar to UM 287, this Ly nebula is believed to result from fluorescent emission powered by the ionizing radiation of the QSO. Thus, the comparable SB in these two distinct cases, together with the factor of difference in the quasar ionizing luminosity191919The ionizing luminosities are and 30.4, respectively for UM 287 and SDSSJ0841+3921 (FAB2015; Hennawi2013)., strongly suggests that brightness of the central source does not play an important role, provided it is bright enough to ionize the surrounding gas distribution. This again favors an optically thin emission scenario, since in this regime eqn. (5) shows that the Ly SB scales as SB, which is independent of the luminosity of the quasar.
Having established that the optically thin regime is a good hypothesis, we can now derive the physical properties of the emitting gas. To break the degeneracy between the covering factor , column density , and volume density intrinsic of our SB measurement in the optically thin regime, we need independent constraints. Regarding , it has been shown that the smooth morphology of the emission in giant Ly nebulae implies a covering factor of (FAB2015; FAB2015LAB; Jackpot). For simplicity, in what follows, we will always assume a covering factor of . Note that, as said before we are assuming , and thus this value is also well motivated by the distribution of optically thick absorbers, which show a high covering factor in the quasar’s CGM (Prochaska2013b).
Further, the total hydrogen column density can be constrained via photoionization modeling of absorption systems along the sightlines of background QSOs that pierce through the halo of a foreground QSO (QPQ8 and references therein). In particular, QPQ8 compared photoionization models to a statistical sample of absorbers in the CGM of typical QSOs, finding a median within kpc from the quasars, although they observed a substantial scatter in the distribution of values, i.e. dex (see Fig. 15 of QPQ8). A similar value (at impact parameter of kpc) was estimated from photoionization modeling of SDSSJ0841+3921, exploiting the so far unique case where a giant Ly nebula could be modeled in both emission and absorption (Jackpot). Given that these are the most recent and reliable estimates for in the literature which are also consistent with each other, we henceforth assume a value of throughout our analysis.
Thus plugging in the values , , and into eqn. (5), for the Ly SB in the optically thin scenario, we find the typical volume density expected on scales of kpc 202020 kpc is the average distance of the bin used in our analysis, i.e. kpc kpc., given that SB erg s cm arcsec from our measurement. We thus obtain an estimate of the volume density of the gas illuminated by the quasar and fluorescing in Ly:
Given this estimate for at kpc, we can compute the ionization parameter at these distances around quasars. If we assume that is not affected by illumination on these large scales, its value holds for both the shadowed and illuminated gas, and varies due to the different ionizing sources, either the ultraviolet background (UVB; e.g., HM2012) if the gas is shadowed from the quasar emission, or the quasar if it is illuminated. We find a value of and for shadowed gas and illuminated gas, respectively. For randomly observed absorption systems, we thus expect to see a strongly bi-modal distribution of values. Hence, it is interesting to compare this result with what was found by QPQ8 from photoionization modeling of individual absorbers. Indeed photoionization modeling constrains the ionization parameter , as long as ionic ratios can be measured (e.g. C/C, Si/Si). Note that the average ionizing luminosity of their sample is , i.e. 0.5 dex lower than our sample, resulting in .
In their data there is no sign for a bi-modal distribution in (see Figure 7 in QPQ8), with low ionization parameter at any distance from the quasar and no systems at higher , implying that these absorbers are likely not illuminated by the quasar. It is thus plausible that the Lau et al. sample was biased against sources with very high values, which would result in fewer detectable metal absorption lines needed to constrain the ionization parameter through photoionization models. Indeed, two systems out of 14 studied by QPQ8, i.e. J0409-0411 and J0341+0000, do not have metal line detections, and thus do not have constraints on the ionization state of the gas, i.e. no estimate on . Nevertheless, the situation is more complicated with effects that can wash out the bi-modal distribution expected from the simple picture of half shadowed and half illuminated gas. In particular, star-forming galaxies in the vicinity of the absorber could boost the value of shadowed sources from to a higher . Further, the real distance between QSO and absorber is unknown, i.e. we have information only on the projected distance on the sky. If the real distance is 2-3 times larger than the projected distance, this would result in a reduction of 4-10 in , and hence in . In addition, if the highly ionized gas has similar velocity as the less ionized gas, its absorption features would then be lost in the stronger absorption of the less ionized gas, resulting in a bias against such high-U absorption systems. Clearly, to confirm this scenario and test the relevance of these additional effects, further observations and analysis are needed.
However, given that we have an estimate for , it is then interesting to understand which kind of absorption features are expected from the gas that we detect here in emission and compare them to the two systems in QPQ8, which do not have metal line detections. To estimate this, we use the photoionization code Cloudy (Ferland2013) to calculate the column densities expected for different ions for an optically thin slab of gas at the distance of kpc from the quasar. As input parameters we used the aforementioned quantities that would give rise to the observed Ly emission, i.e. cm, . We have then included the presence of the ultraviolet background (UVB; HM2012), and assumed the gas to have a metallicity of Z, similar to values in QPQ8. This calculation shows that in the Ly emitting gas cm, cm, cm, cm, and both and have very low values ( cm), confirming that would be extremely difficult to detect metal absorption from such systems. Indeed, the two highly ionized systems reported in QPQ8 have column densities comparable to these expectations. In particular, J0409-0411 shows cm, cm, cm, and has an impact parameter of 190 kpc from the quasar, while J0341+0000 similarly shows cm, cm, cm, cm, and has an impact parameter of 235 kpc. Thus, our analysis suggests that these two systems in QPQ8 are good candidates for background sightlines piercing gas which is illuminated by the quasar, and hence optically thin, and characterized by a high ionization parameter.
Lastly, we can compare our empirically derived value with simulations. We are not aware of an average density profile for the cool gas in massive (M M) halos in the literature, although we plan to determine it in future work. Simulations of massive halos (e.g., Rosdahl12) usually predict values in the range ( cm) on these scales ( kpc) within filamentary structures. These values are in rough agreement with our estimate. However, they may not optimally resolve this gas (Cantalupo2014; Crighton2015), and high densities are also probably present (FAB2015; QPQ8).
Cantalupo2014 showed that one can calibrate relations between SB and the total hydrogen column density in cosmological simulations post-processed with radiative transfer using the RADAMESH Adaptive Mesh Refinement code (Cantalupo2011), resulting in a picture which is consistent with the simple analytical relations discussed above, i.e. eqn. (5) and eqn. (5). By comparing pixel by pixel the mock image obtained for the SB and the column density map of cool gas ( K) in the case where the gas is mostly ionized by the quasar, i.e. what we above referred to as the optically thin scenario, Cantalupo2014 obtained the following relation212121Note that in the highly ionized case is basically . from their simulation:
This relation depends on the clumping factor , introduced by Cantalupo2014 to account for dense gas on small scales that are possibly unresolved by the simulation. If we plug in our value of SB erg s cm arcsec, and assume , we find , which is consistent with the simulation values (see Extended Data Fig. 3 of Cantalupo2014), and in agreement, within the uncertainties, to the observational value we used in the foregoing discussion . We thus find agreement between simulations and our simple analytical model in the case of these low SB levels, as opposed to the case of the giant bright Ly nebulae, where very high clumping factors are required by simulations to match the observed SB, with values up to (Cantalupo2014).
In summary, we argue that the hydrogen volume density around typical QSOs on kpc scales is approximately constrained to be cm, in broad agreement with the cool phase densities found in cosmological simulations 222222Intriguingly this value is in agreement with a radiation pressure confinement scenario (RPC; Stern2014), shown to be in place at least up to 10 kpc from an AGN. In this scenario the radiation pressure is the main contributor to the gas pressure so that , and the density profile is then shown to follow a simple relation. Following eqn. (10) in Stern2014, using kpc, the average ionizing photon flux for our sample s, erg , and , we obtain cm, remarkably close to what we found. However, to establish a RPC scenario on such large scales one would need a particular tuning of the time scales as the photons need to reach 300 kpc from the source. Nonetheless, it is worth noticing that any additional source of external pressure that would break the RPC scenario, would make the gas density higher (Stern2014).. However, bright giant Ly nebulae, like UM 287 (Cantalupo2014) and SDSSJ0841+3921 (Jackpot), appear to require much higher densities cm (Cantalupo2014; FAB2015; Jackpot), to explain their higher SB Ly emission. Thus our analysis indicates a scenario in which % of quasars show bright giant Ly nebulae (SB erg s cm arcsec on scales kpc) which require high gas densities, while the remainder of the population is characterized by fainter extended Ly emission and thus by much lower gas densities. At the moment, it is not clear what physical process is responsible for these factor of difference in gas density. It has been argued that the presence of a giant Ly nebula may be physically connected to the location of over-densities of galaxies and AGN which provides an important clue (Jackpot).
5.1. Comparison with Previous Deep Observations around Quasars
Although many studies have targeted extended emission around individual quasars (HuCowie1987; Heckman1991spec; Heckman1991; Christensen2006; Hennawi2009; North2012; Hennawi2013; Roche2014; Herenz2015), they do not achieve the SB sensitivity indicated by our stacked profile of the quasar CGM. This is primarily due to the limited sensitivity of the instruments employed in these studies or because they had a different chief science goal, e.g. focused on studying the much brighter emission from the EELR on smaller scales (Christensen2006).
On the other hand, both Rauch2008 and Cantalupo2012 have targeted quasars down to interesting sensitivity levels that should, according to our average profile, allow one to detect emission on 100 kpc scales. In particular, Rauch2008 achieve a SB limit of SB erg s cm arcsec ( in 1 arcsec) in their 92 h long-slit observations of the quasar DMS 2139-0405 (; Hall1996), while Cantalupo2012 reach SB erg s cm arcsec ( in 1 arcsec) with 20 h of narrow-band imaging targeting the quasar HE0109-3518. Although these studies reached interesting depths, they do not appear to show evidence for any extended Ly emission on large scales. However, as we have shown, probing these extremely low surface brightness levels on 100 kpc scales requires a very careful analysis with proper accounting of systematics (e.g., flat-fielding, sky and continuum subtraction, contamination from nearby sources) and a specifically tailored analysis. It would be interesting to search for the extremely faint CGM Ly emission that we detected in our stack in these very deep observations of individual quasars, which we plan to pursue in future work.
5.2. Comparison with LBGs and LAEs Ly Profiles
With the aim of comparing our newly determined Ly profile with previous observations at comparable redshift, and to gain a better overview on extended Ly halos, we compile extended Ly profiles for different sources from the literature. The left panel of Figure 12 shows our data together with the average profile of 92 continuum-selected LBGs at (dashed magenta line; Steidel2011), the average profile of 24 proto-cluster LBGs at (brown squares; Matsuda2012), the average profile of 130 LAEs at , selected to be in dense regions (blue squares; Matsuda2012), and the median profile of 27 faint line emitters at (cyan triangles; Rauch2008).
The origin of Ly halos around these star-forming galaxies is still a matter of debate, and, as for the case of giant Ly nebulae around quasar, several mechanisms have been discussed in the literature. First, it is important to point out that Ly cooling radiation is not able to account for the high Ly luminosity of these halos, i.e. erg s, giving rise to halos which are an order of magnitude fainter than observed (e.g., Haiman2000; Dijkstra&Loeb2009; Dijkstra2012, but see also Rosdahl12). However, an additional boost to the Ly emission could result from the ionizing photons emitted by the galaxy itself or by nearby galaxies (e.g., Dijkstra&Loeb2009; MasRibas2016). The photoionized gas could then emit enough Ly photons to be consistent with the observations. This would depend, however, on the ionizing photon escape fraction. Further, resonant scattering is not able to reproduce the extent of the Ly halo, giving rise to more compact emission than observed (e.g., Dijkstra2012; Lake2015), unless bipolar outflows in which the clumps decelerate (in the case of LBGs; Dijkstra2012), or additional sources of Ly emission (in the case of LAEs; Lake2015), such as low star-formation spatially distributed inside the host dark matter halo or cooling radiation, are taken into account. In summary, current studies suggest that in the case of LAEs and LBGs, as it is the case for QSO, a strong contribution from photoionization is needed in order to explain the extent and surface brightness of the Ly halo observed, with up to 50-60% of the observed Ly signal due to this powering mechanism (MasRibas2016).
In the right panel of Figure 12, we show the Ly profile of the UM 287 nebula (red; Cantalupo2014) and of the Ly nebula around SDSSJ0841+3921 (blue; Jackpot) that we have computed from the continuum-subtracted NB image by averaging in circular apertures around the QSO after masking all the sources from the -band (i.e. available broad band) and compact sources from the NB, as done for our individual objects in the GMOS observations (see §3.1). Note that the data reduction of UM 287 and SDSSJ0841+3921 has been performed differently. First of all, the Ly image has been obtained through a continuum-subtraction (see Cantalupo2014; Jackpot), and no empirical PSF-subtraction has been performed on these data. Further, no detailed variance images have been computed. For these reasons, the errors shown in this plot for UM 287 and SDSSJ0841+3921 do not include the errors on the continuum and PSF subtraction, but are instead estimated through a bootstrap analysis (). Note that these profiles, as highlighted in our analysis, could be affected by systematics at large radii (i.e. kpc), where imperfections in the data reduction and in the data itself may be dominant (e.g., flat-fielding, sky-subtraction, continuum subtraction). In both panels we indicate with a vertical line the expected virial radius for a QSO at this redshift ( kpc, Prochaska2014). As stressed already in our previous work, the extent of the Ly emission goes beyond the expected size of the dark matter halo hosting quasars.
Figure 12 shows that QSOs and LBGs have a higher emission profile which extends further out than LAEs. This could reflect the physics behind the Ly emission, with objects able to ionize a larger amount of gas out to larger distances and characterized by a denser environment surrounded by higher and more extended Ly signal. With the amount and distribution (density) of gas probably determining the strength of the Ly emission (i.e. optically thin scenario), and with the environment probably determining the extent (Matsuda2012). In other words, this plot seems to suggest that more massive systems, i.e. the quasar population, show higher Ly surface brightness profiles because the stronger ionizing radiation from the central source (or from nearby galaxies)232323Note that seems unlikely that unresolved/undetected galaxies could have a large contribution to the Ly emission that we detect. Indeed, if we use the flux limit of erg s cm (above which we are surely complete), we get a crude estimate of the star-formation rate SFR M yr for undetected galaxies in our NB data, using the formula in Kennicutt1998 and case B. Satellites with SFR of this order seem to have a small impact on the brightness and morphology of extended Ly emission on scales of kpc (Cen2013). is able to ionize the gas at much larger distances than in less massive systems. Indeed, ‘average’ LBGs and LAEs rely only on the ionizing radiation of stars (i.e. no bright QSO or AGN) and populate lower mass halos, i.e. M and M for LBGs and LAEs, respectively. Thus, probably for these reasons, their Ly profile is lower and less extended.
On the other hand, systems like UM 287 and SDSSJ0841+3921 probably represent rare specific cases (i.e. % of quasars) where the gas densities are simply much higher. However, it is interesting to note that at kpc, both our stacked profile and the profiles for UM 287 and SDSSJ0841+3921 show a Ly signal of SB erg s cm arcsec. It is thus plausible that our stacking analysis is consistent with (or detects) the same gas in the CGM or filamentary structures on these large scales around the quasars as for UM 287 and SDSSJ0841+3921 (Cantalupo2014; FAB2015; Jackpot). However, deep observations of a large sample of QSOs are needed to firmly characterize their extended Ly profile with the same precision as for the LBGs and LAEs, and conclusions could be drawn only after disentangling the several mechanisms that can produce low SB Ly emission.
6. Summary and Conclusions
Using NB imaging data collected as part of the FLASHLIGHT-GMOS survey (see also FAB2013), we have performed a stacking analysis to characterize the extended Ly emission around typical bright QSOs. We find that:
The 15 QSOs in our sample do not show giant Ly nebulae similar to UM 287 (Cantalupo2014) or SDSSJ0841+3921 (Jackpot), i.e. SB erg s cm arcsec emission on scales kpc, even though we would have been able to detect such emission. The PSF subtraction reveals Ly emission above SB erg s cm arcsec extending on scales of radius arcsec ( kpc) around 7 out of 15 quasars, i.e. 47% of the sample.
The average radial profile of the 15 QSOs in our sample shows a deviation from the Moffat PSF of our NB images, starting at kpc at around SB erg s cm arcsec. This can be translated to a low significance first radial profile of the Ly emission of the quasar CGM (see Figure 7). Using a Monte Carlo analysis (see §4.1), we ascertain that we have a detection within an annular bin spanning kpc kpc from the QSOs. The Ly emission in this bin, centered at kpc, is estimated to be SB erg s cm arcsec.
Assuming a scenario in which the illuminated gas is highly ionized by the quasar radiation, the detected SB on scales of hundreds of kpc implies gas volume densities of the order cm. This value is much lower than what has been proposed in the case of extended bright Ly nebulae around UM 287 (FAB2015) and SDSSJ0841+3921 (Jackpot), i.e. cm. However, our stacking analysis results in a Ly profile consistent with the Ly profile of UM 287 and SDSSJ0841+3921 on larger scales of kpc (SB erg s cm arcsec). Thus on these larger scales, it plausible that our stacking analysis detect a signal from the CGM or IGM as has been observed around UM 287 and SDSSJ0841+3921.
Future surveys targeting the Ly line around QSOs with new IFU instruments, i.e. MUSE (Bacon2010) and KCWI (Morrissey2012), will be able to investigate the CGM of quasars without the tight restrictions on redshift range and redshift accuracy inherent in narrow-band imaging. This will allow one to observe a much broader range of quasars, enabling studies of quasars with specific attributes such as hyper-luminous quasars, obscured quasars, radio-loud quasars, or quasars in richer environments. Such studies are needed to test our current interpretation and to determine if specific types of quasar preferentially exhibit bright emission from their CGM.
Further, the lack of bright giant Ly nebulae in our data confirm the necessity of large statistical samples of QSOs to uncover the brightest Ly nebulae on the sky. The study of these bright nebulae, such as the UM 287 nebula (Cantalupo2014) or SDSSJ0841+3921 (Jackpot), enable further characterization of the CGM gas in individual systems via the detection and modeling of additional emission lines, e.g. C iv, He ii (FAB2015LAB; FAB2015). The resulting physical properties can then be used as additional constraints on simulations of galaxy-formation. Indeed, we call attention to the lack of high-resolution cosmological zoom simulations of very massive systems ( M) to which one could compare our results (but see Fumagalli2014 and Faucher-Giguere2016 for recent progress on this front). Given the complexity of the radiative transfer of the resonant Ly line, simulations with both ionizing and Ly radiative transfer of massive halos are needed to fully capture the physics in play in giant Ly nebulae, and to compare to the wealth of observational data that are now becoming available.
Appendix A Estimate of the Flux Losses due to the Error on the Redshifts
As already stated in Section §2, the typical uncertainty on the redshift of our sample is (or equivalently km s), which is much smaller than the width of the narrow-band filter used, i.e. (or equivalently km s). Further, to be sure that the Ly emission of our targets fall within the narrow-band filter, we selected only QSOs whose redshift gives a maximum shift of Å from the filter’s center (or equivalently km s).
Nevertheless here, we test the robustness of our results against errors on the systemic redshift of the quasars in our sample. In particular, we estimate the flux losses due to this uncertainty. We proceed as follows. First, starting from the known systemic redshift of our sample, we build 1000 realizations of our 15 quasars by assuming that the error on the redshift is described by a gaussian distribution with km s, km s, km s, or km s. Second, we assume that the signal we are looking for, i.e. the extended Ly emission on scale of the CGM, is characterize by a km s centered at the systemic redshift of the quasar, as in the case of the giant nebula around UM 287 (FAB2015). Knowing the shape of our narrow-band filter, we are then able to compute the amount of flux lost for each sample of the 1000 realizations by taking the average of the loss per quasar within each sample.
Figure 13 shows the distribution of these average flux losses. This test shows that our sample should only be affected by a flux loss of the order of 3% (see blue histogram), while even doubling the error on the redshift would result in a 7% loss (green histogram). Only error distributions with width comparable to the width of our narrow-band filter can result in significant flux losses (see red and cyan histograms), and thus alter remarkably our analysis. Given the small uncertainty on the redshifts of our quasars, and our conservative selection criteria, we conclude that our sample do not suffer from flux losses.