Jet-driven Galaxy-scale Gas Outflows in the Hyperluminous Quasar 3C 273

MUSE data of 3C 273 were obtained under ESO program 097.B-0080(A) on 2016 March 31 under gray moon, clear sky conditions and a seeing of 07. The data cover a 1' × 1' field-of-view (FoV) at 02 sampling and a wavelength range of 4750–9300 Å at a spectral resolution of R ∼ 2500. Given the brightness of 3C 273 (m_V = 12.6 mag), we split the observations into 18 × 250 s exposures with a total integration time of 4500 s. Small subarcsecond dithering and standard 90° rotations between exposures were applied to mitigate residuals of the flat-fielding. The data were fully reduced and calibrated using the standard ESO pipeline (Weilbacher et al. 2014).

The reconstructed gri MUSE color image is shown in Figure 1 (panel a). The bright AGN point source and the optical jet of 3C 273 are both clearly visible. An important processing step is therefore the deblending of the unresolved AGN emission from the spatially resolved emission of the host galaxy. We use the software qdeblend^3D (Husemann et al. 2013, 2014) to estimate the point-spread function (PSF) from the 2D intensity distribution of the Hδ, Hβ, He i, and Hα broad-line wings originating from the unresolved AGN broad-line region. Afterwards, we interpolate the wavelength-dependent PSF with a second-order polynomial on a spaxel-by-spaxel basis. The AGN spectrum from the central spaxels is then convolved with the PSF and subtracted from the original data. Because the Hα line is saturated at the peak of the PSF, we cannot use the central spaxel for the wavelength range around Hα. Hence, we instead use a region northeast of the AGN, which shows a deficit in extended emission to create a "pure" AGN spectrum around Hα, which we convolve again with the corresponding PSF and subtract it from the original data. Given the brightness of the AGN, we still obtain a very high signal-to-noise ratio (S/N) spectrum of the AGN a few arcseconds away from the center, which does not leave strong residuals after subtraction. Nevertheless, the spaxels with <08 around the AGN position remain corrupted in Hα due to the saturation and need to be ignored completely.

Zoom In Zoom Out Reset image size

As shown in Figure 1, we clearly recover spatially extended [O iii] emission with complex kinematics even at an AGN–host galaxy contrast ratio of about 200 in the continuum. In addition, we also recover the stellar continuum of the host galaxy at high S/N within an annulus of 08–6'' circular radii (Figure 2, top panel). We modeled the stellar continuum by fitting stellar templates from the INDO U.S. library (Valdes et al. 2004) using the PyParadise software (Walcher et al. 2015; Weaver et al. 2018). From the stellar continuum fit, we infer an accurate systemic redshift of z = 0.15850 ± 0.00005 and a stellar velocity dispersion within this larger aperture of σ_* = 210 ± 10 km s⁻¹.  The integrated spectrum reveals line ratios of [O iii λ5007]/Hβ = 3.9 ± 0.1, [N ii λ6583]/Hα = 0.64 ± 0.04, and He ii λ4685/Hβ = 0.11 ± 0.01, where He ii λ4685 is detected with 10σ confidence. We show the corresponding emission-line diagnostic diagrams in Figure 2 together with classical demarcation lines, which clearly show that star formation is not the dominant ionization mechanism of the ionized gas. Due to the complex kinematics and heavy blending of the [N ii] and Hα lines, we cannot robustly map this line ratio across the nebula. However, the [O iii λ5007]/Hβ line ratio map can be constructed (Figure 1, panel c) and reveals that the line ratio is even higher in most parts,  confirming that AGN and/or shock ionization is the main ionization mechanism across the nebula.

Zoom In Zoom Out Reset image size

