© The Authors 2024

1. Introduction

Active galactic nuclei (AGNs) residing in the centres of galaxies are prime suspects with respect to the regulation of star formation in galaxies. Silk & Rees (1998) were first to recognise that a small fraction of the immense energy output of AGNs would suffice to unbind most of the interstellar gas from massive host galaxies – provided that the AGN energy can efficiently be injected into the gas of the host galaxy over kiloparsec scales. This would resolve several of the long standing open questions of extragalactic astrophysics, including, for instance, the old stellar populations, low gas content, star formation rates in massive galaxies, the metal enrichment in the intra-cluster gas of massive galaxy clusters, and the tight relationship between the masses of supermassive black holes and the masses or central velocity dispersions of the bulges of their host galaxies (e.g. Gebhardt et al. 2000; Ferrarese & Merritt 2000; Springel et al. 2005; Croton et al. 2006; Forman et al. 2005; Kirkpatrick et al. 2009; Bharadwaj et al. 2015; Schaye et al. 2015; Weinberger et al. 2018; Werner & Mernier 2020).

Many observational and theoretical studies have been addressing the question of how the AGN energy is injected into the surrounding gas and how this affects star formation. Most of these studies have focused on the rapid gas removal through galaxy-wide outflows. However, although many observations have revealed that outflows of molecular, atomic, or ionised gas are common in the host galaxies of powerful AGN and more likely related to the AGN rather than concomitant star formation (e.g. Fiore et al. 2017; Veilleux et al. 2013), uniquely establishing a physical link between the AGN energy injection and the suppression of star formation is still a challenge. For example, many host galaxies of powerful AGN (and even AGN-driven outflows) do not show a concomitant deficit in star formation (Stanley et al. 2017; Ellison et al. 2016).

A remarkable exception is a subclass of (mostly) powerful radio galaxies. In about 30% of nearby (z < 0.2) galaxies from the 3CRR catalogue (Laing et al. 1983) observed with the Infrared Spectrograph (IRS) on board the Spitzer Space Telescope, Ogle et al. (2007) and Ogle et al. (2010) identified bright mid-infrared line emission from luminous, warm molecular hydrogen (molecular hydrogen emission galaxies, MOHEGS; see also Dicken et al. 2012), while typical AGNs and star-formation tracers, such as forbidden lines, polycyclic aromatic hydrocarbons (PAHs), or the mid-infrared dust continuum, are either weak or undetected.

Due to its very simple and highly symmetric structure, the homonuclear H₂ molecule does not have a permanent electric dipole moment; therefore, it radiates only at high gas temperatures, above an order of 100 K. This makes H₂ line emission so elusive in many star forming molecular clouds, which have typical temperatures on the order of 10 K. The galaxies found by Ogle et al. (2010) are remarkable, in that their H₂ line luminosities can reach up to about 20% of their total infrared dust luminosity (Ogle et al. 2007), which requires a continuous gas heating mechanism over long timescales (Nesvadba et al. 2010).

These galaxies show also a marked offset from the Kennicutt–Schmidt relationship between molecular gas mass and star formation rate surface densities, as seen in multiple star-formation tracers, including PAHs, cold dust observed with Herschel, the UV continuum, or Hα (e.g. Nesvadba et al. 2010; Alatalo et al. 2015; Lanz et al. 2016; Nesvadba et al. 2021; Drevet Mulard et al. 2023). In extreme cases, their star formation rates can be one to two orders of magnitude lower at a given gas-mass surface density than in normal star-forming galaxies, which might be a signature of suppression of star formation through the AGN.

Ogle et al. (2010) and Nesvadba et al. (2010) hypothesised that the bright H₂ line emission in these galaxies is produced by shocks that are driven into the molecular gas through interactions with the radio jet, as indicated by several line diagnostics comparing the luminosity of warm H₂ lines with indicators of AGN and star-formation activity, such as PAHs or the mid-infrared continuum (e.g. Nesvadba et al. 2011; Guillard et al. 2012). This is also suggested by the enhanced line widths of the warm H₂ lines in these galaxies (Nesvadba et al. 2010; Guillard et al. 2012). A comparison with CO observations of cold molecular gas suggest that at least in some of these galaxies, the warm molecular gas component dominates the overall gas budget (Nesvadba et al. 2010; Ogle et al. 2010), with warm molecular gas masses of up to few × 10⁹ M_⊙ (Ogle et al. 2010; Nesvadba et al. 2010).

Nearly 15 years after the end of the cold Spitzer mission, the James Webb Space Telescope is now again enabling new mid-infrared observations of AGN, at a higher sensitivity and spectral and spatial resolution, by several orders of magnitude, but with a somewhat more restrictive spectral range (e.g. Pereira-Santaella et al. 2022; Lai et al. 2022; Armus et al. 2023; Álvarez Márquez et al. 2023). Here, we present new observations of one of the best studied MOHEGs, 3C 326 N, at z = 0.09, which we obtained through open-time observations with the integral-field units (IFUs) of the JWST/Near-Infrared Spectrograph (NIRSpec) and Mid-Infrared Instrument (MIRI) in Cycle 1. In brief, 3C 326 N is a giant FRII radio galaxy with a radio power of 2.5 × 10²⁶ W Hz⁻¹ at 327 MHz, and two Mpc-sized lobes. Willis & Strom (1978) estimated a spectral age for the radio source of 200 Myr. We note that 3C 326 N has a nearby neighbour, 3C 326 S, which makes identifying the source of the radio emission challenging; in addition, both have radio cores. Rawlings et al. (1990) argued that 3C 326 N is the better candidate, given its brighter stellar continuum and higher stellar mass of 3 × 10¹¹ M_⊙ (Nesvadba et al. 2010), and since the only galaxy of the two that has [OIII]5007 line emission consistent with the usual relationship between [OIII]5007 luminosity and radio power.

3C 326 N was observed with Spitzer/IRS by Ogle et al. (2007), who identified an unusually bright H₂ line emission from the pure-rotational lines between 5 and 30 μm, whereas many typical AGN lines, and in particular PAHs, have either been very faint or unobserved. The galaxy was further studied by Nesvadba et al. (2010), who obtained IRAM CO(1–0) observations and found that the warm gas is most likely heated by shocks and that the warm molecular gas dominates the overall molecular gas budget, with a warm molecular gas mass of 2 × 10⁹ M_⊙. They also found Na D absorption indicating a (likely jet-driven) outflow of about 30 − 40 M_⊙ yr⁻¹ in neutral gas and a terminal velocity of −1800 km s⁻¹. 3C 326 N has a very low star formation rate ($\mathrm{SFR}=0.{087}_{-0.046}^{+0.106}$ M_⊙ yr⁻¹ from Herschel dust photometry; Lanz et al. 2016) and a radiatively weak AGN (L_X = 10^40.63 erg s⁻¹; Lanz et al. 2016), so that neither star formation nor AGN are powerful enough to heat the molecular gas (Nesvadba et al. 2010). Moreover, the cosmic ray flux required to heat the observed amount of molecular gas would be high enough to destroy the H₂ molecules (Nesvadba et al. 2010), which leaves mechanical heating through shocks driven by the radio source as the only plausible gas heating mechanism. This was subsequently confirmed through remarkably broad line widths (FWHM = 650 km s⁻¹) of ro-vibrational emission lines of H₂ seen with the VLT/SINFONI near-infrared IFU (Nesvadba et al. 2011), as well as through very high line ratios of H₂ to Paα, which are also a good indicator of shock (Puxley et al. 1990). These authors also suggested that the line broadening may be indicative of high gas turbulence, which (if it goes down to the usual scales of giant molecular clouds at ∼100 pc) may suppress the star formation by the molecular gas. By comparing the energy injection that can plausibly be injected through star formation, AGN radiation, or the radio jets, Nesvadba et al. (2010) concluded that only the radio source is powerful enough to explain the observed gas kinematics.

In this work, we aim to provide a more in-depth look, via novel modelling methods and improved data quality and wavelength coverage, into the molecular and ionised gas structures of 3C 326 N that have initially been probed in these previous works. We provide a rich set of modelled spectral features, which we used to deduce a reliable systemic redshift and investigate the kinematics and heating mechanism of the molecular gas. We also investigated the potential of the ro-vibrational lines that are produced by hot > 1000 K gas (comprising only a small fraction of the total gas mass) as tracers of turbulence in place of pure rotational lines which are produced by the bulk of the molecular gas at cooler ∼100 K temperatures. The rotational lines are red-shifted out of the MIRI bandpass at z ≥ 1.5, whereas the ro-vibrational lines are accessible up to the Epoch of Reionisation if they are bright enough.

The paper is organised as follows. In Section 2, we describe our observations and data reduction for the JWST and ancillary datasets obtained with ALMA and the JVLA. We then explain our methodology and Bayesian fitting routine in Section 3, before describing our analysis of the stellar and dust continuum, radio morphology, and kinematics of the molecular and warm ionised gas in Section 4. In Section 5, we discuss our results in the context of multi-phase molecular gas properties and star formation. We then use this source to investigate what our results imply about observations of distant galaxies in the early Universe, where rotational lines are not observable with MIRI, in Section 6. We summarise our results in Section 7.

Throughout the paper, we use the flat Planck Collaboration XIII (2016) cosmology, where H₀ = 67.8 km s⁻¹ Mpc⁻¹ and Ω_M = 0.308. In this cosmology, the luminosity distance to 3C 326 N at z = 0.08979 is D_L = 423.6 Mpc, while the angular size-distance is D_A = 356.8 Mpc. A total of 1.729 kpc is projected onto 1 arcsec.

2. Observations and data reduction

2.1. NIRSpec and MIRI imaging spectroscopy

3C 326 N was observed with the integral-field units (IFUs) of the Near-Infrared Spectrograph (NIRSpec) and the Medium-Resolution Spectrograph (MRS) of the Mid-Infrared Instrument (MIRI) on the James Webb Space Telescope as part of programme GO1-2162 (PI Nesvadba). The NIRSpec data were obtained on 4 March 2023 through the filter F170LP and with the grating G235H. These data cover a near-continuous spectral range between 1.66 μm and 3.17 μm, with a gap between 2.40 μm and 2.56 μm, due to the physical gap between the two NIRSpec detectors in the focal plane of the telescope. The field of view of the NIRSpec IFU is 3.0″ × 3.0″ and the dispersion through the G235H grating is Δλ = 8.7 × 10⁻⁴ μm, which corresponds to a spectral resolution of 85 − 150 km s⁻¹.

The data were obtained during one visit with 2976 s of on-source observing time at four dither positions, using the dither pattern NRSIRS2RAPID with 50 single-integration groups. We also obtained a background exposure with the same readout pattern, 10 single-integration groups, one dither position, and a total exposure time of 160 s.

The MIRI data were obtained on 28 March 2023 with the four channels of the LONG setting. Channels 1–4 cover spectral ranges between 6.53 μm and 7.65 μm, 10.02 μm and 11.7 μm, 15.41 μm and 17.98 μm, and 24.19 μm and 27.9 μm, respectively. The spectral resolving power is R = 1330 − 3610, corresponding to a spectral resolution between 83 and 226 km s⁻¹. The field of view of the MIRI IFU is between 3.2″ × 3.7″ in Channel 1 and 6.6″ × 7.7″ in Channel 4.

