© The Authors 2023

1. Introduction

J1430+1339, known as the Teacup (Keel et al. 2012), is an active galaxy at z ≃ 0.08506 hosting a highly obscured (N_H ∼ 5 × 10²³ cm⁻²; Lansbury et al. 2018) type-2 quasar (QSO; Villar Martín et al. 2014; Harrison et al. 2014), with a bolometric luminosity of L_AGN ∼ 10^{45 − 46} erg s⁻¹ (Harrison et al. 2015; Lansbury et al. 2018). With a mass of log (M_⋆/M_⊙)≃11.15 (Ramos Almeida et al. 2022) and a star formation rate (SFR) of ∼8–12 M_⊙ yr⁻¹ (Jarvis et al. 2020; Ramos Almeida et al. 2022), corresponding to a specific SFR of log ((SFR/M_⋆)/Gyr⁻¹)∼ − 1.2 to –1.1, the galaxy lies ∼0.6–0.8 dex above the “main sequence” of star-forming galaxies for its redshift and stellar mass (e.g. Brinchmann et al. 2004; Speagle et al. 2014; Popesso et al. 2023).

Its nickname, the Teacup, is due to the presence of a closed loop of ionised gas to the E of the nucleus, resembling the handle of a teacup, with the main body of the galaxy in continuum emission representing the cup itself (Fig. 1). The ionised gas “handle” is particularly prominent in [O III], extending up to distances of ∼10 kpc (Keel et al. 2012; Harrison et al. 2015), and it is dominated by ionisation from the AGN photons, as indicated by longslit observations (Keel et al. 2015). The ionised nebula extends up to ∼140 kpc (∼70 kpc per side; Villar Martín et al. 2018, 2021; Moiseev & Ikhsanova 2023). This lies among the largest ionised nebulae around active galaxies known so far, especially at low z (see Villar Martín et al. 2018 and references therein). The galaxy is bulge-dominated and shows disturbed and shell-like features in optical continuum emission, indicative of past merger activity (Keel et al. 2015). The source is classified as radio quiet in the νL_ν(1.4 GHz)–L_[O III] plane, having [O III] and radio luminosities of L_[O III] ∼ 5 × 10⁴² erg s⁻¹ and L_1.4 GHz ∼ 5 × 10²³ W Hz⁻¹, respectively (Harrison et al. 2015; Jarvis et al. 2019). Nevertheless, it exhibits lively radio activity, with bipolar lobes in the E–W direction on scales of ∼10 kpc per side (Harrison et al. 2015; Jarvis et al. 2019). Integral field spectroscopic (IFS) observations allowed for a kiloparsec-scale high-velocity component (v ∼ −800 km s⁻¹) of the ionised gas to be identified in both the optical (Harrison et al. 2014, 2015) and the near-IR (Ramos Almeida et al. 2017) bands, which has been interpreted as an outflow, after a broad component indicative of an outflow was found in [O III] in the SDSS spectrum of the galaxy (FWHM¹ ∼ 1000 km s⁻¹; Mullaney et al. 2013; Villar Martín et al. 2014).

Fig. 1. Optical (emission lines and continuum) versus radio emission of the Teacup galaxy. Left panel: false-colour image from VLT/MUSE. [O III] ionised gas emission is reported in green, Hα in red, and line-free continuum emission averaged between 5500–6700 Å (observed wavelengths) in blue. The colour intensity scale is the same for [O III] and Hα. The main spatial features mentioned in the text are labelled. Right panel: 20″ × 20″ zoomed-in [O III] emission (same as in left panel), with the contours of the VLA 5.12 GHz (white) and highest-resolution 6.22 GHz (magenta) radio images from Harrison et al. (2015) overlaid. For the former, contour levels are 35.0, 66.0, 124, and 234 μJy beam−1; for the latter, they are 40, 102, 261, and 668 μJy beam−1. The black numbered “+” symbols mark the regions from which the spectra reported in Fig. 2 were extracted.

Specifically, Harrison et al. (2015) performed a comparison between the optical and radio emissions of the Teacup galaxy, exploiting IFS observations from the European Southern Observatory’s (ESO) VIsible MultiObject Spectrograph (VIMOS) at the Very Large Telescope (VLT), Hubble Space Telescope (HST) imaging, and multi-band Karl G. Jansky Very Large Array (VLA) radio imaging. The radio lobes appear as bubbles filled with diffuse emission. The eastern radio bubble brightness peaks on its eastern edge, where it traces the shape of the ionised gas arc. The ionised gas is much fainter at the location of the western radio bubble, potentially explained by a radio-optical ionising flux misalignment or by a lower gas density. The ionised gas kinematic map shows a ∼20° misalignment between the kinematic axis and the radio bubbles. The highest resolution radio images (HPBW² ∼0.6 kpc) reveal that the central core is constituted by two unresolved sources, a brighter central one and a fainter one at ∼0.8 kpc to the NE, named HR-A and HR-B by Harrison et al. (2015), respectively. These two unresolved sources have been interpreted as likely being a small-scale radio jet, with HR-B being a hot spot (as commonly observed in compact jets; e.g. Leipski et al. 2006; Morganti et al. 2015). Harrison et al. (2015) interpret the ∼10 kpc per side radio bubbles to be the result of the action of the compact radio jet, based on the analogy with other similar cases, such as Mrk 6 (7.5 kpc radio bubbles and a ∼1 kpc radio jet; Kharb et al. 2006), as predicted by jet-driven bubble models (e.g. Sutherland & Bicknell 2007; Wagner et al. 2012), or to nuclear AGN winds, as in the energy-conserving expanding bubble models (e.g. Faucher-Giguère & Quataert 2012; Zubovas & King 2012; Wagner et al. 2013). High-resolution Chandra X-ray imaging reveals a loop co-spatial with the eastern ionised gas handle and the edges of the radio bubble, whose X-ray spectral properties are consistent with shocked thermal gas (Lansbury et al. 2018), supporting a jet or AGN wind origin for the bubble inflating mechanism.

Near-infrared (NIR) K-band observations with the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI) at VLT (Ramos Almeida et al. 2017) reveal a nuclear blueshifted broad component (FWHM  ∼  1600 − 1800 km s⁻¹) in both the hydrogen recombination lines (Paα, Brδ, and Brγ) and the [Si VI] λ1.963 μm coronal line. The larger line width of the ionised gas NIR lines compared to their optical counterparts in the same nuclear region suggests that the NIR lines, which are less affected by extinction compared to the optical ones, provide information on the outflowing gas closer to the AGN. Ramos Almeida et al. (2017) found the presence of ionised gas from Paα with FWHM > 250 km s⁻¹ across the NE bubble and SW fan, out to 5.6 kpc from the nucleus, indicative of an extended outflow. Here, the above work also reported the tentative detection of an additional very broad (FWHM ∼ 3000 km s⁻¹) component.

Atacama Large Millimeter/submillimeter Array (ALMA) CO(2–1) observations of the central 3″ × 3″ found a 1.3 kpc-scale rotating disc of molecular gas oriented approximately in the N–S direction, with mass M_H2 ∼ 6 × 10⁹ M_⊙ (consistent with the value inferred in Jarvis et al. 2020 from single-dish APEX observations), and a molecular outflow in the E–W direction with velocities up to 250 km s⁻¹, a mass of M_out = 3.12 × 10⁷ M_⊙, and a (de-projected) mass outflow rate of Ṁ_out = 15.8 M_⊙ yr⁻¹ (Ramos Almeida et al. 2022). Combining the same CO(2–1) data with additional new ALMA CO(3–2) observations, Audibert et al. (2023) found evidence for jet-induced molecular gas excitation and turbulence perpendicular to the compact radio jet. They also explored different outflow scenarios, from more to less conservative, and reported corresponding molecular mass outflow rates between 15 and 41 M_⊙ yr⁻¹.

In this work, we present new optical integral field observations from the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) at VLT, which allowed us to simultaneously map the distribution, kinematics, and excitation properties of the ionised gas out to large scales, ∼100 kpc (corresponding to the 1′ × 1′ field of view of MUSE at the galaxy distance, D_L ∼ 404 Mpc³; the physical scale is 1.66 kpc arcsec⁻¹). We compared the MUSE observations with the VLA radio images from Harrison et al. (2015), to study the interplay between the radio-emitting material and the ionised gas in the optical.

In Sect. 2 we describe the reduction and analysis of the VLT/MUSE observations, in Sect. 3 we present and comment on the MUSE maps of ionised gas obtained from the previous step, in Sect. 4 we characterise and discuss the spatially resolved properties of the galactic outflow, in Sect. 5 we discuss the broad velocity dispersions detected perpendicular to the jet and AGN ionisation axes, in Sect. 6 we present possible evidence for “positive” AGN feedback, and finally in Sect. 7 we summarise the conclusions of this work. Throughout this study, a H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7 cosmology is adopted. We adopt the standard astronomical orientation in the maps (N is up, and E is to the left).

2. Data reduction and analysis

The Teacup was observed with VLT/MUSE in its wide field mode (WFM), which covers ≃1′ × 1′ with a spatial sampling of 0.2″ spaxel⁻¹, and extended spectral mode, spanning the range 4600–9350 Å. It was observed as part of two observing programmes, a seeing-limited one (ID 0102.B-0107, PI L. Sartori) and a seeing-enhancer ground layer adaptive optics (GLAO)-assisted one (ID 0103.B-0071, PI C. Harrison). In this work we only employ the observations from the former programme, the seeing-limited one (∼0.5″ − 0.6″ at zenith), since visually assessing that their spatial quality was significantly better than that of the GLAO-assisted ones, even after the GLAO correction, due to the worse seeing conditions of these latter ones (∼1.0″ − 1.1″ at zenith). The employed observations consist of six observations of 950 s on-target each, for a total exposure time of 5700 s. Consecutive observations were rotated by 90° relative to one another.

Regarding the VLA radio observations that we compare with those from VLT/MUSE, we employ the C-band 5.12 GHz image, with a resolution of HPBW = 1.13″ × 1.04″, and the high-resolution (HPBW = 0.37″ × 0.24″) 6.22 GHz image. For the reduction of these radio data, we refer to Harrison et al. (2015) where they were originally presented.

2.1. Data reduction

We performed the data reduction using the ESO MUSE pipeline (Weilbacher et al. 2020; version 2.8.1), by employing the software ESO Reflex (Recipe flexible execution workbench; Freudling et al. 2013), that gives a graphical and automated way to execute the Common Pipeline Library (CPL; Banse et al. 2004; ESO CPL Development Team 2014) reduction recipes with EsoRex (ESO Recipe Execution Tool; ESO CPL Development Team 2015), within the Kepler workflow engine (Altintas et al. 2004). The data were corrected for bias, flat field, illumination, wavelength calibration, flux calibration, geometric reconstruction of the data cube, sky subtraction, and exposure combination. In order to maximise the on-target time, the observations did not include dedicated sky observations, such that the sky background was extracted from empty regions in the science field itself. Notably, the default pipeline sky extraction, which automatically creates a spatial mask for sky extraction based on the flux of the collapsed continuum image, did not give optimal results for the analysis of this complex system. The morphology of the continuum and ionised gas line emission is very different in the Teacup (see Fig. 1), and the automatic masking based on the continuum left regions containing significant gas line emission unmasked, that resulted in an extracted sky spectrum containing gas emission lines, in particular [O III] λλ4959,5007, Hβ, Hα, [N II] λλ6548,6584, and [S II] λλ6716,6731. This resulted in subtraction of these emission lines from the cube together with the sky. We thus followed a different approach, by creating a custom sky mask based on [O III] λ5007 (since this is the strongest and most extended emission line in the Teacup). The mask included the 10% weakest spaxels in [O III] flux, after having smoothed the [O III] image with a Gaussian kernel having a 3-spaxel standard deviation. The final sky mask was obtained with a logical “AND” between this [O III] mask and the mask produced by the pipeline that included 15% weakest spaxels in the white-light continuum image, so as to select only spaxels weak in both [O III] and continuum. The resulting data cube after having employed this custom sky mask for the sky subtraction showed no over-subtraction of emission lines of the source.

Unfortunately, this final data cube contained sky emission line residuals close to the wavelengths of the [S II] doublet which potentially affected its analysis (the other emission lines of interest were instead not affected by sky residuals), so we produced another cube to be used exclusively for the analysis of [S II] (details in Sect. 2.2). This cube was obtained by running the Zurich Atmosphere Purge (ZAP; Soto et al. 2016) code on the reduced data cube just described. ZAP is a Python code aimed at removing residual sky contamination from the MUSE data cubes by performing principal components analysis (PCA). We preferred using the reduced data cube (without ZAP execution) for the rest of the analysis, since PCA introduces further noise in the data and can also create artefacts. We checked that [S II] was not affected by artefacts and that its line profiles were consistent among the two cubes.

2.2. Data analysis

Data analysis was carried out by employing custom Python scripts, with an approach similar to that detailed in Venturi et al. (2018, 2021) and Mingozzi et al. (2019), summarised below.

Given that the MUSE spectra of the Teacup show stellar continuum emission underlying the ionised gas emission lines, we first modelled and subtracted it in order to separate these two components and allow for an easier subsequent modelling of the emission lines. Before doing that, we applied a Voronoi adaptive binning (Cappellari & Copin 2003) to the data cube to reach a minimum signal-to-noise ratio (S/N) per bin on the stellar continuum. We set an average value of 35 per wavelength channel for the minimum S/N to be achieved per bin, considering the rest-frame continuum in the range 4260–5530 Å which includes prominent stellar absorption features, and excluding gas emission lines. We modelled the stellar continuum in each bin across the whole spectral range of MUSE in its extended mode (4600–9350 Å, corresponding to rest-frame stellar templates in the range ∼4240–8570 Å given the redshift of the source, z ≃ 0.08506), by making use of the E-MILES single-stellar population (SSP) model spectra (Koleva & Vazdekis 2012; Vazdekis et al. 2016; La Barbera et al. 2017), which cover the spectral range 1680–50000 Å. The employed models adopt a unimodal standard Salpeter IMF with slope 1.3 (Salpeter 1955) as described in Vazdekis et al. (1996) and scaled-solar Padova+00 isochrones (Girardi et al. 2000), which have M/H metallicities spanning from –1.71 to +0.22 (in log, with respect to the solar value) and ages between 0.063 and 17.8 Gyr. The continuum modelling was done using the penalized pixel-fitting code (PPXF; Cappellari & Emsellem 2004; Cappellari 2017), which convolves the linearly combined stellar templates with a Gaussian function to reproduce the velocity and velocity dispersion of the absorption features. The main gas emission lines, as well as the residuals of not optimally subtracted sky lines, were masked during the stellar fitting, and underlying absorption features were extrapolated over in order to model the full spectral range. We employed both an additive and a multiplicative polynomial of degrees 2 and 5, respectively, to account for deformations of the spectra with respect to the templates (e.g. because of reddening).

