11institutetext: 1Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir,
km 4, 28850 Torrejón de Ardoz, Madrid, Spain
2Institute of Astronomy and Astrophysics, Academia Sinica, No. 1, Section 4,
Roosevelt Road, Taipei 10617, Taiwan
3Instituto de Astrofísica de Andalucía, CSIC, Apartado de correos 3004,
E-18080 Granada, Spain
4GRANTECAN, Cuesta de San José s/n, 38712 Breña Baja, La Palma, Spain
5Instituto de Astrofísica de Canarias, Vía Láctea s/n, 38200 La Laguna,
Tenerife, Spain
11email<EMAIL_ADDRESS>
# AGN feedback can produce metal enrichment on galaxy scales

M. Villar Martín1 C. López Cobá2 S. Cazzoli3 E. Pérez Montero3 A. Cabrera
Lavers4,5 Based on observations made at Paranal Observatory under programme ID
0102.B-0107

###### Abstract

Context. Giant ($>$100 kpc) nebulae associated with active galaxies provide
rich information about the circumgalactic medium (CGM) around galaxies, its
link with the interstellar medium (ISM) of the hosts and the mechanisms
involved in their evolution.

Aims. We have studied the giant nebula associated with the Teacup ($z$=0.085)
quasar based on VLT MUSE integral field spectroscopy with two main scientific
aims: to investigate whether the well known giant ($\sim$10 kpc) AGN induced
outflow has an impact on the distribution of heavy elements in and outside the
host galaxy, and to establish whether there is a physical and/or evolutionary
link between the CGM and the ISM.

Methods. We have mapped the O/H and N/O gas abundances in two spatial
dimensions across the giant nebula and within the galaxy by means of comparing
emission line ratios with photoionization model predictions. We have also
compared the kinematics of the giant nebula with that of the gas and the stars
in the host galaxy.

Results. The widely studied AGN driven outflow responsible for the $\sim$10
kpc ionized bubble is enhancing the gas metal abundance up to $\sim$10 kpc
from the AGN. O/H is solar or slightly higher in the bubble edges, in
comparison with the subsolar abundances across the rest of the nebula median
(O/H$\sim$0.63 (O/H)⊙). The giant nebula shows complex kinematics. In addition
to non circular motions, we propose that there is a rotational component that
extends up to the nebular edge. It seems to be kinematically coupled at least
partially with the stellar component and the ISM of the host galaxy.

Conclusions. AGN feedback can produce metal enrichment at large extranuclear
distances in galaxies ($\geq$10 kpc). The dynamics of the giant nebula is
consistent with CGM studies based on absorption lines that show evidence for a
substantial rotation component in the CGM of numerous galaxies which,
moreover, is often kinematically coupled with the ISM and stars of the
galaxies. This study shows the great potential of using giant nebulae to
investigate the CGM around active galaxies at all redshifts from its line
emission.

###### Key Words.:

galaxies: evolution – galaxies: active – galaxies: abundances – quasars:
individual: Teacup

##  1 Introduction

The circumgalactic medium (CGM) is the bound gas halo surrounding galaxies
outside their interstellar medium (ISM) and inside their virial radius,
extending out to a few hundred kpc. The CGM may be the key regulator of the
galactic gas supply. Gas flows occurring between the CGM and the interstellar
medium (ISM) are thought to shape galaxies and drive their evolution, via
feedback, accretion and recycling of gas. Thus, investigating the metallicity,
structure and kinematics of the different gas phases will help us to
understand how galaxies gain, eject and recycle the gas during their existence
(see Tumlinson et al. 2017 for a review).

Because the CGM is very diffuse and therefore almost invisible in emission,
its physical properties remain largely unconstrained. Its understanding has
come so far mostly from studies based on absorption lines produced by the CGM
around galaxies, imprinted on the spectrum of background objects such as
quasars and radio galaxies. The CGM of a few non-active nearby galaxies has
also been studied in emission (e.g. Hayes et al. 2016), but the procedure is
extremely challenging. The presence of a powerful active galactic nucleus
(AGN) can render the CGM gas observable in emission around galaxies up to many
10s of kpc, well into the CGM. Giant emission line nebulae (size $\ga$60 kpc
and sometimes $>$100 kpc) associated with quasars and radio galaxies at
different redshifts have been studied since the 80’s (e.g.Baum et al. 1988;
McCarthy et al. 1990; McCarthy 1993; van Ojik et al. 1996; Villar-Martín et
al. 2003; Borisova et al. 2016; Villar-Martín et al. 2018; Fossati et al.
2021; Balmaverde et al. 2022; Wang et al. 2023). The discovery of giant
($\geq$70 kpc) widely spread reservoirs of molecular gas (e.g. Emonts et al.
2016; Falkendal et al. 2021) associated with several high $z$ ($z$$\ga$2)
radio galaxies reveals a multiphase CGM that was chemically enriched when the
Universe was 3 Gyr old and supplies gas from which galaxies grow.

Wide field integral field spectroscopic instruments on 8-10 meter telescopes,
such as MUSE on VLT, have opened excellent opportunities to detect and study
in great detail the elusive material from the warm ($T<$105 K) ionized CGM
around powerful active galaxies at different redshifts thanks to the
illumination by the powerful AGN. If it were not for the excitation by the AGN
continuum, this gas might not be detected or would only be observable through
absorption line studies.