ALMA uses 3C 273 as a calibrator so that a large number of short observations of this source are available. We retrieved all data available for 3C 273 covering the redshifted frequency of ¹²CO(1−0) with a spatial resolution between 05 and 10. This included eight different tracks from five different ALMA projects, all of which observed 3C 273 for ≈5 minutes in order to use it as a bandpass calibrator. In all but one case a 2 GHz correlator window (with a raw spectral resolution of 15.625 MHz) was present over the redshifted frequency of CO(1−0). In one case, a 1850 MHz correlator window was present, with a raw spectral resolution of 3.906 MHz. The other spectral windows present in each observation were used to detect continuum emission.

The raw ALMA data for each track were calibrated using the standard ALMA pipeline in the Common Astronomy Software Applications (CASA) package. Additional flagging was carried out where necessary to improve the data quality. Three iterations of self-calibration were used to improve the phase and amplitude calibration of the data on the bright continuum in 3C 273. The data presented here were produced using natural weighting, yielding a synthesized beam of 078 × 075 at a position angle of 37° (a physical resolution of ≈2.1 kpc). We use data with a channel width of 30 km s⁻¹ and pixels of 02 (resulting in approximately 3.5 pixels across the synthesized beam).

Bright point source continuum emission with a flux of 12.9 ± 0.02 Jy was detected over the full line-free bandwidth. Based on the position of this point source, we are able to register the relative AGN positions of the MUSE and ALMA data with a precision of at least one pixel size (02). We subtracted the continuum emission from the data in the uv plane using the CASA task uvcontsub. The continuum-subtracted dirty cubes were cleaned in regions of source emission (identified interactively) to a threshold equal to 1.5 times the rms noise of the dirty channels. The clean components were then added back and reconvolved using a Gaussian beam of FWHM equal to that of the dirty beam. This produced the final, reduced, and fully calibrated ¹²CO(1−0) data cube of 3C 273, with an rms noise level of 0.18 mJy beam⁻¹ in each 30 km s⁻¹ channel.

CO(1−0) emission was detected from this source, with an integrated line width of ≈490 km s⁻¹ (FWHM) and an integrated intensity of 1.82 ± 0.02 Jy km s⁻¹. This corresponds to an H₂ mass of (1.7 ± 0.3) × 10⁹ M_⊙, assuming α_CO = 0.8 M_⊙ (K km s⁻¹ pc²)⁻¹ as typically used for compact luminous systems, such as starburst galaxies, submillimeter galaxies (SMGs) and quasars (e.g., Downes & Solomon 1998; Scoville et al. 2003; Xia et al. 2012). If we were instead to use a Galactic α_CO, then the derived H₂ mass in this system would increase by a factor of 5.5.

Zeroth-moment (integrated intensity), first-moment (mean velocity), and second-moment (velocity dispersion) maps of the detected line emission were created using a masked moment technique. A copy of the clean data cube was first Gaussian-smoothed spatially (with an FWHM equal to that of the synthesized beam), and then Hanning-smoothed in velocity. A three-dimensional mask was then defined by selecting all pixels above a fixed flux threshold of 1.5σ, adjusted to recover as much flux as possible in the moment maps while minimizing the noise. The moment maps were then created using the unsmoothed cubes within the masked regions only. In Figure 1, we show the zeroth- and first-moment maps in panel (d) and (e), respectively. The morphology of the molecular gas is surprisingly asymmetric and appears as an "arc"-like feature about ∼18 away from the nucleus toward the southwest, almost exactly perpendicular to the jet axis. It seems consistent with the dust lane identified in coronagraphic images presented by Martel et al. (2003), but the significantly lower spatial resolution of the ALMA data does not allow a precise spatial comparison. Importantly, the molecular arc is significantly more extended along the southwest direction on both sides of the spiral continuum excess structure seen in the Hubble images. In particular, there is no corresponding feature in Hubble images at the location of the brightest molecular gas knot seen 2'' west of the nucleus. Because the ionized gas morphology is completely different and does not show a similarly prominent structure, it is unclear whether the molecular gas is bound to a stellar overdensity formed in a recent galaxy interaction or is currently forming in situ due to a density wave and enhanced ambient pressure initiated by the expanding outflow.

