The first comprehensive study of a giant nebula around a radio-quiet quasar in the z < 1 Universe

ABSTRACT

We present the first comprehensive study of a giant, ≈70 kpc-scale nebula around a radio-quiet quasar at z < 1. The analysis is based on deep integral field spectroscopy with Multi-Unit Spectroscopic Explorer of the field of HE 0238−1904, a luminous quasar at z = 0.6282. The nebula emits strongly in [O ii], |$\rm H \beta$|⁠, and [O iii], and the quasar resides in an unusually overdense environment for a radio-quiet system. The environment likely consists of two groups which may be merging, and in total have an estimated dynamical mass of M_dyn ≈ 4 × 10¹³ to 10¹⁴ M_⊙. The nebula exhibits largely quiescent kinematics and irregular morphology. The nebula may arise primarily through interaction-related stripping of circumgalactic and interstellar medium (CGM/ISM) of group members, with some potential contributions from quasar outflows. The simultaneous presence of the giant nebula and a radio-quiet quasar in a rich environment suggests a correlation between such circum-quasar nebulae and environmental effects. This possibility can be tested with larger samples. The upper limits on the electron number density implied by the [O ii] doublet ratio range from |$\log (n_{\rm e, [O\,{\small II}]} /\mathrm{cm}^{-3})<1.2$| to 2.8. However, assuming a constant quasar luminosity and negligible projection effects, the densities implied from the measured line ratios between different ions (e.g. [O ii], [O iii], and [Ne v]) and photoionization simulations are often 10−400 times larger. This large discrepancy can be explained by quasar variability on a time-scale of ≈10⁴−10⁵ yr.

1 INTRODUCTION

Galaxy evolution is a complex process that involves gas inflows and outflows thought to control star formation and black hole growth (for a review, see Naab & Ostriker 2017). Observations of interstellar medium (ISM) gas masses and star formation rates suggest that massive star-forming galaxies have an ISM depletion time-scale much smaller than the age of the Universe at z < 3 (Kennicutt & Evans 2012; Tacconi et al. 2013). This can be explained if galaxies accrete gas from external sources to maintain their star-forming activity and black hole growth (though see Leitner & Kravtsov 2011). At the same time, the ISM of galaxies can lose gas through various processes including stellar (for a review, see Zhang 2018) and AGN feedback (for a review, see Fabian 2012), ram pressure stripping (e.g. Hester 2006), and tidal interactions with neighbouring galaxies (e.g. Marasco et al. 2016). Therefore, observations of the physical conditions, kinematics, and distribution of gas around galaxies can provide insights into the mechanisms governing galaxy formation and evolution. For these reasons, observations of the gaseous cosmic ecosystems of galaxies were highlighted as a key long-term priority by the 2020 Decadal Survey for Astronomy and Astrophysics (National Academies of Sciences 2021).

The properties of gas flows around galaxies, including their morphology and kinematics, can be directly traced by observations of giant gas nebulae with state-of-the-art wide-field integral field spectrographs (IFSs) such as the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) and the Keck Cosmic Web Imager (Martin et al. 2010). At z > 2, systematic IFS surveys around radio-quiet quasars discovered ubiquitous giant H i Lyα nebulae (e.g. Cantalupo et al. 2014; Borisova et al. 2016b; Cai et al. 2019; O’Sullivan et al. 2020; Fossati et al. 2021; Mackenzie et al. 2021). More recently, a study of the ionization states of one of these nebulae found that the gas has a surprisingly large density for halo-scale emission or a very broad density distribution (Cantalupo et al. 2019). However, due to redshifting of optical emission lines into the infrared, surface brightness dimming, and the faintness of galaxies at high redshift, more fully characterizing these z > 2 nebulae is time-consuming even with large space- or ground-based telescopes (though see Langen et al. 2023).

At low redshift, on the other hand, non-resonant emission lines such as [O ii], |$\rm H \beta$|⁠, and [O iii] are available at optical wavelengths, and collecting galaxy spectra is less expensive. The power of IFSs enabled the discoveries of giant nebulae around starburst galaxies, galaxy groups, and quasars (e.g. Epinat et al. 2018; Boselli et al. 2019; Chen et al. 2019; Rupke et al. 2019; Burchett et al. 2021; Zabl et al. 2021; Leclercq et al. 2022; Dutta et al. 2023a), arising from outflows, interactions, and filamentary accretion. These low-redshift nebulae provide an opportunity to study the physical conditions and the processes that may produce giant nebulae at higher redshift. Most published studies of giant nebulae around z < 1 quasars have focused on radio-loud systems (Johnson et al. 2018; Helton et al. 2021; Johnson et al. 2022), which represent a small fraction of the general quasar population (e.g. Kellermann et al. 1989). Furthermore, clustering measurements indicate that radio-loud quasars typically reside in massive galaxy groups with halo masses of |$M \sim 10^{13}\ {\rm {\rm M}_{\odot }}$| while the halo masses of more common radio-quiet systems are approximately five times lower on average (e.g. Shen et al. 2009). This mass miss-match and the possibility of radio jet feedback make the comparison between low-redshift giant nebulae around radio-loud quasars and high-redshift radio-quiet ones difficult.

Recently, Chen et al. (2023) demonstrated the existence of giant nebulae around two radio-quiet quasars as part of a study focused on turbulence using the observed velocity structure function. In this paper, we present the first comprehensive characterization of a giant nebula and associated galaxy environment around a radio-quiet quasar at z < 1, HE 0238−1904. Recently, this nebula was independently discovered and reported by Zhao & Wang (2023). However, our interpretation of the system differs substantially from the one presented by Zhao & Wang (2023) due to adoption of a significantly different quasar systemic redshift. In particular, Zhao & Wang (2023) adopted a Mg ii emission-based redshift of z = 0.631 from the Hamburg/ESO Survey of bright Quasars (Wisotzki et al. 2000). On the other hand, we adopt a redshift estimate of z = 0.6282 based on the [O ii] emission-line centroid measured in the spectrum of the quasar extracted from the same MUSE data set used to measure the kinematics of the giant nebula. The paper is organized as follows: In Section 2, we discuss the observations, data reduction, and processing. In Section 3, we describe our measurements and investigate the group environment and giant nebula properties. In Section 4, we investigate the origin of the nebula and the physical conditions of the gas. In Section 5, we summarize our findings and discuss their implications.

Throughout the paper, we adopt a flat Λ cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and |$H_{0} = 70 \, \rm km \, s^{-1} Mpc^{-1}$|⁠. All magnitudes are given in the AB system unless otherwise stated.

2 OBSERVATIONS AND DATA

The z ≈ 0.63 quasar HE 0238−1904 has high-quality archival UV Hubble Space Telescope (HST) absorption spectra used to study the circumgalactic medium (CGM) of the Milky Way (Zheng et al. 2019; Bish et al. 2021) and distant galaxies (Lehner et al. 2018; Muzahid et al. 2018) in addition to a highly ionized, fast outflow from the quasar itself (Muzahid et al. 2012; Arav et al. 2013). To identify faint foreground galaxies in the quasar field, we observed it with MUSE as part of the Quasar-field Blind Emitter Survey (MUSE-QuBES; Muzahid et al. 2020; Dutta et al. 2023b) on the Very Large Telescope (VLT; PI: J. Schaye, PID: 094.A-0131(B), and 096.A-0222(A)). MUSE is an integral-field spectrograph on the UT4 VLT with a field of view (FoV) of 1 × 1 arcmin² and a spaxel size of 0.2 arcsec in wide-field mode (WFM). MUSE covers the spectral range between |$4750\, \rm{\mathring{\rm A}}$| to |$9350 \, \rm{\mathring{\rm A}}$| and a resolution of R ∼ 3000. The MUSE observations are centred near the quasar sightline, and we obtained eleven exposures collected between 2014 November 18 and 2016 February 2 with a total exposure time of 8.75 h with median seeing full width at half-maximum (FWHM) conditions of 0.7 arcsec. At the redshift of HE 0238−1904, the MUSE FoV corresponds to a projected size of ≈400 proper kpc (pkpc) on a side, and the spectral coverage includes emission lines such as [O ii], Hβ, and [O iii]. These emission lines enable sensitive studies of any ionized nebulae and galaxies in the quasar’s environment.

