MUSE view of PDS 456: kpc-scale wind, extended ionized gas and close environment

PDS 456 is the most luminous ($L_{\rm bol}\sim 10^{47}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}}$) radio-quiet quasar at $z<0.3$ and can be regarded as a local counterpart of the powerful quasars shining at Cosmic Noon. It hosts a strong nuclear X-ray ultra-fast ($\sim 0.3c$) outflow, and a massive and clumpy CO (3-2) molecular outflow extending up to $\sim$5 kpc from the nucleus. We analyzed the first MUSE Wide Field Mode (WFM) and Adaptive-Optics Narrow Field Mode (AO-NFM) optical integral field spectroscopic observations of PDS456. The AO-NFM observations provide an unprecedented spatial resolution, reaching up to $\sim$280 pc. Our findings reveal a complex circumgalactic medium around PDS 456, extending up to a maximum projected size of $\approx$46 kpc. This includes a reservoir of gas with a mass of $\sim 10^{7}-10^{8}\rm\leavevmode\nobreak\ M_{\odot}$, along with eight companion galaxies, and a multi-phase outflow. WFM and NFM MUSE data reveal an outflow on a large scale ($\approx$12 kpc from the quasar) in [O III], and on smaller scales (within 3 kpc) with higher resolution (about 280 pc) in H$\alpha$, respectively. The [O III] outflow mass rate is $\rm 2.3\pm 0.2\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\ yr^{-1}$ which is significantly lower than those typically found in other luminous quasars. Remarkably, the H$\alpha$ outflow shows a similar scale, morphology, and kinematics to the CO (3-2) molecular outflow, with the latter dominating in terms of kinetic energy and mass outflow rate by two and one orders of magnitude, respectively. Our results therefore indicate that mergers, powerful AGN activity, and feedback through AGN-driven winds will collectively contribute to shaping the host galaxy evolution of PDS 456, and likely, that of similar objects at the brightest end of the AGN luminosity function across all redshifts. Moreover, the finding that the momentum boost of the total outflow deviates from the expected energy-conserving expansion for large-scale outflows highlights the need of novel AGN-driven outflow models to comprehensively interpret these phenomena.

1 Introduction

Most of the massive galaxies host a supermassive black hole (SMBH) in their center, that accrete material during a phase called Active Galactic Nucleus (AGN), releasing a large amount of energy affecting the host-galaxy itself. We call this phenomenon “AGN feedback”. SMBHs reach their maximum accreting and feedback efficiency at $z\sim 2-3$, which we refer to as “Cosmic Noon”, i.e. the peak epoch of galaxy assembly and SMBH accretion (see Förster Schreiber & Wuyts, 2020, and references therein). Understanding AGN feedback is essential for explaining key observational findings like SMBH-galactic correlations (Magorrian et al., 1998), star formation inefficiency at high-$\rm M_{\star}$ (Binney & Tabor, 1995), and intergalactic medium enrichment (Tytler et al., 1995).

During the AGN phase, the released energy triggers semi-relativistic nuclear winds via radiative pressure-driven (Proga et al., 1998, 2000) or magnetic pressure-driven (Fukumura et al., 2010; Luminari et al., 2021) mechanisms, leading to X-ray (Ultra Fast Outflow; UFO Reeves et al., 2003) and UV broad absorption line (Jannuzi et al., 1996; Vietri et al., 2022). These winds shock the Inter-Stellar Medium (ISM) gas, generating galactic-scale, multi-phase outflows (King, 2005; Lapi et al., 2005; Menci et al., 2008; Faucher-Giguère & Quataert, 2012).

If the energy of the shocked wind is radiated away and the momentum is conserved (momentum-conserving scenario; Fabian, 1999; King, 2003), outflows are expected to extend up to a few hundred parsecs (Feruglio et al., 2017; Zanchettin et al., 2021). An energy-conserving scenario is proposed to explain the expansion of kiloparsec-scale outflows (Silk & Rees, 1998; King, 2005; Zubovas & King, 2012; Nayakshin, 2014; Fiore et al., 2017, hereafter F17). In this context, the nuclear wind undergoes adiabatic expansion, facilitated by a cooling time for the post-shock outflow that surpasses its dynamical time, thus ensuring energy conservation. This scenario remains valid beyond a critical distance of more than $100\leavevmode\nobreak\ \rm pc$ from the quasar resulting in outflows exhibiting a momentum boost of $\sim$10-20 (Zubovas & Nayakshin, 2014a; King & Pounds, 2015). However, in a two-temperature plasma scenario, where ions decouple thermally from electrons due to the shock, the cooling radius is much smaller than $100\leavevmode\nobreak\ \rm pc$ (Faucher-Giguère & Quataert, 2012). Moreover, Zubovas & Nayakshin (2014b) demonstrate that the efficiency in depositing energy into the ambient gas also depends on the clumpiness of the gas.

More recent observational investigations (Bischetti et al., 2019; Smith et al., 2019; Reeves & Braito, 2019; Sirressi et al., 2019; Marasco et al., 2020; Tozzi et al., 2021; Speranza et al., 2022; Bonanomi et al., 2023) have revealed a more complex scenario, supporting a momentum-conserving scenario for both molecular and ionized galactic-scale outflows, at least in some cases. Still other AGN driven large-scale outflows show momentum loading factors that fall below even the momentum-conserving prediction, while others have extremely high factors, above the energy-driven expectations (Marasco et al., 2020). Several hypotheses have been put forward to explain these observations. One scenario predicts that galactic-scale outflows could be primarily driven by the radiative pressure from AGN photons on dust (Fabian, 2012; Ishibashi et al., 2018; Costa et al., 2018). This mechanism can account for galactic-scale outflows maintaining a low momentum boost. It predicts a super-linear correlation (i.e. correlation with a slope greater than one) between the kinetic power of large-scale outflows and AGN luminosity; such a correlation is supported by outflow observations (e.g., F17, Fluetsch et al. 2021, Lutz et al. 2020). On the other hand, this model does not explain the spread of outflow properties in different AGN with the same luminosity. Another possible explanation relies on considering an intermittent AGN luminosity history (see Nardini & Zubovas, 2018; Zubovas & Nardini, 2020; Zubovas et al., 2022). These papers rely entirely on the energy-driven outflow paradigm but show that outflow properties correlate much more strongly with the long-term (several Myr) average AGN luminosity rather than the instantaneous one. The observed momentum and energy loading factors are then determined to a large extent by the ratio of the current AGN luminosity to the long-term average.

In this paper, we use MUSE Adaptive Optics (AO) -assisted Narrow Field Mode (NFM), and Wide Field Mode (WFM) observations to explore the properties of the environment and the ionized outflow in PDS 456, that is the most luminous (bolometric luminosity; $L_{\rm bol}\sim 10^{47}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}}$) radio-quiet quasar at $z<0.3$ ($z_{\rm CO}=0.185$; Bischetti et al., 2019, hereafter B19), discovered by Torres et al. (1997) in the Pico dos Dias survey. The analysis of the Spectral Energy Distribution (SED) reveals both the quasar-like and the ULIRG-like nature of PDS 456, suggesting that it is undergoing a transition from a Luminous IR Galaxy (LIRG) to a quasar (Yun et al., 2004). In B19, the star formation rate (SFR) of the host galaxy is estimated to be $\approx 30-80\leavevmode\nobreak\ \rm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$. In addition, Yun et al. (2004) conducted an analysis using multiple datasets, including a 0.6^′′ resolution K-band image obtained from the Keck Telescope, a VLA continuum image at 1.2 GHz, and the CO(1-0) emission with the Owens Valley Radio Observatory (OVRO) millimeter array. Their study revealed the presence of three compact continuum sources located $\approx 10\leavevmode\nobreak\ \rm kpc$ southwest of the quasar. Furthermore, in their ALMA observations, over a region of $\rm\sim 50\times 50\leavevmode\nobreak\ kpc^{2}$, B19 identified CO (1-0)-emitting companions near PDS 456.

Simpson et al. (1999) were the first to analyze the optical spectrum, revealing the presence of broad line region (BLR) Balmer lines (FWHM $>4000\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$), a weak [O III] emission line, and prominent Fe II transitions. PDS 456 can be seen as the local counterpart of the luminous quasars shining at $z\gtrsim 2$, where we expect a high efficiency in driving radiative winds (e.g. Brusa et al., 2015; Carniani et al., 2015; Bischetti et al., 2017; Förster Schreiber et al., 2018; Kakkad et al., 2020). Indeed, PDS 456 exhibits a relativistic ($\sim 0.3c$) and powerful ultra fast X-ray winds (kinetic power; $\dot{E}_{\rm{kin}}\simeq 0.2L_{\rm{bol}}$), detected through Fe K-shell absorption features (Reeves et al., 2009; Nardini et al., 2015; Luminari et al., 2018). Hamann et al. (2018) reported the possible presence of a fast ($\sim 0.3c$) UV broad absorption line wind in C IV, exhibiting similar velocities to the nuclear X-ray winds. Additionally, O’Brien et al. (2005) measured a $\sim 5000\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$ blueshift of the C IV emission line with respect to the systemic redshift, likely associated with an outflow in the quasar broad line region. The presence of a molecular outflow was detected in ALMA data by B19. They observed a blueshifted ($<-250\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$) CO (3-2) outflow extending eastward from the quasar, with a maximum projected distance of $\sim$5 kpc. Additionally, a more compact ($<1\leavevmode\nobreak\ \rm kpc$) west-oriented outflow was observed, exhibiting low positive velocities ($\sim 80-150\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$). The molecular outflow has a low momentum boost (i.e. $\dot{P}_{\rm{out}}/\dot{P}_{\rm{rad}}\approx 0.36$, where $\dot{P}_{\rm{out}}=v_{\rm{max}}\times\dot{M}_{\rm{out}}$, with $\dot{M}_{\rm{out}}$ the mass rate of the outflow and $v_{\rm{max}}$ the maximum velocity of the outflow, and $\dot{P}_{\rm{rad}}=L_{\rm bol}/c$) which is not consistent with an energy-conserving scenario. Indeed, in $log(L_{\rm bol}/\rm erg\leavevmode\nobreak\ s^{-1})\gtrsim 47$ quasars, the energetic contribution of the ionized outflows is found to be comparable to that of the molecular phase, thus $\dot{M}_{\rm out,ion}\sim\dot{M}_{\rm out,mol}$ for $L_{\rm bol}=10^{47}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}}$ (F17). Nevertheless, although the mass outflow rate measured for the molecular outflow ($\dot{M}_{\rm out}\sim 290\rm\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\ yr^{-1}$) is well below the prediction from the empirical relation by F17, it is sufficient to deplete the molecular gas reservoir in $\rm\sim 10\leavevmode\nobreak\ Myr$, leading to a rapid suppression of the star formation activity.