With MUSE, we recover extended [O iii] emission out to ∼20 kpc from the quasar. The [O iii] emission is asymmetric toward the east side as shown in Figure 1, panel (b). The [O iii] surface brightness distribution confirms the results of the Fabry–Perot imaging (Hippelein et al. 1996), but MUSE achieves a much higher S/N and covers a larger velocity range. In the bottom panels of Figure 1, we present two pseudo-slit spectra (slit width of 08) for two position angles, one along the radio jet axis (PA −45° N to E) and one perpendicular to it (PA 45°). In both cases, we see highly complex ionized gas motions, with amplitudes of ±800 km s⁻¹ as well as several broad and multiply peaked lines over this large velocity range in the central 4 kpc.

Figure 3 highlights the complexity of the [O iii] line profiles, which are difficult to model with a fixed number of Gaussians. It is highly desired to dissect the [O iii] line profile into distinct components, which has been used to study individual physical conditions by correlating line ratios and kinematics (e.g., McElroy et al. 2015). We have tried this process but find that fitting [O iii] line profiles separately for each spaxel violates the spatial correlation function imposed by the PSF as shown in Figure 3. The [O iii] line shape is significantly changing from spaxel to spaxel, which has a size of 02, about a factor of 3 smaller than the FWHM of the seeing. While three Gaussians provide a good fit to the overall line profile in all cases, the spatial variation of the Gaussian parameters are unphysical as the PSF demands that individual spectral components only vary in flux on such small scales. A fully consistent solution would only be achieved when the line shapes are modeled directly in 3D, taking the PSF into account as a constraint on the spatial flux distribution. Such a method, which is able to model the entire ENLR structure at once with a large number of components, still needs to be developed. For practical reasons, we instead map nonparametric line-shape parameters such as the median velocity (v₅₀), the peak velocity (v_peak), the FWHM, and the width containing 80% of the line flux (W₈₀), which has been previously used to characterize ionized gas outflows in the [O iii] line around luminous AGNs (e.g., Liu et al. 2013b; Harrison et al. 2014, 2016; McElroy et al. 2015; Sun et al. 2017). Because the method requires that the cumulative line flux distribution is monotonically increasing, i.e., no negative flux values across the line shape, we still need to generate a noise-free model of the lines for each spaxel. Given the complex line shapes, we use a large spectral library of [O iii] doublet lines to model the line profiles as a nonnegative linear superposition, which is not restricting the model to a fixed number of Gaussians. The library is constructed with a velocity range of ±1200 km s⁻¹ at 10 km s⁻¹ sampling and four line dispersions (40 km s⁻¹, 80 km s⁻¹, 160 km s⁻¹, 320 km s⁻¹), which are convolved with the instrumental resolution as inferred by Bacon et al. (2017) at the [O iii] wavelength. During the fitting, we discard components with a peak flux density S/N < 3 to avoid fitting noise features. We use a Monte Carlo approach to estimate errors for each parameter by re-fitting the noise-modulated data 50 times. In Figure 1, panels (f), (g), (h), and (i), we show the resulting kinematics maps with a v₅₀ velocity error less than 50 km s⁻¹ and S/N > 5 for W₈₀ and FWHM.

Zoom In Zoom Out Reset image size

The kinematics maps inferred from [O iii] as shown in Figure 1 are complex and are therefore best interpreted in conjunction with the molecular gas distribution and kinematics. The molecular gas corresponds to kinematically quiescent gas with a small line dispersion (≤45 km s⁻¹) and a smooth velocity gradient from blue- to redshifted motion of ±160 km s⁻¹ around the systemic velocity of the stellar body. The ionized gas shows the same smooth velocity gradient and radial velocity amplitude in v₅₀ and v_peak of ±150 km s⁻¹ at the location of the molecular gas with a small line width both in W₈₀ and FWHM. Furthermore, the ionized gas at even larger distance from the AGN exhibits the same direction of motion as the molecular gas, which suggests that the large-scale kinematics are dominated by gravitationally driven rotational motion. These quiescent motions at radii >2'' define the global rotational pattern with redshifted motion on the east and blueshifted motion on the west side of the galaxy.