We then subtracted the modelled stellar continuum on a spaxel basis, by re-scaling the model obtained for a given Voronoi bin to the spectral median of the observed continuum in each spaxel belonging to that bin. From this continuum-subtracted cube, we modelled the ionised gas optical emission lines of interest, namely Hγ, [O III] λ4363, Hβ, [O III] λλ4959,5007, [O I] λλ6300,6363, Hα, [N II] λλ6548,6584, and [S II] λλ6716,6731. Before proceeding with the line modelling, we produced an alternative cube, by applying a Voronoi binning to the continuum-subtracted cube around Hβ (the signal and noise of the spectral channels in an interval of ∼20 Å encompassing the line have been considered), requesting an average S/N per spectral channel of at least 5, in order to retrieve the Hβ emission from a wider area of the galaxy compared to the unbinned case. We chose Hβ because it is used for most of the diagnostic maps employing emission line ratios. The original continuum-subtracted cube and the continuum-subtracted, Voronoi-binned cube were separately fitted as described below, and the results employed to produce different maps (see Sect. 3).

The emission-line fitting was carried out with MPFIT (Markwardt 2009). We performed two fittings of the emission lines, using one or two Gaussian components per line. To reduce the parameter ranges over which the minimisation is made and avoid being trapped in local minima with unphysical bases, we constrained some of the fit parameters. The line ratios Hα/Hβ and Hγ/Hβ were allowed to vary above a minimum of 2.87 for the former and below a maximum of 0.466 for the latter, which are their theoretical extinction-free values for Case B recombination and a gas temperature of 10⁴ K (Osterbrock & Ferland 2006). The [S II] λ6716/λ6731 ratio was allowed to vary in the range of values 0.46–1.43, outside which it stops being sensitive to electron density (Osterbrock & Ferland 2006). All the emission lines were tied to have the same velocity shift and velocity dispersion. These constraints were applied separately to both the one- and two-component fits. To guide the fit, the starting guess of centroid velocity of the first component, meant to reproduce the bulk of the line, was matched with the wavelength of the [O III] line peak, and the parameter was allowed to vary between [ − 50, +50] km s⁻¹ around that. Based on the spectral shifts of the broad line wings observed in the data, the starting guess of velocity of the second component, meant to reproduce outflowing motions in the form of asymmetric wings or additional line components, was set to +300 km s⁻¹ (–300 km s⁻¹) with respect to the line peak, with the positive (negative) sign in case of a redward (blueward) asymmetry of the line profile, and the parameter was allowed to vary by up to 500 km s⁻¹ with respect to the peak velocity towards the asymmetric side of the line. To assess which of the two fits (with one or two Gaussians) to adopt in each spaxel, depending on the complexity of the line profiles, we employed the improvement of the reduced χ² as a criterion. The two-component fit was selected only in those spaxels whose two-component reduced χ² was smaller than 80% of that resulting from the one-component fit. Moreover, we only selected the two-component fit when the peak S/N of the second component was greater than 3 in the [O III] line, to avoid selecting a second component reproducing the noise. In the nuclear region, the line profiles are very complex and could not be reproduced with two components only (see e.g. spectrum number 1 in Fig. 2, where a highly blueshifted outflowing component requires the inclusion of an additional Gaussian in the fitting). We thus added a third component in a radius of 12 spaxels (∼2.4″) from the nucleus, to account for these additional outflowing components. The same selection criterion for the number of components based on the reduced χ² of the three- versus two-component fits described above was applied, as well as the S/N = 3 cut on the third component of [O III].

Fig. 2. Example of modelled spectra from different regions, covering Hβ and [O III] (upper sub-panels) and [N II], Hα, and [S II] (lower sub-panels). From top to bottom, the spectra come from close to the nucleus, to the NE of it (number 1), from the eastern lobe or Teacup handle (number 2), and from the western lobe (number 3); these regions are labelled accordingly in Fig. 1, right panel. For each emission line, blue, green, and orange curves mark the first (narrower), the second, and the third (broader) Gaussian components employed in the fitting, respectively.

The [S II] doublet was separately modelled from the data cube we produced after having removed the sky residuals affecting its spectral region by running ZAP (Sect. 2.1). We first subtracted the stellar continuum previously obtained for the other data cube, again by re-scaling the model obtained for a given bin to the flux of each spaxel belonging to it. As done before, we obtained a Voronoi-binned cube from the continuum-subtracted one, by using the same binning employed for the original data cube (the one on which we did not apply ZAP). We then modelled the [S II] doublet on both the binned and the unbinned cube. To do so, in each spaxel we used the same number of Gaussian components used to model the other emission lines in the original cubes (binned and unbinned) and fixed the velocity and velocity dispersion of each Gaussian component of [S II] to the values obtained for the fit of other emission lines, leaving only the flux free to vary.

The results of all the fitting procedures were then used to produce different maps for the ionised gas, as described in the following in Sect. 3. The stellar kinematics, that is not among the main goals of this work, is reported and briefly discussed in Appendix A. Some examples of modelled spectra are shown in Fig. 2.

3. MUSE ionised gas maps

In this section we present the maps that we have obtained from our analysis of the reduced MUSE data cube of the Teacup galaxy. We applied a cut of S/N = 3 on the emission-line peak to all the maps reported in the following. In case more lines are involved in a given map (e.g. due to a line ratio), the cut is applied to all the lines involved, and the spatial extension of the map is thus limited by the weakest line. The side of the field of view (FOV) of the MUSE observations corresponds to ∼100 kpc at the luminosity distance of the Teacup, D_L ∼ 404 Mpc (1″ ∼ 1.66 kpc).

3.1. Flux maps

In Fig. 1 we present a false-colour image of the galaxy obtained from the full FOV of MUSE. The continuum (blue) shows the irregular shape of the galaxy stellar content, resulting from its past merging activity. The [O III] and Hα emission, stemming from ionised gas, are reported in green and red, respectively, and the former generally dominates over the latter (they are shown using the same flux scale). Close to the nucleus, at ∼10 kpc on its eastern and western sides, the handle of the Teacup and its counter-fan, respectively, are recognisable. Hereafter, we call these two ionised gas structures handle or loop when referring to the eastern one, and counter-fan the western one, in accordance with previous works (Keel et al. 2012; Harrison et al. 2015; Ramos Almeida et al. 2017), or simply lobes when referring to both of them. We note that the western optical lobe or counter-fan also shows shell-like emission, though on smaller scales and with a more poorly defined shape compared to the eastern handle. The sensitivity and large FOV of MUSE allow for ionised gas filaments to be revealed and mapped to much larger distances on both sides, up to distances of ∼50 kpc per side. As part of this extended filamentary emission, we find a second weaker loop reaching about 15″(∼25 kpc) to the E–NE of the nucleus, in the same direction as the Teacup handle, having a more irregular shape than the handle. This second loop might in principle have the same nature as the handle, that is, a bubble of gas inflated by the action of the AGN jet and outflows, but corresponding to a previous expulsion episode. The large-scale filaments seem to follow an X-shaped geometry centred on the nucleus. This ∼100 kpc size structure corresponds to the giant nebula that was spectroscopically identified in Villar Martín et al. (2018). More recent work found arc-shaped ionised gas emission on even larger scales (just outside the FOV of MUSE), to both NE and SW (more prominent and definite in shape to the NE), for a total size of the nebula of ∼120–130 kpc (Villar Martín et al. 2021; Moiseev & Ikhsanova 2023). This appears to be the continuation of the X-shaped filaments detected in the MUSE observations, and would represent an additional large-scale bubble, for a total of three consecutive bubbles (each corresponding to a past AGN outburst), one on ∼10 kpc scales (the handle), one on ∼25 kpc scales, and one on ∼60 kpc scales. These multiple gas shells at increasing distances resulting from repetitive AGN ejections are very similar to those predicted by the TNG50 cosmological simulation as a result of AGN feedback (Nelson et al. 2019) and seen from, for instance, NGC 1275 at the centre of the Perseus cluster (e.g. Fabian et al. 2003, 2006; Sanders & Fabian 2007).

Figure 3 shows maps of the [O III] emission line flux, for the total modelled line profile (left; same reported in Fig. 1) as well as for the first, narrower (centre) and second, broader Gaussian component (right) used for the line modelling (Sect. 2.2). The first component follows the same morphology as the total profile. The second component is concentrated in the central regions, including the eastern handle and western counter fan, although it broadly follows the morphology of the total line flux and first component.

Fig. 3. Maps of the [O III] (top panels) and Hα (bottom panels) emission line flux distribution in the Teacup galaxy. The fluxes of the total modelled line profile (left panels; same reported in Fig. 1) and of the first, narrower (centre panels) and the outflow (second plus third), broader Gaussian components (right panels) making up the total profile are reported. The colour intensity scale is the same for all the three maps of each line, as reported in the colour bar to the right. The contours in upper-left panel mark the VLA radio emission at 5.12 GHz from Harrison et al. (2015), same as in Fig. 1, right panel.

The VLA radio contours at 5.12 GHz (same as in Harrison et al. 2015) are superimposed on the [O III] map in Fig. 3. This matching shows that the eastern radio bubble peaks at its eastern edge, along the [O III] loop, consistent with the results from Harrison et al. (2015). The MUSE-VLA comparison also clearly shows that there is a misalignment between the axis of the radio emission, oriented in the E–W direction, and that of the optical emission, tilted in the NE–SW direction with respect to the radio one. This can be seen in both the SW counter-fan and the larger scale structure of the gas on both sides.

3.2. Kinematics

The kinematic maps of the ionised gas in the Teacup are presented in Fig. 4, with the [O III] velocity (upper panels) and velocity dispersion maps (middle panels). These are respectively given by the centroid and standard deviation, in case of a single Gaussian component, or by the first- and second-order moments of the velocity of the total modelled line profile, in case of multiple components. The maps in the leftmost of these panels are obtained from the total modelled profile of [O III]. The central and right panels show the above two quantities for the first and the second (plus third, when present in the central region) modelled Gaussian components, respectively. They are zoomed in the central 20″ × 20″, outside which only Gaussian component was used in the line fitting, making the maps for the first component identical to those for the total line profile. We also show the zoomed versions of the [O III] velocity and velocity dispersion for the total line profile (same reported in top-left and mid-left panels for the full FOV) in lower panels, together with the [O III] flux for comparison. We considered a systemic velocity of 24 486 km s⁻¹ for the maps, consistent with Moiseev & Ikhsanova (2023). We do not report the velocity and velocity dispersion maps for Hα since they are virtually identical to those of [O III].

Fig. 4. Ionised gas kinematic maps of the Teacup. [O III] velocity (upper panels) and velocity dispersion (middle panels), obtained as first- and second-order moments of the velocity, respectively, of the total modelled profile of [O III] (left panels), as well as for the first (central panels) and the second plus third Gaussian components (first- and second-order moments of velocity of their summed profile; right panels) employed for the line profile modelling. Central and right panels are zoomed in the central 20″ × 20″ (indicated by the black box in left panels). Bottom panels: zoomed versions (also in the central 20″ × 20″) of the [O III] flux (left; same as in Fig. 3, top-left panel), velocity (centre; same as in top-left panel) and velocity dispersion (right; same as in mid-left panel) maps for the total modelled emission line profile. The contours indicate the same VLA 5.12 GHz and highest-resolution 6.22 GHz radio images from Harrison et al. (2015) reported in Fig. 1, right panel, with the same contour levels.

The velocity maps show, on the scale of the optical and radio bubbles, a kinematic structure with opposite velocity, that is present in both the first (narrower; σ ≲ 150–200 km s⁻¹) and second (broader; σ ≳ 200 km s⁻¹) components, with the first possibly tracing a rotational feature (as modelling of regular circular rotation in Moiseev & Ikhsanova 2023 seems to suggest; see also Villar Martín et al. 2018) while the second the ionised outflow, known to be present in this direction (Mullaney et al. 2013; Villar Martín et al. 2014; Harrison et al. 2015 in [O III]; Ramos Almeida et al. 2017 in NIR hydrogen recombination lines and [Si VI]). In the direction and scales (∼10 kpc) of the ionised and radio bubbles there is no evidence for stellar rotation (see Fig. A.1, left panel) or for the presence of a gas disc (which could be responsible for the gas rotation) from the ionised gas morphology (Fig. 3). Moreover, rotating H₂ (Ramos Almeida et al. 2017) and CO (Ramos Almeida et al. 2022; Audibert et al. 2023) molecular gas is observed in a different direction (N–S, i.e. almost perpendicular) and on much smaller scales (∼2 arcsec, i.e. ∼3 kpc). Based on these aspects, another possibility is that all the bipolar ionised gas, including the narrower component, is being evacuated as part of the bubbles. However, we prefer to adopt a more conservative approach in this work and consider only the broader component(s) as part of the outflow, when we estimate its properties (Sect. 4).

On scales larger than the optical and radio bubbles, there are two sudden discontinuities in velocity, which drops to smaller values (in absolute value) at the outer edge of the bubbles (∼10 kpc) to rise again at ∼20–25 kpc (∼15″) from the nucleus. This is more visible on the eastern (redshifted) side. This large-scale pattern can be appreciated even better in the velocity maps obtained from the stellar continuum-subtracted, Voronoi-binned cube, reported in Fig. B.2. This structure, on the scale of the ∼100 kpc giant bubble, has opposite-sign velocities (receding to the E, approaching to the W) and is detected only in the first component, with small velocity dispersions (≲100 km s⁻¹). It might thus either represent a separate rotating component in addition to that observed on the scale of the bubbles, possibly a leftover of the past galaxy merger event, or the distinct kinematics (either rotating and/or still expanding) of the largest-scale bubble expelled in a past AGN outburst (see Sect. 3.1; Villar Martín et al. 2021; Moiseev & Ikhsanova 2023).

The velocity dispersion map of the total profile (mid-left and bottom-right panels of Fig. 4) shows that this is maximal within a radius of about 10 kpc around the nucleus, with values from ∼200 km s⁻¹ up to ∼500 km s⁻¹ close to the nucleus. Such high values stem from the second plus third components (mid-right panel), though also the first component (mid-central panel) has some residual enhancements (∼150–200 km s⁻¹) in the central region compared to the average velocity dispersion elsewhere. Interestingly, the highest values of the velocity dispersion (≳250 km s⁻¹) in the total profile and second (plus third) component velocity dispersion maps are not isotropically distributed around the nucleus, but mainly appear in the direction perpendicular to the ionised and radio bubbles, rather than in the direction of the bipolar outflow as normally expected in AGN (e.g. Venturi et al. 2018; López-Cobá et al. 2020; Juneau et al. 2022; Kakkad et al. 2023). Specifically, they reach the highest values (∼500 km s⁻¹) at the location of the jet head, the eastern hotspot HR-B in the highest-resolution radio image at 6.22 GHz. We discuss this aspect in more detail in Sect. 5.

