© The Authors 2024

1. Introduction

A pressing question in extragalactic astronomy refers to the role of active galactic nuclei (AGN) in the evolution of their host galaxies. This can only be addressed by quantifying the impact of different types of AGN on their hosts and the evolution of this impact over the life-cycle of an active supermassive black hole. In the case of radio-loud AGN, where radio plasma is emitted in the form of collimated jets, evidence is slowly building up that the young radio jets entering the interstellar medium (ISM) significantly impact the surrounding gas (e.g. Holt et al. 2008; Aditya & Kanekar 2018; Morganti & Oosterloo 2018; Glowacki et al. 2017; Maccagni et al. 2017; Murthy et al. 2021; Shih & Stockton 2014; Kukreti et al. 2022, and references therein).

In order to better understand the contribution of radio AGN in shaping their host galaxies and the evolution of this impact with that of the radio source, it is essential to study radio AGN with different properties and in different evolutionary stages using various probes of the multi-phase ISM at different spatial scales. Although the best way to advance our understanding in this respect is with statistical studies of a large number of sources, these studies are immensely time consuming. Instead, detailed investigations of individual sources – tracing the kinematics, physical condition, and location of the affected gas – in combination with theoretical models of the jet–ISM interaction readily help in obtaining an initial overview of the impact radio AGN have on the surrounding medium. (e.g. Mukherjee et al. 2018a; Audibert et al. 2023; Murthy et al. 2019, 2022; Morganti et al. 2013, 2015, 2021; Dasyra et al. 2015; Cicone et al. 2014; Alatalo et al. 2011; Zovaro et al. 2019). While this has led to significant progress in understanding the impact of radio jets on the ISM, further studies are needed in order to cover a broader parameter space, for example in terms of morphology, radio luminosity, age, and size of the radio sources.

In this context, the 21 cm transition of neutral atomic hydrogen (H I) detected in absorption is a useful probe of the gas in the vicinity of radio jets (see Morganti & Oosterloo 2018, for a detailed review). Apart from the opacity of the gas itself, the detection of atomic hydrogen in absorption depends primarily on the strength and shape of the background radio continuum source. Thus, in the case of radio sources that are bright and extended at parsec (pc) scales, H I 21 cm absorption can be used to study jet–ISM interactions at very high spatial resolution using very long baseline interferometry (VLBI) techniques (e.g. Schulz et al. 2018, 2021).

Here we present new results from a detailed study of atomic gas in the central region of a well-known young radio source, 4C 31.04, which is found to be interacting with the ambient ISM (Zovaro et al. 2019). This is the first part of a broader study of the cold gas in this source, and a study of the molecular gas will be presented in a companion paper (Murthy et al., in prep.).

4C 31.04 is a young radio AGN hosted by the gas-rich galaxy MCG 5-4-18, making it an interesting source for studying the possible interaction between the radio jets and the ISM. It is one of the few sources for which the age has been estimated both from the expansion of the hot spots and a spectral index analysis (Giroletti et al. 2003). This age estimation is in the range of 500−5000 years, making it a genuinely young radio source.

The radio source is ∼140 pc in size¹ with two lobes and a core. It is reported to have a radio luminosity of 1.14 × 10²⁵ W Hz⁻¹ at 5 GHz (Giroletti et al. 2003), making it one of the weakest known compact symmetric objects (CSOs; Polatidis et al. 2002). Although at pc scales it consists of two lobes and a core, making it a CSO, the two lobes are of different morphology: the eastern lobe is the brighter of the two, includes a bright hot spot, and has a well-defined shape, while the western lobe is fainter overall, with a fainter hot spot, and has a relaxed, extended morphology (see Fig. 1, left).

Fig. 1. Left: VLBI radio continuum map of 4C 31.04. The contours start at 4σ and increase by a factor of $\sqrt{2}$. The beam is shown in the bottom right corner of the image. The peak brightness is 196 mJy beam−1, measured at the hot spot along the eastern lobe. Right: HST WFPC2/PC zoomed-in image of the inner part of MCG 5-4-18 in colour overlaid with the radio source in the very centre surrounded by black contours. The HST image shows the warped dust lane perpendicular to the radio axis. The HST image is taken from the Mikulski Archive for Space Telescopes (MAST).

The host galaxy, MCG 5-4-18, is part of a small group, being the dominant member thereof. Hubble Space Telescope (HST) observations (Perlman et al. 2001) show that both the nucleus and the outer regions of the galaxy are partly obscured by dust. In the nuclear region, cone-like features are present aligned with the radio lobes (Fig. 1, right), with the western side being redder than the eastern side. The origin of these features is not discussed by Perlman et al. (2001) and they could simply be due to strong extinction from the molecular disc running north–south. Furthermore, there is another, elongated, warped obscuring feature perpendicular to the radio axis, extending 500 pc to the north and 1 kpc to the south.