Finally, although PDS 456 is classified as a radio-quiet object, Yang et al. (2021) recently reported the detection of a faint ($L_{1.66\leavevmode\nobreak\ \rm{GHz}}<10^{40.2}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}}$), complex radio structure. This consists of collimated jets, and an extended (up to 360 pc) non-thermal radio emission, which is explained as shock emission due to the interaction between nuclear winds and high-density regions of the ISM. The energy associated with these radio structures ($L_{R}/L_{\rm{bol}}=10^{-7}$) is significantly lower compared to that of the nuclear winds (Nardini et al., 2015; Luminari et al., 2018).

In conclusion, PDS 456 is an excellent laboratory to study the environment and the feeding and feedback processes involved in the typical scenario of efficient AGN phase observed at Cosmic Noon around hyper-luminous quasars. In particular, the large FoV of the MUSE WFM Integral Field Spectroscopy (IFS) data is crucial to study the environment of this system traced by the ionized phase, i.e. the Narrow Line Region (NLR), host-galaxy, companions and diffuse emission. On the other hand, the MUSE AO-NFM data of the PDS 456 center provides observations with unprecedented spatial-resolution ($280\leavevmode\nobreak\ \rm pc$) for this type of system, with the aim to study the mechanisms of expansion of any possible ionized-phase outflow associated with the molecular one detected with ALMA in B19.

Sect. 2 describes the MUSE WFM and NFM observations, while the analysis of these data and the corresponding results are reported in Sects. 3 and 4, respectively. In these sections, we investigate the environment of PDS 456, as well as the morphology and kinematics of the extended emitting gas and the presence of companion galaxies.

Sect. 5 explores the physical and energetic properties of the outflow detected in both observation modes across different scales. Finally, in Sect. 6, we discuss possible scenarios for the expansion of these outflows based on their properties. A summary of the main results is provided in Sect. 7.

To be consistent with B19, we consider the systemic redshift of PDS 456 of $z=0.185$, as well as a physical scale of 3.126 kpc/arcsec at the following cosmological parameters: $H_{0}=69.6\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}$, $\Omega_{m}=0.286$, and $\Omega_{\Lambda}=0.714$.

Figure 1: RGB wide filters images of the WFM (left) and NFM (right) MUSE datacube, obtained by collapsing the MUSE datacubed in the same spectral regions of the HST/ACS-HCR filters F475W, F625W and F775W for the WFM data, and F459M, F660N and F892N for the NFM data. For the NFM data, spectral regions contaminated by instrument features fallen in the filters are excluded. In the WFM image, we highlight the region studied in this work, including PDS 456 and companions, with a dashed white box measuring $\rm 80\times 60\leavevmode\nobreak\ kpc^{2}$. The cyan box indicates the size ($\approx 23\times 23\leavevmode\nobreak\ \rm kpc^{2}$) of the NFM FOV centered on PDS 456. The positions of the nucleus of PDS 456 are marked with black crosses.

2 Data Reduction

In this paper we analyze MUSE IFS observations in WFM and in AO-NFM taken on April to June 2019 (PI: E. Piconcelli). The MUSE WFM observations consisted of four Observing Blocks (OBs) of 600 sec each, for a total exposure time of 1 h. Each exposure was rotated by 90 deg with the addition of a small dithering. Instead, the MUSE AO-NFM observations consist of three OBs for a total observing time of $\sim$3 h, and each OB consists of 8 exposures including sky pointings. The observations were acquired with seeing and airmass $\sim$0.4”-1” and $\sim$1.1, respectively.

Data were reduced using the ESO MUSE standard pipeline with EsoRex v. 3.12.3 (Weilbacher et al., 2014). The sky frame for MUSE WFM raw data was generated directly from the “OBJECT” exposures using the standard ESO pipeline. For AO-NFM raw data, the sky frame was created from the pixel tables of exposures taken in empty sky regions, i.e. the “SKY” exposures.

The final WFM (NFM) datacube is characterized by a FoV of $\approx 1^{\prime}\times 1^{\prime}$ ($\approx$7$\times$7 $\rm arcsec^{2}$), and a pixel size of 0.2 (0.025) arcsec, that is equivalent to $\sim$625 pc ($\sim$78 pc) at $z=0.185$. The FWHM of the final datacubes’ PSF are $\sim 1^{\prime\prime}$ (i.e. $\sim 3\leavevmode\nobreak\ \rm kpc$) and $\sim 0.09^{\prime\prime}$ (i.e. $\sim 280\leavevmode\nobreak\ \rm pc$) for WFM and NFM observations, respectively. These values are estimated based on the individual point-sources within the fields. The spectral range covered is from 4750 Å to 9350 Å, with a spectral bin of 1.25 Å. However, the presence of the Na ID notch filter, which is required for the laser guide star in the AO mode, significantly contaminates the spectral region ranging from approximately 5500 to 6000 Å in the MUSE NFM data. Consequently, we did not utilize this contaminated portion of the data. We verified the absolute wavelength calibration in our final datacube by checking the positions of the most intense sky lines. Finally, we checked the calibration of the flux in NFM-MUSE data comparing the total PSF of PDS 456 with that in WFM-MUSE data as reference.

Fig 1 shows the RGB images of PDS 456. In the left panel, we present an image of the entire FoV from the WFM-MUSE data. To create this image, we overlay mimic-wide filters images that are generated by collapsing the spectrum of the datacube into three spectral regions simulating the filters in HST/ACS-HCR: F475W (B), F625W (G) and F775W (R). This allows us to obtain a rough approximation of the color of each source within the FoV. The right panel displays the NFM MUSE observation, presented as an RGB image overlaid with mimic-wide filters images collapsed in the spectral regions of the HST/ACS-HCR filters F459M (B), F660N (G) and F892N (R).

We use the positions of a handful of point sources identified by GAIA within the WFM MUSE observation field to refine the astrometry of our data. The astrometry of NFM MUSE data is automatically adjusted through alignment with the peak of the continuum emission from PDS 456.

Figure 2: The extended emission surrounding PDS 456, within the region of the WFM-MUSE FoV indicated by a dashed white box in the left panel of Fig. 1. This RGB composite image is generated from the following images: optimally-extracted narrow band images of H$\alpha$ (R) and [O III]$\lambda 5008$Å (B) obtained using the CubEx tool, and the White-Light Image (WLI; G). The black cross symbols mark the line-emitting companions detected with the Keck telescope (Yun et al., 2004), ALMA (B19), and in this work with MUSE (see Table 1 for details). MUSE continuum sources are visible in green. The white contours show the surface brightness levels of the [O III] nebula at 10, 50, 100 and 500 $10^{-18}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ cm^{-2}\leavevmode\nobreak\ arcsec^{-2}}$.

3 WFM observations

3.1 Subtraction of the Nuclear and Continuum Emission

We analyze the MUSE WFM observations with the purpose to investigate the properties of the kpc-scale diffuse gas around PDS 456 by subtracting the BLR emission, the iron transitions and the continuum emission (hereafter we refer to the combination of these emission components as ”the nuclear component”) from the region associated with the quasar PSF. For the WFM observations this task is complicated by the presence of several continuum sources blended with the quasar. For this reason we do not employ for the WFM the standard PSF subtraction tools used in literature (e.g., CubePSFSub see later for the discussion) and we develop a custom procedure as described in detailed below.

The initial step consists of extracting the spectrum from a small circular region centered on PDS 456 in the assumption that this area is completely dominated by the nuclear component. In order to derive a spectral model for this emission, we adopt an empirical approach since a purely analytical model results in large residuals, particularly from iron templates. The empirically calibrated model is obtained by performing a spline interpolation of the spectrum using third-degree polynomial after carefully masking the spectral regions containing narrow emission lines. In order to pinpoint these spectral regions, we rely on a iterative procedure using visual inspection. The model spectral shape is fixed while we vary the normalization for each spaxel through spectral fitting based on $\chi^{2}$ minimization after masking the spectral regions associated with narrow emission lines. Second order effects due to, e.g. wavelength PSF variations are accounted for with an additional spectral component modeled with a 7-degree polynomial. We stress that this additional component has minimal effects on our final result. Finally, the obtained spectral models are subtracted spaxel by spaxel.