To ensure robustness of results, we analyzed the MUSE data reduced through three independent pipelines including CubEx (Cantalupo et al. 2019), the MUSE GTO team pipeline (Weilbacher et al. 2014), and the ESO reduction pipeline (Weilbacher et al. 2012) and found consistent results with all three. All three pipelines include bias subtraction, flat fielding, wavelength calibration, geometric calibration, sky subtraction, flux calibration, and stacking of exposures. For the ESO reductions, we obtained the final, stacked datacube from the ESO Science Archive and performed additional post-processed sky subtraction with the Zurich Atmosphere Purge package (Soto et al. 2016). For simplicity, we converted the air wavelengths delivered by the three pipelines to vacuum.

To enable more sensitive and higher angular resolution photometric measurements of galaxies in the quasar field, we also obtained an image from the Advanced Camera for Surveys (ACS) on the HST with the F814W filter (PI: L. Straka, PID: 14660) with a total exposure time of 2182 s split between four dithered exposures. We obtained the reduced, stacked image from the Barbara A. Mikulski Archive for Space Telescopes (MAST). In addition, to measure the UV luminosity of the quasar, we obtained the archival UV spectrum from the Cosmic Origins Spectrograph (COS; Green et al. 2012) from MAST. The spectrum consists of a total exposure time of 14 400 and 7496 s in the G130M and G160M gratings, respectively (PI: J. Green and S. Penton, PID: 11 541 and 12505). We reduced and coadded the COS spectrum following procedures outlined in Johnson, Chen & Mulchaey (2015) and Chen et al. (2020).

2.1 Quasar light subtraction

HE 0238−1904 has a Gaia (Gaia Collaboration et al. 2018) G-band magnitude of m_G = 15.2, and this brightness combined with the broad wings of the MUSE point spread function (PSF) causes contamination of nearby galaxy spectra with quasar light. This contamination includes both continuum and line emission due to the unresolved narrow-line region in the nucleus. To study faint extended emission, we removed the contamination by performing quasar light subtraction as described in Helton et al. (2021). In summary, our method of quasar light subtraction does not rely on PSF measurements. Instead, it uses spectral information and the fact that quasars and galaxies have different spectral energy distributions (see also Rupke, Gültekin & Veilleux 2017; Chen et al. 2023).

In ground-based observations, the Earth’s atmosphere scatters bluer photons more than redder ones so that the PSF is wider at bluer wavelengths. The differential scattering makes the spectral slope observed in a spaxel depend on the angular separation from the quasar with steeper (shallower) slopes further from (closer to) the quasar centroid. To account for this, we used a two-component non-negative matrix factorization (NMF; Blanton & Roweis 2007; Ren et al. 2018) of the quasar light, with one component having a shallow slope and a second having a steep slope. Adding additional a third or fourth NMF component(s) did not noticeably improve the results. In general, the spectrum for each spaxel near the quasar has some light from the quasar and potentially nearby galaxies as well. To subtract quasar light while avoiding subtraction of galaxy light, we fit each spaxel with a linear combination of the two quasar non-negative components and the first two Sloan Digital Sky Survey-Baryon Oscillation Spectroscopic Survey (SDSS-BOSS) galaxy eigenspectra (Bolton et al. 2012) and then subtracted the quasar component of the model. Unlike with some other systems (e.g. Johnson et al. 2018), the host of HE 0238−1904 does not exhibit bright, extended starlight, so the contribution inferred by the galaxy model was not significant.

3 MEASUREMENTS AND ENVIRONMENT

3.1 Quasar properties

HE 0238−1904 is a luminous, radio-quiet quasar (Véron-Cetty & Véron 2006; Arav et al. 2013). To ensure self-consistent measurements of the quasar properties, we estimated its redshift, luminosity, and black hole mass using the MUSE spectrum extracted via MPDAF (Bacon et al. 2016) with a r = 3 arcsec aperture. To measure the systemic redshift of the quasar, we fit the [O ii]|$\lambda \lambda 3727, 3729$| doublet with a Gaussian profile following Hewett & Wild (2010) and found z = 0.6282 ± 0.0002, where the uncertainty represents the scatter between the [O ii] centroid and stellar absorption lines of SDSS quasars at similar redshift. This redshift is from a previously reported Mg ii based estimate from Wisotzki et al. (2000). Even so, a more recent Mg ii based redshift of z = 0.628 from Monroe et al. (2016) confirms our [O ii]-based redshift estimate. In general, quasar redshifts measured from the [O ii] doublet are more accurate than those measured from broad-lines like Mg ii, as we argue in Section 4.1.

In addition, we estimated the bolometric luminosity and the black hole mass of HE 0238−1904 by fitting the extracted MUSE spectrum with the Python QSO fitting code (PyQSOFit; Guo et al. 2019). PyQSOFit fits a quasar’s spectrum with a combination of a power-law continuum, Fe  ii template, and sets of Gaussian line profiles for both the broad and narrow lines. We modelled the |$\rm H\beta$| and [O  iii] spectral region with the continuum components, three Gaussian profiles for the broad Hβ, and two for the narrow Hβ and [O iii]. From the fit, we computed a monochromatic luminosity at 5100Å of |$\lambda L_{5100} \approx 1.6 \times 10^{46} \rm \ erg \, s^{-1}$| and a bolometric luminosity of |$L_{\rm bol} \approx 1.7 \times 10^{47} \, \rm erg \, s^{-1}$| using the bolometric correction factor from Richards et al. (2006). Finally, we inferred a black hole mass of |$M_{\rm BH}\approx 10^{9.8}\ {\rm {\rm M}_{\odot }}$| using the single-epoch virial theorem-based approach from Vestergaard & Peterson (2006). Following Kormendy & Ho (2013), this black hole mass corresponds to a stellar mass of |$M_{*}\approx 10^{12.0}\ {\rm {\rm M}_{\odot }}$| for the host galaxy, but we caution this stellar mass may be significantly overestimated due to uncertainty in single-epoch virial theorem-based black hole masses and observed scatter in the black hole mass–stellar mass relation. For example, if the true black hole mass is 1σ below the mean single-epoch virial theorem estimate, and the stellar mass is 1σ below the estimate from the black hole mass–stellar mass relation, the inferred stellar mass would be |$M_{*}\approx 10^{11.4}\ {\rm {\rm M}_{\odot }}$|⁠. Furthermore, the single-epoch virial theorem-based relation used here is not calibrated for quasars as luminous as HE 0238−1904, which may drive disc wind, erroneously inflating the black hole mass estimate. The fitted quasar spectrum is shown in Fig. 1.

3.2 Galaxy measurements and properties

To study the environment of HE 0238−1904, we conducted a galaxy survey by first identifying all continuum sources in MUSE and the ACS+F814W image. We identified continuum sources by running Source Extractor (SE; Bertin & Arnouts 1996) on a median MUSE white light image and the HST image separately. To ensure completeness, we also added sources based on visual inspection. Typically, sources are missing from MUSE due to biased background estimation caused by bright objects in the field or due to blending. Based on the background sky standard deviation and source counts in the ACS+F814W image, the imaging catalogue is complete for objects brighter than m_F814W ≈ 26−27, depending on angular size.

