© The Authors 2024

1. 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⁷ yr), observations do indeed indicate that it can still significantly impact the host galaxy (e.g., total energy output by the jet ≳10⁶⁰ 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⁻¹, 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¹¹ 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. 2007, 2017a,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, 2007; Nesvadba et al. 2006, 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, 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. 2023, 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₀ = 70 km s⁻¹ Mpc⁻¹ 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.16 GHz), 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.

2. 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⁻¹ (Falkendal et al. 2019). Finally, there is significant molecular gas in the host (2.5 × 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⁻¹. 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 (IRS²) read-out pattern. The total on-source 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¹ (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.

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.

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.3 GHz 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.

3. 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.

3.1. Morphology of the extended ionized nebula

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 NIRSpec IFU) has an elongated morphology along the north-south 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⁻¹), systemic redshift ([−100, 100] km s⁻¹), and red wing ([500, 700] km s⁻¹), 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.71 has 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.

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)². 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⁻¹ 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⁻¹. 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 jet-gas 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⁻¹).

3.2. 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 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)³.

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.

3.3. 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⁻¹. 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. 2004b, 2006; Kewley et al. 2001, 2006, 2013; Wylezalek et al. 2017).

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).

We extracted a spectrum from a 0.2″ × 0.2″ or 1.4 × 1.4 kpc² (4 × 4 pix²) 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⁻¹) 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.

4. Discussion and summary

4.1. 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²⁸ W Hz⁻¹, 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² and is similar to the aperture used for 4C+19.71 (1.4 × 1.4 kpc²). 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.

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.71 with S/N ≳ 10 as in Fig. 2. The inset provides a further zoom-in around Hβ–[O III] complex.

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 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.71 and Cygnus A at least from [O II]3726, 3729 to [O III]. We conclude that, to first order, AGN photoionization mechanisms dominate the inner parts of these two galaxies. The 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.

4.2. 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. 2004b, 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 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. These lines are fit simultaneously with the same velocity shift and line width for the corresponding components (“kinematically tied”, e.g., Zakamska et al. 2016). Following Vayner 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.

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. 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 radiation-pressure 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 ≤ vshock ≤ 1000 km s−1 in steps of 100 km s−1 and magnetic parameters of 0.01 ≤ b = B/n1/2 ≤ 100 μG cm−3/2 in steps of 1 dex, Allen et al. 2008). The directions of increasing vshock and b are indicated with purple arrows. We include the empirical classification from Kewley et al. (2006) in gray lines in panels a–c.

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 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 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 radiation-pressure 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_Hc), 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 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 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⁻³ (e.g., De Breuck et al. 2000). For the shock velocities, we use higher values (v_shock ≥ 200 km s⁻¹) 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⁻³ are also physical for AGN (Liu et al. 2020) which occupies a subset of the line ratio space spanned by n = 100 cm⁻³ 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).

4.3. Inefficient quasar-driven outflow at nucleus

Our NIRSpec IFU observations enable studying outflows on scales ≲20 kpc. VLT/SINFONI observations captured more extended outflowing gas (Collet et al. 2016; Nesvadba et al. 2017a,b).

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⁻¹) 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⁻¹ and Δv_blue = −386 km s⁻¹). We follow Vayner et al. (2024) and use the following equation to calculate the ionized outflow gas mass (Osterbrock & Ferland 2006):

$$\begin{array}{c}\hfill M_{\mathrm{ion}}=1.4\left(\frac{m_{\mathrm{p}}{L}_{\mathrm{H}\alpha }}{j_{{\mathrm{H}}_{\alpha }}{n}_{\mathrm{e}}}\right),\end{array}$$(1)

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⁻³ from Nesvadba et al. (2017a). Given that the emissivity, j_Hα, is not sensitive to n_e, we use the value of 2 × 10⁻²⁵ erg cm³ s⁻¹ for T_e ≃ 15 000 K (Vayner et al. 2024). Based on Eq. (1), we find M_ion = 10⁵ and 10⁶ M_⊙ for the red and blue components, respectively. We then estimate the outflow rate using the following equation (Vayner et al. 2024):

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{ion}}=\frac{M_{\mathrm{ion}}{v}_{\mathrm{ion}}}{\mathrm{\Delta }R}·\end{array}$$(2)

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⁻¹. The outflow momentum flux Ṗ_ion and outflow kinetic luminosity $L_{\mathrm{ion}}^{\mathrm{kinetic}}$ are calculated from