We used the MRSLONG and MRSSHORT detectors with the LONG(C) grating, FASTR1 readout pattern, 35 single-integration groups, and four dither positions that were optimised for extended sources. The on-source exposure time was 389 s per detector. A second visit was used to observe a background frame with the same setup directly following the on-source observations.

Both the NIRSpec and MIRI data were reduced with version 1.12.4 of the official pipeline PYTHON package: JWST. We used the CRDS context collection JWST_1146 for static calibration file association. Both datasets require no additional data cleaning beyond the pipeline’s performance. In summary, the raw files were processed with Detector1Pipeline to correct for detector level effects. The resulting ramp files were passed to the Spec2Pipeline, which performs the WCS correction, flat field correction, background correction, and flux calibration. Furthermore, the pipeline applies an MSA imprint and pathloss correction for the NIRSpec data as well as a fringe correction for the MIRI MRS data. We experimentally applied the fringe correction to the NIRSpec data as well but found no improvement. Finally, the data for NIRSpec and MIRI were resampled and combined with the Spec3Pipeline. We used the drizzle combination method to combine and rotate the NIRSpec data and each MIRI channel separately. This pipeline also performed outlier detection and a master background subtraction for the MIRI MRS data.

2.2. ALMA interferometry

We also obtained archival ALMA data in band 3 from programme 2015.1.01120.S, covering the expected wavelength of ¹²CO(1–0) at z = 0.0898, 105.8 GHz. Continuum emission from the nucleus is aptly detected, with a flux density at 115.3 GHz in the rest-frame of S_115 GHz = 1.6 ± 0.01 mJy bm⁻¹. The CO(1–0) line emission is not detected. We used this dataset to derive an upper limit on the line flux, which has RMS = 0.2 mJy bm⁻¹ in 15.4 MHz-wide spectral channels, for a beam size of FWHM = 0.46″ × 0.26″ along the major and minor axes, respectively. The position angle is PA = −13.5°, measured from north towards east.

We followed the prescription of Sage et al. (2007) and Young et al. (2011) by setting σ_line² = Δv²σ²N_l (1.+N_l/N_b), where σ_line is the line-integrated standard deviation of the dataset, Δv, the line width in km s⁻¹, N_l is the sampling of the line, and N_b = 118 is the number of spectral channels from which σ was estimated. Assuming a line width of FWHM = 100 km s⁻¹, a sampling of 43.26 km s⁻¹ (as in our data), and RMS = 0.2 mJy bm⁻¹, we find σ_line = 30.7 mJy km s⁻¹ bm⁻¹. Using Equation (3) from Solomon et al. (1997), we can translate this flux into a line luminosity, by setting $L{\prime }_{\mathrm{CO}}=3.25\times {10}^7{S}_{\mathrm{CO}}\mathrm{\Delta }V{\nu }_{\mathrm{obs}}^{-2}{D}_{\mathrm{L}}^2{(1+z)}^{-3}$, where S_COΔV is the integrated line flux in Jy km s⁻¹, ν_obs is the observed frequency in GHz, D_L is the luminosity distance in Mpc, and z is the redshift. To estimate a molecular gas mass from the luminosity, we adopted the standard conversion factor appropriate for the Milky Way, α_CO = 4.3 M_⊙/ [K km s⁻¹ pc²] (e.g. Bolatto et al. 2013). We found a 3σ upper limit of 1.5 × 10⁸ M_⊙ bm⁻¹. Assuming that CO line emission would extend over the same area also seen with NIRSpec in warm molecular gas, we would expect that the total molecular gas mass would be a factor 5.3 larger, namely, corresponding to a total of M_H2, CO = 8 × 10⁸ M_⊙. This is about a factor 3 lower than the value of M_H2 = 2 × 10⁹ M_⊙ of warm molecular gas previously found with Spitzer (Ogle et al. 2007). This outcome suggests that the majority of the molecular gas in 3C 326 N is indeed warm.

Nesvadba et al. (2010) previously reported weak CO(1–0) line emission detected at 5σ with the IRAM Plateau de Bure interferometer in the D configuration at 5″ beam size. Given the new, deeper observations available now from ALMA with a much larger bandwidth, we consider this previous detection to be spurious. Although the spectrum was carefully extracted from the data cube, the large beam size and narrow bandwidth made it very challenging to isolate the continuum from the line emission, which may explain the faint putative line detection obtained by the best instrument available at that time.

2.3. JVLA interferometry

We also obtained deep radio continuum imaging with the Karl G. Jansky Very Large Array (JVLA) in the A-array in the L, C, and X bands at 1.52, 7.25, and 9.0 GHz through programme 18A−285 (PI Nesvadba). We observed the source in four sessions between 2 March 2018 and 9 June 2018 with 24–27 antennae and a total allocation of 6.4 h. Our main goal was to probe the morphology of the radio core of 3C 326 N and to search for extended emission on scales comparable to the 3″ size of the molecular disk. We used 3C 286 as flux and bandpass calibrator and J1609+2641 as the gain calibrator. In the L-band, we centred the two basebands on 1.264 and 1.776 GHz to obtain 1 GHz of nominal bandwidth. In the C and X-bands, we tuned each baseband to 5.5 and 6.5 GHz, and to 9.5 and 10.5 GHz, respectively, reaching 2 GHz of nominal bandwidth.

Data were reduced in the standard way using CASA Version 5.4.2-5. First, data were flagged, antennae positions corrected, and an initial flux density scaling was derived from our primary calibrator 3C 286 in each band. We then did an initial phase calibration for all bands before solving for antenna-based delays and performing the bandpass calibration in each band, again, using 3C 286. Complex gains were derived from J1609+2641 before we used the flux of our primary calibrator to derive the amplitude gains of J1609+2641 and 3C 326 N in each band. We eventually applied the calibration before reconstructing and cleaning the images, using briggs = 0.5. The final RMS obtained in this way is about 5 mJy bm⁻¹ in the C and X-bands, and 90 mJy bm⁻¹ in the L band, which is strongly affected by RFI. Beam sizes are around 0.2″–0.3″ in the C and X bands, and 1.1″ in L. The central tuning, beam size, position angle, and RMS in each band can be found in Table 1.

2.4. Relative alignment

The JWST/NIRSpec and each channel of the MIRI observations show a small spatial offset with each other and with the ancillary datasets, mainly along right ascension. Relative to the NIRSpec data, MIRI data are offset in RA by −0.7″, 0.58″, and 0.32″ in channel 1, channel 2, and channel 3, respectively. The offsets in Dec are subpixel in MIRI and we do not list the offset for channel 4 due to an extremely tentative continuum detection. The positional offsets are significantly larger than the given 0.1″ pointing accuracy of JWST¹. Therefore, we did not simply use the position on the sky as indicated in the file headers; instead, we used the stellar continuum peak in the NIRSpec data cube to align the JWST data relative to the radio data. In aligning the datasets, we made the simple astrophysical assumption that the radio core falls on top of the peak in the stellar continuum emission, which is also the location of the highest Paα, [FeII]λ1.599, and [FeII]λ1.6440 surface brightness, highlighting the presence of a faint AGN at that position.

3. Emission line and stellar continuum fitting

To model the spectral features of the JWST data, we used a Markov chain Monte Carlo (MCMC) based fitting method. The sampling was performed with the python package EMCEE (Goodman & Weare 2010) and the data file handling was performed by ASTROPY (Astropy Collaboration 2013).

3.1. Emission lines

Unless otherwise stated, the line modelling was performed on one line at a time with the aim to probe differences in line kinematics or morphology (if present). When modelling a line, we cut the cube in the spectral direction such that the modelled wavelength range is λ₀(1 ± 0.01)(1 + z_in), where λ₀ is the lines rest frame wavelength and z_in is our initial redshift guess set to 0.09 based on the measurement of Nesvadba et al. (2011) and visual inspection of the NIRSpec lines in the central 0.5″ of the galaxy.

We reproduced the emission line features using a Gaussian line profile coupled with a linear continuum in frequency. We find that the line shape is often complex and requires multiple Gaussians to explain. Therefore, the line profile can be expressed as:

$$\begin{array}{c}\hfill F_{\nu }(\nu )=m\nu +b+{\displaystyle \sum _{n=1}^N}(1+z_n)A_{n}exp(-0.5{\left(\frac{(1+z_n)\nu -{\nu }_0}{{\sigma }_n}\right)}^2),\end{array}$$(1)

where A_n is the amplitude of Gaussian, n; ν is the frequency; ν₀ is the line rest frame frequency; z_n is the redshift of Gaussian, n; σ_n is the standard deviation of the line, n; m is the continuum line gradient; and b is the continuum intercept. We note that m and b are not adjusted for redshift when modelled. We allowed up to two Gaussians and selected the number used by the Bayesian information criterion (BIC).

When modelling, we measure the goodness of fit from a parameter set using the maximum likelihood + prior probability distribution. Our maximum likelihood is defined as:

$$\begin{array}{cc}\hfill 2ln\left(L\right)& =-{\displaystyle \sum [\frac{{(F_{\nu ,\mathrm{ob}}\left(\nu \right)-F_{\nu }(\nu ,\alpha ))}^2}{\delta }+ln\left(2\pi \delta \right)]},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \delta & ={\mathrm{err}}^2+f^2{F}_{\nu }{(\nu ,\alpha )}^2\hfill \end{array}$$(3)

where err is the data uncertainty; f is the fraction by which the uncertainty is underestimated; and F_ν, ob is the observed flux. We define the prior probability function for all parameters to be uniform between the boundaries and −inf otherwise. We select nonrestrictive boundaries of −5 to 5, −1 to 1 Jy, 0–10 mJy, 0.1–300 GHz, and 0.088–0.0915 for m, b, A_n, σ_n, and z_n, respectively. For z, we defined an additional prior probability, P_z. P_z is a normal distribution such that:

$$\begin{array}{c}\hfill P_z={\displaystyle \prod _{n=1}^N}(\frac{1}{{\sigma }_z\sqrt{2\pi }}exp(-0.5{\left(\frac{z_n-<z>}{{\sigma }_z}\right)}^2)),\end{array}$$(4)

where ⟨z⟩ is the mean z for all N Gaussians in the model and the standard deviation of the z separation probability of σ_z = 0.001. We base this probability distribution on the fact that the separation between Gaussians for one line can differ, but greater separations are less likely to be part of the same line instead of two unrelated features. We chose the σ_z value to be 0.001 to set a 95% chance that a line emitting source is not separated from the average velocity of emitting structures in one spaxel by more than ∼550 km s⁻¹ which is similar to the maximum velocity offset between components for any single line in 3C 326 N from Nesvadba et al. (2011). We find this significantly improves the multi-line fits without being overly constraining on the velocity separation of coincidental projected structures. We also set the rule that z_n < z_n + 1 to prevent the degenerate solution of components switching parameters. These additional rules for z have no effect on the one-component case.