Tentative detection of a very broad (FWHM ∼ 3000 km s⁻¹) Paα component to the NE and SW of the nucleus was reported by Ramos Almeida et al. (2017), which we do not detect in the (higher S/N) MUSE data. We briefly discuss this in Appendix D.

3.3. Radio-optical misalignment

Interestingly, the comparison of the lower-resolution radio emission at 5.12 GHz with the [O III] velocity maps of the total profile and first-component (bottom-central and top-central panels in Fig. 4, respectively) shows that the optical kinematic axis is clearly misaligned with respect to the axis of the radio bubbles (by ∼20°, as found by Harrison et al. 2015). The optical-radio misalignment is even more evident in the velocity field than in the flux map (Fig. 4, bottom-left panel and Fig. 3, top-left panel).

On the other hand, the highest-resolution (beam HPBW = 0.37″ × 0.24″) VLA radio image at 6.22 GHz follows more closely the optical flux and kinematic axes (Fig. 4, bottom-central, top-central, and top-right panels) and is thus misaligned with respect to the large-scale radio structure. On scales of ∼2.5″ from the nucleus, however, the velocity map of the first (narrower) component shows a kinematic structure with opposite velocities (redshifted to the NE and blueshifted to the SW, exceeding 300 km s⁻¹) that is misaligned from both the radio bubbles at 5.12 GHz and the highest resolution radio hotspots at 6.22 GHz and whose centre appears to be offset from the nucleus. This might represent an additional fast rotating structure, possibly resulting from the past merging activity. The kinematics of the narrower component is very complex and understanding it would require accurate kinematic modelling, which goes beyond the scope of this work.

As mentioned before, the two unresolved knots constituting the high-resolution radio structure are likely a compact radio jet. The jet is very well aligned with the velocity of the second (plus third) broader components (top-right panel), which trace a bipolar outflow. Being aligned with the ionised outflow, the jet might contribute to the outflow driving process. Since also the [O III] flux axis is aligned with that of the jet, this means that the accretion disc (to which the jet is expected to be perpendicular; e.g. Blandford et al. 2019) and the obscuring toroidal structure expected to surround the nucleus down to ∼pc scales (defining the opening angle of the AGN ionising field on galactic scales; e.g. Bianchi et al. 2012; Ramos Almeida & Ricci 2017) are co-planar. This is indeed not always the case, since a misalignment between the radio jet and the [O III] emission, associated with the AGN ionisation field, is sometimes observed (e.g. Goosmann & Matt 2011; Fischer et al. 2013, 2014; May et al. 2020). Possible causes for the observed misalignment between the large-scale radio emission and the small-scale radio jet, the optical lobes, and optical kinematics are jet precession (e.g. Kharb et al. 2014; Kharb 2018) or jet-bending due to jet-cloud interactions (as e.g. in the case of NGC 1068; Gallimore et al. 1996a,b).

No clear evidence for the presence of a gas cloud or complex of clouds emerges from the MUSE observations on the scale of the radio knots (∼0.5″ − 1″), nor of a sudden bend of the jet trajectory, that could favour the jet bending scenario due to interaction with gas clouds. Unfortunately, the resolution of MUSE observations (FWHM ∼ 1″) does not allow for this possibility to be tested further. The higher resolution HST images show no evidence for dust lanes in the jet trajectory (dust lanes are only observed to the W of the nucleus; Keel et al. 2015) or for an impact on clouds (Harrison et al. 2015) which could bend the jet. We note however that the scales of the jet-cloud interaction in NGC 1068 and of the high-resolution jet in the Teacup are very different, being ∼20 pc for the former case and ∼800 pc for the latter. In any case, the jet bending scenario seems unlikely, given that the large-scale radio bubbles are misaligned with respect to the small-scale jet on both the E and W sides, that would require a bending by the same angle due to a jet-cloud interaction on both sides of the nucleus.

On the other hand, jet precession could in principle be a viable explanation for the misalignment. However, it seems very unlikely that the jet during its precession is currently pointing exactly in the direction of the optical flux and kinematic axes, whose direction is defined on tens of kpc. Therefore, all in all, the origin of the misalignment between the large-scale radio emission and the small-scale radio jet, the optical lobes, and optical kinematics is not clear.

3.4. Extinction, density, and temperature

In Fig. 5 we show the maps for some of the physical properties of the gas, namely the implied extinction A_V (from both Hα/Hβ and Hγ/Hβ), the electron density (from [S II] λ6716/λ6731; Osterbrock & Ferland 2006), and the electron temperature (from [O III] (λ5007+λ4959)/λ4363). We adopted a Calzetti et al. (2000) attenuation law for galactic diffuse ISM (R_V = 3.12) and intrinsic ratios of Hα/Hβ = 2.86 and Hγ/Hβ = 0.466 (for an electron temperature of T_e = 10⁴ K; Osterbrock & Ferland 2006) to evaluate the extinction A_V. For the electron temperature, T_e, we employed the relation from Proxauf et al. (2014), updated from that in Osterbrock & Ferland (2006).

Fig. 5. Maps of physical properties of the ionised gas in the Teacup. Extinction AV from Hα/Hβ (top-left) and Hγ/Hβ (top-right), electron density (from [S II] λ6716/λ6731; bottom-left), and electron temperature (from [O III] (λ5007+λ4959)/λ4363; bottom-right). The line fluxes involved in the line ratios employed for the maps are those of the total modelled emission line profiles (with one, two, or three Gaussians). The maps are zoomed in the central 20″ × 20″.

All these maps have been obtained from the total modelled emission line profiles. We do not report the maps for the single components as the lines involved (Hβ, [S II] doublet, and [O III] λ4363) are limited only to the very central regions in their second component due to low S/N.

As can be seen in Fig. 5, top panels, the extinction is maximal around the nucleus, in an elongated area in the NW–SE direction, where it reaches values above A_V ∼ 1.2. This is in agreement with the r–i colour map of the central 4″ × 4″ shown in Fig. E.1 of Ramos Almeida et al. (2022), and with the presence of strong dust lanes westwards of the nucleus (Keel et al. 2015). On larger scales, lower values of the extinction, around 0.6, are also found in some regions in the direction of the [O III] E loop and W counter-fan. The values of extinction obtained from Hγ/Hβ tend to be higher around the centre compared to those from Hα/Hβ by up to a factor of ∼2. This can possibly be due to a discrepancy between the adopted attenuation law (Calzetti et al. 2000) and the actual one experienced by the gas in the galaxy, or to deviations from pure Case B ionisation. However, a comparative analysis of different dust attenuation laws goes beyond the scope of this work.

The electron density (Fig. 5, bottom-left panel) peaks in a limited region close to the nucleus, where it reaches values ≳500 cm⁻³. It quickly drops away from the nucleus, down to values close to or below the threshold of 10 cm⁻³ (below which the [S II] doublet ratio becomes virtually insensitive to density variations) at 1″ − 2″ from the nucleus. Interestingly, the density rises again at the eastern edge of the [O III] loop, with values > 100 cm⁻³.

Finally, the electron temperature (Fig. 5, bottom-right panel) is fairly constant (neglecting the noisy pixels with high values at the edges), being comprised in the range T_e ∼ 1.3 − 1.4 × 10⁴ K. The temperature T_e = 10⁴ K assumed for the intrinsic values of the line ratios of Hα/Hβ and Hγ/Hβ is thus a fairly good approximation of the actual temperature.

3.5. Gas excitation

In Fig. 6 we report the spatially resolved emission-line ratio diagnostic diagrams (Baldwin et al. 1981; Veilleux & Osterbrock 1987) for the Teacup, used to trace the dominant gas ionisation mechanism throughout and around the galaxy. The diagrams exploit ratios between emission lines close in wavelength, making these ratios virtually unaffected by dust extinction.

Fig. 6. Spatially resolved diagnostic diagrams of gas excitation for the Teacup, obtained from the star-subtracted Voronoi-binned data cube. [O III] λ5007/Hβ versus [N II] λ6584/Hα (top-left), versus [S II] (λ6716+λ6731)/Hα (top-centre), and versus [O I] λ6300/Hα (top-right) diagrams, together with their respective spatial distributions (bottom panels). The fluxes of the total modelled emission line profiles have been employed for the diagrams. The solid curve defines the theoretical upper bound for pure star formation from Kewley et al. (2001), while the dashed one in the [N II]-diagram is the Kauffmann et al. (2003) empirical classification. The dot-dashed line represents the demarcation line between Seyfert galaxies and shocks or LI(N)ERs from Schawinski et al. (2007) in the [N II]-diagram and from Kewley et al. (2006) in the [S II]- and [O I]-diagrams. Regions dominated by AGN ionisation are then marked in red, composite regions in the [N II]-diagram in purple, shock-ionised or LI(N)ER regions in green, and star formation-dominated regions would be marked in blue (when present).

The diagrams show that the gas ionisation from AGN photons dominates over almost the whole system (red colour), from the nucleus to the optical lobes and up to the outskirts of the 100 kpc giant gas nebula. However, shock-like ionisation is present in a region close to the nucleus perpendicular to the [O III] lobes (green colour), in the same direction where enhanced velocity dispersion is observed (≳300 km s⁻¹; see Fig. 4, bottom-right and mid-right panels). This is also observed in the single-component diagrams reported in Fig. 7, though only in the first (narrower) component and not in the second (broader) one. In the [N II]-diagram this region does not fall within the area of the diagram associated with shocks, but in the composite region (purple colour), however, the [N II]/Hα map in Fig. B.1 shows that this region is characterised by higher values of the line ratio than in the direction of the optical lobes, which may be ascribed to a shock contribution to the gas ionisation. [S II]/Hα and [O I]/Hα (also reported in Fig. B.1) are more sensitive than [N II]/Hα to the presence of shocks (see e.g. Allen et al. 2008; Sharp & Bland-Hawthorn 2010), which can explain why there is this discrepancy between the [N II]-diagram and the other diagrams. Moreover, the curves in these diagnostic diagrams are in fact not strict separations between the different ionisation regimes. Shocks with velocities between ∼200–400 km s⁻¹, comparable with the enhanced velocity dispersion observed, can indeed reproduce the observed ratios in the composite region of the [N II]-diagram as well (e.g. Allen et al. 2008; Ho et al. 2014). We further discuss these aspects in Sect. 5.

Fig. 7. Same as in Fig. 6, but separately for narrower (first) and outflow (second plus third) modelled components.

4. Characterisation of the extended ionised outflow

The spatially resolved nature of MUSE observations, together with its wide FOV (1′×1′), allows for the properties of the ionised outflow in the Teacup to be characterised on a wide range of scales (from ∼1 to 100 kpc) as a function of distance from the active nucleus. We consider the second and third, broader (σ ≳ 150 − 200 km s⁻¹) Gaussian components from our line modelling as pertaining to the outflowing gas and we adopt the physical properties found for these two components to calculate the properties of the outflow. As discussed in Sect. 3.2, while all the ionised gas (including the first, narrower component) might be in outflow as part of the bubbles, we adopt a more conservative approach and consider only the broader components to estimate the outflow properties.

The summed flux of the two broader components is used to compute properties like outflow density and extinction, while velocity and velocity dispersion are obtained as the first- and second-order moments of the velocity of the summed profile of the two broader components, respectively. We employ the maps produced from the stellar continuum-subtracted Voronoi-binned cube, in order to have a sufficiently high S/N also in the fainter regions and thus extend this analysis to larger distances from the nucleus.

First, we estimate the mass of gas contained in the outflow by converting the extinction-corrected flux of the outflow components of Hα and [O III], using the two following relations (Carniani et al. 2015; Fiore et al. 2017):

$$\begin{array}{c}\hfill M_{\mathrm{out},\mathrm{H}\alpha }/M_{\odot }=0.6\times {10}^9\left(\frac{L_{\mathrm{H}\alpha }}{{10}^{44}\mathrm{erg}{\mathrm{s}}^{-1}}\right){\left(\frac{n_{\mathrm{e}}}{500{\mathrm{cm}}^{-3}}\right)}^{-1}\end{array}$$(1)

and

$$M_{\mathrm{o}\mathrm{u}\mathrm{t},[\mathrm{O}\text{III}]}/M_{\odot }=0.8\times {10}^8\left({\displaystyle \frac{L_{[\mathrm{O}\text{III}]}}{{10}^{44}{\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{s}}^{-1}}}\right){\left({\displaystyle \frac{n_{\mathrm{e}}}{500{\mathrm{c}\mathrm{m}}^{-3}}}\right)}^{-1}$$(2)

which assume fully ionised gas with electron temperature T_e = 10⁴ K (fairly in agreement with the one estimated from the [O III] diagnostic ratio in Sect. 3.4), ${\langle n_e\rangle }^2/\langle n_e^2\rangle =1$ and, for the latter relation, a solar oxygen abundance.

We stress that to estimate extinction and density in each bin we also adopted the outflow (second plus third) modelled components, not the total line fluxes. We used Hα/Hβ for tracing extinction instead of Hγ/Hβ or a combination of the two because the former gives more reliable ratios between the faint outflow components due to the higher S/N of Hβ compared to Hγ. The maps of extinction and density of the outflow components (reported in the top-left and bottom-left panels of Fig. B.3, respectively) have some problems. Some deviant values are present due to the fact that the wings of Hβ and [S II] can be faint and the fitting procedure may have not been optimal in all cases. Therefore, in order to obtain more robust spatially resolved estimates of extinction and density for the second component, as needed to calculate the mass of ionised gas, we masked the most deviant values (> 3σ from the mean at a given distance from the nucleus). We also masked isolated external bins (usually at the edge of the S/N threshold), for which the validity of the inferred value could not be assessed from neighbouring bin values. Finally, we applied a median filter with a kernel size of 3 spaxels to the Hα/Hβ and [S II] ratio maps, before calculating extinction and density, in order to mitigate high spatial frequency variations below the spatial resolution of the observations.

To get robust estimates of the outflow quantities in the central region, where the strongest outflows are present, as well as to obtain reliable uncertainties on them, we performed Monte Carlo simulations with 10 000 iterations of the three-component emission line fitting (as described in Sect. 2.2) in each spaxel within a radius of 10 spaxel from the nucleus; in each iteration, the flux in each spectral channel was varied by the error on the flux times a random normal distribution having standard deviation equal to 1. In some spaxels, some parameters accumulated at their extremes allowed by the fitting procedure for most of the iterations, especially the [S II] ratio, accumulating at its upper limit (lower limit for the density). In these cases, we loosened the constraints in the Monte Carlo simulations by allowing the spectra to vary by an amount equal to the error on the flux times a random normal distribution having standard deviation equal to 3. We could not run a Monte Carlo simulation in each spaxel of the data cube because the computational time would have been excessive.