A number of studies have shown that this galaxy is rich in cold gas. van Gorkom et al. (1989) first detected a broad H I absorption feature (FWHM > 100 km s⁻¹) towards 4C 31.04, which was later confirmed by Mirabel (1990), who also detected an additional, narrow absorption feature (FWHM ∼ 6.6 km s⁻¹) redshifted with respect to the broad absorption, which they proposed was an infalling high-velocity cloud. Conway (1999) studied the H I absorption towards this source at pc scales using the Very Long Baseline Array (VLBA) and proposed that the broad absorption arises from a disc with an inclination angle of < 20°, a thickness of 100 pc, and a radius of around 500 pc–1000 pc. Struve & Conway (2012), analysing this further, only report a tentative velocity gradient perpendicular to the radio axis and conclude that the absorption arises from the atomic-gas counterpart of the molecular gas disc reported in the literature.

Molecular gas has been detected from this source by various studies. Ocaña Flaquer et al. (2010) report a molecular gas mass of 6.16 × 10⁸ M_⊙ based on the detection of strong CO(1−0) emission and absorption with the IRAM 30 m telescope. García-Burillo et al. (2007) used the Plateau de Bure interferometer to detect a 1.4 kpc molecular gas disc seen in HCO⁺ emission and absorption and suggested that the disc may not be dynamically relaxed owing to a jet–ISM interaction.

A study of warm molecular gas by Zovaro et al. (2019) further supports the jet–ISM interaction scenario. These authors report the detection of shock-excited, warm molecular H₂ out to ∼1 kpc, which they attribute to the far-reaching impact of the jet–ISM interaction, drawing inspiration from the numerical simulations that predict such an effect in the case of radio jets expanding in a clumpy ISM (e.g. Sutherland & Bicknell 2007; Bicknell et al. 2018; Mukherjee et al. 2016). Thus, 4C 31.04 is particularly interesting for high-spatial-resolution studies of cold gas, where this jet–ISM interaction can be traced directly.

We observed 4C 31.04 with the wide band of the Westerbork Synthesis Radio Telescope (WSRT) in 2004 and detected a shallow, broad, blueshifted absorption feature in addition to the broad and the narrow features previously identified. This new feature – an indicator of a fast outflow – triggered more sensitive VLBI observations to localise the outflow and to study its properties and the possible impact of the radio jets on cold gas at pc scales.

We use a redshift of z = 0.0598² for 4C 31.04 based on our molecular gas observations, which will be presented in a companion paper (Murthy et al., in prep.). We detail the observations and data reduction in Sect. 2, report our results in Sect. 3, discuss their implications for the jet–ISM interaction in Sect. 4, and finally conclude in Sect. 5.

2. Observations and data reduction

2.1. WSRT observations

The WSRT observations of 4C 31.04 were carried out on April 11, 2004. The on-source time was five hours and the setup included a bandwidth of 20 MHz subdivided into 256 channels centred at 1341 MHz. The data reduction was carried out using Miriad (Sault et al. 1995) following standard procedures (e.g. Morganti et al. 2016). Antennas C and D, giving the longest baselines in the array, had to be flagged while reducing the data. Although this limits the spatial resolution of the final images and cubes – as 4C 31.04 is unresolved at arcsecond scales – there is no spatial information lost as a consequence of the exclusion of those antennas.

We made the final continuum image and the cube using uniform weighting. The final continuum image has an RMS noise of 1.5 mJy beam⁻¹. The final cube was made at a velocity resolution of 18 km s⁻¹, with an RMS noise of 1.1 mJy beam⁻¹ channel⁻¹ and a beam size of 82″ × 14″ with PA = 4.9°.

2.2. Global VLBI observations

The phase-referenced Global VLBI observations were carried out on 24 May, 2019 with 17 stations of the EVN and VLBA arrays (Project code: GM045B). The two arrays observed this project for 10.3 h each, with an overlap of three hours between the EVN and VLBA antennas. 4C 31.04 was observed for 12.86 h in total with interleaved scans on the phase-reference calibrator, J0123+0344. We used 3C 84 as a fringe finder and bandpass calibrator. The experiment was correlated into two data sets: a continuum pass with a total bandwidth of 64 MHz starting from 1315.75 MHz subdivided into four sub-bands with 32 spectral channels each, and a spectral-line pass of 16 MHz centred at 1339.75 MHz and subdivided into 512 spectral channels.