For each identified object, we extracted a MUSE spectrum with MPDAF with a circular aperture of r = 0.7 arcsec, which is roughly the size of the MUSE seeing FWHM. The choice of this modest aperture may result in some wavelength-dependent aperture losses but helps increase the signal-to-noise ratio (S/N) for redshift estimation. We then fit each spectrum as a linear combination of SDSS galaxy eigenspectra as described in Helton et al. (2021) to measure the source redshift. In summary, we computed the best-fitting linear combination on a grid from z = 0 to 1 with a step size of Δz = 0.0001 and recorded the goodness-of-fit statistic (χ²) over the entire grid. We adopted the redshift with the minimum global χ² as our initial solution. We then visually inspected each best-fitting model to ensure robustness and assigned the redshift quality. For galaxies with both emission and absorption lines, we masked out strong emission lines and measured the redshift based on stellar absorption features when possible to avoid a potential bias in redshift from large-scale nebulae in the field (which may not be closely associated with the galaxies in question). Finally, we classified our confidence in the redshift measurements based on the number of the detected spectral features. All of the galaxies in the quasar environment have two or more spectral features except for G11 and G18. According to Helton et al. (2021), the uncertainty in galaxy redshifts measured in MUSE spectra with these techniques is |$\sigma \approx \rm 20 \, km \, s^{-1}$|⁠. Comparing the continuum source catalogue and the corresponding redshift measurements, the redshift survey is approximately 100 per cent complete for sources brighter than m_F814W ≈ 24 and approximately 95 per cent complete for those brighter than m_F814W ≈ 25. For comparison, an L_* galaxy at z ≈ 0.6 has m_F814W ≈ 20.6 assuming the luminosity function from Faber et al. (2007). The high completeness of the galaxy survey at faint magnitudes enables us to study the origins of nebulae, even if they arise from interactions involving relatively faint dwarf galaxies.

To examine properties of the quasar host environment, we identified candidate group members based on their line-of-sight (LOS) velocities relative to the quasar (Δv = v − v_QSO). In particular, we selected galaxies with |$|\Delta v| \lt \rm 2000\ km \, s^{-1}$|⁠. We inferred the physical properties of the selected galaxies with Bagpipes (Carnall et al. 2018, 2019). Bagpipes performs stellar population synthesis (SPS) with a stellar evolution model from Bruzual & Charlot (2003), an initial mass function from Kroupa (2001), and the Bayesian inference package Multinest (Feroz, Hobson & Bridges 2009; Buchner et al. 2014; Feroz et al. 2019). We fit both spectroscopic and photometric data simultaneously with Bagpipes. Many of the galaxies in our sample only have one photometric data point available, necessitating the use of the spectra to further inform the SPS. In our fitting procedure, we assumed an exponential star formation history with e-folding time-scale of |$0.01 \lt \rm \tau /Gyr \lt 8.00$|⁠, solar stellar metallicity, and dust attenuation model from Calzetti et al. (2000) with |$0 \lt A_V/\rm mag \lt 2$|⁠. The choice of exponentially declining star formation histories enables more direct comparison with surveys such as MUSE-Wide (Urrutia et al. 2019) and the MUSE Ultra DEEP Field (Fossati et al. 2019). We introduced a second-order multiplicative polynomial to reconcile the potential artificial differences between SED measured in photometry and spectra. This polynomial accounts for systematic uncertainty in the MUSE flux due to wavelength-dependent aperture losses and uncertainty in the flux calibration (Weilbacher et al. 2020). We also used Bagpipes spectrum noise scaling to allow the relative weighting of the photometry and spectrum to be a nuisance parameter. We note that the results are not sensitive to this scaling in our case (see Carnall et al. 2019). In addition to the ACS+F814W photometry, we also included grizY photometric data from the Dark Energy Survey (DES; Abbott et al. 2021) available for 16 galaxies. The resulting stellar mass estimates and dust attenuation A_V values are reported in Table 1. The stellar masses have associated systematic uncertainties of ≈0.2 dex. Galaxies close to the quasar (G1–G7) are contaminated by the quasar light, and we used the quasar-light subtracted spectra for Bagpipes fitting when possible. Galaxies G1, G3, G11, G13, G18, G20, and G31 do not have a stellar mass estimate because their continua are too faint or are too badly contaminated by the quasar continuum. To further characterize these galaxies, we also report 4000 Å  break strength (D4000; Gallazzi et al. 2005) and rest-frame B-band absolute magnitude with K-corrections calculated using templates from Coleman, Wu & Weedman (1980) chosen based on the strength of the 4000 Å break. The IDs, galaxy coordinates (RA, Dec.), redshifts, ACS+F814W apparent magnitudes, absolute B-band magnitudes, adopted K-correction templates (S0, Scd, or Irregular), and D4000 measurements are reported in Table 1, along with the angular distances, projected distances, and LOS velocity differences from the quasar sightline. The locations of these galaxies are shown in Fig. 2 and several example MUSE spectra are overplotted with their best-fitting PCA spectral models in Fig. 3. An interactive view of the galaxy environment and spectra is available online.¹

3.3 The galactic environment

In the MUSE field of HE 0238−1904, we identified 35 galaxies, including the quasar host, with LOS velocities |$|\Delta v| \lt \rm 2000 \, km \, s^{-1}$| of the quasar systemic velocity, which is sufficient to encompass most members of even massive galaxy clusters. Fig. 2 shows a 1.5 × 1.5 arcmin² FoV image from the ACS |$+$| F814W observations of the field where we marked the quasar with a grey star and labelled galaxies with circles as well as their ID. The colour of the circle represents the LOS velocity of each galaxy relative to the quasar. Additionally, we display the 1 × 1 arcmin² MUSE FoV, and a smaller 30 × 30 arcsec² region which is the focus of later figures in this work.

Among the 35 galaxies in the environment of HE 0238−1904, four (two) exhibit stellar masses of |$\log (M_*/{\rm {\rm M}_{\odot }})\gt 10.5\ (\gt \!11)$| (excluding the quasar), indicating a significant overdensity and likely a massive group. To further characterize the environment, we show the distribution of galaxies’ LOS velocities relative to the quasar (Δv = v − v_QSO) in the bottom right-hand panel of Fig. 2. The LOS velocity distribution peaks around |$\rm -100 \, km \, s^{-1}$| but exhibits a non-Gaussian tail toward higher velocity of |$\rm +100\, km \, s^{-1}$| to |$\rm +1400 \, km \, s^{-1}$|⁠. There is a clear trend between LOS velocity and location on the sky visible in Fig. 2 with galaxies with |$\Delta v \gt 0 \rm \, km \, s^{-1}$| largely falling north-east of the quasar and those with |$\Delta v \lt 0 \, \rm km \, s^{-1}$| falling near the quasar or south-west of it. To better visualize the location−velocity trend, we divided the field into two regions, one NE of the quasar and one SW of it. The NE (SW) one is marked by an orange (purple) trapezoid in Fig. 2. We also show the LOS velocity distribution of the galaxies in each trapezoidal region by the corresponding histograms in the inset panel in Fig. 2. The peak and the tail in the histogram correspond closely to these two regions, respectively. The non-Gaussian LOS velocity distribution and correlation with spatial location suggest that the overdensity near the quasar host may consist of two distinct, but possibly interacting, galaxy groups.