To evaluate any potential overestimation or underestimation of the narrow emission lines resulting from the subtraction of the best-fit nuclear component of PDS 456, we examine the [N II]$\lambda 6583$Å/[N II]$\lambda 6548$Å and [O III]$\lambda 5008$Å/[O III]$\lambda 4960$Å line ratios. Our analysis indicates that, within the $1\sigma$ uncertainties, these ratios match with theoretical predictions ($\approx 1/3$; Storey & Zeippen, 2000). As a further check we verify that the spatial extension of the nuclear component emission matches with the PSF estimated from some point sources in the WFM MUSE field. As a final step, we follow the same procedure described in Sect 3.2 from Borisova et al. (2016) (see also Cantalupo et al. 2019 for further details) in order to eliminate both any possible residual emission resulting from the subtraction of the nuclear component and any continuum sources across the entire WFM field of view.

3.2 Morphology and extension of the ionized gas emission

Table 1: List of companions detected around PDS 456. Col. (1) source name; col. (2) and (3): coordinates in degrees; col. (4) and (5): projected distances in kpc from PDS 456 and redshift, respectively; col. (6): tracers, such as continuous lines or bands, were employed to identify the redshift of the sources.

Src RA DEC D(kpc) $z$ tracers (1) (2) (3) (4) (5) (6) K1 262.0818 -14.26599 9 0.1847^a Continuum K-band, Emission line CO(3-2), Absorption lines Na ID K2 262.0820 -14.26629 10 0.1845^a Continuum K-band, Emission line CO(3-2), Absorption lines Na ID K3 262.0821 -14.26640 11 0.1837^a Continuum K-band, Emission line CO(3-2), Absorption lines Na ID M1 262.0826 -14.2665 11 0.1847 Rest-frame optical emission lines and continuum in MUSE M2 262.0809 -14.2666 21 0.1829/0.1836^b Emission line CO(3-2), rest-frame optical emission lines and continuum in MUSE M3 262.0802 -14.2679 38 0.1838 Rest-frame optical emission lines and continuum in MUSE M4 262.0802 -14.2653 25 0.1843 Rest-frame optical emission lines and continuum in MUSE M5 262.0800 -14.2647 28 0.1840 Rest-frame optical emission lines and continuum in MUSE The sources have been detected with ALMA observations in B19, with the Near-Infrared Camera on the W. M. Keck Telescope in Yun et al. (2004), and with the IFS MUSE-WFM data in this work. ^a The redshift is obtained by fitting the Na ID absorption feature. ^b The first (second) value indicates the CO (1-0) (H$\alpha$) based redshift;

To obtain a 3D map of the ionized gas emission, we use the CubExtractor tool on the datacube after subtracting the nuclear component and the continuum. This tool generates a 3D-mask comprising a minimum number of connected voxels (i.e. spatial and spectral elements; MinNVox=10000) above a user-defined threshold, and whose line emission is above a minimum signal-to-noise ratio (SN_Threshold=3). These values have been widely used in the literature (e.g., Borisova et al., 2016; Arrigoni Battaia et al., 2019; Cantalupo et al., 2019) to extract extended emission. This method allows us to trace the circumgalactic medium (CGM) in PDS 456 through the extended emission in the H$\beta$, H$\alpha$, [N II]$\lambda 6548$Å, [N II]$\lambda 6583$Å  [O III] and [S II] optical transitions. We produce the optimally-extracted Narrow Band (NB) images for each detected emission line. These are images obtained by collapsing the only voxels in which the signal-to-noise ratio (SNR) of the emission line is above SN_Threshold we set.

Fig. 2 shows an RGB image produced by including optimally-extracted NB images of the extended H$\alpha$ (red) and [O III] (blue) emission and the White Light Image (WLI) highlighting the continuum (green) of the original datacube. This image reveals a complex morphology of the CGM around PDS 456, extending up to a maximum projected size of $\approx$46 kpc.

Figure 3: The top panel displays the spectra and their best-fit models, while the bottom panel presents the [N II]-BPT diagram, for the galaxy companions of PDS 456. These spectra were extracted from circular regions with a radius of 2 pixels in the datacube, where the nuclear and continuum emission components of PDS 456 have been subtracted. Transparent gray dots show the values from the SDSS survey (credits: Jake Vanderplas $\&$ AstroML Developers), The black lines dividing different types of emissions in the BPT diagram were taken from Kewley et al. (2001).

3.3 Multiple Companions of PDS 456

The CGM emission shown in Fig. 2 also reveals the existence of components resembling bridges that connect various sources marked with black crosses. These sources share the same redshifts as PDS 456, confirmed through the detection of emission and absorption lines with ALMA (B19), and with MUSE data in this study. A comprehensive list of companion sources, including their coordinates, distances from PDS 456, redshifts, and associated tracers, is provided in Table 1.

Prominent continuum emission is observed from three sources near PDS 456, i.e. K1, K2, and K3, which were also detected in the K-band continuum using the Near-Infrared Camera on the Keck Telescope (Yun et al., 2004). MUSE data reveal Na ID absorption features, serving as proxies for neutral gas. The redshifts of these three sources are estimated by fitting the Na ID absorption line using the model adopted in Sato et al. (2009) and Perna et al. (2020). The best-fit models are presented in Fig. 13 in Appendix A. The redshifts derived in this way for K1 and K3 match with those obtained from CO (3-2) in B19. The absence of these Na ID absorption features at the center of PDS456 is consistent with other luminous unobscured AGN (Rupke et al., 2005; Villar Martín et al., 2014; Perna et al., 2017, 2019). Regarding the other continuum sources highlighted in green in Fig. 2, no distinctive features can be discerned to determine their redshifts. A more detailed analysis for this purpose is deferred to future investigations. The MUSE spectra of five sources, i.e. M1–M5, show emission lines and and a less prominent continuum emission. Notably, the continuum emission for these sources is not visible in Fig. 2, but it is observed in the extracted spectra. Therefore, all eight sources listed in Table 1 exhibit both emission/absorption lines and continuum, identifying them as galaxy companions. Notably, the majority of companions are situated to the west of the central quasar at distances ranging from 9 to 38 kpc.

To determine the nature and properties of most of the PDS 456 companions detected with MUSE, we employ the Baldwin, Phillips, and Terlevich (BPT Baldwin et al., 1981) diagram. As done above, we extract the same spectra from the datacube nuclear- and continuum-subtracted (Fig. 3, top panel). The emission lines are fit with Gaussian profiles, and the emission line ratios are derived and plotted on the BPT diagram (Fig. 3, bottom panel), with the color bar showing the distance in kpc from PDS 456. The BPT diagram indicates that M1, M4, and M5 fall within the AGN-dominated region, suggesting a significant contribution of AGN illumination. M2 is classified as HII regions, indicating ongoing star formation. Unfortunately, for the companion M3, information on the y-axis is unavailable, and we are unable to assign a specific classification. For the sources M1, M4, and M5, we employ the formula from Keel et al. (2012) to calculate the lower limit of the incident ionizing flux necessary for the observed H$\beta$ line emission, taking into account the distance of these companions from PDS 456. These estimates suggest that the ionizing flux originating from PDS 456 can explain the high ionization observed in the emission lines of all these companions. Therefore, their position on the BPT diagram is not attributed to the presence of an AGN.

3.4 Luminosity and Mass of the diffuse ionized gas

Figure 4: Maps of the extinction-corrected mass estimated from the H$\alpha$ (top) and [O III] (bottom) emission, divided into Voronoi regions derived to reach a SNR$>3$ of the H$\beta$ line emission. The black thick contours in each map represent the SNR levels at 3, 5, 10, 30, and 50 of the [O III] line emission. The black cross marks the position of PDS 456.

In this section, our primary objective is to conduct a comprehensive examination of the properties of the large-scale ionized gas surrounding PDS 456. To achieve this, we exclude the influence of galaxy companions by utilizing circular apertures with a 5-pixel diameter (i.e., approximately the PSF) centered on their positions.

To estimate the gas properties, we utilize all available optical emission lines within Voronoi regions (see Venturi et al., 2018; Molina et al., 2022). These regions are selected to have a SNR of H$\beta$ larger than 3, and we use the python tool vorbin for this purpose (Cappellari & Copin, 2003). The H$\beta$ is crucial to estimate the extinction of the emission based on the H$\alpha$/H$\beta$ line ratio.

The maps of the ionized gas mass estimated from the extended H$\alpha$ and [O III] emission in Voronoi regions (i.e. $M_{\ce{H$\alpha$}}$ and $M_{[{O\textsc{III}}]}$, respectively) can be seen in Fig. 4. Within each Voronoi region, we extract the spectrum and perform a fit. The fitting procedure involves modeling all emission lines using individual Gaussian profiles with identical widths, except for the [O III] doublet transition. In a few cases, we may observe a broader width in [O III] compared to the other emission lines. This discrepancy is attributed to the fact that [O III] is the most reliable tracer for ionized winds and shows a higher SNR, along with H$\alpha$, which is more affected by stellar absorption features. The centroid position of each Gaussian is set to have the same redshift, with an adjustment margin equal to the MUSE spectral resolution. The widths obtained from the best-fit models are deconvolved for the MUSE spectral resolution. Adopting these constraints in the modeling helps to minimize the impact of noise, instrumental and sky features, as well as any absorption that might affect the shape of individual emission lines (see also Veilleux et al., 2023).