The density could not be constrained in a handful of spaxels (seven, for a total of ∼0.8″ in diameter) around the nucleus, even with the Monte Carlo simulations, due to the [S II] doublet line profiles being very broad and blended, making the fit highly degenerate. We therefore assigned a typical density value of the central regions to these degenerate spaxels, to allow for the calculation of outflow mass and derived quantities. This value has been estimated by first calculating the median density in each 1-spaxel wide radial annulus comprised between 3 and 8 spaxels from the nucleus, and then obtaining the 50th-percentile of these annular median densities. We stress however that in these seven nuclear spaxels the density is basically unconstrained and this has a deep impact on the total nuclear mass since they host the highest outflow line fluxes. We return to this aspect later in this section.

The resulting maps of extinction and density of the outflow are reported in top-right and bottom-right panels of Fig. B.3, respectively. Here, we set to 10 cm⁻³ those densities falling below this value, where the [S II] ratio gets insensitive to density. Outflow extinction and electron density range between A_V  ∼  0 − 1.6 mag and n_e  ≲  10 to ∼5000 cm⁻³, respectively. In Appendix C we explore alternative independent methods to estimate the electron density.

In order to get a value of extinction and density in the outer regions, where no values are available, and be able to calculate the gas mass in those regions, we assigned them the median values calculated over the available areas (excluding the central region, where these quantities peak and is therefore not representative of the more external areas). The medians were calculated at the level of the Hα/Hβ and [S II] line ratios, and then converted to their relative physical quantities. This resulted in a median extinction and density of A_V ≃ 0 and n_e ≃ 12 cm⁻³.

We thus employed these polished maps of extinction and density of the outflow components to obtain the mass of ionised gas in each spaxel using Eqs. (1) and (2). We stress that electron densities of ∼10 cm⁻³ are to be considered as upper limits, since this is the threshold below which the [S II] ratio is not sensitive to density anymore. Therefore, mass outflow rates resulting from densities of ∼10 cm⁻³ are to be considered as lower limits. Given that electron densities of the ISM can be as low as ∼1 cm⁻³ (e.g. Osterbrock & Ferland 2006), the difference can be at most by a factor of ∼10.

For the mass outflow rate calculation, we follow the same approach adopted in Venturi et al. (2018; see also e.g. Maiolino et al. 2012; González-Alfonso et al. 2017; Harrison et al. 2018; Revalski et al. 2021 for a thorough discussion on mass outflow rate calculation). The mass outflow rate through a spherical surface of radius r subtended by a solid angle Ω is Ṁ(r) = Ωr²ρ(r)v(r). Assuming that, within each shell with radial thickness ΔR, density and outflow velocity are constant, the radial average of the outflow rate within ΔR is

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{out}}=\frac{M_{\mathrm{out}}{v}_{\mathrm{out}}}{\mathrm{\Delta }R},\end{array}$$(3)

where M_out is the outflow mass (obtained from either Hα or [O III]), v_out the outflow velocity, and ΔR the physical width of each spaxel (0.2″, i.e. ∼0.3 kpc). For the outflow velocity, we adopted a combination of velocity and width of the outflow Gaussian component(s) from the emission-line modelling, in accordance with many other works (e.g. Rupke et al. 2005; Fiore et al. 2017; Fluetsch et al. 2019). The idea behind this is that the observed line profile of the outflowing gas is given by the combination of velocities directed at different angles with respect to the line of sight, and therefore, due to projection effects, only the highest observed velocities represent the actual velocity of the outflowing material, when this is directed in the observer’s line of sight. Specifically, we adopt the definition used for instance in Rupke et al. (2005) and Fluetsch et al. (2019), that is,

$$\begin{array}{c}\hfill v_{\mathrm{out}}=v_2+{\mathit{FWHM}}_2/2\simeq v_2+1.18{\sigma }_2,\end{array}$$(4)

where v₂, FWHM₂, and σ₂ are the centroid velocity, the FWHM, and the velocity dispersion of the second modelled component, respectively. In the central region, where also a third Gaussian component was used, the first- and second-order moments of velocity calculated over the summed line profile of second plus third components (v_2 + 3 and σ_2 + 3, respectively), have to be substituted in the equations above.

The map of the mass outflow rate thus obtained from the outflow (i.e. second, plus third when present) modelled Gaussian component(s) of Hα is reported in Fig. 8, top-left panel. The analogous map obtained from [O III] outflow component(s) is very similar (except for having smaller values overall, as we discuss below), therefore we do not show it. We then produced radial profiles of the mass outflow rate as a function of distance from the nucleus, that are displayed in the top-right panels of Fig. 8. The profiles are sampled in bins of 1 kpc. These radial profiles have been produced by summing, in each radial bin, the single-spaxel mass outflow rates from the maps in left panels, and re-scaling by ΔR/ΔR_radbin, being ΔR_radbin = 1 kpc.

Fig. 8. Spatially resolved properties of the galactic ionised outflow in the Teacup. Top-left: mass outflow rate map, obtained from the outflow (second plus third) modelled component(s) of Hα. Top-right: radial profiles of mass outflow rate (top), kinetic rate (middle), and momentum rate (bottom) as a function of distance from the nucleus, from Hα (red and orange) and [O III] (blue and light blue) outflow components. These quantities were obtained by summing up the single-spaxel values (calculated through Eqs. (3), (5), and (6)) contained in each radial bin, and re-scaling by ΔR/ΔRradbin, being ΔRradbin the radial bin width (=1 kpc) and ΔR the spaxel size. Orange and light blue points mark those external radial bins for which we adopted a single median gas density, since no density measurement was available (see Sect. 4 for details). A different symbol is used for the innermost radial bin, which comprises a few spaxels whose outflow density could not be constrained (see Sect. 4). Lighter error bars are used when the assigned statistical uncertainties were estimated from representative regions instead of on a spaxel basis (see Sect. 4). The vertical dashed lines mark the distance of the handle. On the rightmost y axis, we also report mass-loading factor (top), coupling efficiency (middle), and momentum boost (bottom) of the outflow, for which we consider SFR ∼ 10 M⊙ and Lbol ∼ 5 × 1045 erg s−1 (see text). Mid-left: radial profiles of ionised outflow mass, from the flux of the outflow component(s) of Hα and [O III]. Bottom: radial profiles of outflow velocity (left; vout = v2 + 3 + FWHM2 + 3/2 ≃ 1.18σ2 + 3), combination of velocity, v2 + 3 (centre), and velocity dispersion, σ2 + 3 (right), of second plus third modelled components, for both [O III] and Hα. The reported values are the flux-weighted (by the outflow component(s)) averages in each radial bin. We stress that the reported error bars are not statistical errors, but are a combination of the statistical errors and the standard deviation of the values, to show the extent of variation of the velocities in each radial bin.

We first note that the mass outflow rate obtained from Hα is larger than that obtained from [O III] by about a factor of 3–4. This is not surprising, as also Carniani et al. (2015) and Fiore et al. (2017) report a higher outflow mass from Hβ than [O III], by a factor of ∼2 and ∼3, respectively. Carniani et al. (2015) ascribe it as due to the different volumes from which the two lines are emitted, differing exactly by a factor ∼2 in their case study. The discrepancy is a bit larger in the case of the Teacup. The reason for this discrepancy is unclear, one possibility could be that the radiation field in the Teacup is less energetic, giving a weaker [O III] emission for a given Balmer emission line flux (the latter instead directly relating to the number of ionising photons and not to the hardness of the radiation field; e.g. Osterbrock & Ferland 2006). In the rest of the work, we focus on the outflow quantities obtained from Hα, which is a more robust tracer of the outflow mass compared to [O III], not subject to metallicity variations and to the uncertainties mentioned above. We leave the radial profiles of outflow mass (and derived quantities) obtained from [O III] for comparison.

We stress that the mass outflow rate in the innermost radial bin is mostly unconstrained, due to density being unconstrained in seven nuclear spaxels, as mentioned earlier in this section. As said, the outflow line fluxes in these spaxels are among the highest, therefore, even if they are only very few, their impact on the mass outflow rate of the innermost radial bin is dominant. We took this into account by assigning errors spanning the whole density range (∼1–50 000 cm⁻³) to the nuclear problematic spaxels and highlighted this aspect by using a different symbol for the innermost radial bin in the figure.

The mass outflow rate radial profiles exhibit a decreasing trend with distance, dropping by about three orders of magnitude over 30 kpc, from ∼10² M_⊙ yr⁻¹ to ≲10⁻¹ M_⊙ yr⁻¹. Radially decreasing mass outflow rates are observed in the case of other ionised outflows in AGN (e.g. Venturi et al. 2018; Revalski et al. 2018, 2021) and in principle could be both ascribed to a slowing down of the outflow with distance and to the fact that the AGN ionising flux drops as r⁻².

We note that in the case of the cited studies the mass outflow rate peaks at distances in the range ∼0.1–1 kpc, not on the nucleus. In the case of the Teacup, it is unclear whether the mass outflow rate peaks on the nucleus or at larger distances, since it is highly unconstrained in the innermost radial bin, as discussed above. Moreover, the spatial resolution of the MUSE WFM observations analysed in this work, ∼1″, corresponding to ∼1.7 kpc, is coarser than the scales of the peaks reported in the above works. Models of outflows driven by AGN radiation pressure predict that the mass outflow rate is expected to peak at about 1–2 kpc from the nucleus and then decrease with distance (the latter aspect is observed in this study and in those mentioned above) as the outflow propagates (Costa et al. 2018; Meena et al. 2023). Higher spatial resolution observations would therefore be required to investigate whether the mass outflow rate of the Teacup peaks on the nucleus or not. At larger distances, the mass outflow rate exhibits a small secondary peak at ∼10 kpc from the nucleus, corresponding to the size of the handle of the Teacup, indicated as a vertical dashed line in the radial profiles.

