Discovery of Optically Emitting Circumgalactic Nebulae around the Majority of UV-luminous Quasars at Intermediate Redshift

Surveys of the star formation rates and available interstellar gas reservoirs (Kennicutt & Evans 2012; Tacconi et al. 2013) as well as the metallicity distribution of stars (Hayden et al. 2015; Greener et al. 2021) in the Milky Way and distant, massive galaxies are inconsistent with simple "closed-box" models (Tinsley 1974) of galaxy evolution. These discrepancies can be resolved if galaxies accrete fresh material from external reservoirs as they evolve (for a review, see Putman 2017). At the same time, outflows driven by supernovae (e.g., Dekel & Silk 1986) and active galactic nuclei (AGNs; e.g., Silk & Rees 1998) are required in state-of-the-art cosmological simulations to reproduce the observed stellar masses of galaxies (e.g., Vogelsberger et al. 2014; Schaye et al. 2015). These same feedback mechanisms eject heavy elements produced by stars and supernovae out of the interstellar medium (ISM) and into intergalactic space to reproduce the mass–metallicity relation (e.g., Tremonti et al. 2004; Ma et al. 2016). Consequently, improving our empirical understanding of the rich set of physical processes that regulate gas exchange between galaxies and their surrounding cosmic ecosystems (for reviews, see Donahue & Voit 2022 and Faucher-Giguere & Oh 2023) represents a key priority of extragalactic astrophysics (National Academies of Sciences 2021).

Direct observations of the intergalactic and circumgalactic medium (IGM/CGM) are critical to developing a more complete understanding of galaxy evolution. However, the diffuse nature and cool–warm temperature range (T ≈10⁴–10⁶ K) of the IGM/CGM around typical galaxies make emission observations difficult with current facilities except in the very local Universe (e.g., Chynoweth et al. 2008; Lokhorst et al. 2022) and rare cases with unusually high surface brightness. Consequently, emission observations of the IGM/CGM often rely on deep, coadded observations averaging over hundreds to thousands of galaxies (e.g., Wisotzki et al. 2016; Dutta et al. 2023; Guo et al. 2023a, 2023b).

Much of our knowledge of the physical nature and extent of the CGM relies on sensitive absorption spectroscopy of UV-bright background sightlines that pass through the halos of foreground galaxies (for a review, see Tumlinson et al. 2017). These observations enable constraints on the densities, temperatures, and metallicities of the CGM around a wide variety of galaxies and across cosmic time (e.g., Rudie et al. 2019; Zahedy et al. 2019; Cooper et al. 2021; Zahedy et al. 2021; Qu et al. 2022). However, the lack of morphological information in absorption-line studies of the IGM/CGM necessitates nontrivial model assumptions when attributing the gas to inflows, outflows, or other origins (see Gauthier & Chen 2012; Ho et al. 2017; Schroetter et al. 2019; Zabl et al. 2019) except in rare cases where multiple sightlines enable velocity shear measurements (e.g., Chen et al. 2014; Zahedy et al. 2016; Lopez et al. 2018, 2020).

While direct emission studies of the CGM and IGM are not generally feasible for individual galaxies beyond the local Universe, IGM/CGM emission can be detected in some rare or extreme cases. In particular, systems like quasars that produce a significantly elevated ionizing radiation background increase the equilibrium recombination rate and temperature of the surrounding IGM/CGM, resulting in increased emission in recombination and collisionally excited lines (e.g., Chelouche et al. 2008). Consequently, observers have invested time on large ground-based telescopes toward identifying and following up systems where the IGM/CGM is expected to be observable in emission with state-of-the-art wide-field integral field spectrographs (IFSs; Bacon et al. 2010; Martin et al. 2010; Mateo et al. 2022). Despite cosmological surface-brightness dimming, the first successful surveys of individual CGM emission were conducted at z = 2–4 in H i Lyα emission due to the inherent strength of the line and its redshifting into the optical window. These surveys demonstrated that giant, ≈100 kpc CGM nebulae are nearly ubiquitous around luminous quasars at z = 2–4 (Borisova et al. 2016; Cai et al. 2019; O'Sullivan et al. 2020; Fossati et al. 2021; Mackenzie et al. 2021), and some extend to hundreds of kiloparsecs in projection (e.g., Cantalupo et al. 2014; Cai et al. 2018).

Recent IFS-enabled discoveries at lower redshifts of z < 1.5 revealed the potential for studies of ≈30–100 kpc-scale CGM nebulae observed in emission in nonresonant, rest-frame optical lines (e.g., [O ii], Hβ, and [O iii]) and resonant near-UV (NUV) lines (e.g., Mg ii) around quasars (Johnson et al. 2018; Helton et al. 2021; Johnson et al. 2022; Dutta et al. 2023; Epinat et al. 2024; Liu et al. 2024), galaxy groups (Epinat et al. 2018; Chen et al. 2019; Leclercq et al. 2022; Dutta et al. 2023; Epinat et al. 2024), and galaxies with evidence of recent bursts of star formation (Rupke et al. 2019; Burchett et al. 2021; Zabl et al. 2021). The joint morphological and kinematic analysis of the nebulae combined with deep galaxy surveys enabled by IFSs provide direct insights into the origins of the gas, including interactions, accretion, and outflows. Furthermore, statistical analysis of velocity structure functions around the largest of these nebulae provide new insights into halo-scale turbulence (M. C. Chen et al. 2023, 2024). However, while such giant, relatively high-surface-brightness, optically emitting nebulae have led to significant insights in individual case studies, their incidence rate and properties around quasars have yet to be quantified in statistical samples.

Here, we report the first IFS survey constraining the incidence of optically emitting circumgalactic nebulae around 30 UV-luminous, unobscured quasars at z = 0.4–1.4 using deep observations from the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) collected by the MUSE Quasar Blind Emitters Survey (MUSEQuBES; e.g., Dutta et al. 2024) and the Cosmic Ultraviolet Baryon Surveys (CUBS; Chen et al. 2020). The paper proceeds as follows. In Section 2, we describe the observations and data-reduction steps. In Section 3, we characterize the properties of the MUSEQuBES and CUBS quasars. In Section 4, we present the discovery of large nebulae in the quasar fields. In Section 5, we discuss the implication of our findings, search for correlations between the nebulae and quasar properties, and compare with circum-quasar nebulae at higher redshift. Finally, we summarize our results and make concluding remarks in Section 6.

Throughout, we adopt a flat Λ cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and a Hubble constant of H₀ = 70 km s⁻¹. Unless otherwise stated, all magnitudes are given in the AB system (Oke & Gunn 1983), all wavelengths are given in vacuum, and all distances are proper.

Both the CUBS and MUSEQuBES surveys were designed to study the CGM and IGM around galaxies at z ≈ 0.1–1.4 using deep and highly complete galaxy redshift surveys with MUSE in the fields of UV-bright background quasars with high-quality UV absorption spectra from the Cosmic Origins Spectrograph (Green et al. 2012) aboard the Hubble Space Telescope (HST). The broad spectral coverage (4700–9350 Å), wide field of view ($1^{\prime} \times 1^{\prime}$), and high throughput (peaking at 35% at 7000 Å) of MUSE (Bacon et al. 2010) collectively enable searches for low-surface-brightness nebulae in the environments of the CUBS and MUSEQuBES background quasars at z = 0.4–1.4 (Johnson et al. 2018, 2022; M. C. Chen et al. 2023; Liu et al. 2024) at the same time. The MUSEQuBES MUSE observations (PI: J. Schaye) were conducted in the MUSE wide-field mode under natural seeing conditions characterized by a FWHM = 05–10 with 2 to 10 hr of integration. The CUBS MUSE observations (PI: H.-W. Chen) have integration times ranging from 1.4 to 3.4 hr and were conducted in wide-field mode with ground-layer adaptive optics (Kolb et al. 2016; Madec et al. 2018). As a result, the CUBS data cubes exhibit improved an image quality of FWHM = 05–07.

