JWST discovers an AGN ionization cone but only weak radiative-driven feedback in a powerful z ≈ 3.5 radio-loud AGN

We present the first results from a JWST program studying the role played by powerful radio jets in the evolution of the most massive galaxies at the onset of Cosmic Noon. Using NIRSpec integral field spectroscopy, we detect 24 rest-frame optical emission lines from the z = 3 . 5892 radio galaxy 4C + 19.71. 4C + 19.71 contains one of the most energetic radio jets known, making it perfect for testing radio-mode feedback on the interstellar medium (ISM) of a M ⋆ ∼ 10 11 M ⊙ galaxy. The rich spectrum enables line ratio diagnostics showing that the radiation from the active galactic nucleus (AGN) dominates the ionization of the entire ISM out to at least 25 kpc, the edge of the detection. Sub-kpc resolution reveals filamentary structures and emission blobs in the warm ionized ISM distributed on scales of ∼ 5 to ∼ 20kpc. A large fraction of the extended gaseous nebula is located near the systemic velocity. This nebula may thus be the patchy ISM which is illuminated by the AGN after the passage of the jet. A radiatively-driven outflow is observed within ∼ 5kpc from the nucleus. The ine ffi cient coupling ( ≲ 10 − 4 ) between this outflow and the quasar and the lack of extreme gas motions on galactic scales are inconsistent with other high-z powerful quasars. Combining our data with ground-based studies, we conclude that only a minor fraction of the feedback processes is happening on < 25kpc scales.


Introduction
Active galactic nucleus (AGN) feedback can be categorized into either quasar mode, where the output from the supermassive black hole (SMBH) is coupled radiatively to the gas, or radio mode, where the feedback is due to the kinetic energy of the jet (e.g., Fabian 2012;Harrison 2017).The quasar mode is known to be prevalent (e.g., Harrison et al. 2014).The radio mode also plays an important role in AGN-gas interactions, leading to quenching.This is especially true in massive galaxies in the early Universe, as shown by both numerical simulations and observations (e.g., Heckman & Best 2014;Mukherjee et al. 2018;Kondapally et al. 2023).Despite the relatively short lifetime of a powerful jet (∼10 7 yr), observations do indeed indicate that it can still significantly impact the host galaxy (e.g., total energy output by the jet 10 60 erg, Nesvadba et al. 2006;Miley & De Breuck 2008).Relativistic jets trigger shocks and drive outflows through an expanding overpressured bubble (e.g., Begelman & Cioffi 1989).While hydrodynamic simulations demonstrated that jets have the ability to interact with the interstellar medium (ISM) out to several kiloparsecs (e.g., creating turbulence, driving outflows, and compressing the gas, Dugan et al. 2017;Mukherjee et al. 2021), observations paint a more complicated picture.Specially, jet-ISM interactions are found to be different between AGN with intermediate radio power (L 1.4 GHz = 10 23−25 W Hz −1 , e.g., Mullaney et al. 2013) and those with a higher power (e.g., Mukherjee et al. 2016).Moreover, it is difficult to study the feedback mechanisms of the most powerful jetted AGN due to observational limitations and the simultaneous presence of energetic, quasar mode feedback.
High-redshift radio galaxies (HzRGs) are the best targets to test AGN feedback in massive host galaxies (∼10 11 M , De Breuck et al. 2010) because the gas, dust, and stellar populations in their host galaxies can be observed without contamination from the bright, point-like quasar light.In contrast to low-z, where radio mode AGN are believed to have low black hole accretion rates with inefficient radiative luminosity, HzRGs are observed to have both a vigorous radio mode and quasar mode energy output (e.g., Vernet et al. 2001;Nesvadba et al. 2007Nesvadba et al. , 2017a,b),b).There are numerous studies focused on the outflows on scales from several to tens of kiloparsecs for both low-and high-redshift radio galaxies, which find evidence of them being jet driven (Tadhunter et al. 2001(Tadhunter et al. , 2007;;Nesvadba et al. 2006Nesvadba et al. , 2007)).A closer look (scales of tens of parsecs) at the feedback in the vicinity of radio AGN is possible at low-z (e.g., Tadhunter et al. 2003).However, spotting the warm, ionized quasar outflow around cosmic noon (Madau & Dickinson 2014) down to subkiloparsec scales is challenging.Therefore, it is still unclear how energetic jet mode and quasar mode feedback couple with the ISM, how efficient these mechanisms are in driving outflows, and on which scales the different feedback mechanisms dominate.All of these are critical issues to be addressed before we can achieve a deeper understanding of the quenching of star formation in these massive host galaxies (Falkendal et al. 2019).
We can now finally take a step forward thanks to the Near Infrared Spectrograph (NIRSpec; Jakobsen et al. 2022) on board JWST, which provides an order of magnitude improvement in sensitivity.At least for HzRGs, there is evidence that the morphology and/or kinematics of the ionized ISM and even circumgalactic medium (CGM) gas are impacted by jets (Nesvadba et al. 2008(Nesvadba et al. , 2017a;;Falkendal et al. 2021;Wang et al. 2023), which may indicate their role in feedback.However, the situation may be different when using deeper observations (e.g., Wylezalek et al. 2017).For example, how does the ISM of HzRGs look in the vicinity of and inside the host galaxies (e.g., ∼10 to 20 kpc).NIRSpec integral field unit (IFU) observations of powerful quasars at cosmic noon have already revealed the detailed ionization mechanisms and quantified the outflow rates on galactic scales (e.g., Wylezalek et al. 2022;Vayner et al. 2023Vayner et al. , 2024)).NIRSpec will revolutionize our view of the inner ISM of HzRGs where the detailed physics are still unclear.A resolution-matched comparison can finally be done for high-and low-z radio AGN.NIRSpec IFU will also enable observational tests of simulated jet feedback (e.g., Mukherjee et al. 2016) in the early Universe.
In this work, we present the first JWST view of the ionized ISM of a z ∼ 3.5 radio-loud AGN.In Sect.2, we describe the observations and data processing.We then show in Sect. 3 the morphology of the warm ionized gas and identification of the observed optical emission lines from the extracted spectra.Finally, we discuss the impact of the jet on the ISM using line-ratio diagnostics and summarize in Sect. 4. Throughout this paper, we assume a flat Λ cosmology with H 0 = 70 km s −1 Mpc −1 and Ω m = 0.3.Following this cosmology, 1 arcsec = 7.25 kpc at the redshift of 4C+19.71.Throughout this work, we use z = 3.5892 ± 0.0004 as determined from atomic carbon, [C i](1−0) (ν rest = 492.16GHz), as the systemic redshift (Falkendal et al. 2021;Kolwa et al. 2023).Due to the relatively large synthesized beam (1.8 ×1.9 , Falkendal et al. 2021) of the Atacama Large Millimeter/submillimeter Array (ALMA) [C i](1−0) observation, it is uncertain whether the molecular gas is located in the host galaxy (Appendix D).Though JWST observes in a vacuum, we label the emission lines with their air wavelengths in Å following the convention, for example [O iii]5007.

