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

Zhuoqi (Will) Liu1, Sean D. Johnson1, Jennifer I-Hsiu Li1,2, Gwen C. Rudie3,
Joop Schaye4, Hsiao-Wen Chen5, Jarle Brinchmann6, Sebastiano Cantalupo7, Mandy
C. Chen5, Wolfram Kollatschny8, Michael V. Maseda9, Nishant Mishra1, Sowgat
Muzahid10
1Department of Astronomy, University of Michigan, 1085 S. University, Ann
Arbor, MI 48109, USA
2Michigan Institute for Data Science, University of Michigan, Ann Arbor, MI,
48109, USA
3The Observatories of the Carnegie Institution for Science, 813 Santa Barbara
Street, Pasadena, CA 91101, USA
4Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, The
Netherlands
5Department of Astronomy & Astrophysics, The University of Chicago, Chicago,
IL 60637, USA
6Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP,
Rua das Estrelas, PT4150-762 Porto, Portugal
7Department of Physics, University of Milan Bicocca, Piazza della Scienza 3,
I-20126 Milano, Italy
8Institut für Astrophysik und Geophysik, Universität Göttingen, Friedrich-Hund
Platz 1, D-37077 Göttingen, Germany
9Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St.,
Madison, WI 53706 USA
10Inter-University Centre for Astronomy and Astrophysics (IUCAA), Post Bag 4,
Ganeshkhind, Pune 411 007, India E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present the first comprehensive study of a giant, $\approx\\!\\!70$ kpc-
scale nebula around a radio-quiet quasar at $z<1$. The analysis is based on
deep integral field spectroscopy with MUSE of the field of HE 0238$-$1904, a
luminous quasar at $z=0.6282$. The nebula emits strongly in
$\mathrm{[O\,II]}$, $\rm H\beta$, and $\mathrm{[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_{\rm dyn}\approx 4\times 10^{13}$ to
$10^{14}\ {\rm M_{\odot}}$. 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\,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 timescale of $\approx 10^{4}{-}10^{5}$ years.

###### keywords:

quasars: supermassive black holes – galaxies: groups – intergalactic medium

††pubyear: 2023††pagerange: The first comprehensive study of a giant nebula
around a radio-quiet quasar in the $z<1$ Universe–The first comprehensive
study of a giant nebula around a radio-quiet quasar in the $z<1$ Universe

##  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 timescale 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 neighboring 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 (KCWI; Martin et al. 2010). At $z>2$, systematic IFS surveys
around radio-quiet quasars discovered ubiquitous giant H I Ly$\alpha$ 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
$\rm[O\,II]$, $\rm H\beta$, and $\rm[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; Zabl et al., 2021; Burchett 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 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 dataset 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 $\Lambda$ cosmology with $\Omega_{\rm
m}=0.3$, $\Omega_{\rm\Lambda}=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\approx 0.63$ quasar HE 0238$-$1904 has high-quality archival UV HST
absorption spectra used to study the CGM of the Milky Way (Zheng et al., 2019;
Bish et al., 2021) and distant galaxies (Muzahid et al., 2018; Lehner 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) &
096.A-0222(A)). MUSE is an integral-field spectrograph on the UT4 VLT with a
field of view (FoV) of $1^{\prime}\times 1^{\prime}$ and a spaxel size of
$0.2^{\prime\prime}$ in wide-field mode (WFM). MUSE covers the spectral range
between $4750\,\text{\AA}$ to $9350\,\text{\AA}$ and a resolution of $R\sim
3000$. The MUSE observations are centered near the quasar sightline, and we
obtained eleven exposures collected between November 18th, 2014 and February
2nd, 2016 with a total exposure time of 8.75 hr with median seeing full-width-
at-half-maximum (FWHM) conditions of $0.7^{\prime\prime}$. At the redshift of
HE 0238$-$1904, the MUSE FoV corresponds to a projected size of $\approx 400$
proper kpc (pkpc) on a side, and the spectral coverage includes emission lines
such as [O II], H$\beta$, and [O III]. These emission lines enable sensitive
studies of any ionized nebulae and galaxies in the quasar’s environment.

Figure 1: MUSE spectrum of HE 0238$-$1904 overplotted with best-fit models.
The MUSE spectrum is shown as a solid black line, the power-law continuum
model is shown as a dashed purple line, and the iron template model is shown
using a solid blue line. The bottom left inset panel shows the
$\mathrm{[O\,II]}$ line emission with the best-fit continuum+line model shown
in red. The top right inset panel shows the $\rm H\beta$ and [O III] emission
with the best-fit shown in red.

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 (ZAP; 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 Hubble Space Telescope (HST) with
the F814W filter (PI: L. Straka, PID: 14660) with a total exposure time of
$2182$ seconds 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 14400
seconds and 7496 seconds in the G130M and G160M gratings, respectively (PI: J.
Green and S. Penton, PID: 11541 and 12505). We reduced and coadded the COS
spectrum following procedures outlined in Johnson et al. (2015); 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 et al., 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^{\prime\prime}$ aperture. To measure the systemic redshift of the quasar,
we fit the $\rm[O\,II]\lambda\lambda 3727,3729$ doublet with a Gaussian
profile following Hewett & Wild (2010) and found $z=0.6282\pm 0.0002$, where
the uncertainty represents the scatter between the $\rm[O\,II]$ centroid and
stellar absorption lines of SDSS quasars at similar redshift. This redshift is
$\approx+500\ {\rm km\,s^{-1}}$ 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, $\mathrm{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$\beta$, and two for the
narrow H$\beta$ 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 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 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\sigma$ below the
mean single-epoch virial theorem estimate, and the stellar mass is $1\sigma$
below the estimate from the black hole mass-stellar mass relation, the
inferred stellar mass would be $M_{*}\approx 10^{11.4}\ {\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 disk
wind, erroneously inflating the black hole mass estimate. The fitted quasar
spectrum is shown in Figure 1.

Figure 2: HST ACS+F814W image of the field of HE 0238-1904. The full image has
a FoV of $1.5^{\prime}\times 1.5^{\prime}$. The larger dashed box shows the
$1^{\prime}\times 1^{\prime}$ MUSE FoV. The smaller dashed box marks the
$30^{\prime\prime}\times 30^{\prime\prime}$ region displayed in Figure 4. The
LOS velocities of galaxies relative to the quasar are denoted with outlining
colors and the corresponding colorbar is shown on the bottom left. The
histogram in the bottom right inset panel shows the velocity distribution of
galaxies where galaxies in both orange and purple outlined regions are plotted
separately. We note that the orange and purple regions and corresponding
histograms are only for visualization. The two-Gaussian fitting of the
velocity distribution does not rely on any spatial information. Galaxies in
the quasar host environment are labeled with black circles and labeled by
their IDs. The approximate stellar mass weighted group center is marked with a
white asterisk while the weighted centers of the richer, redshifted group and
less rich, blueshifted group are marked with red and blue asterisks,
respectively.

###  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 catalog is complete for objects brighter than $m_{\rm
F814W}\approx 26\\!-\\!27$, depending on angular size.

Figure 3: MUSE galaxy spectra with the best-fit spectral models. The MUSE
spectrum is shown by a solid black line. The uncertainty is shown by a solid
grey line. The best-fit model used for redshift measurement is shown by a
solid red line.

For each identified object, we extracted a MUSE spectrum with MPDAF with a
circular aperture of $r=0.7^{\prime\prime}$, 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 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-fit linear combination on a grid
from $z=0$ to $z=1$ with a step size of $\Delta z=0.0001$ and recorded the
goodness-of-fit statistic ($\chi^{2}$) over the entire grid. We adopted the
redshift with the minimum global $\chi^{2}$ as our initial solution. We then
visually inspected each best-fit 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 catalog and the corresponding
redshift measurements, the redshift survey is approximately $100\%$ complete
for sources brighter than $m_{\rm F814W}\\!\approx\\!24$ and approximately
$95\%$ complete for those brighter than $m_{\rm F814W}\\!\approx\\!25$. For
comparison, an $L_{*}$ galaxy at $z\approx 0.6$ has $m_{\rm
F814W}\\!\approx\\!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.

Table 1: Summary of Galaxies in the Field of HE 0238$-$1904 at $z\approx
z_{\rm{QSO}}$.

ID | R.A.a | Decl.b | zc |  $m_{\mathrm{F814W}}$d |  $M_{B}$e |  K-correction | D4000 | $A_{V}$ |  $\log(M_{*}/{\rm M_{\odot}})$f |  $\Delta\theta$g | dh |  $\Delta$vi
---|---|---|---|---|---|---|---|---|---|---|---|---
| (J2000) | (J2000) |  | (AB) | (AB) | template |  | (mag) |  | (′′) | (pkpc) | ($\mathrm{km\ s^{-1}}$)
Host | 02:40:32.58 |  $-$18:51:51.4 | 0.6282 | … | … | … | … | … | … | 0.0 | 0.0 | 0
G1 | 02:40:32.63 |  $-$18:51:55.8 | 0.6278 | 24.3 |  $-$17.5 | S0 | $1.26\pm 0.57$ | … | 9.3j | 4.4 | 30.4 | -76
G2 | 02:40:32.73 |  $-$18:51:47.1 | 0.6270 | 23.3 |  $-$18.5 | S0 | $1.56\pm 0.08$ | 0.1 | 9.5 | 4.8 | 32.7 | -224
G3k | 02:40:32.74 |  $-$18:51:55.9 | 0.6280 | 23.8 |  $-$18.3 | Irr | … | … | 9.6j | 5.0 | 34.3 | -40
G4 | 02:40:32.57 |  $-$18:51:56.7 | 0.6284 | 24.9 |  $-$17.3 | Irr | $1.05\pm 0.07$ | 0.2 | 8.3 | 5.4 | 36.7 | +34
G5 | 02:40:32.71 |  $-$18:51:57.0 | 0.6280 | 25.2 |  $-$17.0 | Irr | $0.64\pm 0.08$ | 0.1 | 7.4 | 5.9 | 40.1 | -40
G6 | 02:40:32.96 |  $-$18:51:54.4 | 0.6295 | 22.4 |  $-$19.4 | S0 | $1.35\pm 0.02$ | 0.1 | 10.1 | 6.1 | 41.5 | +237
G7 | 02:40:33.04 |  $-$18:51:53.8 | 0.6275 | 23.8 |  $-$18.0 | S0 | $1.30\pm 0.04$ | 0.0 | 9.3 | 6.9 | 46.9 | -132
G8 | 02:40:32.21 |  $-$18:51:58.7 | 0.6284 | 21.8 |  $-$20.0 | S0 | $1.62\pm 0.02$ | 0.2 | 10.4 | 9.1 | 61.9 | +34
G9 | 02:40:33.44 |  $-$18:51:50.7 | 0.6330 | 23.8 |  $-$18.1 | S0 | $1.49\pm 0.05$ | 0.2 | 9.7 | 12.2 | 82.2 | +882
G10 | 02:40:33.53 |  $-$18:51:48.4 | 0.6323 | 20.0 |  $-$21.9 | S0 | $1.71\pm 0.01$ | 0.8 | 11.5 | 13.8 | 94.3 | +753
G11 | 02:40:32.37 |  $-$18:51:37.6 | 0.6302 | … | … | … | … | … | … | 14.1 | 96.3 | +360
G12 | 02:40:32.00 |  $-$18:51:39.9 | 0.6297 | 21.4 |  $-$20.4 | S0 | $1.64\pm 0.02$ | 0.2 | 10.6 | 14.1 | 96.5 | +274
G13 | 02:40:32.28 |  $-$18:52:04.9 | 0.6272 | … | … | … | … | … | … | 14.2 | 97.0 | -187
G14 | 02:40:33.17 |  $-$18:51:37.9 | 0.6310 | 22.6 |  $-$19.2 | S0 | $1.37\pm 0.03$ | 0.7 | 10.0 | 15.8 | 108.0 | +513
G15 | 02:40:33.62 |  $-$18:51:43.2 | 0.6253 | 24.8 |  $-$17.0 | S0 | $1.99\pm 0.22$ | 0.4 | 9.0 | 16.8 | 115.0 | -537
G16 | 02:40:31.85 |  $-$18:52:05.5 | 0.6279 | 23.8 |  $-$18.0 | S0 | $1.98\pm 0.16$ | 1.1 | 9.5 | 17.5 | 119.8 | -58
G17 | 02:40:33.75 |  $-$18:51:45.5 | 0.6332 | 22.7 |  $-$19.1 | S0 | $1.57\pm 0.03$ | 0.6 | 10.1 | 17.6 | 120.3 | +919
G18 | 02:40:33.53 |  $-$18:51:39.6 | 0.6332 | … | … | … | … | … | … | 17.9 | 121.9 | +922
G19 | 02:40:33.69 |  $-$18:52:00.1 | 0.6358 | 22.2 |  $-$19.7 | S0 | $1.60\pm 0.02$ | 0.4 | 10.3 | 18.0 | 122.9 | +1398
G20 | 02:40:31.97 |  $-$18:52:07.9 | 0.6271 | … | … | … | … | … | … | 18.8 | 128.1 | -205
G21 | 02:40:33.48 |  $-$18:51:36.9 | 0.6341 | 22.1 |  $-$19.7 | S0 | $1.26\pm 0.02$ | 1.4 | 10.3 | 19.3 | 131.8 | +1084
G22 | 02:40:31.34 |  $-$18:52:02.5 | 0.6268 | 23.0 |  $-$18.9 | S0 | $1.66\pm 0.05$ | 0.5 | 10.1 | 20.9 | 142.8 | -261
G23 | 02:40:33.76 |  $-$18:51:38.2 | 0.6319 | 24.4 |  $-$17.6 | S0 | $1.62\pm 0.11$ | 0.6 | 9.5 | 21.3 | 145.5 | +679
G24 | 02:40:33.87 |  $-$18:51:36.1 | 0.6333 | 23.6 |  $-$18.5 | Scd | $1.07\pm 0.04$ | 1.8 | 9.8 | 23.8 | 162.4 | +937
G25 | 02:40:33.26 |  $-$18:52:13.9 | 0.6277 | 25.5 |  $-$16.3 | S0 | $1.46\pm 0.17$ | 1.8 | 8.8 | 24.5 | 167.5 | -95
G26 | 02:40:30.93 |  $-$18:51:43.7 | 0.6272 | 23.0 |  $-$18.8 | S0 | $1.66\pm 0.05$ | 0.6 | 9.9 | 24.7 | 168.3 | -187
G27 | 02:40:34.29 |  $-$18:51:46.3 | 0.6297 | 23.7 |  $-$18.1 | S0 | $1.30\pm 0.06$ | 0.5 | 9.5 | 24.8 | 169.0 | +274
G28 | 02:40:32.96 |  $-$18:52:17.2 | 0.6282 | 23.0 |  $-$19.1 | Scd | $1.07\pm 0.02$ | 1.0 | 9.5 | 26.4 | 180.0 | -3
G29 | 02:40:32.32 |  $-$18:51:24.4 | 0.6357 | 24.5 |  $-$17.8 | Irr | $0.83\pm 0.06$ | 1.4 | 8.3 | 27.2 | 185.6 | +1379
G30 | 02:40:34.59 |  $-$18:51:45.3 | 0.6323 | 24.5 |  $-$17.6 | Scd | $0.96\pm 0.06$ | 0.1 | 8.9 | 29.2 | 199.2 | +753
G31 | 02:40:34.57 |  $-$18:52:00.4 | 0.6312 | … | … | … | … | … | … | 29.6 | 201.9 | +550
G32 | 02:40:34.83 |  $-$18:51:55.7 | 0.6354 | 24.8 |  $-$17.1 | S0 | $1.25\pm 0.08$ | 0.0 | 9.0 | 32.2 | 219.9 | +1324
G33 | 02:40:34.55 |  $-$18:51:34.9 | 0.6313 | 20.0 |  $-$22.1 | Scd | $1.17\pm 0.01$ | 0.4 | 10.8 | 32.4 | 220.9 | +569
G34 | 02:40:34.88 |  $-$18:51:53.0 | 0.6349 | 20.6 |  $-$21.2 | S0 | $1.55\pm 0.01$ | 0.4 | 11.2 | 32.7 | 222.9 | +1232

  * •

Notes.

  * a

Right ascension.

  * b

Declination.

  * c

Best-fit redshift, from principal component analysis of SDSS galaxy
eigenspectra from BOSS. G11/G18 have only one spectral feature.

  * d

Apparent HST ACS+F814W magnitude.

  * e

Absolute B-band magnitude.

  * f

Stellar mass from stellar population fits to the MUSE spectrum and DES & HST
photometry.

  * g

Angular distance from the quasar.

  * h

Projected physical distance from the quasar.

  * i

LOS velocity from the quasar.

  * j

Stellar mass estimated from the median $M_{*}/L$ ratio of the group resulting
in large systematic uncertainties.

  * k

The uncertainty in the position of G3 is larger than other galaxies due to the
diffraction spike in the HST ACS+F814W image.

To examine properties of the quasar host environment, we identified candidate
group members based on their LOS velocities relative to the quasar ($\Delta
v=v-v_{\rm QSO}$). In particular, we selected galaxies with $|\Delta v|<\rm
2000\ km\,s^{-1}$. We inferred the physical properties of the selected
galaxies with Bagpipes (Carnall et al., 2018; Carnall et al., 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 (Buchner et al., 2014; Feroz et
al., 2009; 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 datapoint available, necessitating the use of the spectra to
further inform the stellar population synthesis. In our fitting procedure, we
assumed an exponential star formation history with e-folding time scale of
$0.01<\rm\tau/Gyr<8.00$, solar stellar metallicity, and dust attenuation model
from Calzetti et al. (2000) with $0<A_{V}/\rm mag<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 2nd 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 $\approx 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 et al. (1980) chosen based on the strength of the 4000 Å break. The
IDs, galaxy coordinates (R.A., Decl.), 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 Figure 2
and several example MUSE spectra are overplotted with their best-fit PCA
spectral models in Figure 3. An interactive view of the galaxy environment and
spectra is available online111http://zhuoqiliu.com/HE0238-1904.html.

###  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|<\rm 2000\,km\,s^{-1}$ of the
quasar systemic velocity, which is sufficient to encompass most members of
even massive galaxy clusters. Figure 2 shows a $1.5^{\prime}\times
1.5^{\prime}$ 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 color of the circle represents the LOS velocity of each
galaxy relative to the quasar. Additionally, we display the $1^{\prime}\times
1^{\prime}$ MUSE FoV, and a smaller $30^{\prime\prime}\times
30^{\prime\prime}$ 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 M_{\odot}})>10.5\ (>\\!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 ($\Delta v=v-v_{\rm QSO}$) in the bottom
right panel of Figure 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 Figure 2 with galaxies
with $\Delta v>0\rm\,km\,s^{-1}$ largely falling North East of the quasar and
those with $\Delta v<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 Figure 2. We also show the LOS
velocity distribution of the galaxies in each trapezoidal region by the
corresponding histograms in the inset panel in Figure 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
suggests 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\%$ 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\%$ 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 Figure 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
M_{\odot}}$ and $M_{\rm dyn}\sim 5.7\times 10^{11}\ {\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>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 M_{\odot}}$. The mass range of
$M_{\rm dyn}\approx 4\times 10^{13}-10^{14}\ {\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 Figure 2, we show the stellar mass
weighted group center as a white asterisk, and membership weighted
($\frac{P_{\rm blue/red}}{P_{\rm blue}+P_{\rm red}}$) centers 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
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_{*}/\mathrm{M}_{\odot})=11.5$
and $11.2$, respectively. Furthermore, the richer, redshifted group contains a
massive disc galaxy, G33, with $\log(M_{*}/{\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.
Nonetheless, 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.

Figure 4: Visualizations of the nebula discovered around HE 0238$-$1904\.
Panel (a): HST ACS+F814W image of the field. Galaxies are circled in black and
labelled with their IDs. Panel (b): map of the nebular LOS velocity relative
to the quasar systemic velocity. Galaxies are circled in black and colored
with their velocities. Panel (c): map of nebular photoionization shown as the
line ratio $\rm[O\,III]\lambda 5008/[O\,II]\lambda\lambda 3727+3729$.
Panel(d)-(f) and Panel (g)-(i): narrow-band $\rm[O\,II]$ and $\rm[O\,III]$
surface brightness maps extracted from the MUSE datacube over the velocity
intervals labelled in each panel. The inset panel in Panel (h) shows a zoomed,
unsmoothed map around G3 and G5 to emphasize the possible existence of a tidal
tail. These maps are overlaid with $\rm[O\,II]$ and $\rm[O\,III]$ surface
brightness contours at levels of $0.08$ and $0.3\times
10^{-17}\rm\,erg\,cm^{-2}\,s^{-1}\,arcsec^{-2}$. The contours shown in panel
(e) and panel (h) are overlaid on the HST image in blue and red respectively.
We note that surface brightness maps and contours are smoothed with $3$ pixel
kernels. A version of this figure with the region circles marked in every
velocity panel is available online1. Table 2: Summary of emission-line
measurements for extracted regions in the nebula around HE 0238$-$1904.

ID | Distancea | Extraction | $\rm[O\,II]$ | $\rm H\beta$ | $\rm[O\,III]$ | $\rm[Ne\,V]$ | $\rm[O\,III]$ | $\rm He\,II$ |  $\Delta$vb |  $\sigma$c
---|---|---|---|---|---|---|---|---|---|---
| ($\rm pkpc$) | radius | $\lambda\lambda 3727+3729$ |  | $\lambda 5008$ | $\lambda 3346$ | $\lambda 4364$ | $\lambda 4687$ | ($\mathrm{km\ s^{-1}}$) | ($\mathrm{km\ s^{-1}})$
|  | (′′) | ($\mathrm{10^{-17}\,erg}$ | ($\mathrm{10^{-17}\,erg}$ | ($\mathrm{10^{-17}\,erg}$ | ($\mathrm{10^{-17}\,erg}$ | ($\mathrm{10^{-17}\,erg}$ | ($\mathrm{10^{-17}\,erg}$ |  |
|  |  |  $\mathrm{s^{-1}cm^{-2}}$) |  $\mathrm{s^{-1}cm^{-2}}$) |  $\mathrm{s^{-1}cm^{-2}}$) |  $\mathrm{s^{-1}cm^{-2}}$) |  $\mathrm{s^{-1}cm^{-2}}$) |  $\mathrm{s^{-1}cm^{-2}}$) |  |
S1 | 45 | 0.7 | $1.73\pm 0.05$ | $0.69\pm 0.06$ | $9.17\pm 0.05$ | $0.15\pm 0.03$ | $0.21\pm 0.02$ | $<0.21$ | $-11\pm 3$ | $62\pm 4$
S2 | 36 | 0.7 | $3.55\pm 0.08$ | $1.14\pm 0.14$ | $23.48\pm 0.10$ | $0.37\pm 0.05$ | $0.40\pm 0.04$ | $0.35\pm 0.11$ | $-55\pm 3$ | $43\pm 4$
S3 | 25 | 0.7 | $<0.30$ | $<0.27$ | $6.27\pm 0.22$ | $<0.15$ | $<0.09$ | $<0.18$ | $-107\pm 3$ | $61\pm 4$
$\rm S3_{wing}$ | 25 | 0.7 | $2.90\pm 0.10$ | $0.73\pm 0.09$ | $2.44\pm 0.22$ | $<0.18$ | $<0.12$ | $<0.21$ | $-14\pm 9$ | $104\pm 5$
S4 | 17 | 0.7 | $1.34\pm 0.18$ | $0.28\pm 0.08$ | $3.39\pm 0.10$ | $<0.18$ | $<0.09$ | $<0.15$ | $-114\pm 3$ | $45\pm 4$
$\rm S4_{wing}$ | 17 | 0.7 | $4.17\pm 0.20$ | $0.52\pm 0.09$ | $3.14\pm 0.12$ | $<0.27$ | $<0.15$ | $<0.27$ | $\hskip 0.85358pt+12\pm 8$ | $169\pm 6$
S5 | 9 | 0.7 | $5.96\pm 0.28$ | $0.77\pm 0.26$ | $2.51\pm 0.22$ | $<0.84$ | $<0.51$ | $<0.78$ | $\hskip 8.82036pt+8\pm 11$ | $140\pm 11$
S6 | 20 | 0.7 | $5.04\pm 0.07$ | $1.47\pm 0.12$ | $14.03\pm 0.07$ | $0.15\pm 0.05$ | $0.22\pm 0.04$ | $0.34\pm 0.09$ | $-62\pm 3$ | $68\pm 4$
S7 | 29 | 0.7 | $0.99\pm 0.04$ | $0.18\pm 0.06$ | $0.63\pm 0.04$ | $<0.09$ | $<0.06$ | $<0.18$ | $-72\pm 8$ | $111\pm 8$
S8 | 18 | 0.7 | $2.33\pm 0.04$ | $0.52\pm 0.06$ | $1.98\pm 0.04$ | $<0.09$ | $<0.06$ | $<0.15$ | $-119\pm 4$ | $89\pm 4$
S9 | 11 | 0.7 | $3.71\pm 0.16$ | $1.10\pm 0.15$ | $2.56\pm 0.13$ | $<0.45$ | $<0.27$ | $<0.39$ | $+173\pm 7$ | $110\pm 7$
S10 | 15 | 0.7 | $1.96\pm 0.05$ | $0.47\pm 0.05$ | $1.58\pm 0.04$ | $<0.12$ | $<0.09$ | $<0.15$ | $\hskip 0.85358pt+58\pm 4$ | $79\pm 5$
B1 | 49 | 1.4 | $1.14\pm 0.08$ | $0.89\pm 0.12$ | $2.21\pm 0.08$ | $<0.21$ | $<0.15$ | $<0.33$ | $\hskip 0.85358pt+50\pm 6$ | $128\pm 7$
B2 | 47 | 1.8 | $4.32\pm 0.13$ | $<0.57$ | $1.96\pm 0.11$ | $<0.30$ | $<0.21$ | $<0.57$ | $-36\pm 8$ | $119\pm 8$
B3 | 76 | 1.7 | $1.03\pm 0.09$ | $<0.60$ | $1.37\pm 0.07$ | $<0.21$ | $<0.15$ | $<0.42$ | $-69\pm 5$ | $50\pm 6$
B4 | 79 | 1.4 | $0.31\pm 0.11$ | $<0.24$ | $0.83\pm 0.06$ | $<0.18$ | $<0.12$ | $<0.33$ | $\hskip 0.85358pt+30\pm 4$ | $30\pm 6$
$\rm B4_{wing}$ | 79 | 1.4 | $0.99\pm 0.16$ | $<0.24$ | $0.40\pm 0.11$ | $<0.18$ | $<0.12$ | $<0.33$ | $\hskip 3.69885pt-83\pm 42$ | $201\pm 36$

  * a

Projected physical distance from the quasar.

  * b

LOS velocity relative to the quasar with uncertainty, assuming a systematic
uncertainty of $3\rm\,km\,s^{-1}$ (Weilbacher et al., 2020).

  * c

LOS velocity dispersion with uncertainty, assuming a systematic uncertainty of
$4\rm\,km\,s^{-1}$ (Kamann et al., 2016).

###  3.4 Nebular Environment

Due to ionizing radiation from the accretion disk, wide-field IFS observations
of quasar fields often find large nebulae (Johnson et al., in prep). To search
for the nebulae around HE 0238$-$1904, we conducted continuum subtraction of
the datacube locally for the $\rm[O\,II]$, $\rm H\beta$, and $\rm[O\,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\,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-fit model to
complete the continuum subtraction.

This continuum-subtracted MUSE datacube enabled the discovery of a giant
ionized nebula in $\rm[O\,II]$, $\rm H\beta$, and $\rm[O\,III]$ around HE
0238$-$1904 with a total area of $\approx 5000\ {\rm kpc}^{2}$ which is
visualized in Figure 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\,II]$ and $\rm[O\,III]$ emission contours over the
HST image in panel (a) of Figure 4. We also display a nebular LOS velocity map
in panel (b) and a $\rm[O\,III]/[O\,II]$ line ratio map in panel (c). We
constructed these two maps by jointly fitting Gaussian line profiles to the
continuum-subtracted $\rm[O\,II]$, $\rm H\beta$, and $\rm[O\,III]$ datacubes.
Instead of fitting the spectrum of each individual spaxel, we averaged over
circular apertures of $r=1^{\prime\prime}$ 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 narrowband images generated from the continuum subtracted
datacubes for each of $\rm[O\,II]$ and $\rm[O\,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 a 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 line width 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[O\,III]/[O\,II]$ line ratio. The majority of the nebula is $\rm[O\,II]$
dominated but the region South East of the quasar has greater $\rm[O\,III]$
emission, particularly, at a few $\rm[O\,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 ($\rm[O\,II]$, $\rm H\beta$, and
$\rm[O\,III]$) and weak lines ($\rm[Ne\,V]\lambda 3427$, $\rm H\delta$, $\rm
H\gamma$, $\rm[O\,III]\lambda 4364$, and $\rm He\,II\lambda 4687$222Other weak
lines such as $\rm[Ne\,III]\lambda 3869$, $\rm He\,I\lambda 3889$ & $\rm H8$,
and $\rm H\epsilon$ are covered by MUSE but we do not use them in this work
because of contaminating sky lines or blending with other lines.). 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 Figure 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^{\prime\prime}$ while B# labels low surface
brightness regions which required a larger extraction radius
($>1^{\prime\prime}$) 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., $\rm[O\,III]\lambda 4960$
and $\rm[O\,III]\lambda 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 1-D velocity dispersions, in Table 2. We display strong and
weak line spectra as well as their best-fit models in Figure 5 and Figure 6
respectively for a representative subset of the regions.

Figure 5: Examples of nebular spectra (stronger lines) and best-fit spectral
models for multiple regions. The locations of these regions are shown as
circles and labelled by their IDs in Figure 4. The extracted spectrum is shown
as solid black lines and the error array is shown as grey lines. The best-fit
models are shown as red solid lines. Figure 6: Examples of nebular spectra
(fainter lines) and best-fit spectral models for multiple regions. The
locations of these regions are shown as circles and labelled by their IDs in
Figure 4. The plotting style is as described in Figure 5.

##  4 Discussion

Figure 7: Emission line surface brightness profile for the nebula around HE
0238$-$1904\. The [O II] and [O III] profiles are extracted over a velocity
interval of $-600$ to $600\,\rm km\,s^{-1}$, and are circularly averaged at
different distances from the quasar centroid.

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 a 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 timescales of decades from the SDSS DR14 quasar catalog (Pâris et al.,
2018) and the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010)
R90 quasar catalog (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 timescale of the nebula is approximately $1\%$. 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 $\approx 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_{\rm e}/{\rm
cm}^{-3})\approx-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<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 Figure 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, 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. Second, 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. The first scenario requires a
coincident alignment between the nebula stripped from the dwarfs and the
lightcone of the quasar. Alternatively, the second scenario would naturally
result in the illumination of the nebula by the quasar since the nebula
surrounds it.

To investigate the location of the nebula relative to the quasar and the
emission geometry of HE 0238$-$1904, we show the surface brightness profiles
of [O II] and [O III] made with Photutils (Bradley, 2023) in Figure 7. The
profile of [O II] declines smoothly as a function of radius, while the [O III]
exhibits shallower drop due to the bright knots seen in the narrow-band
images. The centroids of narrow-band [O II] and [O III] surface brightness
maps are 10 and 19 pkpc away from the quasar, respectively. This
correspondence in position can be explained if the nebula surrounds the
quasar. Alternatively, the nebula could be more distant from the quasar and
reside within its ionization cone by chance. This is more probable if the
opening angle is wide, as suggested by Trainor & Steidel (2013); Borisova et
al. (2016a); Schmidt et al. (2018); den Brok et al. (2020). However, a large
opening angle also makes chance alignment of the [O II] centroid and the
quasar centroid within $15\%$ of the size of the nebula less likely.

The giant nebulae around HE 0238$-$1904 was 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 dataset 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$\beta$ 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.

Table 3: Summary of nebula regions in the Field of HE 0238$-$1904.

ID |  $\log(n_{\rm e,[O\,II]}/\mathrm{cm}^{-3})$a |  $\log(n_{\rm H,Cloudy}/\mathrm{cm}^{-3})$b |  ${\rm log}(U_{\rm Cloudy})$c
---|---|---|---
S1 | $<1.6$ | $1.6^{\hskip 0.56905pt+\hskip 0.28453pt0.1}_{-0.1}$ | $-2.2^{-0.1}_{\hskip 0.56905pt+\hskip 0.42677pt0.1}$
S2 | $<1.7$ | $1.7^{\hskip 0.56905pt+\hskip 0.28453pt0.1}_{-0.1}$ | $-2.1^{-0.1}_{\hskip 0.56905pt+\hskip 0.42677pt0.1}$
S3 | … | … | …
S4 | … | … | …
S5 | $<1.6$ | $4.2^{\hskip 0.56905pt+\hskip 0.28453pt0.2}_{-0.3}$ | $-3.0^{-0.2}_{\hskip 0.56905pt+\hskip 0.42677pt0.3}$
S6 | $\hskip 20.77051pt1.8^{\hskip 0.56905pt+\hskip 0.28453pt0.1}_{-0.1}$ | $2.7^{\hskip 0.56905pt+\hskip 0.42677pt0.1}_{-0.1}$ | $-2.5^{-0.1}_{\hskip 0.56905pt+\hskip 0.42677pt0.1}$
S7 | $<1.9$ | $3.0^{\hskip 0.56905pt+\hskip 0.28453pt0.3}_{-0.3}$ | $-3.2^{-0.3}_{\hskip 0.56905pt+\hskip 0.42677pt0.3}$
S8 | $<1.3$ | $3.5^{\hskip 0.56905pt+\hskip 0.28453pt0.1}_{-0.2}$ | $-3.3^{-0.1}_{\hskip 0.56905pt+\hskip 0.42677pt0.2}$
S9 | $<2.3$ | $4.1^{\hskip 0.56905pt+\hskip 0.28453pt0.2}_{-0.3}$ | $-3.5^{-0.2}_{\hskip 0.56905pt+\hskip 0.42677pt0.3}$
S10 | $<1.4$ | $3.6^{\hskip 0.56905pt+\hskip 0.28453pt0.2}_{-0.2}$ | $-3.3^{-0.2}_{\hskip 0.56905pt+\hskip 0.42677pt0.2}$
B1 | $<2.8$ | $2.1^{\hskip 0.56905pt+\hskip 0.28453pt0.1}_{-0.2}$ | $-2.7^{-0.1}_{\hskip 0.56905pt+\hskip 0.42677pt0.2}$
B2 | $<1.2$ | $2.9^{\hskip 0.56905pt+\hskip 0.28453pt0.1}_{-0.3}$ | $-3.4^{-0.1}_{\hskip 0.56905pt+\hskip 0.42677pt0.3}$
B3 | $<2.5$ | $1.9^{\hskip 0.56905pt+\hskip 0.28453pt0.1}_{-0.2}$ | $-2.8^{-0.1}_{\hskip 0.56905pt+\hskip 0.42677pt0.2}$
B4 | … | … | …

  * •

Notes.

  * a

Number density measurement or $95$% upper limit measured from
$\rm[O\,II]\lambda 3729/[O\,II]\lambda 3727$.

  * b

Number density inferred from Cloudy simulation described in the text.

  * c

Best-fit ionization parameter computed by Cloudy simulation.

###  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; Rupke et al., 2019; Chen et al., 2019; Helton et al.,
2021). The presence of the quasar suggests the source of ionization is AGN-
related. To study the physical conditions of the gas, we measured the the
density- and temperature-sensitive $\rm[O\,II]\lambda 3729/[O\,II]\lambda
3727$ and $\rm[O\,III]\lambda 4364/[O\,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 Figure 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\,II]\lambda 3727$, $\rm[O\,II]\lambda 3729$,
$\rm[O\,III]\lambda 4364$, and $\rm[O\,III]\lambda 5008$, we can directly
measure electron density ($n_{\rm e}$) and temperature ($T_{\rm e}$), as
discussed in Osterbrock & Ferland (2006). The $\rm[O\,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[O\,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\,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 et
al., 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\,II]$ doublet is resolved, the line ratio is
in the low density limit except S6. We therefore report 95$\%$ upper limits in
density for all but S6. The inferred electron number density upper limits
range from $1.2<\log(n_{\rm e,[O\,II]}/\mathrm{cm^{-3}})<2.8$, with a median
of $\log(n_{\rm e,[O\,II]}/\mathrm{cm^{-3}})<1.6$. These density upper limits
are consistent with gas arising from ionized ISM (Draine, 2011) or CGM. We
detected $\rm[O\,III]\lambda 4364$ in only three luminous regions, S1, S2, and
S6. The inferred temperatures for S1, S2, and S6 are
$\log(T/\mathrm{K})\approx 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_{\rm H}/\mathrm{cm^{-3}})\approx 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\,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\propto\nu^{\alpha}$) between 0.37 and 73.5 $\rm Ryd$, with $\alpha$
between $-1.8<\alpha<0$ following Groves et al. (2004) but extending to a
higher $\alpha$. We set the modeled 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_{\rm H}/\mathrm{cm}^{-3})<4.6$ and metallicity of
$-1.5<\log(Z/Z_{\odot})<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
$\alpha$ 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_{\odot})>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 $\alpha$. 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 Figure 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-fit $\alpha$ values for most regions are within
$-1.0<\alpha<-0.6$, somewhat greater than ones given in Groves et al. (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_{\odot})<0.2$. The densities
inferred from these photoionization simulations range from $\log(n_{\rm
H,Cloudy}/\mathrm{cm}^{-3})=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<\log(n_{\rm H,\,Cloudy}/\mathrm{cm}^{-3})<4.2$. However, the more direct
$\rm[O\,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, Figure 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\,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 mis-match, we modeled the emitting gas as a multi-
phase 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\,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 log-normal density
distribution (see Cantalupo et al., 2019). A log-normal distribution is
defined as

${\rm PDF}(n)\text{d}n=\frac{1}{\sqrt{2\pi}\sigma}{\rm
exp}\Big{[}-\frac{[\text{ln}(n)-\text{ln}(\mu)]^{2}}{2\sigma^{2}}\Big{]}\text{dln}(n)$
(1)

where $\sigma$ is the dispersion and $\mu$ 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 log-normal density distribution by interpolating Cloudy
models and integrating over the PDF. Our results show that a log-normal
distribution with a large $\sigma$ can reproduce the ionization-sensitive line
ratios, but the log-normal 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})>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 ionized by the quasar due to obscuring dust blocking its light from
reaching this gas. Filtering the quasar’s radiation through dust would soften
the incident ionizing radiation field. However, the best-fit $\alpha$ values
from our photoionization analysis suggests a hard ionizing spectrum for almost
all regions. Alternatively, the ionization of the nebulae could be due to
young stellar populations (Morisset et al., 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 Figure 8. Producing the observed [O III]/[O II] and
[Ne V]/[O II]333We note that [Ne V]/[Ne III] as a better shock tracer cannot
be used due to [Ne III]$\lambda$3869 is severely contaminated by skylines.
ratios requires shock velocities of $v_{\rm shock}>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 $>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 favor 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.

Figure 8: The emission line diagnostic diagrams log([O III]$\lambda 5008$/[O
II]$\lambda\lambda 3727,3729$) versus log([O III]$\lambda 5008$/H$\beta$) and
log([Ne V]$\lambda 3427$/[O II]$\lambda\lambda 3727,3729$) versus log([O
III]$\lambda 5008$/[O II]$\lambda\lambda 3727,3729$) for nebular regions. Line
ratios are shown as orange points with error bars, or shown as $3\sigma$ upper
limits for non-detections. For S3, S4, and B4, total line ratios (main+wing)
are shown as orange diamonds with error bars, or with upper limits. We note
that S3, S4, and B4 might have large uncertainty due to multiple components
detected within $150\rm\,km\,s^{-1}$. We compare these line ratios with the
fast radiative shock models (shock plus precursor) from Allen et al. (2008).
Emission-line ratio grids with solar metallicity, a preshock density of
$n_{\rm e}=100\,\rm cm^{-3}$, and a magnetic field strength of
$B=0.001{-}100\rm\,\mu G$ are shown in orange and grey for a shock velocity of
$100{-}250\rm\,km\,s^{-1}$ and $250{-}1000\rm\,km\,s^{-1}$ respectively.

Perhaps more likely, the difference in the density estimates could be due to
quasar variability (Richstone & Oke, 1977). Quasar variability is directly
observed on timescales 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^{4}$ to $10^{7}$ years
and AGN episodes may be highly clustered (e.g., Schirber et al., 2004;
Gonçalves et al., 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\,II]$ ratio by re-
running Cloudy models and MCMC analysis after significantly decreasing the
quasar luminosity.

Figure 9: Examples of violin plots and density posteriors for nebular regions.
For each subfigure, the left panel shows the density- and temperature-
sensitive line ratios from the flux measurements as red points with error bars
and the photoionization model posteriors are shown in orange. The middle panel
shows the ionization-sensitive line ratios from the flux measurements as red
points with error bars and the photoionization model posteriors in orange. The
right panel shows the density posterior from both direct (red histogram) and
indirect density estimates (orange filled histogram). The density posterior
inferred from the [O II] doublet extends below the plotted range as indicated
by the red arrows.

##  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.

  1. 1.

We found that HE 0238$-$1904 resides in an overdense environment containing
two potentially merging galaxy groups based 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 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 M_{\odot}}$ (Shen et al., 2009).

  2. 2.

We identified a giant nebula covering a projected area of $\approx 5000\ {\rm
kpc}^{2}$ around HE 0238$-$1904 emitting strongly in $\rm[O\,II]$,
$\mathrm{H}\beta$, and $\rm[O\,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.

  3. 3.

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.

  4. 4.

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\,II]$ doublet range from $\log(n_{\rm
e,[O\,II]}/\mathrm{cm^{-3}})<1.2$ to $2.8$, with a median of $\log(n_{\rm
e,[O\,II]}/\mathrm{cm^{-3}})<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\,II]$-based density upper
limits, reaching values of up to $400$ times higher.

  5. 5.

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 timescales of $10^{4}$ to $10^{5}$ years. 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 timescales of $10^{4}$ to
$10^{5}$ years has implications far beyond accretion disk physics in the
central engine. In particular, significant fluctuations on these timescales
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-meter 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 the ESO and HST data
archives.

## References

  * Abbott et al. (2021) Abbott T. M. C., et al., 2021, ApJS, 255, 20
  * Allen et al. (2008) Allen M. G., Groves B. A., Dopita M. A., Sutherland R. S., Kewley L. J., 2008, ApJS, 178, 20
  * Arav et al. (2013) Arav N., Borguet B., Chamberlain C., Edmonds D., Danforth C., 2013, MNRAS, 436, 3286
  * Assef et al. (2018) Assef R. J., Stern D., Noirot G., Jun H. D., Cutri R. M., Eisenhardt P. R. M., 2018, ApJS, 234, 23
  * Astropy Collaboration et al. (2022) Astropy Collaboration et al., 2022, ApJ, 935, 167
  * Bacon et al. (2010) Bacon R., et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. p. 773508 (arXiv:2211.16795), doi:10.1117/12.856027
  * Bacon et al. (2016) Bacon R., Piqueras L., Conseil S., Richard J., Shepherd M., 2016, MPDAF: MUSE Python Data Analysis Framework, Astrophysics Source Code Library, record ascl:1611.003 (ascl:1611.003)
  * Bertin & Arnouts (1996) Bertin E., Arnouts S., 1996, A&AS, 117, 393
  * Bish et al. (2021) Bish H. V., Werk J. K., Peek J., Zheng Y., Putman M., 2021, ApJ, 912, 8
  * Blanton & Roweis (2007) Blanton M. R., Roweis S., 2007, AJ, 133, 734
  * Bolton et al. (2012) Bolton A. S., et al., 2012, AJ, 144, 144
  * Borisova et al. (2016a) Borisova E., Lilly S. J., Cantalupo S., Prochaska J. X., Rakic O., Worseck G., 2016a, ApJ, 830, 120
  * Borisova et al. (2016b) Borisova E., et al., 2016b, ApJ, 831, 39
  * Boselli et al. (2019) Boselli A., et al., 2019, A&A, 631, A114
  * Bradley (2023) Bradley L., 2023, astropy/photutils: 1.8.0, doi:10.5281/zenodo.7946442, https://doi.org/10.5281/zenodo.7946442
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Buchner et al. (2014) Buchner J., et al., 2014, A&A, 564, A125
  * Burchett et al. (2021) Burchett J. N., Rubin K. H. R., Prochaska J. X., Coil A. L., Vaught R. R., Hennawi J. F., 2021, ApJ, 909, 151
  * Cai et al. (2019) Cai Z., et al., 2019, ApJS, 245, 23
  * Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, ApJ, 533, 682
  * Cantalupo et al. (2014) Cantalupo S., Arrigoni-Battaia F., Prochaska J. X., Hennawi J. F., Madau P., 2014, Nature, 506, 63
  * Cantalupo et al. (2019) Cantalupo S., et al., 2019, MNRAS, 483, 5188
  * Carnall et al. (2018) Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379
  * Carnall et al. (2019) Carnall A. C., et al., 2019, MNRAS, 490, 417
  * Chen et al. (2019) Chen H.-W., Boettcher E., Johnson S. D., Zahedy F. S., Rudie G. C., Cooksey K. L., Rauch M., Mulchaey J. S., 2019, ApJ, 878, L33
  * Chen et al. (2020) Chen H.-W., et al., 2020, MNRAS, 497, 498
  * Chen et al. (2023) Chen M. C., et al., 2023, MNRAS, 518, 2354
  * Coleman et al. (1980) Coleman G. D., Wu C. C., Weedman D. W., 1980, ApJS, 43, 393
  * Comerford et al. (2017) Comerford J. M., Barrows R. S., Müller-Sánchez F., Nevin R., Greene J. E., Pooley D., Stern D., Harrison F. A., 2017, ApJ, 849, 102
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium
  * Dutta et al. (2023a) Dutta R., et al., 2023a, arXiv e-prints, p. arXiv:2302.09087
  * Dutta et al. (2023b) Dutta S., Muzahid S., Schaye J., Mishra S., Chen H.-W., Johnson S., Wisotzki L., Cantalupo S., 2023b, arXiv e-prints, p. arXiv:2303.16933
  * Epinat et al. (2018) Epinat B., et al., 2018, A&A, 609, A40
  * Faber et al. (2007) Faber S. M., et al., 2007, ApJ, 665, 265
  * Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  * Ferland et al. (2017) Ferland G. J., et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385
  * Feroz et al. (2009) Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601
  * Feroz et al. (2019) Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2019, The Open Journal of Astrophysics, 2, 10
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fossati et al. (2019) Fossati M., et al., 2019, MNRAS, 490, 1451
  * Fossati et al. (2021) Fossati M., et al., 2021, MNRAS, 503, 3044
  * Furesz et al. (2020) Furesz G., et al., 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. p. 114470A, doi:10.1117/12.2562803
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Gallazzi et al. (2005) Gallazzi A., Charlot S., Brinchmann J., White S. D. M., Tremonti C. A., 2005, MNRAS, 362, 41
  * Gonçalves et al. (2008) Gonçalves T. S., Steidel C. C., Pettini M., 2008, ApJ, 676, 816
  * Green et al. (2012) Green J. C., et al., 2012, ApJ, 744, 60
  * Groves et al. (2004) Groves B. A., Dopita M. A., Sutherland R. S., 2004, ApJS, 153, 9
  * Guo et al. (2019) Guo H., Liu X., Shen Y., Loeb A., Monroe T., Prochaska J. X., 2019, MNRAS, 482, 3288
  * Heckman et al. (1981) Heckman T. M., Miley G. K., van Breugel W. J. M., Butcher H. R., 1981, ApJ, 247, 403
  * Helton et al. (2021) Helton J. M., Johnson S. D., Greene J. E., Chen H.-W., 2021, MNRAS, 505, 5497
  * Hess et al. (2017) Hess K. M., Cluver M. E., Yahya S., Leisman L., Serra P., Lucero D. M., Passmoor S. S., Carignan C., 2017, MNRAS, 464, 957
  * Hester (2006) Hester J. A., 2006, ApJ, 647, 910
  * Hewett & Wild (2010) Hewett P. C., Wild V., 2010, MNRAS, 405, 2302
  * Johnson et al. (2015) Johnson S. D., Chen H.-W., Mulchaey J. S., 2015, MNRAS, 449, 3263
  * Johnson et al. (2018) Johnson S. D., et al., 2018, ApJ, 869, L1
  * Johnson et al. (2022) Johnson S. D., et al., 2022, ApJ, 940, L40
  * Kamann et al. (2016) Kamann S., et al., 2016, A&A, 588, A149
  * Kellermann et al. (1989) Kellermann K. I., Sramek R., Schmidt M., Shaffer D. B., Green R., 1989, AJ, 98, 1195
  * Kennicutt & Evans (2012) Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531
  * Kirkman & Tytler (2008) Kirkman D., Tytler D., 2008, MNRAS, 391, 1457
  * Konidaris et al. (2020) Konidaris N. P., Rudie G. C., Newman A. B., Hare T. S., Williams J. E., Lanz A. E., Kelson D. D., Crane J., 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. p. 114471E, doi:10.1117/12.2561171
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Langen et al. (2023) Langen V., Cantalupo S., Steidel C. C., Chen Y., Pezzulli G., Gallego S. G., 2023, MNRAS, 519, 5099
  * Leclercq et al. (2022) Leclercq F., et al., 2022, A&A, 663, A11
  * Lehner et al. (2018) Lehner N., Wotta C. B., Howk J. C., O’Meara J. M., Oppenheimer B. D., Cooksey K. L., 2018, ApJ, 866, 33
  * Leitner & Kravtsov (2011) Leitner S. N., Kravtsov A. V., 2011, ApJ, 734, 48
  * Liu et al. (2013a) Liu G., Zakamska N. L., Greene J. E., Nesvadba N. P. H., Liu X., 2013a, MNRAS, 430, 2327
  * Liu et al. (2013b) Liu G., Zakamska N. L., Greene J. E., Nesvadba N. P. H., Liu X., 2013b, MNRAS, 436, 2576
  * Luridiana et al. (2015) Luridiana V., Morisset C., Shaw R. A., 2015, A&A, 573, A42
  * Mackenzie et al. (2021) Mackenzie R., et al., 2021, MNRAS, 502, 494
  * Marasco et al. (2016) Marasco A., Crain R. A., Schaye J., Bahé Y. M., van der Hulst T., Theuns T., Bower R. G., 2016, MNRAS, 461, 2630
  * Martin et al. (2010) Martin C., Moore A., Morrissey P., Matuszewski M., Rahman S., Adkins S., Epps H., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. p. 77350M, doi:10.1117/12.858227
  * Mateo et al. (2022) Mateo M., Bailey J. I., Song Y., Crane J., Hull C., Shectman S., Birk C., 2022, in Evans C. J., Bryant J. J., Motohara K., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 12184, Ground-based and Airborne Instrumentation for Astronomy IX. p. 121845P, doi:10.1117/12.2629506
  * Monroe et al. (2016) Monroe T. R., Prochaska J. X., Tejos N., Worseck G., Hennawi J. F., Schmidt T., Tumlinson J., Shen Y., 2016, AJ, 152, 25
  * Morisset et al. (2015) Morisset C., Delgado-Inglada G., Flores-Fajardo N., 2015, Rev. Mex. Astron. Astrofis., 51, 103
  * Munari et al. (2013) Munari E., Biviano A., Borgani S., Murante G., Fabjan D., 2013, MNRAS, 430, 2638
  * Muzahid et al. (2012) Muzahid S., Srianand R., Savage B. D., Narayanan A., Mohan V., Dewangan G. C., 2012, MNRAS, 424, L59
  * Muzahid et al. (2018) Muzahid S., Fonseca G., Roberts A., Rosenwasser B., Richter P., Narayanan A., Churchill C., Charlton J., 2018, MNRAS, 476, 4965
  * Muzahid et al. (2020) Muzahid S., et al., 2020, MNRAS, 496, 1013
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * National Academies of Sciences (2021) National Academies of Sciences 2021, Pathways to Discovery in Astronomy and Astrophysics for the 2020s, doi:10.17226/26141.
  * Newville et al. (2014) Newville M., Stensitzki T., Allen D. B., Ingargiola A., 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, Zenodo, doi:10.5281/zenodo.11813
  * Nyland et al. (2020) Nyland K., et al., 2020, ApJ, 905, 74
  * O’Sullivan et al. (2020) O’Sullivan D. B., Martin C., Matuszewski M., Hoadley K., Hamden E., Neill J. D., Lin Z., Parihar P., 2020, ApJ, 894, 3
  * Oppenheimer & Schaye (2013) Oppenheimer B. D., Schaye J., 2013, MNRAS, 434, 1063
  * Oppenheimer et al. (2018) Oppenheimer B. D., Segers M., Schaye J., Richings A. J., Crain R. A., 2018, MNRAS, 474, 4740
  * Osterbrock & Ferland (2006) Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Pâris et al. (2018) Pâris I., et al., 2018, A&A, 613, A51
  * Poggianti et al. (2016) Poggianti B. M., et al., 2016, AJ, 151, 78
  * Ren et al. (2018) Ren B., Pueyo L., Zhu G. B., Debes J., Duchêne G., 2018, ApJ, 852, 104
  * Richard (2019) Richard J., 2019, in The Very Large Telescope in 2030. p. 24, doi:10.5281/zenodo.3356260
  * Richards et al. (2006) Richards G. T., et al., 2006, ApJS, 166, 470
  * Richstone & Oke (1977) Richstone D. O., Oke J. B., 1977, ApJ, 213, 8
  * Robitaille & Bressert (2012) Robitaille T., Bressert E., 2012, APLpy: Astronomical Plotting Library in Python, Astrophysics Source Code Library (ascl:1208.017)
  * Rupke et al. (2017) Rupke D. S. N., Gültekin K., Veilleux S., 2017, ApJ, 850, 40
  * Rupke et al. (2019) Rupke D. S. N., et al., 2019, Nature, 574, 643
  * Schawinski et al. (2015) Schawinski K., Koss M., Berney S., Sartori L. F., 2015, MNRAS, 451, 2517
  * Schirber et al. (2004) Schirber M., Miralda-Escudé J., McDonald P., 2004, ApJ, 610, 105
  * Schmidt et al. (2018) Schmidt T. M., Hennawi J. F., Worseck G., Davies F. B., Lukić Z., Oñorbe J., 2018, ApJ, 861, 122
  * Segers et al. (2017) Segers M. C., Oppenheimer B. D., Schaye J., Richings A. J., 2017, MNRAS, 471, 1026
  * Shen (2021) Shen Y., 2021, ApJ, 921, 70
  * Shen et al. (2009) Shen Y., et al., 2009, ApJ, 697, 1656
  * Soto et al. (2016) Soto K. T., Lilly S. J., Bacon R., Richard J., Conseil S., 2016, MNRAS, 458, 3210
  * Stone et al. (2022) Stone Z., et al., 2022, MNRAS, 514, 164
  * Syphers & Shull (2014) Syphers D., Shull J. M., 2014, ApJ, 784, 42
  * Tacconi et al. (2013) Tacconi L. J., et al., 2013, ApJ, 768, 74
  * Thatte et al. (2022) Thatte N. A., et al., 2022, in Evans C. J., Bryant J. J., Motohara K., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 12184, Ground-based and Airborne Instrumentation for Astronomy IX. p. 1218420, doi:10.1117/12.2628834
  * Trainor & Steidel (2013) Trainor R., Steidel C. C., 2013, ApJ, 775, L3
  * Urrutia et al. (2019) Urrutia T., et al., 2019, A&A, 624, A141
  * Véron-Cetty & Véron (2006) Véron-Cetty M. P., Véron P., 2006, A&A, 455, 773
  * Vestergaard & Peterson (2006) Vestergaard M., Peterson B. M., 2006, ApJ, 641, 689
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Weilbacher et al. (2012) Weilbacher P. M., Streicher O., Urrutia T., Jarno A., Pécontal-Rousset A., Bacon R., Böhm P., 2012, in Radziwill N. M., Chiozzi G., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 8451, Software and Cyberinfrastructure for Astronomy II. p. 84510B, doi:10.1117/12.925114
  * Weilbacher et al. (2014) Weilbacher P. M., Streicher O., Urrutia T., Pécontal-Rousset A., Jarno A., Bacon R., 2014, in Manset N., Forshay P., eds, Astronomical Society of the Pacific Conference Series Vol. 485, Astronomical Data Analysis Software and Systems XXIII. p. 451 (arXiv:1507.00034), doi:10.48550/arXiv.1507.00034
  * Weilbacher et al. (2020) Weilbacher P. M., et al., 2020, A&A, 641, A28
  * Wisotzki et al. (2000) Wisotzki L., Christlieb N., Bade N., Beckmann V., Köhler T., Vanelle C., Reimers D., 2000, A&A, 358, 77
  * Wright et al. (2010) Wright E. L., et al., 2010, AJ, 140, 1868
  * Zabl et al. (2021) Zabl J., et al., 2021, MNRAS, 507, 4294
  * Zahedy et al. (2019) Zahedy F. S., Chen H.-W., Johnson S. D., Pierce R. M., Rauch M., Huang Y.-H., Weiner B. J., Gauthier J.-R., 2019, MNRAS, 484, 2257
  * Zhang (2018) Zhang D., 2018, Galaxies, 6, 114
  * Zhao & Wang (2023) Zhao Q., Wang J., 2023, ApJ, 943, L25
  * Zheng et al. (2019) Zheng Y., Peek J. E. G., Putman M. E., Werk J. K., 2019, ApJ, 871, 35
  * den Brok et al. (2020) den Brok J. S., et al., 2020, MNRAS, 495, 1874