We reduced the MUSEQuBES data with the MUSE Guaranteed Time Observations (GTO) team pipeline as described in Bacon et al. (2017), which includes standard reduction steps in the MUSE ESO pipeline (Weilbacher et al. 2020) but with postprocessing that improves image quality via self-calibration and principal-component-based sky subtraction (Soto et al. 2016). We also reduced the data using the CubEx postprocessing pipeline described in Borisova et al. (2016) and Cantalupo et al. (2019). For MUSEQuBES fields, we performed the analysis in this paper on both the GTO and CubEx reduction products and found consistent results in all cases.

We reduced the CUBS data cubes with the ESO pipeline both with and without CubEx postprocessing as described in Chen et al. (2020) and M. C. Chen et al. (2023). For all of the CUBS and MUSEQuBES fields, we also retrieved the reduced, coadded MUSE-DEEP data cubes from the ESO archive. Compared to the GTO and CubEx data products, MUSE-DEEP reductions exhibit more artifacts in continuum images and larger spectral residuals near sky lines. However, for the narrowband, continuum-subtracted data cubes and surface-brightness maps used in this work, we obtained similar results with the ESO-DEEP products except when near a night-sky emission line, in which case the CubEx and GTO reductions exhibit significantly improved residuals. After reducing the data, we converted the air wavelengths delivered by the MUSE pipelines to vacuum for convenience. The total MUSE exposure time available and seeing FWHM measured in the final stacked data cubes are listed for each quasar field in Table 1.

Table 1. Summary of Quasar Properties and MUSE Observations, Sorted by Quasar Redshift

	Quasar Properties									MUSE Data
Quasar Name	R.A. a	Decl. b	zc	Line d	$\mathrm{log}\lambda {L}_{\lambda }/\mathrm{erg}\,{{\rm{s}}}^{-1}$ e	$\mathrm{log}{L}_{\mathrm{bol}}/\mathrm{erg}\,{{\rm{s}}}^{-1}$ f	$\mathrm{log}{M}_{\mathrm{BH}}/{M}_{\odot }$ g	$\mathrm{log}{RL}$ h	Wr([O II]) i	${t}_{\exp }$ j	Seeing k
	(deg)	(deg)							(Å)	(hr)	(arcsec)
HE 0435−5304	69.2118	−52.9800	0.4279	[O ii]	44.9	45.8	8.4	<1.0	−0.78 ± 0.04	2.0	1.0
HE 0153−4520	28.8050	−45.1032	0.4532	Hβ	45.5	46.5	8.9	<0.9	−0.01 ± 0.01	2.0	0.7
HE 0226−4110	37.0635	−40.9541	0.4936	[O ii]	46.0	47.0	9.0	<0.1	−0.17 ± 0.01	8.0	0.5
PKS 0405−123	61.9518	−12.1935	0.5731	[O ii]	46.4	47.3	9.4	3.1	−0.36 ± 0.01	9.8	1.0
HE 0238−1904	40.1357	−18.8642	0.6282	[O ii]	46.2	47.2	9.8	<0.3	−0.34 ± 0.01	10.0	0.7
3C 57	30.4882	−11.5426	0.6718	[O ii]	45.8	46.8	8.9	3.7	−0.55 ± 0.02	2.0	0.7
PKS 0552−640	88.1020	−64.0363	0.6824	[O ii]	46.4	47.4	9.7	2.3	−0.32 ± 0.02	2.0	0.7
J0110−1648	17.6479	−16.8077	0.7822	[O ii]	46.2	46.9	9.1	2.4	−0.23 ± 0.02	1.4	0.6
J0454−6116	73.5664	−61.2740	0.7864	[O ii]	46.2	46.9	9.4	<0.9	−0.84 ± 0.02	1.4	0.7
J2135−5316	323.9716	−53.2821	0.8123	[O ii]	46.6	47.3	9.4	<0.6	−0.22 ± 0.01	1.9	0.5
J0119−2010	19.9837	−20.1729	0.8160	[O ii]	46.4	47.1	9.3	<0.9	−0.27 ± 0.01	1.4	0.5
HE 0246−4101	42.0261	−40.8093	0.8840	[O ii]	46.7	47.4	9.7	<0.7	−0.22 ± 0.03	2.1	0.7
J0028−3305	7.1266	−33.0970	0.8915	[O ii]	46.3	47.0	9.4	<1.0	−0.13 ± 0.01	1.4	0.5
HE 0419−5657	65.2246	−56.8455	0.9481	[O ii]	46.1	46.8	9.2	<1.2	−0.27 ± 0.05	1.4	0.6
Q0107−025	17.5677	−2.3141	0.9545	[O ii]	46.2	46.9	9.4	<1.1	−0.67 ± 0.10	2.0	0.9
Q0107−0235	17.5548	−2.3313	0.9574	[O ii]	46.1	46.8	9.6	3.2	−1.36 ± 0.05	2.0	0.8
PKS 2242−498	341.2508	−49.5301	1.0011	[O ii]	46.3	47.0	9.6	3.4	−0.83 ± 0.03	1.4	0.6
PKS 0355−483	59.3413	−48.2042	1.0128	[O ii]	46.3	47.0	9.6	2.9	−0.84 ± 0.02	2.6	0.5
HE 0112−4145	18.5921	−41.4964	1.0238	[O ii]	46.2	46.9	9.5	1.2	−2.04 ± 0.05	1.4	0.7
HE 0439−5254	70.0502	−52.8048	1.0530	Mg ii	46.3	47.0	9.3	<1.1	−0.03 ± 0.04	2.5	0.7
HE 2305−5315	347.1574	−52.9802	1.0733	Hβ	46.5	47.2	9.5	<1.0	−0.03 ± 0.02	1.9	0.6
HE 1003+0149	151.3968	+1.5793	1.0807	[O ii]	46.3	47.0	9.3	1.5	−0.62 ± 0.03	2.0	0.9
HE 0331−4112	53.2794	−41.0336	1.1153	[O iii]	46.8	47.5	9.9	<0.8	−0.11 ± 0.03	2.6	0.6
TXS 0206−048	32.3781	−4.6406	1.1317	[O ii]	46.4	47.2	9.8	3.0	−1.65 ± 0.03	8.0	0.6
Q1354+048	209.3594	+4.5948	1.2335	[O ii]	46.5	47.2	9.5	<1.1	−0.67 ± 0.04	2.0	0.5
J0154−0712	28.7278	−7.2061	1.2957	[O ii]	46.8	47.5	9.5	<0.9	−0.11 ± 0.01	3.4	0.5
Q1435−0134	219.4511	−1.7863	1.3117	[O ii]	47.1	47.8	10.0	1.8	−0.41 ± 0.02	6.1	0.6
PG 1522+101	231.1021	+9.9747	1.3302	[O ii]	47.0	47.7	10.1	<0.7	−0.18 ± 0.02	2.0	0.5
HE 2336−5540	354.8050	−55.3974	1.3531	[O ii]	47.1	47.8	10.0	2.5	−0.23 ± 0.01	2.6	0.6
PKS 0232−04	38.7805	−4.0348	1.4450	[O ii]	46.7	47.4	9.7	4.0	−0.70 ± 0.12	10.0	0.7