4C+19.71, observation, and data processing
4C+19.71 (MG 2144+1928) is the first target observed as part of the JWST Cycle 1 General Observers program (ID 1970, PI: W. Wang) "Zooming into the monster's mouth: tracing feedback from their hosts to circumgalactic medium in z = 3.5 radio-loud AGN".The program includes four HzRGs with a diversity of star formation rates and radio jet morphologies (van Ojik et al. 1996;Carilli et al. 1997;Nesvadba et al. 2007;Falkendal et al. 2019).We selected them to have similar redshifts (z 3.5) to maximize the number of emission lines in a single observational setting.
4C+19.71 has a FR II-type radio jet with a double-sided lobe extending over 60 kpc (Fanaroff & Riley 1974;Carilli et al. 1997).Besides the radio data, there are multiwavelength studies available for 4C+19.71(e.g., Pentericci et al. 1999;Maxfield et al. 2002;Seymour et al. 2007;De Breuck et al. 2010).The HzRG is surrounded by an X-ray halo and [O iii]5007 nebula with similar extension.Both are elongated in the same direction as the jet (∼60 kpc; Armus et al. 1998;Smail et al. 2012;Nesvadba et al. 2017b).A more extended, ∼100 kpc Lyα nebula is also detected around the radio AGN.It also shows an elongated morphology along the jet axis, even beyond the radio lobes (Falkendal et al. 2021;Wang et al. 2023).
The stellar mass of the host galaxy of 4C+19.71 is constrained to be 10 11.1 M .It forms stars at the rate of ∼84 M yr −1 (Falkendal et al. 2019).Finally, there is significant molecular gas in the host (2.5 × 10 10 M , Kolwa et al. 2023).
The observation of 4C+19.71 was conducted on UT 2022 October 30 with NIRSpec in IFU mode (Böker et al. 2022;Jakobsen et al. 2022).We selected G235H/F170LP as the disperser and filter setup, which covers 1.70−3.15µm in the observed frame or ∼3700−7620 Å in the rest frame.The combination of filter and grating offers the spectral resolution of 85−150 km s −1 .This ensures major optical emission lines seen in type-2 AGN ([O ii]3726, 3729 to Hα, including [S ii]6716, 6731) can be observed at once (e.g., Zakamska et al. 2003).We adopted a 9-point dither pattern and the improved reference sampling and subtraction (IRS2 ) read-out pattern.The total onsource exposure time was 3.7 h.We designed the leakage exposure, 0.4 h, to identify light leaking through the failed open shutters on the microshutter array (MSA) at the first dither position.
We downloaded the raw data from the Mikulski Archive for Space Telescopes (MAST).The data were processed with the JWST Science Calibration pipeline 1 (v1.9.4) with the Calibration Reference Data System (CRDS) context file jwst_1041.pmap.Our procedure is similar to Vayner et al. (2023), executing the standard first and second stages of the pipeline with one modification.Specifically, we used the "emsm" method instead of "drizzle" during the data cube construction step to reduce the low-frequency ripples due to undersampling.The third and final stage of the processing is to combine cubes from different exposures into the final cube.We instead used the script from Vayner et al. (2023), who applied the Python package reproject because the third stage of the pipeline falsely rejects the bright emission peak while keeping some noise spikes.During this step, sigma clipping (2σ) was also included to reject outliers.The data cube was resampled to have a pixel scale of 0.05 .We processed the raw NIRSpec IFU data of the standard star TYC 4433-1800-1 (PID 1128) for flux calibration.Finally, we performed an additional background subtraction to the flux-calibrated science data cube to alleviate negative background continuum seen in source-free regions.The size of the point spread function (PSF) of the NIRSpec IFU depends on wavelength and spatial position.We mark the full width at half maximum (FWHM) of the PSF constructed in Vayner et al. (2024) in Fig. 1a.
Several works based on NIRSpec IFU observations report a World Coordinate System (WCS) offset in the final data cube (e.g., Wylezalek et al. 2022;Perna et al. 2023).We aligned the continuum emission position of the west foreground galaxy in our NIRSpec cube to the Hubble Space Telescope (HST) WFPC2 image and used this new WCS in the following analysis (see Fig. 1a and Pentericci et al. 1999).The astrometry of the HST image was corrected by cross-matching with Gaia DR3 targets (Gaia Collaboration 2023).Our final WCS solution resulted in a shift of ∆RA = 0.43 and ∆Dec = −0.22 compared to the pipeline output WCS.We discuss the astrometry correction in Appendix A.
We also present the ALMA Band 8 observation of 4C+19.71.This observation was carried out under program ID 2021.1.00576.S (PI: W. Wang) on UT 2022 June 17.In this work, we use only the archived image at ν obs = 400.3GHz to show the position of the cold dust emission.A detailed analysis of the ALMA observations is beyond the scope of this paper and will be presented in a forthcoming paper.