On the contrary, the ionized gas kinematics enclosed within the radius of the molecular gas arc is highly complex as demonstrated in Figure 3. Nevertheless, the light-weighted bulk motion of the gas shows a bipolar velocity gradient along the radio jet axis with kinematics predominately redshifted on the approaching jet side and blueshifted on the opposite side. While the FWHM remains narrow with <400 km s⁻¹ across most of the ENLR, it appears larger at the transition to the rotation-dominated region. This is likely a superposition effect due to the PSF smearing at the transition region. Likewise, the large W₈₀ > 800 km s⁻¹ width is not caused by individual clouds with high-velocity dispersion but rather the spatial blending of many narrow emission lines with a broad distribution in radial velocities. While this complexity limits our ability to understand the interaction with ambient medium in detail, a key observation is that the radio jet axis is clearly aligned with the bipolar velocity gradient in the ionized gas,  indicative of a strong coupling between the expanding radio jet and the ambient gas in 3C 273. Although it is strange that the approaching side of the jet appears to push the gas backwards, we will outline a potential geometrical scenario in Section 3.2 based on the directly inferred relative orientation of the jet and gas disk.

A lot of different definitions have been adopted for the ENLR size, such as an effective radius (e.g., Husemann et al. 2014; Kang & Woo 2018), the maximum detectable size (Bennert et al. 2002; Villar-Martín et al. 2010, 2018; Greene et al. 2012; Storchi-Bergmann et al. 2018), or the size until an intrinsic surface brightness limit corrected for cosmological dimming (e.g., Hainline et al. 2013, 2014; Liu et al. 2013a; Sun et al. 2017). By design, all these definitions are prone to different observational effects, and their application is limited to specific science goals. The maximum detectable [O iii]-emitting region size is strongly dependent on the depth of observations and whether the ionization mechanisms can be ambiguously attributed to AGN ionization. It is therefore not well suited for studying scaling relations between AGN luminosity and ENLR size. Additionally, the luminosity-weighted radius depends on whether the dominating NLR is included or excluded from the radius estimations.

Several studies therefore adopted a characteristic ENLR size using a fixed intrinsic [O iii] line flux surface brightness limit of 10⁻¹⁵ erg s⁻¹ cm⁻² arcsec⁻²(1 + z)⁻⁴ to allow a comparable ENLR size–luminosity relation to be constructed at various redshifts. For 3C 273, we determine a characteristic ENLR size of R_ENLR = 12.1 ± 0.2 kpc, which is robust given the good spatial resolution and the application of the PSF subtraction technique. In Figure 4 (left panel), we compare this ENLR to literature measurements. Previously, a break in the L_AGN–R_ENLR scaling relation has been reported to occur somewhere at ∼6–10 kpc (e.g., Hainline et al. 2013, 2014; Liu et al. 2013a). In order to mitigate the effect of beam smearing for their high-redshift targets, Hainline et al. (2014) applied a surface brightness profile fitting convolved with an approximate PSF to recover the intrinsic size of the ENLR. Considering that we were able to subtract the dominating AGN point source directly and that we can exclude a top-hat function as a radial ENLR light profile, we compare our inferred ENLR size for 3C 273 with the Sérsic profile measurements from Hainline et al. (2014). We find that the ENLR size of 3C 273 is 2.7σ larger than the proposed break radius and sample distribution reported by Hainline et al. (2014). This suggests that either the ENLR size–luminosity relation has no break radius at all, that the break in the ENLR size occurs at significantly larger radii,  or that the ENLR is generically asymmetric at the chosen surface brightness limit so that the sizes are systematically underestimated in randomly oriented long-slit observations. Given that 3C 273 is an unobscured AGN, the intrinsic ENLR size may be even larger compared to the obscured AGN used in previous studies because the unified model predicts that the ionization cone of unobscured AGNs should be more aligned toward the line-of-sight (LOS) than that of obscured AGNs. This would mean that statistically, a large correction factor needs to be applied to turn the measured projected ENLR sizes into intrinsic sizes. We can rule out the ionization cone being directly aligned with our LOS to explain the relatively round morphology, because the kinematic substructure and the known direction of the radio jet implies a significant inclination of the cone. In any case, those inclination effects would only strengthen our result as the intrinsic size of the ENLR would be even higher for 3C 273. Interestingly, we note that adopting a 3× higher intrinsic [O iii] surface brightness limit as proposed by Sun et al. (2018) would imply a size of R_NLR = 4.15 ± 0.1 kpc, which is a factor of 4 smaller than the corresponding L_AGN–R_ENLR scaling relation reported by the authors. Hence, the choice of the surface brightness limit to measure the size has a great impact on the signatures of an ENLR break radius. Additionally, unaccounted systematic effects in the relation of Sun et al. (2018) may contribute to this discrepancy, due to the use of broadband imaging that still requires spectroscopic confirmation.