The data reduction was carried out using standard VLBI data reduction techniques (e.g. Murthy et al. 2019; Schulz et al. 2018, 2021) in “classic” AIPS (Greisen 2003). For both the continuum and line passes, we first applied the parallactic-angle and a priori gain corrections, and initial flagging using the calibration and flag tables obtained from the EVN pipeline and then corrected for the ionospheric dispersive delays. For the continuum pass data, we additionally corrected for the instrumental delays. We subsequently flagged bad data followed by a global fringe fit to correct for phase delay and rate, and then estimated the bandpass solutions. After applying these solutions to the target uv data, we used only the line-free channels to iteratively improve the gain solutions via self-calibration involving imaging and phase-only self-calibration cycles with the solution interval for the phase-only self-calibration being progressively reduced from 30 min to 1 min. We then applied a single round of amplitude and phase self-calibration. Next, we subtracted the continuum from the visibilities and flagged the bad data in the residuals, copied the flag table, and then repeated the self-calibration cycle mentioned above. We continued this until the continuum image showed no further improvement. The final continuum image presented in this paper was made with natural weighting using the continuum-pass data. In the case of the line pass, we subtracted the final continuum model from the visibilities and further subtracted a first-order polynomial from each visibility spectrum to remove any residual continuum emission. We note that due to the fact that the complex structure of the radio source limits accurate continuum subtraction, the noise on the spectra are not perfectly Gaussian. Finally, we deredshifted the uv data using z = 0.0598. Next we imaged each channel using natural weighting to obtain the final spectral cube using the same restoring beam as the continuum image.

The continuum image has an RMS noise of 330 μJy beam⁻¹, and the RMS noise on the cube is 460 μJy beam⁻¹ per 6.98 km s⁻¹ channel. Both have a restoring beam of 12.68 mas × 3.18 mas with PA = −13.25°. The details of observations, images, and cubes are listed in Table 1.

The detected absorption is well resolved against the bright radio source. Hence, we used the Source Finding Application (SoFiA; Serra et al. 2015; Westmeier et al. 2021)³ to extract moment maps from only the broad absorption feature in order to further study the kinematics of the gas. We extracted the moment maps using a mask to include all absorption beyond 4σ in the line channels from the region co-spatial with the radio source. Table 1 lists the observation details as well as the image and cube parameters.

3. Results

3.1. WSRT at kpc scales

The radio continuum is unresolved at the kpc-scale resolution of WSRT and has a peak brightness of 2.67 Jy beam⁻¹, which agrees with the value reported by Mirabel (1990). The H I absorption profile is shown in Fig. 2. For z = 0.0598, the profile includes the broad absorption centred at the systemic velocity and the narrow redshifted absorption at 220 km s⁻¹, which are known from earlier studies. Additionally, we see a shallow, broad blueshifted wing not detected by earlier studies and extending up to −400 km s⁻¹. The broad absorption profile and the shallow blueshifted wing combined have a full width at zero intensity (FWZI) of ∼600 km s⁻¹. We estimated the H I column density using the standard relation: N_H I = 1.82 × 10¹⁸ (T_spin/c_f) ∫ τdv cm⁻², where we assume a gas spin temperature of T_spin = 100 K, the typical temperature of the cold neutral medium (see Morganti & Oosterloo 2018), and the covering factor, c_f = 1. We estimated the total optical depth by integrating over the entire absorption profile, and using a peak brightness of 2.67 Jy beam⁻¹. We obtain N_H I = (1.04 ± 0.09) × 10²¹ cm⁻² for the broad absorption centred at the systemic velocity. This is consistent with the earlier estimations in the literature (e.g. Mirabel 1990).

Fig. 2. Comparison of the WSRT and the global VLBI integrated H I absorption profiles. The WSRT spectrum has a velocity resolution of 18 km s−1 and an RMS noise of 1.1 mJy while the VLBI spectrum has a velocity resolution of 28 km s−1 and an RMS noise of 10 mJy.

For the narrow absorption profile, we obtain N_HI = (5.6 ± 1.1)×10¹⁹ cm⁻². This too is consistent with earlier studies by Mirabel (1990) and Conway (1999) who concluded that the absorber is a foreground gas cloud far away from the nuclear region. As our focus is on the jet–ISM interaction, hereafter we only consider the broad absorption, which, as we discuss below, arises from the nuclear gas co-spatial with the radio source.

3.2. Global VLBI at parsec scales

The pc-scale VLBI continuum map is shown in the left panel of Fig. 1. The radio source has a projected size of ∼140 pc with two lobes and a core, which is in agreement with earlier studies (e.g. Conway 1996; Giroletti et al. 2003). Assuming an uncertainty on the flux density scale of 10%, we measure an integrated flux density of 1.9 ± 0.2 Jy and a peak brightness of 196 ± 20 mJy beam⁻¹. Cotton et al. (1995) report a higher total flux measurement of 2.46 Jy at slightly lower spatial resolution (with a beam size of 17 mas × 10 mas, PA = 72°) with the VLBA at 1.8 GHz. Struve & Conway (2012) however report a significantly lower value of 1.16 Jy. As no variability is observed in 4C 31.04 (Giroletti et al. 2003), we conclude that the lower value of Struve & Conway (2012) is likely due to a combination of the lower sensitivity of their observations and calibration artefacts. The core has a total flux density of 13.5 ± 1.4 mJy. The total flux density within the eastern lobe is 1.0 ± 0.1 Jy and that within the western lobe is 900 ± 90 mJy. We detect both the redshifted narrow feature, and the broad feature centred at the systemic velocity against most of the radio source. The integrated H I spectrum from our observations is shown in Fig. 2. We measure FWZI = 360 km s⁻¹ for the absorption feature at the systemic velocity. We recover all the flux of the absorption detected with the WSRT.