All Gaussian components were then modelled simultaneously for the line and each spaxel was modelled independently. We used a uniform initial distribution between the boundaries of 0–1 mJy and 1–100 GHz for the amplitude and standard deviation, respectively, and a normal distribution around 0.09 ± 0.001 for z. We used 5000 steps and 200 walkers. We then performed a second run, starting from the best fit parameters by likelihood of the first run with a normal distribution with a 0.1% standard deviation for the initial positions. The second run also used 5000 steps and 200 walkers. An example result for a two-Gaussian fit to the H₂1–0 S(3) line for different regions of the galaxy is given in Fig. 1.

Fig. 1. H21–0 S(3) Gaussian line model results compared to data for different spaxels and resulting integrated line emission map. Middle: integrated model flux for H21–0 S(3). Edges: Gaussian model fit to the H21–0 S(3) line for different spaxels of the central image compared to data. The coloured squares in the image represent the position of the plotted spaxel of the same colour. Dashed lines are individual Gaussian components and the solid line is the summed Gaussians. Shaded regions represent the 1σ errors. These spectra area also a good illustration of the high quality of the date and the automatic fitting routine.

When reporting Gaussian line component amplitudes and widths, we account for the instrumental line spread function. These vary slightly between NIRSpec and the different MIRI channels. The resolution for the NIRSpec IFU and MIRI MRS long channels 1, 2, 3, and 4 is: 2700, 3100–3610, 2860–3300, 1980–2790, and 1630–1330, respectively².

3.2. CO band head absorption

To model the CO band head absorption features, we used a template provided by Winge et al. (2009) instead of a Gaussian line profile. The templates are of single stars; to adjust for a stellar population we assume the distribution of velocities is normal. We convolved the template spectrum with a Gaussian with a standard deviation of the stellar velocity dispersion, σ_SVD, and resampled it to the spectral resolution of the NIRSpec spectrum using the python package SPECTRES (Carnall 2017). We then multiply the spectrum by a second-order polynomial continuum to correct for the normalisation applied by Winge et al. (2009). We compared the adjusted template to the data between the observed wavelengths of 2.485 μm and 2.62 μm which contains only CO band head features and no emission features (Fig. 2). Similarly to the emission lines, the instrumental line spread function (R = 2700) is accounted for in the given velocity dispersions.

Fig. 2. Results of CO band head template modelling. Top-left: CO band head velocity map from our template fit. For clarity, while we model all spaxels, here we only plot the velocities for spaxels with an average signal-to-noise ratio of S/N > 10. Top right: FWHM of the velocity dispersion for the CO features, the spaxels with an average S/N < 10 are omitted. Overplotted as coloured squares: Positions of the four spaxels displayed in the bottom four panels. The spaxel positions and labels match those of Fig. 1. Overplotted on both top plots in purple: Contours of the near-IR continuum emission in Fig. 6. Bottom: model fit to the CO features in four selected spaxels. The filled regions are the 1σ uncertainties.

To find the best template to reproduce our data, we compared every template to the highest signal-to-noise spaxel in the NIRSpec dataset. We found that the template star of HD 9138, a K4III star, is a very good reproduction of our data without additional spectral types. We proceed to fit only the HD 9138 template to every spaxel. We generated an initial position for the walkers using a normal distribution around the best fit parameters of the highest S/N spaxel fit with a standard deviation of 1%. We used 2000 steps and 100 walkers when performing the fitting. An example of the model fitting can be seen for four selected spaxels in Fig. 2. It can be seen that the one template can reproduce our observed data very well.

4. Analysis

4.1. Properties of the radio core

The GHz radio morphology of 3C 326 N is compact on scales of a few arcsec in all three bands we observed. As an example, the 10 GHz radio morphology is shown as circumnuclear contours in Fig. 3. Radio synchrotron emission follows a simple power law, where the flux, f, is set by the frequency, ν, and a spectral index, α, as f ∝ ν^−α. With the values for the radio flux and central frequencies given in Sect. 2.3, we found a result of α = 0.21 between 1.52 GHz and 10.0 GHz, a typical value for the flat-spectrum radio core of a radio AGN. Although we obtained fairly deep data, with RMS = 90 μJy at 1.52 GHz, and RMS ∼ 5μJy at 7.5 and 10 GHz, we did not detect diffuse emission around the nucleus in any band.

Fig. 3. H2 1–0 S(3) emission-line morphology (top), velocities (centre), and FWHM line widths (bottom) with the JVLA 8.4 GHz radio continuum shown as thick contours near the centre. Thin contours show the emission-line morphology in the top left panel to ease orientation. Left to right: single-component fit, systemic line component and blue-shifted line component. The emission-line morphology and single-component and systemic velocity maps show an hourglass gas morphology with a strongly blue- and red-shifted bipolar velocity field, consistent with a pair of bubbles that are being inflated by the radio source at the centre. The yellow ellipse and blue circle in this and subsequent figures indicate the beam size of the JVLA and the PSF of the JWST dataset.

To investigate whether we may have missed a diffuse component, we compare the flux measured at 1.52 GHz, S_1.52 = 2.10 ± 0.09 mJy, with that of FIRST (Becker et al. 1995; White et al. 1997), S_1.4 = 3.94 ± 0.14 mJy, finding that 47% of the flux measured with FIRST at 5.4″ beam size are missing in our data obtained with a 1.2″ beam. Assuming that such emission would be of uniform surface brightness and beam-filling in the FIRST data, we find that we would need to probe down to a surface brightness of 95 μJy in the L-band, which corresponds to our RMS, so that it is well possible that faint, diffuse emission is present but not detected. In the higher frequency bands, due to the much smaller beam size, and likely steep spectral index, such an extended component would also have been missed. The 5.4″ beam of FIRST is slightly larger than the molecular disk in 3C 326 N.

A similar situation was encountered by Zovaro et al. (2019), who identified ro-vibrational line emission from molecular hydrogen in the z = 0.06 radio galaxy 4C31.04 using GEMINI/NIFS, which extended over much larger radii than the radio source. Their detailed comparison with the relativistic hydrodynamic simulations of Mukherjee et al. (2018) showed a significant component of diffuse radio emission which extends well beyond the brightest regions, but would require a dynamic range much greater than the factor 100 or so possible with current radio facilities to be probed. The same arguments are applicable to our present case.

4.2. Emission line and CO band head detections

We show the source-integrated spectrum obtained with NIRSpec in Fig. 4 and channel 1 of MIRI MRS in Fig. 5. As already pointed out by Nesvadba et al. (2011), the spectrum is strongly dominated by the ro-vibrational emission lines of molecular hydrogen. In NIRSpec we cover the 1–0 S and the 1–0 Q series from S(0) to S(7) and Q(1) to Q(9). In NIRSpec, we also identify the 1–0 O(2) and O(3) lines, and potentially the 2–1 S(3) and Q(3) lines. With MIRI we detect the rotational lines, H₂ 0–0 S(3), S(5), and S(6), previously observed by Ogle et al. (2007), as well as [ArII]λ6.985 and [NeIII]λ15.56. Ogle et al. (2007) detected additional forbidden lines of [FeII]λ5.34, [SIV]λ10.51, [NeII]λ12.81, [SIII]λ18.71, and [OIV]λ25.89, that are not covered by our data, but to our knowledge, these are the first detections of [ArII]λ6.985 and [NeIII]λ15.56. Other lines of the 0–0 S series are not covered by our setup.

Fig. 4. Source-integrated spectrum of 3C 326 N from the NIRSpec observation from 1.5 to 2.9 μm showing most of the lines. The flux was extracted with an elliptical aperture centred on the peak continuum flux with a major axis along Dec of 1.1″ and a minor axis of 0.9″. Empty pixels were not included and the pixels within the aperture were sigma clipped with a tolerance of 5σ. The spectral gap at 2.2 μm is caused by the NIRSpec IFU detector gap.

Fig. 5. Integrated source spectra from the MIRI MRS Channel 1 observation. The flux was extracted with an elliptical aperture centred on the peak continuum flux with a major axis along Dec of 1.1″ and a minor axis of 0.9″. Empty pixels were not included and the pixels within the aperture were sigma clipped with a tolerance of 5σ.

Moreover, with NIRSpec we identify the warm ionised gas line Paα as well as the iron line [FeII]λ1.644. [FeII]λ1.599 is tentatively detected but falls within the strong CO 5–2 band head absorption feature. Brγ is covered, but not detected, probably because the interstellar line emission is filling in the stellar Brγ absorption line. At R ∼ 2700, the spectral resolving power of NIRSpec does not allow us to robustly isolate these two components.

In addition to these emission lines originating from interstellar gas, we also detect the Δv = 2 and Δv = 3 stellar absorption-line series of ¹²C¹⁶O, namely, eight lines between v = 2 − 0 and v = 9 − 7 at 2.29 to 2.51 μm, and six absorption lines between v = 3 − 0 to v = 8 − 5 at 1.56–1.66 μm. We subsequently use the absorption lines between 2.3 μm and 2.4 μm to determine the systemic redshift and search for evidence of rotation of the stellar component. We do not use the redward lines due to contamination from the H₂ 1–0 Q series of emission lines and our templates do not cover the bluer absorption lines.

While a thorough analysis of all detected emission lines is beyond the scope of this paper, we provide the source integrated fluxes from our spatially resolved one- or two-Gaussian component model fitting described Sect. 3 in Table 4. In Table 3, we give the 5σ upper limits to Brγ and the 6.2 μm PAH feature, which were covered by our data but not detected. As a courtesy to the reader, we also performed single Gaussian fits to the source integrated spectrum and provide the source-integrated redshift and FWHM line widths of all emission lines except H₂ 1–0 Q(2) and H₂ 1–0 Q(8) in Table 5. The lines are more complex than a single Gaussian and therefore the given values are only to give the reader an estimate of the integrated line properties. For fluxes, we recommend only those from Table 4 be used. The two omitted Q series lines become blurred with other spectral features in the source integrated spectrum. However, they are distinct in the spatially resolved spectrum.

4.3. Continuum morphology and stellar absorption bands

3C 326 N has a well extended continuum morphology clearly observable in the near-IR. We plot the continuum morphology at 1.8 μm in Fig. 6. The extended emission is elliptical with a major axis position angle of ∼155°.

Fig. 6. 1.8 μm continuum with contours overplotted.

Since NIRSpec covers the Δv = 2 and Δv = 3 series of CO band heads, we can use these absorption features to derive the systemic redshift and stellar velocity dispersion. Stellar absorption lines are the most robust tracers of the systemic redshift, because the stellar component dominates the mass budget within the central few kpc of early-type galaxies (3 × 10¹¹ M_⊙ in 3C 326 N compared to 2 × 10⁹ M_⊙ in molecular gas Nesvadba et al. 2010), and because it is not affected by hydrodynamic processes like AGN or starburst-driven winds, which may significantly impact the gas kinematics. This is a particular worry in case of radio galaxies like 3C 326 N.

When we model the CO band heads, we find no evidence of a velocity gradient across the galaxy (see Fig. 2). This suggests that the galaxy stellar population is not rotationally dominated. Thus, 3C 326 N resembles a classical “slow rotator” as is typical for a massive galaxy of a few 10¹¹ M_⊙ of stellar mass (e.g. Emsellem et al. 2011). We also find no significant change in the FWHM of the model stellar velocity distribution (SVD) with radius. Therefore, we can find the systemic redshift by merging all posterior probability distributions of each spaxel fit. To reduce uncertainty, we merge the parameter probability distributions for all spaxels with an average wavelength bin S/N > 10 in the modelled data to find the systemic redshift and its uncertainty.