Results
In this paper, we inspect both the 2D morphology and 1D spectrum of the warm ionized gas.Given the signal-to-noise ratio (S/N) of the optical emission lines covered in our setting, we choose the [O iii]5007 ([O iii] hereafter if not specified) as the proxy for the morphology.It is extensively used in the studies of AGN narrow line regions (e.g., Husemann et al. 2013;Harrison et al. 2014;Wylezalek et al. 2016).This is the first time that the distribution of ionized gas around a z > 3.5 radio loud AGN is seen with 1 kpc resolution.100, 100] km s −1 ), and red wing ([500, 700] km s −1 ), respectively.On each of the images, there are two distinct peaks located south and north which roughly align with the jet axis (dark purple contours in Fig. 1a, Carilli et al. 1997).Based on [O iii] gas kinematics probed by VLT/SINFONI near the jet hot spot positions (Fig. A.7 of Nesvadba et al. 2017a), we assume the southern jet is the approaching one.There is no doubt that 4C+19.71has a hidden quasar (e.g., Seymour et al. 2007;De Breuck et al. 2010;Falkendal et al. 2019, and Sect. 4.1).Assuming the ionization cone of the quasar is aligned with the jet axis (Drouart et al. 2012), it is likely that we are observing, in both the south and north regions, the foreground and background parts of the cone.

Morphology of the extended ionized nebula
To assist in studying the morphology, we show example [O iii]4959, 5007 spectra extracted at four different locations (Figs.1c1-c4).We note that the analysis of the extended emission-line region kinematics is not the focus of this paper.
The per-spatial-pixel (spaxel) spectra are fit using q3dfit (Rupke et al. 2023), which is a Python tool for analyzing JWST IFU data based on the software IFSFIT (Rupke 2014; Rupke et al. 2021) 2 .It is clear that at least three different kinematic components are present on galactic scales ( 10 kpc).This is especially obvious for the gas in the south, with three distinct peaks present (e.g., Fig. 1c2).We stress that we refer to the component closest to 0 km s −1 as "systemic" and the others as "blue" or "red" depending on their relative velocity shift with respect to the systemic component.At this point, it remains unclear whether they are the blueshifted or redshifted outflowing gas clouds.
The most striking morphological features are the spatially extended emission blobs and filamentary structures around the systemic velocity (Fig. 1b1).For example, both components detected in the east filament (Fig. 1c1) and the gas of the northern spot (Fig. 1c3) all have |∆v| < 300 km s −1 .Neither the previous narrowband imaging nor ground-based IFU could resolve these components spatially and spectrally (Armus et al. 1998;Nesvadba et al. 2017a).With the relatively small field of view (FoV) of the NIRSpec IFU, the jet hot spots are not captured which are presumed to be the regions with the most energetic jetgas interactions.Given the fact that these extended filamentary features sit at the systemic redshift, we arrive at the conclusion that they are not associated with the outflow but may be high density gas clumps.Their morphologies may have been disrupted by the jet-induced bubble (Mukherjee et al. 2016;Dugan et al. 2017).Indeed, at least in the example spectra extracted at the east filament, we observe a broad component indicative of disruption by the jet (FWHM = 1052 km s −1 ).

Continuum emission
We overlay in Fig. 1b the continuum emission detected in our ALMA Band 8 observation, which likely indicates the location of the cold dust emission heated by the newly formed stellar population (e.g., Herrera-Camus et al. 2021).In Fig. 1d we show that the AGN position (NIRSpec continuum, Sect.4.1) agrees with the ALMA dust location within their uncertainties.The 1σ uncertainty, ∼0.02 , of the ALMA dust position is estimated 2 https://q3dfit.readthedocs.io/en/latest/using pos acc = beam FWHP /(S /N)/0.9 with beam FWHP 0.2 and S /N 10 (see Appendix A for the 1σ uncertainty, ∼0.06 , of the NIRSpec continuum) 3 .
A previous HST image (WFPC2/F702W) indicates other two objects in the NIRSpec IFU FoV (Pentericci et al. 1999).They are also seen in our NIRSpec observations (Fig. 1a).We label them G1 and G2 for the galaxy in the west and east, respectively.G1 is used as the WCS alignment reference in Sect. 2. They are identified as foreground galaxies with z 1.786 for G1 and z 1.643 for G2.Their spectra and emission line identifications are shown in Appendix B.

Emission line identification
We present the identification of the emission lines in Fig. 2.This is the first time that optical nebular emission lines from [O ii]3726, 3729 to [S ii]6716, 6731 are observed for z 3.5 HzRGs with the spectral resolution of 85 to 150 km s −1 .This makes it possible to investigate the ionization mechanisms and metallicities (e.g., line ratio diagnostic diagrams, Baldwin et al. 1981;Veilleux & Osterbrock 1987;Kauffmann et al. 2003;Groves et al. 2004bGroves et al. , 2006;;Kewley et al. 2001Kewley et al. , 2006Kewley et al. , 2013;;Wylezalek et al. 2017).
We extracted a spectrum from a 0.2 × 0.2 or 1.4 × 1.4 kpc 2 (4 × 4 pix 2 ) square aperture near the position of the continuum emission peak.The aperture and size were chosen to also include the [O iii] emission peak which maximizes the S/N for emission line identification.The aperture was aligned with the extension of the brightest [O iii] emission blob.We show this aperture in the black box in Fig. 1d.The spectrum is presented in Fig. 2. We summed over five wavelength pixels around the systemic zero (velocity range of ∼270 to ∼150 km s −1 ) to calculate the line S/N for identification.We identify ∼24 emission lines (doublets) with S /N > 3. The ones with S /N 10 are marked by red lines.We implemented a visual check by examining the spatial distribution of lines with 3 < S /N < 10.The spatially resolved ones are marked by black dashed lines.
We fit the spectrum with q3dfit.We used three kinematic components.The line center and width of the same component are connected for different lines.We present these kinematic results in Table 1.We refer to them as central, red, and blue components based on their velocity shifts.In Appendix C, we report fit line fluxes from the 1D spectrum.

