The nearby extreme accretion and feedback system PDS 456: finding a complex radio-emitting nucleus

When a black hole accretes close to the Eddington limit, the astrophysical jet is often accompanied by radiatively driven, wide-aperture and mildly relativistic winds. Powerful winds can produce significant non-thermal radio emission via shocks. Among the nearby critical accretion quasars, PDS 456 has a very massive black hole (about one billion solar masses), shows a significant star-forming activity (about seventy solar masses per year) and hosts exceptionally energetic X-ray winds (power up to twenty per cent of the Eddington luminosity). To probe the radio activity in this extreme accretion and feedback system, we performed very-long-baseline interferometric (VLBI) observations of PDS 456 at 1.66 GHz with the European VLBI Network (EVN) and the enhanced Multi-Element Remotely Linked Interferometry Network (e-MERLIN). We find a rarely-seen complex radio-emitting nucleus consisting of a collimated jet and an extended non-thermal radio emission region. The diffuse emission region has a size of about 360 pc and a radio luminosity about three times higher than the nearby extreme starburst galaxy Arp 220. The powerful nuclear radio activity could result from either a relic jet with a peculiar geometry (nearly along the line of sight) or more likely from diffuse shocks formed naturally by the existing high-speed winds impacting on high-density star-forming regions.


INTRODUCTION
Accreting supermassive black holes (SMBHs) can provide mechanical feedback on their host galaxies via launching jets and winds. Jets are collimated relativistic outflows emitting synchrotron emission (e.g. Blandford, Meier & Readhead 2019). Winds are non-jetted, relatively wider angled outflows with a relatively low speed of < ∼ 0.3 c (e.g. Tombesi 2016) and are mainly traced by optical and X-ray spectroscopic as well as radio continuum observations (e.g. ⋆ E-mail: jun.yang@chalmers.se Panessa et al. 2019). When the SMBHs have accretion rates approaching the Eddington limit, they may produce not only jets (e.g. Giroletti et al. 2017;Yang et al. 2019Yang et al. , 2020a) but also extremely powerful winds (e.g. Nardini et al. 2015;Tombesi et al. 2015;Fiore et al. 2017). These winds can sweep out the interstellar medium (e.g. Zakamska & Greene 2014), produce radio-emitting shocks (e.g. Nims, Quataert & Faucher-Giguère 2015), contribute the extragalactic gamma-ray background (e.g. Lamastra et al. 2017) and quench star formation (e.g. Kormendy & Ho 2013;Fiore et al. 2017). Observing this complex nuclear radio activity with the very-long-baseline interferometric (VLBI) imaging technique can provide information to constrain physical properties in the extreme accreting environment and help understand feedback with the host galaxy.
As a critical accretion system, PDS 456 is of great interest for probing powerful multi-phase winds and outflows. It hosts a quasi-spherical mildly relativistic ( < ∼ 0.3 c) X-ray wind with a very high kinetic energy ∼0.2L Edd (Nardini et al. 2015). The power of X-ray winds positively correlates with the X-ray luminosity, and thus these winds are likely radiation pressure driven (Matzeu et al. 2017). From soft X-ray (Boissay-Malaquin et al. 2019;Reeves et al. 2020) to ultra-violet (O'Brien et al. 2005;Hamann et al. 2018) bands, there are also more reports of highly blue-shifted absorption lines at velocities consistent with the hard Xray winds. Moreover, the CO (3-2) emission line observations at ∼1 mm show that there exists not only some spatially extended molecular outflows up to ∼5 kpc but also high-velocity (∼800 km s −1 ) compact outflows in the nucleus (Bischetti et al. 2019).
Based on the multi-wavelength spectral energy distribution (SED), PDS 456 is a radio-quiet analogue of the well-known radio-loud blazar 3C 273 (Yun et al. 2004). The radio counterpart of PDS 456 is a steep-spectrum source (Yun et al. 2004) with a flux density of 25 ± 4 mJy at 1 GHz and a spectral index of α = −0.84 ± 0.11 (Yang et al. 2019). Besides the intense active galactic nucleus (AGN), its SED in the far IR shows evidence for significant star-forming activity in its host galaxy (Yun et al. 2004). The radio emission is likely dominated by the AGN activity (Yun et al. 2004). The pilot high-resolution VLBI observations at 5 GHz found a faint jet structure consisting of a few components on the deca-pc scale in the nuclear region, while failed to image some very diffuse radio structure (∼70 per cent of its total flux density, Yang et al. 2019). To get a complete view on the potential large scale radio activity which may owe to an episodic jet, and understand nuclear starbursts and/or wind shocks, we performed new VLBI imaging observations at a suitable frequency of 1.66 GHz with the European VLBI Network (EVN) plus the enhanced Multi-Element Remotely Linked Interferometry Network (e-MERLIN).
This paper is organised as follows. We describe the VLBI observations and the data reduction in Section 2 and present the imaging results of PDS 456 in Section 3. We discuss star-formation rate based on the IR SED, physical characteristics of the observed complex nuclear radio activity and potential identification with the jet, starburst and shock activity, and the accretion rate in Section 4. Throughout the paper, a standard ΛCDM cosmological model with The high-dynamical range map of its complex jet structure. The full width at half maximum (FWHM) of the synthesised beam is 29.6 mas × 6.6 mas at 73.9 • . The contours are 0.05 × (−1, 1, 2, ..., 64) mJy beam −1 . The image peak is 229 mJy beam −1 . With respect to the image noise level 0.017 mJy beam −1 , the dynamic range reaches 13 740. The inset plots the non-optimal (u, v) coverage. (b) The noise distribution in the residual map. The display range is between −0.083 to +0.070 mJy beam −1 . Both observations were carried out in the e-VLBI mode. The data transfer speeds were 1024 Mbps (16 sub-bands in dual polarisation, 16 MHz per sub-band, 2-bit quantisation) at the EVN stations and 512 Mbps (2 sub-bands in dual polarisation, 64 MHz per sub-band, 2-bit quantisation) at the e-MERLIN stations (Kn, Da, Pi, Cm and De). The data correlation was done in real time by the EVN software correlator (SFXC, Keimpema et al. 2015) at JIVE (Joint Institute for VLBI, ERIC) using an integration time of 1 s and a frequency resolution of 1 MHz.
The observations of the faint source PDS 456 applied the phase-referencing observing technique. The pcscale compact source J1724−1443 (Petrov et al. 2006), about 59 arcmin apart from our target, was used as the phase-referencing calibrator. Its J2000 position is RA = 17 h 24 m 46. s 96654 (σ ra = 0.1 mas), Dec. = −14 • 43 ′ 59. ′′ 7609 (σ dec = 0.2 mas) in the source catalogue 2016a from the Goddard Space Flight Centre (GSFC) VLBI group. The cycle time of the periodic nodding observations was about four minutes (∼0.5 min for J1724−1443, ∼2.5 min for PDS 456, ∼1.0 min for two gaps). A very bright flat-spectrum radio quasar NRAO 530 (J1733−1304, e.g. An et al. 2013) was also observed to determine the instrumental phases and bandpass shapes.
The data were calibrated using the National Radio Astronomy Observatory (NRAO) software package Astronomical Image Processing System (aips, Greisen 2003). First, we reviewed the data carefully and flagged out offsource data and some very low-sensitivity data. Second, we ran a-priori amplitude calibration with properly smoothed antenna monitoring data (system temperatures and gain curves) or nominal system equivalent flux densities when these monitoring data were not available. Third, we updated the position of the antenna Sr with the aips task clcor (options ANAX and ANTP) according to the first geodetic VLBI measurements (project code EL054, reported by Leonid Petrov on the web 1 ): J2000 epoch, position: (X, Y, Z) = (+4865183.542, +791922.255, +4035136.024) m, velocity: (V x , V y , V z ) = (−12.15, +19.29, +10.86) mm yr −1 and axis offset 0.031 m. The corrections are (∆X, ∆Y, ∆Z) = (+0.6116, −0.1553, −1.0446) m for the antenna position and +0.031 m for the antenna axis offset. Fourth, the general bootstrap phase calibrations were performed. We corrected the ionospheric dispersive delays according to the maps of total electron content provided by Global Positioning System (GPS) satellite observations, removed phase errors due to the antenna parallactic angle variations, aligned the phases across the subbands via running fringe-fitting with a two-minute scan of the NRAO 530 data. After these calibrations, we combined all the subbands, ran the fringe-fitting and applied the solutions to PDS 456 by interpolation. Finally, the bandpass calibration was performed. All the above steps were scripted in the ParselTongue interface (Kettenis et al. 2006).
We first imaged the phase-referencing source J1724−1443. The imaging procedure was performed through a number of iterations of model fitting with a group of delta functions, i.e., point source models, and the self-calibration in difmap (Shepherd, Pearson & Taylor 1994). This is similar to the non-negative least-squares algorithm (NNLS, Briggs 1995) in the image plane and thus can also allow users to achieve a high-dynamic-range map (Yang et al. 2020b). We re-ran the fringe-fitting and the amplitude and phase self-calibration in aips with the input source model made in difmap. All these solutions were also transferred to the data of PDS 456 by the linear interpolation.
The total intensity and residual maps from the direct model fitting are shown in Fig. 1. The calibrator J1724−1443 has a single-sided core-jet structure with a total flux density of 0.27 ± 0.02 Jy. Its radio core has a flux density of 0.23 ± 0.01 Jy and is used as the reference point in the phase-referencing calibration. There were 79 point source models used in the Stokes I map. Both the Stoke I and zeroflux density V maps have almost the same off-source noise level, σ rms = 0.017 mJy beam −1 . The noise distribution in the residual map of Stokes I is quite uniform in particular in the on-source region, although there is only one station (Hh) on the long baselines. Because the method tried to use the minimum number of point sources and the (u, v) coverage is poor, the faint jet structure does not looks very smooth. The e-MERLIN data were excluded in the final image because they gave some low-level noise peaks (∼0.15 mJy beam −1 , ∼0.07 per cent of the peak brightness) in the residual map.
We imaged the faint target PDS 456 without any self-calibration in difmap. First, we made a low-resolution clean map with the e-MERLIN and Jb. The target PDS 456 has a peak brightness of 7.6±0.2 mJy beam −1 using natural weighting (beam FWHM: 506 × 111 mas 2 in position angle 14 • ), and indicates a slightly resolved structure with a total flux density of 12.3±0.3 mJy and a de-convolved size of ∼123 mas. After the deconvolution of the major component with a proper window, there is a hint of two-sided jet-like structure at position angle (PA) of about −25 • and +145 • in the residual map but they are not bright enough (<5σ) for a confirmation as real features. This is followed up by making a high-resolution map with the EVN plus e-MERLIN data. We used σ −2 (σ, data error) as the visibility data weight and the purely natural grid weighting to clean the diffuse emission (peak brightness: 0.23 mJy beam −1 ) in the nuclear region. Owing to residual phase errors of the phase-referencing calibration, our final image has a noise level about a factor of two higher than the statistical value estimated in the zeroflux Stokes V map.

EVN PLUS E-MERLIN IMAGING RESULTS
The total intensity maps of PDS 456 observed with the EVN plus the e-MERLIN on 2020 January 22 are shown in Fig. 2. The beam synthesised with natural weighting has an FWHM of 23 × 5.2 mas 2 in PA = +73 • . To avoid resolving out diffuse radio features, we artificially increased the beam area and used a circular beam with a FWHM = 23 mas in the map. The nucleus of PDS 456 consists of two relatively compact knots, denoted as the components C and NW, and a diffuse emission region extending mainly toward West and South. According to the position of its centroid, we labelled the complex structure as the component W. The unshown image from the first epoch observation has a significantly lower quality owing to the limited telescopes and sensitivity, yet it independently confirms the existence of the components C, NW and W. The relatively bright spot at ∼12σ in the South, Fig. 2 was only marginally (∼3σ) detected in the early epoch. The circular Gaussian model fitting results in difmap (version 2.5e) are given in Table 2. All the errors are formal uncertainties measured at the reduced χ 2 = 1. Empirically, the flux density measurements have a systematic error of around ten per cent. Compared to the flux density (16 ± 2 mJy at 1.66 GHz) predicted by its radio spectrum (Yang et al. 2019), the VLBI map only restores ∼70 per cent. With respect to J1724−1443, the differential astrometry gives a J2000 position of RA = 17 h 28 m 19. s 78958 and Dec. = −14 • 15 ′ 55. ′′ 8520 for the peak component C. This is in agreement with the average position of components C1 and C2 at 5 GHz reported by (Yang et al. 2019) if we consider significant systematic error due to the extended source structure and the different (u, v) coverage, ∼0.1θ size (∼1 mas). With respect to the component C, we reported the relative offsets of components NW and W in Table 2.
In the last two columns of Table 2, we also report radio luminosity L R = νL ν (ν, observing frequency) and brightness temperature T b , estimated (e.g., Condon et al. 1982) as where S int is the integrated flux density in mJy, ν obs is the observation frequency in GHz, θ size is the FWHM of the circular Gaussian model in mas, and z is the redshift. The peak component C has a brightness temperature of ∼1 × 10 7 K. We use the Bayesian SED modelling and interpreting code BAYESED (Han & Han 2012 to decompose the IR SED by using an AGN torus model and a simple modified blackbody model to represent the contribution from cold dust emission heated by a young stellar population (for more details, see Fan et al. 2016Fan et al. , 2017Fan et al. , 2019. The derived IR luminosity of cold dust is ∼ 6.9 × 10 11 L ⊙ , which is lower than that of Yun et al. (2004) by a factor of about two. The difference is mainly due to the following two factors. First, our SED included the more data points at far IR wavelengths than that of Yun et al. (2004). Second, we excluded the contribution of the AGN torus emission in the IR luminosity estimate. We convert the derived IR luminosity of cold dust into a SFR of about 69 M ⊙ yr −1 , assuming a Chabrier initial mass function (Chabrier 2003). We note that this SFR estimate should be taken as an upper limit due to source blending in the IR bands. The SFR is also in agreement with the results reported by Bischetti et al. (2019), suggesting significant but no extreme star forming activity in the host galaxy.

Complex nuclear radio activity
The sub-kpc scale radio nucleus of PDS 456 shows a complex radio morphology. This is not unusual in ultraviolet to sub-mm luminous but radio weak/quiet sources because their radio emission originate from multiple physical mechanisms including low-radio-power jets, starbursts and shocks (e.g. Panessa et al. 2019). Owning to the absence of bright radio cores and the missing of sensitive short baselines, VLBI observations of the radio counterparts of these galaxies often show nondetections or significantly over-resolved radio structures. These cases involve extreme UV-selected starburst galaxies (Alexandroff et al. 2012), ULIRG IRAS 23365+3604 Table 2. Summary of characteristic parameters of the radio components detected in PDS 456. Columns give (1) component name, (2) total flux density, (3-4) relative offsets in Right Ascension and Declination with respect to the peak component C, (5) best-fit size θ size of the Circular Gaussian model, (6) brightness temperature T b and (7) radio luminosity L R .  (Frey et al. 2016) and submillimetre-selected galaxies (Chen et al. 2020).

Location of the radio core
The

Evidence of a low-radio-power jet
There exists some faint emission smoothly connecting the components C and NW. The component NW could represent the head of a jet along PA ∼ −52. • 1. With respect to the component C, the apparent jet opening angle for the component NW is ∼10. • 7. This is about a factor of two smaller than the median value of 21. • 5 found in significantly beamed AGN jets (Pushkarev et al. 2017). Thus, the jet is unlikely to be very close to the line of sight. Assuming that the jet is close to the kinematic axis of the molecular disk (Bischetti et al. 2019), the jet viewing angle is θ view = 25 • ± 10 • . This indicates a de-projected length of 570 ± 20 pc for the component NW. Since there is no highly relativistic jet observed among extreme accretion systems (Yang et al. 2020a), the jet in PDS 456 might not be significantly Doppler boosted. Based on its low radio luminosity, it can be identified as a low radio power jet (e.g. Kunert-Bajraszewska et al. 2010;An & Baan 2012).
The component NW has also been marginally detected as the component X2 in the previous EVN observations at 5 GHz (Yang et al. 2019). The two observations give a spectral index of α = −0.8 ± 0.1. Owing to the near face-on disk geometry (Bischetti et al. 2019) and the southern diffuse emission, only the jet on the approaching side is clearly identified. In the VLBI image at 1.6 GHz, we failed to confirm the existence of the component X1 (relative offsets, δ ra ≈ 37.4 mas, δ dec ≈ −99.4 mas), which was tentatively detected based on the shortest baseline of Ef-Wb on the receding side in the early observations at 5 GHz. Thus, the component X1 is probably an instrumental noise peak.

Diffuse component W: a relic jet or a composite of starbursts and shocks
PDS 456 has an optically thin integrated spectrum up to >15 GHz (Yang et al. 2019). As the component W contributes about half of the total radio emission of PDS 456 and the components C has a relatively hard spectrum, the component W may have a similar optically thin spectrum originating from non-thermal radio activity. It is hard to naturally identify the component W as a jet component extending further from the component NW, as W has a diffuse radio morphology and a rather different extension direction from the existing jet, and there is no hint of a bending point. With respect to the component C, the component W shows a very wide-angle (>180 • ) extension with a size of 118 ± 2 mas (366 ± 6 pc) and a centroid in PA ∼ − 103 • .
The component W may be an extended relic jet component close to the line of sight. Due to the peculiar jet geometry, the projected structure does not resemble an elongated jet. Compared to the inferred jet component NW, the identification necessitates a significant change in the jet direction. Currently, this information is unavailable.
Because of its high luminosity and extended morphology, it can not be simply explained as a single supernova or a supernova remnant, which have typical peak monochromatic luminosities only up to L R ∼10 38.7 erg s −1 (Weiler et al. 2002) and might only explain the brightened spot in the South. It is also hard to explain the component W as a group of supernovae and supernova remnants, i.e. starbursts. The nearby extreme starburst galaxy Arp 220 has a starformation rate (SFR) of ∼220 M ⊙ yr −1 and a radio luminosity of 4 × 10 39 erg s −1 (Varenius et al. 2016(Varenius et al. , 2019. Compared to Arp 220, PDS 456 has at least three times lower SFR (cf. Sect. 4.1). Assuming that fifty per cent of radio emission is from the starburst activity, the SFR derived based on the radio luminosity (Eq. 6, Bell 2003) would be 690 ± 70 M ⊙ yr −1 . This is inconsistent with the SED-based SFR and thus requires an additional source of energy injection.
Powerful high-speed AGN winds and outflows can naturally produce the additional synchrotron radio emission via shocks on scales > ∼ 100 pc in the nuclear highdensity star-forming region (e.g. Zakamska & Greene 2014;Nims, Quataert & Faucher-Giguère 2015). Mildly relativistic and wide-opening-angle winds have been identified in X-ray spectroscopic observations (e.g. Nardini et al. 2015;Matzeu et al. 2017;Reeves et al. 2020), with reports of multiple velocity components, up to 0.46 c . A possible ultraviolet outflow at 0.3 c has also been reported by Hamann et al. (2018). Thus, mildly relativistic winds can extend farther out to pc scale (Reeves et al. 2020), as also inferred from Section 4.3 and thus might have a large impact on the whole nuclear region. Two CO(3-2) molecular outflow components were identified by Bischetti et al. (2019) in the compact nuclear region. One is a blue-shifted, highvelocity (v ∼ −800 km s −1 outflow component. The other is a low-velocity ( < ∼ 500 km s −1 ), high-velocity dispersion (peak: 360 km s −1 ) outflow component, which is ∼50 mas west from the quasar and spatially coincident with the centroid of the radio component W. The coincidence agrees with the shock model.
The latter scenario involving shocks due to wind interaction is now probed further to estimate the SFR. This involves driving of galactic winds and consequent feedback with the host galaxy through star formation. The relevant forces in action on cold gas clouds (temperature ≤ 10 4 K) at pc-kpc scales include ram pressure from hot outflowing gas (in which the cold clouds are embedded) and radiation pressure from the galactic disk which are in opposition to the gravitational force due to the self-gravitating region of the disk (e.g. Sharma & Nath 2012); the corresponding SFR is approximated in terms of the wind velocity v w and the circular velocity of the galactic disk v c For the above estimated v w of 360-800 km s −1 and a fiducial value of v c = 120 km s −1 , the SFR is 13.9-102.1 M ⊙ yr −1 . As this range is consistent with the estimated 69 M ⊙ yr −1 from the IR observations, the latter is preferred. If the component W is associated with the nuclear molecular outflows, and formed by the AGN winds sweeping up the high-density star-forming region, the energy conversion efficiency is L R L bol ∼10 −7 . This is two orders of magnitude lower than predicted by the model of spherical wind shocks (Nims, Quataert & Faucher-Giguère 2015). However, this might be acceptable for PDS 456 as the nuclear winds can preferentially propagate along the low-density polar directions. Our estimate is also consistent with the study of Zakamska & Greene (2014) which finds a strong association between powerful outflows and the radio luminosity in radio quiet quasars, with a statistical median conversion efficiency of L R L bol ∼ 10 −6 .

Wind velocity and the Eddington ratio
If the powerful wide-angle X-ray winds (Nardini et al. 2015) are radiatively driven (Matzeu et al. 2017), PDS 456 might have a super-Eddington accretion rate. The Eddington ratio Γ e is the ratio of the accretion disk luminosity (L disk ) to the Eddington luminosity (L Edd ). For an optically thick geometrically thick disk (Shakura & Sunyaev 1973) in the vicinity of the SMBH, where M bh is the SMBH mass, M is the accretion rate and m is the accretion rate scaled in terms of the Eddington rate M Edd , L Edd = ǫ M Edd c 2 where ǫ is the efficiency factor, r bh = 6 r G = 6 GM bh /c 2 (r G is the gravitational radius) is the assumed distance from which the radiation is assumed to peak and corresponds to the innermost stable circular orbit for a Schwarzschild black hole (non-rotating). The evolution of velocity of a free particle driven by radiation pressure from the inner disk is given by ( where x = r/r G , ξ = (1 − 1/x), ξ s = ξ(x = 6) and we have assumed a radial particle trajectory which is shown to be the case for the large distance regime (Mohan & Mangalam 2015). The formulation corresponds to a wide-angle wind outflow launched from the innermost region subject to a radiative deceleration (an effective drag force), aided by the gravitational potential of the black hole. This can be used to estimate the Eddington ratio corresponding to an observed asymptotic wind velocity at the sub-pc -pc-scale. The equation (4) is solved assuming ǫ = 0.1 and Γ e = 0.1 − 10.0 (sub -super Eddington luminosity) which corresponds to m = 0.12 − 12.0 using equation (3). For a Keplerial angular velocity Ω bh = GM bh /r 3 bh associated with the launch radius r bh , the rotational velocity v φ = r bh Ω bh = c/6 1/2 ≈ 0.41 c. This could represent the minimum velocity of disk material that is advected into the winds and possibly into the collimated jet; for the boundary condition associated with equation (4), we assume an initial launch velocity β(r bh ) = 0.41. A launch velocity exceeding 0.41 c could indicate a poloidal component away from the disk plane and hint at mechanisms driving the material into trajectories away from the canonical Keplerian disk orbits, including a truncated disk scenario with a magnetically arrested disk accretion in the innermost region (e.g. Tchekhovskoy, Narayan & McKinney 2011) and energy extraction from the spin of the black hole (Blandford & Znajek 1977;Blandford & Payne 1982).
The above formulation can be generalized to estimate the velocity (or Lorentz factor) in the launching region for sources which indicate accretion rates near or exceeding the Eddington rate. The contours of the asymptotic velocity β(Γ e , x) are plotted in Fig. 3. The particle velocity tends to a constant by ≈ 1000 r G ≈ 0.1 pc (in the AGN rest frame). With an increase in Γ e , the velocity saturates to a larger value for the assumed β(r bh ). Wind velocities of ≤ 0.3 c, as suggested by the radio and X-ray observations can originate from the inner disk. These can be achieved well in advance and stay at the saturated value up to the pc-scale and possibly beyond. The Eddington ratio corresponding to β ≤ 0.3 c is Γ e ≤ 2.6 ( M ≤ 71.3 M ⊙ yr −1 ), indicating that PDS 456 likely accretes at rates exceeding the Eddington rate. The estimated Γ e ≤ 2.6 suggests that the bolometric luminosity (Reeves et al. 2000) may have been slightly underestimated (upto 3 × 10 47 erg s −1 ) or the SMBH mass (Nardini et al. 2015) may have been overestimated (M bh as low as 3 × 10 8 M ⊙ ).
The formulation can act an independent manner of constraining the Eddington ratio (and the associated accretion rate) in this and similar super-Eddington sources using the observed wind velocity as an input.  (4) as a function of the Eddington ratio and the distance from the SMBH r/r G (up to 1000, about 0.1 pc). The β values tend to a constant before about 0.1 pc for an initial launch velocity of 0.41 c. The radiatively driven and mildly relativistic winds with β = 0.1-0.3 c indicate the Eddington ratio Γ = 0.2-2.6.

CONCLUSIONS
The luminous radio-quiet quasar PDS 456 is a nearby extreme accretion and feedback system. With the EVN plus e-MERLIN observations at 1.66 GHz, we achieved a more complete view of its complex radio morphology. Our VLBI images revealed a complex radio-emitting nucleus consisting of a collimated jet and a very extended structure. The latter has a size of about 120 mas and a radio luminosity about three times higher than the nearby extreme starburst galaxy Arp 220. The powerful nuclear radio activity could be explained as a relic jet with an unusual jet geometry, or a more natural consequence of star-forming activity plus AGN-driven wind shocks. In addition, fainter nuclear radio activity and a sub-mJy flat-spectrum jet base may exist. The observed relativistic X-ray winds of ∼0.3 c in the inner sub-pc-scale region suggests that PDS 456 likely accretes exceeding the Eddington rate. Future multi-frequency deep VLBI observations of PDS 456 can potentially provide more details on the complex nuclear activity.