The left panel of Fig. 3 shows the N_HI map of the gas. The column density of the gas along the eastern lobe is higher compared to that along the western lobe with N_{H I,average} = 2.5 × 10²¹ cm⁻² for the former and N_{H I,average} = 7 × 10²⁰ cm⁻² for the latter. At the location of the core, the column density is N_HI = 1.6 × 10²¹ cm⁻². These estimations are for T_spin = 100 K. Struve & Conway (2012) reported a similar difference in the column density between the two lobes. Their values of N_HI for the core and the eastern lobe agree with our estimates, although these authors do not report an estimation for the western lobe.

Fig. 3. Left: H I column density map of the gas seen in absorption in our VLBI data estimated assuming a gas spin temperature of 100 K. Right: velocity field of the gas detected in absorption, shown in colour. There is a slight velocity gradient along the western lobe and also along the parts of the eastern lobe closer to the core in the north–south direction while such a gradient is entirely missing along the rest of the eastern lobe. The contours in both images represent the radio continuum with the same contour levels as those in Fig. 1.

The velocity field of the gas is shown in the right panel of Fig. 3. The gas along the western lobe and parts of the eastern lobe close to the core shows a velocity gradient along PA ∼5° −10°, which is consistent with the molecular gas disc presented by García-Burillo et al. (2007) in terms of the sense of rotation and velocity gradient. On the other hand, the gas in the brightest part of the eastern lobe – including the hot spot – exhibits no significant gradient. We use these features to discuss the morphology of the absorbing gas in Sect. 4.1.

The velocity dispersion of the absorbing gas is shown in the left panel of Fig. 4. The gas seen against the eastern lobe exhibits a higher velocity dispersion compared to that along the western lobe. Interestingly, we detect the shallow, blueshifted absorption wing against only a part of the eastern lobe, ∼31 mas from the core, corresponding to a projected distance of ∼35 pc. The region from which this absorption arises is marked in the left panel of Fig. 4. The width of this broad blueshifted wing is lower than that of the wing detected with the WSRT (see Sect. 4.2). This is due to the lower sensitivity of the VLBI H I spectrum, which prevents us from detecting the faintest part of the blueshifted wing. We further note that, since the shallow absorption feature is detected over only a small region, it is not visible in the integrated spectrum (Fig. 2). We fitted Gaussian components to the shallow and the deep absorption features in the VLBI spectrum extracted from the region containing the blueshifted wing. To quantify the significance of the shallow absorption feature, we used the Bayesian information criterion (BIC). We use BIC = χ² + k ln(n), where χ² is the chi-square statistic of the model, n is the sample size, and k is the number of model parameters. The two-component Gaussian model gives BIC = −1.35, while for the single-component Gaussian model BIC = 7.57, with the absolute BIC difference between the two models being 8.9, which constitutes strong evidence in favour of the two-component Gaussian model (Kass & Raftery 1995). The fits to the shallow and the deep absorption features are shown in the bottom-right panel of Fig. 4. We use this information in Sect. 4.2.

Fig. 4. Left: velocity dispersion map of the H I seen in absorption with our VLBI data, shown in colour. The grey contours represent the radio continuum with the same contour levels as those in Fig. 1. The region where the shallow blueshifted wing is detected is surrounded with a black contour; this was ascertained by visual inspection. Top right: VLBI spectrum extracted from the region marked with the black contour in the velocity dispersion map – where the blueshifted, shallow wing is detected – compared with the WSRT spectrum. The spectra are normalised to emphasise the shallow wing in the VLBI spectrum. Bottom right: Gaussian fits to the deep and the shallow absorption features in the spectrum extracted from the region marked in black contour in the velocity dispersion map. The RMS noise on the VLBI spectrum is 0.7 mJy. See Sect. 4.2 for details.

4. Discussion

We detect multiple H I absorption features in our observations: a narrow feature redshifted by 220 km s⁻¹ with respect to the systemic velocity, a deep, broad absorption centred at the systemic velocity for z = 0.0598, and a shallow, broad, blueshifted absorption wing. The broad shallow wing is a new detection and the high spatial resolution achievable by the Global VLBI array further allows us to localise it. As mentioned above, the focus of this paper is on the gas giving rise to the broad, shallow absorption feature and the feature centred at the systemic velocity. We first discuss the morphology of the gas giving rise to the latter absorption feature in Sect. 4.1 and then discuss the former in Sect. 4.2.

4.1. A disturbed circumnuclear disc

The absorption feature centred at the systemic velocity has a high FWZI = 360 km s⁻¹, and is detected across the entire radio source (see Figs. 2 and 4). These together imply that the gas causing this absorption is close to the radio source and screens it in a smooth distribution at the resolution of our observations (∼4 pc). However, the shallow blueshifted absorption wing, which is indicative of outflowing gas, is detected only towards a small part of the eastern lobe, suggesting the presence of clumps of gas as well.