Comparison to Cygnus A
Cygnus A (3C405) is the most powerful radio AGN in the local Universe with comparably high radio power as 4C+19.71(z = 0.0562, P 178 MHz ∼ 5.5 × 10 28 W Hz −1 , Carilli & Barthel 1996).This makes it the perfect target to compare with the properties of the warm ionized ISM gas using the common optical tracers (Fosbury et al. 1999;Vernet 2001).Using the spectropolarimeter on the Keck telescope, Ogle et al. (1997) observed the nebular emission lines of Cygnus A. We overlay our 1D spectrum (Fig. 2) with the spectrum of Cygnus A extracted at its nucleus in Fig. 3.The extraction aperture of Cygnus A is ∼1 × 1.1 which corresponds to 1.1 × 1.2 kpc 2 and is similar to the aperture used for 4C+19.71(1.4 × 1.4 kpc 2 ).This is the first time that a rest-frame optical spectrum of a z 1 radio AGN can be studied in comparable detail as a local example.1d).We mark the emission lines detected with S /N 10 in red lines.Lines detected with 3 < S /N < 10 and visually checked to be extended are marked with gray lines.We note that the feature at ∼2.945 µm is likely instrumental due to the leakage flux from the MSA.The central gray shaded region is the detector gap.The top panel is the full view of the spectrum while the zoom-in view of the faint lines are shown in the bottom panels.The magenta curves show the best fit of the spectrum (Appendix C).
To first order, the continuum shape of 4C+19.71 is consistent with Cygnus A. We note that the Cygnus A spectrum is corrected for Galactic reddening and had its host galaxy continuum subtracted (Ogle et al. 1997).Though there might be contamination of continuum (i.e., due to the leak of from the MSA) for our NIRSpec data, we conclude that the emission from evolved stellar population in the host of 4C+19.71 is faint, that is, the We conclude that, to first order, AGN photoionization mechanisms dominate the inner parts of these two galaxies.The A169, page 5 of 15 similarity of the Balmer lines, Hβ, Hγ, and Hδ, suggests the ISM gas conditions are similar (Osterbrock & Ferland 2006).
The dissimilarity of the flux at redder wavelength (e.g., [O i], [N ii], and [S ii]) may suggest differences in gas enrichment.
van Bemmel et al. (2003) proposed that the quasar wind-driven outflowing dust clouds scatter the optical line emitting photons in Cygnus A. We know there are at least three kinematic components present in the nucleus of 4C+19.71(Sect.3.1).We defer to a future publication the detailed study of the gas kinematics to unveil the spatially resolved gas motions and constrain the ISM distribution.