We measure the electron density ($n_{e}$), color excess ($\rm E(B-V)$), ionized gas luminosity ($L_{\rm{ion}}$) and mass ($M_{\rm{ion}}$) derived from the [O III] and H$\alpha$ emission lines. To estimate $n_{e}$, we utilize the doublet [S II] emission line ratio, as described in Osterbrock & Ferland (2006). The errors in $n_{e}$ are determined by applying a bootstrap algorithm for a sample of $n_{e}$ values obtained by randomly varying the intensities of the [S II] emission lines within 1$\sigma$ error. In cases in which the SNR of the [S II] doublet lines is less than 2.5, which tends to occur at larger distances from the quasar, we adopt a fixed value of $n_{e}=150\leavevmode\nobreak\ \rm cm^{-3}$, which represents the minimum value calculated through the [S II] line ratio in the outer regions. In this way, we find that $n_{e}$ is $\approx 600\leavevmode\nobreak\ \rm cm^{-3}$ at the center, decreasing to $\approx 150\leavevmode\nobreak\ \rm cm^{-3}$ in the outer regions. To correct for dust extinction in the observed flux, we estimate $\rm E(B-V)$ using the Balmer decrement H$\alpha$/H$\beta$ flux ratio. We apply the Calzetti et al. (2000) attenuation law, assuming an intrinsic Balmer ratio H$\alpha$/H$\beta$ of 2.86 for gas at an electron temperature of $T_{\rm{e}}=10^{4}\rm\leavevmode\nobreak\ K$ (Osterbrock & Ferland, 2006), and utilizing $\rm R_{V}=3.12$. Near PDS 456, we obtain $\rm E(B-V)$ ranging from 0.2 to 1.1$\leavevmode\nobreak\ \rm{mag}$, decreasing at larger distances. Then, we estimate the total luminosity, corrected for dust extinction, in H$\alpha$ and [O III] to be $log(L_{\ce{H$\alpha$}}/\rm{erg\leavevmode\nobreak\ s^{-1}})=42.47_{-0.31}^{+0.25}$ and $log(L_{[{O\textsc{III}}]}/\rm{erg\leavevmode\nobreak\ s^{-1}})=43.43_{-0.30}^{+0.26}$.

To estimate $M_{\rm{ion}}$, we use both the H$\alpha$ and [O III] transition luminosity, corrected for dust extinction, according to the following formulas::

$$M_{[{O\textsc{III}}]}=8\times 10^{7}\rm{M_{\odot}}\biggl{(}\frac{C}{10^{[O/H]-[O/H]_{\odot}}}\Biggr{)}\biggl{(}\frac{L_{[OIII]}}{10^{44}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}}}\Biggr{)}\biggl{(}\frac{n_{e}}{500\leavevmode\nobreak\ \rm{cm^{-3}}}\Biggr{)}^{-1}$$ (1)

by Osterbrock & Ferland (2006) (or in Carniani et al., 2015; Bischetti et al., 2017), where $C$ is the clumpiness of the gas (i.e., $\rm\langle n_{e}^{2}\rangle/\langle n_{e}\rangle^{2}$), $\rm[O/H]$ is the oxygen abundance relative to hydrogen, and

$$M_{\ce{H$\alpha$}}=3.3\times 10^{8}\rm{M_{\odot}}\biggl{(}\frac{L_{H\alpha}}{10^{43}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}}}\Biggr{)}\biggl{(}\frac{n_{e}}{100\leavevmode\nobreak\ \rm{cm^{-3}}}\Biggr{)}^{-1}$$ (2)

by Nesvadba et al. (2017) and Leung et al. (2019). In both cases, we assume a $T_{e}=10^{4}\leavevmode\nobreak\ \rm K$ and solar metallicities. We correct $M_{[{O\textsc{III}}]}$ by multiplying it by a factor of three, as determined by F17 due to the discrepancies between mass estimates obtained using equations 1 and 2. The total mass, as derived through the H$\alpha$ and [O III] emission lines, is $2.75_{-2.26}^{+2.96}\times 10^{7}\leavevmode\nobreak\ \rm M_{\odot}$ and $8.69_{-6.74}^{+9.41}\times 10^{7}\leavevmode\nobreak\ \rm M_{\odot}$, respectively. These are two-three orders of magnitude lower than the virial dynamical mass of $\sim 10^{10}\leavevmode\nobreak\ \rm M_{\odot}$ derived within 1.3 kpc by B19.

Figure 5: Total observed [O III] luminosity ($L_{[{O\textsc{III}}]}$) of PDS 456 (red star) compared with those measured for other samples of Type 1 quasars in the redshift range $z\sim$0.1 to $\sim$4. The quasars at $z<$1 are from Husemann et al. (2013) (green dots), Husemann et al. (2014) (yellow dots), Liu et al. (2014) (purple dots) and Perna et al. (2017) (orange dots). Quasars at $z>$1 are from Coatman et al. (2019) (blue dots) and Bischetti et al. (2017) (red dots).

In Fig. 5 the total $L_{[{O\textsc{III}}]}$ measured for PDS 456 is compared to those derived for other samples of Type 1 quasars over a large redshift range ($z\sim$0.1-4). PDS 456 clearly lies at the brightest end of the $L_{[{O\textsc{III}}]}$ distribution for sources at $z\lesssim 1$. Moreover, its $L_{[{O\textsc{III}}]}$ is also consistent with that observed for the bulk of the quasar population at Cosmic Noon (e.g., Coatman et al., 2019) supporting the role of PDS 456 as a local counterpart of high-z quasars (and it is smaller by a factor of tens only when compared to that of hyper-luminous quasars like WISSH ones).

3.5 The large-scale [O III] ionized outflow

Figure 6: Maps of $w_{80}$ (top left), $v_{50}$ (top right), $v_{10}$ (bottom left), and $v_{90}$ (bottom right), of the [O III] emission line. These maps provide an overview of the widths and velocity shifts from the entire [O III] profile. In each panel, the black thin contours represent the SNR levels shown in Fig. 4. The black thick contours in the center of the figure depict the molecular outflow, detected in B19. These contours correspond to the emission levels of 1, 4, 8, and 10 mJy. The white cross symbol marks the position of the quasar. The white lines in each map represent 2$\rm\sigma$ Gaussian-smoothed contours that enclose the regions where $w_{80}$ exceeds 400$\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$.

In this section, we aim at revealing the presence of an ionized outflow based on the kinematics of the gas. To achieve an adequate SNR in the spectra, we bin the final nuclear- and continuum-subtracted datacube into $3\times 3$ pixel regions, equivalent to $1.87\leavevmode\nobreak\ \times 1.87\leavevmode\nobreak\ \rm kpc^{2}$. Each spectrum extracted from this binned datacube is analyzed to model only the [O III] line emission, which is used as best proxy for the outflow. We adopt a single or double Gaussian model according with the Bayesian Information Criterion (BIC; Schwarz, 1978). From the best-fit models we derive the velocities maps: $v_{50}$, $v_{90}$, and $v_{10}$, representing the velocities containing a specific percentage of the total integrated flux of the emission line. Additionally, we calculate $w_{80}$, defined as $v_{90}-v_{10}$ (see Harrison et al., 2014). Commonly, the latter is used to trace kinematic properties associated with outflows (Vega Beltrán et al., 2001; Collet et al., 2016; Harrison et al., 2016).

In Fig. 6, the $w_{80}$ map reveals velocities exceeding $500\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ in the vicinity of the PDS 456 core. These values gradually decrease to $\approx 300\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ at larger radii and decrease further near the location of the companions. The $v_{50}$ map reveals negative velocities with $v_{50}<-100\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ to the east of the quasar’s position (Fig. 6). On the west side, velocities are in the range $|v_{50}|<50\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$, consistently with the maximum velocities observed in the CO (3-2) disk probed in ALMA data (B19). Velocities for both the blue and red wings are determined using $v_{10}$ and $v_{90}$, respectively. These measurements indicate an east-west velocity gradient, with the blue wings extending to $-450\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ in the east and the red wings reaching $250\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ in the west.

The kinematic maps highlight the following features: a notable velocity gradient in the east-west direction with negative to positive $v_{50}$ values, and perpendicular to the gradient of the compact molecular disk; a central $w_{80}$ exceeding $\rm\approx 500\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$, more than twice the central velocity dispersion of the molecular gas (i.e. $\sigma_{\rm{vel}}$$=w_{80}/1.09/2.355$ for a Gaussian profile); and $v_{10}<-300\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$, aligning with the extended CO (3-2) outflow represented by thick black contours in Fig. 6, which is also blueshifted. These findings collectively suggest the presence of outflowing gas, which is further supported by the analysis of NFM data presented in Sect. 4.2. The latter reveal an unmistakable blueshifted outflow detected in H$\alpha$ within $\approx$3 kpc from the center, aligned with the direction of the CO (3-2) outflow.

This information can be used for pointing out the region where the outflow dominates in WFM MUSE data, essential for establishing the maximum projected distance of the [O III] outflow. As a rough estimate, we calculate an average $\sigma_{\rm{vel}}$ of the H$\alpha$ outflow in NFM MUSE data of $\rm\sim 150\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ (see Fig.9). This corresponds to a $w_{80}=400\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$ and can be used as a reliable tracer for the [O III] outflow. Applying this criteria, we observe that the region where the ionized outflow dominates has a maximum projected size of $\approx 20\leavevmode\nobreak\ \rm kpc$, as shown by the white contours in each map in Fig. 6.

In conclusion, our findings reveal compelling evidence of an outflow with a potential maximum projected distance of 20 kpc. This ionized outflow is characterized by $w_{80}\rm\geq 500\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ in the center, and also exhibits blueshifted and redshifted velocities up to $\rm-450\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ and $\rm 250\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ to the east and west of the quasar, respectively.

3.6 No Rotational pattern in the ionized disk

ALMA data on PDS 456 reveal a compact (effective radius of 1.3 kpc) and nearly face-on (i$\sim$25 deg) molecular disk in CO (3-2) with a peak-to-peak velocity of $\rm 100\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ along the NS direction (see Fig. 4 in B19). Under the assumption that the molecular and ionized phases of the gas disk share the same velocity gradient (Levy et al., 2018), the kinematics of the ionized gas, probed with both [O III] and H$\alpha$, do not reveal an ionized counterpart to this disk. Instead, there is a gradient in velocities along the E-W direction (see the $v_{50}$ map in Fig. 6), perpendicular to the velocity gradient of the molecular disk and indicative of the outflow kinematics. Therefore, we suggest that the ionized disk is not detected, likely due to the rotational velocity of the disk being comparable to the MUSE spectral resolution, and the kinematics being dominated by the outflow.