Zoom In Zoom Out Reset image size

A break in the size–luminosity relation would necessarily imply a transition from an ionization-bounded to a matter-bounded region,  e.g., running out of gas, as the radiation field smoothly decreased as R⁻². Hence, a transition may occur when the gas density drops faster than the radiation field at the edge of galaxies (Netzer et al. 2004; Hainline et al. 2013). The stellar body of 3C 273 can be traced out to a radius of 7'' (19 kpc) with an effective radius of r_eff = 26 (Martel et al. 2003). Hence, the ENLR does not reach beyond the host galaxy and is well covered within 3 × r_eff of the stellar distribution. A break in the ENLR size–luminosity relation may not occur at the chosen surface brightness limit simply because of the large size of 3C 273's host galaxy. Alternatively, Dempsey & Zakamska (2018) proposed that the transition from being ionization- to being matter-bounded occurs because optically thick clouds lead to a steeper decline of the radiation field than R⁻², causing a break in the relation at small ENLR sizes. For 3C 273, we find the dense molecular gas to be  located inside the characteristic ENLR radius so that dense gas shielding and filtering of AGN radiation may be responsible for the asymmetry of the ENLR, possibly reducing the ENLR brightness along the southwest direction but not toward the east side. Another scenario relevant for 3C 273 is the cooling of gas from the hot halo gas. According to the chaotic cold accretion (CCA) scenario, ionized gas filaments and clouds are condensing out of the hot phase within ∼10 kpc (Gaspari et al. 2018, their Figure 5). Due to CCA, the transition from ionization- to matter-bounded ENLR clouds would be related to the cooling radius where a large portion of gas becomes neutral again to be ionized by the AGN radiation. While a fraction of the condensed gas will be funneled toward the SMBH, a residual part will settle onto a clumpy rotating structure within the central 10 kpc scale, which is in agreement with our observations.

Overall, none of the processes that have been proposed to reduce the ENLR break radius seems relevant for 3C 273. This is in agreement with our result that 3C 273 is inconsistent with the reported break radius in the ENLR size–luminosity relation. However, very large EELRs have been reported around luminous AGNs at low surface brightness levels (e.g., Villar-Martín et al. 2010; Husemann et al. 2011; Kreimeyer & Veilleux 2013; Villar-Martín et al. 2018), which are an order of magnitude below the intrinsic surface brightness assumed for the characteristic ENLR radius here. Usually, those EELR structures are in a matter-bounded regime far beyond the transition radius (e.g., Fu & Stockton 2007; Kreimeyer & Veilleux 2013; Storchi-Bergmann et al. 2018) and are therefore usually associated with gas density enhancements in the intergalactic medium due to tidal tails, ongoing galaxy interactions, or ram pressure from an expanding radio jet. For 3C 273, we do not detect any [O iii] emission beyond 23 kpc at our surface brightness detection limit of Σ_{[O iii]} > 4 × 10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻² at the 3σ level per 02 pixel. This suggests that 3C 273 has not been subject to a gas-rich merger very recently, and the expanding radio jet has not expelled gas beyond the 20 kpc scale.