Line ratio diagnostics and ionization mechanisms
To study the extended irregular ISM structures around the hidden powerful radio AGN (Sect.3.1), we construct line ratio maps and compare them to emission line diagnostics (e.g., Groves et al. 2004bGroves et al. , 2006;;Kewley et al. 2006).We use q3dfit to fit the cube to systematically separate the blended emission lines and account for doublets.The fit region was selected based on the [O iii] narrowband image with a S/N cut of >5 (white contour in Fig. 4).We limit the spatial study in this paper to the six The Balmer line ratios are sensitive to dust attenuation (e.g., Veilleux & Osterbrock 1987;Osterbrock & Ferland 2006).We present the color excess, E(B − V), maps based on the fit Hα/Hβ flux ratios in Fig. 4a.The map is constructed using an intrinsic flux ratio of 3.1 assuming Case B recombination in AGN ionization and a Calzetti et al. (2000) extinction law following Vayner et al. (2023).We find that there is a dusty region with high E(B − V) ∼ 1 around the AGN position, coincident with our ALMA dust continuum (pink contours in Fig. 4a), confirming the importance of dust causing extinction near the AGN.This spatial coincidence also validates our manual astrometry correction based on foreground galaxy (Appendix A).Fosbury et al. (1999), Vernet (2001) proposed that a dust ring is present around the nucleus of Cygnus A. Given the spectral similarity of 4C+19.71 to Cygnus A (Fig. 3), it is possible that we also detect a dust lane around the center of 4C+19.71.The line ratios shown in Figs.4b and c  they fall in the line ratio space spanned by the northern regions (green squares) in all diagrams.Despite having a relatively large dispersion, the median line ratios of the northern region is consistent with the spaxels from the nucleus region (black triangles).We report the median line ratios of each region in Table 2.We overlay the Groves et al. (2004a) dusty radiationpressure dominated model grids on the diagnostic diagrams.
We chose the 2 Z models (black) as indicated by the [N ii]/Hα ratios (e.g., metallicity from Groves et al. 2006;Nesvadba et al. 2017a).The models are constructed with various combinations of the power-law index, α, of the ionizing source (i.e., F ν ∝ ν α ) and ionization parameter, U = S /(n H c), where c is the speed of light, n H is the hydrogen number density and S is the flux of ionizing photons entering the cloud.The goal is not to perform a quantitative comparison, but rather to test which parameter(s) could be responsible for the line ratio differences.Although other interpretations are possible, a higher U parameter could explain the lower [O i]/Hα and a higher [O iii]/Hβ (e.g., Villar-Martín et al. 1999).Our comparison indicates a scenario where the north, south, and nucleus regions have higher U parameters than the EW region.This is expected given that the jet axis is along the north-south direction which is also the direction of the AGN ionization cone (Drouart et al. 2012).Hence, the ionizing photons of the AGN can reach the ISM with less extinction along the north-south direction.In this configuration, the EW region is outside of the ionization cone where the photons will encounter more obstacles before ionizing the gas (e.g., Fosbury et al. 1999;Vernet et al. 2001).This indicates larger particle density n H and less ionizing photon flux (smaller S ) in the EW which leads to lower observed U parameters.
The south region has lower [N ii]/Hα, [O i]/Hα, [S ii]/Hα, and [O ii]/[O iii].The U parameter appears at similar levels in the south as in the nucleus.Comparing with the model grids, this could indicate that the ratio difference is due to lower α (steeper power-law ionizing spectrum).If this is the case, then there may be additional extinction between the ionizing photons and the southern clouds.The attenuation map (Fig. 4a) indeed A169, page 6 of 15   1d).The dashed white lines and circle in panel b and c divide the FoV into five regions for line ratio diagnostic analysis (Fig. 5).FWHM of the JWST PSF is marked at the left corner.
shows a higher value toward the south.If this is due to the geometry of the dust torus, it may be inconsistent with the indicated orientation from Nesvadba et al. (2017a).We also overlay the similar radiation model with 0.5 Z in cyan in Fig. 5c and find that there is an alternative possibility that the southern region is more metal poor than the nuclear region.Additionally, the complex kinematics in the southern region (Fig. 1c3), suggesting a more complex scenario in the south (e.g., a late-stage merger, see Sect.4.5).The detailed analysis of this southern region requires more lines which is beyond the scope of this paper (e.g., He lines, Groves et al. 2004a;Holden et al. 2023).
To further test the ionization mechanisms, we also overlay the empirical classification from Kewley et al. (2006) in Figs.5a-c.Although these are calibrated using low-z data (i.e., may not be relaible for high-z, see Kewley et al. 2019), the results provide evidence of the AGN radiation dominating ionization of the ISM within 30 kpc of the hidden quasar.Some spaxels in the north and east regions are located in the low A169, page 7 of 15 ionization or H ii region.This suggests a different ionization mechanism or mix of different mechanisms at these spatial locations.Indeed, we find a part of them are roughly consistent with the "shock+precursor" models (light purple grids Fig. 5b, Allen et al. 2008).The shock models can cover a large fraction of the parameter space of the line ratio diagnostics.The identification of shock ionization requires combination of flux ratios with kinematics, electron temperature and Hα fluxes (e.g., Alatalo et al. 2016;Baron et al. 2017;Alarie & Morisset 2019).Our discussion is not meant to be decisive but to be inclusive with the possibility of shock.We choose the models with pre-shock number density, n = 100 cm −3 (e.g., De Breuck et al. 2000).For the shock velocities, we use higher values (v shock ≥ 200 km s −1 ) following De Breuck et al. (2000) studying the emission lines of HzRGs.We note that shock+precursor model with pre-shock density n = 1000 cm −3 are also physical for AGN (Liu et al. 2020) which occupies a subset of the line ratio space spanned by n = 100 cm −3 model using the same set of v shock and magnetic parameters.The driver of the shock is unclear.The EW region is roughly in the direction perpendicular to the ionization cone (jet axis) which suggests a less collimated shock driver or the ionization cone is not perpendicular to the dusty obscuring disk (e.g., Dugan et al. 2017;Tanner & Weaver 2022;Vayner et al. 2023).As for the more extended parts ( 10 kpc) in the north, the passage of the powerful radio jet could also have left remaining shock effects.Our results show that the inner part of the ionized ISM in the z ∼ 3.5 radio-loud AGN is dominated by photons from the hidden AGN with potential evidence for shock ionization.Its structure agrees with the classical ionization cone.Although the analysis is based on knowledge gained at low-z and there are more unsolved complexities, JWST will be the key to address this (e.g., Nakajima & Maiolino 2022;Sanders et al. 2024).
A169, page 8 of 15 In Sect.3.1, we show that the "high-velocity" clouds are only detected within the ∼5 kpc from the AGN and their morphology is aligned with the jet axis (Fig. 1).This is inconsistent with the results from observations of quasars at cosmic noon where fast (|∆v| ∼ 500 km s −1 ) outflows are seen at 10 kpc and are less directionally confined (e.g., Vayner et al. 2024).However, Cresci et al. (2023), Veilleux et al. (2023) also showed one contradicting case.Hence, this discussion is yet settled and required more observations.Through our discussions in Sects.4.1 and 4.2, these may be the outflows inside the ionization cone.The double-sided radio lobes have a projected extent of ∼60 kpc, which is beyond the JWST FoV.Hence, this potential outflow maybe radiatively driven.
The systematic per-spaxel kinematic study is beyond the scope of this work.Additionally, as shown in Sect.4.2 (also in Sect.4.5), there may be more complex scenarios happening in the southern cone.Hence, we use the results from the integrated 1D spectrum (Figs.1d and 2, Tables 1 and C.1) at the nucleus region as a proxy to estimate the kinetic energy of the potential outflow near the AGN.For this calculation, we focus on the red and blue components (Tables 1 and C.1) and assume them to be the redshifted and blueshifted outflow clouds seen on the back and front side of the northern part of the ionization cone, respectively.Then, we assume their relatively velocity shifts with respect to the central component (Appendix D) as their outflow velocities (∆v red = 334 km s −1 and ∆v blue = −386 km s −1 ).We follow Vayner et al. (2024) and use the following equation to calculate the ionized outflow gas mass (Osterbrock & Ferland 2006): where L Hα and n e are the Hα luminosity and electron density, respectively.We assume an electron temperature T e 15 000 K (Nesvadba et al. 2017a;Vayner et al. 2024).We use the fiducial value n e = 500 cm −3 from Nesvadba et al. (2017a).Given that the emissivity, j H α , is not sensitive to n e , we use the value of 2 × 10 −25 erg cm 3 s −1 for T e 15 000 K (Vayner et al. 2024).Based on Eq. (1), we find M ion = 10 5 and 10 6 M for the red and blue components, respectively.We then estimate the outflow rate using the following equation (Vayner et al. 2024): For the distance ∆R, we use half of the physical size of the extraction aperture, 0.7 kpc.For the velocity, we use the v ion = |∆v|+σ v (Rupke & Veilleux 2013;Vayner et al. 2024).We report the summation of the red and blue components as the estimated total outflow rate.This gives a total Ṁion of 2 M yr −1 .The outflow momentum flux Ṗion and outflow kinetic luminosity L kinetic ion are calculated from We calculate values for Ṗoutflow ion and L kinetic ion of 1.5 × 10 34 dyne and 8.0 × 10 41 erg s −1 , respectively.These are much smaller (∼2 dex) than the values derived on larger scales by SINFONI (Nesvadba et al. 2017a,b).If we use the infrared AGN luminosity, L IR AGN = 10 10.91 L , and assume a conversion factor, κ bol AGN = 6, the bolometric luminosity of 4C+19.71 is estimated to be L bol ∼ 2.5 × 10 47 erg s −1 (Drouart et al. 2014;Falkendal et al. 2019).Then the coupling efficiency of this potential radiatively driven outflow, L kinetic ion /L bol , is ∼10 −5 .Using NIRSpec IFU, Perna et al. (2023) and Vayner et al. (2024) reported higher coupling efficiency for z ∼ 3 quasars (0.02 and 10 −3 , respectively).Veilleux et al. (2023) found a relatively low coupling efficiency using the same instrument, 1.8×10 −4 , for a z ∼ 1.6 type-1 quasar which is still 1 dex higher than 4C+19.71.Harrison et al. (2018) summarized the coupling efficiency for a large number of AGN (e.g., their Fig. 2).Compared to their statistics, the value from this work is located at the lower end.Our estimation based on the JWST observation is a ∼2 to 3 dex lower than the coupling efficiency of the jet kinetic energy probed on larger scales for a sample of ∼50 HzRGs (10 −3 to 10 −2 , Nesvadba et al. 2008Nesvadba et al. , 2007Nesvadba et al. , 2017a,b),b).
Our NIRSpec IFU data is only sensitive to the central clumps, 20 kpc, but not to extended diffuse emission.If we assume there are fast outflows (∼500 km s −1 ) and use the distance from the center quasar to the edge of our FoV (∼20 kpc), the crossing time of this outflow would be ∼40 Myr.This is shorter than the jet age, ∼60 Myr, of 4C+19.71(Nesvadba et al. 2017a).Hence, it is plausible that we do not capture the fast and (presumably) energetic outflowing gas on large scales (e.g., size of the radio lobes).Nevertheless, if the radiative coupling efficiency is 10 −5 at the nucleus where the radiation is the strongest, it will be less efficient farther outside.Simulations (e.g., Costa et al. 2018) find that the efficiency for the radiation-driven outflows can be ∼0.1% in luminous quasars (L bol 10 47 erg s −1 ) within galactic scales ( 10 kpc).The efficiency drops after ∼10 Myr which may resemble the case of our observation.Though the ISM and CGM gas of 4C+19.71seem more quiescent than in other HzRGs (Nesvadba et al. 2017a;Wang et al. 2023), a gas kinetic energy luminosity of 10 45 erg s −1 was found for this source by including the ISM at the jet lobes.Comparing to the results from Nesvadba et al. (2017a,b), we conclude that, at least for this HzRG, the mechanical feedback from the radio jet takes the leading role for driving the outflow, and it is happening in the radio lobes beyond the galactic scale (e.g., ∼30 kpc from the AGN).
One of the main uncertainties is the estimation of M ion , which is inversely proportional to n e .We find that M ion (and thus L kinetic ion ) could be ∼2 dex lower if we use the n e estimated based on our observation (Appendix C).Therefore, the radiative coupling could be even more inefficient.Another uncertainty is the outflow speed.We take the unknown orientation of the ionization cone into account by using v ion = |∆v| + σ v .Even if the estimated speed is one of the three components in 3D space (i.e., the intrinsic v ion is higher by another factor of three), the radiative coupling efficiency will increase by a factor of nine, ∼10 −4 , which is still ∼1 dex lower than the kinetic efficiency from Nesvadba et al. (2017a).