In the remaining top-right panels of Fig. 8 we report the radial profiles of kinetic rate (middle) and momentum rate (bottom) of the ionised outflow. These are obtained as follows:

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{kin}}={\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^2/2\end{array}$$(5)

and

$$\begin{array}{c}\hfill {\dot{p}}_{\mathrm{out}}={\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}\end{array}$$(6)

for kinetic and momentum rate, respectively. They also show a decreasing trend with distance from the nucleus and a secondary peak at ∼10 kpc, at the scale of the handle, as in the case of the mass outflow rate.

By inspecting in Fig. 8 the radial profiles of outflow mass (mid-left), outflow velocity (bottom-left), as well as velocity, v₂₍₊₃₎ (bottom-centre), and velocity dispersion, σ₂₍₊₃₎ (bottom-right), of the second (plus third, when adopted) modelled spectral components, which define the outflow velocity (v_out ≃ v₂₍₊₃₎ + 1.18σ₂₍₊₃₎), it is possible to understand which of these quantities determines the shape of the radial profiles of mass outflow rate, kinetic rate, and momentum rate. We thus see that the outflow mass and the velocity dispersion of the second plus third component(s), rather than their velocity shift, drive the shape of the radial profiles of the above outflow properties.

If we adopt instead the more extreme definition from Rupke & Veilleux (2013), Feruglio et al. (2015), Fiore et al. (2017), and many others, v_out = v₂ + 2σ₂, the resulting mass outflow, kinetic, and momentum rates are a factor of ∼1.7, 5, and 3 higher, respectively, than those obtained with the definition we adopted. Another different approach would have been to employ the centroid velocity of the second (plus third) component(s) for the outflow velocity, without any combination with the velocity dispersion, v_out = v₂₍₊₃₎. The assumption behind this would be that the broadening of the emission line profiles is given by a real distribution of velocities of the outflow and not by projection effects due to the outflow propagating along different angles with respect to the observer’s line of sight. However, considering that the observed v₂₍₊₃₎ are quite small (≲100 km s⁻¹) where the line profiles are broadest, around the centre (up to σ₂₍₊₃₎  ∼  400 − 500 km s⁻¹; see Fig. 4, top-right and mid-right panels, and Fig. 8, bottom-centre and bottom-right panels), we find this alternative assumption unlikely and we have thus not adopted this approach. We point out that, in this case, the values of mass outflow, kinetic, and momentum rate would be a factor of ∼3–20, 10–3000, and 5–200 times lower, respectively, compared to those resulting from the adopted approach, with the largest differences occurring at the smaller radii from the nucleus, where the two approaches give very different outflow velocities (≳500 km s⁻¹ in case of the adopted definition versus ≲100 km s⁻¹ in case of the alternative v_out = v₂₍₊₃₎ assumption; Fig. 8, bottom panels).

4.1. Outflow energetics: Jet versus AGN acceleration

The comparison between the outflow and the radio jet energetics allows us to obtain insights on the outflow driving mechanisms. Harrison et al. (2015) infer a jet kinetic luminosity of Ṗ_jet ∼ 2.5 × 10⁴² erg s⁻¹ from the radio luminosity of the two high-resolution knots HR-A and HR-B, assuming that 1% of the total jet energy is converted into radio luminosity, while Audibert et al. (2023) find Ṗ_jet ∼ 1 − 3 × 10⁴³ erg s⁻¹ from the empirical relations of Bîrzan et al. (2008) and Cavagnolo et al. (2010). Therefore, we adopt a range of jet powers, Ṗ_jet ∼ 10^42.5−43.5 erg s⁻¹ for our comparisons.

To check whether the jet is capable of accelerating the ionised outflow, its power should be compared with the outflow kinetic rate on similar scales (∼1 kpc), which correspond to those of the innermost radial bin (Fig. 8). However, this is mostly unconstrained since it is affected by the degeneracy in estimating the density in a few nuclear spaxels. We therefore compare the jet power with the kinetic rate in the second innermost radial bin. This ranges between 5 × 10⁴² − 2 × 10⁴³ erg s⁻¹ considering the uncertainties. This would require a transfer of ∼15% up to 100% of the jet kinetic power to the ISM. Transfer efficiencies of the jet energy to the kinetic energy of the ISM of ≲30% are expected from simulations (e.g. Wagner & Bicknell 2011; Mukherjee et al. 2016), which only partially overlap with the above range. Therefore, purely from the energetic point of view, the jet alone is able to drive the outflow only when simultaneously considering the lower end of the inferred outflow kinetic rates and the upper end of the jet powers.

We can also investigate if the AGN alone is able to accelerate the outflow. Different outflow acceleration modes are predicted by theory. In the thermal feedback scenario, a fast wind accelerated by the AGN on accretion-disc scales drives the expansion of a hot, shocked bubble in the galaxy ISM (e.g. King 2003; King & Pounds 2015; Costa et al. 2014). Depending on the cooling efficiency, the bubble may either radiate away most of its energy and push the ambient ISM through its ram pressure (“momentum-conserving” outflow) or retain its energy while it expands, efficiently sweeping up the ISM it encounters (“energy-conserving” outflow). In the radiation pressure-driven scenario, the AGN is able to directly accelerate gas located on galactic scales through radiation pressure on dust, providing large momentum boosts of the outflow, ṗ/(L_AGN/c) (e.g. Thompson et al. 2015; Ishibashi & Fabian 2015; Ishibashi et al. 2018; Costa et al. 2018).

Therefore, we compare the above ionised outflow kinetic rate of the Teacup with its AGN bolometric luminosity, L_AGN. We estimate L_AGN by using the 2–10 keV X-ray luminosity from Lansbury et al. (2018; L_X ∼ 0.8 − 1.4 × 10⁴⁴ erg s⁻¹) and the X-ray dependant bolometric correction k_bol(L_X) from Duras et al. (2020) for type-2 AGN, including its intrinsic spread of 0.27 dex. We thus obtain L_AGN ∼ 7.2  ×  10⁴⁴ − 5.6  ×  10⁴⁵ erg s⁻¹. We adopt an approximate mean value of L_AGN ∼ 2  ×  10⁴⁵ erg s⁻¹, which matches the estimate from Harrison et al. (2014) from mid-to-far-IR SED fitting. We also adopt this value for the kinetic coupling efficiency and momentum boost radial profiles in Fig. 8 (rightmost y axes). We thus obtain a kinetic coupling efficiency of Ė_kin/L_AGN ∼ 0.003–0.009, where, as said above, we have considered the range of kinetic rates in the second innermost radial bin defined by the uncertainties. When considering the lower and upper estimates of the AGN bolometric luminosity, Ė_kin/L_AGN varies from ∼0.007–0.03 (for L_AGN ∼ 7.2 × 10⁴⁴ erg s⁻¹) to Ė_kin/L_AGN ∼ 0.0009 − 0.003 (for L_AGN ∼ 5.6 × 10⁴⁵ erg s⁻¹).

Despite the large uncertainty on L_AGN which prevents firm conclusions to be drawn, we can see that these ranges of kinetic coupling efficiencies are only partially consistent with the theoretical predictions for an energy-conserving outflow, ∼0.005–0.05 (e.g. Hopkins & Elvis 2010; Zubovas & King 2012; Costa et al. 2014; King & Pounds 2015), while they are more generally compatible with a radiation-pressure driven one, for which kinetic coupling efficiencies of ∼0.001–0.01 are expected (e.g. Ishibashi et al. 2018; Costa et al. 2018).

The momentum boost, considering again the ionised outflow momentum rate in the second innermost radial bin, which ranges between ṗ ∼ 2 − 6 × 10³⁵ g cm s⁻² considering the uncertainties (Fig. 8), is ṗ/(L_AGN/c) ∼ 2–8 (for L_AGN  ∼  2 × 10⁴⁵ erg s⁻¹), ∼7 − 20 (∼7.2 × 10⁴⁴ erg s⁻¹), and ∼0.9 − 3 (∼5.6 × 10⁴⁵ erg s⁻¹). These are generally too low compared to expectations for energy-conserving outflows (∼20) and too high for a momentum-conserving one (∼1). It is instead generally consistent with the range of momentum boosts of ∼1–10 predicted by models of outflows driven by direct radiation pressure on dust within scales of a few kpc (the first ∼5 Myr of propagation; Ishibashi et al. 2018; Costa et al. 2018), which are compatible with the innermost kpc scales we are considering here.

All in all, the large uncertainties involved in the above comparisons prevent firm conclusions to be drawn on the outflow driving mechanism. It seems however that direct AGN radiation pressure on dust is more generally compatible with the observations compared to a thermal acceleration mechanism, either energy- or momentum-conserving. From purely energetic arguments, the jet alone is capable of driving the outflow only when considering the lowermost kinetic rates and the uppermost jet powers. Another possibility is that the ionised outflow is driven by the combined action of the jet and the AGN, through either direct radiation pressure or the action of an energy- or momentum-conserving bubble.

We note that the above estimates of the ionised outflow properties for the Teacup are around the high end of those found for ionised outflows in AGN of similar luminosity (L_AGN ≲ 5 × 10⁴⁵ erg s⁻¹), based on current observational literature (Carniani et al. 2015; Fiore et al. 2017; Bischetti et al. 2019). This may be due to the fact that these literature estimates are based on integrated spectra or barely resolved observations, for which an outflow mass integrated over several kpc and a single outflow radius are employed in Eq. (3). Such a radius may correspond to the largest distance reached by the outflow, when available, or by the ionised gas emission from imaging, or to the entire spectral aperture. This may consequently lead to a lower mass outflow rate than in the spatially resolved case, when considering the inner regions from the nucleus (corresponding to a small outflow radius), if the outflow mass peaks there as in the case of the Teacup.

We then attempted to obtain an integrated mass outflow rate by modelling the spectrum extracted from a circular aperture centred on the nucleus. However, due to the mixing of multiple velocity components stemming from different regions resulting in very broad, blended [S II] line profiles, the modelling is degenerate and cannot constrain the outflow density and thus its mass, as in the case of the handful of nuclear spaxels discussed in Sect. 4. This issue persisted with different aperture sizes. Hence, we adopted a different approach and obtained a pseudo-integrated mass outflow rate by summing up the single-spaxel mass outflow rates within a circular aperture re-scaling the result by ΔR/R_ap, with ΔR the width of each spaxel and R_ap the radius of the aperture. By using apertures with radii of 1.5 (SDSS-like), 5, 10, and 20 arcsec (corresponding to about 2.5, 8.3, 16.6, and 33.2 kpc, respectively), we obtained Ṁ_out ∼ 43, 21, 11, and 6 M_⊙ yr⁻¹, respectively. The value obtained from the SDSS-like aperture is consistent with those measured on similar scales for the second and third innermost radial bins in the mass outflow rate radial profile (Fig. 8, top-right panel) and the aperture size is compatible with the AGN outflow radii reported in the literature (Carniani et al. 2015; Fiore et al. 2017; Fluetsch et al. 2019, 2021). In light of this, we adopt the mass outflow rate in the second innermost radial bin to carry out quantitative comparisons since it is more robust than the SDSS-like pseudo-integrated one, which includes the nuclear spaxels where the outflow electron density, and thus the outflow mass and mass outflow rate, could not be constrained.

The higher ionised mass outflow rates of the Teacup compared to those observed in other AGN might be related to the fact that it was originally selected from SDSS images precisely because of its spatially resolved line emission, combining high surface brightness with large angular extent (Keel et al. 2012). This selection bias may have enhanced the probability of detecting a more powerful and extended ionised outflow.

We also note that the momentum boosts we find for the ionised outflow in the Teacup are much lower than those reported for the same object and gas phase by Harrison et al. (2014) of ṗ/(L_AGN/c) ∼ 40. However, in their work, the mass outflow rate and related quantities are the result of the mean between two different methods, one employing the same approach used in this work (i.e. a relation similar to Eq. (3)), the other considering an energy conserving bubble in a uniform medium (e.g. Heckman et al. 1990; Nesvadba et al. 2006; Veilleux et al. 2005), which gives values larger by a factor of ≳100–1000 compared to the former method.

4.2. Multi-phase outflow

The recent study of the molecular phase of the outflow in Ramos Almeida et al. (2022) allows for a comparison with the ionised outflow component presented in this work. They detect a molecular outflow in the inner ∼0.5″(∼800 pc) from the nucleus, with maximum velocities of ∼250 and −180 km s⁻¹ for the approaching and the receding parts of the outflow, respectively. These are much lower (by a factor ≳2–3) than the velocities up to more than 600 km s⁻¹ found for the ionised component. A lower outflow velocity for the molecular phase is not surprising (see e.g. Carniani et al. 2015; Fluetsch et al. 2019) and can be ascribed to the stronger resistance to acceleration by the dense molecular gas. Ramos Almeida et al. (2022) report a de-projected mass outflow rate of Ṁ_out = 15.8 M_⊙ yr⁻¹ for the cold molecular phase from ALMA CO observations on the same spatial scales. In a more recent work, Audibert et al. (2023) considered four different outflow scenarios to calculate the molecular outflow mass of the Teacup, using the same CO(2–1) dataset as Ramos Almeida et al. (2022), and they reported mass outflow rates ranging between 15 and 41 M_⊙ yr⁻¹.

The mass outflow rate of the ionised gas that we obtain in the second innermost radial bin from Fig. 8 (being the innermost one basically unconstrained; see Sect. 4), ranging between ∼40–130 M_⊙ yr⁻¹ considering the uncertainties, is larger than that of the molecular gas, by a factor of ∼3–8 or ∼1–3 when considering the lowest (15 M_⊙ yr⁻¹) or the uppermost (41 M_⊙ yr⁻¹) values for the molecular phase, respectively. Therefore, the addition of the molecular budget to the ionised one does not substantially change the picture that we have obtained from the ionised phase alone, in terms of outflow acceleration mechanisms and impact on the host galaxy.

The discrepancy found in the Teacup between the two gas phases is surprising, given that the molecular phase is usually found to dominate the mass outflow rate budget, exceeding the ionised outflow one by a factor of ∼10–10³ at L_AGN ≲ 10⁴⁶ erg s⁻¹ (Carniani et al. 2015; Fiore et al. 2017; Bischetti et al. 2019; Fluetsch et al. 2019, 2021). We stress, however, that the limited number of molecular outflows so far quoted in literature may be biased to the most extreme values, because they are the easiest to detect in CO-bright galaxies. We point out that, in the case of the Teacup, the ionised and molecular outflow properties are obtained with different methods, since Ramos Almeida et al. (2022) extract the outflowing flux by integrating only the high-velocity gas in a slit of 0.2″ oriented along the CO disc minor axis (E–W direction), while we employ the flux of the second, broad component of the fit (calculated all over the regions at the same distance from the nucleus, not only in a 0.2″-wide slit). This might explain the large discrepancy in favour of the ionised outflow. Moreover, the molecular outflow resides in the inner ∼0.5″ (∼0.8 kpc), while the ionised mass outflow rate of the second innermost radial bin we have considered for the comparison resides on scales of ∼2 kpc. Higher-resolution observations for the ionised gas are therefore needed for a matching with the molecular outflow on the same spatial scales.

We can then evaluate the effect of the outflow on the star formation activity and molecular gas reservoir in the Teacup. Considering the SFR of the Teacup (between ∼8–12 M_⊙ yr⁻¹, from FIR emission; Jarvis et al. 2020; Ramos Almeida et al. 2022), the mass-loading factor (η = Ṁ_out/SFR) of the ionised gas is ∼3–30 (∼5–30 when also considering the molecular outflow) at the second innermost radial bin. This means that the outflow consumes the gas reservoir much faster than star formation processes. Mass-loading factors above 1, and even exceeding 10, are often found in AGN (e.g. Cicone et al. 2014; Fiore et al. 2017; Fluetsch et al. 2019).

We can calculate the molecular gas depletion time due to outflows, that is, the time needed to remove from the galaxy all the molecular gas reservoir available for star formation given the current mass outflow rate (assuming that the cold gas content is not replenished with fresh supply), as τ_depl = M_gas/Ṁ_out, where Ṁ_out is the mass outflow rate of the multi-phase (ionised plus molecular) outflow. Here we consider M_gas ∼ M_H2, for which we employ the total molecular gas mass of M_H2 ∼ 6 × 10⁹ M_⊙ from Jarvis et al. (2020) and Ramos Almeida et al. (2022). This gives τ_depl ∼ 4 × 10⁷–10⁸ yr, consistent with the values found for AGN of similar luminosity (L_AGN  ∼  10^{45 − 46} erg s⁻¹; Cicone et al. 2014; Fluetsch et al. 2019), while the depletion time due to star formation (M_gas/SFR) is larger, ∼5 − 8 × 10⁸ yr (consistent with the values expected for galaxies above the star-forming main sequence at the redshift of the Teacup; e.g. Schinnerer et al. 2016; Tacchella et al. 2016; Tacconi et al. 2018). Such short depletion timescales indicate that the outflow may be effective in depleting the molecular gas reservoir. We stress, however, that these are the actual mass-loading factor and depletion timescale only if the outflow continues at the present rate (actually, at the rate given by the second innermost bin considered, at a scale of ∼2 kpc).

The ionised gas emission observed on large scales actually indicates that the AGN in the Teacup was more luminous in the past (∼10⁵ yr timescales) and has faded since then, by a factor of ∼10 considering a present-day bolometric luminosity of L_AGN ∼ 2 × 10⁴⁵ erg s⁻¹ (see Gagne et al. 2014; Keel et al. 2017; Villar Martín et al. 2018 and discussion in Lansbury et al. 2018). Therefore, the nuclear outflow would have also likely been higher in the past than today and consequently the past gas depletion times may have been even lower than that estimated above from present-day nuclear outflow rate, suggesting that the outflow in the Teacup could have been even more effective than today in depleting the molecular gas reservoir in the galaxy.

To estimate to which extent the outflow has actually affected the host galaxy, we compare the mass of gas participating in the outflow with the molecular gas reservoir available from star formation. We obtain the total ionised mass in outflow by summing together the outflow mass values found from our spatially resolved analysis. This gives an ionised outflow mass of ∼3.6 × 10⁸ M_⊙, or a total outflow mass of ∼4 × 10⁸ M_⊙ when also considering the molecular mass (∼3 × 10⁷ M_⊙; Ramos Almeida et al. 2022). The outflow mass is then about 5% the total molecular gas mass in the galaxy, ∼ 6×10⁹ M_⊙.

We can evaluate what fraction of the outflow is actually able to escape the galaxy gravitational potential. We calculate the escape velocity by considering a Hernquist potential (Hernquist 1990) for the stellar component and a Navarro–Frenk–White (NFW) potential (Navarro et al. 1997) for the dark matter (DM) halo, given by the following density profiles, respectively:

$$\begin{array}{c}\hfill \rho (r)=\frac{M_{\star }}{2\pi }\frac{a}{r}\frac{a}{{(r+a)}^3}\end{array}$$(7)

and

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_{\mathrm{crit}}{\delta }_{\mathrm{c}}}{(r/r_{\mathrm{s}}){(1+r/r_{\mathrm{s}})}^2}·\end{array}$$(8)

In the former, M_⋆ is the stellar mass of the Teacup, for which we adopt log (M_⋆/M_⊙)≃11.15 (Ramos Almeida et al. 2022), and r_eff ∼ 3.7 kpc its effective (half-light) radius from the SDSS g-band de Vaucouleurs profile fit (Abdurro’uf 2022), where r_eff ≃ 1.8135a. In the latter, ρ_crit is the critical density of the Universe, δ_c the characteristic overdensity for the halo, and r_s = r₂₀₀/c the characteristic radius, c being the concentration parameter. The virial DM halo mass (M₂₀₀) is inferred from the stellar mass via the stellar-to-halo mass relation (SHMR) of Moster et al. (2013), inverted following Posti et al. (2019), and the concentration c is obtained through the M₂₀₀ − c relation of Dutton & Macciò (2014). We thus obtain log (M₂₀₀/M_⊙)∼13.2. We calculate the escape velocity for each spaxel in the FOV, each having a given distance r from the centre of the galaxy. This results to be in the range ∼1200–1500 km s⁻¹ for the radial distances considered (< 30 kpc).