HST observations by Perlman et al. (2001) show the presence of a dust lane almost perpendicular to the radio axis (i.e. PA ∼ 200°; see Fig. 1). Molecular gas observations further suggest the presence of a warped gas disc that is almost edge-on with a position angle of 5° − 10°, at ∼1.4 kpc (García-Burillo et al. 2007), changing to PA = 350° at ∼10 kpc (Morganti et al. 2023; Murthy et al., in prep.). Thus, it is likely that the H I too belongs to this circumnuclear gas disc.

Conway (1996) and Struve & Conway (2012) previously proposed this scenario based on the difference in H I column density along the two lobes, and the detection of a tentative velocity gradient in H I in their VLBA data along the eastern lobe at PA = 100°, which is almost along the radio axis. Struve & Conway (2012) further estimated a scale height of ∼50 pc for this disc.

With more sensitive observations, however, we find that while the velocity gradient does indeed exist, it is along an axis perpendicular to the radio axis. This can be seen in Fig. 3, the velocity field of H I. The velocity gradient is present along the western lobe and parts of the eastern lobe close to the core at PA ∼ 5° −10°, which shows similarities to the kpc-scale disc reported by García-Burillo et al. (2007). However, this gradient is absent along the rest of the eastern lobe, including the radio hot spot, suggesting that the regular rotation of the gas is significantly disrupted there due to the eastern lobe interacting with the gas disc. Furthermore, the velocity dispersion is high – at ≥40 km s⁻¹, as compared to a typical value of 20 km s⁻¹ in the central regions of early-type galaxies – for most of the gas along both lobes, implying that the entirety of the circumnuclear gas is highly turbulent.

Based on our findings, we suggest that the radio source is expanding into the circumnuclear disc. The eastern lobe is interacting directly with one or more gas clouds, disrupting the regular rotation of the gas locally. The western lobe, although diffuse in morphology, consists of a weaker hot spot, indicating the presence of a milder interaction with the gas. Overall, the gas has a high velocity dispersion everywhere (≥40 km s⁻¹; Fig. 4, left), suggesting that the jet–ISM interaction has affected the disc significantly. Furthermore, the asymmetric morphology of the radio source could be due to the two radio lobes encountering different environments in the process of expansion. The direct jet–cloud interaction occurring along the eastern lobe is also driving an outflow, which we discuss in the following subsection.

4.2. Parsec-scale H I outflow

Given the strong jet–ISM interaction taking place along the eastern lobe, it is very likely that the shallow blueshifted absorption wing arises from an outflowing gas cloud. In the WSRT spectrum, this feature extends up to ∼400 km s⁻¹, while at VLBI scales we only recover the feature extending up to ∼300 km s⁻¹, and the signal is also weaker. As the outflow is resolved even at parsec scales (see Fig. 4), it is likely that the reduced optical depth sensitivity of our VLBI observations is the cause of this partial detection.

This outflow is detected only in a small region along the eastern lobe ∼35 pc from the core and is not detected against the radio hot spot where the optical depth sensitivity is the highest (Fig. 4). This suggests that the outflowing gas is clumpy. Similar clumpy gas outflows at parsec scales have been reported in other young or restarted radio sources (e.g. Morganti et al. 2013; Schulz et al. 2018, 2021).

Following Mahony et al. (2013) and Heckman et al. (2002), we estimate a mass-outflow rate (Ṁ) using the expression:

$$\begin{array}{c}\hfill \dot{M}=30\frac{\mathrm{\Omega }}{4\pi }\times \frac{r_{\ast }}{1\mathrm{kpc}}\times \frac{N_{\mathrm{H}i}}{{10}^{21}{\mathrm{cm}}^{-2}}\times \frac{v}{300\mathrm{km}{\mathrm{s}}^{-1}}{M}_{\odot }{\mathrm{year}}^{-1},\end{array}$$(1)

where Ω is the solid angle that the outflow subtends which we assume to be equal to π, r_* is the radius of the outflow which we measure from the core to be 35 pc, N_HI is the H I column density, and v is the velocity of the outflow.