AGN illuminates the filamentary ISM after jet passage
In Sect.3.1, we found that ionized ISM emission with 30 kpc of the AGN is dominated by gas close to the systemic redshift (Fig. 1c).Due to the limitation of the FoV, we only captured the northern extended region.In general, the gas morphology is elongated along the north-south direction, aligned with the jet axis.When we focus on the smaller scale structures (∼5 kpc), these gaseous clouds are filamentary and patchy and do not follow the collimated jet but have larger angle separations (∼30 deg).As discussed in Sect.4.2, we are observing the ionization cone-like structures along the jet.This provides enough ionizing photons to illuminate the extended gaseous nebula beyond A169, page 9 of 15 the galactic scale.The fact that we do not observe the extended ISM in directions outside of the ionization cone, for example in the west where we have the spatial coverage (Fig. 1), may be due to the shortage of ionizing photons.Hence, it may be that the patchy filaments in the north are the "intrinsic" structures, which are not seen in other directions as they are not emitting or simply outside the NIRSpec FoV.
The HzRGs we observed may have entered their fluxlimited parent radio surveys because their radio lobes are brightened when they encounter a higher density environment which provides a larger "working surface" for the jet (Eales 1992;Wang et al. 2023).Combining the discussion in Sect. 4.3 and Nesvadba et al. (2017a), the situation of 4C+19.71 may be that the jet was launched ∼60 Myr ago, and in its younger years could have pushed some gas outside of the host galaxy.Once the radio lobes have escaped the host galaxy, the radio jet is very collimated and is no longer interacting with the ISM; the observed systemic gas in Fig. 1 is thus no longer kinematically affected by the radio jets, but the gas is still photoionized by the AGN photons within the ionization cone along the jet axis.Fu & Stockton (2009) found similar quiescent nebulae whose morphology do not fully agree with their radio maps in low-z quasars on comparable physical scales.They proposed that the expanding bubbles following the jet disturb the ISM.4C+19.71has one of the most powerful radio jets, P jet ∼ 10 47 erg s −1 (using log(P 1.4 GHz /W Hz −1 ) = 28.6 ± 0.1, Nesvadba et al. 2017a;Cavagnolo et al. 2010).The situation may be different than the low-z low radio power quasars.Nevertheless, we find evidence of shock ionization on the extended northern nebula of which the jet could be the driver (Sect.4.2).Further analysis (e.g., T e mapping) will be helpful to determine the mechanisms.

Complex southern ISM
In Sect.4.3, our discussion of the outflow is based on the northern region near the quasar as the nature of the southern "highvelocity" region is more complicated.In Fig. 1c2, we clearly find that there are three velocity components with the "highvelocity" ones (blue and red) dominating the flux.The morphology of the blue and red parts are not spatially overlapping (Figs. 1a,b1,and b3).Such kinematics resemble that of a rotating disk.If we simply take the velocity shifts of the blue and red components from the one spaxel spectrum (Fig. 1c2), ∆v blue = −585 km s −1 and ∆v red = 535 km s −1 .This is very high for a rotating disk; using the half size of ∼4 kpc, this would imply a dynamical mass of ∼2 × 10 11 M .We do not detect any continuum emitting source at this position, though we note that the partially overlapping foreground G2 may obscure some of the emission (Fig. 1a).Extracting from a r = 0.1 aperture, the observed flux density upper limit is ∼10 −20 erg s −1 cm −2 Å −1 at 2.5 µm (i.e., rest frame V band).This then gives an upper limit to the absolute V band magnitude of ∼−17 which could be a galaxy with M ∼ 10 8 to 10 9 M (e.g., Weaver et al. 2022).This corresponds to a M /M dyn.∼ 10 −3 which is very small for a galaxy at z ∼ 3.5.The southern region shows lower [O i]/Hα and [S ii]/Hα than the rest of the nebula.It also has the potential to have subsolar metallicity (Sect.4.2).Though our estimates have large uncertainties, we cannot exclude the possibility of it being a late stage merger.A more detailed kinematic study is required to unveil the nature of this region with complex kinematics.