The object of this work is the well known Teacup radio quiet type 2 quasar at
$z=$0.085. It shows a $\sim$10 kpc loop of ionized gas resembling a handle of
a teacup (hereof its nickname), which was discovered by volunteers of the
Galaxy Zoo project (Keel et al. 2012) and has been widely studied in the
context of AGN feedback and its potential impact on galaxy evolution. The
system has been proposed to be the scenario of a giant outflow generated
either by an AGN wind or induced by a 1 kpc radio jet whose effects are
noticed up to at least $\sim$10 kpc from the AGN and might be responsible for
the bubble-like morphology (Gagne et al. 2014; Harrison et al. 2015; Keel et
al. 2015; Ramos Almeida et al. 2017; Villar-Martín et al. 2018; Moiseev &
Ikhsanova 2023; Venturi et al. (2023, hereafter V23 ).

Villar-Martín et al. (2018) discovered a $>$100 kpc ionized nebula associated
with this object (see also Villar Martín et al. 2021, which could be a product
of a merger that occurred 1-2 Gyr ago (Keel et al. 2015). This rich gas
reservoir, which extends into the CGM, has been rendered visible due to the
activity of the quasar nucleus. AGN photoionization dominates the excitation
of the spatially resolved gas emission up to its boundary (Gagne et al. 2014;
Villar-Martín et al. 2018; V23 ; Moiseev & Ikhsanova 2023), except at some
locations outside the putative quasar ionization cones where evidence for
shock excitation has been found (V23 ). Stellar photoionization could also
contribute to ionize the gas locally in some tidal features. The large scale
kinematics are strongly reminiscent of rotation (Villar-Martín et al. 2018;
Moiseev & Ikhsanova 2023) and tentative results suggest subsolar nebular
abundances ($\sim$0.5$Z_{\rm\odot}$, Villar-Martín et al. 2018). The well
known bubble appears to be expanding from the nucleus and out into the nebula.

We present in this paper a detailed optical spectroscopic study of the Teacup
nebula based on VLT-MUSE archival data with two main goals: mapping the gas
abundances in two spatial dimensions and compare the kinematics of the giant
nebula on the largest scales with that of the ISM and specially the stars in
the host galaxy. The ultimate science goals are to establish whether there is
a physical and/or evolutionary link between the CGM and the ISM and whether
the giant outflow has an impact on the distribution of heavy elements from the
nucleus on large spatial scales. We will discuss the results in the context of
other studies of the CGM and its role in the evolution of galaxies.

The paper is organized as follows: the data and analysis method are presented
in Sect. 2 and 3 respectively. The latter includes the description of the
methods used to derive the gas chemical abundances and physical properties and
to model the gas and stellar kinematics. The results are presented in Sect. 4
and discussed in Sect. 5. The main conclusions are summarised in Sect. 6.

Throughout this paper, we assume flat $\Lambda$CDM cosmology following Planck
Collaboration et al. (2020), with H0 = 67.4 km s-1 Mpc-1, $\Omega_{\rm m}$ =
0.31. This gives a spatial scale of 1.65 kpc arcsec-1 at $z=0.085$.

##  2 Data

The data were collected for the 0102.B-0107 program (principal investigator,
PI: L. Sartori; see V23  for details) with the European Southern Observatory
(ESO) Very Large Telescope (VLT) and the Multi Unit Spectroscopic Explorer
(MUSE, Bacon et al. 2010). This instrument covers a 1$\arcmin\times$1$\arcmin$
field of view (FoV) in the Wide Field Mode (WFM), with a spatial sampling of
0.2$\arcsec$ pix-1. The wavelength coverage is $\sim$4650-9300 Å with a 1.25 Å
pix-1 spectral sampling and a resolving power
$R=\lambda$/$\Delta\lambda\sim$1700-3400 ($\Delta V\sim$176-88 km s-1).

The observations were performed in March 2019. The processed archive data cube
was used for this study. The Teacup ionized nebula fills a large fraction of
the MUSE FoV. Separate sky cubes are not available from these observations.
Sky oversubtraction due to contamination of the sky spectrum by object
emission was identified at a few spatial locations of very low surface
brightness which appeared as artefacts that mimic absorption features adjacent
to the strongest emission lines (specially [OIII]$\lambda$5007). Because our
attempts to improve the sky subtraction did not achieve significant
improvements, we finally used the archive cube. To further evaluate the
potential impact of this effect, a comparison with the processed archive
datacube from program 0103.B-0071 (PI: C. Harrison) was also performed.
Although separate sky cubes were obtained for this program, the resulting sky
subtraction was not significantly better. The comparison was in any case
valuable to characterize the impact of the sky subtraction on the data
quality. The artefacts are in general faint in comparison with the emission
lines and are shifted in wavelength so that they could be efficiently
identified and isolated. We confirm that imperfections on the sky subtraction
do not affect our results and conclusions. These are all based on the
0102.B-0107 datacube because of the significantly better seeing
(FWHM=0.74$\arcsec$ vs. 1.3$\arcsec$ for the 0103.B-0071 program) and of
somewhat better signal-to-noise ratio (SN).

##  3 Analysis

###  3.1 Spatially resolved emission line flux measurements

One of our main goals is to map the electron temperature, $T_{\rm e}$, and the
oxygen abundance O/H in two spatial dimensions to investigate whether the
giant outflow has an impact on the distribution of heavy elements across the
galaxy and out into the CGM. For the $T_{\rm e}$ determination, it is
essential to measure the flux ratio of a nebular to an auroral emission line,
such as the [OIII]$\lambda$5007/$\lambda$4363 (or
[NII]$\lambda$6583/$\lambda$5755). For this, we defined spatial apertures at
different locations through the nebula and extracted the integrated spectrum
from each one.

The apertures were selected based on the visual inspection of the
[OIII]$\lambda$5007 morphology at different wavelengths (i.e. velocities)
scanned through the line profile. These scans reveal striking morphological
changes with velocity and are specially useful to identify some faint
structures, such as tidal tails, arcs, knots, filaments, which have similar
velocities (Fig. 2). By restricting the central $\lambda$ and the width,
$\Delta\lambda$, of the narrow band images, the contrast of the lowest surface
brightness features is enhanced and the spatial apertures covering them can be
defined with greater accuracy to maximize the SN of the integrated spectra.
This was useful to detect faint emission lines through the largest possible
extension, while preserving at the same time the morphological information on
the diversity of nebular structural elements.

Figure 1: [OIII] continuum subtracted images covering different spectral
windows that were selected to highlight the diversity of nebular morphological
features. Each image covers a different velocity (i.e. spectral) range
relative to the nuclear systemic velocity as indicated on top. The nebular
morphology strongly varies with velocity. The left panel shows the total
[OIII] flux narrow band image. The well known $\sim$10 kpc ionized bubble is
marked with a tiny yellow star in the left panel. The green lines indicate the
position angles of the radio axis to the NE and to the W (Harrison et al.
2015). To guide the reader, the letters ‘A’ to ‘F’ mark some emission line
features that can also be identified in the mask map of Fig. 2.

Figure 2: Map of the masks used in our analysis (left) and [OIII] kinematic
maps (middle and right panels). The colors in the first map have no particular
meaning but help differentiating the apertures. A 1D spectrum was extracted
from each one, so that a single $W_{\rm 80}$ and $V_{\rm s}$ (middle and left
panels) values are associated to each aperture. The maps cover the total MUSE
FoV ($\sim$1$\arcmin\times$1$\arcmin$). $W_{\rm 80}$ and $V_{\rm s}$ are in km
s-1. Letters ‘A’ to ‘F’ , yellow star and solid lines have the same meaning as
in Fig. 2.

The map of the positions of the resulting 64 selected apertures is shown if
Fig. 2 (left panel). We also show in adjacent panels of the same figure the
$W_{\rm 80}$ and $V_{\rm s}$ maps based on [OIII]$\lambda$5007\. $W_{\rm 80}$
is the velocity width that encloses 80% of the total line flux and $V_{\rm s}$
is the velocity shift relative to the narrow core of the nuclear [OIII] line.
These maps can be directly compared with those of V23 , that were based on a
spaxel by spaxel analysis. The huge dimensions of the nebula both along and
perpendicular to the radio axis are apparent. Although some prominent features
stand out (knots, filaments, etc), gas emission seems to fill the entire area
within the nebular outer boundaries. The large scale rotation pattern (Gagne
et al. 2014; Villar-Martín et al. 2018; Moiseev & Ikhsanova 2023) and the
enhanced width, $W_{\rm 80}$, of the lines in the direction perpendicular to
the radio axis are clear (V23 ).

The 1D spectra extracted from the individual apertures were used to measure
the fluxes of the lines involved in our analysis:
[OIII]$\lambda\lambda$4959,5007 and $\lambda$4363, H$\gamma$, H$\beta$,
H$\alpha$, [NII]$\lambda\lambda$6548,6583, and $\lambda$5755 and
[SII]$\lambda$6716,6731.

For a subset of spectra the stellar continuum was relatively strong compared
with the emission line fluxes, that is, the equivalent width of the lines was
low, so that it was first necessary to fit it and subtract it. The effects of
stellar absorption can affect the Balmer emission lines and need to be
accounted for to obtain accurate reddening estimates. It was also necessary in
several spectra to reconstruct the H$\gamma$ and [OIII]$\lambda$4363 baseline
and recover both line fluxes. The continuum was fitted for each spectrum with
the pypipe3d tool (e.g., Lacerda et al., 2022), that will be described in more
detail in 3.3. As a test, we checked the reddening inferred from
H$\gamma$/H$\beta$ and H$\alpha$/H$\beta$ in the final, corrected spectra.
Case B H$\alpha$/H$\beta$=2.87 and H$\gamma$/H$\beta$=0.47 were assumed.
Reddening estimations from both ratios are in general consistent within the
errors. Changing these values within the range H$\alpha$/H$\beta$=2.76-3.05
and H$\gamma$/H$\beta$=0.45-0.47 (Osterbrock 1989) has no significant impact.
The $E_{B-V}$ values implied by H$\alpha$/H$\beta$ were used for extinction
correction since they have smaller errors. Using instead the $E_{B-V}$ implied
by H$\gamma$/H$\beta$ has also a negligible impact on the results.

###  3.2 Derivation of temperatures and gas chemical abundances

We used the python code PyNeb (Luridiana et al. 2015) to estimate $T_{\rm e}$
of [OIII] for spectra with detected [OIII]$\lambda$4363 using
[OIII]$\lambda$5007/$\lambda$4363, and, when possible, of [NII] using
[NII]$\lambda$6583/$\lambda$5755\. This code was also used to infer the
electron density $n_{\rm e}$ using [SII]$\lambda\lambda$6717,6731. We
calculated the corresponding errors by applying to the calculations a Monte
Carlo iteration using the nominal fluxes perturbed with the observational
errors.

For the derivation of abundances, the direct method (i.e. the derivation of
total chemical abundances from their ionic fractions and the measured electron
temperature) cannot be satisfactorily applied for the case of the Narrow Line
Regions (NLR) in Active Galactic Nuclei (AGN) (e.g. Dors et al. 2015).
Instead, we used the code HII-CHI-mistry (hereafter HCm, Pérez-Montero 2014).
The approach of this code consists of a bayesian-like comparison between
certain reddening-corrected ratios of emission line fluxes relative to a
recombination H line (H$\beta$ in our case), that are sensitive to total
oxygen abundance (O/H), nitrogen-to-oxygen ratio (N/O), and ionization
parameter ($U$), with the predictions from a large grid of photoionization
models. In particular, we used version 5.3 of HCm, which considers models
calculated under the most usual conditions in the NLR of active galaxies. As a
result, the most probable values and their uncertainties for O/H, N/O and
$log(U)$ are obtained.

The grid of models is explained in detail in Pérez-Montero et al. (2019). The
gas is assumed to be distributed homogeneously with a filling factor of 0.1
and a constant electron density of $n_{\rm e}$=500 cm-3. The densities across
the Teacup nebula are however significantly lower ($n_{\rm e}\la$150 cm-3, see
Table 2 and V23 ). Since collisional de-excitation effects are not relevant
for the line ratios under consideration in this density regime, the inferred
abundances are not affected by the comparatively high density used in the
models (Pérez-Montero et al. 2019).

Figure 3: Stellar and H$\alpha$ velocity maps for the $28\arcsec\times
19\arcsec$ resized data cube (see text). Top panels from left to right: $gri$
continuum image reconstructed from the MUSE cube; stellar velocity map from
the SSP analysis; best kinematic model from $\mathtt{XS}$; residual of the
model (data minus model); rotational curve, with $V_{t}$ being the circular
velocity. Bottom panels: Similar to the stellar velocity but for the H$\alpha$
velocity map. The shaded region in the rotational curves represent $1\sigma$
error bars. Black contours in the first and second panels represent iso-
velocities spaced each $\pm 50$ km s-1. The yellow contours highlight the
V-band light distribution, with contours starting from 0.003 to 0.4$\times
10^{-16}$ erg s-1 cm-2 Å-1.

The Spectral Energy Distribution (SED) was considered to be composed by two
components: one representing the Big Blue Bump peaking at 1 Ryd, and the other
a power law with spectral index $\alpha_{\rm OX}$=-1.2 representing the non-
thermal X-rays radiation (see Pérez-Montero et al. 2019 for additional
details). The stopping criterion to measure the resulting emergent spectrum is
that the proportion of free electrons in the ionized gas is lower than 2%. The
models consider a Mathis et al. (1977) grain size distribution and a dust-to-
gas mass ratio of 7.5$\times$10-3 (Rémy-Ruyer et al. 2014).

The grid of models vary the oxygen and nitrogen abundances within
12+log(O/H)=6.9-9.1 in bins of 0.1 dex and log(N/O)=-2.0-0.0 in bins of 0.125
dex, respectively. The rest of chemical abundances are scaled in the models
with respect to oxygen following the solar proportions,
12+log(O/H)=8.69$\pm$0.04, given by (Asplund et al. 2021). In accordance with
these authors, we also assumed 12+log(N/H)=7.83$\pm$0.07 for the Sun.

The original grid of models in Pérez-Montero et al. (2019) covers values of
$log(U)$ from -4.0 to -0.5 in bins of 0.25 dex, although we only assumed
$log(U)>$-2.5 to break the degeneracy with $U$ for the ratios of high- to low-
excitation lines in the optical range predicted by the models. This
restriction does not imply any significant change in the derivation of the
chemical abundances (see Pérez-Montero et al. 2019 for a better
clarification). The code also provides uncertainties to the derived quantities
calculated as the quadratical addition of the standard deviation of the
obtained bayesian distribution with the uncertainty obtained from a Monte
Carlo iteration through the given input nominal flux of each line perturbed
with its corresponding observational error.

###  3.3 Analysis and modelling of the stellar and gas kinematics

The kinematics of the Teacup ionized gas has been studied in numerous works
(Gagne et al. 2014; Harrison et al. 2015; Ramos Almeida et al. 2017; Villar-
Martín et al. 2018; Moiseev & Ikhsanova 2023; V23 ). Our main purpose here is
to compare the latter with the stellar kinematics in the host galaxy. To
achieve this, the analysis of the cube was performed using the pypipe3d tool
(e.g., Lacerda et al., 2022); this pipeline has been extensively used for the
analysis of IFS data and MUSE data in particular (e.g., López-Cobá et al.,
2020). pypipe3d performs a decomposition of the observed stellar spectra into
multiple simple stellar populations (SSPs), each of different age and
metallicity. The result is a set of maps comprising information of the stellar
populations (namely, velocity, velocity dispersion, metallicity among others),
and the ionized gas (namely, flux, velocity, velocity dispersion and their
errors). In addition to integrated properties such as the stellar mass. We
refer to Lacerda et al. (2022) for a thorough description of this analysis.

In this work we adopted the stellar libraries based on the extended MILES
(E-MILES, Vazdekis et al., 2016), FWHM=2.5 Å. We used 78 different templates
corresponding to 13 different ages (Age/Gyr = 0.03, 0.05, 0.08, 0.15, 0.25,
0.4, 0.6, 1, 1.75, 3, 5, 8.5, 13.5), and 6 metallicities
($\mathrm{[Z/H]}=-1.49,-0.96,-0.35,+0.06,+0.26,+0.40$). In regions where
E-MILES has a superior spectral resolution than the observed spectra, we
convolved the libraries to that of the MUSE line-spread function (e.g., Bacon
et al., 2017); in this way both templates and data share the same spectral
resolution within the range 4800 Å-7000 Å.

Unlike the ionized gas properties, the stellar kinematics analysis requires
spatial coadding of large bins to recover the stellar populations properties
at the outskirts of the FoV. While this has the advantage of increasing the
SN, including these pixels in the analysis and imposing a high SN would
compromise the spatial resolution of the data (particularly at the distance of
this object that results on an angular scale of 1.65 kpc arcsec-1) and affect
the final arrangement of bins as well as their individual sizes. All this
could distort the reconstruction of the kinematic maps and the modelling. We
thus have made a compromise between preserving the light distribution of the
object and maintaining a reasonable SN of the stellar spectra.

With this goal, before proceeding with the analysis we performed a SN
inspection of the datacube. In general, objects show a SN that decreases as a
function of the galactocentric distance and after certain radius the SN of the
continuum becomes too small to get reliable estimations of the stellar
population properties. In order to optimize the analysis, we selected a
rectangular region that contains only spaxels with a signal in the V-band
larger than 3. This procedure reduces the FoV of the observation to
$28\arcsec\times 19\arcsec$ in each wavelength channel. To increase the SN of
the stellar continuum in each spaxel from the resized cube, we performed a
Voronoi binning segmentation (Cappellari & Copin 2003) over a 2D-slice
centered around the V-band rest frame $\sim$5500$\pm$50 Å, with a SN goal of
30 per binned region.

This spatial binning is only used for the analysis of the stellar populations
but not for extracting the ionized gas kinematic properties, which were
obtained in a spaxel-wise sense for the purpose of comparison with the stars.
The H$\alpha$ flux and velocity maps are a subproduct of the SSP analysis with
pypipe3d. A velocity map is first created with single Gaussian fits. Having
the H$\alpha$ velocity recorded in each spaxel the line flux is then extracted
by creating a mask centred in H$\alpha\pm$200 km s-1. This method filters high
velocity gas. It isolates gas that traces more closely the systemic component
and, in principle, is less affected by extreme motions.

Our procedure does not consider the presence of multiple components.
Understanding that the assumption of single Gaussian fits is an
oversimplification (Villar Martín et al. 2014; Villar-Martín et al. 2018;
Harrison et al. 2015; V23 ), it helps towards tracing the dominant kinematic
component and filtering fainter, broad components, such as the outflow
emission. The potential impact of this will be discussed Sect. 4.2.1 in
comparison with V23 , who performed a multi-component analysis.

We also mapped the gas kinematics using [OIII]$\lambda$5007 following a
similar strategy. Because the results are basically identical, we only present
the H$\alpha$ kinematics and refer the reader to Harrison et al. (2015) and
V23  for a thorough analysis based on [OIII]$\lambda$5007.

Figures 3 shows the stellar and H$\alpha$ velocity maps from the SSP analysis
(Sect. 3.3). The kinematic modelling was performed using the $\mathtt{XS}$
tool (e.g. López-Cobá et al., 2024). This code implements Bayesian methods to
estimate the best disk geometry, namely the inclination $i$, the kinematic
position angle $\phi^{\prime}$ and the kinematic centre, as well as the
rotational velocities that better describe the observed velocity field. It
assumes the flat disk approximation to create a two dimensional interpolated
model. We refer to López-Cobá et al. (2024) for a thorough description of this
analysis. The SSP analysis gives a stellar mass of $\mathrm{\log
M_{\star}/M_{\odot}=11.1}$, which is consistent with Ramos Almeida et al.
(2022).

Figure 4: Maps of (O/H)/(O/H)⊙ ratio and $T_{\rm[OIII]}$. The masks used in
this analysis are shown in the left panel. They correspond to the apertures
for which [OIII]$\lambda$4363 is detected. Notice that the FoV is smaller than
in Fig. 2 and 2. The exact values of O/H and $T_{\rm[OIII]}$ are in Table 2.
$T_{\rm[OIII]}$ is in units of 104 K.  Figure 5: Map of (N/O)/(N/O)⊙. The FoV
is the same as in Fig. 4. The nuclear value is 1.99, which is outside the
color bar range to enhance the contrast for visualization purposes. The
(N/O)/(N/O)⊙ values are in Table 2.

##  4 Results

###  4.1 [OIII] and [NII] electron temperatures and gas chemical abundances

We show in Table 2 the O/H and N/O abudances predicted by photoionization
models for different apertures, as well as the [OIII] and [NII] electron
temperatures inferred with pyNeb for the densities estimated with
[SII]$\lambda\lambda$6716,6731 (see Sect. 3.2). The observed extinction
corrected line ratios relative to H$\beta$ and measured with the 1dim spectra
are also shown. This analysis was performed for a set of 23 (including the
nucleus) apertures where [OIII]$\lambda$4363 has been detected (see left panel
in Fig. 4). The top (lower) part of the table comprises the apertures in which
[NII]$\lambda$5755 is detected (undetected). This does not affect the
estimated chemical abundances, but simply indicates whether $T_{\rm[NII]}$
could be inferred or not. The O/H, N/O and $T_{\rm[OIII]}$ information in
Table 2 is also shown as 2dim maps in Fig. 4 and 5.

The main result which stands out in Fig. 4 is that the giant ionized bubble
shows significantly higher O/H and lower $T_{\rm[OIII]}$ than the gas at most
locations across the giant nebula. This is confirmed in all apertures along
the bubble edge, where $12+log$(O/H) is in the range
8.68$\pm$0.09-8.77$\pm$0.12. As in the nucleus (8.66$\pm$0.11), these
abundances are consistent with the solar value or slightly higher.

For the rest of the nebula (this is, the extended gas outside and beyond the
bubble), except in Ap. 13, the oxygen abundances are subsolar everywhere (as
tentatively found by Villar-Martín et al. 2018). Considering all apertures
across the nebula, the $log$ of the median is 12+$log$(O/H)=8.49 (or 63%
solar), with values as low as 8.37$\pm$0.08-8.44$\pm$0.13 (48%-56% solar) in
Ap. 15, 16, 19, 21, 22. Interestingly, the gas encircled by the bubble edges
(Ap. 11), with 12+$log$(O/H)=8.53$\pm$0.17 (69% solar), shows an intermediate
abundance between the most enriched gas (the bubble edges) and the more metal
poor gas in the rest of the nebula.

It is important to highlight that the abundance predictions are based on the
assumption that the gas is photoionized by the AGN everywhere. This is
reasonable for most apertures except, possibly, at a few positions located in
the direction perpendicular to the main axis of the nebula, where shocks might
be contributing to the excitation of the gas (V23 ). These apertures are
labelled Ap. 12, 18, 19 in Fig. 4, left panel. In spite of this uncertainty,
the fact that (except for Ap. 13) this gas follows the general behaviour of
the O/H and $T_{\rm[OIII]}$ maps suggests that the derived abundances are
reliable.

In general, the $T_{\rm[OIII]}$ map (right panel in Fig. 4) mimics the
behaviour of the O/H maps in an anticorrelation. The nucleus shows the minimum
temperature, $T_{4}$=$T_{\rm[OIII]}$/104=1.27$\pm$0.08 K, in comparison with
the rest of the gas, which shows $T_{4}>$1.5 everywhere. As before, the bubble
stands out in this map, with the edges being significantly colder
($T_{4}\sim$1.5-1.6 in Ap. 3 to 9) than the giant nebula, which is very hot
with $T_{4}\sim$1.7-1.9 almost everywhere. V23  presented a 2dim
$T_{\rm[OIII]}$ map tracing the nucleus and the ionized bubble. They inferred
a narrow range of significantly lower $T_{\rm[OIII]}\sim$(1.3-1.4)$\sim$104 K
compared to our estimations. The reason for this discrepancy remains unknown.

In regard to the N/O abundance (Fig. 5 and Table 2), the most obvious result
is the much higher nuclear value ($log$(N/O)=-0.56$\pm$0.08, i.e. roughly
twice the solar abundance ratio of -0.86$\pm$0.07), compared with the solar or
somewhat below solar N/O everywhere else (-1.03 $\pm$0.09$\leq
log{\rm(N/O)}\leq$-0.84$\pm$0.20), including the bubble edge. There is not
such a striking trend of N/O with location as found for O/H. The bubble edge
shows among the lowest N/O ratio ($\sim$68-85% N/O⊙), although not unique in
comparison with other regions of the giant extended nebula. The bubble edge
N/O values are lower than the gas encircled by it (ap. 11), the nucleus and
the gas in between (ap. 1 and 2).

$T_{\rm[NII]}$ could be measured only for the nucleus and up to the bubble
edge (Table 2). The temperature difference between these two regions is less
pronounced than for $T_{\rm[OIII]}$. $T_{4\leavevmode\nobreak\
\rm[NII]}\sim$0.983$\pm$0.066 K in the nucleus and somewhat higher in the
bubble ($\sim$1.01-1.12 K).

###  4.2 Kinematics

####  4.2.1 Stars vs. ISM gas within the inner $r\leq 7.5\arcsec$

The H$\alpha$ and stellar velocity maps inferred from the SSP analysis are
shown in Figure 3. The rotation pattern of the giant nebula discussed in
previous works is apparent. The rotation velocity $V_{\rm rot}$ reaches a
maximum within the inner few arcsec and drops outwards (Gagne et al. 2014;
Harrison et al. 2015; Ramos Almeida et al. 2017; Villar-Martín et al. 2018;
Moiseev & Ikhsanova 2023; V23 ; see also Fig. 2, left panel).

Some distortions are also apparent. While the gas in the inner regions
($r\la$5$\arcsec$ from the nucleus) follows an ordered rotation as seen by the
iso-velocity contours, at increasing galactocentric distances the contours
become more elongated, revealing a disturbed gaseous disk and/or the presence
of multiple kinematic components along the line-of-sight. The diversity of
structural components apparent in the narrow band images (e.g. filaments,
clumps, tidal tails; Fig. 2), some of which show distinct kinematics, as well
as the effects of feedback along and perpendicular to the radio axis can be at
the origin of such deviations from perfect circular rotation at each
galactocentric radii. To account for these additional motions, multi- (rather
than single-) Gaussian fits could be applied to the H$\alpha$ profile (see
Sect. 3.3). V23  performed a thorough multi-component analysis. Remarkably,
the velocity maps of the different kinematic components (including that
associated with the outflow) follow similar velocity fields, consistent with
the one traced by the single-Gaussian fits. Although the broad width of the
lines is due to the impact of the feedback induced by the nuclear activity,
the overall velocity field shows a pattern that is typical of rotation.

A novel result from this work is the identification of a rotation pattern for
the stellar kinematics in the region where it could be recovered
($r\la$7.5$\arcsec$ or $\la$12 kpc from the AGN). This pattern is, moreover,
similar to the gaseous kinematic pattern observed within the same region (Fig.
3, ”LoS vel.” panels). The models further support the existence of a rotation
pattern for the gas and the stars, and the similarity between both: gas and
stars show similar $\phi^{\prime}$ and $V_{\rm rot}$ (Fig. 3 and Table 1), and
a small difference in inclination $i\sim$15$\degr$. In both cases, $V_{\rm
rot}$ rises sharply up to a $\sim$450 km s-1 at $r\sim$1.5$\arcsec$ and falls
less steeply down to $\sim$300 km s-1 outwards, up to the maximum distance
traced by the modelled data at $r\sim$7.5$\arcsec$. This suggests that both
components are kinematically coupled and are corotating.

Including large bins | $\phi^{\prime}$ | $i$ | $\mathrm{V_{sys}}$
---|---|---|---
Stars | 51$\pm 1^{\circ}$ | $26\pm 1^{\circ}$ | 25551$\pm$1
H$\alpha$ | 58$\pm 1^{\circ}$ | $38\pm 1^{\circ}$ | 25500$\pm$2
Masking large bins | $\phi^{\prime}$ | $i$ | $\mathrm{V_{sys}}$
Stars | 46$\pm 1^{\circ}$ | $27\pm 1^{\circ}$ | 25567$\pm$3
H$\alpha$ | 46$\pm 1^{\circ}$ | $42\pm 1^{\circ}$ | 25491$\pm$2
Table 1: Best parameters describing the ionized gas and stellar disks
geometry. $\phi^{\prime}$ is the position angle of the major axis line of
nodes measured from N to E, $i$ the inclination and $\mathrm{V_{sys}}$ is the
inferred systemic velocity, which is assumed to be constant throughout the
disk.

Previous studies derived different results (Moiseev & Ikhsanova 2023; V23 ).
Based on the same VLT MUSE integral field spectroscopy data set analysed in
the current work, V23  found no clear pattern for the stellar motions, but a
general tentative trend of more redshifted velocities towards the outskirts of
the stellar emission, and lower velocities in the inner, most luminous part of
the galaxy (see Fig. 1 in their Appendix).

We believe that this discrepancy is due to the adopted analysis method. As
explained in Sect. 3.3, unlike the ionized gas properties, the stellar
kinematics analysis is affected by the spatial coadding. We tried to alleviate
this problem by resizing the datacube in which we excluded regions with poor
SN (in our analysis the high SN bins have sizes $\lesssim 1\arcsec$). For
further confirmation, we compared our results before and after removing the
largest Voronoi bins with a mask. In both cases, two symmetric blue shifted
and redshifted halves sides are clearly observed in the stellar velocity map,
both aligned at a roughly constant position angle in the sky. This general
behaviour corresponds to that of a rotating disk. The results of both methods
are in good agreement (Table 1).

Moiseev & Ikhsanova (2023) similarly studied the [OIII]$\lambda$5007 and
stellar kinematics using long slit spectroscopy and 3D spectroscopic
observations at the 6-m telescope of the Special Astrophysical Observatory of
the Russian Academy of Sciences (SAO RAS) and the SCORPIO-2 instrument. While
their results for the gas are in reasonable agreement with ours (except for
$V_{\rm rot}$), those concerning the stars are very different (they are also
different in comparison with V23 ). They described the observed LoS ionized
gas velocity distribution within $r<25\arcsec$ using a model that assumes a
regular circular rotation of a tilted ring. They inferred
$\phi^{\prime}$=62$\pm 4^{\circ}$ and $i$=43$\pm 7^{\circ}$, in reasonable
agreement with our models and a radial dependent $V_{\rm rot}$ that reaches a
maximum of 250 km s-1 at $r\sim$5$\arcsec$ (well below our maximum $V_{\rm
rot}\sim$400 km s-1) at $\sim$1.5$\arcsec$) and drops down to $\sim$100 km -1
further out (see their Fig. 5).

They measured a near zero radial gradient of stellar LoS velocities at $r\la
4\arcsec$, with a hint of slow counter-rotation. Assuming the stars lie in a
disk, they infer an orientation from the stellar continuum isophotes
$\phi^{\prime}$=157$\pm 3^{\circ}$ at $r=18-20\arcsec$ (compared to
$\phi^{\prime}\sim$46-51∘ from our models). The stellar disk is therefore
significantly inclined, even orthogonal, to the gaseous disk, suggesting that
they are kinematically decoupled. The authors measured a near zero radial
gradient of stellar LoS velocities at $r\la 4\arcsec$, The reasons for the
discrepancy with respect to our results can be diverse, including their use of
a long slit spectrum along PA 60∘ to trace the stars (versus the 2dim spatial
coverage provided by the MUSE data).

Rotation seems clearly present in the MUSE data for both the stars and part of
the ISM (clearly, this does not exclude non radial motions). Furthermore, both
components seem to trace each other following a similar kinematic patterns at
$r\la$7.5$\arcsec$. However, the stellar rotational speed implied by the
models is unusually large. According to the Tully-Fisher relation for rotating
disk galaxies, for a galaxy with a stellar mass similar to the Teacup
(log($M_{\rm*}/M_{\odot}$=11.15$\pm$0.05, Ramos Almeida et al. 2022; see also
Sec. 3.3) a value $V_{\rm rot}\sim$260 km s-1 is expected (Simons et al.
2015). The past merger event underwent by the Teacup 1-2 Gyr ago (Keel et al.
2015) is unlikely to be at the origin of the large rotation velocity, since
for a given $M_{\rm*}$, mergers would most likely shift $V_{\rm rot}$ towards
lower values (Covington et al. 2010).

Uncertainties on the value and on the radial dependence of both $i$ and
$\phi^{\prime}$ can affect $V_{\rm rot}$. Our method assumes the disk is
described by constant values in these angles. If the inclination or the
position angle depended strongly on the galactocentric distance (Moiseev &
Ikhsanova 2023) for instance a warped disk, or if the dynamics were dominated
by non-circular motions, different results are expected.

Since independently of the precise nature of the stellar and gas dynamics,
both components appear to follow a similar kinematic pattern at
$r\la$7.5$\arcsec$, and this is better traced up to larger distances and on a
spaxel to spaxel basis, we run several additional tests with the gas. On one
hand, rather than leaving the angles $i$ and $\phi^{\prime}$ as free
parameters, we fixed them: first to the values suggested by the H$\alpha$
morphology ($i$=50$\degr$ and $\phi^{\prime}$=68$\degr$) and subsequently to
the values inferred by Moiseev & Ikhsanova (2023) for their tilted ring model
($i$=43$\degr$ and $\phi^{\prime}$=62$\degr$). In both cases $V_{\rm rot}$
rises up to $\sim$400 km s-1.

An additional test consisted in modelling the gas kinematics traced with
H$\alpha$ by using 3DFIT, the main routine of BBarolo’s tool, which can fit a
3D tilted-ring model to an emission-line data-cube (Di Teodoro & Fraternali
2015). The results are shown in Fig. 6. The models imply $i$= 48$\degr$ (not
too different from previous results), $V_{\rm sys}$= 25554 km s-1 (the maximum
difference for $V_{\rm sys}$ for all the models we considered is 72 km s-1)
and $\phi^{\prime}$ increasing by up to $\sim$60$\degr$. Most importantly, the
rotation curve (red points in the figure) shows a similar amplitude and
behaviour as our previous models: a sharp increase in the center of $V_{\rm
rot}$ up to almost 400 km s-1 that declines down to $\sim$200 km s-1 at larger
distances. The results do not change much if we use a wider velocity cut from
-600 to +600 km s-1.

Figure 6: Tilted ring model (red) of the H$\alpha$ kinematics (grey circles)
using the BBarolo routine 3DFIT. The uncertainties on the derived parameters
are shown in shadowed colors. See text for details.

To summarize, we propose that, independently of the models, both stars and a
fraction of the gas within $r\la$7.5$\arcsec$ follow a similar rotation
kinematic pattern. Under this hypothesis and with the goal of isolating the
rotating component from the non circular-motions, the data were fitted with
circular rotation models. An unrealistically high $V_{\rm rot}$ value is
obtained in the inner few kpc. This suggests potential problems in the
derivation of the model parameters (such as the inclination and geometry)
and/or an inadequate isolation from non circular motions. As explained by Di
Teodoro & Fraternali (2015), who discussed in depth the robustness and the
limits of BBarolo, the goodness of a fit is mainly determined by the
combination of four factors: the inclination of the galaxy disc, the spatial
resolution, the spectral resolution and the S/N. When the inclination is not
known (as is the case in our study), large uncertainties affect the results.
For $i\la$45$\degr$ the fitted rotation curves can become progressively more
uncertain due to the smaller rotational component along the line of sight and
the larger impact of inclination errors. Clearly, deeper wide FoV 2-dim
spectroscopy would be of great value to map and model the motions with higher
precision specially for the stars.

###  4.3 Large scale rotation

The results of modelling the nebular kinematics up to $\sim$20$\arcsec$ with
BBarolo’s 3DFIT are shown in Fig. 7. The large scale rotation is very clear,
as already mapped by Moiseev & Ikhsanova 2023; V23 . The novelty introduced
here is the modelling of both the line of sight velocities and the velocity
dispersion, $\sigma$. Our main goals are on the one hand, to compare the main
rotation pattern on large spatial scales with that at $\sim$7.5$\arcsec$ (see
previous section) and on the other, to study the residuals and consider the
potential implications on the nature of the non radial motions.

The main kinematic axis (grey lines in the middle panels) has
$\phi^{\prime}\sim$67$\degr$ (dashed lines). The wiggles at increasing $r$
might be the result of distortions due to the more patchy H$\alpha$ emission.
For comparison, at $\la$7.5$\arcsec$, within the stellar body, the range of
$\phi^{\prime}$ derived of the different model assumptions and authors is
$\sim$(46-62)$\degr$. The conclusion is that, in spite of the strong non
radial motions existing in the central region, this gas shares the same sense
of rotation as the giant nebula and a kinematic axis that differs at most by
$\la$20$\degr$. It thus appears that the kinematics of the giant nebula is
linked to that of the gas (and the stars, as a consequence) within the central
region. The implications will be discussed in Sect. 5.2.

Rich additional information is provided in the maps of residuals (data minus
model, right panels in Fig. 7). The spatial distribution and morphology of the
residuals are remarkably similar in the $V_{\rm LOS}$ and $\sigma$ maps. They
deviate from the perfect disk morphology and the general rotation pattern
expected from the model. They are filamentary and patchy and could well be
outside the dominant nebular disk.

Figure 7: Moment maps obtained from the cube analysis of the H$\alpha$ line
with BBarolo, using a mask centered in H$\alpha\pm$200 km s-1(same model as in
Fig. 6). From top to bottom: zero moment map, velocity and velocity dispersion
of H$\alpha$. Each column represents the data, the corresponding moment map
and the residual map, respectively. The black dashed line represents the
average position angle from the tilted ring model, namely $\phi=67^{\circ}$
while the ellipse oriented at the latter PA and inclination $i=38^{\circ}$.
The grey line in the moment 1 map shows the radial variation of the position
angle shown in Fig. 6. The kinematic center is shown with an “x” in the first
panel.

V23  studied the spatially resolved properties of the ionized outflow and
giant nebula in the Teacup by modeling the emission lines with up to three
Gaussians in two spatial dimensions. The second and third broader
($\ga$150-200 km s-1) components are morphologically associated with the radio
bubbles and with the outflow. The authors proposed that this gas is moving
radially away from the galaxy and is being evacuated as part of the bubbles
(no models were presented). They hint at the possibility that the whole nebula
has been ejected as a consequence of such process.

We propose instead that a significant fraction of the gas, out to the outer
edge of the nebula follows a large scale velocity rotation pattern. This may
include also part of the gas associated with the otuflow. Although the broad
lines indicate turbulence induced by such mechanism and non-circular motions
are clearly present, the gas may still trace a global rotation pattern that
follows the gravitational potential (Gagne et al. 2014; Villar Martín et al.
2014; Moiseev & Ikhsanova 2023). This is supported by our finding that the
stellar component appears to be rotating, with an axis that is similar to that
of the gas. It is also supported by the overall similarity between the
velocity fields shown by the different gas kinematic components (see Fig. 4 in
V23 ). Although the non-outflow (first, narrow Gaussian component) and the
outflow components show very different line widths, they follow rather similar
velocity fields.

##  5 Discussion

###  5.1 Positive feedback and metal enrichment

We have shown that the Teacup bubble, known to be driven by an AGN wind or the
small nuclear radio jet, is associated with obvious changes in the gas
abundance. The bubble edge shows significantly higher O/H (solar or slightly
solar, similar to the nucleus) in comparison with the subsolar O/H across the
rest of the nebula. Most likely as a consequence, the bubble edge is also
significantly less hot ($T_{4}\sim$1.5-1.6) than the rest of the nebula
(1.7$\la T_{4}\la$1.9). Therefore, the outflow appears to be causing a change
in the gas metal content from the nucleus up to $\sim$10 kpc.

This mechanism may also explain the nuclear deficit of O/H. For the Teacup
values of log($M_{\rm*}/M_{\odot}$)=11.15$\pm$0.05 and star forming rate
$SFR\sim$10 M⊙ yr-1 (Jarvis et al. 2020; Ramos Almeida et al. 2022),
12+log(O/H)$\sim$8.85 is expected, according to the mass-metallicity-SFR
($M_{\rm*}$-Z-SFR) relation by Pérez-Montero et al. (2013) or $\sim$8.8 for
the extrapolation of the $M_{\rm*}$-Z-SFR relation by Andrews & Martini (2013)
to high $M_{\rm*}$. Seyfert 2 galaxies with similar $M_{\rm*}$ tend to show
O/H consistent with these predictions (Pérez-Díaz et al. 2021). The Teacup
nucleus, on the contrary, shows lower than expected O/H (8.66$\pm$0.05), while
the nuclear N/O (-0.56$\pm$0.08, Table 2) is consistent with that expected for
its $M_{\rm*}$ ($log$(N/O)=-0.54, Pérez-Montero et al. 2013; Andrews & Martini
2013).

The dilution produced by inward flows of low-metallicity gas (for instance,
from the outskirts of the two merging galaxies) could explain the nuclear O/H
and N/O (Edmunds 1990; Rupke et al. 2010). The outflow is also an interesting
possibility. This scenario is supported by the similar O/H of the bubble edge
and the nuclear gas, which is moreover enhanced in comparison with the rest of
the nebula. On one hand, radiative outflows can couple more efficiently with
metals via resonance line scattering (e.g. Pauldrach et al. 1994; Arav et al.
1994; Higginbottom et al. 2024; see also Edmunds 1990). On the other, galaxies
show metallicity gradients such that the inner regions have higher abundances
(Searle 1971; Kewley et al. 2019; Maiolino & Mannucci 2019). Based on this, a
gradient is expected to exist within the $\sim$4 kpc diameter nuclear aperture
used in the Teacup analysis. The outflow has been generated in the inner
regions ($\la$1 kpc) close to the AGN (Harrison et al. 2015; Ramos Almeida et
al. 2017; V23 ), where the gas is expected to be more metal rich. If it drags
gas out to large distances, as proposed in this scenario, the global
metallicity of the residual gas within the nuclear aperture would be lower as
a consequence.

An implication of this scenario is that the outflow has been capable of
ejecting gas from the galaxy center and dragged it up to $\sim$10 kpc. If the
outflow expanded without displacing significant amounts of gas to large
distances, the bubble abundance would be similar to the rest of the nebula.
The implications are important. It supports that metal-enriched galactic
outflows (driven by an AGN in this particular case) shape the mass-metallicity
relationship, by removing metals from galaxy potential wells and ejecting them
to large distances, possibly out into the CGM (Tremonti et al. 2004; Peeples
et al. 2014; Chisholm et al. 2018; Tortora et al. 2022).

The behaviour of N/O remains to be explained. An outflow could preserve N/O
(Edmunds 1990 (specific modeling would be valuable for AGN generated winds).
On the contrary, this ratio is depleted in the bubble edge in comparison with
the nucleus and its gradient shows no obvious spatial correlation with the
bubble, but just a tentative trend to show among the lowest N/O. This is not
necessarily a discrepancy. Given the complexity of the N/O behaviour in terms
of secondary and primary stellar production processes, and the fact that the
ejected gas would mix with gas across the nebula with a non-uniform N/O
distribution, it is difficult to predict how this ratio would behave as the
bubble expands and mixes with the pre-existing reservoir.

An alternative scenario to gas ejection from the center is that local chemical
enrichment has been produced by young stars. This is supported by the
detection of blue-coloured continuum emission co-spatial with the bubble edge
due to a population of stars that are younger ($\la$100-150 Myr) than in the
rest of the galaxy ($\ga$0.5-1 Gyr, V23 ). According to these authors,
widespread star formation has been triggered at the edge of the bubble due to
the compressing action of the jet and outflow (positive feedback). The time
scale could be long enough to enrich the local gas with O, but not with
secondary N (Mollá et al. 2006; Kumari et al. 2018). This would explain the
enhanced O/H in comparison with the rest of the nebula, while N/O is not
clearly different, and being at the same time significantly lower than in the
nucleus.

Yet another possible scenario to explain the abundance values in the bubble
relates to the depletion by dust of metals from the gas-phase ISM. If shocks
destroyed dust as the bubble expands (e.g. Dopita et al. 2018), metals could
be released to the gas. Since oxygen is more sensitive to depletion than
nitrogen (Jones & Ysard 2019), this could explain the higher O/H of the bubble
in comparison with the rest of the nebula, while having a tentatively lower
N/O.

Whether due to nuclear gas ejection, to local star formation or dust
destruction, in all three scenarios the AGN induced outflow is responsible for
the metal enrichment of the gas at distances as large as $\sim$10 kpc. The
implications are different. In the ejection scenario, the behaviour of O/H
provides observational evidence of how AGN induced outflows can deprive the
central regions of galaxies from metals and transfer them up to very large
distances, possibly out of the galaxy and into the CGM. In the second
scenario, the behaviour of O/H provides evidence of how AGN induced outflows
can produce local metal enrichment (this is also the case in the third, the
dust depletion scenario) at very large distances from the nucleus, with a
delay between the quasar onset and the induced metal enrichment of
$\sim$100-150 Myr.

###  5.2 The link between the giant nebula and the circumgalactic medium

With a maximum extension as traced by the MUSE data of $\sim$126 kpc (this is
a lower limit since the gas fills the FoV in this direction), the Teacup giant
nebula traces part of the CGM. Its properties are strongly conditioned by the
AGN, but it still provides valuable information about the CGM. If it was not
for the nuclear activity, most (if not all) of this gas would remain
invisible.

The giant nebula shows subsolar abundances, with
O/H$\sim$(48%-84%)$\times$(O/H)⊙ almost everywhere, and median 63%, well below
the nuclear, roughly solar abundance. For comparison, different works based on
absorption line studies have shown that the dense gas in the CGM of $z\la$1
galaxies has a bimodal metallicity distribution function, with an equal number
of absorbers in the low-metallicity ($Z\la$0.03$Z_{\odot}$) and high-
metallicity ($Z\sim$0.4$Z_{\odot}$) branches (e.g. Lehner et al. 2013; Wotta
et al. 2016). The abundance of the Teacup ionized nebula falls in the latter
group (clearly, this does not discard the presence of lower metallicity gas).
The high-metallicity branch has been proposed to trace galactic winds,
recycled outflows, and tidally stripped gas (Lehner et al. 2013).

Lower CGM metallicities compared with the ISM have been found for different
galaxy types, based on studies of absorption line systems at $\la$200 kpc from
their host galaxies (e.g. Kacprzak et al. 2019). These authors found a large
range of abundances in the CGM of isolated star forming galaxies, 0.01$\la$
Z/Z${}_{\odot}<1$ and an offset of $log(dZ)$=-1.17$\pm$0.11 between the CGM
and ISM, which shows no dependence with stellar mass. The relation has a large
scatter of 1$\sigma$=0.72. The offset for the Teacup is smaller for O/H
($log(dZ)\sim$-0.17 using the median nebular O/H), although still within the
scatter. Different processes may be at work in this system, related to the
nuclear activity and/or to its merger history.

We have proposed that part of the giant nebula is rotating (see also Villar-
Martín et al. 2018; Moiseev & Ikhsanova 2023, out to its outer edge. We
consider that this interpretation is independent of the fact that emission
lines trace only part of the total gas reservoir, or the fact that prominent
non radial motions are also present. Because the quasar preferentially
illuminates the gas within the ionization cones, the true geometry of the
giant gas reservoir is unknown. In spite of this uncertainty, as already
mentioned, much of the line emission, in and outside the cones, follows a
general velocity pattern strongly reminiscent of rotation (see our Fig. 2 and
specially the more detailed study by V23 ).

Evidence for giant rotating nebulae ($r\ga$several$\times$10 kpc), sometimes
defined as “superdisks” (Gopal-Krishna & Wiita 2000), has been found at
different $z$ (up to z$\sim$4). With a few excpetions, most are associated
with powerful radio galaxies (van Ojik et al. 1996; Gopal-Krishna & Wiita
2000; Villar-Martín et al. 2006; Martin et al. 2015; Balmaverde et al. 2022;
see also Prescott et al. 2015). This is consistent with observations that
suggest a substantial rotational component in the CGM of numerous systems with
$V_{\rm rot}$ sometimes as large as $\sim$few$\times$100 km s-1, as the Teacup
nebula (Ho et al. 2017; French & Wakker 2020; Zabl et al. 2019; Nateghi et al.
2023).

We have seen that the kinematics of the giant nebula on the largest scale
appears to be linked to that of the gas and the stars within $r\la$12 kpc.
This suggests that the CGM and the ISM are kinematically coupled, at least
partially. Kinematic coupling (at least partial) between the CGM and the ISM
and the stellar disk has been inferred in numerous galaxies, with corotation
being identified sometimes at least out to 0.5 $r_{\rm vir}$, specially for
the low ionization CGM state gas (Ho et al. 2017; French & Wakker 2020; Zabl
et al. 2019; Nateghi et al. 2023). Simulations also find that the cool
($T<$105 K) CGM forms giant corotating structures in different types of
galaxies (e.g. DeFelippis et al. 2020; Ho et al. 2020; Huscher et al. 2021;
Zabl et al. 2019).

The Teacup giant nebula may be a remnant of the formation process (Villar-
Martín et al. 2003; Martin et al. 2015). According to the current lambda cold
dark matter ($\Lambda$CDM) cosmological model, galaxies form hierarchically
out of overdensities in the dark matter distribution. As the surrounding
intergalactic medium is channeled toward a growing galaxy, simulations predict
that the angular momentum of the infalling gas is redistributed onto the disk
and seeds the rotation of the galaxy. As this infalling gas continues to feed
the galaxies throughout their lifetimes, it is expected that the extended
gaseous halos will rotate in the same sense as both the galactic disks (as we
see for the Teacup) and the dark matter halos (e.g. Chen et al. 2003; Pichon
et al. 2011; Stewart et al. 2011; Ho et al. 2020; French & Wakker 2020).

Alternatively, the nebula could be a giant rotating disk formed ) as a
byproduct of the past merger (Gopal-Krishna & Wiita 2000; Villar-Martín et al.
2018; Balmaverde et al. 2022). Large scale (up to $\sim$200 kpc) HI disks have
been found around several low $z$ radio galaxies and have been interpreted as
the result of merger events that occurred $>$108 yr ago (e.g. Morganti et al.
2002a, b; Emonts et al. 2007, 2010). The rotational velocities of these large
HI disks can be very large, ranging from $\sim$200 km s-1(Emonts et al. 2006;
Struve et al. 2010 up to $\sim$400 km s-1(Emonts et al. 2008). We may be
witnessing a similar scenario for the Teacup quasar, in which the giant
rotating gaseous disk is ionized (rather than neutral).

In both scenarios, the gaseous curved filaments identified in Fig. 7 may be a
result of different processes, including spiralling filamentary flows that
could still be feeding the central galaxy (Martin et al. 2019), tidal features
from the merger event underwent by this system (Keel et al. 2012) and outflow
products associated with repeated AGN outbursts (Nelson et al. 2019; V23 ).

##  6 Conclusions

The giant ($\ga$126 kpc) nebula associated with the Teacup QSO2 at $z=$0.085
traces part of its CGM. Its properties are strongly influenced by the nuclear
activity up to the outer detected emission line regions, where it still
provides valuable information about the CGM surrounding the quasar host
galaxy. If it was not for the nuclear activity, most (if not all) of this gas
would remain invisible. This study is an example of the great potential of
studying giant nebulae to investigate in emission the CGM around active
galaxies at all redshifts.

The widely studied AGN driven outflow responsible for the well known ionized
bubble is enhancing the gas metal content (O/H) up to $\sim$10 kpc from the
AGN. This could be a consequence of the ejection of metal rich gas from the
nucleus. In such scenario, the Teacup provides observational evidence for how
AGN feedback can deprive the central regions of galaxies from gas and displace
metals out to very large distances, possibly out of the galaxy. It supports
that metal-enriched AGN outflows can shape the massâ metallicity relationship.
Alternatively, the O/H enrichment could have been produced locally by the
young stellar population formed in the bubble edge, possibly formed as a
consequence of positive feedback $\sim$100-150 Myr (V23 ). A third possibility
is the release of oxygen to the gas phase as a consequence of dust destruction
in the bubble by shocks triggered by the expanding outflow. In any of the
scenarios considered, the nuclear activity is the ultimate mechanism
responsible for the metal enrichment of the gas at large extranuclear
distances ($\sim$10 kpc).

The giant nebula shows subsolar metallicity almost everywhere, except the
bubble. We propose that the nebula has a rotational component that extends up
to its outer edge. Prominent non-circular motions are also present, probably
originated by a diversity of mechanisms, including the outflow and the merging
process. This is consistent with CGM studies based on absorption lines that
show evidence for a substantial rotation component in the CGM of numerous
galaxies, in addition to non-rotational motions. The rotational component of
the nebula appears to be kinematically coupled, at least partially, with the
galaxy ISM and the stellar dynamics. CGM-ISM coupling has also been inferred
from observations and simulations of the CGM surrounding galaxies. We propose
that we are witnessing a similar behaviour, with the CGM traced in emission,
rather than absorption.

## Acknowledgments

We are grateful to Luc Binette and Bjorn Emonts for very valuable feedback.
MVM and ACL research has been funded by grant Nr. PID2021-124665NB-I00 by the
Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/
10.13039/501100011033 and by ”ERDF A way of making Europe”. SC acknowledges
financial support from the Severo Ochoa grant CEX2021-001131-S and the
Ministry of Science, Innovation and Universities (MCIU) under grants
PID2019-106027GB-C41. EPM acknowledges financial support by project
Estallidos8 PID2022-136598NB-C32 (Spanish Ministry of Science and Innovation).
The Cosmology calculator by Wright (2006) has been used.

Table 2: Photoionization model predictions of the O/H and N/O relative abundances, the electron temperatures, $T_{\rm[OIII]}$ and $T_{\rm[NII]}$ and electron densities $n_{\rm e}$ inferred from the [SII]$\frac{\lambda 6716}{\lambda 6731}$ ratio. The observed, extinction corrected line fluxes used in the models are quoted relative to H$\beta$. E(B-V) values were inferred from H$\alpha$/H$\beta$. The spatial location of the apertures can be seen in Fig. 4. Apertures 3 to 9 cover the edge of the ionized bubble. The top (lower) part of the table comprises the apertures in which [NII]$\lambda$5755 is detected (undetected). N/A in the [SII]$\frac{\lambda 6716}{\lambda 6731}$ column means that the doublet ratio could not be measured accurately due to residual sky artefacts. The solar abundances are 12+$log$(O/H)=8.69$\pm$0.04 and $log$(N/O)=-0.86$\pm$0.07 for the Sun (Asplund et al. 2021). Aperture | E(B-V) |  [OIII]$\lambda$4363 |  [OIII]$\lambda$5007 |  [NII]$\lambda$5755 |  [NII]$\lambda$6584 |  [SII]$\lambda$6725 |  [SII]$\frac{\lambda 6716}{\lambda 6731}$ |  12+$log$(O/H) |  $log$(N/O) | $T_{\rm[OIII]}$ | $T_{\rm[NII]}$ | $n_{\rm e}$
---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  |  |  |  | K | K |  cm-3
1 |  0.074$\pm$0.012 |  0.144$\pm$0.014 |  6.000$\pm$0.215 |  0.019$\pm$0.001 |  1.437$\pm$0.114 |  1.598$\pm$0.047 |  1.316$\pm$0.056 |  8.68$\pm$0.15 |  -0.86$\pm$0.11 |  16565$\pm$875 |  9508$\pm$312 |  149$\pm$79
2 |  0.155$\pm$0.015 |  0.141$\pm$0.012 |  5.802$\pm$0.588 |  0.017$\pm$0.001 |  1.396$\pm$0.044 |  1.507$\pm$0.051 |  1.396$\pm$0.026 |  8.61$\pm$0.13 |  -0.85$\pm$0.11 |  16739$\pm$1196 |  9236$\pm$225 |  34$\pm$36
3 |  0.054$\pm$0.013 |  0.135$\pm$0.015 |  6.043$\pm$0.251 |  0.020$\pm$0.002 |  1.044$\pm$0.042 |  1.467$\pm$0.060 |  1.349$\pm$0.048 |  8.76$\pm$0.11 |  -1.03$\pm$0.09 |  16058$\pm$1071 |  11185$\pm$448 |  100$\pm$63
4 |  0.034$\pm$0.013 |  0.158$\pm$0.012 |  7.290$\pm$0.251 |  0.018$\pm$0.002 |  1.012$\pm$0.048 |  1.359$\pm$0.071 |  1.382$\pm$0.027 |  8.69$\pm$0.10 |  -1.00$\pm$0.08 |  15801$\pm$699 |  10822$\pm$543 |  51$\pm$34
5 |  0.063$\pm$0.018 |  0.142$\pm$0.008 |  7.199$\pm$0.322 |  0.016$\pm$0.002 |  0.991$\pm$0.062 |  1.196$\pm$0.056 |  1.343$\pm$0.033 |  8.77$\pm$0.12 |  -0.99$\pm$0.08 |  15064$\pm$493 |  10325$\pm$565 |  109$\pm$41
6 |  0.046$\pm$0.017 |  0.143$\pm$0.010 |  6.885$\pm$0.319 |  0.019$\pm$0.001 |  1.086$\pm$0.063 |  1.326$\pm$0.078 |  1.321$\pm$0.036 |  8.69$\pm$0.12 |  -0.98$\pm$0.11 |  15433$\pm$641 |  10656$\pm$344 |  141$\pm$52
7 |  0.101$\pm$0.014 |  0.162$\pm$0.009 |  7.342$\pm$0.323 |  0.019$\pm$0.001 |  1.070$\pm$0.067 |  1.419$\pm$0.059 |  1.341$\pm$0.036 |  8.68$\pm$0.09 |  -1.00$\pm$0.08 |  15938$\pm$501 |  10743$\pm$353 |  114$\pm$50
8 |  0.087$\pm$0.021 |  0.128$\pm$0.014 |  5.665$\pm$0.326 |  0.020$\pm$0.001 |  1.291$\pm$0.112 |  1.589$\pm$0.088 |  1.369$\pm$0.046 |  8.76$\pm$0.15 |  -0.93$\pm$0.11 |  16181$\pm$1037 |  10134$\pm$391 |  74$\pm$60
9 |  0.031$\pm$0.013 |  0.145$\pm$0.006 |  6.900$\pm$0.316 |  0.020$\pm$0.002 |  1.195$\pm$0.069 |  1.423$\pm$0.072 |  1.394$\pm$0.021 |  8.72$\pm$0.10 |  -0.94$\pm$0.10 |  15453$\pm$442 |  10495$\pm$516 |  40$\pm$33
10 |  0.105$\pm$0.016 |  0.115$\pm$0.019 |  4.557$\pm$0.208 |  0.020$\pm$0.002 |  1.169$\pm$0.079 |  1.632$\pm$0.074 |  1.411$\pm$0.016 |  8.56$\pm$0.20 |  -0.98$\pm$0.11 |  16888$\pm$1484 |  10580$\pm$530 |  $\la$35
Nucleus |  0.186$\pm$0.033 |  0.086$\pm$0.014 |  6.550$\pm$0.310 |  0.021$\pm$0.003 |  1.497$\pm$0.157 |  0.920$\pm$0.059 |  1.242$\pm$0.138 |  8.66$\pm$0.11 |  -0.56$\pm$0.08 |  12772$\pm$785 |  9825$\pm$658 |  301$\pm$270
11 |  0.000$\pm$0.046 |  0.146$\pm$0.021 |  5.941$\pm$0.387 | N/A |  1.271$\pm$0.183 |  1.455$\pm$0.082 |  1.335$\pm$0.017 |  8.53$\pm$0.17 |  -0.85$\pm$0.15 |  16899$\pm$1339 | N/A |  115$\pm$22
12 |  0.132$\pm$0.015 |  0.113$\pm$0.014 |  3.570$\pm$0.170 | N/A |  1.350$\pm$0.061 |  1.865$\pm$0.085 |  1.450$\pm$0.050 |  8.52$\pm$0.14 |  -0.87$\pm$0.11 |  19540$\pm$1235 | N/A |  16$\pm$34
13 |  0.000$\pm$0.038 |  0.083$\pm$0.023 |  3.197$\pm$0.325 | N/A |  1.378$\pm$0.166 |  1.973$\pm$0.228 |  1.417$\pm$0.054 |  8.71$\pm$0.24 |  -0.89$\pm$0.14 |  17766$\pm$2661 | N/A |  34$\pm$47
14 |  0.031$\pm$0.025 |  0.175$\pm$0.035 |  6.901$\pm$0.531 | N/A |  1.244$\pm$0.104 |  1.684$\pm$0.126 |  1.330$\pm$0.029 |  8.53$\pm$0.15 |  -0.94$\pm$0.11 |  17197$\pm$2124 | N/A |  129$\pm$42
15 |  0.077$\pm$0.041 |  0.198$\pm$0.024 |  7.222$\pm$0.808 | N/A |  0.914$\pm$0.115 |  1.480$\pm$0.183 |  1.316$\pm$0.025 |  8.40$\pm$0.12 |  -0.94$\pm$0.13 |  17989$\pm$1652 | N/A |  145$\pm$34
16 |  0.023$\pm$0.015 |  0.271$\pm$0.040 |  8.844$\pm$0.387 | N/A |  0.647$\pm$0.072 |  1.382$\pm$0.319 | N/A |  8.37$\pm$0.08 |  -1.00$\pm$0.11 |  18829$\pm$1648 | N/A | N/A
17 |  0.029$\pm$0.020 |  0.222$\pm$0.017 |  8.941$\pm$0.404 | N/A |  1.123$\pm$0.102 |  1.339$\pm$0.146 |  1.429$\pm$0.061 |  8.49$\pm$0.05 |  -0.87$\pm$0.12 |  16910$\pm$754 | N/A |  33$\pm$48
18 |  0.121$\pm$0.012 |  0.102$\pm$0.023 |  3.917$\pm$0.210 | N/A |  1.097$\pm$0.058 |  1.602$\pm$0.191 |  1.466$\pm$0.030 |  8.56$\pm$0.24 |  -0.93$\pm$0.13 |  17299$\pm$2055 | N/A |  $\la$75
19 |  0.027$\pm$0.044 |  0.090$\pm$0.028 |  2.626$\pm$0.481 | N/A |  1.245$\pm$0.283 |  1.125$\pm$0.378 |  1.386$\pm$0.039 |  8.39$\pm$0.30 |  -0.86$\pm$0.28 |  18925$\pm$3439 | N/A |  48$\pm$45
20 |  0.077$\pm$0.031 |  0.087$\pm$0.016 |  3.075$\pm$0.420 | N/A |  1.494$\pm$0.282 |  1.888$\pm$0.284 |  1.396$\pm$0.039 |  8.60$\pm$0.25 |  -0.84$\pm$0.20 |  18398$\pm$2835 | N/A |  41$\pm$45
21 |  0.126$\pm$0.059 |  0.201$\pm$0.033 |  8.234$\pm$1.289 | N/A |  0.716$\pm$0.126 |  0.947$\pm$0.164 |  1.372$\pm$0.028 |  8.44$\pm$0.13 |  -0.94$\pm$0.16 |  16972$\pm$1929 | N/A |  70$\pm$38
22 |  0.000$\pm$0.007 |  0.154$\pm$0.012 |  7.099$\pm$0.143 | N/A |  0.776$\pm$0.038 |  1.164$\pm$0.338 | N/A |  8.44$\pm$0.18 |  -0.97$\pm$0.17 |  15798$\pm$608 | N/A | N/A

## References

  * Andrews & Martini (2013) Andrews, B. H. & Martini, P. 2013, ApJ, 765, 140
  * Arav et al. (1994) Arav, N., Li, Z.-Y., & Begelman, M. C. 1994, ApJ, 432, 62
  * Asplund et al. (2021) Asplund, M., Amarsi, A. M., & Grevesse, N. 2021, A&A, 653, A141
  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, Proc. SPIE, 7735, 773508
  * Bacon et al. (2017) Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1
  * Balmaverde et al. (2022) Balmaverde, B., Capetti, A., Baldi, R. D., et al. 2022, A&A, 662, A23
  * Baum et al. (1988) Baum, S. A., Heckman, T. M., Bridle, A., et al. 1988, ApJS, 68, 643
  * Borisova et al. (2016) Borisova, E., Cantalupo, S., Lilly, S. J., et al. 2016, ApJ, 831, 39
  * Cappellari & Copin (2003) Cappellari, M. & Copin, Y. 2003, MNRAS, 342, 345
  * Chen et al. (2003) Chen, D. N., Jing, Y. P., & Yoshikaw, K. 2003, ApJ, 597, 35.
  * Chisholm et al. (2018) Chisholm, J., Tremonti, C., & Leitherer, C. 2018, MNRAS, 481, 1690
  * Covington et al. (2010) Covington, M. D., Kassin, S. A., Dutton, A. A., et al. 2010, ApJ, 710, 279
  * Cox et al. (2006) Cox, T. J., Di Matteo, T., Hernquist, L., et al. 2006, ApJ, 643, 692
  * DeFelippis et al. (2020) DeFelippis, D., Genel, S., Bryan, G. L., et al. 2020, ApJ, 895, 17
  * Di Teodoro & Fraternali (2015) Di Teodoro, E. M. & Fraternali, F. 2015, MNRAS, 451, 3021.
  * Dopita et al. (2018) Dopita, M. A., Vogt, F. P. A., Sutherland, R. S., et al. 2018, ApJS, 237, 10
  * Dors et al. (2015) Dors, O. L., Cardaci, M. V., Hägele, G. F., et al. 2015, MNRAS, 453, 4102
  * Edmunds (1990) Edmunds, M. G. 1990, MNRAS, 246, 678
  * Emonts et al. (2006) Emonts, B. H. C., Morganti, R., Tadhunter, C. N., et al. 2006, A&A, 454, 125
  * Emonts et al. (2007) Emonts, B. H. C., Morganti, R., Oosterloo, T. A., et al. 2007, A&A, 464, L1
  * Emonts et al. (2008) Emonts, B. H. C., Morganti, R., Oosterloo, T. A., et al. 2008, MNRAS, 387, 197
  * Emonts et al. (2010) Emonts, B. H. C., Morganti, R., Struve, C., et al. 2010, MNRAS, 406, 987
  * Emonts et al. (2016) Emonts, B. H. C., Lehnert, M. D., Villar-Martín, M., et al. 2016, Science, 354, 1128
  * Falkendal et al. (2021) Falkendal, T., Lehnert, M. D., Vernet, J., et al. 2021, A&A, 645, A120
  * Fossati et al. (2021) Fossati, M., Fumagalli, M., Lofthouse, E. K., et al. 2021, MNRAS, 503, 3044
  * French & Wakker (2020) French, D. M. & Wakker, B. P. 2020, ApJ, 897, 151
  * Gagne et al. (2014) Gagne, J. P., Crenshaw, D. M., Kraemer, S. B., et al. 2014, ApJ, 792, 72
  * Gopal-Krishna & Wiita (2000) Gopal-Krishna & Wiita, P. J. 2000, ApJ, 529, 189
  * Hani et al. (2018) Hani, M. H., Sparre, M., Ellison, S. L., et al. 2018, MNRAS, 475, 1160
  * Harrison et al. (2015) Harrison C. M., Thomson A.P., Alexander D. M., Bauer F.E., Edge A.C., Hogan M.T., Mullaney J. R., Swinbank A. M., 2015, ApJ, 800, 45
  * Hayes et al. (2016) Hayes, M., Melinder, J., Östlin, G., et al. 2016, ApJ, 828, 49
  * Higginbottom et al. (2024) Higginbottom, N., Scepi, N., Knigge, C., et al. 2024, MNRAS, 527, 9236
  * Ho et al. (2017) Ho, S. H., Martin, C. L., Kacprzak, G. G., et al. 2017, ApJ, 835, 267
  * Ho et al. (2020) Ho, S. H., Martin, C. L., & Schaye, J. 2020a, ApJ, 904, 76
  * Huang et al. (2016) Huang, Y., Liu, X.-W., Yuan, H.-B., et al. 2016, MNRAS, 463, 2623
  * Huscher et al. (2021) Huscher, E., Oppenheimer, B. D., Lonardi, A., et al. 2021, MNRAS, 500, 1476
  * Jarvis et al. (2020) Jarvis, M. E., Harrison, C. M., Mainieri, V., et al. 2020, MNRAS, 498, 1560
  * Jones & Ysard (2019) Jones, A. P. & Ysard, N. 2019, A&A, 627, A38.
  * Kacprzak et al. (2019) Kacprzak, G. G., Pointon, S. K., Nielsen, N. M., et al. 2019, ApJ, 886, 91
  * Keel et al. (2012) Keel, W. C., Chojnowski, S. D., Bennert, V. N., et al. 2012, MNRAS, 420, 878
  * Keel et al. (2015) Keel, W. C., Maksym, W. P., Bennert, V. N., et al. 2015, AJ, 149, 155
  * Kewley et al. (2019) Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511
  * Kumari et al. (2018) Kumari, N., James, B. L., Irwin, M. J., et al. 2018, MNRAS, 476, 3793
  * Lacerda et al. (2022) Lacerda, E. A. D., Sánchez, S. F., Mejía-Narváez, A., et al. 2022, New A, 97, 101895
  * Lehner et al. (2013) Lehner, N., Howk, J. C., Tripp, T. M., et al. 2013, ApJ, 770, 138
  * López-Cobá et al. (2020) López-Cobá, C., Sánchez, S. F., Anderson, J. P., et al. 2020, AJ, 159, 167
  * López-Cobá et al. (2024) López-Cobá, C., Lihwai, L., & Sánchez, S. F. 2024, Rev. Mexicana Astron. Astrofis., 60, 19
  * Luridiana et al. (2015) Luridiana, V., Morisset, C., & Shaw, R. A. 2015, A&A, 573, A42
  * Maiolino & Mannucci (2019) Maiolino, R. & Mannucci, F. 2019, A&A Rev., 27, 3
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
  * McCarthy et al. (1990) McCarthy, P. J., Spinrad, H., van Breugel, W., et al. 1990, ApJ, 365, 487
  * McCarthy (1993) McCarthy, P. J. 1993, ARA&A, 31, 639
  * Mannucci et al. (2010) Mannucci, F., Cresci, G., Maiolino, R., et al. 2010, MNRAS, 408, 2115
  * Martin et al. (2015) Martin, D. C., Matuszewski, M., Morrissey, P., et al. 2015, Nature, 524, 192
  * Martin et al. (2019) Martin, D. C., O’Sullivan, D., Matuszewski, M., et al. 2019, Nature Astronomy, 3, 822
  * Mollá et al. (2006) Mollá, M., Vílchez, J. M., Gavilán, M., et al. 2006, MNRAS, 372, 1069
  * Morganti et al. (2002a) Morganti, R., Oosterloo, T. A., Tinti, S., et al. 2002a, Seeing Through the Dust: The Detection of HI and the Exploration of the ISM in Galaxies, 276, 423
  * Morganti et al. (2002b) Morganti, R., Oosterloo, T. A., Tinti, S., et al. 2002b, A&A, 387, 830
  * Moiseev & Ikhsanova (2023) Moiseev, A. V. & Ikhsanova, A. I. 2023, Universe, 9, 66
  * Mukherjee et al. (2018) Mukherjee, D., Bicknell, G. V., Wagner, A. Y., et al. 2018, MNRAS, 479, 5544
  * Nateghi et al. (2023) Nateghi, H., Kacprzak, G. G., Nielsen, N. M., et al. 2023, arXiv:2311.05165
  * Nelson et al. (2019) Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234
  * Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei, by Donald E. Osterbrock. Published by University Science Books, ISBN 0-935702-22-9
  * Pauldrach et al. (1994) Pauldrach, A. W. A., Kudritzki, R. P., Puls, J., et al. 1994, A&A, 283, 525
  * Peeples et al. (2014) Peeples, M. S., Werk, J. K., Tumlinson, J., et al. 2014, ApJ, 786, 54
  * Pérez-Díaz et al. (2021) Pérez-Díaz, B., Masegosa, J., Márquez, I., et al. 2021, MNRAS, 505, 4289
  * Pérez-Montero et al. (2013) Pérez-Montero, E., Contini, T., Lamareille, F., et al. 2013, A&A, 549, A25
  * Pérez-Montero (2014) Pérez-Montero, E. 2014, MNRAS, 441, 2663
  * Pérez-Montero et al. (2019) Pérez-Montero, E., Dors, O. L., Vílchez, J. M., et al. 2019, MNRAS, 489, 2652
  * Pichon et al. (2011) Pichon, C., Pogosyan, D., Kimm, T., et al. 2011, MNRAS, 418, 2493
  * Pointon et al. (2019) Pointon, S. K., Kacprzak, G. G., Nielsen, N. M., et al. 2019, ApJ, 883, 78
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6.
  * Prescott et al. (2015) Prescott, M. K. M., Martin, C. L., & Dey, A. 2015, ApJ, 799, 62
  * Ramos Almeida et al. (2017) Ramos Almeida, C., Piqueras López, J., Villar-Martín, M., et al. 2017, MNRAS, 470, 964
  * Ramos Almeida et al. (2022) Ramos Almeida, C., Bischetti, M., García-Burillo, S., et al. 2022, A&A, 658, A155.
  * Rémy-Ruyer et al. (2014) Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31
  * Rupke et al. (2010) Rupke, D. S. N., Kewley, L. J., & Barnes, J. E. 2010, ApJ, 710, L156.
  * Searle (1971) Searle, L. 1971, ApJ, 168, 327
  * Sellwood & Spekkens (2015) Sellwood, J. A. & Spekkens, K. 2015, arXiv e-prints, arXiv:1509.07120
  * Simons et al. (2015) Simons, R. C., Kassin, S. A., Weiner, B. J., et al. 2015, MNRAS, 452, 986
  * Stewart et al. (2011) Stewart, K. R., Kaufmann, T., Bullock, J. S., et al. 2011, ApJ, 738, 39
  * Struve et al. (2010) Struve, C., Oosterloo, T., Sancisi, R., et al. 2010, A&A, 523, A75.
  * Tortora et al. (2022) Tortora, C., Hunt, L. K., & Ginolfi, M. 2022, A&A, 657, A19
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898
  * Tumlinson et al. (2017) Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389
  * Valdes et al. (2004) Valdes, F., Gupta, R., Rose, J. A., et al. 2004, ApJS, 152, 251
  * van Ojik et al. (1996) van Ojik, R., Roettgering, H. J. A., Carilli, C. L., et al. 1996, A&A, 313, 25
  * Vazdekis et al. (2016) Vazdekis, A., Koleva, M., Ricciardelli, E., et al. 2016, MNRAS, 463, 3409
  * (89) Venturi, G., Treister, E., Finlez, C., et al. 2023, A&A, 678, A127
  * Villar-Martín et al. (2003) Villar-Martín, M., Vernet, J., di Serego Alighieri, S., et al. 2003, MNRAS, 346, 273
  * Villar-Martín et al. (2006) Villar-Martín, M., Sánchez, S. F., De Breuck, C., et al. 2006, MNRAS, 366, L1
  * Villar Martín et al. (2014) Villar Martín, M., Emonts, B., Humphrey, A., et al. 2014, MNRAS, 440, 3202
  * Villar-Martín et al. (2018) Villar-Martín, M., Cabrera-Lavers, A., Humphrey, A., et al. 2018, MNRAS, 474, 2302
  * Villar Martín et al. (2021) Villar Martín, M., Emonts, B. H. C., Cabrera Lavers, A., et al. 2021, A&A, 650, A84
  * Wang et al. (2023) Wang, W., Wylezalek, D., Vernet, J., et al. 2023, A&A, 680, A70
  * Wotta et al. (2016) Wotta, C. B., Lehner, N., Howk, J. C., et al. 2016, ApJ, 831, 95
  * Wright (2006) Wright, E.L. 2006 PASP, 118, 1711
  * Zabl et al. (2019) Zabl, J., Bouché, N. F., Schroetter, I., et al. 2019, MNRAS, 485, 1961
  * Zinchenko (2023) Zinchenko, I. A. 2023, A&A, 674, L7