To quantify the velocity dispersions of these two potential groups, we fit two Gaussians to the entire LOS velocity distribution. This results in one narrow, blueshifted Gaussian and one broader, redshifted one. The blueshifted Gaussian has a mean LOS velocity of |$\Delta v_{\rm group} = \rm -99 \pm 25 \, km \, s^{-1}$| and a 1D velocity dispersion of |$\sigma _{\rm group} = \rm 92 \pm 50 \, km \, s^{-1}$| and includes |$\approx 35~{{\ \rm per\ cent}}$| of the galaxies near HE 0238−1904. The redshifted Gaussian has |$\Delta v_{\rm group} = \rm 629 \pm 140 \, km \, s^{-1}$| and |$\sigma _{\rm group} = \rm 506 \pm 90 \, km \, s^{-1}$| and includes |$\approx 65~{{\ \rm per\ cent}}$| of the galaxies. In both cases, the uncertainty estimates are based on bootstrap resampling. While the Gaussian fitting did not include any spatial information, the two Gaussians closely match the purple and orange velocity histograms formed from a spatial separation (see Fig. 2). These fitting results suggest that the environment around the quasar includes one massive group at |$\Delta v_{\rm group} \approx \rm 600 \, km \, s^{-1}$| and one less massive group closer to the quasar velocity. Assuming each group is virialized, we estimate dynamical masses of |$M_{\rm dyn} \sim 9.8 \times 10^{13}\ {\rm {\rm M}_{\odot }}$| and |$M_{\rm dyn}\sim 5.7 \times 10^{11}\ {\rm {\rm M}_{\odot }}$| (Munari et al. 2013) for the richer, redshifted group and less rich, blueshifted group, respectively. To place a lower limit on the mass estimate, we fit a single Gaussian to galaxies with |$\Delta v \gt 200 \, \rm km \, s^{-1}$|⁠. We found the velocity dispersion is |$\approx 400 \, \rm km \, s^{-1}$|⁠, corresponding to a mass of |$M_{\rm dyn} \sim 3.8 \times 10^{13}\ {\rm {\rm M}_{\odot }}$|⁠. The mass range of |$M_{\rm dyn} \approx 4 \times 10^{13}{-}10^{14}\ {\rm {\rm M}_{\odot }}$| is consistent with massive group or modest mass cluster. However, we caution that the assumption that the groups are virialized introduces additional uncertainty given the complex environment. Finally, in Fig. 2, we show the stellar mass weighted group centre as a white asterisk, and membership weighted (⁠|$\frac{P_{\rm blue/red}}{P_{\rm blue} + P_{\rm red}}$|⁠) centres as red and blue asterisks for the richer, redshifted group and less rich, blueshifted group, respectively.

To test the expectation that dynamically more massive groups will contain more massive galaxies, we investigate the most massive galaxies in each group. G8 and G22 are the most massive galaxies with a stellar mass of |$\log (M_*/{\rm {\rm M}_{\odot }})=10.4$| and 10.1, respectively, in the less rich, blueshifted group. On the other hand, the richer, redshifted group includes two massive elliptical galaxies, G10 and G34, with log (M_*/M_⊙) = 11.5 and 11.2, respectively. Furthermore, the richer, redshifted group contains a massive disc galaxy, G33, with |$\log (M_*/{\rm {\rm M}_{\odot }})=10.8$|⁠. This is consistent with HE 0238−1904 residing in an overdense region likely made of two groups with the redshifted one being richer and more massive. However, the quasar redshift falls between the centroids of the two groups indicating that it could arise in either or truly be located between them. Despite the large uncertainty in the stellar mass of the quasar host galaxy (see Section 3.1), the large black hole mass suggests it is a massive galaxy, possibly the largest in the overdensity around HE 0238−1904. It is therefore more probable that HE 0238−1904 resides in the richer, redshifted group. None the less, we cannot completely rule out the possibility that HE 0238−1904 originates from the less rich, blueshifted group. In either case, the dynamically rich and likely unrelaxed environment could result in galaxy interactions that can produce giant nebulae via ram pressure and tidal stripping.

3.4 Nebular environment

Due to ionizing radiation from the accretion disc, wide-field IFS observations of quasar fields often find large nebulae (Johnson et al., in preparation). To search for the nebulae around HE 0238−1904, we conducted continuum subtraction of the datacube locally for the |$\rm [O\,{\small II}]$|⁠, |$\rm H\beta$|⁠, and |$\rm [\rm O\,{\small III}]$| emission lines around the quasar. For continuum fitting near each of the three lines, we masked the spectral region within |$\pm 500{-}1000 \, \rm km\, s^{-1}$| of the expected observed wavelength at the quasar’s redshift. We fine-tuned the masked region individually for each of the three lines to avoid skyline contamination and to account for the larger width |$\rm [O\,{\small II}]$| doublet. For each spaxel in the masked datacube, we then fit a third-order polynomial to the continuum regions around each line and subtracted the best-fitting model to complete the continuum subtraction.

This continuum-subtracted MUSE datacube enabled the discovery of a giant ionized nebula in |$\rm [O\,{\small II}]$|⁠, |$\rm H\beta$|⁠, and |$\rm [\rm O\,{\small III}]$| around HE 0238−1904 with a total area of ≈5000 kpc² which is visualized in Fig. 4. This nebula surrounds the quasar with projected radii of |$d \rm \approx \! 30 \ to \ 50 \, pkpc$| and with LOS velocities of |$\Delta v \approx -250 \ \rm to \ {+}250 \ km \, s^{-1}$| from the quasar. The nebula is more extended to the south-east and the south-west of the quasar. The south-east extension of the nebula is spatially coincident with galaxies G1, G3, G4, and G5. Additionally, the tail extending south-west of the quasar is distinct from but approximately in the direction of G8.

To examine the nebula and any relationship with galaxies in the quasar environment, we show |$\rm [O\,{\small II}]$| and |$\rm [\rm O\,{\small III}]$| emission contours over the HST image in panel (a) of Fig. 4. We also display a nebular LOS velocity map in panel (b) and a |$\rm [\rm O\,{\small III}]/[O\,{\small II}]$| line ratio map in panel (c). We constructed these two maps by jointly fitting Gaussian line profiles to the continuum-subtracted |$\rm [O\,{\small II}]$|⁠, |$\rm H\beta$|⁠, and |$\rm [\rm O\,{\small III}]$| datacubes. Instead of fitting the spectrum of each individual spaxel, we averaged over circular apertures of r = 1 arcsec to enhance S/N. We chose this aperture radius based on experimentation to visualize even faint parts of the nebula. These two maps provide an opportunity to study the spatial dependence of the kinematics and the ionization state of the gas. In addition, we show three panels of narrow-band images generated from the continuum-subtracted datacubes for each of |$\rm [O\,{\small II}]$| and |$\rm [\rm O\,{\small III}]$| in velocity ranges of |$-300 \rm \ to \ {-}100 \ km \, s^{-1}$|⁠, |$-100 \rm \ to \ {+}100 \ km \, s^{-1}$|⁠, and |$+100 \rm \ to\ {+}300 \ km \, s^{-1}$| in panel (d)–(f) and (g)–(i), respectively.

The nebula exhibits an irregular morphology but with a spatial trend in kinematics. In particular, the region north of the quasar is redshifted relative to the quasar and has an LOS velocity of |$\Delta v = 0{-}250\rm \, km \, s^{-1}$|⁠. The region south of the quasar including the tail to the West is mainly blueshifted relative to the quasar but with a small redshifted region in the most southern points. This southern region is spatially coincident and potentially kinematically coincident with G1, G3, G4, and G5. However, the continua of these galaxies are too faint to measure stellar absorption-based redshifts. This raises the possibility that their nebular spectra may be contaminated by the surrounding nebulae, resulting in a biased redshift measurement. In the case of G3 and G4, the linewidth of the nebular emission near the galaxies is significantly narrower than the more extended emission from nearby parts of the giant nebula, indicating that the galaxy line emission likely arises in the ISM of the two dwarfs.