Notes.

^a Right ascension of the quasar in J2000 coordinates from the Gaia survey Data Release 3 (DR3; Gaia Collaboration et al. 2023). ^b Decl. of the quasar in J2000 coordinates from the Gaia survey DR3 (Gaia Collaboration et al. 2023). ^c Quasar systemic redshift measured as described in Section 3. ^d Emission line used to measure the quasar systemic redshift. Systematic redshift uncertainties correspond to ≈20, 50, and 130 km s⁻¹ for [O ii]-, [O iii]-, and Hβ/Mg ii-based redshifts, respectively. ^e Monochromatic luminosity measured at rest-frame λ = 5100 and 3000 Å for quasars at z < 0.7 and z > 0.7, respectively. Systematic uncertainties in the monochromatic luminosities are dominated by ≈10% flux-calibration errors. ^f Bolometric luminosity of the quasar estimated as described in Section 3. Systematic uncertainty in bolometric luminosities are dominated by population scatter in the bolometric corrections of ≈0.1 dex (Richards et al. 2006). ^g Supermassive black hole mass estimated with a single-epoch virial theorem approach described in Section 3. Systematic uncertainties in these black hole mass estimates are ≈0.4 dex (Shen et al. 2023). ^h Radio loudness based on flux at rest-frame wavelengths of 6 cm and 2500 Å. Systematic uncertainties in $\mathrm{log}{RL}$ are ≈0.1 dex. ⁱ Rest-frame [O ii] equivalent widths measured in the MUSE quasar spectra within the seeing disk corresponding to ≲5 kpc. ^j Total exposure time of the MUSE observations. ^k FWHM seeing measured in the MUSE white-light image of the quasar.

Download table as:  ASCIITypeset image

The ESO, GTO, and CubEx pipelines all share common wavelength and flux-calibration approaches. MUSE wavelength calibration frames are taken with three built-in lamps (HgCd, Ne, and Xe) by ESO Science Operations staff in the morning, after the science exposures are acquired. The observed locations of the arc lines on the detectors are fit with outlier-resistant two-dimensional polynomials to simultaneously solve for the dispersion, tilt, and curvature of the spectral lines in the arc frames. Comparisons of arc-lamp spectra taken over the course of a night and observations of telluric absorption features in MUSE spectra of bright stars indicate that systematic uncertainties in the wavelength calibration are typically less than 1 km s⁻¹ (Kamann et al. 2018; Weilbacher et al. 2020).

MUSE flux calibrations are taken with spectrophotometric standards nightly. Systematic uncertainties in the spectrophotometry as a function of position in the data cube and wavelength are typically less than 5% for data taken under photometric conditions (Weilbacher et al. 2020), though others have reported systematic errors in the flux calibration of up to 25% (e.g., Péroux et al. 2019). To quantify the level of systematic errors in the CUBS and MUSEQuBES flux calibration, we identified bright stars and galaxies in each MUSE data cube to serve as references and calculated synthetic photometry in the ACS+F814W and Dark Energy Camera i bands. We then compared these synthetic magnitudes with those measured in archival HST ACS+F814W images of the MUSEQuBES fields and public Dark Energy Survey (DES) Data Release 2 (Abbott et al. 2021) photometry for the CUBS fields. Based on comparisons of the synthetic MUSE photometry and reference photometry from the HST and DES, we estimate a 1σ uncertainty of 10%.

One of the primary differences in the ESO, GTO, and CubEx pipelines is their approach to illumination corrections. The ESO pipeline performs illumination corrections using a combination of flats taken with the calibration unit (Kelz et al. 2012) and twilight flats. However, subtle differences in the true illumination pattern during the science observations and those taken with the calibration unit and twilight sky correlate with IFU slice, resulting in artificial grid patterns in narrowband images and emission-line maps (see discussion in Bacon et al. 2015 and Borisova et al. 2016). To improve the illumination correction and remove these artificial patterns from MUSE data, both the GTO and CubEx pipelines perform "self-calibration" by producing illumination corrections based on the night-sky continuum and emission lines observed in science frames after masking continuum sources. These self-calibration steps effectively remove the grid pattern present in ESO cubes and reduce residuals in synthetic images formed from the data cubes by a factor of approximately 3 (Lofthouse et al. 2020).

The other major differences in the ESO, GTO, and CubEx pipeline products results from their night-sky-subtraction approaches. All three pipelines begin by masking continuum sources to identify spaxels where the detected flux is dominated by sky emission, but their modeling and subtraction of the night-sky emission differ significantly. The ESO sky-subtraction pipeline adopts a model of the night-sky line emission based on a fixed set of 5100 cataloged emission lines (Cosby et al. 2006) in the MUSE wavelength range. These are categorized into 52 groups, with lines in each group sharing the same excited state. The night-sky line-emission model fixes the relative intensity of lines in each of these groups and fits the observed global sky spectrum, allowing the total intensity of each group to vary and convolving with the line-spread function (LSF). After this, the best-fitting global sky model is subtracted and the residuals in source-free spaxels are used to estimate the continuous component of the night-sky emission for each slice (Streicher et al. 2011).

To reduce night-sky residuals present after the ESO sky subtraction, the GTO pipeline employs postprocessed sky subtraction with principal component analysis using the Zurich Atmosphere Purge package (Soto et al. 2016). On the other hand, the CubEx pipeline operates with postprocessing on ESO data products with the ESO sky-subtraction module turned off. It then uses the CubeSharp (Borisova et al. 2016; Cantalupo et al. 2019) routine to produce an empirical model of the LSF on a line-by-line basis (or groups of lines if they are close in wavelength). One key advantage of the CubEx sky-subtraction approach is that the LSF can be finely adjusted to improve sky subtraction while also conserving flux. Both the GTO and CubEx pipelines improve sky line residuals significantly, resulting in a factor of approximately 3 reduction in the standard deviation of flux in source-free spaxels (Lofthouse et al. 2020).

To study faint line emission near the bright quasars, we performed quasar light subtraction as described in Johnson et al. (2018) and Helton et al. (2021), which effectively removes the spatially unresolved continuum, broad-line emission, and narrow-line emission from the nucleus. Finally, to search for line-emitting nebulae, we produced continuum-subtracted emission-line maps around the expected wavelength of the [O ii] doublet at the redshift of each quasar. To do so, we first identified suitable continuum regions on the blue and red sides of the potential emission line. The continuum regions are typically 3000–4000 km s⁻¹ in width but tailored for each quasar based on proximity to night-sky emission lines. We then fit a low-order polynomial to the continuum regions in each spaxel of the data cube and subtracted the best-fit model to produce continuum-subtracted data cubes using the MUSE Python Data Analysis Framework (MPDAF; Bacon et al. 2016).