The sizes of outflows have been measured by defining a kinematically disturbed region (KDR) either based on a fixed cut in line width, W₈₀ > 600 km s⁻¹ (Sun et al. 2017), or in relation to the stellar velocity dispersion, σ_{[O iii]}/σ_* > 1 (Kang & Woo 2018). For 3C 273, we measure σ_* = 210 ± 10 km s⁻¹, which is close to the expected value of σ_* ∼ 250 km s⁻¹ from the M_BH–σ_* relation (Kormendy & Ho 2013) adopting M_BH ∼ 5 × 10⁹ M_⊙ (Kaspi et al. 2000), so that both KDR definitions are comparable. Indeed, previous measurements reveal a tight relation between L_bol and R_KDR as shown in Figure 4. For 3C 273, we measure R_KDR = 6.1 ± 1.5 kpc, adopting W₈₀ > 600 km s⁻¹, which is significantly smaller than the extrapolated relations proposed by Sun et al. (2017), but close to those by Kang & Woo (2018). Interestingly, the relation by Sun et al. (2017) for R_KDR predicts sizes larger than the R_ENLR at the high-luminosity end. As discussed in Husemann et al. (2016), it is possible that beam smearing leads to a systematic overprediction of KDR sizes. This effect is minimized for 3C 273, due to the high physical resolution and application of the AGN-host deblending.