The nebula also shows a spatial trend in the ionization-state-sensitive |$\rm [\rm O\,{\small III}]/[O\,{\small II}]$| line ratio. The majority of the nebula is |$\rm [O\,{\small II}]$| dominated but the region south-east of the quasar has greater |$\rm [\rm O\,{\small III}]$| emission, particularly, at a few |$\rm [\rm O\,{\small III}]$| knots near G1, G3, and G5. The knots near G3 and G5 have the highest surface brightness in the nebula. Furthermore, the bright region extending to the south of the brightest knot near G3 is reminiscent of a tidal tail.

To better explore the properties of the nebula, we selected several representative regions in it and extracted their full spectra to infer physical conditions from both strong ([O ii], |$\rm H\beta$|⁠, and [O iii]) and weak lines ([Ne v]λ3427, |$\rm H\delta$|⁠, |$\rm H\gamma$|⁠, [O iii]λ4364, and He iiλ4687²). We picked the locations of these regions to cover a wide range in line ratios, surface brightness, and projected locations relative to the quasar. These regions are shown in Panel (g) of Fig. 4 and labelled with letters and numbers where S# refers to regions with higher surface brightness for which we used an extraction radius of 0.7 arcsec while B# labels low surface brightness regions which required a larger extraction radius (>1 arcsec) to achieve sufficient S/N.

To measure the emission properties for each region, we jointly fit the strong and weak emission lines described above with Gaussian profiles using LMFIT (Newville et al. 2014). For each region, all fitted lines share the same redshift and velocity width, but line fluxes are free parameters except for cases with line ratios set by atomic physics (e.g. [O iii]λ4960 and [O iii]λ5008). In most cases, a single set of Gaussians is enough to describe the emission-line profiles, except for S3, S4, and B4 which require a second set of Gaussians to account for broader (⁠|$\sigma \approx 100{-}170 \, \rm km \, s^{-1}$|⁠) emission wings. Such emission wings are often seen around luminous quasars due to quasar-driven outflows (Heckman et al. 1981; Liu et al. 2013a, b), but the wings on S3, S4, and B4 may also be due to projection effects. We summarize the measurements for these regions, including their distances from the quasar, extraction radii, line fluxes, LOS velocities, and 1D velocity dispersions, in Table 2. We display strong and weak-line spectra as well as their best-fitting models in Figs 5 and 6, respectively, for a representative subset of the regions.

4 DISCUSSION

As discussed in Section 3.3, the environment of HE 0238−1904 is overdense and includes a massive galaxy group or cluster. Based on clustering studies, this environment is richer than those of most radio-quiet systems, but consistent with expectation for radio-loud ones. This demonstrates that radio-quiet systems like HE 0238−1904 are diverse in terms of their host environment. Nevertheless, the lack of detected radio emission and amorphous morphology of the nebula suggests that it is not jet related. Considering that most published giant nebulae at z < 1 are in rich environments, the presence of giant nebulae might be correlated with group properties. A larger sample size of quasars with wide IFS observations is required to investigate this possibility.

Alternatively, such a rich environment can be explained by variable radio quasars. Quasars are capable of changing from radio-quiet to radio-loud or vice versa. Nyland et al. (2020) found 26 sources showing radio variability over time-scales of decades from the SDSS DR14 quasar catalogue (Pâris et al. 2018) and the Wide-field Infrared Survey Explorer (Wright et al. 2010) R90 quasar catalogue (Assef et al. 2018). These sources, once considered radio-quiet quasars, now meet the criteria for radio-loud ones. It implies that the probability that any particular radio-quiet quasar becomes radio-loud on the light-crossing time-scale of the nebula is approximately 1 per cent. However, the presence of a massive group and nebula mean that HE 0238−1904 is not a representative quasar and so may be more likely to transition to radio-loud relatively soon. On the other hand, the possibility that HE 0238−1904 was previously radio-loud and is now radio-quiet is harder to address since such transitions are not well studied.

In the following subsections, we discuss insights into the physical origins and state of the giant nebula which includes analyses of density and ionization-state sensitive diagnostic emission lines. Several of these analyses require priors on the dust content and density of the gas. To investigate dust content, we estimate Balmer-line ratios, and find |$\rm H\delta /H\gamma$| ratios of ≈0.55. These ratios are consistent with Case B recombination (Osterbrock & Ferland 2006) in the absence of dust. To obtain density estimates, we infer emission measure of the nebula from the surface brightness of |$\rm H \beta$| following Chen et al. (2019). Assuming |$\rm H\alpha /H\beta \approx 3$|⁠, a clumping factor of 1, and length-scale |$30 \rm \, pkpc$|⁠, we found an electron density of log (n_e/cm⁻³) ≈ −1. However, this density estimate has a large uncertainty and is effectively a lower limit due to the assumption of a unity clumping factor.

4.1 Origin of the nebular gas

Giant nebulae can be produced via ram pressure and tidal stripping, AGN and stellar feedback, or filamentary accretion. The nebula around HE 0238−1904 is unlikely to arise from a jet-driven outflow given the fact that the quasar is radio-quiet and exhibits no detectable radio jet. While S3 and S4 exhibit broad emission wings, most regions are well characterized by a single Gaussian profile with narrow velocity dispersion (⁠|$\sigma \lt 120 \, \rm km \, s^{-1}$|⁠; see Table 2). These quiescent kinematics are inconsistent with the broad velocity dispersion expected from radio-quiet AGN and stellar feedback (Liu et al. 2013b; Rupke et al. 2019). In addition, the morphology is inconsistent with expectations for filamentary accretion (Johnson et al. 2022). On the other hand, the nebula is spatially and kinematically coincident with likely interacting galaxies in the field of HE 0238−1904, suggesting that stripping from interactions is likely responsible for most of the nebula with possible subdominant contributions from outflows.

The nebula spatially surrounds the Host, G1, G3, G4, and G5, and extends to the south-west of the quasar to a projected distance of |$d \sim \rm 70 \, pkpc$|⁠. This spatial coincidence suggests that the nebula likely arises from interaction-related stripping. The dwarf galaxies G3 and G5 show a possible tidal-tail-like structure as shown in panels (e) and (h) of Fig. 4, suggesting that this part of the nebula might be created from tidal stripping. In addition to this, the emission maps on larger scales resemble a head-tail morphology with the head around the quasar and with the tail extending to the south-west of the quasar. Head-tail morphologies are commonly seen in nebulae originated from ram pressure-stripped ISM (e.g. Poggianti et al. 2016; Boselli et al. 2019; Chen et al. 2019). Interestingly, while the nebula exhibits a head-tail morphology, it does not exhibit multiple filaments like some ‘jellyfish’ galaxies observed in the optical line emission. Instead, it resembles the smoother emission profile sometimes seen in ram pressure debris observed in H i 21 cm (Hess et al. 2017). There are two plausible explanations for ram pressure stripping in the environment of HE 0238−1904. First, the nebula may arise from stripping of the quasar host’s ISM and CGM if it is falling into the richer, redshifted group and passing through the associated hot halo. Secondly, dwarf galaxies may have travelled through the hot halo of the massive group from west to east, leaving their ram pressure-stripped ISM and CGM behind along their path.