To characterize the properties of the CUBS and MUSEQuBES quasars, we extracted their optical MUSE spectra with MPDAF (Bacon et al. 2016; Piqueras et al. 2017) prior to quasar and continuum subtraction with apertures chosen to include >95% of the quasar light. We measured precise, narrow-line redshifts for the quasars using the [O ii] λ λ3727, 3729 doublet whenever possible, assuming an effective centroid of λ3728.6 expected for the 0.8:1 to 0.9:1 doublet ratio typical of quasar narrow-line emission (Hewett & Wild 2010). At the same time, we also measured the rest-frame equivalent width of the [O ii] doublet in the quasar spectra extracted within a seeing disk, W_r([O ii]). The seeing disks correspond to ≲5 kpc from the quasars, ensuring that the [O ii] equivalent widths capture emission from the narrow-line region rather than more extended nebulae. Quasar redshifts based on [O ii] emission are characterized by typical uncertainties of 20–30 km s⁻¹ (Boroson 2005; Hewett & Wild 2010). Four quasars do not exhibit detectable [O ii] despite the high signal-to-noise ratio (S/N) available in the MUSE spectra, necessitating the use of other emission lines observed in the near-IR by Sulentic et al. (2004). For HE 0331−4112, we estimated the systemic redshift using the [O iii] λ5008 emission-line peak. Quasar redshifts measured with the peak of the [O iii] line are characterized by typical systematic uncertainties of ≈50 km s⁻¹ when compared to stellar absorption-based redshifts (see discussion in Shen 2016). For the remaining three quasars, narrow emission lines are not detectable, necessitating the use of broad Hβ (HE 0153−4520 and HE 2305−5315) or Mg ii (HE 0439−5254). We estimated the systemic redshifts of these three quasars by cross-correlating the observed quasar spectra with the template from Hewett & Wild (2010) using methods described in that work. As discussed in Hewett & Wild (2010), the resulting broad Hβ and Mg ii emission redshifts exhibit uncertainties of ≈130 km s⁻¹ when compared to more precise narrow-line redshifts. The quasar redshifts and the emission lines driving them are listed in Table 1.

To estimate the luminosities and black hole masses of the quasars, we used PyQSOFit (Guo et al. 2018, 2019) to fit the quasar continuum with a power-law and Fe ii templates (Boroson & Green 1992; Vestergaard & Wilkes 2001) and multiple Gaussians for broad and narrow emission lines. To enable estimates of the bolometric luminosities and black hole masses of the quasars, we measured the monochromatic luminosity at 5100 Å and line width of broad Hβ for systems at z < 0.7, and the monochromatic luminosity at 3000 Å and line width of broad Mg ii for those at z > 0.7. For four CUBS quasars (HE 0331−4112, J2135−5316, PKS 2242−498, and HE 2305−5315), the Mg ii emission line falls near or in the MUSE sodium laser gap necessitated by the adaptive-optics system. For these, we obtained supplementary spectra with the Magellan Echellete (MagE) spectrograph (Marshall et al. 2008) to enable Mg ii line-width measurements (for details of the MagE observations and reductions, see Li et al. 2024). We then combined the monochromatic luminosities and line widths to infer supermassive black hole masses (M_BH) using the single-epoch virial theorem relations described in Shen et al. (2011). Typical uncertainties in single-epoch black hole mass estimates are 0.4 dex, though uncertainties may be higher for the CUBS and MUSEQuBES quasars because they are higher in luminosity than typical reverberation mapping quasar samples used to calibrate the relations (see discussion in Shen et al. 2023).

To enable comparisons of the luminosities of the quasars across redshifts, we converted the monochromatic luminosity measurements to bolometric luminosity estimates (L_bol) using bolometric corrections from Richards et al. (2006) of 10.3× and 5.6× for monochromatic luminosities at 5100 Å and 3000 Å, respectively. Uncertainties in the bolometric luminosities are dominated by ≈0.1 dex intrinsic scatter in bolometric corrections (Richards et al. 2006). Finally, we estimated the radio loudness, defined as the flux density ratio at rest-frame 6 cm versus 2500 Å, of each quasar using radio observations from the Rapid Australian Square Kilometre Array Pathfinder (ASKAP) Continuum Survey (RACS) Data Release 1 (Hale et al. 2021) as described in Li et al. (2024). In the case of nondetections, we calculated 3σ upper limits on the radio flux and converted these to upper limits on radio loudness. The RACS radio observations and MUSE rest-NUV observations are separated by a time period of several years, resulting in a nonnegligible systematic uncertainty of 0.1 dex due to potential NUV variability on these timescales (Welsh et al. 2011), though we caution that rare but significant radio variability could result in significantly larger errors in the radio loudness (Nyland et al. 2020). The monochromatic luminosities, bolometric luminosities, supermassive black hole mass estimates, radio-loudness measurements, and [O ii] rest-frame equivalent widths for each quasar are reported in Table 1.

To contextualize the properties of the CUBS and MUSEQuBES quasars, we plot their luminosities, black hole masses, and radio loudnesses versus redshift and compare with quasars in the Sloan Digital Sky Survey (SDSS) sample from Shen et al. (2011) in Figure 1. The CUBS and MUSEQuBES quasars were selected to be sufficiently UV bright for high-quality absorption spectroscopy with COS. As a result, the CUBS and MUSEQuBES quasars are among the most luminous quasars in the z < 1.4 Universe, with luminosities typically ≈1.5 dex higher than the median SDSS quasar at the same redshift. Similarly, the inferred SMBH masses of the CUBS and MUSEQuBES quasars are ≈1 dex above the SDSS median at similar redshifts (Shen et al. 2011). The 30 quasars are approximately evenly split between radio loud and radio quiet, with a division between the two classes of $\mathrm{log}{RL}\approx 1$.

Zoom In Zoom Out Reset image size

To identify large nebulae around the CUBS and MUSEQuBES quasars, we first visually inspected the quasar light and continuum-subtracted MUSE data cubes by searching for extended emission within ∣Δv∣ < 2000 km s⁻¹ of the expected wavelength of the [O ii] doublet. We chose this large velocity interval to be sufficient to include emission from bound gas in even massive galaxy clusters. Furthermore, the velocity interval is sufficient to account for potential systematic uncertainty in the quasar systemic redshifts, even for those with redshifts based on Mg ii or broad Hβ. We focused the search on the [O ii] doublet because the MUSE spectral range covers it for all of the CUBS and MUSEQuBES quasars. Moreover, detailed studies of individual nebulae around quasars at z ≈ 0.5–0.6, where a more comprehensive suite of emission lines is available, suggest that [O ii]-dominated regions of circum-quasar nebulae are often more extended than [O iii]-dominated areas (e.g., M. C. Chen et al. 2023; Liu et al. 2024). Our initial inspection revealed that most of the CUBS and MUSEQuBES quasars exhibit extended [O ii] emission on ≳20–30 kpc scales.

To enable quantitative detection and analysis of the nebulae, we performed thresholded detection and segmentation in three dimensions (3D; two spatial and one spectral) following the approach often used in studies of Lyα nebulae at z = 2–3 (e.g., Borisova et al. 2016; Arrigoni Battaia et al. 2019b; Sanderson et al. 2021). In particular, we formed a S/N data cube and then smoothed it in both the spatial and spectral dimensions using a Gaussian kernel with σ ≈ 1.5 pixels, though with a somewhat larger σ in cases with elevated noise due to night-sky lines. We then performed 3D source detection and segmentation on the smoothed S/N maps with a threshold of S/N ≈ 1.5, though with a slightly higher threshold in cases where the night-sky line residuals are elevated relative to the variance map. We then connected adjacent pixels above this threshold and required a minimum of 10 connected pixels and a total S/N ≳ 5 to define a detection. To uniformly visualize and characterize these nebulae, we produced [O ii] surface-brightness maps by integrating the unsmoothed flux in each spaxel over the wavelength interval defined by the 3D segmentation. For spaxels that do not include any detected nebulae, the background and noise properties correspond to those expected of three spectral pixels at the wavelength where the mean S/N per pixel of the nebula is highest. The resulting [O ii] surface-brightness maps for each quasar are shown in Figure 2.

Zoom In Zoom Out Reset image size