We find a systemic redshift of z_syst = 0.08979 ± 0.0003, consistent with the value of z = 0.09 from Guillard et al. (2015) measured from [CII]. Ogle et al. (2007) previously measured z = 0.089 ± 0.001 with Spitzer from the pure-rotational lines of warm H₂. Nesvadba et al. (2011) found z = 0.0903 ± 0.0002 from ro-vibrational lines and Paα measured with VLT/SINFONI. We use our value as the systemic redshift for all other velocity plots. From our CO modelling, we also find that the FWHM of the CO band head velocity dispersion is 5.1 ± 0.8 nm, equivalent to a velocity of 660 ± 100 km s⁻¹.

4.4. Emission-line morphology and kinematics

4.4.1. Ro-vibrational line emission of H₂

We start our discussion of the emission-line morphology and kinematics with the properties of the ro-vibrational H₂ 1–0 lines observed with NIRSpec, using the S(3) line as reference for this set of lines. S(3) is amongst the brightest lines we detected, not noticeably contaminated by other nearby lines, and in the blue part of the spectral range covered by our data, thus maximising spatial resolution. H₂ 1–0 S(3) can be contaminated by [SiVI]λ1.962; however, we note that we find similar modelling results in all high bright H₂ lines, including in MIRI, and that we find no obvious contamination in this line. Fig. 12 shows comparisons between different H₂ lines and shows no significant variations in line profile for H₂ 1–0 S(3). Other H₂ 1–0 lines also show the same spatial distribution and kinematics within the measurement uncertainties.

We find that the molecular hydrogen lines show a clear north-south morphology shown in Fig. 3 roughly centred on the galaxy nucleus as indicated by the peak in surface brightness of the stellar continuum emission. Bright line emission is probed over 2.0″ in north-south direction, and over 0.5″ and 1.2″, in east-west direction in the southern and northern part of the source, respectively. north of the nucleus we find a cavity-like structure delineated by several bright, unresolved knots of line emission at the far (northern) end. This cavity and the structure of the components themselves are not mirrored in the continuum supporting the idea that it is not caused by stellar emission. South of the nucleus, the emission-line morphology is more regular and featureless, and extends to about 0.8″ south from the nucleus.

We fit the Gaussian line profiles to the line emission in all spaxels as detailed in Sect. 3. Examining the line profiles (for example, those shown in Fig. 1), we find that fits with two Gaussians are preferred in almost all regions of the galaxy over single-Gaussian fits. The exception to this is the lowest surface-brightness regions where the signal-to-noise is not sufficient to fit two lines components. Some areas close to the nucleus can also be explained by a single component where the systemic redder component becomes dominant. For simplicity, we focus most of our discussion on the two-component modelling. In this discussion, we distinguish the two components as a redder (higher z) systemic (or “primary”) component and a bluer (lower z, or “secondary”), often broader, component. The line profiles, quality of the spectra and Gaussian fits can be seen in Fig. 1.

Generally speaking, the gas velocities in both Gaussian components show a monotonous gradient from blue-shifted in the south to more red-shifted in the north (Fig. 3). The systemic southern component shows red-shifted velocities along its western edge. As can be seen from the left panel of Fig. 3, the same trend can be seen in the single-component fit, although the gradient is smaller and less pronounced, because of the lower quality of the fit stemming from fitting a double-peaked profile with a single Gaussian.

The secondary component, which is always blue-shifted relative to the first, shows the same trend between blue-shifted gas southwards and red-shifted gas northwards; however, it is blue-shifted by a value of 100–200 km s⁻¹ relative to the first component. Typical line FWHM values are ∼400 − 1000 km s⁻¹ for the bluer component and ∼200 − 400 km s⁻¹ for the redder. Other H₂ 1–0 lines have similar morphologies and kinematics within the measurement uncertainties.

4.4.2. Rotational line emission of H₂

The pure rotational lines of H₂ were observed with MIRI at about 3 − 4× less spatial resolution than the ro-vibrational lines with NIRSpec, namely, with a PSF of 0.3–0.42″. We will also discuss these lines on the example of the H₂ 0 − 0 S(3) line, which is at a rest-frame wavelength of λ = 9.66 μm. These lines are interesting to study in parallel to the H₂ 1–0 lines, because they probe the dominant molecular gas component in this galaxy (Sect. 4.5).

The integrated emission-line morphology is shown in Fig. 7. H₂ 0 − 0 S(3) line emission is also clearly elongated along north-south direction with a projected size of 2.0″. In the northern and southern part of the source, the east-west extent is 1.0″, and 0.5″, respectively. Most of the line emission is originating from the blob south of the nucleus and a bright emission-line region north of it. There is a region of much fainter line emission between the nucleus and the northern emission-line region, which corresponds to the cavity already seen in H₂ 1–0 S(3). Overall, the morphology and gas kinematics closely follow those already seen in the H₂ 1–0 S(3) line, although the blending due to the lower spatial resolution of the data tend to lower the velocity offsets, enhance the line widths, and generally, smooth out small-scale features. We will do a more quantitative comparison between the two lines, including the line kinematics, in Sect. 6, that takes into account the effect of the different beam sizes. We do not observe significant differences between the morphology of the H₂ 0–0 S(3) line described here, and the 0–0 S(5) and S(6) lines, which we also covered.

Fig. 7. H2 0–0 S(3) emission-line morphology (top), velocities (centre), and FWHM line widths (bottom) with the JVLA 8.4 GHz radio continuum shown as thick contours near the centre. Thin contours show the emission-line morphology in the top left panel to ease orientation. Left to right: single-component fit, systemic line component, and blue-shifted line component. The emission-line morphology and single-component and systemic velocity maps show an hour-glass gas morphology with a strongly blue and red-shifted bipolar velocity field, consistent with a pair of bubbles that are being inflated by the radio source at the centre.

4.4.3. Warm ionised gas lines

In addition, NIRSpec and MIRI have also observed several warm ionised gas lines. Here, we discuss the lines of Paα, [FeII]λ1.644, and [NeIII]λ15.56. While it has already been shown conclusively by Ogle et al. (2007), Nesvadba et al. (2010), and Nesvadba et al. (2011) that the emission lines of warm H₂ are mainly heated by shocks driven by the radio source into the gas, the ionisation mechanism of the other lines is not clear. Paα can be either heated by shocks, AGN photoionisation, or UV photons from young stellar populations (for a review, see e.g. Kewley et al. 2019), [FeII]λ1.644 can be either heated by shocks or AGN photoionisation (e.g. Rodríguez-Ardila et al. 2004)], and [NeIII]λ15.56, which has a high ionisation potential of E_ion = 40.96 eV (Levesque & Richardson 2014), either by AGN photoionisation, or, in extreme cases, star formation (Ho & Keto 2007). Spoon & Holt (2009) and Spoon et al. (2009) found that [NeIII]λ15.56 line emission likely traces circumnuclear gas in nearby ULIRGs that is strongly perturbed through interactions with the AGN.

We can use the Paα and [FeII]λ1.644 morphology to conclude that both lines are most likely heated by shocks. Their morphology follows very closely that of the H₂ 1–0 and H₂ 0–0 lines. A contribution from AGN photoionisation near the nucleus, where Paα and [FeII]λ1.644 are brighter than the H₂ 1–0 S(3) line, is possible. They follow the north and south components seen in H₂ 1–0 S(3), but the north component is comparatively fainter relative to the south and the north component is brighter on the western side. This can be seen in the single component fits or the systemic component of the two-Gaussian fit in Figs. 3, 7, 8, and 10; we discuss this further in Sect. 5.

Fig. 8. Paα emission-line morphology (top), velocities (centre) and FWHM line widths (bottom) with the JVLA 8.4 GHz radio continuum shown as thick contours near the centre. Thin contours show the emission-line morphology in the top left panel to ease orientation. Left to right: single-component fit, systemic line component, and blue-shifted line component. The emission-line morphology and single-component and systemic velocity maps show an hour-glass gas morphology with a strongly blue- and red-shifted bipolar velocity field, consistent with a pair of bubbles that are being inflated by the radio source at the centre.

The [NeIII]λ15.56 emission, as shown in Fig. 9, does extend north-south but it does not show the same northern bright spot, and its morphology is more symmetric around the nucleus than, for instance, H₂ 0–0 S(3), which shows a depression between the nucleus and the northern clumps co-spatial with the cavity already identified in the H₂ 1–0 S(3) line.

Fig. 9. Single-component [NeIII]λ15.56 emission-line morphology (left), velocities (centre), and FWHM line widths (bottom) with the JVLA 8.4 GHz radio continuum shown as thick contours. Due to strong fringing in MIRI channel 3, the blue-shifted line, while present, is noisy and difficult show clearly. Therefore, we do not plot the two-component fit.

4.5. Mass budgets

3C 326 N has already several detailed mass estimates of different phases of warm and cold gas. Ogle et al. (2007) and Nesvadba et al. (2010) estimated a mass of warm H₂ from the pure rotational line emission finding M_H2, rot = 1.1 × 10⁹ M_⊙ and M_H2, rot = 1.3 − 2.7 × 10⁹ M_⊙, respectively. We will in the following use the latter value, which is appropriate for gas at temperatures T ≥ 150 K, and has taken explicitly into account that the gas is heated by shocks. Nesvadba et al. (2010) also found M_ion = 2 × 10⁷ M_⊙ for warm ionised gas probed by Hα. Guillard et al. (2015) estimated M_CNM = 9.5 × 10⁸ M_⊙, for gas with T = 100 K and n_H = 1000 cm⁻³. Lanz et al. (2016) estimated a dust mass of $M_{\mathrm{dust}}=6.{05}_{-0.218}^{+0.376}\times {10}^7$ M_⊙. To complete these different mass estimates, we now also provide an estimate of the warmest H₂ observed through the ro-vibrational lines.

We use the approach of Scoville et al. (1982) and Mazzalay et al. (2013) to estimate the mass of very warm H₂ from the H₂1–0 S(1) line. They set M_H2,ro-vib = 5.0875 × 10¹³D²F_{1 − 0S(1)} 10^0.4A_2.2, where D is the luminosity distance in Mpc, F_{1 − 0S(1)} the flux of the H₂ 1–0 S(1) line in erg s⁻¹ cm s⁻², and A_2.2 the extinction at 2.2 μm. Their estimate is appropriate for gas with T = 2000 K, a transition probability of A_S(1) = 3.47 × 10⁻⁷ s⁻¹, and a population fraction in the upper level of 1.22%.

We have no estimate of A_2.2, as Paα is the only line of warm ionised gas that is detected in our datasets; however, Nesvadba et al. (2010) obtained a global estimate of A_V = 0.17 mag from stellar population modelling in the optical, assuming a Calzetti et al. (2000) extinction law. This suggests that extinction in the near- to mid-infrared can be safely ignored, at least for source-integrated estimates like the one we are doing here. We will therefore adopt A_2.2 = 0 mag.