The discovery of a giant nebula requires both the presence of gas and its positioning within quasar’s ionization cone. However, due to projection effects, the relative position between the quasar and the nebula remains uncertain. The two previously mentioned hypotheses provide potential frameworks. (1) If the gas results from stripping of the quasar host’s ISM, the nebula is likely to surround the quasar. In this case, it will naturally be illuminated by the quasar. Alternatively (2) if the nebula arises from the stripped CGM/ISM of other galaxies in the overdensity, the gas will be widely distributed throughout the groups and more distant from the quasar. Only a fraction of this gas might coincidentally fall within the quasar’s ionization cone, consistent with the large opening angle suggested by Trainor & Steidel (2013), Borisova et al. (2016a), Schmidt et al. (2018), and den Brok et al. (2020).

To examine between these scenarios, we show the surface brightness profiles of [O ii] and [O  iii] made with Photutils (Bradley 2023) in Fig. 7. The profile of [O ii] declines smoothly as a function of radius, and plateaus at |$\approx \rm 50 \, pkpc$|⁠. In contrast, the [O iii] profile exhibits shallower drop due to the bright knots seen in the narrow-band images. The plateau in the [O ii] profile corresponds to the head-tail morphology of the nebula, and the bright knots hints at a dwarf-related origin for part of the nebula. Collectively, the [O ii] and [O iii] profiles suggest a complex scenario. The centroids of narrow-band [O iii] and [O iii] surface brightness maps are 10 and 19 pkpc away from the quasar, respectively, an alignment to within 15 per cent of the size of the nebula. This coincidence could be explained if the gas surrounds the quasar or if the quasars ionization cone is fairly well centred on our LOS. However, the significant contributions of individual dwarf galaxies to the [O iii] surface brightness profile underscore the challenge in precisely determining the nebula’s position relative to the quasar. Consequently, it is plausible that both scenarios (1) and (2) contribute to the nebula.

The giant nebulae around HE 0238−1904 were independently discovered and reported by Zhao & Wang (2023). They attributed the gas to a superbubble driven by the quasar based on an apparent large velocity shift between the nebula and the quasar redshift and as well as broad-line widths reported near the quasar. However, the large velocity shift is due to the reliance on an older, Mg ii-based redshift of z = 0.631, which is |$\approx +500\ {\rm km\, s^{-1}}$| from our [O ii]-based redshift of z = 0.6282. Rather than relying on a redshift estimate from the literature, we measured the quasar redshift and kinematics of the giant nebula from the same MUSE data set to avoid any systematic uncertainty due to wavelength calibration errors. Moreover, quasar redshifts based on [O ii] are generally more accurate than those measured from Mg ii due to the narrowness of the line and lack of blueshifted wings on [O ii]. In particular, quasar redshifts measured from [O ii] trace the underlying quasar host redshifts measured in stellar absorption to within |$\approx \pm 20\ {\rm km\, s^{-1}}$| (Hewett & Wild 2010). Finally, our redshift estimate of z = 0.6282 is more consistent with the centroid of the broad Hβ line, aligns with the peak of the quasar’s [O iii] emission line, and matches a more recent Mg ii-based redshift of z = 0.628 from the UV-bright Quasar Survey (Monroe et al. 2016). Furthermore, we measured significantly narrower line widths near the quasar. This is likely due to our removal of [O iii] and [O ii] emission from the unresolved narrow-line emission region of the quasar while Zhao & Wang (2023) only removed emission from the broad-line region. In summary, the modest velocity shifts and largely narrow emission-line widths are consistent with much of the gas originating from interactions with more minor possible contributions from an outflow. When using the updated quasar redshift and quasar-light subtracted datacube, we find no evidence for a fast, quasar-driven superbubble in the system.

4.2 Physical conditions of the emitting gas

Previous studies of giant nebulae have attributed the ionization of the gas to ionizing photons from AGN, shocks, and young stellar populations (e.g. Johnson et al. 2018; Chen et al. 2019; Rupke et al. 2019; Helton et al. 2021; Zhang et al. 2023). The presence of the quasar suggests the source of ionization is AGN related. To study the physical conditions of the gas, we measured the density- and temperature-sensitive |$\rm [O\,{\small II}]\lambda 3729/[O\,{\small II}]\lambda 3727$| and |$\rm [\rm O\,{\small III}]\lambda 4364/[\rm O\,{\small III}]\lambda 5008$| line ratios as well as ionization state-sensitive strong and weak line ratios in each region. These line ratio measurements are reported in Table 2 and a [O iii]/[O ii] map is shown in panel (c) of Fig. 4. We discuss these measurements and their implications in the following three subsections.

4.2.1 Direct density and temperature estimates

With spectral coverage of |$\rm [O\,{\small II}]\lambda 3727$|⁠, |$\rm [O\,{\small II}]\lambda 3729$|⁠, |$\rm [\rm O\,{\small III}]\lambda 4364$|⁠, and |$\rm [\rm O\,{\small III}]\lambda 5008$|⁠, we can directly measure electron density (n_e) and temperature (T_e), as discussed in Osterbrock & Ferland (2006). The |$\rm [O\,{\small II}]$| doublet is a good density estimator because the difference in excitation energy between these two upper states is small so that the relative population in the two states is determined by electron density and is insensitive to temperature. In contrast, the |$\rm [\rm O\,{\small III}]$| doublet upper states have a larger excitation energy difference, making the populations of these states mainly sensitive to electron temperature and insensitive to electron density. Electron number densities from the |$\rm [O\,{\small II}]$| doublet are reasonable proxies for the overall densities of ionized nebulae because |$\rm H$| and |$\rm O$| share similar ionization energies of |$13.6 \, \rm eV$|⁠.

To translate line ratios into physical conditions, we used Pyneb (Luridiana, Morisset & Shaw 2015) which predicts the [O ii] and [O iii] line ratios at a given density and temperature by solving the detailed balance equation for an n-level atom. We fit the measured line ratios with Pyneb models by performing Markov chain Monte Carlo (MCMC) analysis with emcee (Foreman-Mackey et al. 2013), and inferred physical conditions from the resulting posteriors. We report the densities in Table 3, though we omit measurements in cases where the S/N or broad-line width results in poorly constrained conditions.

For all regions, where the |$\rm [O\,{\small II}]$| doublet is resolved, the line ratio is in the low-density limit except S6. We therefore report 95  per cent upper limits in density for all but S6. The inferred electron number density upper limits range from |$1.2 \lt \log (n_{\rm e, [O\,{\small II}]}/\mathrm{cm^{-3}}) \lt 2.8$|⁠, with a median of |$\log (n_{\rm e, [O\,{\small II}]}/\mathrm{cm^{-3}}) \lt 1.6$|⁠. These density upper limits are consistent with gas arising from ionized ISM (Draine 2011) or CGM. We detected |$\rm [\rm O\,{\small III}]\lambda 4364$| in only three luminous regions, S1, S2, and S6. The inferred temperatures for S1, S2, and S6 are log (T/K) ≈ 4.2, 4.2, and 4.1, respectively.

4.2.2 Indirect density estimates from photoionization simulations

Under the assumption that the nebula is ionized by the quasar, its ionization states are set by the luminosity of the quasar, density of the gas, and distance from the quasar, with secondary effects from metallicity and ionizing spectral shape. With an estimate of the quasar’s luminosity and assuming projection effects are negligible, the density structure of the gas can be inferred from measured line ratios (see Cantalupo et al. 2019). Studies of high-redshift quasar nebulae found ionization states can only be explained by a density of log (n_H/cm⁻³) ≈ 1.9, significantly higher than expected CGM/IGM densities, or alternatively by a broad density distribution (see Cantalupo et al. 2019). At low redshift, this kind of scenario can be further explored with insight from rest-optical lines to compare ionization-based densities with more direct density estimates from the |$\rm [O\,{\small II}]$| doublet.