Summary and future
In this work, we present the first NIRSpec IFU view of the ionized ISM in the vicinity of the z 3.59 radio-loud AGN 4C+19.71.With unprecedented resolution and sensitivity, we study the gas morphology, ionization states and preliminary kinematics within 10 kpc of the AGN.This makes it possible to study the role of the jet and radiative feedback near cosmic noon and to compare to low-z AGN and theoretical models.We find that the radiation from the hidden quasar dominates the ionization of the ISM of 4C+19.71which resembles the scenario of other powerful quasars at cosmic coon.The radiatively driven outflow is only found within 5 kpc of the AGN and is inefficiently (∼10 −5 ) coupled to the central L bol ∼ 10 47 erg s −1 quasar even at the nucleus.Combining with ground-based studies, we conclude that the jet kinetic energy takes a leading role in the feedback of this HzRG at z ∼ 3.5.Line ratio diagnostics infer the existence of a ionization cone along the jet axis.The AGN ionizing photons illuminate the filamentary extended ISM along the cone.The observed morphology may be the "intrinsic" shape of the dense gas around this relatively massive galaxy.These conclusions are based on the first source from our sample.The picture of the feedback from the most powerful radio sources will be clearer with our full sample which will allow kinematics studies of the ionized ISM of radio-loud AGN with very different jet morphologies.For the manual correction of the astrometry, we matched G1 positions (Section 3) measured from the ∼ 2.7µm continuum image collapsed from our NIRSpec IFU cube and the HST/WFPC2 (Whitmore et al. 2016) image.This was done by matching the centroids of the fit 2D Gaussian profile.We caution that although we updated the astrometry of the HST image using Gaia DR3 (Gaia Collaboration 2023), the uncertainty of the WCS position is ∼ 0.06 .This is estimated based on the standard deviation of the offsets of the cross matched Gaia targets in the HST FoV.This uncertainty (∼ 1.2 given the unknown direction) is marked by the size of the plus marker in Fig. 1d.We further note that the HST image was observed in 1997.The proper monition of the stars were taken into account during the correction by adopting the reported epoch=J2000 coordinates in the Gaia DR3 catalog.Another source of uncertainty could be the intrinsic continuum emission offset from different bands.The HST/WFPC2 image was taken with the F702W filter, which corresponds to 2520 Å in the rest frame of G1.The NIRSpec IFU data covers ∼ 0.57 − 1.15µm for G1.We examined the continuum centroid offset of G1 by extracting narrowband images at the blue and red ends of the data cube.A ∼ 0.09 shift is seen.
We present the offsets of these corrections with respect to the manual shift (aligned with HST) used in the analysis in Table A .1. Fig. A.1 presents the visual inspection of the WCS correction by comparing the results from different methods.We mark the HST positions of the three continuum emission objects with green circles.The radius of each circle is 0.12 which can be used as the spatial uncertainty of the HST object.We see that Correction2+− and Correction2−+ have ∼ 0.1 offsets from the manual shift (or HST WCS).Given that there is a ∼ 0.1 uncertainty in the HST-NIRSpec alignment, we decide that the IFU WCS offset can be solved by correcting the systematic error.Since the direction of the error is not well-understood, we use the manual match in this paper.