With these assumptions, and for an integrated H₂ 1–0 S(1) line flux of 1.5 × 10⁻¹⁸ W m⁻², we find M_H2,ro-vib = 1.4 × 10⁴ M_⊙. This is thus a very small mass component compared to the other estimates, in spite of the bright line luminosity. The hot gas is expected to only be a small fraction of the total gas mass because it needs to be heated to relatively high temperatures, for example by an AGN or by shocks, and the gas cools rapidly.

It is also interesting to compare the warm and cold molecular gas mass in 3C 326 N with that of other galaxies, to highlight the unusual environment that this galaxy provides. Mazzalay et al. (2013) additionally propose an empirical calibration of warm to cold gas probed by CO(1–0), which was derived empirically by comparing both masses in about 50 low-redshift galaxies spanning a wide range of types. Based on their relationship, we would expect to find a cold molecular gas mass of 9.87 × 10⁹ M_⊙. In Sect. 2 we have estimated a 3σ upper limit on the cold molecular gas mass in 3C 326 N of M_H2, CO = 8 × 10⁸ M_⊙, assuming that cold molecular gas extends over the same area as the warm molecular gas. This is more than an order of magnitude less than what would be expected in normally star-forming galaxies.

Likewise, from the dust mass estimated by Lanz et al. (2016) from Herschel SPIRE observations, $M_{\mathrm{dust}}=6.{05}_{-0.218}^{+0.376}\times {10}^7$ M_⊙, we would expect a cold molecular gas mass of about 6 × 10⁹ M_⊙, assuming a typical gas-to-dust mass ratio of G/D = 100, as would be typical for massive early-type galaxies with solar or super-solar metallicity (Rémy-Ruyer et al. 2014). This is again about a factor 10 greater than what is observed from CO line emission for cold molecular gas. Both estimates fall however within factors of a few from the warm H₂ mass observed in 3C 326 N.

5. Discussion

5.1. Nature of the emission-line regions in 3C 326 N

In the previous section, we describe the properties of the molecular disk in 3C 326 N, which is roughly centred on the nucleus and extends to a radius of about 1.5 kpc. This disk is strongly kinematically disturbed with dual-component line profiles throughout the disk, and large line widths, typically 200–400 km s⁻¹ in the systemic, and up to 1000 km s⁻¹ in the secondary, bluer component. Both components show a velocity gradient that increases from south to north, and which resembles disk rotation when seen at low spatial resolution (Nesvadba et al. 2011). The two components are typically offset from each other by about 80 − 100 km s⁻¹.

The excellent spatial resolution of NIRSpec reveals a large cavity in the northern part of the disk, that is 0.6″ × 0.8″ across along its minor and major axis, respectively, corresponding to 1.0 kpc ×1.4 kpc at the redshift of 3C 326 N (Fig. 3). This cavity is delineated by three bright clumps in the far north. Inside the cavity the gas is much fainter (but still detected) and reaches the highest velocities relative to the nucleus (< 380 km s⁻¹), which are best seen in the central panel of Fig. 3.

Such a cavity may either be produced by an expanding bubble driven by AGN or star-formation feedback (e.g. Capetti et al. 1999; Mayya et al. 2023); alternatively, it may indicate the presence of a massive dust lane across the galaxy. Star-formation can be ruled out given the low star formation rate estimated with Herschel, $\mathrm{SFR}=0.{087}_{-0.046}^{+0.106}$ M_⊙ (Lanz et al. 2016). We do not identify any bright, warm dust or PAH emission at this location, and in the optical, the global extinction is very low, A_V = 0.17 mag (Nesvadba et al. 2010), so that the presence of a massive dust lane able to obscure near and mid-infrared emission in the H₂ 1–0 S(3) and 0–0 S(3) lines, seems unlikely. We thus conclude that we are seeing an inflating, AGN-driven bubble expanding through the molecular gas, reaching an expansion velocity of up to 380 km s⁻¹.

The single-line fits plotted in the left panel of Fig. 3 do not show the true projected gas velocities, because the complex line profiles are not well represented with a single Gaussian, but they probe deeper than the dual-Gaussian fit. They show a clear hour-glass shaped velocity distribution, which in the north is associated with the cavity, and surrounded by gas at more moderate velocities. The same feature is identified in the systemic (redder) component in the dual-Gaussian fits.

In the south, we did not identify an analogous blue-shifted bubble, but the global gas morphology and velocity field of the systemic component are clearly hourglass shaped, implying a symmetric underlying structure. The onset of this structure can be seen in the right panel of Fig. 3, which shows the blue-shifted gas component of H₂ 1–0 S(3), with velocities that are offset by up to about −300 km s⁻¹ from the systemic redshift. The structure is better seen in H₂ 0–0 S(3) (Fig. 7), and also in warm ionised gas, in particular Paα (Fig. 8).

The gas velocities are more than previously estimated from the seeing-limited VLT/SINFONI data (Nesvadba et al. 2011) at 1″ spatial resolution; they are also significantly more than the rotation velocity expected for a massive elliptical galaxy. The turnover velocities of rotating gas disks in elliptical galaxies with 3 × 10¹¹ M_⊙ in stellar mass, similar to 3C 326 N, are ∼150 km s⁻¹, with turnover radii around 1.5 kpc, comparable to the size of the molecular disk in 3C 326 N (Yoon et al. 2021).

The most likely explanation of these features is that we are seeing a bubble expanding within a massive, rotating, molecular disk, where the observed velocity gradient is most likely a combination of disk rotation and outflowing motion. The excellent 0.1″ resolution provided by NIRSpec is critical to identify the complex kinematic features of such a configuration, and to measure gas velocities robustly. Nonetheless, the overall complexity of the gas kinematics and morphology make it difficult to clearly separate a kinematically perturbed from the quiescently rotating part of the disk, and to derive, for instance, a robust estimate of the gas energetics in a simple, unique, and analytical way. In a companion paper, we will use hydrodynamic simulations to characterise the interactions between AGN and gas disk in 3C326 N in a more robust way (Shende et al., in prep.) Given that star formation rates and the bolometric luminosity of the AGN are very low, and are not sufficient to either power the observed emission-line fluxes or drive outflows or major gas motions (Nesvadba et al. 2010, 2011), the only possible culprit that could inflate this bubble is the radio jet.

The radio core is compact in all frequency bands in which we observed 3C 326 N with the VLA (Sect. 2.3), corresponding to beam sizes between 0.2″ and 1.2″, respectively, at 9 GHz and 1.4 GHz, respectively. However, we also observe 47% of missing flux at 1.4 GHz compared to the flux measured at 5.4″ beam size with FIRST (Becker et al. 1995, Sect. 4.1). This suggests that significant amounts of radio emission must be present at scales between about 1 and 4.6 kpc around the nucleus. This is likely the radio plasma that is permeating the gas disk in 3C 326 N, stirring up the gas, creating turbulence, and inflating the cavity seen in the northern hemisphere. A similar situation has previously been found by Zovaro et al. (2019) in the radio galaxy 4C 31.04.

Qualitatively similar, bipolar gas motions have already been predicted by Meenakshi et al. (2022) from hydrodynamic simulations of the passage of radio jets through a surrounding gas disk (their Fig. 5). The jet energy of their simulation, P_j = 10⁴⁴ erg s⁻¹, is somewhat less than that estimated by Nesvadba et al. (2010) for 3C 326 N, P_jet, 326 = 2 − 13 × 10⁴⁴ erg s⁻¹, obtained with different methods ranging from X-ray cavity inflation (Bîrzan et al. 2008) to the 5 GHz core luminosity (Merloni & Heinz 2007). Nonetheless, the simulated bubbles of Meenakshi et al. (2022) reach gas velocities of up to ∼300 km s⁻¹, slightly lower than the 380 km s⁻¹ that we observed (albeit not correcting for inclination), illustrating that the jet in 3C 326 N may well cause the outstanding gas properties that we observe.

We also note that warm ionised gas in the hourglass component has a general redshift of about 100 − 200 km s⁻¹ relative to H₂. This effect is only seen in the systemic gas component, which shows this hourglass structure, not in the blue-shifted component. We are not aware of any systematic effect that could cause this behaviour in only one line component, and which is also seen in other warm ionised gas lines (Fig. 12). We might be seeing an effect of stratification, where the warm ionised gas is accelerated to somewhat different velocities, or where different parts of the same structure are particularly bright in the warm ionised and molecular gas (namely, different components of the same structure.).

5.2. An extended, turbulent, multi-phase interstellar medium throughout the galaxy

The spatially resolved line ratios in 3C 326 N offer interesting information about the gas conditions and heating mechanism of the gas in this galaxy. While a detailed analysis of the physical conditions of the warm molecular gas is beyond the scope of this paper, we can draw several interesting conclusions about these mechanisms from comparing the line ratios in the different components of the galaxy.

Nesvadba et al. (2010) and Nesvadba et al. (2011) already showed that the source-integrated rotational and ro-vibrational line emission of this galaxy is heated by shocks; however, their data did not have the spatial resolution necessary to resolve the different components of this galaxy, and to search for pockets of star formation that could be embedded in a globally shock heated interstellar medium.

As first argued by Puxley et al. (1990), ro-vibrational lines in the NIR are an excellent shock tracer, when comparing their line fluxes with those of warm ionised gas lines like Brγ or Paα. Puxley et al. (1990) used 44 HII regions observed in 30 galaxies to infer that line ratios between H₂ S(1) and Brγ are typically between 0.1 and 1.5 in star-forming regions, with an average of 0.5 ± 0.3. We note that this empirical ratio does not take into account potential effects from shocks present in these star-forming regions. Outliers in their sample with much brighter H₂ emission include in particular NGC 6240 and NGC 3079, known to have large amounts of shocked gas (e.g. Keel 1990; Veilleux et al. 1994; Meijerink et al. 2013; Medling et al. 2021).

Instead of Brγ and S(1), in the following we use Paα and S(3) for our comparison. Paα is significantly brighter than Brγ, but falls in between the atmospheric windows for nearby galaxies observed from the ground. For a galaxy at z ∼ 0.1 like 3C 326 N and observations from space like with JWST, this is however not a concern. Moreover, Paα is less affected by underlying absorption lines from stellar photospheres, which can artificially lower the measured Brγ flux in medium-resolution spectra like those provided by the NIRSpec IFU.

The decrement between Brγ and Paα for gas with an electron temperature, T_e = 10⁴ K, is F_Paα/F_Brγ = 11.2 (e.g. Dopita et al. 2003), and the line ratio between H₂ 1–0 S(3) and S(1) in 3C 326 N is on average F_S(3)/F_(S1) = 1.2 (Table 4). This implies that the Puxley et al. (1990) criterion of F_Brγ/F_{H21 − 0S(3)} ∼ 1 for line emission from shocked gas translates into a criterion of F_{H21 − 0S(3)} > 0.1 F_Paα.

Line emission throughout 3C 326 N is well above this threshold. When we take the ratio of total source Paα to H₂ 1–0 S(3) we find the same low value 0.5 as in Nesvadba et al. (2011). If we take the ratios within the apertures defined in Fig. 11, we find a ratio of 0.5 for aperture S, of 0.6 in aperture NW, 0.45 in N, and 0.4 in NE.