To infer the physical conditions from the line ratios in Table 2, we ran photoionization simulations for each region with Cloudy version C17.03 (Ferland et al. 2017). We modelled the quasar’s radiation field using a power law (I ∝ ν^α) between 0.37 and 73.5 |$\rm Ryd$|⁠, with α between −1.8 < α < 0 following Groves, Dopita & Sutherland (2004) but extending to a higher α. We set the modelled quasar luminosity at |$1\, \rm Ryd$| using direct measurement of the monochromatic UV luminosity from COS. For the gas, we adopted single density and single metallicity models, with density of −2 < log (n_H/cm⁻³) < 4.6 and metallicity of −1.5 < log (Z/Z_⊙) < 0.5. We chose this metallicity range to cover the characteristic metallicities of the cool CGM around massive elliptical galaxies (Zahedy et al. 2019) but extended it to higher metallicity in case some gas has ISM origins. Due to limited ion coverage, metallicity, and α are degenerate in some cases, so we treated them as nuisance parameters and focused on inferred densities. We note that there is relatively little degeneracy between density and metallicity except at high metallicities of log (Z/Z_⊙) > 0.2 when increased cooling from metal lines begins to substantially change the equilibrium temperature.

For each region, we conducted these models in grids with a step of 0.2 dex in density and metallicity, and 0.2 in α. We then interpolated these models with the RegularGridInterpolator function from scipy.interpolate (Virtanen et al. 2020) within these ranges after checking for convergence. Finally, we ran emcee to estimate posteriors given the measured line ratios and uncertainties. We verified the quality of the fits by comparing the posteriors of the model line ratios with the measured line ratios using violin plots shown in Fig. 9. The violin plots verify that the ionization-state-sensitive line ratios (shown in the middle panels) are consistent with the measured line ratios. The best-fitting α values for most regions are within −1.0 < α < −0.6, somewhat greater than ones given in Groves, Dopita & Sutherland (2004). Inferred metallicities for S1, S2, and S6, with He ii and [Ne v] detections, are well constrained to be −0.2 < log (Z/Z_⊙) < 0.2. The densities inferred from these photoionization simulations range from log (n_{H, Cloudy}/cm⁻³) = 1.6 to 4.2 and are reported in the right column of Table 3, though we stress that these densities neglect potential quasar variability and projection effects.

4.2.3 Comparison of the density estimates

Previous photoionization-based estimates of the density of quasar nebulae at high redshift found unexpectedly high densities, close to or exceeding typical densities for the ISM, despite being measured on CGM/IGM scale (Cantalupo et al. 2019). The ionization sensitive line ratios of the nebula around HE 0238−1904 also imply high photoionization-based densities of |$1.6\lt \log (n_{\rm H, \, Cloudy}/\mathrm{cm}^{-3})\lt 4.2$|⁠. However, the more direct |$\rm [O\,{\small II}]$|-based densities are inconsistent with and significantly smaller than the photoionization-based densities for most regions, as shown in Table 3. To better demonstrate this inconsistency, Fig. 9 shows both the measured line ratios and the posteriors inferred from the photoionization models for S2, S6, and S9. The ionization-state-sensitive line ratios are consistent with the model posteriors for all three regions, while the |$\rm [O\,{\small II}]$| line ratios are highly discrepant for S6 and S9. The right panel of each subfigure shows the density posteriors from both direct and indirect density estimates.

As shown in Table 3, we found that all regions with photoionization-based density estimates except S1, S2, B1, and B3 have a large (1−2 dex) discrepancy when compared to the [O ii] doublet-based densities. In the most extreme case, S5, the two density estimates are off by 2.6 dex or a factor of 400. In principle, the inferred density mismatch could be explained by a non-uniform density distribution if the [O ii] arises from less dense gas than the other emission lines. To test whether a more complicated density structure could explain the density mismatch, we modelled the emitting gas as a multiphase system consisting of one low-density component and one high-density component with the relative contribution of each treated as an additional free parameter. This model successfully reproduces the observed emission-line ratios, and the density inferred for the high-density component matches the single-phase model results. Furthermore, the posteriors of the two-component model indicate that the high-density component dominates the |$\rm [O\,{\small II}]$| emission. Therefore, a two-phase model cannot explain the density discrepancy between the direct [O ii]-based density measurements and the ionization-state-based density estimates.

To test if a broad, continuous density distribution can explain the discrepancy, we modelled the emitting gas with a lognormal density distribution (see Cantalupo et al. 2019). A lognormal distribution is defined as

where σ is the dispersion and μ is the mean density. We started with calculating emission-line emissivity in an extended Cloudy model grid, similar to ones discussed in Section 4.2.2. We then computed the predicted line ratios for a lognormal density distribution by interpolating Cloudy models and integrating over the PDF. Our results show that a lognormal distribution with a large σ can reproduce the ionization-sensitive line ratios, but the lognormal models predict that the [O  ii] emission arises from dense gas, resulting in [O ii] line ratios of |$\log (\frac{\lambda 3729}{\lambda 3727}) = {-}0.4 \, {\rm to} \, {-}0.1$|⁠, inconsistent with the observed ratios of |$\log (\frac{\lambda 3729}{\lambda 3727}) \gt 0.1$|⁠. Therefore, a broad density distribution is unlikely to reconcile the density discrepancy.

Alternatively, projection effects can also result in disagreement between the two density estimates. However, assuming that the gas is randomly and approximately spherically distributed around the quasar, the projected distance is unlikely to be much smaller than the radial distance between the quasar and the nebula. For example, producing a factor of 400 mismatch in density requires the radial distance to be 20 times larger than the projected distance. While such projection effects are possible, in principle, the required contrived geometry is unlikely.

In principle, the discrepancy in density could be explained if the nebula is not directly ionized by the quasar due to obscuring dust or translucent clouds blocking its light from reaching this gas. Filtering the quasar’s radiation through dust would soften the incident ionizing radiation field. However, the best-fitting α values from our photoionization analysis suggest a hard ionizing spectrum for almost all regions. The hard inferred ionizing spectrum is inconsistent with expectations from a quasar SED filtered through dust clouds.

Alternatively, translucent clouds of moderate optical thickness to ionizing photons can also filter the quasar’s radiation. Depending on the density and the physical size, these clouds could produce distinct line ratios as a function of depth into the cloud (Liu et al. 2013a). Typically, the outer parts of the cloud produce no significant [O ii] or [O iii] emission because oxygen is highly ionized. However, |$\rm H\beta$| is a recombination line and so a non-negligible fraction of the |$\rm H\beta$| emission arises from outer parts of the cloud that do not emit in [O ii] or [O iii]. As a result, translucent regions are expected to have stronger |$\rm H\beta$| emission than |$\rm [O\,{\small II}]$| and |$\rm [\rm O\,{\small III}]$|⁠. Yet, none of the nebular regions have such |$\rm [\rm O\,{\small III}]/H\beta$| ratio. If these translucent clouds exist around HE 0238−1904, they therefore must be blended with optically thick clouds due to seeing conditions and projection effects. The presence of unresolved translucent clouds could be investigated by observing the nebula with higher spatial resolution instruments such as NIRSpec on the JWST or with adaptive optics from the ground. Nevertheless, while translucent clouds may help reconcile the density discrepancy in some cases, moderate optical depth clouds can only absorb a modest portion of the quasar’s radiation. Therefore, it is unlikely to explain the largest density discrepancies.