These 3D extractions have the advantage of visualizing the full extent of detected emission while reducing noise for kinematically quiescent or undetected regions compared to surface-brightness maps defined by fixed velocity intervals. However, they also exhibit complicated, spatially varying noise characteristics (see discussion in Borisova et al. 2016, Arrigoni Battaia et al. 2019a, and Mackenzie et al. 2021). To quantify systematic uncertainty in the measured properties of the nebulae, we also produced surface-brightness maps integrated over fixed velocity intervals chosen for each quasar field to include 95% of the detected [O ii] nebular emission in the quasar light and continuum-subtracted data cubes. The resulting narrowband surface-brightness maps are broadly consistent with those obtained from 3D segmentation. However, the S/N in the 3D segmentation-based maps is noticeably improved compared to synthetic narrowband images in cases where the nebular emission is spread over a wide velocity range. To ensure the robustness of the results, we performed all of the analysis in this paper on the surface-brightness maps produced by 3D segmentation and synthetic narrow bands. We confirmed that the conclusions of the paper hold with both approaches.

To characterize the size and flux of the nebulae around the CUBS and MUSEQuBES quasars relative to the seeing and noise level, we used the source measurement functions of Photutils (Bradley et al. 2016) to measure the FWHM of the nebulae in the surface-brightness maps. We integrated the 3D segmentation to estimate their total [O ii] flux. We define extended nebulae to be detections if their FWHM exceeds the FWHM seeing measured in the data cube and if their total S/N is greater than 5. With these criteria, all but three of the quasars have detected [O ii]-emitting, extended nebulae surrounding them. One of these, HE 2305−5315, does exhibit a small nebula offset from the quasar centroid, but its FWHM is comparable to the seeing, so we do not consider it further in this paper. In all 27 cases with detected, extended nebulae, the measured FWHM of the nebulae exceeds the seeing FWHM by a factor of >2.4. This ensures that the effect of seeing on the measured properties of the nebulae are minimal. For example, the contributions of the seeing to the measured FWHM of the nebulae are <10% when considered in quadrature.

To estimate the size of the nebulae, we measured the major axis length scale determined by the maximum projected separation between pixels contained within the S/N-defined segmentation map corresponding to each nebula. We adopted this definition of the size of the nebulae to enable comparison with studies of z ≈ 3 Lyα nebulae around quasars such as Borisova et al. (2016) and Mackenzie et al. (2021), which used a similar approach. In two cases, PKS 0405−123 (Johnson et al. 2018) and TXS 0206−048 (Johnson et al. 2022), multiple large nebulae are detected in the quasar environment, some of which are well offset from the quasars themselves. For these, we report the projected size for the largest contiguous nebula as defined in the 3D segmentation.

Several of the nebulae are quite elongated along one axis, meaning that the major axis size does not fully encapsulate the extent of the systems. For this reason, we also report the area, in square kiloparsecs, defined by the projection of the 3D segmentation of each nebula. In cases with multiple nebulae, we report the sum of their areas. Finally, to quantify the surface brightness of the nebulae, we computed the redshift-dimming-corrected surface brightness measured in an annulus centered on the quasar with an inner radius corresponding to 10 kpc and an outer radius corresponding to 30 kpc, which we denote as SB_10–30(1 + z)⁴. The outer radius of 30 kpc is chosen to be half of the median major axis full size of the nebulae around the CUBS and MUSEQuBES quasars. The inner radius is chosen to be large enough to avoid residuals from quasar light subtraction and to exceed the effective radii expected in stellar continuum images of late-type and early-type massive galaxies at intermediate redshift (e.g., van der Wel et al. 2014). The FWHM, major axis size, projected area, total S/N, total [O ii] luminosity, and SB_10–30(1 + z)⁴ measured for the nebulae are reported in Table 2. To quantify the systematic uncertainty in these quantities, we remeasured them using the synthetic narrowband images of the [O ii] nebulae. While the systematic differences inferred from the median ratios between the quantities measured with the 3D segmentation approach versus the synthetic narrowband images are minimal, there is a nonnegligible scatter corresponding to 10% for FWHM, 15% for size, 20% for area, and 25% for both luminosity and SB_10–30(1 + z)⁴.

Table 2. Summary of the [O ii]-emitting Nebulae around the CUBS and MUSEQuBES Quasars, Ordered by Redshift

Quasar	za	Flag b	FWHM c	Size d	Area e	S/N f	[O ii] Luminosity g	${\mathrm{SB}}_{10-30}{\left(1+z\right)}^{4}$ h
			(arcsec)	(kpc)	(kpc2)		(erg s−1)	($\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{arcsec}}^{-2}$)
HE 0435−5304	0.4279	=	2.6	55	1170	35.3	9.6 × 1041	5.6 × 10−17
HE 0153−4520	0.4532	<	0.6	7	35	2.4	1.9 × 1040	1.4 × 10−18
HE 0226−4110	0.4936	=	3.9	84	2740	156.3	1.1 × 1042	3.9 × 10−17
PKS 0405−123 i	0.5731	=	9.1	129	8180	134.8	2.7 × 1042	2.2 × 10−17
HE 0238−1904 j	0.6282	=	4.9	103	4820	187.9	2.3 × 1042	9.0 × 10−17
3C 57	0.6718	=	3.6	71	1930	71.3	1.2 × 1042	5.6 × 10−17
PKS 0552−640	0.6824	=	6.2	153	5180	119.6	3.4 × 1042	1.1 × 10−16
J0110−1648	0.7822	=	1.8	29	530	20.3	4.1 × 1041	1.3 × 10−17
J0454−6116 k	0.7864	=	3.5	102	4350	234.0	4.8 × 1042	1.3 × 10−16
J2135−5316 k	0.8123	=	3.8	83	2250	47.9	2.1 × 1042	9.4 × 10−17
J0119−2010	0.8160	=	4.2	96	3090	60.1	2.5 × 1042	1.4 × 10−16
HE 0246−4101	0.8840	=	3.8	74	2280	44.3	1.0 × 1042	4.6 × 10−17
J0028−3305	0.8915	=	2.1	42	500	23.2	2.9 × 1041	1.1 × 10−17
HE 0419−5657	0.9481	=	2.1	35	620	23.2	5.6 × 1041	2.5 × 10−17
Q0107−025	0.9545	=	1.9	28	410	13.9	2.0 × 1041	1.4 × 10−17
Q0107−0235	0.9574	=	3.1	90	2130	144.8	3.7 × 1042	1.3 × 10−16
PKS 2242−498	1.0011	=	3.4	71	2260	54.6	2.3 × 1042	1.2 × 10−16
PKS 0355−483	1.0128	=	2.8	50	1310	70.5	1.0 × 1042	6.2 × 10−17
HE 0112−4145	1.0238	=	1.8	38	890	70.8	2.7 × 1042	7.8 × 10−17
HE 0439−5254	1.0530	=	2.5	47	970	34.3	9.1 × 1041	5.6 × 10−17
HE 2305−5315	1.0733	<	0.6	15	100	3.9	9.6 × 1040	1.1 × 10−17
HE 1003+0149	1.0807	=	2.8	53	1660	94.8	7.0 × 1042	3.9 × 10−16
HE 0331−4112	1.1153	=	1.8	32	520	25.9	6.9 × 1041	3.1 × 10−17
TXS 0206−048 l	1.1317	=	5.0	200	10800	257.1	9.7 × 1042	3.0 × 10−16
Q1354+048	1.2335	=	6.4	126	2860	66.5	2.9 × 1042	9.9 × 10−17
J0154−0712	1.2957	=	2.8	63	1340	41.9	2.2 × 1042	1.2 × 10−16
Q1435−0134	1.3117	=	2.7	63	2010	76.0	7.3 × 1042	2.6 × 10−16
PG 1522+101	1.3302	=	2.9	50	1220	28.3	2.5 × 1042	1.6 × 10−16
HE 2336−5540	1.3531	<	0.4	7	30	3.6	5.0 × 1040	4.3 × 10−17
PKS 0232−04	1.4450	=	4.0	116	3670	97.1	7.8 × 1042	3.6 × 10−16