Interactions of jets with the ambient medium have been extensively studied in a wide variety of jet properties from compact to extended systems (e.g., Tadhunter et al. 1989, 2000; O'Dea et al. 2002; Holt et al. 2008; Morganti et al. 2013; Mahony et al. 2016; Villar-Martín et al. 2017; Tremblay et al. 2018; Jarvis et al. 2019; Kolwa et al. 2019). In the vast majority of these cases, the most extreme gas motions have been associated with the head of the jets as they push through the ambient gas, leaving turbulent gas behind their path. The impact of those outflows have been proposed to either suppress star formation (negative feedback) or even promote star formation (positive feedback). In the first case, the outflow is assumed to disrupt dense cold gas in the host galaxies and drag it outside the host so that less cold gas is able to form stars (e.g., Nesvadba et al. 2006; Schawinski et al. 2009; Fabian 2012). In the second case, a density wave caused by the expanding shock front may allow molecular gas to rapidly form through enhanced cooling in a compressed dense gas phase. The second process has mainly been theoretically predicted (Silk 2005; Gaibler et al. 2012; Ishibashi & Fabian 2012; Zubovas et al. 2013), but some observational evidence for enhancement in star formation has been observed in the filament of Centaurus A (Crockett et al. 2012), in companion galaxies impacted by jets (e.g., Croft et al. 2006; Molnár et al. 2017), and nearby galaxies with strong outflows (e.g., Cresci et al. 2015; Maiolino et al. 2017).

3C 273 hosts a powerful jet and a bright optical nucleus, so that enough energy is released to provide strong feedback on the ambient gas. Given that there is no ionized gas detected on large scales associated with the head of the jet suggests that the jet axis is not going directly through the gas disk in the host galaxy. Still, there is a clear alignment of the kinematically disturbed ionized gas on kiloparsec scales. Here, we refrain from computing gas outflow energetics to distinguish between a mechanically or radiatively driven outflow, because neither the electron densities nor the actual velocities of individual clouds can be robustly measured with our current data. However, the well-constrained jet axis about ∼20° away from our LOS (e.g., Stawarz 2004) helps us draw a potential scenario for the evolution of the outflow. The jet axis provides an independent estimate of the 3D orientation of the primary AGN ionization cone axis, which can otherwise be estimated from the ENLR morphology and kinematics (e.g., Fischer et al. 2013).

Surprisingly, we see preferentially receding gas motion on the approaching jet side and approaching gas motions on the opposite side. This is inconsistent with a purely radiatively driven outflow, which should lead to radial gas acceleration if the ionization cone is roughly aligned with the jet axis. However, the expanding cocoon around the radio jet will lead to shocks propagating in all directions and create turbulent and bulk motion in the surrounding interstellar medium as it expands with time. Due to the inclination of the gas disk, the expanding jet-driven hot gas cocoon will accelerate the gas in transverse or even backflow motion. Indeed, the presence of molecular gas with an arc-like structure just outside the KDR of the outflow may be interpreted as the consequence of a gas density wave generated by the expanding shock perpendicular to the jet. While the molecular gas may also be bound to the remnant of an infalling galaxy, we do not see such a well-defined structure in the ionized gas phase. The jet axis breaks the symmetry and goes exactly through the center of the molecular arc as shown in Figure 1, panel (e), which is unlikely to be a pure coincidence. The motion of the cold gas is not yet affected too much by the outflow and predominantly traces the ordered gravitational motion like the kinematically quiescent ionized gas at the same location. Hence, we can estimate the inclination of the undisturbed rotation-dominated gas disk by comparing the dynamical mass with the expected enclosed mass. The stellar mass of 3C 273 is roughly M_* ∼ 4 × 10¹² M_⊙ according to the M_BH–M_* relation of Häring & Rix (2004). This leads to an intrinsic rotational velocity of v_rot = 400 km s⁻¹ at 8 kpc distance, whereas the observed maximum radial velocity in the ionized and cold gas phase is about v_rad = 150 km s⁻¹. This implies an inclination of ∼20° of the gas disk.

Our proposed scenario is illustrated as a cartoon in Figure 5. As the jet has expanded close to the speed of light for 10⁵ yr in 3C 273 (Stawarz 2004), it has already reached beyond the galaxy and is no longer directly interacting with the gas. However, the continuously expanding cocoon of the jet still affects the gas disk and pushes it in the transverse direction over time. Those transverse and backflow motions from the expanding cocoon have been observed in hydrodynamical simulations of expanding jets in a dense gas medium (Mukherjee et al. 2018). This transverse motion and shock front generates a symmetric density wave in the preshock phase, which allows the formation of a significant reservoir of cold molecular gas through enhanced cooling. While a counter-jet was recently detected in 3C 273 (Punsly & Kharb 2016) and may also explain the disturbed kinematics on the other side of the AGN, it may be intrinsically less powerful and may fail to produce a shock front in which molecular gas can form. Because the ionized emission-line ratios do not show an enhancement of Hβ with respect to [O iii] at the location of the cold gas, this suggests that the cold gas is not actively forming stars. However, the α_CO factor is unknown and diffuse star formation may not be easily recognized, due to the surrounding AGN-ionized region. All of these measurements highlight the complex interplay between outflows and the multiphase medium, which depends on a lot of parameters.

Zoom In Zoom Out Reset image size

While the kinematics may also be interpreted as inflow motion, the velocities of the gas clouds of up to ±800 km s⁻¹ are too high to be driven purely by gravity. Nevertheless, the accelerated gas clouds will contribute to the turbulence in the system and likely enhance the inflow of some material toward the nucleus. Such a chaotic inflow of gas may account for the flicker-noise AGN variability, which is predicted by the CCA model (Gaspari et al. 2017) and was discovered in 3C 273 over a 80 yr light-curve monitoring (Press 1978). A more detailed model of the ionized gas distribution and kinematics to overcome the current limitations of the seeing-limited observations will be performed with future high-resolution observations allowed by the new MUSE narrow-field mode.

It is well known that there is an anticorrelation between the Eddington ratio of an AGN and the strength of the [O iii] line in the NLR (e.g., Boroson & Green 1992; Boroson 2002; Marziani et al. 2003). Due to selection effects, high Eddington ratios dominate the luminous AGN population at high redshift, and it has been proposed that the NLR is disappearing in these AGNs (e.g., Netzer et al. 2004). Indeed, the ENLR on kiloparsec scales remain undetected preferentially around high Eddington-ratio AGNs in a sample of luminous AGNs at z < 0.3 (Husemann et al. 2008). 3C 273 belongs to the high Eddington-ratio regime with L_bol/L_Edd ∼ 0.6 (Husemann et al. 2013) and shows a weak [O iii] line in the AGN spectrum. The ENLR was indeed reported to be undetected in initial IFU observations of this source (Husemann et al. 2013). We conclude from this work that the high contrast between the luminous AGN and the ENLR has apparently caused previous nondetections. It is simply more difficult to detect the ENLR underneath a beam-smeared point-like AGN signal and requires much higher S/N observations to detect the [O iii] line for a given surface brightness limit. Several AGN surveys at high redshift, such as KASHz (Harrison et al. 2016), WISSH (Bischetti et al. 2017), or SUPER (Circosta et al. 2018), have attempted to detect ionized gas around AGNs. Our MUSE data provide an ideal data set to empirically constrain the distribution, kinematics, and surface brightness of the emission lines around a hyperluminous quasar. It is therefore instructive to simulate observations of 3C 273 at z = 2.3 where [O iii] is shifted into the H band.

The angular scale at z = 2.3 is 8.2 kpc/arcsec, which is nearly a factor of 3 larger compared to the physical resolution at the redshift of 3C 273. In order to observe such a source at the same spatial resolution as provided by the MUSE data for 3C 273, the angular resolution would need to be three times better than our seeing (i.e., 023). This matches the typical angular resolutions achieved with adaptive optics observations in the H band. Furthermore, the spectral resolution of MUSE is almost the same as the one of the SINFONI H-band grating or the high-resolution gratings of NIRspec aboard the James Webb Space Telescope. Hence, we have to scale the integrated flux of the MUSE cube after redshifting the spectra by

$$\begin{eqnarray}&&f=\displaystyle \frac{{D}_{L}{\left(z=0.158\right)}^{2}}{{D}_{L}{\left(z=2.3\right)}^{2}}\displaystyle \frac{(0.158+1)}{(2.3+1)}=5.9\times {10}^{-4}\end{eqnarray} \tag{ 1 }$$

and reassigning the pixel scale from 02 to 0066 or 005 for practical reasons to avoid subpixel resampling.

We then simulate realistic SINFONI observations at z = 2.3 by binning 2 × 2 pixels in the MUSE cube, scaling the fluxes according to Equation (1) and redshifting the spectra. Empirical noise (without sky lines) is added for SINFONI for 1 hr and 9 hr of on-source time as directly inferred from real observations with a 01 plate scale and a 3'' × 3'' FoV. We also simulate data for the NIRspec IFU with the same on-source exposure times adopting the noise model of the official ETC (Pontoppidan et al. 2016). In Figure 6, we show the narrowband images after subtracting the QSO contribution with qdeblend^3D in the same way as for the original MUSE data. Ground-based NIR spectrographs such as SINFONI are hardly able to detect the very brightest inner region of the ENLR even after 9 hr of integration time, which may explain the low detection rate of kiloparsec-scale [O iii] outflows around optically bright AGNs (Vietri et al. 2018). NIRSpec will be able to recover the brightest part of the ENLR, but likely miss the full ENLR extent. High-redshift QSO studies therefore need to be aware that a much higher S/N is needed to detect the diffuse extended outflows underneath highly luminous AGNs when designing observations.

Zoom In Zoom Out Reset image size