4 NFM observations

The AO-NFM MUSE data provide us an unprecedented zoom-in of the central $\sim 24\times 24\leavevmode\nobreak\ \rm kpc^{2}$ region of PDS 456 with a high spatial resolution ($\sim 280\leavevmode\nobreak\ \rm pc$), reaching a surface brightness limit of $\approx 3\times 10^{-19}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ cm^{-2}\leavevmode\nobreak\ arcsec^{-2}}$ ($3\sigma$). This allows us to study in detail the morphology and kinematics of the ionized gas on the same scales of the molecular outflow (B19).

Figure 7: MUSE NFM Surface Brightness map of the H$\alpha$ emission. The black contours show the SNR levels of the H$\alpha$ emission of 3, 5, and 10. Grey contours show the morphology of the blueshifted CO (3-2) molecular outflow (B19), at the emission levels of 1, 4, 8 and 10 mJy. The black circles mark regions affected by significant residuals in the PSF subtraction due to the quasar and a nearby continuum source.

4.1 PSF and Continuum Subtraction

We use CubExtractor tools to subtract the PSF contribution of the quasar (including emission non-resolved from host-galaxy or NLR), and the continuum of the sources in the datacube. In particular, the CubePSFSub tool uses an empirical approach for the PSF subtraction at each wavelength based on pseudo-broad band images produced from the cube itself as described in detail in Cantalupo et al. (2019). This tool is specifically designed to extract faint and extended line emission on large scales from the bright continuum PSF of an isolated quasar assuming that it is the dominant source of continuum emission in the PSF area as it is in our NFM observations (e.g. Borisova et al., 2016; Arrigoni Battaia et al., 2019; Farina et al., 2019; Cantalupo et al., 2019; Travascio et al., 2020).

As anticipated in Sect. 2, the [O III] and H$\beta$ transitions are not covered as they fall in the Na ID notch filter from the laser guide star, that contaminate that spectral region. Hence, we perform the PSF and continuum subtraction on a datacube in which we mask the spectral regions including H$\alpha$, [N II] doublet and [S II] doublet emission lines. The analysis to extract with CubExtractor the extended emission in AO-NFM MUSE data is the same used for WFM MUSE data described in Sect. 3.2. We extract the extended emission (in the H$\alpha$, [N II], and [S II] transitions) by setting, in the CubEx function, a SNR threshold of 3 and a minimum number of connected voxels of 3000.

4.2 High-Resolution view of the Ionized versus Molecular Spatially-Resolved Outflows

The morphology of the H$\alpha$ emission shown in Fig. 7 (obtained from the AO-NFM observations) extends in a north-east direction with respect to the quasar, along the same direction of both the large-scale ionized and extended molecular outflows. The emission structure shows two peaks emission connected by a diffuse extended component. One peak is located close to the nucleus, while the other is situated 2 kpc east of the former, exhibiting a thick shell-like geometry. The morphology of this structure closely overlaps with that of the molecular outflow (B19; grey contours), suggesting a multi-phase composition of the outflowing gas. Due to the lower sensitivity of AO-NFM MUSE data, compared to that in the WFM MUSE data presented in the previous sections, we do not detect a similar elongated structure in the south direction.

Figure 8: First (top) and Second (bottom) Moment of the flux distribution of the H$\alpha$ emission obtained with the Cube2Im function in CubExtractor, representing the $v_{\rm{shift}}$ and $\sigma_{\rm{vel}}$ maps. Gray contours show the morphology of the extended (up to $\approx 5\rm\leavevmode\nobreak\ kpc$), blueshifted CO (3-2) molecular outflow, detected by B19, at the emission levels of 1, 4, 8 and 10 mJy. The black circles and contours are the same as in Fig. 7.

Fig.8 shows the maps of the 1st (top panel; $v_{\rm{shift}}$) and 2nd (bottom panel; $\sigma_{\rm{vel}}$) moment of the H$\alpha$ flux distribution obtained with the Cube2Im tool (see Borisova et al., 2016). The thick gray contours overlaid on the maps represent the morphology of the molecular outflow detected by B19. In the $\approx 1.4\rm\leavevmode\nobreak\ kpc$ thick shell-like region traced by the H$\alpha$ emission, the $v_{\rm{shift}}$ values are lower than $-150\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$ and $\sigma_{\rm{vel}}$ values $\geq 120\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$. These velocities differ from those present in the remaining regions, showing $v_{\rm{shift}}$ between $-50$ and $50\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$ and $\sigma_{\rm{vel}}$ around $<100\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$. The spatial correlation between this shell and the molecular and large-scale ionized outflows suggests that this perturbed emission, up to $\sigma_{\rm{vel}}$ of $200\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$, is likely tracing an outflowing gas front.

Figure 9: Spectra extracted from nuclear- and continuum-subtracted MUSE WFM (red; H$\alpha$) and NFM (black; H$\alpha$) datacubes, and from the ALMA data (green; CO (3-2)), in the region where the H$\alpha$ outflow is observed in NFM MUSE data.

Fig.9 displays the spectrum of the H$\alpha$ transition, zoomed in on the black H$\alpha$ transition, extracted from the thick ”shell”-like region with $\sigma_{\rm{vel}}$ higher than 90$\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ (maximum $\sigma_{\rm{vel}}$ in the center of the compact molecular disk). The spectrum of the H$\alpha$ line emission (in red) from the WFM MUSE data, extracted from the same region, and the spectrum of the CO (3-2) outflow (in green) from B19, are overlaid. The H$\alpha$-NFM and the CO (3-2) spectrum are normalized at the peak of the emission line, while the H$\alpha$-WFM spectrum is normalized to match its blueshifted wing with the H$\alpha$-NFM. These spectra are plotted as a function of $v_{\rm{shift}}$ with respect to the wavelength of the relative transition at the systemic redshift of the quasar. By comparing the maximum velocities reached by the wings of the emission lines, we observe identical $v_{\rm{shift}}$ values for the H$\alpha$ line emission in NFM and WFM MUSE data. We conclude that the molecular and ionized outflows share similar kinematics. These similarities (i.e. extent, morphology, kinematics) strongly suggests the co-existence of different gas phases within these outflows. Moreover, the morphology of the emission of the CO (3-2) molecular and the H$\alpha$ ionized outflows may depend on the spatial variation of density, implying a clumpy structure of the gas (Baron & Netzer, 2019), although it could also be related to spatially varying emissivity or intrinsic dust absorption.

Figure 10: Top panel: [N II]/H$\alpha$ emission line ratio map, at $\rm SNR_{\ce{H$\alpha$}}>2.5$, whose pixel size is 5 times larger than the native pixel in NFM MUSE data, that corresponds to $\sim 1.3\rm\leavevmode\nobreak\ FWHM_{PSF}$. The black dot marks the position of the central quasar, and has the FWHM size of the $\rm PSF_{NFM}$. Bottom panel: map of the [O III]/H$\beta$ line ratio estimated through the [O III] and H$\beta$ Optimally-Extracted NB images extracted from the WFM MUSE data with CubExtractor. The black contours represent the SNR levels of the H$\alpha$ emission extracted from NFM MUSE data at 3, 5, and 10.

4.3 Ionization mechanisms of the outflowing gas

In this section, we discuss the ionization mechanism of the outflowing gas in H$\alpha$, based on the spatially-resolved [N II]/H$\alpha$ and [O III]/H$\beta$ line ratios measured in NFM and WFM MUSE data, respectively, as key components in the diagnostic BPT diagram.

Fig. 10 (top panel) shows a smoothed map of the [N II]/H$\alpha$ emission line ratio obtained from the respective optimally-extracted images derived with Cube2Im. Most of the emission associated with the outflow region exhibits a drop in the log([N II]/H$\alpha$) ($<-0.3$) with respect to the values in the entire map. Unfortunately, no other relevant emission lines fall in the spectral range covered by NFM MUSE data, while WFM MUSE data cannot provide line ratio estimates at the same spatial resolution of NFM data. The low [N II]/H$\alpha$ line ratio weakens the hypothesis of shock excitation, which is typically associated with an increase in the [N II]/H$\alpha$ and [S II]/H$\alpha$ line ratios in the regions with highest values of $\sigma_{\rm{vel}}$ (e.g. Leung et al., 2019). This could be an indicator of star-formation occurring in the outflowing gas (e.g., Maiolino et al., 2017) or AGN ionization followed by the recombination of hydrogen atoms in a photon-dominated scenario, depending on the [O III]/H$\beta$ line ratio (e.g., O’Dell et al., 2009). To discern among these envisaged scenarios, the bottom panel in Fig. 10 presents a map of the [O III]/H$\beta$ line ratio, derived from the optimally-extracted NB images obtained from WFM MUSE data with the CubExtractor tools by setting a SNR threshold of 3. We observe that the [O III]/H$\beta$ line ratio tends to increase in the region of the outflow detected in H$\alpha$ emission in NFM MUSE data.

This, combined with the [N II]/H$\alpha$ line ratio, agrees with results of Hinkle et al. (2019), who suggest that an outflowing gas showing a larger [O III]/H$\beta$ line ratio and a smaller low-ionization line ratio (i.e. [N II]/ H$\alpha$, [S II]/ H$\alpha$, [O I]/ H$\alpha$) than those of the host galaxy ISM implies a difference in the ionization parameter. In addition, Mingozzi et al. (2019) explain that outflowing gas in the NLR of Seyfert galaxies shows a low-[N II]/H$\alpha$ line ratio as it is directly illuminated by the AGN continuum. This may result from the H$\alpha$ emission tracing matter-bound clouds within the outflow, which dominate the emission in the ionization cone (see Mingozzi et al., 2019, for a discussion).