We use the Gaussian components of the VLBI absorption spectrum shown in the bottom-right panel of Fig. 4, from which we obtain an FWHM for the outflow of 134 km s⁻¹ with the velocity centroid being −300 km s⁻¹. We obtain a peak absorption of 1.0 mJy beam⁻¹ from the fit. The total flux from the region where the outflow is seen is 210 mJy. We estimate an integrated optical depth for the outflow using the approximation: ∫τdv = 1.06 × FWHM × τ_peak = 0.68 km s⁻¹. Closer to the AGN and with disturbed kinematics, as in this case, the spin temperature of the gas is generally higher than the typically assumed 100 K (see Morganti & Oosterloo 2018). Thus, assuming T_spin = 1000 K for the outflow gives a total H I column density of N_HI = 1.2 × 10²¹ cm⁻². This gives us a mass-outflow rate of ∼0.3 M_⊙ year⁻¹. As we do not detect the entire outflow, our estimate of Ṁ is only a lower limit. If we assume that the entire blueshifted wing seen in the WSRT spectrum arises from the outflowing gas in the same region then we obtain an upper limit to the mass-outflow rate. Fitting Gaussian components to the WSRT spectrum gives a peak absorption of 4.0 mJy beam⁻¹, FWHM = 218 km s⁻¹, and a velocity centroid of −205 km s⁻¹ for the blueshifted wing. Using these, and assuming T_spin = 1000 K, we obtain Ṁ = 1.4 M_⊙ year⁻¹. Ocaña Flaquer et al. (2010) report a star formation rate (SFR) of 4.9 M_⊙ yr⁻¹ based on far-infrared measurements. Thus if the H I mass-outflow rate is close to the derived upper limit of 1.4 M_⊙ yr⁻¹, it could have an impact on quenching the star formation in the galaxy. However, the reported SFR is an average over the entire galaxy. Our observations, in comparison, cover only a small region of the centre of the galaxy. Thus, given the different scales, a fair comparison between the two quantities is not possible. We can also obtain an order-of-magnitude estimate of the mass of the gas constituting the outflow, assuming a spherical geometry for the gas cloud. The radius of the gas cloud measured based on the size of the region where the outflow is detected is ∼10 pc giving us a mass of ∼10⁴ M_⊙. The total gas mass in the circumnuclear region of this galaxy is ∼6 × 10⁸ M_⊙ (Ocaña Flaquer et al. 2010). Thus, the outflow does not constitute the bulk of the gas in the circumnuclear region in this case. This has been shown to be case in other sources with pc-scale VLBI H I absorption studies (e.g. 4C 12.50, 3C 236), where the detected gas clouds have masses of ∼10⁵ M_⊙ and the pc-scale outflows (at most tens of M_⊙ year⁻¹) constitute only a small fraction of the circumnuclear gas (e.g. Morganti et al. 2013; Schulz et al. 2021).

4.3. Jet-induced feedback

H I outflows have only been studied from kpc to pc scales in four sources so far, namely 4C 12.50, 3C 236, 4C 52.37, and 3C 293 (Morganti et al. 2013; Schulz et al. 2018, 2021). The estimated ages of these sources range between ∼10⁴ and 10⁵ years, and in all these sources, the H I is found to be clumpy at pc scales. In the younger sources, the entire outflow is concentrated in the central few pc, while in relatively older sources, only part of the outflow is recovered at pc scales, with the rest of the outflow being spread out over hundreds of pc, resolved out at VLBI scales. 4C 31.04 is at least an order of magnitude younger than these sources, where – while gas clumps are indeed present – most of the H I is in a continuous distribution even at a resolution of 4 pc. In all these cases, the pc-scale outflows do not appear to be the dominant form of negative feedback from radio jets. Instead, the main impact of the radio jets on the ambient ISM is the increased turbulence of the gas present in the circumnuclear region.

These findings, albeit drawn from only a handful of sources, agree with numerical simulations of jets expanding into gas discs, which predict an evolution in the nature of the jet–ISM interaction over the age of a radio source (e.g. Wagner & Bicknell 2011; Wagner et al. 2012; Mukherjee et al. 2018b; Dutta et al. 2023). Moreover, Sutherland & Bicknell (2007) claim that the radio morphology of 4C 31.04 matches that of a radio source in one of the earliest phases of evolution in their simulations. These authors predict that radio jets in this case are interacting strongly with an asymmetric ISM, creating an expanding bubble of plasma and gas. The western jet is about to break out of this bubble after interacting with the ISM, while the interaction is still strong along the eastern lobe where the radio jet is still within the created bubble.

Our results strongly support this scenario. We find that in 4C 31.04 the pc-scale radio lobes are expanding into circumnuclear gas and are interacting directly with an asymmetric but predominantly smoothly distributed ISM along the two lobes. Along the eastern lobe, the jet has encountered gas clouds and the resulting interaction has disturbed the gas and is also driving a cold gas outflow at ∼35 pc from the core, in projection. The optical spectrum of 4C 31.04 exhibits only a few weak emission lines (Marcha et al. 1996). Moreover, Zovaro et al. (2019) found that the line ratios of the near-infrared lines are consistent with shock excitation rather than photoionisation. Thus, the effect of radiation on the ambient gas is not likely to be significant, making the radio jets the main cause of the disturbed kinematics of the gas. Also, we find no strong jet–ISM interaction along the western lobe, although it has a hot spot and the gas in this region is highly turbulent (σ ≥ 40 km s⁻¹), suggestive of a strong interaction in the past.