We then obtain the portion of flux of the outflow component (the second and third we have used in the line fitting; see Sect. 4) having a velocity, either blue- or red-shifted, larger than the escape velocity, |v| > v_esc. From the sum over all the spaxels of the fluxes of this high-velocity part of the outflow component (F_{|v|  >  vesc}) and of the fluxes of the entire outflow component (F_out), we derive the outflow escape fraction as F_{|v|  >  vesc}/F_out ∼ 0.4%. This is the fraction of the outflowing material that is fast enough to be able to escape the galaxy gravitational potential. Considering that the outflow mass is about 5% the molecular gas mass in the galaxy (see above), this means that the outflowing gas that is actually able to escape the DM halo amounts to about 0.02 % of the molecular gas mass.

We stress however that the actual escape fraction may be higher, since the outflow velocities measured from the line profiles are subject to projection effects and the fraction of outflowing gas having intrinsic velocities higher than the escape velocity may then be larger. Moreover, the above calculations of the escape velocity only hold for ballistic motions. If instead the outflow keeps on being accelerated (as it is in the case of the driving mechanisms probed in Sect. 4.1), a larger fraction of it would be able to escape the galaxy and DM halo potential. On the other hand, the ejected gas could encounter mechanical resistance opposed by the gaseous halo, which would have the effect of slowing down the outflow. It is not obvious which of these effects would dominate and we refrain from further speculations.

Summarising, while on one hand the ionised outflow is observed at distances of kpc to tens of kpc, meaning that it is actually able to escape the galaxy, and its depletion times are quite small (≲10⁸ yr), on the other hand only a negligible fraction of the galaxy gas reservoir appears to be actually able to escape the DM halo. The gas that does not leave the DM halo potential may possibly be re-cycled and re-accreted on the galaxy at later times, as invoked by theory (e.g. Tumlinson et al. 2017 and references therein). Even if the gas is unable to escape the DM halo potential, the outflow could still have a significant feedback effect without this being completely ejective. The material expelled from the galaxy could progressively inject energy into the gaseous halo which could hinder, over time, the gas cooling and its (re-)accretion on the galaxy (e.g. Costa et al. 2020).

5. Enhanced line velocity widths perpendicular to the radio structure

The ionised gas velocity dispersion presents a strong enhancement in an elongated region to the NW and SE of the nucleus, along the optical minor axis, perpendicular to the [O III] lobes and radio jet, with values ≳300 km s⁻¹, visible in both the [O III] full-profile and second-component maps (Fig. 4, bottom-right and mid-right panels, respectively). The enhancement region extends over ∼5″, corresponding to ∼8 kpc. Such line width enhancement can also be seen in the VIMOS IFU data in Harrison et al. (2015) and is also present in the molecular gas from CO observations (Ramos Almeida et al. 2022; Audibert et al. 2023), though in the molecular phase it only extends on ∼1.6 kpc and has lower velocity dispersions (∼100 − 140 km s⁻¹) than in the ionised phase.

Regarding the dominant ionisation mechanism, the gas in this velocity dispersion-enhanced region has line ratios compatible with the presence of shocks according to the emission-line ratio diagnostic diagrams (Figs. 6 and 7), especially when considering only the first, narrower component (upper panels of Fig. 7).

We can take advantage of the findings from recent works for interpreting this phenomenon, that is, a velocity dispersion enhancement elongated perpendicular to the radio jets and the AGN ionisation cones observed in radio-quiet local AGN (e.g. Riffel et al. 2014; Shin et al. 2019; Feruglio et al. 2020; Cazzoli et al. 2022; Venturi et al. 2021, and references therein). Venturi et al. (2021), based on the fact that all the sources in which the phenomenon is observed so far host a low-power (≲10⁴⁴ erg s⁻¹), compact (≲1 kpc) AGN jet having low inclinations with respect to their host galaxy discs, concluded that the most likely origin for the observed phenomenon may be the action of the low-power jets interacting with the ISM in the galaxy disc. Supporting this interpretation, simulation of jet-ISM interaction (e.g. Wagner & Bicknell 2011; Mukherjee et al. 2016, 2018a,b) predict that low-power (≲10⁴⁴ erg s⁻¹) jets having low inclinations with respect to the galaxy disc (≲45°) will have maximal interaction with the disc material and, during their slow propagation through the disc, will perturb and shock the ISM and drive turbulence perpendicular to their direction of motion and to the disc plane, along the direction of minor resistance.

In the case of the Teacup, the same phenomenon may be at work, operated by the kiloparsec-scale jet hosted in the galaxy centre. The jet is found to subtend a small inclination angle relative to the molecular gas disc (∼32°), leading to maximal jet-ISM interaction, and the connection between the jet and the perpendicular molecular gas with large velocity dispersions is supported by tailored simulations (Audibert et al. 2023). From the ionised gas velocity dispersion maps in Fig. 4, bottom-right and mid-right panels, we see that the maximum values (∼500 km s⁻¹) are co-spatial with the head of the jet (the knot HR-B). This strongly points to a causal link between the jet and the large velocity dispersions perpendicular to it.

The detection of this phenomenon in the Teacup indicates that it occurs not only in Seyferts, but also in more radiatively powerful (L_AGN ≳ 10⁴⁵ erg s⁻¹) QSOs, as also reported in Girdhar et al. (2022) and Ulivi et al. (2023) in QSOs of similar redshift and AGN luminosity as the Teacup. Moreover, if the large-scale radio lobes are the result of past nuclear activity and the small-scale radio jet consequently represent a re-triggering of jet activity, it follows that the perpendicular velocity dispersion enhancement produced by a compact radio jet in its early phases may occur multiple times during an AGN lifetime, as a result of re-triggering of jet activity. This would help to mix the ISM and, if the turbulent material is lifted perpendicular to the galaxy disc, even the circumgalactic medium, therefore chemically enriching it.

The fact that the molecular gas is less affected by this phenomenon compared to the ionised gas in terms of magnitude and spatial extension of the velocity dispersion enhancement, as also observed in the case of other sources exhibiting the same phenomenon (Diniz et al. 2015; Feruglio et al. 2020; Riffel et al. 2015; Shimizu et al. 2019; Girdhar et al. 2022), may be due to the higher density and clumpiness of the molecular gas, making it harder to perturb and disperse it. On the other hand, by comparing the CO(2–1) and CO(3–2) emission, Audibert et al. (2023) also reported enhanced gas temperature in the direction perpendicular to the jet in the Teacup, as also predicted by the simulations.

6. Stars in the Teacup handle

In Fig. 9, we report an image of the stellar continuum emission, in the (observed) spectral range 5500–6700 Å, free from the strongest gas emission lines. This shows an excess co-spatial with the prominent eastern ionised gas loop, or Teacup handle. However, the handle is a place of intense gas ionisation where even usually faint emission lines can be quite intense with respect to the continuum. Therefore, the continuum emission reported in Fig. 9, despite the exclusion of the strongest emission lines (such as [O III], Hβ, [O I], and Hα), may be contaminated by other emission lines, that could be responsible for the continuum excess in the handle.

Fig. 9. Map of continuum emission, obtained by averaging the (observed) spectral channels of the data cube between 5500–6700 Å, free from the strongest gas emission lines (same as in Fig. 1, left panel, here zoomed in the central 40″ × 40″).

For this reason, we performed a refined extraction of the continuum by carefully selecting different spectral intervals which are completely free even of very faint emission lines in the handle, that we checked visually. To further minimise possible contributions from faint emission lines, the continuum is calculated as the median, instead of the mean, of the signal in the selected spectral ranges. The continuum images obtained in this way are shown in Fig. 10, together with the spectral ranges selected for continuum extraction (grey shaded areas). The spectra, extracted from the aperture labelled as number 1 in Fig. 11, top-left panel, show that no faint line emission is present in the selected spectral ranges at the noise level. Moreover, since the spectra are obtained by collapsing an extended region of 362 spaxels and show no significant presence of weak lines, we can rule out that they contribute on a spaxel level at a much lower S/N. The continuum images from all the six different extraction spectral ranges still show the presence of an excess of continuum emission stemming from the Teacup handle. A higher signal-to-noise image, produced by averaging together the continuum from all the six spectral ranges, is shown in Fig. 11, top-left panel.

Fig. 10. Zoomed-in maps of median stellar continuum in selected spectral ranges free of gas emission lines. The spectral range in which the median has been calculated are indicated by the light grey shaded area in the spectra below. The spectra are extracted by summing all the spaxels in the handle, as traced by [O III] emission, whose contours are reported in top-left image for comparison; the extraction aperture of the spectra is shown in Fig. 11, top-left panel, labelled as number 1.

The continuum excess indicates that stars are actually present in the handle, in addition to the underlying large-scale weak stellar continuum of the galaxy. The presence of stars is confirmed by the sharp stellar absorption lines observed in the spectrum extracted from the handle shown in Fig. 11, bottom panel, whose y-axis is stretched to highlight the absorption features.

Fig. 11. Properties of the stellar populations, showing evidence for young stars in the Teacup handle. Top-left: zoomed-in map of the stellar continuum obtained by performing the median in all the six spectral ranges reported in Fig. 10 together (light grey shaded areas). Mid-left: map of continuum flux ratio (i.e. colour) between the images corresponding to the bluest and the reddest spectral ranges reported in Fig. 10 (top-left and bottom-right panels, respectively). Prior to performing the ratio, the continuum images were smoothed with a Gaussian kernel having σ = 1.5 spaxels, for a better visual output. The patch having low continuum flux ratio in the southern part of the handle is due to a background red galaxy. Top-right: luminosity-weighted fractions of stellar populations in the stellar metallicity (in log scale, with respect to solar metallicity) versus age plane resulting from the stellar population modelling (details in text) of each numbered aperture shown in top-left panel. Bottom: spectrum (black) extracted from the aperture covering the handle (number 1 in top-left panel) with the modelled stellar continuum overlaid (red); shaded grey vertical bands indicate spectral ranges which were masked in the fitting, and the spectrum within these ranges is marked in blue.

Two are the possible scenarios explaining the presence of stars in the handle. The first is that the handle is a tidal feature due to the past merger occurred in this system, in which both stars and gas are present, the gas being illuminated by the AGN. The second, in which the handle is a bubble of gas inflated by the action of radio jet and outflowing gas, as also suggested by the handle being in the same direction of jet and outflow, by its co-spatiality with the eastern radio bubble (see Sect. 3 and Harrison et al. 2015), and by having its vertex on the nucleus (as expected for AGN-inflated bubbles from simulations, see e.g. Nelson et al. 2019), is that the stars have formed directly in the handle due to the action of the jet and outflow, that through the compression of gas clouds triggered the formation of stars (so-called positive feedback).

We checked the colour of the continuum to investigate if a bluer colour is associated with the handle, that would suggest the presence of a younger stellar population and support the second scenario above. To do so, we performed the ratio between the bluest and the reddest line-free continuum images from Fig. 10. The map of this blue-vs-red continuum ratio is presented in Fig. 11, bottom-left panel. The map shows that the handle is characterised by a higher ratio (i.e a bluer continuum) than the rest of the galaxy. The region within ∼2 − 3″ to the E of the nucleus, constituting the base of the handle, is also characterised by a bluer continuum.

To check for the actual presence of a younger stellar population in the handle than in the rest of the galaxy, we modelled the stellar continuum in the handle in terms of its stellar population and compared with that obtained in other regions of the galaxy (dashed apertures in Fig. 11, top-left panel). To do so, we made use of a python-based graphical user interface (GUI) which allows us to extract spectra from single spaxels or multi-spaxel regions (details in D’Ago et al. 2021). The interface allowed us to interactively probe a variety of regions and to test different fitting parameters to model the spectra with the software PPXF (Cappellari & Emsellem 2004; Cappellari 2017). As done for the stellar continuum fitting of the cube, we employed the E-MILES single-stellar population (SSP) model spectra (details in Sect. 2.2). In this case, since we are interested in the properties of stellar populations (rather than in simply modelling and subtracting the stellar continuum as in Sect. 2.2), we adopt a more realistic Chabrier (2003) IMF and BaSTI isochrones (Pietrinferni et al. 2004, 2006), which provide a finer sampling of the parameters of the stellar populations (age and metallicity) compared to the Padova+00 isochrones. Gas emission lines were masked and the automatic sigma-clipping included in PPXF was also employed to mask additional fainter emission lines. Since our goal in this case is the modelling of the stellar population and not of its kinematics, we did not employ an additive polynomial but only a multiplicative one, to avoid changes in the strength of the absorption line features in the templates, as indicated in Cappellari (2017). We tested different multiplicative polynomial degrees, between 5 and 30. The variation of the degree did not result in a significant change in the dominant stellar population of each region. We thus employed an intermediate value of 15 for the multiplicative polynomial degree. We used a regularisation parameter of 100 to obtain the smoothest solution among the many degenerate ones that equally reproduce the data (details in Cappellari 2017).

The resulting luminosity-weighted fractions of stellar populations in the stellar metallicity versus population age plane are reported in Fig. 11, top-right panels. We see that while other parts of the galaxy (regions number 3, 4, and 5) are dominated by a relatively older stellar population with an age between ≳0.5–1 Gyr, the handle (region number 1) and the central part of the galaxy (region number 2) show a significant presence of a younger stellar population, with ages ≲100–150 Myr, typical of OB stars. In particular, in the handle (as well as in the centre) this younger component by far dominates the stellar populations. Moiseev & Ikhsanova (2023) confirm that a younger stellar population is present in the central regions of the Teacup, on scales comparable with those of the handle (≲5 arcsec).

This result indicates that star formation activity in the last ∼150 Myr gave rise to the population of young stars, which in the case of the handle may have been triggered by the action of jet, outflow, and expanding radio lobe compressing the gas at the edge of the bubble. Star formation triggered by AGN jets and outflows is predicted by models, due to the compression of the gas clouds, either by direct impact or at the edge of an over-pressurised expanding shocked bubble, increasing their density, leading to their collapse and eventually to star formation (e.g. Rees 1989; Begelman & Cioffi 1989; De Young 1989; Daly 1990; Mellema et al. 2002; Silk & Norman 2009; Gaibler et al. 2012; Nayakshin & Zubovas 2012; Dugan et al. 2014; Silk 2013; Zubovas et al. 2013; Nayakshin 2014; Bieri et al. 2016; Fragile et al. 2004, 2017).

The Teacup handle could therefore represent the expanding bow shock of gas and stars formed in-situ predicted by models of AGN jet-triggered star formation (e.g. Begelman & Cioffi 1989; Gaibler et al. 2012; Dugan et al. 2014; Bieri et al. 2016; Fragile et al. 2004, 2017). In this context, the enhancement of gas density observed at the edge of the handle in the direction of the jet, radio lobes, and outflow may favour the scenario of gas in the handle being compressed, invoked as a necessary condition to form stars by the above models. The timescale of ∼100–150 Myr is compatible with the time needed by the bubble to expand to its actual size (∼10 kpc) considering an expansion velocity of ∼100 km s⁻¹.