$$\begin{array}{cc}\hfill {\dot{P}}_{\mathrm{ion}}^{\mathrm{outflow}}& ={\dot{M}}_{\mathrm{ion}}\times v_{\mathrm{ion}},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill L_{\mathrm{ion}}^{\mathrm{kinetic}}& =\frac{1}{2}{\dot{M}}_{\mathrm{ion}}\times v_{\mathrm{ion}}^2.\hfill \end{array}$$(4)

We calculate values for ${\dot{P}}_{\mathrm{ion}}^{\mathrm{outflow}}$ and $L_{\mathrm{ion}}^{\mathrm{kinetic}}$ of 1.5 × 10³⁴ dyne and 8.0 × 10⁴¹ erg s⁻¹, 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_{\mathrm{AGN}}^{\mathrm{IR}}={10}^{10.91}{L}_{\odot }$, and assume a conversion factor, ${\kappa }_{\mathrm{AGN}}^{\mathrm{bol}}=6$, the bolometric luminosity of 4C+19.71 is estimated to be L_bol ∼ 2.5 × 10⁴⁷ erg s⁻¹ (Drouart et al. 2014; Falkendal et al. 2019). Then the coupling efficiency of this potential radiatively driven outflow, $L_{\mathrm{ion}}^{\mathrm{kinetic}}/L_{\mathrm{bol}}$, is ∼10⁻⁵. Using NIRSpec IFU, Perna et al. (2023) and Vayner et al. (2024) reported higher coupling efficiency for z ∼ 3 quasars (0.02 and 10⁻³, respectively). Veilleux et al. (2023) found a relatively low coupling efficiency using the same instrument, 1.8 × 10⁻⁴, 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⁻³ to 10⁻², Nesvadba et al. 2008, 2007, 2017a,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⁻¹) 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⁻⁵ 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⁴⁷ erg s⁻¹) 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.71 seem more quiescent than in other HzRGs (Nesvadba et al. 2017a; Wang et al. 2023), a gas kinetic energy luminosity of 10⁴⁵ erg s⁻¹ 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_{\mathrm{ion}}^{\mathrm{kinetic}}$) 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⁻⁴, which is still ∼1 dex lower than the kinetic efficiency from Nesvadba et al. (2017a).

4.4. 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 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 flux-limited 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.71 has one of the most powerful radio jets, P_jet ∼ 10⁴⁷ erg s⁻¹ (using log(P_1.4 GHz/W Hz⁻¹) = 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.

4.5. 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 “high-velocity” region is more complicated. In Fig. 1c2, we clearly find that there are three velocity components with the “high-velocity” 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⁻¹ and Δv_red = 535 km s⁻¹. This is very high for a rotating disk; using the half size of ∼4 kpc, this would imply a dynamical mass of ∼2 × 10¹¹ 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⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹ 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⁸ to 10⁹ M_⊙ (e.g., Weaver et al. 2022). This corresponds to a M_⋆/M_dyn. ∼ 10⁻³ 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.

4.6. 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.71 which 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⁻⁵) coupled to the central L_bol ∼ 10⁴⁷ erg s⁻¹ 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.

Acknowledgments

This work uses the NASA’s Astrophysics Data System and a number of open source software other than the aforementioned ones such as Jupyter notebook (Kluyver et al. 2016); matplotlib (Hunter 2007); SciPy (Virtanen et al. 2020); NumPy (Harris et al. 2020); Astropy (Astropy Collaboration 2018); LMFIT (Newville et al. 2016). We thank Daizhong Liu and Markus Hundertmark for the useful discussion of the WCS alignment. This work is based in part on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program JWST-GO-01970. Support for program JWST-GO-01970 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. D.W. and W.W. acknowledge support through an Emmy Noether Grant of the German Research Foundation, a stipend by the Daimler and Benz Foundation and the MERAC foundation and a Verbundforschung grant by the German Space Agency. N.L.Z. and A.V. are supported in part by NASA through STScI grants JWST-ERS-01335 and JWST-GO-01970. The work of D.S. was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2021.1.00576.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

References

Appendix A: Astrometric correction