This implies that shocks dominate the gas heating of molecular gas throughout the galaxy not only in a global sense, but also for individual features seen in the line maps of Fig. 3. This includes in particular the compact clumps in the north of the galaxy, which from their morphology are reminiscent of young star-forming clusters, as are commonly found in other galaxies. The line ratios however, suggest that they are not heated by UV photons from young stars, but from shocks like the ambient gas, otherwise Paα should be an order of magnitude brighter compared to the S(3) line than what is observed. Thus, they are more likely to be gas that has fragmented along the rim of the jet-inflated bubble.

Shocks are also the dominant heating mechanism of the pure-rotational line emission of H₂ in 3C 326 N, as already shown by Nesvadba et al. (2010). They constructed a diagnostic diagram based on PAH, H₂ 0–0 S(1), CO(1–0), and 24 μm dust emission, and compared with expectations from theoretical models of star-forming, photon-dominated regions (PDRs), to find that the source-integrated S(1) flux of 3C 326 N exceeded the gas heating possible with UV photons from young stars by more than an order of magnitude. They also disfavoured cosmic ray heating, arguing that the required cosmic ray flux densities would be so high that the H₂ molecule would be destroyed, leaving shocks as only plausible gas heating source.

Nesvadba et al. (2010) also compared the source-integrated flux measurements of H₂ 0–0 S(0) to S(7), and fitted the shock models of Flower & Pineau Des Forêts (2010), finding that at least three different shock velocities between v_s = 4 and 40 km s⁻¹ were needed for C-type shocks to explain the excitation diagram, for two assumptions on gas density, respectively (n_H2 = 10³ and 10⁴ cm⁻³, respectively). From their Fig. 3, we can see that the H₂ 0–0 S(3) line is most sensitive to gas heated by intermediate-velocity shocks, whereas the S(5) line is dominated by the highest shock velocities, and representative of the least important mass component. This is one of the reasons why we focus the discussion of the warm molecular gas in 3C 326 N onto the S(3) component.

More recently, Villa-Vélez et al. (2024) showed that the H₂ lines from our work agree extremely well with predicted fluxes they derive by performing an in-depth modelling of previously observed ionised and H₂ emission lines in 3C 326 N with a photo-dissociation and shock heating model. Their results confirm that the global H₂ emission lines of this work are dominantly shock excited.

We also note that, on top of the similar gas morphologies in different lines of molecular and warm ionised gas highlighted in Sect. 4, the line profiles and velocities of the H₂ 1–0 S(3), H₂ 1–0 S(5), H₂ 0–0 S(3), and H₂ 0–0 S(5) all agree very well. [FeII]λ1.644 and Paα also agree with each other and relatively well with the molecular line profiles. This can be seen by comparing the line maps shown in Figs. 3, 7, 8, and 10, and also by comparing the line profiles extracted from several apertures within 3C 326 N (shown in Fig. 12).

Fig. 10. [FeII]λ1.644 emission-line morphology (top), velocities (centre), and FWHM line widths (bottom) with the JVLA 8.4 GHz radio continuum shown as thick contours near the centre. Thin contours show the emission-line morphology in the top left panel to ease orientation. Left to right: single-component fit, systemic line component and blue-shifted line component. The emission-line morphology and single-component and systemic velocity maps show an hour-glass gas morphology with a strongly blue- and red-shifted bipolar velocity field, consistent with a pair of bubbles that are being inflated by the radio source at the centre.

Figure 12 again shows the offset already discussed in Sect. 5.1 between the velocities of the ionised and molecular emission. This offset is most apparent in the systemic line component which is red-shifted with respect to the molecular lines by 100–200 km s⁻¹ in the southern component, the galaxy centre, and the centre of the northern component (aperture N in Fig. 11). The systemic component is then in agreement in velocity with the molecular lines in the western wing of the northern component. However, in the eastern wing of the northern component (aperture NW) both ionised emission lines are double peaked with both peaks being comparable in amplitude. The second peak is blue-shifted with respect to the molecular lines by 200 km s⁻¹.

Fig. 11. Comparison of rotational and ro-vibrational line emission for several apertures. Top: integrated flux for the H20–0 S(5) line and the H21–0 S(5) line as titled. The H21–0 S(5) flux image has been convolved to the PSF at the H20–0 S(5) line. The original PSF for H21–0 S(5) and the post-convolution PSF are the white and black filled circles in the H21–0 S(5) plot, respectively. Bottom: ratio of 1–0 to 0–0 line fluxes for the apertures shown in the top image.

Fig. 12. Extracted line profiles for a selection of observed lines. The zero velocity is set by the CO band head absorption lines. The lines are extracted using the apertures in Fig. 11 as indicated by the coloured marker in the corner of each subplot.

The one line with a strongly different morphology is [NeIII]λ15.56, which does not follow the same morphology seen in the other lines, but is fairly symmetric about the nucleus. [NeIII]λ15.56 can be photoionised by AGN, star formation in extreme cases, or shocks. Nesvadba et al. (2010) estimated a line ratio of [NeII]λ12.81/[NeIII]λ15.56 ∼ 2.2 for shock conditions favoured by the molecular gas lines. We can use the [NeII] flux measured by Ogle et al. (2007) of F([NeII]λ12.82) = 4 × 10⁻¹⁹ W m⁻², and the source-integrated [NeIII]λ15.56 flux in Table 4, to measure a source-integrated line flux of [NeII]λ12.81/[NeIII]λ15.56 = 4.3, which is about twice as large as the value expected for shocked gas. This would be consistent with having a second line component present, presumably heated by AGN photoionisation.

In conclusion, the ionised and molecular lines follow globally the same approximate bulk velocity, although local differences do exist. Similar gas kinematics in different gas phases are expected in multi-phase media, where the different gas phases do not probe different components of a galaxy, but are co-spatial and undergoing a constant mass and energy exchange (Guillard et al. 2009). This is characteristic for gas in non-equilibrium, rapidly changing conditions, for instance, in rapidly cooling post-shock gas following the passage of a shock front, coming from an expanding radio source, for example (Mukherjee et al. 2018).

5.3. 3C 326 N in the context of turbulence-driven star formation

5.3.1. Kennicutt–Schmidt diagram

One of the outstanding properties of radio galaxies with large amounts of warm molecular hydrogen are their low star formation rates compared to their molecular gas content. Schmidt (1968) and Kennicutt (1989) showed that gas-mass surface densities are a major determinant of the star formation rate density in galaxies. A subset of MOHEGs show a significant offset from the typical relationship between gas mass and star formation rate surface densities in normal star-forming galaxies by factors 10–100. This offset is seen even in spiral galaxies with normal stellar mass surface densities and moderately bright radio jets (Nesvadba et al. 2021).

The offset of radio galaxies in the Kennicutt–Schmidt diagram was first noticed by Nesvadba et al. (2010), who compared the molecular gas mass surface densities of these galaxies with star formation traced by PAH emission. PAHs are an imperfect tracer of star formation, as PAHs are relatively easily destroyed by shocks (Micelotta et al. 2010; Egorov et al. 2023). Subsequently, this observation was confirmed with Herschel/SPIRE photometry (Lanz et al. 2016). The star formation rate of 3C 326 N that these measurements found is $\mathrm{SFR}=0.{087}_{-0.046}^{+0.106}$ M_⊙ yr⁻¹.

These new ALMA observations (Sect. 2.2) have revealed that the previous CO(1–0) detection of Nesvadba et al. (2010), based on shallower data observed with the IRAM Plateau de Bure Interferometer in 2009, could not be confirmed. These data had much smaller bandwidth, and a large, 5″ beam, making it impossible to separate the central continuum emission from putative more extended line emission. The new ALMA data place an upper limit of 8 × 10⁸ M_⊙ on the cold molecular gas mass derived from the CO(1–0) line flux for a Milky-Way like H₂-to-CO conversion factor of α_CO = 4.3M_⊙/ (K km s⁻¹ pc²).

In Fig. 13, we show the Kennicutt-Schmidt diagram with the updated position of 3C 326 N, taking these new measurements into account. We estimated a gas mass surface density of cold molecular gas by using the upper limit on the molecular gas derived in Sect. 2.2, and assumed that the cold molecular gas morphology follows that of H₂ 0–0 S(3). We also assume that star formation is distributed over the same area. This is a conservative estimate, as we do not observe any signatures of star formation in our data. With these estimates and assumptions, 3C 326 N falls a factor 23 below the usual Kennicutt-Schmidt relationship for normally star-forming galaxies. Villar Martín et al. (2023) show that cold gas morphology can be smaller than the H₂; however, a smaller morphology would only increase the discrepancy.

Fig. 13. Kennicutt–Schmidt law between star formation rate and molecular gas mass surface density. The large blue and red dot show the position of 3C 326 N when using the mass surface density of cold molecular gas probed by CO(1–0) and H2 0–0 lines, respectively. For the latter estimate, the black line shows the range of surface-brightness observed with MIRI. Small empty circles show several sets of normally star-forming galaxies in the nearby Universe taken from Fischer et al. (2013), Kennicutt (2020), and Saintonge et al. (2013). Green and grey squares show the MOHEGs of Ogle et al. (2010) and Alatalo et al. (2015).

For comparison, we also show where the galaxy falls when using the observed H₂ 0–0 S(3) mass surface densities. To estimate these, we adopted the total mass estimate of warm molecular hydrogen of 2 × 10⁹ M_⊙ (Nesvadba et al. 2010) and assumed that changes in gas excitation are not large enough to strongly alter the overall mass-to-light ratio of the gas. This allowed us to estimate a range of gas-mass surface densities, Σ = 200 − 680 M_⊙ pc⁻², from the surface-brightness distribution of the line shown in Fig. 7. Showing the position of 3C 326 N for warm molecular hydrogen in this diagram bears of course a large risk of systematic errors, but allows us also to estimate where this galaxy would approximately fall if the gas heating that is required to keep this gas warm, was absent. It can be seen that the galaxy falls about two orders of magnitude below normal star-forming galaxies in this case.

5.3.2. Turbulent energy in the warm molecular gas

The large line widths observed with NIRSpec and MIRI MRS show clearly that the gas in 3C 326 N is strongly stirred up and likely made turbulent through the interactions with the radio jet (Nesvadba et al. 2010, 2011). The 0.11″ size of the PSF corresponds to a projected size scale of 190 pc, which demonstrates that the large line widths that have been previously observed are not due to velocity offsets between different gas components on kpc scales, but are present even at size scales typical of giant molecular clouds.

This is an important observation, because interstellar turbulence is now recognised as a major regulation mechanism of star formation (Krumholz & McKee 2005; Federrath et al. 2010; Padoan & Nordlund 2011; Hennebelle & Chabrier 2011; Hennebelle & Falgarone 2012). In 3C 326 N, the only energy source that is powerful enough to power velocity offsets and line widths as observed, is the radio jet. We know also from hydrodynamic simulations that radio jets can make surrounding gas turbulent (e.g. Mukherjee et al. 2018), and impact star formation, by locally compressing, but globally dispersing the gas (Mandal et al. 2021; Mandal et al., in prep.).