These simulations also predict that the radio plasma that percolates through the disc during the expansion of the jets drives shocks into the ISM over spatial scales that are larger than the apparent extent of the radio source. Although H I absorption observations are very useful for studying the kinematics of the gas at very high spatial resolution, the spatial extent of the gas that can be studied is limited to the size of the background radio source. Thus, to understand the impact of the radio jets at larger spatial scales, this needs to be combined with complementary emission studies that do not suffer from the same limitation. In the case of 4C 31.04, warm molecular gas observations by Zovaro et al. (2019) have already shown that the gas exhibits clear signatures of shock excitation even at kpc scales owing to the impact of the radio jets, which is in agreement with the aforementioned numerical simulations.

4.4. Survival of atomic gas in subkpc-scale outflows

Although these simulations are able to reproduce the observed kinematics of the circumnuclear gas remarkably well (e.g. Mukherjee et al. 2018a,b; Zovaro et al. 2019; Murthy et al. 2022, and references therein), they only consider gas at very high temperatures (> 10⁵ K), while the atomic and molecular gas that is a part of the ISM directly interacting with the radio jets is expected to be at lower temperatures (< 10⁵ K; e.g. Perucho 2024). Thus, it still needs to be verified whether the conditions defined in these simulations are conducive for the survival of cold gas and, as such, are a good representation of the actual situation.

Perucho (2024) provide an analytical solution to this difficult problem with a focus on atomic gas outflows. These authors predict that the shocks driven by low-power radio jets (≤10⁴³ erg s⁻¹) are not strong enough to ionise entire atomic gas clouds at any spatial scale at any given time. Moreover, the cooling time for the ionised gas in the clouds is almost instantaneous compared to the lifetime of the source. Thus, atomic gas should exist in outflows even within the central subkpc region affected by the expanding radio jets either via direct interaction or shocks.

On the other hand, the cold gas may or may not be detected at subkpc scales in higher power sources (> 10⁴³ erg s⁻¹) depending on the density of the gas. Such sources are capable of ionising a large fraction of the gas and also cause the cooling times to be longer due to the mixing of the clouds with the hotter environment. At low gas densities (< 10³ cm⁻³), this would result in the atomic gas outflows being detected only at kpc scales. However, when the density of the gas is higher (> 10³ cm⁻³), the post-shock temperature is lower and the cooling times are very short compared to the lifetime of the source (see Fig. 4 in Perucho 2024). This results in rapid recombination allowing the atomic gas to exist even in the extreme conditions at subkpc scales.

We estimate the jet power of 4C 31.04 to be in the range of 6 × 10⁴³ erg s⁻¹ to 10⁴⁴ erg s⁻¹ using the relations presented in Cavagnolo et al. (2010) and Wu (2009). As these estimations could at best be a lower limit (e.g. Mukherjee et al. 2018b), we note that 4C 31.04 is an intermediate or high-power radio source. Thus, in the scenario presented by Perucho (2024) the presence of a subpc outflow in 4C 31.04 indicates that the gas is sufficiently dense for the cooling times to be very short and the recombination rate to be high, while the mixing between the cold and hotter gas via instabilities is also likely to be inefficient.

Thus, overall, these results add to the growing body of evidence that radio jets have a significant impact on cold gas in the early stages of their evolution. In the very central region, radio jets primarily make the gas highly turbulent, while also driving outflows, albeit only constituting a small fraction of the total circumnuclear gas. Furthermore, the overall impact on the gas spans a much larger spatial scale – even an order of magnitude larger in the case of 4C 31.04 – than the apparent extent of the radio source itself. Similar studies of radio sources spanning a wide range in age, radio luminosity, and size are needed to build a clearer picture of the nature of the jet–ISM interaction and the impact it has on the host galaxies.

5. Summary

We present deep H I absorption observations at kpc and pc scales using the WSRT and the global VLBI array of a young radio AGN, 4C 31.04, a well-known CSO with radio lobes extending up to ∼140 pc. At kpc scales, we detect a narrow absorption (FWZI = 6.6 km s⁻¹) redshifted by 220 km s⁻¹ with respect to the systemic velocity estimated using z = 0.0598, a broad absorption centred at the systemic velocity, and a shallow blueshifted wing with a combined FWZI of 525 km s⁻¹. This study is focused on the broad absorption at the systemic velocity and the blueshifted wing, as these correspond to the gas in close proximity to the radio lobes. We measure a radio continuum flux density of (2.67 ± 0.267) Jy and a total H I column density of the gas – excluding that giving rise to the narrow absorption profile – of (1.04 ± 0.09)×10²¹ cm⁻².