Figure 11: Spectra, in velocity units, extracted from the regions in which the outflow dominates the emission (see text for details), in WFM (top panel) and NFM (bottom panel) MUSE data. We report a zoom-in of the [O III] (top panel) and H$\alpha$ (bottom panel) transitions. The red lines show the best-fit model consisting of the combination of two Gaussian components: a narrow one (green) and a broader one (blue), tracing the outflow.

Table 2: List of the properties of the ionized (H$\alpha$, [O III]; top subtable) and molecular (CO (3-2); bottom subtable in B19) outflows in PDS 456. The numbers within the square brackets represent the minimum and maximum values within 1$\sigma$, respectively (see text for details). See text for details on the meaning of estimated properties.^a For the ionized outflow, the reported values include $R_{\rm out}$ for the region traced by the [O III] and $\Delta R$ for the region traced by the H$\alpha$ (see the text for details), while for the CO (3-2) outflow, the presented information includes the range of distances in which the molecular outflowing blobs are detected. ^b $v_{\rm out}$ for the the ionized and molecular outflows indicate the $v_{\rm max}$ (i.e. $v_{\rm{shift}}+2\sigma_{\rm vel}$) or the interval of velocities, respectively. The same ranges of radii and velocities for the molecular outflow are reported in Table 1 of B19. They provide ranges for the radius and velocity of the outflow, as the properties of both the total and the extended outflow were derived by combining the characteristics of individual blobs. We consider only the molecular blobs spatially connected to the H$\alpha$ outflow in the NFM MUSE observations when assessing the properties of the 1-3 kpc-scale molecular outflow.

Total Outflow (extended + unresolved)
	$L_{\rm out}$	$R_{\rm out}$ a	$v_{\rm out}$ b	$\langle n_{e}\rangle$	$M_{\rm out}$	$\dot{M}_{\rm out}$	$\dot{E}_{\rm kin}$	$\dot{P}_{\rm{out}}$
	$[10^{41}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}}]$	$[\rm{kpc}]$	$[\rm km\leavevmode\nobreak\ s^{-1}]$	$\rm[cm^{-3}]$	$[10^{7}\leavevmode\nobreak\ \rm M_{\odot}]$	$[\rm M_{\odot}\leavevmode\nobreak\ yr^{-1}]$	$[10^{42}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}}]$	$[10^{34}\leavevmode\nobreak\ \rm dyne]$
[O III] (WFM)	38 [35-42]	11.2 [6.5-11.4]	$593\leavevmode\nobreak\ [518-668]$	$200$	2.3 [2.1-2.5]	3.8 [3.1-4.9]	0.44 [0.35-0.54]	1.5 [1.2-1.7]
				410 [334-497]	1.1 [0.9-1.3]	1.84 [1.47-2.37]	0.21 [0.17-0.26]	0.71 [0.60-0.84]
CO (3-2)		$\sim$1.8-5, $<$1.2	$\rm[-1000,650]$		25[16-28]	290[180-760]	40[7-58]	120[78-320]
1-3 kpc-scale Outflow (without unresolved component)
H$\alpha$ (NFM)	8 [6-9]	1.4 [0.9-1.9]	690 [540-840]	500 [300-700]	0.47 [0.35-0.67]	2.5 [1.9 - 3.1]	0.38 [0.28-0.51]	1.1 [0.9-1.3]
CO (3-2)		1.8-5	$\rm[-1000,-250]$		5.1 [4.5-5.7]	42 [37-47]	7.3 [1.3-10.6]	19 [17-21]

5 Energetics of the Ionized Outflow in PDS 456

In this section, we estimate the properties ($M_{\rm out}$, $\dot{M}_{\rm{out}}$, $\dot{E}_{\rm kin}$ and $\dot{P}_{\rm out}$) of the ionized outflows in PDS 456 probed with H$\alpha$ emission in NFM MUSE data (up to $\approx 3\rm\leavevmode\nobreak\ kpc$ scales) and [O III] emission in the WFM MUSE data (up to $\approx 12\rm\leavevmode\nobreak\ kpc$ scales). First, we compare the properties of the extended ionized outflow probed with WFM MUSE data with those found in other AGNs. Most of them are reported in F17 and updated in B19, while others values can be found in Fluetsch et al. (2019) and Speranza et al. (2024). F17 assembled a sample of cold molecular outflows (CO, OH) and warm ionized outflows ([O III], H$\alpha$, H$\beta$) from AGN spanning a broad range of redshifts ($z<3$) and bolometric luminosities. Most of these AGN were studied using integral field unit (IFU) data. Subsequently, given that the morphology and kinematics of the 1-3 kpc-scale (excluding the blob B, as detected by B19, located 5 kpc south of PDS 456; see Fig. 7) blueshifted CO (3-2) and H$\alpha$ outflows observed in the MUSE-NFM data suggest that these components are part of the same multi-phase outflow and are likely driven by the same past AGN feedback event (Sect. 4.2), we proceed to compare the properties and energy of these two outflow phases. This is crucial to explore the distribution of energy in different outflow phases and its correlation with distinct mechanisms that determine energy efficiency coupling. We therefore calculate the integrated properties of the ionized outflow. To achieve this, we extract a spectrum from the WFM and NFM MUSE data in regions dominated by the outflow, defined as the region where $w_{80}\geq 400\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$ (see Sect. 3.5), and the region where $\sigma_{\rm{vel}}$ is larger than $90\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ (see Sect. 4.2), respectively. In both cases, we define the outflow based on a fitting model with two Gaussians: a narrow one ($\sigma_{\rm{vel}}$$\approx 76\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$) representing the systemic, non-outflowing gas, and a broader one ($\sigma_{\rm{vel}}$$\approx 240\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$) tracing the outflow component. From the best-fit model of the outflowing component, we derive some parameters for estimating additional outflow properties, and their associated 1$\sigma$ uncertainties are calculated using a bootstrap algorithm with $10^{4}$ iterations. This process involves propagating the errors obtained from the best-fit models of the initial spectra.

We derive the maximum velocity (i.e. $v_{\rm max}\equiv$ $v_{\rm{shift}}$ $+2\times$ $\sigma_{\rm{vel}}$; F17) and the observed integrated fluxes, that is corrected for dust extinction. The latter, characterized by $\rm E(B-V)\approx 0.52\leavevmode\nobreak\ \rm{mag}$, is estimated from the H$\alpha$/H$\beta$ line ratio observed in the central $5\times 5\rm\leavevmode\nobreak\ kpc^{2}$ region in the WFM MUSE data, with the same method used in Sect. 3.4. This value is consistent with the results obtained by Reeves et al. (2021) for PDS 456 and is also used for the extinction correction of the H$\alpha$ outflow.

The energetic properties that we derive for the outflow are dependent on $n_{e}$. Most of the studies of spatially-resolved kpc-scale outflows in AGN assume a uniform $n_{e}$ value, e.g. $<500\rm\leavevmode\nobreak\ cm^{-3}$ (Kakkad et al. 2016, Rupke et al. 2017, F17, Carniani et al. 2015, Bischetti et al. 2017) or $>1000\rm\leavevmode\nobreak\ cm^{-3}$ (Müller-Sánchez et al., 2011; Santoro et al., 2018; Förster Schreiber et al., 2019; Baron & Netzer, 2019). On the other hand, several studies show that $n_{e}$ of the outflows in AGNs can vary between $\sim 50\rm\leavevmode\nobreak\ cm^{-3}$ and $10^{6}\rm\leavevmode\nobreak\ cm^{-3}$ (Perna et al., 2017; Kakkad et al., 2018; Rose et al., 2018; Baron & Netzer, 2019; Holden et al., 2023; Venturi et al., 2023). In addition, uncertainties on $n_{e}$ can be also large due to the limited data quality and the strong assumptions adopted in each method (Baron & Netzer, 2019; Davies et al., 2020; Revalski et al., 2022).

Figure 12: $v_{\rm max}$ (top panel), $\dot{M}_{\rm{out}}$ (central panel) and $\dot{E}_{\rm kin}$ (bottom panel) of the molecular and ionized outflow in PDS 456 using star symbols, versus $L_{\rm bol}$. We compare these values with those of ionized and molecular outflows reported in F17, B19, Fluetsch et al. (2019) and Speranza et al. (2024). For the estimation of the ionized outflow’s properties, we use an $n_{e}$ of $200\leavevmode\nobreak\ \rm cm^{-3}$ from F17.

For a proper comparison of the [O III] outflow properties in PDS 456, as observed in WFM MUSE data, with those reported in F17, we use the same electron density value, $n_{e}=200\rm\leavevmode\nobreak\ cm^{-3}$, as employed in their study. The mass of the [O III] outflow is derived using the equations (1) and, based on F17, the outflow mass estimated with the [O III] line is systematically smaller by a factor of 3 than that derived with H$\alpha$ or H$\beta$ transitions, which are robust estimators of the mass. Therefore, we multiply the $M_{[{O\textsc{III}}]}$ in eq. (1) by a factor 3. To calculate the mass outflow rate of the [O III] outflow we use the formula in F17 assuming a cone geometry:

$$\dot{M}_{\rm out}=3\times\Biggl{(}\frac{M_{[OIII]}\times v_{\rm{max}}}{R_{\rm{out}}}\Biggr{)}$$ (3)