On the other hand, the ionization of the nebulae could be due to young stellar populations (Morisset, Delgado-Inglada & Flores-Fajardo 2015) or fast shocks (Allen et al. 2008). However, there is no evidence of extended star-formation in rest-frame u-band images of the system formed from the MUSE datacube. To investigate the possibility of fast shocks, we show two emission-line diagnostic diagrams overlaid with shock models in a grid of shock velocity and magnetic field strength in Fig. 8. Producing the observed [O iii]/[O ii] and [Ne v]/[O ii]³ ratios requires shock velocities of |$v_{\rm shock}\gt 250 \, \rm km \, s^{-1}$| (Allen et al. 2008). These shock velocities are greater than the LOS velocity and velocity dispersion of the nebula in nearly all locations, even after accounting for projection effects. For example, some regions (S1 and S2) would require shock velocities exceeding |$1000 \, \rm km \, s^{-1}$| and most regions (S3, S4, S6, S8, S10, B1, B2, B3, and B4) would require |$\gt 300\!-\!400 \, \rm km \, s^{-1}$|⁠, making them unlikely to be ionized by shocks. On the other hand, while the observed line ratios of S5, S7, and S9 favour AGN photoionization, large uncertainties in their |$\rm H\beta$| flux can accommodate shocks with velocities as low as |$200 \, \rm km \, s^{-1}$|⁠. This would alleviate the density discrepancy in these three regions. However, for most regions, the shock velocity required to reproduce the observed line ratios exceeds velocities observed in the system. Shocks are therefore unlikely to explain the density discrepancy in most cases.

Perhaps more likely, the difference in the density estimates could be due to quasar variability (Richstone & Oke 1977). Quasar variability is directly observed on time-scales of decades (Stone et al. 2022). Observations of ‘changing-look’ AGN, light echoes, and quasar proximity zones suggest the average episodic lifetime of quasars may range from 10⁴ to 10⁷ yr and AGN episodes may be highly clustered (e.g. Schirber, Miralda-Escudé & McDonald 2004; Gonçalves, Steidel & Pettini 2008; Kirkman & Tytler 2008; Trainor & Steidel 2013; Syphers & Shull 2014; Schawinski et al. 2015; Comerford et al. 2017; Schmidt et al. 2018; Shen 2021). Therefore, each region of the nebula around HE 0238−1904 may experience a drastically different radiation field from the quasar, depending on the light travel time. For example, S5 and S6 are at a projected distance of 10 to |$\rm 20 \, pkpc$| from the quasar, respectively, and their line ratios can be explained if the quasar was 400 and 10 times less luminous than currently observed. In contrast, S1 and S2 are at a projected distance of |$\approx \rm 40 \, pkpc$| from the quasar, and their properties can be explained if they received ionizing radiation consistent with the current luminosity of the quasar. We confirmed that quasar variability could explain the ionization state and |$\rm [O\,{\small II}]$| ratio by re-running Cloudy models and MCMC analysis after significantly decreasing the quasar luminosity.

5 SUMMARY AND CONCLUSIONS

In this paper, we presented the first comprehensive analysis of a giant nebula around a radio-quiet quasar at z < 1 based on MUSE observations of the field of HE 0238−1904. The wide FoV, high spatial sampling, and wide wavelength coverage enabled us to investigate the origin and the physical condition of the group and gaseous environment with a spatially resolved analysis of the morphologies, kinematics, and nebular photoionization properties. Our finding can be summarized as follows:

• We found that HE 0238−1904 resides in an overdense environment containing two potentially merging galaxy groups based on spatial distribution and kinematics. This includes a less rich, blueshifted group with 12 galaxies and a richer, redshifted group with 22 galaxies. Assuming the more massive group is virialized, its dynamical mass is |$M_{\rm dyn} \sim 4 \times 10^{13}{-}10^{14}\ {\rm {\rm M}_{\odot }}$|⁠. Such a massive, rich environment is unusual for a radio-quiet quasar, which typically resides in a halo with a mass of |$\sim 3 \times 10^{12}\ {\rm {\rm M}_{\odot }}$| (Shen et al. 2009).

• We identified a giant nebula covering a projected area of ≈5000 kpc² around HE 0238−1904 emitting strongly in |$\rm [O\,{\small II}]$|⁠, Hβ, and |$\rm [\rm O\,{\small III}]$|⁠. The nebula has an irregular morphology with a spatial trend in kinematics where the region north of the quasar is redshifted and the region south of the quasar is mainly blueshifted relative to the quasar. The southern region is spatially coincident with four dwarf galaxies.

• The coincidence with nearby galaxies suggests that it arises from stripping of ISM or CGM, which is consistent with its morphology and largely narrow LOS velocity dispersion. In addition, the nebula shows a head-tail morphology with the head near the quasar and with the tail extending toward south-west of the quasar. The head-tail structure may originate from ram pressure if the quasar and the surrounding nebula are infalling toward the massive galaxy group to the north-east. However, we note there are some small regions at |$d \approx 20 \, \rm pkpc$| from the quasar that have broader emission wings, perhaps suggesting an outflow origin.

• To better characterize the physical conditions of the nebula, we measured the fluxes of strong and weak emission-line fluxes. The inferred electron number density upper limits from the |$\rm [O\,{\small II}]$| doublet range from |$\log (n_{\rm e, [O\,{\small II}]}/\mathrm{cm^{-3}}) \lt 1.2$| to 2.8, with a median of |$\log (n_{\rm e, [O\,{\small II}]}/\mathrm{cm^{-3}}) \lt 1.6$|⁠. These density upper limits are consistent with ISM or CGM origin. However, densities inferred from photoionization models are often inconsistent with the |$\rm [O\,{\small II}]$|-based density upper limits, reaching values of up to 400 times higher.

• The disagreement in density estimates is unlikely to be due to density inhomogeneities, but can be explained by quasar variability, if the quasar varied significantly on time-scales of 10⁴ to 10⁵ yr. This finding suggest that long-term quasar variability should be included when considering ionization-based inferences into the physical conditions of giant nebulae around quasars.

The possibility of significant quasar variability on time-scales of 10⁴ to 10⁵ yr has implications far beyond accretion disc physics in the central engine. In particular, significant fluctuations on these time-scales can result in out-of-equilibrium conditions in the low-density circumgalactic medium due to the long recombination time of low-density gas (Oppenheimer & Schaye 2013; Segers et al. 2017). Indeed, such AGN ‘flickering’ may be responsible for strong O vi absorption observed around Milky Way-like galaxies at low redshift (Oppenheimer et al. 2018). The recent and upcoming commissioning of new IFSs on large telescopes, such as LLAMAS (Furesz et al. 2020), IFUM (Mateo et al. 2022), Blue MUSE (Richard 2019), and MIRMOS (Konidaris et al. 2020), will continue to drive further discoveries of giant nebulae which could be followed up with IFS like HARMONI (Thatte et al. 2022) on future, 30-m class telescopes, extending similar insights to higher redshifts and fainter systems.

ACKNOWLEDGEMENTS

SDJ and ZQL acknowledge partial support from HST-GO-15280.009-A, HST-GO-15298.007-A, HST-GO-15655.018-A, and HST-GO-15935.021-A. JIL is supported by the Eric and Wendy Schmidt AI in Science Postdoctoral Fellowship, a Schmidt Futures program. SC gratefully acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation programme grant agreement No 864361. This paper is based on observations from the European Organization for Astronomical Research in the Southern Hemisphere under ESO (PI: J. Schaye, PID: 094.A-0131(B) & 096.A-0222(A)), and the NASA/ESA Hubble Space Telescope (PI: L. Straka, PID: 14660; PI: J. Green, 11541; PI: S. Penton, PID: 12505). Additionally, this paper made use of the NASA/IPAC Extragalactic Database, the NASA Astrophysics Data System, Astropy (Astropy Collaboration et al. 2022), Aplpy (Robitaille & Bressert 2012), and Photutils (Bradley 2023).

DATA AVAILABILITY

The data used in this paper are available from the ESO and HST data archives.

Footnotes

References