The JWST Fine Guidance Sensor (FGS) is responsible for guiding the pointing, and there is no known issue in the guiding process. However, there is a known systematic WCS offset of the NIRSpec IFU⁴. This is due to poorly registered IFU detector coordinates on the JWST focal panel with respect to the reference V2 and V3 axes. This corresponds to 0.173″ and 0.178″ in the V2 and V3 direction⁵, respectively (Help Desk, private communication through). This offset projection to sky (RA, Dec) coordinates could then be calculated using pysiaf by assuming the targeted RA and Dec of the IFU pointing is correct. The projection matrix for converting V2, V3 coordinates to the sky plane can be set by pysiaf.utils.rotations.attitude_matrix using the old (V2, V3) reference coordinates, target RA and Dec, and angle between the north and V3. This angle is measured from north to east with a value of 83.0652845 deg for our observation. It is calculated using the angle between the IFU with respect to the North (from the header) and the fixed angle between the IFU and V3. There is an additional ∼3″ systematic rotational angle error between the IFU plane and V3 axis. Since this will only result in a secondary (∼0.007″) offset, we neglect it here. The final systemic offsets of the pointing, from updated position to the original position, are ΔRA = 0.21″ and ΔDec= − 0.15″. We refer to this correction as “Correction1". This error has been reported as partially corrected using the up-to-date context file. We test the new pipeline assigned WCS using jwst_1093.pmap and find it matches closely with our “Correction1". We report an offset of |ΔRA| = 0.09″ and |ΔDec| = 0.07″ using the fit 2D Gaussian centroids of G1 at ∼2.7μm. Hence, we consider this systematic WCS error has been fixed. However, other source of uncertainty is still presented and still required manual alignment.

In addition, there is a further offset due to the zero-point shift of the NIRSpec detector and the FGS which adds another term to the aforementioned values in the V2V3 plane. The shifts are -0.016″ and 0.111″ in V2 and V3 coordinates, respectively (Help Desk, private communication). However, a caution is noted that the direction of this shift is not well-understood. We assumed this direction and followed the aforementioned procedure using pysiaf. This resulted in the ΔRA = 0.09″ and ΔDec= − 0.14″. We refer to this correction as “Correction2++". Since the direction is uncertain, we produced further corrections “Correction2+−", “Correction2−+", and “Correction2−−". The “+" and “−" signs indicate which direction we adopted for the offset with respect to the original direction suggested by, the Help Desk.

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. 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).

Appendix B: Foreground galaxies

We report in Sect. 3 that two foreground galaxies have been identified with continuum and emission lines: G1 to the west of the radio galaxy and G2 to the east. There we present the spectra of G1 and G2. Their positions are marked with blue dots in Fig. 1 which are the centers of their extraction apertures. For both galaxies, we used square apertures with the size of 0.25″ × 0.25″ (∼5 × 5 pix²) and summed each spectrum from the spaxels inside the aperture. For partially covered spaxels, we multiply the area by the fractional overlap with the extraction aperture. The aperture for G1 was centered at its continuum peak position. As for G2 which is only partially captured in our FoV, we placed the aperture to be as close as possible to its continuum emission peak while minimizing the impact from the detector edge. In Fig. B.1 and B.2, we show the spectra extracted for G1 and G2, respectively. The redshifts of G1 and G2 are determined based on their Hα lines and have values of z_G1 ≃ 1.786 and z_G2 ≃ 1.643. For each foreground galaxy, at least two more emission lines are identified based on their redshifts. We use thick black vertical lines to mark the spatially and spectrally detected emission lines. The expected positions of emission lines near Hα based on the derived redshift of the targets are marked in gray lines if not clearly detected. We note that the [S II]6716,6731 doublet of G2 overlaps with the [Ne III]3869 line of 4C+19.71 in wavelength (orange line in top panel of Fig. B.2). The spatial check confirms that there are distinct emission peaks separated by ∼0.4″.

Fig. B.1. Spectrum extracted at the position of foreground galaxy, G1, without continuum subtraction. Upper panel: Zoom-in view of the spectrum around Hα. Middle panel: Zoom-in view of the spectrum around [S III]9071. Bottom panel: Zoom-in view of the spectrum around [S III]9533. We mark the lines detected both spectrally and spatially in thick black lines. The expected positions of the lines around Hα are marked in gray lines.

Fig. B.2. Spectrum extracted at the position of foreground galaxy, G2, without continuum subtraction. Upper panel: Zoom-in view of the spectrum around Hα. Middle panel: Zoom-in view of the spectrum around [S III]9533. Bottom panel: Zoom-in view of the spectrum around He I10833. Same line-marking conventions are adopted as Fig. B.1. We further mark in orange the position of [Ne III]3869 emission from 4C+19.71 which is near the [S II]6716,6731 doublet of G2 (upper panel).

Appendix C: Emission line list