The virial parameter, α_vir allows us to quantify the amount of turbulent energy compared to gravitational binding energy in molecular clouds (Bertoldi & McKee 1992; Sun et al. 2018). We set α_vir = 5σ_cl²(πGR_clΣ_gas)⁻¹, following Bertoldi & McKee (1992). σ_cl is the line width dominated by turbulent motion, Σ_gas the gas-mass surface density, R_cl the size of the cloud, and G the gravitational constant. Giant molecular clouds in the Milky Way, as well as nearby spiral galaxies are found to be star-forming for α_vir ∼ 1. In gas with much higher turbulent velocity, star formation is suppressed. Nesvadba et al. (2011) showed with low-resolution SINFONI data, that the high interstellar turbulence in 3C 326 N caused by the radio jet can suppress star formation, provided that high turbulent velocities prevail down to the typical scales of giant molecular clouds and cloud complexes.

With our NIRSpec and MRS data, we have been able to test this prediction, by measuring the line broadening at the spatial resolution limit of our data, namely: 0.11″ for H₂ 1–0 S(3), corresponding to 190 pc, and 0.3″ for H₂ 0–0 S(3), corresponding to 530 pc. In particular, we identified the narrowest widths of both lines, FWHM = 170 km s⁻¹ for H₂ 1–0 S(3) and FWHM = 147 km s⁻¹ for H₂ 0–0 S(3).

Both line widths are shown in Fig. 14 as large, upward-pointing arrows to indicate that these values are the lowest of a set of measurements, where the highest values are > 1000 km s⁻¹. The position of these measurements is compared with several expected distributions. The red dotted line indicates where clouds fall with α_vir = 1, for a gas mass surface density of Σ = 250 M_⊙ pc⁻² appropriate for 3C 326 N. Molecular clouds in the Milky Way fall slightly below that line because the Galactic value is somewhat lower (Heyer et al. 2009). The yellow band indicates the range implied by Larson’s scaling relationship between radius and turbulent velocity for turbulent molecular clouds (Larson 1981), extrapolating from the measured values on 1 kpc scales downward to the scales probed by our new data, and beyond. We allowed for a variation range of about a factor 8. The blue line supposes that the H₂ line emission observed in 3C 326 N is powered by the dissipation of turbulent energy on sub-parsec scales, which implies that L_H2/M_H2 = 3/2f_H2σ_cl³/R_cl (McKee & Ostriker 2007). Also, L_H2 and M_H2 are the luminosity and mass of warm H₂, respectively, with L_H2/M_H2 = 0.06 L_⊙/M_⊙ (Nesvadba et al. 2010, 2011). We allowed for a fiducial fraction of about 50% to the total radiative gas cooling from other gas phases, parameterised by the correction factor f_H2 = 0.5.

Fig. 14. Gaussian line width, σ, of warm molecular gas in 3C 326 N as a function of radius. The minimal line widths observed in H2 1–0 S(3) and H2 0–0 S(3) at spatial resolutions corresponding to 190 pc and 530 pc, respectively, are shown as large, upward-pointing arrows. The broadest line widths observed in both lines are > 1000 km s−1, and fall well outside the plot. Black dots show previous SINFONI observations of H2 1–0 S(3) at kpc-scale resolution. The red dotted line shows the position of gas clouds with αvir = 1, for the gas-mass surface densities observed in 3C 326 N. Light blue circles show molecular clouds in the Milky Way. The yellow band indicates the range where the data should fall if the Larson (1981) scaling relations of turbulent clouds apply. The blue line indicates where the lines should fall for typical mass-to-light ratios of gas cooling through turbulent dissipation. See text for details.

We find that the new observations fall well within the range predicted by the previous VLT/SINFONI observations obtained on kpc scales, which reinforces the previous statement of Nesvadba et al. (2011) that the absence of star formation in 3C 326 N is well explained by the high turbulent gas velocities. One could argue that the high temperature rather than the mechanical energy could be the cause, however, since mechanical interactions of the radio jet with the gas are the primary heating mechanism in this galaxy, the gas heating is a concomitant effect of the same process.

The good agreement of our new and old data with the Larson scaling relationship and radiative efficiency of H₂ line emission produced through the dissipation of turbulent energy adds further credence that the interactions between radio jet and gas are indeed raising the turbulent energy in the gas of 3C 326 N. Thanks to the high spatial resolution of the new data, we can also rule out that small pockets of less turbulent gas exist in some regions of 3C 326 N. Very broad line widths are found to be ubiquitous in the molecular gas throughout the galaxy, at least down to scales of 200 pc.

6. On using the H₂ 1–0 lines as proxy to the H₂ 0–0 lines in distant galaxies

The rotational H₂ lines probe the by far dominant mass component of warm, shocked hydrogen, and are thus very interesting tracers of interstellar turbulence. This makes them important lines to study in particular in distant galaxies at intermediate to high redshifts, where massive galaxies were forming stars at high rates from interstellar gas that had significantly higher mass surface densities and turbulent velocities than today (e.g. Guillard et al. 2015). Unfortunately, the rotational lines are red-shifted out of the MIRI bandpass at z = 1.89 for H₂ 0–0 S(3), and at z = 3.04 for H₂ 0–0 S(5) for the nominal spectral coverage of the MRS. The redshift limit becomes even more constraining when taking into account the loss in count rate in Channel 4 observed during the first year of JWST operations. For example, without including Channel 4C, the limiting redshifts are z = 0.83 for H₂ 0–0 S(3), and z = 1.56 for H₂ 0–0 S(5).

The situation is much better for the ro-vibrational lines which fall into the MRS bandpass out to, for instance, a redshift of z = 12.1 and z = 7.3, respectively, for observations of 1–0 H₂ S(3) with and without Channel 4C. However, it is not clear to what degree observations of the ro-vibrational lines can be used as a proxy of the pure-rotational lines. The temperature range they probe is much higher than for the rotational lines (T ∼ 1000 K instead to T ∼ 100 K), and in some environments, where the mid-infrared spectrum is strongly dominated by rotational line emission from warm shocked molecular hydrogen, like Stephan’s Quintet (Appleton et al. 2006), the ro-vibrational lines are very faint.

Our MIRI and NIRSpec observations provide us with an excellent opportunity to study in a very clean environment not contaminated by star formation or strong AGN emission, to what degree information from ro-vibrational lines heated by shocks can be trusted to infer the basic properties of the gas probed by pure-rotational lines. For different regions of the galaxy we compare the 1–0 and 0–0 S(3) and S(5) lines to demonstrate whether the 1–0 lines are good proxies for the 0–0 lines in the case of galaxies with strong shocks. To this end, we convolved the data cubes of H₂ 1–0 S(3) and S(5) with a Gaussian PSF to produce data cubes that have the same spatial resolution as the data cubes of the mid-infrared lines.

We blurred the 1–0 lines with a Gaussian using a FWHM based on the quadratic difference between the 1–0 and 0–0 PSFs. Both 1–0 lines are covered by NIRSpec and therefore can be explained by a PSF of 0.115″ in our observing mode (D’Eugenio et al. 2023). The PSF for the 0–0 lines can be derived by setting FWHM_PSF = 0.033λ(μm)+0.106 from Law et al. (2023). When extracting line fluxes, we used the apertures in Fig. 11, while discarding pixels that have a flux of more than 5σ away from the average of the 4 × 4 pixel median filtered image where σ is the standard deviation of the difference between the model flux image and the median filtered image. The spectral resolving power is similar with both instruments, so that no further correction is required in spectral direction. This allows us to compare directly whether both sets of lines probe similar morphologies and gas kinematics.

This can be seen in the images of Fig. 11. For example, there is a bright blob in the northern part of the galaxy that is less pronounced in the 1–0 line flux than that of the 0–0 line. Overall, the ratio for both the 0–0 and 1–0 S(3) and S(5) lines differ slightly between the northern and southern components, but remain consistent within 15%. The average 0–0/1–0 ratio for all apertures in Fig. 11, not including the cavity or total source (blue and green apertures), is 5.5 ± 0.6. The southern component has a lower line ratio than the northern component for both S(3) and S(5) lines.

These variations indicate that the shock conditions are not exactly the same everywhere in the galaxy, but that the uncertainties that this implies, for instance, for determining global emission-line morphologies and to isolate individual structures such as shocks and clumps within the galaxy, are not larger than other systematic uncertainties in this kind of work, for instance, due to the stellar initial mass function or the CO-to-H₂ conversion factor. Of course, for a more detailed analysis of the gas conditions, measurements of the 0–0 lines are still required.

The gas kinematics probed by the rotational and ro-vibrational lines are also very similar throughout the galaxy. Figure 12 shows a comparison of the line profiles of H₂ 1–0 S(3) and S(5), H₂ 0–0 S(3) and S(5), of Paα, [FeII]λ1.644, and [NeIII]λ15.56. This comparison is made for the six apertures shown in Fig. 11.

The systemic components of the 1–0 and 0–0 lines are in good agreement with each other in terms of the velocity and line widths for both S(3) and (5). The blue-shifted component is more enhanced in the 1–0 lines than in the 0–0 lines in certain apertures, however, these differences are small. This suggests that turbulent line broadening and global outflow properties can be estimated from the 1–0 lines when the 0–0 lines are not available. We conclude that the shocked molecular hydrogen in galaxies dominated by turbulent gas motion shows globally similar kinematics and morphologies in the 1–0 and 0–0 lines. Differences present in blue wings and line ratios however suggest that dedicated observations of the 0–0 lines are required to derive robust estimates, for instance, for the outflow masses or gas excitation conditions.

7. Summary

We have obtained new Cycle-1 open time observations with the JWST NIRSpec and MIRI integral-field units to study the interactions of the radio jet and molecular gas in the nearby (z = 0.09) radio galaxy 3C 326 N. This galaxy is one of the best-studied examples of a subset of radio galaxies characterised by bright line emission from warm molecular hydrogen that was previously seen with the Spitzer Space Telescope (Ogle et al. 2007, 2010); meanwhile, typical star-formation tracers like PAHs or a dust continuum are either weak or undetected. Subsequently, Nesvadba et al. (2010) showed that some of these galaxies fall well below the usual star-forming sequence in the Kennicutt-Schmidt diagram, suggesting that star formation may proceed at lower rates than in galaxies without radio jets. The main results of our new study are as follows:

• We find bright line emission for 19 ro-vibrational lines of the S, Q, and O series of warm molecular hydrogen, along with three lines of pure-rotational line emission, H₂ 0–0 S(3), S(5), and S(6). We also obtained a robust systemic redshift estimate of z = 0.08979 from the CO band heads observed with NIRSpec, using a K4III giant as well matching template.

• The gas probed by these emission lines is in a 3 kpc-sized disk, already identified by Nesvadba et al. (2011) with seeing-limited data. At 10× higher spatial resolution with NIRSpec, we identify a kpc-sized cavity in the northern hemisphere of the disk, delineated by three bright knots of line emission at the far end. We interpret this feature as the signature of an expanding, jet-inflated bubble within a rotating gas disk. Similar morphologies are seen in ro-vibrational and rotational lines, as well as the warm ionised gas lines [FeII]λ1.644, Paα and [ArII]λ6.98. Only [NeIII]λ15.56 has a distinct morphology that is more centred on the nucleus.