Observational evidence for jet- or outflow-induced star formation has been found in a number of cases so far, for example within the filaments along the jet of Centaurus A (Mould et al. 2000; Rejkuba et al. 2002; Crockett et al. 2012; Santoro et al. 2016) and of 4C 41.17 (Bicknell et al. 2000), in companions aligned with jets from another galaxy, as in the case of the Minkowski’s Object around NGC 541 (Croft et al. 2006; Salomé et al. 2015) and the companions of PKS2250-41 (Inskip et al. 2008) and 3C 285 (Salomé et al. 2015), at the jet termination in NGC 5643 (Cresci et al. 2015), in the gas filaments surrounding the giant radio lobes in Coma A (Capetti et al. 2022), and at the edge of the outflow in the type-2 QSO Mrk 34 (Bessiere & Ramos Almeida 2022). See also the reviews in Miley & De Breuck (2008), Wagner et al. (2016), and O’Dea & Saikia (2021).

To our knowledge, the Teacup represents the first case that evidence for widespread stars formed around a galactic bubble inflated by the action of an AGN is found, confirming the predictions from theory.

7. Conclusions

In this work we have presented a spatially resolved study of the ionised gas in the Teacup QSO (z ≃ 0.08506; L_AGN ∼ 2 × 10⁴⁵ erg s⁻¹) employing optical integral field observations from MUSE at VLT. This complex galaxy constitutes a great laboratory to study AGN feedback in all of its manifestations, since it hosts ∼10 kpc-scale co-spatial optical ionised-gas and radio bubbles, a compact (∼1 kpc), low-power (Ṗ_jet ∼ 10^42.5−43.5 erg s⁻¹) jet, multi-phase (ionised plus molecular) outflow activity, and a giant ionised nebula (≳100 kpc size in total; Fig. 1).

We produced maps of ionised gas line-emission, kinematics, excitation, and physical properties (density, temperature, and dust extinction) spanning a total physical size of about 100 kpc. We confirm the presence of an ionised outflow aligned with the optical AGN-ionised lobes (one of which resembles the handle of a teacup, giving the name to the galaxy), with velocities up to ≳500 km s⁻¹ (Figs. 2, 4, and 8). We measured dust extinction from the Balmer decrement, Hα/Hβ, reaching values above A_V ∼ 1.2 around the nucleus, where it peaks, and with values of ∼0.6 or lower in correspondence with the optical emission-line lobes (Fig. 5). We found electron densities reaching ≳500 cm⁻³ at the nucleus and above 100 cm⁻³ at the edge of the Teacup handle, while in the rest of the galaxy it is generally a few tens of cm⁻³ or around or below the threshold value of 10 cm⁻³ under which the [S II] doublet diagnostic is not sensitive to electron density anymore (Fig. 5). We found electron temperatures around 1.3 − 1.4 × 10⁴ K across the galaxy (Fig. 5). The gas ionisation is dominated by the AGN photons (Figs. 6 and 7).

We detected strongly enhanced ionised gas velocity dispersions (≳300 km s⁻¹) elongated over ∼17 kpc perpendicular to the kiloparsec-size radio jet, AGN ionisation lobes, and fast outflows (Fig. 4). The same is observed in the cold molecular phase, though with smaller velocity dispersions (by a factor of ∼4–5) and spatial extension (by a factor of ∼5). The largest ionised gas velocity dispersion (∼500 km s⁻¹) is observed co-spatial with the head of the jet, strongly pointing to a causal connection between the two. This phenomenon is similar to what was found in nearby AGN (Seyferts) also hosting compact, low-power jets, interpreted as due to the action of the jets interacting with the galaxy ISM and driving turbulent material perpendicular to their propagation direction. If the same process is occurring here, this work demonstrates that the phenomenon of line velocity width enhancement in the direction perpendicular to jets can also occur in sources that are more radiatively powerful (QSOs) and more evolved from the radio point of view, such as the Teacup (whose radio bubbles reach ∼10 kpc per side).

Three ionised gas bubbles at increasing distances from the nucleus can be identified in the Teacup in the NE direction (with less definite counterparts to the SW), each likely corresponding to a past AGN outburst, at ∼10 kpc (the handle), ∼25 kpc, and ∼60 kpc, respectively (Figs. 1 and 3). These multiple gas shells resulting from repetitive AGN ejections are very similar to what the most recent simulations of AGN feedback predict.

The ∼10 kpc-scale radio bubbles are misaligned by ∼20° with respect to the small-scale (∼1 kpc) jet and the optical kinematic axes, which are instead aligned (Figs. 1, 3, and 4). The origin of the misalignment, whether due to jet precession or jet-bending due to jet-cloud interactions for instance, is unclear.

We determined the properties of the ionised outflow as a function of distance from the AGN (Fig. 8). We observed a decreasing radial trend of the mass outflow rate, kinetic rate, and momentum rate. The mass outflow rate drops from around 40–130 M_⊙ yr⁻¹ in the inner 1–2 kpc to ≲0.1 M_⊙ yr⁻¹ at 30 kpc (i.e. by up to about three orders of magnitude), with a small secondary peak at the distance of the eastern ionised-gas loop (∼10 kpc; the handle). If using a more extreme approach to estimate the outflow velocity adopted in many outflow studies, the mass outflow rate, kinetic rate, and momentum rate would be a factor ∼1.7, 5, and 3 higher, respectively, than those we obtained.

We compared the ionised outflow energetics with the power of the nuclear radio jet (Ṗ_jet ∼ 10^42.5−43.5 erg s⁻¹) and with the AGN radiative output (L_AGN ∼ 7 × 10⁴⁴ − 6 × 10⁴⁵ erg s⁻¹) to assess the dominant outflow driving mechanism. Due to the large uncertainties in the jet and AGN luminosities, as well as in the outflow energetics, we could not robustly constrain the outflow driving mechanism from energetic considerations. Nevertheless, we find that direct AGN radiation pressure on dust is more generally compatible with the kinetic coupling efficiencies (Ė_kin/L_AGN) and momentum boosts (ṗ/(L_AGN/c)) which we estimated for the outflow, compared to an energy- or momentum-conserving outflow scenario due to an expanding shocked bubble driven by an accretion-disc wind. From purely energetic arguments, the jet alone is able to drive the outflow only when simultaneously considering the lower end of the inferred kinetic rates and the upper end of the jet powers. Another, more likely possibility is that the outflow is driven by the combined action of the jet and AGN radiation. The fact that the head of the kiloparsec-scale jet is co-spatial with the largest ionised gas velocity dispersions (Fig. 4) suggests that the jet has a significant role in driving the outflow.

The properties found for the ionised outflow in the Teacup are at the high end for ionised outflows in AGN with a similar luminosity. Moreover, the mass outflow rate in the ionised phase (∼40–130 M_⊙ yr⁻¹) is significantly higher than that in the molecular phase (∼15–40 M_⊙ yr⁻¹; only detected in the inner ∼800 pc from the nucleus) on similar scales. This is in contrast to what is generally observed in other AGN, where the molecular mass outflow rate exceeds the ionised one by a factor of ∼10–10³. The addition of the molecular budget to the ionised one does not substantially change the picture obtained from the ionised phase alone, in terms of the outflow acceleration mechanisms.

The mass-loading factor (Ṁ_out/SFR) of the ionised (plus molecular) gas is ∼3–16 (∼5–20), indicating that the outflow consumes the gas reservoir much faster than star formation processes. The molecular gas depletion time due to the multi-phase (ionised plus molecular) outflows, ${\tau }_{\mathrm{depl}}=M_{{\mathrm{H}}_2}/{\dot{M}}_{\mathrm{out}}\lesssim {10}^8$ yr, is short enough for the outflow to be effective in depleting the molecular gas reservoir. However, we find that the fraction of the ionised outflow that is able to escape the dark matter halo gravitational potential is only ∼0.4%. Therefore, while the outflow depletion time is short and the expelled gas is located at distances of tens of kiloparsecs, meaning that the outflow can escape the galaxy, only a negligible fraction of it is actually able to leave the DM halo without being possibly re-accreted at later times. On the other hand, the expelled material could progressively inject energy into the halo and thus hamper the gas cooling and its (re-)accretion on the galaxy.

Finally, we detected blue-coloured continuum emission co-spatial with the eastern ionised gas loop (the handle; scales of ∼10 kpc), whose modelled stellar populations are younger (≲100–150 Myr) than in the rest of the galaxy (≳0.5–1 Gyr; Fig. 11). This constitutes possible evidence for widespread star formation triggered at the edge of the bubble due to the compressing action of the jet and outflow (positive feedback), as predicted by theory. The presence of an enhanced gas electron density at the eastern edge of the ionised gas loop, in the direction of the jet and outflow (Fig. 5), is consistent with the compression scenario.

Acknowledgments

We acknowledge funding from ANID programmes FONDECYT Postdoctorado 3200802 (G.V.), FONDECYT Regular – 1190818 (E.T., F.E.B.) and 1200495 (F.E.B., E.T.), CATA-BASAL AFB170002 (G.V., E.T.), ACE210002, and FB210003 (E.T., F.E.B.), and Millennium Science Initiative Program – NCN19_058 (E.T.), European Union’s HE ERC Starting Grant No. 101040227 – WINGS (G.V.), a UK Research and Innovation grant, code: MR/V022830/1 (C.M.H.), DLR grant FKZ 50 OR 2203 (D.T.), and the projects “Feeding and feedback in active galaxies”, with reference PID2019-106027GB-C42, funded by MICINN-AEI/10.13039/501100011033, and “Quantifying the impact of quasar feedback on galaxy evolution”, with reference EUR2020-112266, funded by MICINN-AEI/10.13039/501100011033 and the European Union NextGenerationEU/PRTR (C.R.A.). For the purpose of open access, the authors have applied a Creative Commons Attribution (CC-BY) license to any author accepted version arising. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. This research has made use of the Python packages SciPy (Virtanen et al. 2020), NumPy (Harris et al. 2020), IPython (Perez & Granger 2007), Matplotlib (Hunter 2007), galpy (Bovy 2015), and Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). This research has made use of the FITS files visualisation tools QFitsView, https://www.mpe.mpg.de/~ott/QFitsView/) and DS9 (Joye & Mandel 2003). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. We thank Antonino Marasco and Giulia Tozzi for the useful conversation on the implementation of DM halo potentials. We thank Anna Gallazzi for the insightful discussion on stellar metallicities in galaxies.

References

Appendix A: Stellar kinematics

The stellar velocity and velocity dispersion, obtained from the stellar continuum modelling of the Voronoi-binned cube (described in Sect. 2.2), are reported in Fig. A.1 (left and right panels, respectively). The stellar velocity does not seem to show any clear pattern that could suggest rotation. It shows, however, a general pattern of more redshifted velocities at the outskirts of the stellar emission of the galaxy and lower velocities in the inner, most luminous part of the galaxy. This might be related to the post-merger nature of the source, with the disturbed stellar features of the galaxy (see continuum image in Fig. 9) having different velocities since the galaxy still has to completely relax. The stellar velocity dispersion is also disturbed. It is not radially symmetric and shows an elongation in the same SE-NW direction along which morphological features are seen in the continuum image, in particular extended tails in the SE and S direction and a shell-like feature to the NW.

Fig. A.1. Maps of stellar velocity (left) and velocity dispersion (right), obtained from the stellar-continuum fitting of the Voronoi-binned data cube (see Sect. 2.2).

Appendix B: Additional figures for ionised gas

In Fig. B.1, we show the maps for the [O III]/Hβ, [N II]/Hα, [S II]/Hα, and [O I]/Hα line ratios, from the stellar continuum-subtracted, Voronoi-binned data cube (see Sect. 2.2). [O III]/Hβ is strongest in the direction of the [O III] lobes, while the ratios involving low-ionisation lines ([N II], [S II], and [O I]) are maximal in the perpendicular direction, where the gas velocity dispersion enhancement is observed (Sect. 5).

Fig. B.1. Emission-line ratio maps obtained from the fit of the stellar continuum-subtracted, Voronoi-binned data cube (see Sect. 2.2). [O III]/Hβ (top-left), [N II]/Hα (top-right), [S II]/Hα (bottom-left), and [O I]/Hα (bottom-right) maps for the total modelled line profiles.

Fig. B.2 displays the maps of the [O III] velocity for the total line profile, as well as for the first and the second plus third modelled Gaussian components, from the stellar continuum-subtracted, Voronoi-binned data cube.

Fig. B.2. Ionised gas velocity maps of the Teacup from the fit of the stellar continuum-subtracted, Voronoi-binned data cube (see Sect. 2.2). [O III] velocity, obtained as the first-order moment of velocity of the total modelled line profile (left panel), as well as velocity of the first (middle panel) and second plus third Gaussian components (first-order moment of velocity of their summed profile; right panel) employed for the line profile modelling. In the left panel we report the contours of [O III] flux (from Fig. 3, top-left panel) to allow for an easier visual comparison between the Voronoi-binned velocity field and the morphology of the ionised gas emission.

Fig. B.3 shows the maps of extinction and electron density of the second plus third modelled Gaussian components from the stellar continuum-subtracted, Voronoi-binned data cube, together with their polished versions used to derive the mass outflow rate and related quantities (Sect. 4).

Fig. B.3. Maps of extinction from Hα/Hβ (top-left panel) and electron density from [S II] ratio (bottom-left panel) of the outflow (second plus third) components from the stellar continuum-subtracted, Voronoi-binned data cube, and polished versions of these maps (top-right and bottom-right panels, respectively) used to derive the mass outflow rate and related quantities (see Sect. 4).

Appendix C: Outflow electron density radial profile

In Fig. C.1 we show the radial profile of the outflow electron density, calculated as the mean of the values in each radial bin, weighted by the respective fluxes of [S II] (λ6716+λ6731) of the outflow component(s). We stress that this radial profile does not fully capture the spread of density values at each radial bin, and should therefore not be strictly taken as a quantitative measurement of the outflow density at each radial bin, but more like a qualitative indication of its variation with distance from the nucleus.

Fig. C.1. Radial profile of the ionised outflow electron density, obtained as the mean of the outflow densities (from the [S II] flux ratio of second plus third modelled components), weighted by the respective [S II] (λ6716+λ6731) fluxes of the outflow component(s), in each radial bin. We stress that the reported bars are not statistical errors but a combination of the statistical errors and the standard deviation of the values to show the extent of variation of the densities in each radial bin. Symbols and bar shading have the same meaning as in Fig. 8. Solid and dashed grey lines show the electron density profiles inferred from solving Eq. C.1, by employing the Q(H) radial profiles from Gagne et al. (2014) and Keel et al. (2017), respectively, with a typical uncertainty of ±0.3 dex (shown as light-grey shaded regions).