Notes.

^a Quasar systemic redshift repeated from Table 1 for convenience. ^b Nebula detection flag, with "=" indicating a detection and "<" indicating a nondetection, in which case all nebula measurements represent upper limits. ^c FWHM of the [O ii] surface-brightness profiles, which have systematic uncertainties of ≈10%. ^d Projected linear sizes of the [O ii] nebula defined as the major axis diameter, which have systematic uncertainties of 15%. ^e Projected areas of the [O ii] nebula, which have systematic uncertainties of ≈20%. ^f Total signal-to-noise ratios of the [O ii] nebula. ^g Total luminosities of the [O ii] nebula, which have systematic uncertainties of ≈25%. ^h Mean cosmological-dimming-corrected surface brightness measured in an annulus with inner and outer radii of 10 and 30 kpc, which have systematic uncertainties of ≈25%. ⁱ For more details on this nebula and its environment, see Johnson et al. (2018). ^j For more details on this nebula and its environment, see Zhao & Wang (2023) and Liu et al. 2024. ^k For more details on these nebulae and their environments, see M. C. Chen et al. (2023). ^l For more details on this nebula and its environment, see Johnson et al. (2022).

Download table as:  ASCIITypeset image

While half of the 30 quasars in this work have HST images available in the archive, the images were acquired with long exposures to enable study of the morphologies of faint foreground galaxies and to cross-correlate with CGM/IGM absorption. As a result, the bright quasars are saturated, preventing studies of stellar continuum from the host galaxies. However, the mean effective radius of the stellar component of massive late- and early-type galaxies at these redshifts is R_eff ≈ 9–10 kpc (e.g., van der Wel et al. 2014), corresponding to a diameter of ≈20 kpc. If quasar hosts represent 2σ outliers in the mass–size relation, we expect their effective radii to be ≈15 kpc, corresponding to a diameter of 30 kpc.

To gain insights into the extent of the nebulae in comparison to the expected extents of the stellar components of their host galaxies, we compute the fraction of systems with nebulae larger than 100, 50, 30, and 20 kpc in the following. Among the 30 quasars, seven (${23}_{-6}^{+9} \%$) have detected [O ii] nebulae with major axis sizes greater than 100 kpc, 20 (${67}_{-9}^{+7} \%$) greater than 50 kpc, 25 (${83}_{-9}^{+5} \%$) greater than 30 kpc, and 27 (${90}_{-8}^{+3} \%$) greater than 20 kpc. The median size of the [O ii] nebulae detected around the CUBS and MUSEQuBES quasars is 60 kpc. The sizes of the [O ii] nebulae are therefore larger than the expected sizes of the stellar component of their host galaxies, indicating that the majority of UV-luminous quasars at z = 0.4–1.4 are surrounded by [O ii]-emitting, CGM-scale nebulae.

Interestingly, the smaller 20–60 kpc nebulae are generally well centered on the quasars and relatively symmetric, but the larger ones often exhibit irregular morphology and may not be well centered (e.g., PKS 0405−123 and Q1354+048). In some cases, the morphologies of the nebulae are coincident with nearby galaxies with redshifts similar to the quasars, as previously reported for PKS 0405−123 (Johnson et al. 2018), TXS 0206−048 (Johnson et al. 2022), and HE 0238−1904 (Liu et al. 2024). This can be explained if the larger nebulae often result from ram pressure and tidal stripping experienced during galaxy interactions, as previously suggested by Stockton & MacKenty (1987). The nebulae do not exhibit obvious coincidences with the radio lobe and jet orientation of the radio-loud quasars, disfavoring jet-driven outflow origin (though see Fu & Stockton 2009). Further investigation of the coincidences between the extended nebulae and interacting galaxies in the quasar-host environment will require detailed investigation of the group environments like those in Johnson et al. (2018, 2022), Helton et al. (2021), and Liu et al. (2024).