where $R_{\rm out}$ is the maximum radius up to which high velocity gas (i.e. $w_{80}>400\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$) is detected, following the approach used in the F17 sample. Specifically, to take into account the uncertainties in the estimates of the radius of the bulk of the outflow, we compute the 50th percentile (the median value), as well as the 16th and 84th percentiles of the distances, weighted by the fluxes. The main properties of the total outflow traced from the [O III] emission in WFM MUSE data are listed in Table 2 (top raw).

For the [O III] outflow detected in WFM MUSE data, we estimate $M_{\rm{out}}=(2.3\pm 0.2)\times 10^{7}\leavevmode\nobreak\ \rm M_{\odot}$, $\dot{M}_{\rm{out}}=3.6_{-0.8}^{+1.3}\leavevmode\nobreak\ \rm M_{\odot}\leavevmode\nobreak\ yr^{-1}$ and $\dot{E}_{\rm kin}=(4.2\pm 1)\times 10^{41}\leavevmode\nobreak\ \rm{erg\leavevmode\nobreak\ s^{-1}}$. In Fig. 12 we report the plots $L_{\rm bol}$ versus $v_{\rm max}$ (top panel), $\dot{M}_{\rm{out}}$ (central panel) and $\dot{E}_{\rm kin}$ (bottom panel). The ionized outflow in PDS 456 exhibits a $v_{\rm max}$ that is consistent with the scatter of the relation in F17, at fixed $L_{\rm bol}$. The $v_{\rm max}$ of the total molecular outflow in PDS 456 (B19) reaches values around $1000\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$, that is consistent with the values observed in the ionized outflow. Furthermore, it falls within the lower limit of the scatter observed in the population of molecular outflows. $\dot{M}_{\rm{out}}$ and $\dot{E}_{\rm kin}$ of the ionized outflow in PDS 456 are $\approx$2-4 orders of magnitude lower than those estimated in the sample of F17 at a similar $L_{\rm bol}$. The properties of the molecular outflow in PDS 456 are slightly below, at the lower end, the extrapolation of the high-$L_{\rm{bol}}$ points for the molecular outflow in F17. $\dot{M}_{\rm{out}}$ and $\dot{E}_{\rm kin}$ of the total ionized outflow observed in PDS 456 are consistent with those of the well-studied outflow in the quasar XID2028 at $z\sim 1.59$ by using high-resolution JWST data (Cresci et al., 2023; Veilleux et al., 2023). The ionized outflow in XID2028 exhibits larger velocities ($\approx 1100\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$) compared to that in PDS 456, but shows comparable extent ($\approx 10\leavevmode\nobreak\ \rm kpc$), $M_{\rm out}$ ($10^{7}-10^{8}\leavevmode\nobreak\ M_{\odot}$), and $\dot{M}_{\rm{out}}$ ($3-6\leavevmode\nobreak\ \rm M_{\odot}\leavevmode\nobreak\ yr^{-1}$). However, the morphology of the ionized outflow and radio lobes in XID2028 suggests that they are connected, possibly indicating a jet-driven outflow. On the contrary, the orientation and size of the ionized outflow and radio jet (up to 200 pc in size) in PDS 456 are uncorrelated, suggesting a radiative-driven outflow mechanism. On the other hand, if we consider the [S II] line ratio to estimate $n_{e}$ of the [O III] outflow (Osterbrock & Ferland, 2006), we obtain a value of approximately $400\leavevmode\nobreak\ \rm cm^{-3}$. Assuming this as the most accurate value, both $\dot{M}_{\rm{out}}$ and $\dot{E}_{\rm kin}$ must be rescaled by a factor of two compared to those reported in Fig. 12 (see Tab. 2).

Thanks to the superb resolution of the NFM MUSE we are able to directly compare the ionized outflow traced by the H$\alpha$, and the CO (3-2) molecular outflows discovered by B19, within $\approx$1-3 kpc from PDS 456 (see Fig. 7). For the H$\alpha$ outflow, we derive $n_{e}=(500\pm 200)\rm\leavevmode\nobreak\ cm^{-3}$ using the [S II] lines ratios (Osterbrock & Ferland, 2006; Jiang et al., 2019), which is consistent with $n_{e}$ derived for the [O III] outflow observed in WFM MUSE data. The mass of this H$\alpha$ outflow is estimated using the equation (2). To estimate the mass outflow rate, we adopt the formula from (Husemann et al., 2019):

$$\dot{M}_{\rm out}=\Biggl{(}\frac{v_{\rm{max}}}{100\leavevmode\nobreak\ \rm{km\leavevmode\nobreak\ s^{-1}}}\Biggr{)}\leavevmode\nobreak\ \Biggl{(}\frac{M_{H\alpha}}{10^{6}\leavevmode\nobreak\ \rm{M_{\odot}}}\Biggr{)}\leavevmode\nobreak\ \Biggl{(}\frac{100\leavevmode\nobreak\ \rm{pc}}{\Delta R}\Biggr{)}\rm\leavevmode\nobreak\ \rm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$$ (4)

which provides the mass rate within an outflowing shell of thickness $\Delta R$. The key properties of the H$\alpha$ ionized outflow in the NFM MUSE data are summarized in Table 2 (bottom row). $M_{\rm out}$, $\dot{M}_{\rm{out}}$ and $\dot{P}_{\rm out}$ for the H$\alpha$ ionized outflow detected in NFM MUSE data, are about 1 order of magnitude lower than those for the molecular outflow. The distinct gap between the properties carried by the ionized and molecular outflows is also evident for the total outflows as discussed in the previous section. These discrepancies suggest that the ionized component of the outflow makes a negligible contribution in terms of energy. This finding is inconsistent with the expected energy equivalence based on the high-$L_{\rm bol}$ relation presented in F17.

Menci et al. (2019) stressed that the mass outflow rate of outflows depends strongly on the direction of the outflow with respect to the plane of the disk. In particular, the mass rate of outflows is highest in the direction of the disk and decreases as the outflow moves away from the disk. This is due to the fact that the density of the gas in the disk is higher than that in the outflowing gas, and the mass outflow rate is proportional to the density of the gas. Based on this scenario, the observed difference in mass outflow rates between ionized and molecular outflows at the $\approx$1-3 kpc-scales, as well as the overall outflow, may be attributed to the different directions of these outflows in relation to the disk’s plane. In this scenario, the molecular outflow could be triggered by nuclear winds impacting the ISM gas, while the ionized outflow might result from the radiative pressure of AGN photons on NLR clouds. However, this hypothesis needs to be confirmed through further detailed analysis of the multi-phase outflow observed at the $\approx$1-3 kpc-scales.

6 Scenarios for the expansion of multi-phase galactic-scale outflows

We discuss here a variety of possible scenarios that could explain the properties of the multi-phase outflow in PDS 456. B19 suggested that in a high-$L_{\rm bol}$ regime the ionized outflow may represent a significant fraction of the outflow mass, reducing the discrepancy between the low momentum boost measured from the molecular phase alone with the expectations for an energy conserving expansion of the nuclear UFO. Nonetheless, the momentum boost of the multiphase, galaxy-scale outflow of $(\dot{P}_{\rm mol}+\dot{P}_{\rm ion})/\dot{P}_{\rm rad}\approx 0.37$ is not consistent with an energy-conserving (i.e. non-radiative) expansion powered by the quasar, which predicts values in the range 5-10 (Zubovas & King, 2012).

6.1 Intermittent AGN phases scenario

A possibility to reconcile the observations with an energy conserving expansion relies on the fact that the outflow might have been inflated during a previous AGN phase. Indeed, considering a velocity $\approx 300\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$, it takes the outflow $\sim$1 Myr to expand to the 0.3 kpc scale of the “compact” outflow, and $\sim$10 Myr to reach the 3 kpc size of the extended outflow. Even 1 Myr is longer than the expected typical duration of a single AGN episode, which is estimated by both statistical (Schawinski et al., 2015) and analytical arguments (King & Nixon, 2015) to be $<$0.1 Myr. So it appears likely that the outflows that are observed now were inflated by multiple previous AGN episodes. This scenario would also reconcile the kinetic power of the galactic-scale outflow in PDS (Fig. 12) with that expected from the F17 relations, which reflect the overall quasar population and, thus, the long-term average of $L_{\rm bol}$. The luminosity during a single AGN episode is not constant; rather, it probably decreases as the accretion disc is consumed, while feedback prevents more material from reaching the disc efficiently (King & Pringle, 2007). Zubovas & Nardini (2020) showed that in such cases, outflow properties tend to correlate better with the long-term average AGN luminosity than the instantaneous one. Given that the current luminosity of PDS 456 is $L_{\rm bol}\simeq 0.7\leavevmode\nobreak\ L_{\rm Edd}$ (Reeves et al., 2000; Nardini et al., 2015) it is very likely that the long-term average luminosity is much lower. The UFO revealed in X-rays, on the other hand, probably evolves on sub-year timescales and so closely tracks the present-day luminosity. Assuming that the nuclear wind always had kinetic power equal to $\sim 20\%\leavevmode\nobreak\ L_{\rm bol}$, as today, a long-term average luminosity 10-20 times lower than the present-day value would bring the outflow properties in line with the expectation of an energy-driven scenario. Such a low average luminosity is easy to achieve: it requires the product of AGN duty cycle over the past 1-10 Myr and the average Eddington factor during an episode to be of order 0.035-0.07. The AGN duty cycle is defined in Zubovas & Nardini (2020) as $\rm\delta_{\rm{AGN}}=t_{\rm{ep}}/t_{r}$, where $\rm t_{\rm{ep}}$ represents the duration of the AGN episode, which is a measure of how long the AGN remains in an active state, while $\rm t_{r}$ refers to the recurrence time scale, which determines the frequency of AGN activity. By definition, the average luminosity during an episode is at least a factor few below Eddington (i.e. $L_{\rm{AGN}}>0.01\times L_{\rm{Edd}}$), so the duty cycle has to be somewhat higher than the population average of $\sim$0.07 to achieve this. Moreover, Zubovas et al. (2022) suggested that AGN episodes may be clustered hierarchically in time, with longer phases lasting $10^{7}\leavevmode\nobreak\ \rm yr$ (as suggested by Hopkins et al., 2005) during which multiple AGN events, each lasting $\sim 10^{5}\leavevmode\nobreak\ \rm yr$, occur.