At pc scales with our Global VLBI observations, we detect the redshifted absorption and the broad absorption at the systemic velocity. We measure a total continuum flux density of (1.19 ± 0.2) Jy with our VLBI data. The absorption at the systemic velocity is well resolved at the resolution of our VLBI observation and detected across the entire radio source. We also detect a part of the blueshifted wing but only against a small region of the eastern lobe, ∼35 pc from the core, just prior to the radio hot spot. The gas has high velocity dispersion everywhere (≥40 km s⁻¹) and the region of the highest velocity dispersion (∼90 km s⁻¹) includes the location of the outflow and the region containing the radio hot spot along the eastern lobe. We detect a velocity gradient along the western lobe and parts of the eastern lobe close to the core at PA ∼ 5° −10°, while such a gradient is considerably weaker along the rest of the eastern lobe. We also find that the H I column density is higher along the eastern lobe, with an average N_HI = 2.5 × 10²¹ cm⁻² compared to the western lobe where the average N_HI is 7 × 10²⁰ cm⁻².

Combining the above findings, we arrive at the conclusion that the lobes of 4C 31.04 are expanding into a circumnuclear disc, partially disrupting it along the eastern lobe, and overall, making the gas highly turbulent. While the gas is predominantly in a smooth distribution at the resolution of our data (∼4 pc), being detected all across the radio source, we also find the presence of clumps at least along the eastern side. Here, the eastern lobe is strongly interacting with gas clouds and also driving an outflow of an H I cloud. We estimate a mass-outflow rate of 0.3 ≤ Ṁ ≤ 1.4 M_⊙ year⁻¹. We detect only a part of the outflow compared to that detected by the WSRT. We suggest that this is due to lower optical depth sensitivity of our observations on account of the outflow being resolved.

Our results support the scenario proposed by Sutherland & Bicknell (2007) where 4C 31.04 is in a very early stage of expansion into an asymmetric ISM. More broadly, our results also agree with the predictions of similar theoretical models regarding the nature of the ISM and the nature of jet–ISM interaction in young radio sources expanding into the circumnuclear gas.

Furthermore, we compare our observations with the predictions of a recent analytical model of the impact of radio jets on atomic gas by Perucho (2024). We find that the existence of a subkpc-scale outflow in this case could be due to inefficient mixing of the cold gas with the hot medium, and high gas density leading to an almost instantaneous recombination of the ionised gas.

Overall, our study demonstrates the importance of probing the jet–ISM interaction at multiple spatial scales and supports the findings presented of earlier studies that radio AGN in the early stages of their evolution significantly impact the gas into which they are expanding, especially that in cold phase.

Acknowledgments

We thank the referee for insightful comments that improved the presentation of the paper. We also thank Manel Perucho for an interesting and useful discussion on sub-kpc H I outflows. S.M. thanks the organisers of the Lorentz Centre workshop, “The importance of jet-induced feedback on galaxy scales” where interesting discussions regarding some of the results presented in this paper were held. The European VLBI Network is a joint facility of independent European, African, Asian, and North American radio astronomy institutes. Scientific results from data presented in this publication are derived from the following EVN project code: GS045. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Westerbork Synthesis Radio Telescope was operated by ASTRON (Netherlands Institute for Radio Astronomy) with support from the Netherlands Foundation for Scientific Research (NWO).

References

All Tables

All Figures

Fig. 1. Left: VLBI radio continuum map of 4C 31.04. The contours start at 4σ and increase by a factor of $\sqrt{2}$. The beam is shown in the bottom right corner of the image. The peak brightness is 196 mJy beam−1, measured at the hot spot along the eastern lobe. Right: HST WFPC2/PC zoomed-in image of the inner part of MCG 5-4-18 in colour overlaid with the radio source in the very centre surrounded by black contours. The HST image shows the warped dust lane perpendicular to the radio axis. The HST image is taken from the Mikulski Archive for Space Telescopes (MAST).

In the text

	Fig. 2. Comparison of the WSRT and the global VLBI integrated H I absorption profiles. The WSRT spectrum has a velocity resolution of 18 km s−1 and an RMS noise of 1.1 mJy while the VLBI spectrum has a velocity resolution of 28 km s−1 and an RMS noise of 10 mJy.
In the text

Fig. 3. Left: H I column density map of the gas seen in absorption in our VLBI data estimated assuming a gas spin temperature of 100 K. Right: velocity field of the gas detected in absorption, shown in colour. There is a slight velocity gradient along the western lobe and also along the parts of the eastern lobe closer to the core in the north–south direction while such a gradient is entirely missing along the rest of the eastern lobe. The contours in both images represent the radio continuum with the same contour levels as those in Fig. 1.

In the text

Fig. 4. Left: velocity dispersion map of the H I seen in absorption with our VLBI data, shown in colour. The grey contours represent the radio continuum with the same contour levels as those in Fig. 1. The region where the shallow blueshifted wing is detected is surrounded with a black contour; this was ascertained by visual inspection. Top right: VLBI spectrum extracted from the region marked with the black contour in the velocity dispersion map – where the blueshifted, shallow wing is detected – compared with the WSRT spectrum. The spectra are normalised to emphasise the shallow wing in the VLBI spectrum. Bottom right: Gaussian fits to the deep and the shallow absorption features in the spectrum extracted from the region marked in black contour in the velocity dispersion map. The RMS noise on the VLBI spectrum is 0.7 mJy. See Sect. 4.2 for details.

In the text