We fit the 1D spectrum presented in Section 3.3 (Fig. 2) using q3dfit. The continuum is modeled by a third order polynomial and the emission lines (S/N > 10) are fit with Gaussian profiles. We set the maximum number of the Gaussian components to be three for the fit. The centroid (width) of each corresponding Gaussian component is set to be the same (Vayner et al. 2023). The fit line fluxes are shown in Table C.1. For the doublet [O II]3726,3729, we set the flux ratio to be 1:1 and report the summation since they are spectrally unresolved. For the low S/N lines (Sect. 3.3), we report 3σ upper limits.

Appendix D: Systemic redshift based on ionized gas

The systemic redshift (velocity zero) used in this work is based on the [C I](1-0) (z = 3.5892, Kolwa et al. 2023). Given its relatively faint flux and large observed beam size, it is unclear whether this traces the “authentic” redshift of the host galaxy. The ideal tracer would be the stellar absorption lines (e.g., Ca II HK). However, the continuum emission near the nucleus is AGN dominated (Sect. 4.1) for our target. Nevertheless, we report the [O III] flux-weighted redshift based on the fit “low-velocity” components from the entire FoV (Sect. 4.2). We define the “low-velocity” components as the ones with |Δv|< 300 km s⁻¹ (e.g., Fu & Stockton 2009). The redshift is z_[O III] = 3.5908 ± 0.0017 which corresponds to 108 ± 114 km s⁻¹ with respect to the [C I](1-0) redshift. The z_[O III] and z_[C I] are consistent within the relatively large uncertainty.

In our calculation in Sect. 4.3, we assume the central component is the systemic component. Its Δv_central = 238 ± 10 km s⁻¹ agrees with the systemic velocity estimated here from the entire FoV.

Appendix E: Estimation of electron density

We use the fiducial value of n_e in the calculation in Sect. 4.3 (Nesvadba et al. 2017a). The n_e for the red and blue components can also be estimated based on the fit [S II]6716 and [S II]6731 lines ratios (0.50 and 0.44 for red and blue components respectively; Table C.1). This result in a n_e ∼ 8 × 10³ cm⁻³ and ∼8 × 10⁴ cm⁻³ for the red and blue components, respectively (e.g., Sanders et al. 2016). These are much larger than the fiducial value, n_e = 500 cm⁻³, from Nesvadba et al. (2017a). However, the fit line ratios for these two kinematic components are in the range where they are not sensitive to n_e. If we nevertheless use n_e based on our estimation, the radiative coupling efficiency becomes ∼10⁻⁷, and continues to support the conclusion is valid that the radiative-driven outflow is inefficient in 4C+19.71.

We report that the n_e of the central component is ∼800 cm⁻³ based on its [S II]6716 and [S II]6731 ratio (0.89, Table C.1; Sanders et al. 2016) which is the same order of magnitude as the value reported by Nesvadba et al. (2017a). Further analysis is required to understand the electron temperature and density at the center.

All Tables

All Figures

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.

In the text

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).

In the text

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.71 with S/N ≳ 10 as in Fig. 2. The inset provides a further zoom-in around Hβ–[O III] complex.

In the text

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.

In the text

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 radiation-pressure 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 ≤ vshock ≤ 1000 km s−1 in steps of 100 km s−1 and magnetic parameters of 0.01 ≤ b = B/n1/2 ≤ 100 μG cm−3/2 in steps of 1 dex, Allen et al. 2008). The directions of increasing vshock and b are indicated with purple arrows. We include the empirical classification from Kewley et al. (2006) in gray lines in panels a–c.

In the text

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).

In the text

Fig. B.1. Spectrum extracted at the position of foreground galaxy, G1, without continuum subtraction. Upper panel: Zoom-in view of the spectrum around Hα. Middle panel: Zoom-in view of the spectrum around [S III]9071. Bottom panel: Zoom-in view of the spectrum around [S III]9533. We mark the lines detected both spectrally and spatially in thick black lines. The expected positions of the lines around Hα are marked in gray lines.

In the text

Fig. B.2. Spectrum extracted at the position of foreground galaxy, G2, without continuum subtraction. Upper panel: Zoom-in view of the spectrum around Hα. Middle panel: Zoom-in view of the spectrum around [S III]9533. Bottom panel: Zoom-in view of the spectrum around He I10833. Same line-marking conventions are adopted as Fig. B.1. We further mark in orange the position of [Ne III]3869 emission from 4C+19.71 which is near the [S II]6716,6731 doublet of G2 (upper panel).

In the text