6.2 Radiation-pressure driven outflow scenario

Yet another possibility, as suggested by B19, is that radiation pressure on dust has triggered these outflows in PDS 456, implying a value of the momentum load of galactic-scale outflow around unity (Ishibashi & Fabian, 2014; Costa et al., 2018). This scenario might be supported by the detection of outflowing emission traced by blueshifted high-ionization emission lines in the mid-IR (i.e. $15.56\leavevmode\nobreak\ \rm\mu m$ [Ne III] and $14.32\leavevmode\nobreak\ \rm\mu m$ [Ne V]) with a velocity larger than those found for [O III] and H$\alpha$ and, therefore, related to an inner, highly-extincted ionized wind component (e.g., Spoon & Holt, 2009); or detection of dust emission co-spatial with the multi-phase outflow.

6.3 Star-formation and radio-jet driven scenario

Both an origin of the multi-phase outflow in PDS 456 due to star formation activity or one linked to the power of a radio jet seem to be unlikely. Specifically, B19 excluded a dominant starburst contribution to the outflow acceleration since the kinetic power of the molecular outflow is significantly larger than that expected for a wind triggered by a star formation activity with a $SFR\simeq 50-80\leavevmode\nobreak\ \rm M_{\odot}$ as observed in PDS 456.

Yang et al. (2021) detected a complex radio structure comprising a radio core, a radio jet, and a diffuse component. These radio components extend to distances of $\sim 370\leavevmode\nobreak\ pc$ from the quasar, with luminosities of $log(L_{\nu}\nu/\rm{erg\leavevmode\nobreak\ s^{-1}})=40$ ($\rm\nu=1.66\leavevmode\nobreak\ GHz$), and $log(L_{\nu}\nu/\rm{erg\leavevmode\nobreak\ s^{-1}})=39.3$ for the radio core and jet, respectively. A strong correlation between the galactic-scale outflows and the radio jets seems unlikely because the radio luminosity is more than two orders of magnitude lower than the kinetic power associated with the galactic-scale outflows, although based only on the luminosity at $\rm\nu=1.66\leavevmode\nobreak\ GHz$. Moreover, the scales and orientation of the radio emission (see Fig. 2 in Yang et al. (2021)) are uncorrelated with those of the ionized and molecular outflows (Fig. 6 and 8), suggesting that they are not related.

7 Summary and Conclusions

We carried out the analysis of the first VLT/MUSE WFM and AO-NFM observation of the nearby ($z=0.185$), luminous ($L_{\rm bol}\sim 10^{47}\rm erg\leavevmode\nobreak\ s^{-1}$) quasar PDS 456. These data provided us with the opportunity of mapping with unprecedented detail the ionized emission from a luminous quasar with a multi-scale, multi-phase outflow.

Specifically, with the WFM MUSE data, we were able to investigate the environment and a kpc-scale outflow. The AO-NFM MUSE data offered an unmatched spatial resolution of $\sim$280 pc, enabling us to accurately study in detail the morphology and kinematics of the ionized outflow and compare it with those of the molecular outflow, with equally resolution and at the same scales. Our analysis provided the following key results:

1. PDS 456 resides in a complex environment characterized by the diffuse emission of ionized gas extending up to a maximum projected distance of $\sim 46\leavevmode\nobreak\ \rm kpc$, which traces the CGM and the presence of multiple line-emitting companion galaxies within $\sim$30-40 kpc (See Fig. 2). We estimate the mass of the emitting ionized gas, through the [O III] and H$\alpha$ emission lines, $\sim$3 and $\sim$9$\times 10^{7}\leavevmode\nobreak\ \rm M_{\odot}$, respectively (see Sect. 3.4). These findings match with theoretical and observational expectations suggesting that hyper-luminous quasars live in over-dense regions in terms of companion galaxies (Wagg et al., 2012; Decarli et al., 2017; Trakhtenbrot et al., 2017; Bischetti et al., 2021; Nguyen et al., 2020; Perna et al., 2023) and reservoir of gas (Bowen et al., 2006; Prochaska & Hennawi, 2009; Farina et al., 2013).

2. We also discover the existence of an ionized outflow ($v_{max}\sim\rm 600\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$), mainly detected in [O III], and characterized by a maximum projected size of $\approx 20\leavevmode\nobreak\ \rm kpc$ and a velocity gradient along the east-west direction (see Sect. 4.2; Fig. 6).

3. $\dot{M}_{\rm{out}}$ and $\dot{E}_{\rm kin}$ of the [O III] ionized outflow in PDS 456 are 2-4 orders of magnitude less than those derived in objects with similar $L_{\rm bol}$ in literature (; ; Fluetsch et al., 2019; Speranza et al., 2024), while the measured value of $v_{\rm max}\approx 600\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ is consistent with those typically found in other quasars (see Fig. 12).

4. The total multiphase outflow expands in a regime not consistent with energy conservation, contrary to what might be expected given its scale (King & Pounds, 2015). In Sect. 6.1, we discuss in detail the intermitted AGN phase scenario, that could provide an explanation for this finding.

5. The spatially resolved H$\alpha$ emission in the AO-NFM MUSE data partially traces the same [O III] outflow detected with WFM MUSE data down to $\approx 1-3\leavevmode\nobreak\$kpc-scales from the quasar. The high resolution enables us to reveal a thick shell-like geometry (see Fig. 7).

6. The remarkable similarity in morphology, direction, and kinematics between the H$\alpha$ outflow detected in NFM MUSE data and the extended CO (3-2) molecular outflow strongly suggests that these two components belong to the same multi-phase outflow and could be driven from the same past AGN feedback episode (see Figs. 8 and 9).

7. The molecular phase dominates the kinetic power of the multi-phase, galaxy-scale outflow, with a $\dot{M}_{\rm mol}$ and a $\dot{P}_{\rm out}$ that is about 1 order of magnitude larger than that of the ionized outflow.

In conclusion, our analysis of the MUSE data unveils novel and important details on the complex interplay between the different phases of AGN-driven outflows extending up to $\sim 20\rm\leavevmode\nobreak\ kpc$ (i.e. interacting with the CGM), and the richness of circum-galactic environment of PDS 456. This may help in planning and performing future multi-band investigations of the properties of the luminous quasars shining at cosmic noon.

We highlight that deeper MUSE observations in both WFM and NFM modes would allow us to enhance and broaden our analyses. Specifically, these WFM observations could reveal a more extensive CGM structure associated with PDS 456 across all optical transitions. Data with an larger SNR may help in identifying a fainter wing of the [O III], suggesting the presence of a less dominant or more obscured outflow at higher velocities. In addition, deeper MUSE NFM observations would be crucial for investigating the presence of more extended quiescent and outflowing ionized gas. This is very important for comparing the differences in clumpiness between the ionized and molecular phases of the gas and improving our ability to produce spatially resolved $n_{e}$ maps of the outflow, thanks to the high SNR maps of [S II] emission lines.

Furthermore, JWST/NIRSpec-IFU data (filter G235H/F170LP) and MIRI-IFU data (filter F770W) have been recently acquired. These data are crucial to investigate both hot molecular and highly-extincted ionized gas, potentially revealing additional phases of the outflow in PDS 456 (Spoon & Holt, 2009; Bianchin et al., 2022). Additionally, the rest-frame near-IR $\rm H_{2}$ emission line can be used as a valuable tracer for shocks and contributing to a clearer understanding of the main expansion mechanism in the shell.

Appendix A Best-fit Na ID absorption transition

Fig. 13 shows the best-fit models of the Na ID absorption lines observed at the position of the K1, K2 and K3 sources in WFM MUSE data. As detailed in Sect. 3.3, we extract the spectra at the sources position in 2-pixel radius circles and we apply the fitting model to follow:

$$I(\lambda)=I_{\rm{em}}\times(1-C_{f}\times[1-e^{-\tau_{0}^{-(\lambda-\lambda_{K})^{2}/(\lambda_{K}b/c)}-2\tau e^{-(\lambda-\lambda_{H})^{2}/(\lambda_{H}b/c)}}])$$ (5)

where $I_{\rm{em}}$ is the normalization parameter, $\lambda_{K}$ and $\lambda_{H}$ are the wavelengths of the sodium at 5891$\leavevmode\nobreak\ \AA$ and 5896$\leavevmode\nobreak\ \AA$ respectively, $C_{f}$ is the covering factor, $\tau_{0}$ is the optical depth at the central $\lambda$, $b$ is the Doppler parameter, and $c$ is the speed of light. This model is widely used in literature (e.g., Rupke et al. 2005; Sato et al. 2009; Perna et al. 2020) and provides a physically-motivated description of the absorption features, which can be used to derive column density and mass of the absorbing gas (Travascio et al. in prep.). The redshift estimated in this way for the sources K1 and K3, are consistent with those derived through molecular transitions with ALMA by B19.

Figure 13: Spectra extracted from 2-pixel radius circles centered on the sources K1, K2, and K3, zoomed in on the Na ID absorption line. The red lines is the best-fit models of the Na ID absorption features.