In order to explore the effects of different density estimates, we compare the electron densities derived from the outflowing portions of the [S II] emission line ratio to the hydrogen densities derived by solving the ionisation parameter equation (e.g. Baron & Netzer 2019; Davies et al. 2020; Revalski et al. 2022). Specifically, we solve for the density at each radius using

$$\begin{array}{c}\hfill n_{\mathrm{H}}=\frac{Q(H)}{4\pi r^{2}cU}\end{array}$$(C.1)

where Q(H) is the number of ionising photons s⁻¹ from the AGN, r the distance of the gas from the AGN, U the ionisation parameter, and c the speed of light (see Revalski et al. 2022 for details). In the case of the Teacup, Q(H) is a function of the distance from the nucleus because the AGN luminosity changes over time (see e.g. Gagne et al. 2014; Keel et al. 2017). We explore using the different Q(H) radial profiles derived by each of the above two studies. In the first case, we adopt Q(H) = C × 9.781 × 10⁵³ s⁻¹, where C = 1 in the nucleus, and increases by steps of 100/16 over the next 16 radial bins, so the last value at r = 16 kpc is Q(H) = 9.781 × 10⁵⁵ s⁻¹ (Gagne et al. 2014). In the second case, we use the variable Q(H) profile derived by Keel et al. (2017), which is estimated based on the Hα emissivity of the gas.

The ionisation parameter is then approximated from our [O III]/Hβ ratio at each radius (each obtained as the mean of the ratios at each radius, weighted by the [O III] flux, considering the outflow components only), which is relatively insensitive to other physical conditions, following Revalski et al. (2022). [O III]/Hβ decreases from ∼8.5 (log(U) = –3.0) in the nucleus to ∼2.5 (log(U) = –3.4) at larger radii. We adopt a smooth monotonic decrease of log(U) by 0.4/16 dex in each of the 16 radial bins from the nucleus up to 16 kpc. This direct estimate of log(U) from [O III]/Hβ, though simplistic, is in general agreement with the models of Gagne et al. (2014).

Finally, we solve Eq. C.1 for the hydrogen gas density using the Q(H) profiles from Gagne et al. (2014) and Keel et al. (2017), and show the resulting radial profiles in Fig. C.1, after reporting them to electron densities, being n_H/n_e ≃ 0.85 (e.g. Crenshaw et al. 2015; Revalski et al. 2018).

We see a general agreement between the various profiles at some radii, with small systematic offsets due to the different assumptions in each technique. Overall, the density profiles obtained from Eq. C.1 tend to be a bit larger (by ∼0.5 dex on average) than that from [S II]. There is good agreement in the inner regions (≲3 kpc), where the outflow is strongest, except in the innermost radial bin, where the outflow electron density from [S II] was poorly constrained and was therefore not used for quantitative purposes (see Sect. 4). We stress that, as Fig. C.1 shows, differences are not systematic in the overall profile. At the same time, uncertainties of this magnitude are expected due to the assumptions made by each individual method.

We emphasise, however, that Fig. C.1 must not be taken as a strict quantitative comparison between the density estimates at each radius from different techniques, but rather as an overall qualitative one, given the uncertainties and different assumptions involved in each of them. For example, even the same method (solving Eq. C.1) gives quite different densities, depending on the employed Q(H) radial profile (from either Gagne et al. 2014 or Keel et al. 2017; solid and dashed lines, respectively). A more detailed comparison of different gas density measurement techniques is beyond the scope of this work.

Appendix D: Very broad Paα line emission

Tentative detection of a very broad (FWHM ∼ 3000 km s⁻¹) Paα component across the base of the NE bubble and the SW fan was reported by Ramos Almeida et al. (2017), who did not detect it in the other near-IR emission lines, though. We do not detect such a broad component in the optical emission lines. Being the flux of its peak about 10% the peak of the main body of the line, this broad component should be visible in the MUSE optical spectra, as it would be much stronger than the noise level, assuming that the above relative proportion between the broad component and the main body of the line holds also in the optical. Nevertheless, the centroid velocity of these very broad components, with respect to that of the peak of the narrower component, have the same sign that we also observe for the broad components in the ionised gas, that is, blueshifted to the NE and redshifted to the SW of the nucleus (see panels 1 and 3 in Fig. 2), which only reach FWHMs of ∼1300 km/s, though. Therefore, the very broad Paα components reported in Ramos Almeida et al. (2017) could be real, rather than being only an instrumental artefact in the SINFONI data, but their extreme broadness could be the result of the high noisiness of the data, affecting especially the line wings, which would have artificially broadened the modelled component.

All Figures

Fig. 1. Optical (emission lines and continuum) versus radio emission of the Teacup galaxy. Left panel: false-colour image from VLT/MUSE. [O III] ionised gas emission is reported in green, Hα in red, and line-free continuum emission averaged between 5500–6700 Å (observed wavelengths) in blue. The colour intensity scale is the same for [O III] and Hα. The main spatial features mentioned in the text are labelled. Right panel: 20″ × 20″ zoomed-in [O III] emission (same as in left panel), with the contours of the VLA 5.12 GHz (white) and highest-resolution 6.22 GHz (magenta) radio images from Harrison et al. (2015) overlaid. For the former, contour levels are 35.0, 66.0, 124, and 234 μJy beam−1; for the latter, they are 40, 102, 261, and 668 μJy beam−1. The black numbered “+” symbols mark the regions from which the spectra reported in Fig. 2 were extracted.

In the text

Fig. 2. Example of modelled spectra from different regions, covering Hβ and [O III] (upper sub-panels) and [N II], Hα, and [S II] (lower sub-panels). From top to bottom, the spectra come from close to the nucleus, to the NE of it (number 1), from the eastern lobe or Teacup handle (number 2), and from the western lobe (number 3); these regions are labelled accordingly in Fig. 1, right panel. For each emission line, blue, green, and orange curves mark the first (narrower), the second, and the third (broader) Gaussian components employed in the fitting, respectively.

In the text

Fig. 3. Maps of the [O III] (top panels) and Hα (bottom panels) emission line flux distribution in the Teacup galaxy. The fluxes of the total modelled line profile (left panels; same reported in Fig. 1) and of the first, narrower (centre panels) and the outflow (second plus third), broader Gaussian components (right panels) making up the total profile are reported. The colour intensity scale is the same for all the three maps of each line, as reported in the colour bar to the right. The contours in upper-left panel mark the VLA radio emission at 5.12 GHz from Harrison et al. (2015), same as in Fig. 1, right panel.

In the text

Fig. 4. Ionised gas kinematic maps of the Teacup. [O III] velocity (upper panels) and velocity dispersion (middle panels), obtained as first- and second-order moments of the velocity, respectively, of the total modelled profile of [O III] (left panels), as well as for the first (central panels) and the second plus third Gaussian components (first- and second-order moments of velocity of their summed profile; right panels) employed for the line profile modelling. Central and right panels are zoomed in the central 20″ × 20″ (indicated by the black box in left panels). Bottom panels: zoomed versions (also in the central 20″ × 20″) of the [O III] flux (left; same as in Fig. 3, top-left panel), velocity (centre; same as in top-left panel) and velocity dispersion (right; same as in mid-left panel) maps for the total modelled emission line profile. The contours indicate the same VLA 5.12 GHz and highest-resolution 6.22 GHz radio images from Harrison et al. (2015) reported in Fig. 1, right panel, with the same contour levels.

In the text

Fig. 5. Maps of physical properties of the ionised gas in the Teacup. Extinction AV from Hα/Hβ (top-left) and Hγ/Hβ (top-right), electron density (from [S II] λ6716/λ6731; bottom-left), and electron temperature (from [O III] (λ5007+λ4959)/λ4363; bottom-right). The line fluxes involved in the line ratios employed for the maps are those of the total modelled emission line profiles (with one, two, or three Gaussians). The maps are zoomed in the central 20″ × 20″.

In the text

Fig. 6. Spatially resolved diagnostic diagrams of gas excitation for the Teacup, obtained from the star-subtracted Voronoi-binned data cube. [O III] λ5007/Hβ versus [N II] λ6584/Hα (top-left), versus [S II] (λ6716+λ6731)/Hα (top-centre), and versus [O I] λ6300/Hα (top-right) diagrams, together with their respective spatial distributions (bottom panels). The fluxes of the total modelled emission line profiles have been employed for the diagrams. The solid curve defines the theoretical upper bound for pure star formation from Kewley et al. (2001), while the dashed one in the [N II]-diagram is the Kauffmann et al. (2003) empirical classification. The dot-dashed line represents the demarcation line between Seyfert galaxies and shocks or LI(N)ERs from Schawinski et al. (2007) in the [N II]-diagram and from Kewley et al. (2006) in the [S II]- and [O I]-diagrams. Regions dominated by AGN ionisation are then marked in red, composite regions in the [N II]-diagram in purple, shock-ionised or LI(N)ER regions in green, and star formation-dominated regions would be marked in blue (when present).

In the text

	Fig. 7. Same as in Fig. 6, but separately for narrower (first) and outflow (second plus third) modelled components.
In the text

Fig. 8. Spatially resolved properties of the galactic ionised outflow in the Teacup. Top-left: mass outflow rate map, obtained from the outflow (second plus third) modelled component(s) of Hα. Top-right: radial profiles of mass outflow rate (top), kinetic rate (middle), and momentum rate (bottom) as a function of distance from the nucleus, from Hα (red and orange) and [O III] (blue and light blue) outflow components. These quantities were obtained by summing up the single-spaxel values (calculated through Eqs. (3), (5), and (6)) contained in each radial bin, and re-scaling by ΔR/ΔRradbin, being ΔRradbin the radial bin width (=1 kpc) and ΔR the spaxel size. Orange and light blue points mark those external radial bins for which we adopted a single median gas density, since no density measurement was available (see Sect. 4 for details). A different symbol is used for the innermost radial bin, which comprises a few spaxels whose outflow density could not be constrained (see Sect. 4). Lighter error bars are used when the assigned statistical uncertainties were estimated from representative regions instead of on a spaxel basis (see Sect. 4). The vertical dashed lines mark the distance of the handle. On the rightmost y axis, we also report mass-loading factor (top), coupling efficiency (middle), and momentum boost (bottom) of the outflow, for which we consider SFR ∼ 10 M⊙ and Lbol ∼ 5 × 1045 erg s−1 (see text). Mid-left: radial profiles of ionised outflow mass, from the flux of the outflow component(s) of Hα and [O III]. Bottom: radial profiles of outflow velocity (left; vout = v2 + 3 + FWHM2 + 3/2 ≃ 1.18σ2 + 3), combination of velocity, v2 + 3 (centre), and velocity dispersion, σ2 + 3 (right), of second plus third modelled components, for both [O III] and Hα. The reported values are the flux-weighted (by the outflow component(s)) averages in each radial bin. We stress that the reported error bars are not statistical errors, but are a combination of the statistical errors and the standard deviation of the values, to show the extent of variation of the velocities in each radial bin.

In the text

	Fig. 9. Map of continuum emission, obtained by averaging the (observed) spectral channels of the data cube between 5500–6700 Å, free from the strongest gas emission lines (same as in Fig. 1, left panel, here zoomed in the central 40″ × 40″).
In the text

Fig. 10. Zoomed-in maps of median stellar continuum in selected spectral ranges free of gas emission lines. The spectral range in which the median has been calculated are indicated by the light grey shaded area in the spectra below. The spectra are extracted by summing all the spaxels in the handle, as traced by [O III] emission, whose contours are reported in top-left image for comparison; the extraction aperture of the spectra is shown in Fig. 11, top-left panel, labelled as number 1.

In the text

Fig. 11. Properties of the stellar populations, showing evidence for young stars in the Teacup handle. Top-left: zoomed-in map of the stellar continuum obtained by performing the median in all the six spectral ranges reported in Fig. 10 together (light grey shaded areas). Mid-left: map of continuum flux ratio (i.e. colour) between the images corresponding to the bluest and the reddest spectral ranges reported in Fig. 10 (top-left and bottom-right panels, respectively). Prior to performing the ratio, the continuum images were smoothed with a Gaussian kernel having σ = 1.5 spaxels, for a better visual output. The patch having low continuum flux ratio in the southern part of the handle is due to a background red galaxy. Top-right: luminosity-weighted fractions of stellar populations in the stellar metallicity (in log scale, with respect to solar metallicity) versus age plane resulting from the stellar population modelling (details in text) of each numbered aperture shown in top-left panel. Bottom: spectrum (black) extracted from the aperture covering the handle (number 1 in top-left panel) with the modelled stellar continuum overlaid (red); shaded grey vertical bands indicate spectral ranges which were masked in the fitting, and the spectrum within these ranges is marked in blue.

In the text

	Fig. A.1. Maps of stellar velocity (left) and velocity dispersion (right), obtained from the stellar-continuum fitting of the Voronoi-binned data cube (see Sect. 2.2).
In the text

	Fig. B.1. Emission-line ratio maps obtained from the fit of the stellar continuum-subtracted, Voronoi-binned data cube (see Sect. 2.2). [O III]/Hβ (top-left), [N II]/Hα (top-right), [S II]/Hα (bottom-left), and [O I]/Hα (bottom-right) maps for the total modelled line profiles.
In the text

Fig. B.2. Ionised gas velocity maps of the Teacup from the fit of the stellar continuum-subtracted, Voronoi-binned data cube (see Sect. 2.2). [O III] velocity, obtained as the first-order moment of velocity of the total modelled line profile (left panel), as well as velocity of the first (middle panel) and second plus third Gaussian components (first-order moment of velocity of their summed profile; right panel) employed for the line profile modelling. In the left panel we report the contours of [O III] flux (from Fig. 3, top-left panel) to allow for an easier visual comparison between the Voronoi-binned velocity field and the morphology of the ionised gas emission.

In the text

	Fig. B.3. Maps of extinction from Hα/Hβ (top-left panel) and electron density from [S II] ratio (bottom-left panel) of the outflow (second plus third) components from the stellar continuum-subtracted, Voronoi-binned data cube, and polished versions of these maps (top-right and bottom-right panels, respectively) used to derive the mass outflow rate and related quantities (see Sect. 4).
In the text

Fig. C.1. Radial profile of the ionised outflow electron density, obtained as the mean of the outflow densities (from the [S II] flux ratio of second plus third modelled components), weighted by the respective [S II] (λ6716+λ6731) fluxes of the outflow component(s), in each radial bin. We stress that the reported bars are not statistical errors but a combination of the statistical errors and the standard deviation of the values to show the extent of variation of the densities in each radial bin. Symbols and bar shading have the same meaning as in Fig. 8. Solid and dashed grey lines show the electron density profiles inferred from solving Eq. C.1, by employing the Q(H) radial profiles from Gagne et al. (2014) and Keel et al. (2017), respectively, with a typical uncertainty of ±0.3 dex (shown as light-grey shaded regions).

In the text