Fig. 1 .
Fig. 1.NIRSpec [O iii] color composite and example spectra of [O iii] doublet.(a) Three-color composite image of narrowband [O iii]5007.We overlay the VLA 4.7 GHz jet hot spots in dark purple contours.Blue dots indicate the position of foreground galaxies (see text).FWHM of the PSF is marked at the left corner.(b) "Monochrome" narrowband images of [O iii]5007 collapsed at blue wing (b1), systemic redshift (b2), and red wing (b3), respectively.Yellow contours show dust continuum emission from the ALMA data.(c) Example spectra of continuum-subtracted [O iii]4959, 5007 extracted at four different spatial locations normalized to their peak flux density.We fit the spectra using q3dfit with up to three kinematic components.The light purple curve indicates the overall fit.The individual components are shown in blue, green, and red with a negative offset from the zero level.(d) Zoom-in view of the central 0.5 × 0.5 region of the narrowband image (gray box in panel a).Dark purple plus + sign marks the position of the continuum emission of the radio galaxy determined from the NIRSpec cube while the yellow triple-spike triangle shows peak position of the ALMA 400.3 GHz dust.The sizes of the markers represent the position uncertainty, 1.2 and 0.04 for radio galaxy and ALMA dust respectively.Black box in the insert shows the aperture where the 1D spectrum (Fig. 2) is extracted.
Armus et al. (1998) observed the [O iii] emission of 4C+19.71 with Keck narrowband imaging.They reported the central [O iii](i.e., a similar region as observed in the field of view of the A169, page 3 of 15 NIRSpec IFU) has an elongated morphology along the northsouth direction with two distinctive parts.We show the first zoom-in view of the [O iii] nebula of 4C+19.71 in Fig. 1.We produce three pseudonarrowband images and show their composite in Fig. 1a.The individual images are shown in Figs.1b1-b3 for the blue wing ([−700, −500] km s −1 ), systemic redshift ([−

Fig. 2 .
Fig. 2. Full spectrum extracted at the AGN position without continuum subtraction (black box in Fig.1d).We mark the emission lines detected with S /N 10 in red lines.Lines detected with 3 < S /N < 10 and visually checked to be extended are marked with gray lines.We note that the feature at ∼2.945 µm is likely instrumental due to the leakage flux from the MSA.The central gray shaded region is the detector gap.The top panel is the full view of the spectrum while the zoom-in view of the faint lines are shown in the bottom panels.The magenta curves show the best fit of the spectrum (Appendix C).
spectrum is dominated by the (possibly scattered) light of the hidden AGN.This indicates that the stellar population has a subdominant contribution in the regions near the AGN.The relative emission line fluxes are comparable for 4C+19.71and Cygnus A at least from [O ii]3726, 3729 to [O iii].
brightest emission lines: [O ii]3726, 3729, Hβ, [O iii]4959, 5007, [O i]6300, [N ii]6548, 6583, Hα, and [S ii]6716, 6731.Theselines are fit simultaneously with the same velocity shift and line width for the corresponding components ("kinematically tied", e.g.,Zakamska et al. 2016).FollowingVayner et al. (2023), we ran the fit with 3, 2, and 1 as the maximum number allowed Gaussian components.We performed a visual inspection to determine the best fit for each spaxel.In the following discussions, we sum the fluxes of the doublet [O ii]3726, 3729 and [S ii]6716, 6731, respectively.We use the fluxes of the individual [O iii]5007 and [N ii]6583 lines, respectively.Hereinafter, we omit the wavelength when referring to the line name.The line ratio maps and diagnostics based on the fit are shown in Figs.4 and 5, respectively.As mentioned before, a kinematical study of the ionized ISM is beyond the scope of this paper; we only derive the line ratios of the integrated line fluxes of all velocity components.
are corrected for the dust attenuation based on this E(B − V) map.In Figs.4b,c, we show respectively the ratio maps of [O i]/Hα and [S ii]/Hα based on the integrated line fluxes of all components.The following line ratios are in logarithmic scale, unless otherwise specified.We group the spaxels in the FoV into five regions (white dashed lines in Figs.4b,c), following the apparent geometry of the AGN ionization cone, with the nucleus added as an individual region.The southern part of the ISM is distinctive as having lower [O i]/Hα and [S ii]/Hα values.In the extended filamentary part of the nebula toward the north, there is a large scatter with the presence of both higher and lower [O i]/Hα and [S ii]/Hα values.In general, the east and west regions have higher [O i]/Hα and [S ii]/Hα.To simplify, we combine the east and west regions and refer to it as EW.We inspect our results in line ratio diagnostic diagrams and present in Fig. 5.In addition to [O i]/Hα (Fig. 5b) and [S ii]/Hα (Fig. 5c), we present the [N ii]/Hα-[O iii]/Hβ (Fig. 5a) and [O ii]/[O iii]-[O iii]/Hβ diagnostics (Fig. 5d).The EW region (red diamonds) occupies different locations on the diagnostic diagrams than the other regions, with higher [O i]/Hα, [S ii]/Hα, and [O ii]/[O iii] values.The EW region is also distinctive with a lower [O iii]/Hβ compared to the other regions.Although the southern spaxels (blue hexagons) are clustered toward lower [N ii]/Hα, [O i]/Hα, [S ii]/Hα, and [O ii]/[O iii] values, we find

Fig. 3 .
Fig. 3. Spectral comparison of 4C+19.71 to Cygnus A (Ogle et al. 1997).Blue histogram shows the same 1D spectrum as in Fig. 2. The spectrum of Cygnus A is shown in black which is extracted at its nucleus position after Galactic reddening correction and subtraction of a elliptical galaxy template.Both spectra are normalized to the peak flux density of their Hγ.The bottom panel is a zoom-in version of the region between the gray dotted lines in the upper panel.We marked the emission lines detected for 4C+19.71with S /N 10 as in Fig. 2. The inset provides a further zoom-in around Hβ-[O iii] complex.

Fig. 4 .
Fig. 4. Dust attenuation (a) and line ratio maps (b), (c) constructed based on the fit fluxes integrated for all components.The white and gray contours indicate the 5σ and 15σ [O iii] surface brightness levels of the un-smoothed systemic narrowband image, respectively (Fig.1).Dark purple lines show the directions toward the radio jet hot spots.Black triangle marks the position of the AGN while magenta contours trace the ALMA continuum (Fig.1d).The dashed white lines and circle in panel b and c divide the FoV into five regions for line ratio diagnostic analysis (Fig.5).FWHM of the JWST PSF is marked at the left corner.

Fig. 5 .
Fig. 5. Line ratio diagnostic diagrams based on fitting and corrected for dust attenuation.Each smaller mark corresponds to one spaxel in the FoV.The larger symbols represent the median.Data points with the same symbol are the same region marked as Fig. 4: Nuc.-nucleus (circle with r = 0.25 ), N -north, S -south, EW -east+west.The typical 1σ error is shown at the left corner.We show the 2 Z dusty radiationpressure dominated photonionization models from Groves et al. (2004a) in black and mark the direction of increasing ionization parameter, log U ∈ [−3.0, −2.0, −1.0, 0.0], and power-law index, α ∈ [−2.0, −1.7, −1.4,−1.2].The similar model with 0.5 Z is shown in panel c in cyan.Purple grids are the "shock+precursor" models with per-shock density 100 cm −3 , (shock velocities of 200 ≤ v shock ≤ 1000 km s −1 in steps of 100 km s −1 and magnetic parameters of 0.01 ≤ b = B/n 1/2 ≤ 100 µG cm −3/2 in steps of 1 dex, Allen et al. 2008).The directions of increasing v shock and b are indicated with purple arrows.We include the empirical classification from Kewley et al. (2006) in gray lines in panels a-c.

Fig. A. 1 .
Fig. A.1.Narrowband continuum images (observed frame 2.721 − 2.778µm) with WCS corrected with different methods.The method is labeled at the top right corner of each panel.The three green circles mark the positions of the host galaxy of 4C+19.71(RG), G1, and G2.The radius of each circle is 0.12 which can be used roughly as the uncertainty of the HST WCS.The white dot in panel (b) marks the centroid of G1 using the updated CRDS (jwst_1093.pmap).

Table 1 .
Fit kinematic results from the 1D spectrum at nucleus.Central comp.Red comp.Blue comp.

Table 2 .
Median line ratio values for different regions in Table A.1.Offset between different astrometry and manual shift.