• The line profiles are complex, with at least two line components throughout the disk seen in the ro-vibrational and rotational lines of warm H₂ as well as the warm ionised gas lines of Paα and [FeII]. We developed a new Bayesian line fitting technique based on a MCMC approach, which allows us to robustly fit both line components throughout the dataset. The cavity is associated with the red-shifted, northern side of a clear hourglass-shaped velocity distribution in the southern and northern part of the disk, with maximal velocities of 380 km s⁻¹. We interpret this velocity pattern as the signature of two expanding bubbles, given the large velocity gradient and line widths, more moderate velocities in fainter surrounding gas, and the clear association with the cavity in the north.

• We also present new JVLA radio observations at 1.5 GHz, 7.3 GHz, and 9 GHz, identifying a compact radio core. A comparison with FIRST shows 47% of missing flux at 1.5 GHz on scales between 1 and 4.6 kpc from the nucleus. This is likely emission from the radio plasma that is inflating the cavity, stirring up the gas, and making it turbulent.

• Line emission of H₂ and most warm ionised lines is heated by shocks. This has already been shown based on source-integrated data (Nesvadba et al. 2010); however, we can now rule out the presence of smaller areas dominated by heating from young stars, down to scales of 190 pc.

• Based on archival ALMA CO(1–0) observations, we found that a previous detection of cold molecular gas done with the IRAM Plateau de Bure Interferometer in 2009 was spurious. We have now placed an upper limit of 8 × 10⁸ M_⊙, which is about 3× less than the cold molecular gas mass of Nesvadba et al. (2010). This highlights the outstanding nature of this galaxy. 3C 326 N falls still below the Kennicutt-Schmidt law if the new upper limit is not very far from the actual gas mass. When taking into account the fact that maintaining the observed high gas mass requires near-constant gas heating through the radio jet (Nesvadba et al. 2010), we also argue that, in the absence of that heating, 3C 326 N would fall about a factor of 100 below the Kennicutt-Schmidt relation of star forming galaxies. This is an argument that results from using the warm instead of the cold molecular gas mass for this comparison.

• The broad line widths and bright line emission can be taken as indicative of enhanced turbulent gas motion and rapid turbulent energy dissipation. We compared the observed line widths and size of the scales with the characteristics of molecular clouds and the line emission produced by turbulent dissipation on sub-pc scales, finding that our data agree with this scenario from kpc scales and down to scales of 190 pc. Star formation theory postulates that enhanced turbulence can suppress star formation for line widths, as observed in 3C 326 N. This reinforces previous arguments that the turbulence injected by radio jets into the interstellar gas of their host galaxies may indeed inhibit star formation, as postulated by cosmological simulations of galaxy evolution.

3C 326 N presents a very clear environment in which shocks dominate the properties of the interstellar gas. We compared the morphology and line widths of the ro-vibrational and rotational lines of warm H₂ to infer whether the former can be used as a proxy to derive the properties of the latter at redshifts where rotational lines have moved out of the MIRI bandpass, for instance, at Cosmic Noon. Rotational lines probe by far the dominant mass component of warm molecular gas. Galaxies in the early Universe often show significant line broadening attributed to enhanced levels of turbulence, and often contain high column densities of molecular gas, which makes observations of warm H₂ particular interesting. Our results suggest that while local differences do exist due to a variation in excitation conditions, the ratio of the S(3) to S(5) line fluxes in both phases agree to within about 60%, while the line widths remain within the measurement uncertainties. This is comparable to other systematic uncertainties related to studies of distant galaxies, such as those due to the stellar initial mass function or the CO-to-H₂ conversion factor. This suggests that ro-vibrational line emission may serve as an interesting tracer for studying the turbulent interstellar medium of distant galaxies, for instance, at Cosmic Noon and all the way up to the Epoch of Reionisation for particularly bright galaxies and gravitational lenses.

Acknowledgments

We would like to thank the referee for the kind and insightful comments and suggestions which helped improve this work. We would further like to thank the staff at STScI, at the JVLA, and at ALMA for their fantastic assistance with preparing the data taking, and with carrying out the observations. JL, NPHN and DM acknowledge partial financial support through CEFIPRA programme 6504-2. This work was supported by the French government through the National Research Agency (ANR) with funding grant ANR AGN_MELBa (ANR-21-CE31-0011). This work is based in large parts 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 programme #GO1-2162. This work is partly based on data obtained with the Karl G. Jansky Very Large Array (JVLA) of the National Radio Astronomy Observatory (NRAO) through programme ID 18A-285. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This paper makes use of the following ALMA data: 2015.1.01120.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC 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

All Tables

All Figures

Fig. 1. H21–0 S(3) Gaussian line model results compared to data for different spaxels and resulting integrated line emission map. Middle: integrated model flux for H21–0 S(3). Edges: Gaussian model fit to the H21–0 S(3) line for different spaxels of the central image compared to data. The coloured squares in the image represent the position of the plotted spaxel of the same colour. Dashed lines are individual Gaussian components and the solid line is the summed Gaussians. Shaded regions represent the 1σ errors. These spectra area also a good illustration of the high quality of the date and the automatic fitting routine.

In the text

Fig. 2. Results of CO band head template modelling. Top-left: CO band head velocity map from our template fit. For clarity, while we model all spaxels, here we only plot the velocities for spaxels with an average signal-to-noise ratio of S/N > 10. Top right: FWHM of the velocity dispersion for the CO features, the spaxels with an average S/N < 10 are omitted. Overplotted as coloured squares: Positions of the four spaxels displayed in the bottom four panels. The spaxel positions and labels match those of Fig. 1. Overplotted on both top plots in purple: Contours of the near-IR continuum emission in Fig. 6. Bottom: model fit to the CO features in four selected spaxels. The filled regions are the 1σ uncertainties.

In the text

Fig. 3. H2 1–0 S(3) emission-line morphology (top), velocities (centre), and FWHM line widths (bottom) with the JVLA 8.4 GHz radio continuum shown as thick contours near the centre. Thin contours show the emission-line morphology in the top left panel to ease orientation. Left to right: single-component fit, systemic line component and blue-shifted line component. The emission-line morphology and single-component and systemic velocity maps show an hourglass gas morphology with a strongly blue- and red-shifted bipolar velocity field, consistent with a pair of bubbles that are being inflated by the radio source at the centre. The yellow ellipse and blue circle in this and subsequent figures indicate the beam size of the JVLA and the PSF of the JWST dataset.

In the text

Fig. 4. Source-integrated spectrum of 3C 326 N from the NIRSpec observation from 1.5 to 2.9 μm showing most of the lines. The flux was extracted with an elliptical aperture centred on the peak continuum flux with a major axis along Dec of 1.1″ and a minor axis of 0.9″. Empty pixels were not included and the pixels within the aperture were sigma clipped with a tolerance of 5σ. The spectral gap at 2.2 μm is caused by the NIRSpec IFU detector gap.

In the text

	Fig. 5. Integrated source spectra from the MIRI MRS Channel 1 observation. The flux was extracted with an elliptical aperture centred on the peak continuum flux with a major axis along Dec of 1.1″ and a minor axis of 0.9″. Empty pixels were not included and the pixels within the aperture were sigma clipped with a tolerance of 5σ.
In the text

	Fig. 6. 1.8 μm continuum with contours overplotted.
In the text

Fig. 7. H2 0–0 S(3) emission-line morphology (top), velocities (centre), and FWHM line widths (bottom) with the JVLA 8.4 GHz radio continuum shown as thick contours near the centre. Thin contours show the emission-line morphology in the top left panel to ease orientation. Left to right: single-component fit, systemic line component, and blue-shifted line component. The emission-line morphology and single-component and systemic velocity maps show an hour-glass gas morphology with a strongly blue and red-shifted bipolar velocity field, consistent with a pair of bubbles that are being inflated by the radio source at the centre.

In the text

Fig. 8. Paα emission-line morphology (top), velocities (centre) and FWHM line widths (bottom) with the JVLA 8.4 GHz radio continuum shown as thick contours near the centre. Thin contours show the emission-line morphology in the top left panel to ease orientation. Left to right: single-component fit, systemic line component, and blue-shifted line component. The emission-line morphology and single-component and systemic velocity maps show an hour-glass gas morphology with a strongly blue- and red-shifted bipolar velocity field, consistent with a pair of bubbles that are being inflated by the radio source at the centre.

In the text

	Fig. 9. Single-component [NeIII]λ15.56 emission-line morphology (left), velocities (centre), and FWHM line widths (bottom) with the JVLA 8.4 GHz radio continuum shown as thick contours. Due to strong fringing in MIRI channel 3, the blue-shifted line, while present, is noisy and difficult show clearly. Therefore, we do not plot the two-component fit.
In the text

Fig. 10. [FeII]λ1.644 emission-line morphology (top), velocities (centre), and FWHM line widths (bottom) with the JVLA 8.4 GHz radio continuum shown as thick contours near the centre. Thin contours show the emission-line morphology in the top left panel to ease orientation. Left to right: single-component fit, systemic line component and blue-shifted line component. The emission-line morphology and single-component and systemic velocity maps show an hour-glass gas morphology with a strongly blue- and red-shifted bipolar velocity field, consistent with a pair of bubbles that are being inflated by the radio source at the centre.

In the text

Fig. 11. Comparison of rotational and ro-vibrational line emission for several apertures. Top: integrated flux for the H20–0 S(5) line and the H21–0 S(5) line as titled. The H21–0 S(5) flux image has been convolved to the PSF at the H20–0 S(5) line. The original PSF for H21–0 S(5) and the post-convolution PSF are the white and black filled circles in the H21–0 S(5) plot, respectively. Bottom: ratio of 1–0 to 0–0 line fluxes for the apertures shown in the top image.

In the text

	Fig. 12. Extracted line profiles for a selection of observed lines. The zero velocity is set by the CO band head absorption lines. The lines are extracted using the apertures in Fig. 11 as indicated by the coloured marker in the corner of each subplot.
In the text

Fig. 13. Kennicutt–Schmidt law between star formation rate and molecular gas mass surface density. The large blue and red dot show the position of 3C 326 N when using the mass surface density of cold molecular gas probed by CO(1–0) and H2 0–0 lines, respectively. For the latter estimate, the black line shows the range of surface-brightness observed with MIRI. Small empty circles show several sets of normally star-forming galaxies in the nearby Universe taken from Fischer et al. (2013), Kennicutt (2020), and Saintonge et al. (2013). Green and grey squares show the MOHEGs of Ogle et al. (2010) and Alatalo et al. (2015).

In the text

Fig. 14. Gaussian line width, σ, of warm molecular gas in 3C 326 N as a function of radius. The minimal line widths observed in H2 1–0 S(3) and H2 0–0 S(3) at spatial resolutions corresponding to 190 pc and 530 pc, respectively, are shown as large, upward-pointing arrows. The broadest line widths observed in both lines are > 1000 km s−1, and fall well outside the plot. Black dots show previous SINFONI observations of H2 1–0 S(3) at kpc-scale resolution. The red dotted line shows the position of gas clouds with αvir = 1, for the gas-mass surface densities observed in 3C 326 N. Light blue circles show molecular clouds in the Milky Way. The yellow band indicates the range where the data should fall if the Larson (1981) scaling relations of turbulent clouds apply. The blue line indicates where the lines should fall for typical mass-to-light ratios of gas cooling through turbulent dissipation. See text for details.

In the text