Empirical characterizations of cool CGM around quasars have advanced significantly in the last decade, thanks to a combination of large ground-based surveys like SDSS that enable studies of Mg ii absorption at z ≲ 2 (e.g., Bowen et al. 2006; Farina et al. 2014; Johnson et al. 2015; Z.-F. Chen et al. 2023), dedicated absorption surveys with large telescopes at z = 2–4 in rest-UV absorption (e.g., Hennawi et al. 2006; Prochaska et al. 2013; Lau et al. 2016), and observations of Lyα emission nebulae at z > 2 (e.g., Borisova et al. 2016; Cai et al. 2019; O'Sullivan et al. 2020; Fossati et al. 2021; Mackenzie et al. 2021). The discovery of rest-frame optical emission from >20 to 30 kpc ionized nebulae around the majority of UV-luminous quasars at z = 0.4–1.4 represents a new opportunity to study the relationship between cool CGM and the galactic systems they surround, while also opening a new set of physical diagnostics based on nonresonant emission lines.

5.1. Possible Relationships between Nebulae and Quasar Properties

Previous studies of the cool, absorbing CGM around quasars at z < 2 demonstrated that strong Mg ii absorbers, which typically trace cool (T ≈ 10⁴ K), high-column-density ($\mathrm{log}N({\rm{H}}\,{\rm\small{I}})/{\mathrm{cm}}^{-2}\gtrsim 17$) enriched gas, are common around quasar-host galaxies (Bowen et al. 2006; Farina et al. 2014; Johnson et al. 2015). The Mg ii covering factors observed around quasar hosts is significantly greater than those found around luminous red galaxies (LRGs; e.g., Huang et al. 2016) and comparable to or higher than those of massive star-forming galaxies (e.g., Huang et al. 2021). Moreover, Johnson et al. (2015) identified a strong correlation between Mg ii absorption covering fractions at d < 200 kpc and quasar luminosity, with luminous quasars ($\mathrm{log}{L}_{\mathrm{bol}}/\mathrm{erg}\,{{\rm{s}}}^{-1}\gt 45.5$) exhibiting Mg ii covering fractions approximately double those of lower-luminosity quasars. Early studies of large-scale optical emission around quasars based on a combination of narrowband [O iii] imaging and follow-up long-slit spectroscopy suggested that emitting nebulae on scales of tens of kiloparsecs are common around radio-loud quasars and uncommon around radio-quiet ones (e.g., Boroson et al. 1985; Stockton & MacKenty 1987). This correlation between extended nebulosity and radio properties led to suggestions that the nebulae arise from radio-mode feedback (e.g., Fu & Stockton 2009).

To search for correlations between the properties of the quasars and surrounding nebulae, Figure 3 displays plots of the sizes and dimming-corrected surface brightnesses (SB₁₀₋₃₀(1 + z)⁴) of the nebulae versus quasar redshift, bolometric luminosity, black hole mass, radio loudness, and [O ii] rest-frame equivalent width. The panels in Figure 3 exhibit substantial scatter with no clear trend between quasar redshifts, bolometric luminosities, and the sizes and areas of the [O ii] nebulae. One the other hand, the [O ii] emission levels in the quasar spectra (with negative values indicating stronger emission relative to the continuum) appear to be weakly correlated with the sizes and surface brightnesses of the nebulae, though with substantial scatter (see the right panel of Figure 3). We quantified the correlations and level of significance with the generalized Kendall τ test, accounting for upper limits and nondetections following Isobe et al. (1986). Tests for correlations between the properties of the nebulae and quasar redshifts, bolometric luminosities, black hole masses, and radio loudness returned no statistically significant results (p-values between 0.2 and 0.8), consistent with the lack of a visual trend seen in Figure 3. On the other hand, the apparent correlation between [O ii] equivalent width and the sizes and surface brightnesses of the nebulae is marginally significant, with a correlation coefficient of −0.3 (indicating a positive correlation between the emission-line strengths in the quasar spectra and large-scale nebulae) and a p-value of 0.01 (2.5σ). To demonstrate that the extended nebulae themselves are not responsible for this correlation, we remeasured the [O ii] equivalent widths in quasar spectra extracted with an aperture radius of 06 (3–5 kpc at z = 0.4–1.4) and found consistent results. The tentative correlation between [O ii] equivalent width and the nebulae suggests a general correlation between circumnuclear and ISM-scale gas observed in the quasar spectra and larger-scale gas supplies in the massive halos that host UV-luminous quasars (see Li et al. 2024). This trend is reminiscent of a previously reported potential correlation between ISM or nuclear [O ii] emission in LRGs and Mg ii absorption in their CGM at projected distances of d ≲ 50 kpc (Huang et al. 2016).

Zoom In Zoom Out Reset image size

The lack of correlations between the properties of the nebulae and the bolometric luminosities of the quasars that they surround suggests that the presence of the emitting gas may be unrelated to quasar outflows in many cases. Alternatively, the lack of correlations between the nebulae and quasar luminosity could be partly due to a lack of dynamic range in the sample. The CUBS and MUSEQuBES surveys do not include any quasars of $\mathrm{log}{L}_{\mathrm{bol}}/\mathrm{erg}\,{{\rm{s}}}^{-1}\lt 45.5$. The correlation between Mg ii absorption at d < 200 kpc and quasar luminosity reported by Johnson et al. (2015) sets in when comparing luminous quasars of $\mathrm{log}{L}_{\mathrm{bol}}/\mathrm{erg}\,{{\rm{s}}}^{-1}\gt 45.5$ with less luminous ones of $\mathrm{log}{L}_{\mathrm{bol}}/\mathrm{erg}\,{{\rm{s}}}^{-1}\lt 45.5$, and this study reported no trend within the more luminous quasar population. Furthermore, previous studies of the extended narrow-line regions around obscured AGNs found a trend between luminosity and [O iii]-detected narrow-line region size, with AGNs of luminosity $\mathrm{log}{L}_{\mathrm{bol}}/\mathrm{erg}\,{{\rm{s}}}^{-1}\lt 45.5$ exhibiting narrow-line region sizes of R_NLR ≲ 6 kpc and more luminous obscured quasars exhibiting larger extended narrow-line regions with R_NLR ≈ 10–15 kpc (Sun et al. 2017). We note that these narrow-line region scales are smaller than most of the nebulae found around the CUBS and MUSEQuBES quasars, which is likely the result of the choice of emission line ([O ii] versus [O iii]) or increased sensitivity of MUSE relative to more traditional long-slit spectrographs. Consequently, determining whether the correlation between cool CGM and quasar luminosity extends to ionized, line-emitting nebulae on circumgalactic scales will require samples of lower-luminosity systems observed with wide-field IFSs.

The lack of correlation between the presence, size, and surface brightness of the [O ii]-emitting nebulae and radio loudness also stands in contrast to previous studies, which found [O iii]-emitting nebulae around radio-loud quasars but not radio-quiet ones (e.g., Boroson et al. 1985; Stockton & MacKenty 1987). While the nebulae presented in this work are observed in [O ii] for the sake of uniformity given our lack of [O iii] coverage for the z > 0.85 systems, the results are not substantially different for [O iii] within the lower-z CUBS and MUSEQuBES quasars.

Our finding of common large nebulae around radio-quiet quasars compared to their absence in previous surveys may be the result of the sensitivity of MUSE. The larger telescope aperture, higher throughput, and higher spectral resolution of MUSE result in dramatically improved surface-brightness limits compared to previous long-slit spectroscopy and narrowband imaging surveys. Indeed, Lyα nebulae around z > 2 quasars were commonly detected around radio-loud systems but were thought to be rare around radio-quiet ones until new, more sensitive observations with MUSE and Keck Cosmic Web Imager became available. The increased sensitivity resulted in the discovery of ubiquitous, 100 kpc-scale nebulae around radio-quiet quasars at high redshift, despite nondetections with less sensitive instruments (see discussion in Borisova et al. 2016 and Cantalupo 2017). However, there is no correlation between surface brightness and radio loudness within the sample of 30 quasars studied here. Reconciling this result with previous studies that found a correlation between radio loudness and the presence of [O iii]-emitting nebulae will likely require larger samples observed with wide-field IFS data covering both [O ii] and [O iii].

5.2. Comparison with Giant Lyα Nebulae Observed around Quasars at z = 2–4

The discovery that half of luminous quasars at z ≈ 1 are surrounded (in projection) by [O ii]-emitting nebulae with sizes greater than 60 kpc indicates that dense, cool CGM is common in quasar-host halos at this epoch. However, this represents a significantly lower incidence and smaller size scale than findings of ubiquitous, ≳100 kpc H i Lyα-emitting nebulae around quasars at z > 2 (Borisova et al. 2016; Cai et al. 2019; O'Sullivan et al. 2020; Fossati et al. 2021; Mackenzie et al. 2021). However, physical interpretation of this apparent difference requires accounting for the expected [O ii]-to-H i Lyα line ratio and the difference in surface brightness dimming between z ≈ 1 and z > 2.

The radiative transfer involved in modeling Lyα is complex, so we opt for an empirical approach and adopt a measured Lyα/Hα ratio observed around quasars and then combine this with model expectations for Hα/[O ii] to determine whether we expect the types of nebulae observed in Lyα at z = 2–4 around quasars to be observable in [O ii] at z ≈ 1. To set expectations for Hα/[O ii], we considered emission models of moderately enriched, dust-free gas photoionized by an AGN using Cloudy v17.03 (Ferland et al. 2017), with model parameters motivated by Groves et al. (2004). In particular, we chose a power-law ionizing spectrum with spectral slope ranging from α = −2.0 to −1.4. We note that the range of predicted Hα/[O ii] ratios for AGN-photoionized gas are similar to predictions for fast shocks (e.g., Allen et al. 2008) and star-forming galaxies (e.g., Kewley et al. 2004), indicating that our conclusions are not sensitive to the ionization source.

One of the dominant drivers of variation in expected Hα/[O ii] ratios is gas metallicity. While the metallicity of the CGM around quasars at z ≈ 1 is not well studied due to observational limitations, Liu et al. (2014) inferred approximately solar metallicity in the nebula around HE 0238−1904 for regions detectable in weak-line diagnostics. Furthermore, the cool absorbing CGM around luminous quasars at z = 2–4 from Lau et al. (2016) and LRGs at z = 0.4 from Zahedy et al. (2019) exhibit typical metallicities of Z ≈ 0.25 Z_⊙, though with substantial scatter in the case of LRGs. We therefore adopt a fiducial metallicity of Z = 0.25 Z_⊙ for the AGN-photoionized gas models but also explore a range of Z = 0.1–1.0 Z_⊙. The presence of substantial dust in large nebulae around quasars is disfavored in cases where Balmer line ratios are measurable (e.g., Helton et al. 2021; Liu et al. 2024). The left axis of Figure 4 shows the Hα/[O ii] ratio versus [O iii]/[O ii] over the range of [O iii]/[O ii] observed in extended nebulae around z ≈ 0.6 quasars (Johnson et al. 2018; Helton et al. 2021; Liu et al. 2024).

Zoom In Zoom Out Reset image size

To translate this expected Hα/[O ii] into expected Lyα/[O ii], we multiply by a factor of 3.7 corresponding to the observed Lyα/Hα ratio of 3.7 ± 0.3 observed in nebulae around three z ≈ 3 quasars from Langen et al. (2023). We note this ratio is a factor of ≈2–3 below that expected from pure recombination radiation. Langen et al. (2023) and previous works (e.g., Cantalupo et al. 2019) attribute this decreased Lyα/Hα to radiative transfer effects combined with aperture losses. Combining this empirical Lyα/[O ii] ratio with the modeled Hα/[O ii], we expect log (Lyα/[O ii]) ≈ 0.25–1.4 assuming AGN photoionization, Lyα emission arising primarily from recombination with moderate scattering, and modest metal enrichment levels. The expected Lyα/[O ii] ratios are shown as a function of [O iii]/[O ii] in the right axis of Figure 4. The presence of dust would only increase the strength of [O ii] relative to Lyα.

Surveys of Lyα emission around quasars at 〈z〉 ≈ 3.3 such as Borisova et al. (2016) and Mackenzie et al. (2021) achieve observed-frame surface-brightness limits similar to the ones enabled by the CUBS and MUSEQuBES observations. However, the difference in surface-brightness dimming at the mean redshift of the Borisova et al. (2016) and Mackenzie et al. (2021) studies versus 〈z〉 = 0.95 in this work introduces a factor of (1 + 3.3)⁴/(1 + 0.95)⁴ ≈ 24 difference in rest-frame surface-brightness sensitivity. This factor of 24 difference in surface-brightness dimming is equivalent to 1.4 dex, compensating for the expectation that [O ii] will be a weaker line than Lyα. Therefore, if the cool, Lyα-emitting CGM observed around quasars primarily arise from recombination radiation with only modest scattering, if the quantities of cool, emitting CGM around quasars do not evolve significantly between z ≈ 3 and z ≈ 1, if CGM metallicities around z ≈ 1 quasars are not substantially lower than expected from absorption surveys, and if quasar episodic lifetimes are not smaller at intermediate redshift compared to z ≈ 3, we expect that the [O ii] nebulae observed in this work would be similar in size and incidence rate to Lyα at higher redshift.

In contrast to this expectation, the [O ii]-emitting nebulae found around CUBS and MUSEQuBES quasars are smaller and less common than their Lyα counterparts at z ≈ 3, as shown in the left panel Figure 5. In particular, while 100 kpc-scale Lyα nebulae are fairly ubiquitous at z = 3, only seven of the 30 (${23}_{-6}^{+9} \%$) of the CUBS and MUSEQuBES quasars exhibit such large [O ii] nebulae. The smaller sizes of the [O ii] nebulae at z ≈ 1 are even starker when considered relative to the estimated virial radii of the quasar-host halos, which are ≈120 kpc at z ≈ 3 (Shen et al. 2007) compared to ≈500 kpc for UV-luminous quasars at z ≈ 1 (Li et al. 2024). The fact that [O ii]-emitting nebulae are smaller and less common around z ≈ 1 quasars than Lyα nebulae around z = 3 ones can be explained if a substantial fraction of Lyα in the outer regions of the z = 3 nebulae arise from scattered Lyα photons in a more neutral phase of the CGM that does not emit in [O ii]. The contrast could also be explained if there is substantial evolution in the quantities or density distribution of the cool, dense phase of halo gas around quasars between z ≈ 3 and z ≈ 1 or if the metallicities of the CGM of z ≈ 1 quasar hosts on ≈50 kpc scales are substantially lower than 0.25 Z_⊙. We note, however, that the Z ≈ 0.25 Z_⊙ metallicity used to estimate the [O ii]-to-Lyα ratio in Figure 4 is based on absorption measurements in massive halos at distances of d ≈ 50–100 kpc (e.g., Lau et al. 2016; Zahedy et al. 2019), similar to the size scale of the Lyα nebulae. Further, Johnson et al. (2015) found little evidence for evolution in the Mg ii covering factor observed in the CGM of quasars over this same cosmic period. Alternatively, the difference in the size and incidence of Lyα nebulae at z ≈ 3 and [O ii]-emitting nebulae at z ≈ 1 could also be explained if luminous quasar episodic lifetimes at z ≈ 1 are shorter than at z ≈ 2–3. Further investigation of the differences in optically emitting [O ii] nebulae versus UV-emitting Lyα nebulae around quasars requires multiwavelength observations at both redshifts.

Zoom In Zoom Out Reset image size

In this paper, we presented the first comprehensive search for nonresonant, optically emitting [O ii] nebulae around 30 UV-luminous quasars at z = 0.4–1.4 using deep IFSs acquired as part of the CUBS and MUSEQuBES. Our findings indicate that >60 kpc-scale circumgalactic nebulae are common around UV-luminous quasars, with approximately half of the sample exhibiting such large-scale emission, while ≈80%–90% exhibit emission on >20–30 kpc scales. Within the CUBS and MUSEQuBES quasar samples, the presence and size of the [O ii] nebulae are not correlated with quasar redshift, luminosity, supermassive black hole mass, or radio loudness. On the other hand, the sizes and surface brightnesses of the [O ii] nebulae are tentatively correlated (2.5σ significance) with the rest-frame [O ii] equivalent width detected in the quasar spectra due to circumnuclear and ISM-scale ionized gas. Further investigation of the relationship between quasars and their surrounding cool gas supplies will require larger samples extending over a broader range of quasar properties, particularly to lower luminosities and black hole masses.

The [O ii]-emitting nebulae around quasars at z ≈ 1 found in this work are less common and smaller on average than Lyα-emitting nebulae found around quasars at z > 2. This difference can be explained, for example, if the dense, cool phase of the CGM around quasars evolves significantly with redshift, if the metallicities in quasar-host CGM on ≈50 scales are lower than expected from absorption surveys, if quasar episodic lifetimes are shorter at z ≈ 1 compared to z ≈ 3, or if the scattering of Lyα photons off of more neutral gas contributes nonnegligibly to the outer regions of Lyα halos around quasars at z = 2–3. These scenarios can be tested by observing the nebulae with broader spectral coverage using upcoming near-IR IFSs such as MIRMOS (Konidaris et al. 2020) and HARMONI (Thatte et al. 2021), or with new UV space-based instrumentation targeting low-z systems in H i Lyα.

We thank the anonymous referee for constructive comments and suggestions that significantly strengthened this paper. We are grateful to Eric Bell and Kayhan Gültekin for insightful discussions which helped shape interpretation of the results. S.D.J. acknowledges partial support from grant Nos. HST-GO-15280.009-A, HST-GO-15298.007-A, HST-GO-15655.018-A, and HST-GO-15935.021-A. J.I.L. is supported by the Eric and Wendy Schmidt AI in Science Postdoctoral Fellowship, a Schmidt Futures program. Based on observations from the European Organization for Astronomical Research in the Southern Hemisphere under ESO. The paper made use of the NASA/IPAC Extragalactic Database, the NASA Astrophysics Data System, Astropy (Astropy Collaboration et al. 2018), and Aplpy (Robitaille & Bressert 2012).

The raw MUSE data, associated calibrations, and ESO reduced data cubes are publicly available in the ESO archive.