© The Authors 2023

1. Introduction

Active galactic nuclei (AGN) release large amounts of energy, which impacts the evolution of the host galaxy. This phenomenon, known as AGN feedback, is important for understanding galaxy evolution (for example, Silk & Rees 1998; Fabian 2012; Harrison et al. 2014; McNamara & Nulsen 2012). Tracing the impact of the various forms of nuclear activity on the host galaxy requires observing its effect on the surrounding gas. One signature of this impact is the presence of kinematically disturbed gas and gas outflows.

Evidence for outflows has been found for many AGN in ionised, neutral atomic, and molecular gas (for example, Cicone et al. 2014; García-Burillo et al. 2014; Harrison et al. 2018; Husemann et al. 2019; Maksym et al. 2022; Morganti & Oosterloo 2018; Müller-Sánchez et al. 2011; Trindade Falcão et al. 2021; Veilleux et al. 2013, 2020; Winkel et al. 2023). Broadly, two mechanisms have been suggested to explain the origin of the outflows. The first mechanism is radiative or quasar mode feedback, which occurs when gas is ejected in winds due to radiation pressure from the nucleus. The second mechanism is mechanical or radio mode feedback, which occurs due to radio jets’ interaction with the surrounding gas (for example, Fabian 2012; King & Pounds 2015). Both mechanisms are likely to play a role in driving feedback in different AGN. Indeed, radiative mode is expected to be dominant in high luminosity sources, close to Eddington, whereas mechanical mode is found more often in sources with luminosity much lower than Eddington. Therefore, it is important to quantify their impact in different AGN. Another important parameter is the changes in the impact of feedback during the AGN lifetime. This would provide a better understanding of AGN feedback, and it is also an important input for numerical simulations of galaxy evolution.

In radio AGN, many studies have found evidence for jet-driven feedback in AGN host galaxies (for example, Husemann et al. 2019; Jarvis et al. 2019; Morganti et al. 2021a; Murthy et al. 2022; Nesvadba et al. 2007, 2021; Ruffa et al. 2022; Schulz et al. 2021; Venturi et al. 2021 for some recent results). Some trends have also been found for the possible impact of jets on the surrounding gas. For example, in a sample of optically selected AGN, Mullaney et al. (2013) found that 1.4 GHz radio luminosity had the most influence on the ionised gas kinematics of the narrow line region, compared to the bolometric luminosity. Therefore, they could disentangle the impact of jets from radiation. A positive link between young radio AGN (and thus a young jet) and disturbed gas kinematics has also been suggested by both ionised and neutral gas studies (Holt et al. 2008; Geréb et al. 2015; Roche et al. 2016; Maccagni et al. 2017; Molyneux et al. 2019; Murthy et al. 2019; Santoro et al. 2020), which is in line with the picture where feedback evolves over the AGN lifetime. This is supported by numerical simulations of an interaction between jets and a clumpy medium, which allowed authors of previous studies to conclude that the impact changes as the jets evolve to larger scales (Sutherland & Bicknell 2007; Mukherjee et al. 2018; Meenakshi et al. 2022).

While Mullaney et al. (2013) could disentangle the effect of radiative and mechanical mode feedback, emphasising the importance of the latter for driving gas outflows, the connection to young radio AGN is still mostly limited to powerful sources or small groups. Furthermore, most of the studies use high frequency (GHz range) radio morphology as a proxy for the evolutionary stage of radio AGN (for example, Holt et al. 2008; Mullaney et al. 2013; Molyneux et al. 2019). However, as we explain below, radio spectral properties can give more accurate information about the evolutionary stage of the AGN.

In this paper, we investigate a sample of radio sources for which we characterise the evolutionary stage using the radio spectral information, extended to low frequency (∼150 MHz) by LOFAR. We combine this with the properties of ionised gas provided by SDSS, thus studying the presence of outflows and quantifying the feedback effect. The sample size allows us to perform a stacking analysis to investigate the presence of outflows in groups of sources selected based on their optical and radio properties, so to disentangle the relevance and effect of various key parameters (for example age of the radio source, Eddington rate, excitation state etc.).

Over the years, radio AGN have been observed in phases of activity and dormancy, which is collectively called the radio AGN life cycle (see Morganti 2017 for a review). Young radio AGN are expected to show a peak in their radio spectrum (usually in the GHz frequency range), due to absorption of the synchrotron emission, and bright jets on small scales (sub-kpc to kpc). The absorption is either due to synchrotron self-absorption (SSA) or free-free absorption (FFA). This phase lasts 10⁵ − 10⁶ Myr and has been detected in peaked spectrum sources (and in particular compact steep spectrum, CSS, and gigahertz peaked spectrum, GPS, sources; Slob et al. 2022 and see O’Dea 1998; O’Dea & Saikia 2021 for a review). As the source evolves and the jets grow, the peak in the radio spectrum shifts to lower frequencies, eventually below the observing frequencies of most radio telescopes. At this stage, the radio spectra is observed to be optically thin. However, with low frequency radio telescopes such as the LOFAR telescope, it is possible to detect the peak down to ∼150 MHz, allowing us to detect peaked sources that would be missed in higher frequency data.

After a phase of activity which can last up to 10⁸ yr, the fuelling of the radio jets stops and there is no more injection of fresh electrons. This can lead to a remnant phase with low surface brightness diffuse emission, and a steep radio spectra, due to preferential radiative cooling at high frequencies (∼GHz, for example, Saripalli et al. 2012; Brienza et al. 2017; Jurlin et al. 2021). Since the lower frequency end of the radio spectrum is affected last, LOFAR is an ideal instrument to detect extended radio emission from AGN, even in the remnant phase. Such low frequency data is also useful to accurately characterise source morphology. In some cases (13−15%), a new phase of activity can be detected even before the radio plasma from the older phase has faded below detection limits (Jurlin et al. 2020). These sources are called restarted radio AGN (for example, Parma et al. 2007; Shulevski et al. 2015; Brienza et al. 2018; Morganti et al. 2021b,c, and see Saikia & Jamrozy 2009 for a review). These sources are expected to have a young phase of activity in the centre, characterised by a peaked spectrum, surrounded by diffuse emission from an older phase of activity (for example in 3C 293 by Kukreti et al. 2022). Therefore, the radio spectral shape can be used as a good tracer for the AGN life cycle.

Tracing the radio spectral shape accurately requires a wide frequency coverage. To achieve this, we used a combination of radio surveys at three frequencies – the LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017) at 144 MHz, the Faint Images of the Radio Sky at Twenty-cm survey (FIRST; Becker et al. 1995) at 1400 MHz and the Very Large Array Sky Survey (VLASS; Lacy et al. 2016) at 3000 MHz. These surveys were also chosen for their high angular resolution (6″ for LoTSS, 5.4″ for FIRST, and ∼2.5″ for VLASS), which ensures that we trace the spectral shape of the central few kpc (≈6−12 kpc at the mean redshift of the sample), where the emission is expected to be dominated by the radio jets. This is needed in order to relate the radio spectral shape to the evolutionary stage of the radio jet. The resolution of these surveys is also comparable to SDSS fibre size (∼3″), which allows us to compare the radio and optical emission over similar spatial scales.

Our goal is to connect the AGN life cycle to feedback. To trace the feedback on the ionised gas and its kinematic properties, we used the [O III]λλ4958,5007 Å doublet emission line. We used the [O III] line because it is sensitive to the impact of radiation and jets and therefore a good tracer for outflows. It is a forbidden transition line, and therefore can only be produced in the low density environments of the narrow line region (NLR) clouds that surround the AGN. The NLR extends up to kpc sizes, and is thus ideal for studying AGN feedback on larger scales than the high density sub-parsec-sized broad line region. [O III] is also a strong emission line, making it easier to detect for a large number of sources.

The paper is organised in the following manner: in Sect. 2 we describe the sample used for the study, in Sect. 3 we describe the methods to obtain the radio and optical properties, and our analysis, in Sect. 4 we present our results, and in Sect. 5 we discuss the ionised gas kinematics as a function of different source properties. Finally, we summarise our results in Sect. 6. In Appendix A, we provide LOFAR 144 MHz continuum maps for sources where we detect extended emission for the first time. Throughout the paper, we have used the ΛCDM cosmological model, with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3 and Ω_vac = 0.7. At the average redshift of our sample (z = 0.11), 1″ corresponds to 2 kpc.

2. Sample description

For our study we used the sample from Maccagni et al. (2017, hereafter M17) and Geréb et al. (2014, hereafter G14). The sample was constructed by cross-correlating the Sloan Digital Sky Survey catalogue (SDSS DR7; York et al. 2000) with the FIRST catalogue. The sample was radio flux selected to have a peak flux density S​_1.4 GHz > 30 mJy in the redshift range 0.02 < z < 0.23 and it had 248 sources. This radio flux density is not the total flux of the source, but the flux of the core coincident with the SDSS counterpart.

We cross-matched this sample with the LoTSS DR2 (Shimwell et al. 2022) and the VLASS (Gordon et al. 2021) catalogue. The cross-match was done using Topcat’s nearest neighbour routine Taylor (2005) with a radius of 7″, comparable to the LoTSS resolution. After the cross-match, our final sample has 135 sources covering a radio luminosity range of L_1.4 GHz = 10^22.6 − 10^26.3 W Hz⁻¹. The sample contains six sources that are known galaxy mergers, with signatures of an ongoing interaction. Since we want to investigate the gas kinematics of radio AGN only, we exclude these six sources from our sample. Finally, we have a sample of 129 radio AGN. The redshift, 1.4 GHz radio luminosity and observed [O III] luminosity distributions (Sect. 3.2.1) of the sample are shown in Fig. 1.

Fig. 1. Sample properties. Left: redshift and (middle) peak 1.4 GHz radio luminosity distributions of the 135 sources in the sample (including mergers). The dip in the number of sources in the interval 0.136 < z < 0.175 is due to the exclusion of sources affected by RFI (G14). Right: [O III] luminosity distribution for all the sources. For non-detections, we used the upper limit of the [O III] luminosity. The hatched area shows the distribution for only the [O III] detections.

3. Methods

3.1. Radio properties

3.1.1. Radio spectral indices

The combination of the high resolution LoTSS, FIRST and VLASS surveys, gives us a unique opportunity to study the radio properties of the same central region for which we have the optical properties from SDSS, thus allowing a direct comparison between the two. To do this, we have used the peak flux densities from the surveys at the central region of the sources, to estimate the spectral indices. Throughout the paper, spectral index refers to this central region spectral index.

The angular resolution of VLASS (∼2.5″) is significantly higher than LoTSS (6″) and FIRST (5.4″). Therefore, we smoothed all the VLASS images for our sources using a 2D Gaussian with the task IMSMOOTH in CASA, to achieve an angular resolution of 5.4″, that is the same as FIRST. This causes a significant change in spectral index for extended sources, since they have emission that is missed if we only use the peak flux density seen by the ∼2.5″ VLASS beam. Finally, Gordon et al. (2021) have also reported that the VLASS flux densities are underestimated by a factor of 0.87. Therefore, we have corrected the VLASS flux densities by this factor. We have then estimated the spectral indices from 144 MHz to 3000 MHz, using the convention S ∝ ν^α. We estimate the spectral index errors using a quadrature combination of rms noise and flux scale errors of the surveys. The flux scale errors used are 10% for LoTSS, 5% for FIRST and 10% for VLASS. The typical error on the spectral indices are 0.05 for ${\alpha }_{1400}^{144}$ and 0.15 for ${\alpha }_{3000}^{1400}$.

3.1.2. Radio spectral shape classification

Our aim is to link the AGN life cycle to feedback. The AGN life cycle can be probed using the radio spectral shape. To characterise the spectral shape, we have plotted our sources on a colour–colour plot of spectral indices, from 144−1400 MHz and 1400−3000 MHz, shown in Fig. 2. Depending on the location of a source on this plot, we can identify its radio spectral shape. As mentioned in the introduction, in the very early stages of evolution, the radio spectra of a radio AGN is peaked at GHz frequencies. As the source evolves, the spectral peak moves to lower frequencies and both optically thick and thin parts of the spectra can be observed. Ultimately, the radio source becomes completely optically thin and the spectral peak cannot be detected anymore. With our wide frequency coverage from 144 MHz to 3000 MHz, we can track the spectral peak over a large range, and identify sources in different stages of evolution.

Fig. 2. Colour–colour plot of spectral indices for the central region of the 129 radio AGN. x-axis shows ${\alpha }_{1400}^{144}$ derived using LoTSS and FIRST peak flux densities and y-axis shows ${\alpha }_{3000}^{1400}$ derived using FIRST and VLASS peak flux densities. Regions in the plot for peaked, non-peaked and convex spectrum sources are marked with blue, orange and grey colour respectively. The spectral shape of a source in each quadrant is drawn at the corners. The error bar shows the typical error on the spectral indices. The hatched design marks the area which we classify as flat-spectrum region. Sources lying in this region are not included in peaked, non-peaked or convex spectrum groups. [O III] detections are marked with filled circles, and non-detections are marked with hollow circles.

Sources with ${\alpha }_{1400}^{144}>0$ and ${\alpha }_{3000}^{1400}>0$ have an absorbed spectra with a spectral peak above 3000 MHz. These sources are in the ‘youngest’ stage of evolution. As the spectral peak moves to lower frequencies, ${\alpha }_{3000}^{1400}$ becomes optically thin (< 0) while ${\alpha }_{1400}^{144}$ is still optically thick (> 0). These are GPS and CSS sources. They lie in the blue region of Fig. 2 and we label them collectively as peaked spectrum sources. Although GPS and CSS sources are usually defined to have a high frequency spectral index α ≤ −0.5 (O’Dea 1998), this definition is arbitrary and the high frequency spectral index has been observed to have a continuous distribution (Fanti et al. 1990; Stanghellini et al. 1998; Callingham et al. 2017).

When the spectral peak moves to frequencies below 144 MHz, that is ${\alpha }_{1400}^{144}<0$ and ${\alpha }_{3000}^{1400}<0$, the radio spectra is completely optically thin in our frequency range. These are the most evolved sources in our sample. They lie in the orange region of Fig. 2 and we label them non-peaked sources. It is possible that there are still some sources in this region that have a peak between 144 MHz and 1400 MHz, and more frequency points are needed to identify them. Sources in the grey region have a convex spectrum, that is ${\alpha }_{1400}^{144}<0$ and ${\alpha }_{3000}^{1400}>0$, and could be sources with multiple epochs of activity, where the optically thick high frequency spectral index is due to the newer phase of activity.

In Fig. 2, sources within the hatched region are not included in any spectral shape group. Sources in this region can have a flat spectra either due to the core dominating the emission, or relativistic beaming due to a small jet viewing angle. These sources are called flat-spectrum radio AGN, and are usually defined in literature to have a spectral index α ≤ 0.5 at all frequencies (O’Dea 1998). In the plot, this region is within spectral indices $-0.5\ge {\alpha }_{1400}^{144}\le 0.05$ and $-0.5\ge {\alpha }_{3000}^{1400}\le 0.15$. This includes sources with spectral indices within 1σ of both axes. It is also possible that there are some genuine peaked spectrum sources in this region, but we cannot identify them with the current frequency points. All spectral shape groups and their properties are summarised in Table 1.

Although peaked spectrum sources are understood to be young radio AGN, another explanation used for these sources is the frustration hypothesis. In this case, peaked sources are not young AGN but confined to small scales by a dense ambient medium (for example, van Breugel et al. 1984; O’Dea et al. 1991; Callingham et al. 2015). Such a dense ambient medium can cause absorption of the radio emission by FFA. Therefore, we also check further in Sect. 4.2.1 if the peaked sources in our sample are genuine young sources, or absorbed due to a high density medium.

3.1.3. Radio morphology and power

The sample has been classified into compact and extended radio sources by G14, using a combination of 1.4 GHz NRAO VLA Sky Survey (NVSS; Condon et al. 1998) major-to-minor axis ratio and FIRST peak-to-integrated flux ratio (see Sect. 2.1 of their paper). However, the low frequency and high sensitivity of LoTSS to diffuse low surface brightness emission means that it can detect low level extended emission, which could be missed at 1.4 GHz. Therefore, we used the criteria of Shimwell et al. (2022) for LoTSS-DR2, to classify a source as resolved or unresolved. Their criteria (described in detail in Sect. 3.1 of their paper) uses a combination of signal-to-noise ratio and integrated-to-peak flux ratio to separate sources. This method attempts to account for calibration inaccuracies and signal-to-noise ratio of the sources, and is more robust than just using the catalogued source sizes. However, it is important to note that the LoTSS criteria was designed for a statistically large sample, and can have field-to-field variations leading to a small percentage of sources being misclassified.

At the average redshift of our sample, an unresolved source with 6″ resolution of LoTSS corresponds to a physical size upper limit of ∼12 kpc. This is close to the typical upper limit of ∼20 kpc for CSS sources (for example, O’Dea 1998; Fanti et al. 1990). Therefore, we consider a source to be compact only if it is classified as unresolved by Shimwell et al. (2022) and compact by G14 criteria. All other sources are classified as extended. A visual inspection of the compact sources revealed that eight sources show clear extended emission in their LoTSS images, on scales of a few tens to hundred kpc. This could be due to the caveat in using the LoTSS criteria mentioned in the paragraph above. However, these source do not show any extended emission in their FIRST or VLASS images. Therefore, these sources could be restarted radio AGN, as mentioned in the introduction, and we discuss them in Sect. 5.4. These sources were also moved to the extended group and are shown in Fig. A.1. In summary, our sample has 64 compact and 65 extended radio sources. Out of the peaked sources discussed in the previous section, 68% are compact and 32% are extended. Comparing this to the non-peaked sources, where 52% are compact and 48% are extended, shows that peaked sources tend to be more compact than non-peaked sources, as would be expected from the AGN life cycle scenario.

We also calculated the 1.4 GHz radio power of the sources using their peak flux density from FIRST. The distribution of the radio powers is shown in Fig. 1. To assess the impact of radio power on feedback, we split the sample along the median radio power. The groups are summarised in Table 1.

3.2. Optical properties

3.2.1. Stellar continuum and emission line modelling

The stellar continuum and emission line modelling (using an F-test approach) presented for the sample was done in Santoro (2018). In this section, we describe the procedure to obtain the emission line spectra of the sources. We have taken the optical spectroscopic data from SDSS DR12 (Alam et al. 2015). The spectra cover the wavelength range 3800−9200 Å with an average spectral resolution of about 160 km s⁻¹. The accuracy on the wavelength calibration is < 5 km s⁻¹. The SDSS pipeline provides properties of the emission lines, but only with a single-component fit to the line profile. A first visual inspection of the spectra already indicated that, when emission lines were present, they often had complex line profiles which required multiple kinematical components for proper modelling. For this reason we re-processed the SDSS spectra with our fitting procedure to extract more reliable information on the kinematics and excitation of the ionised gas.

Considering that in a number of sources the emission lines are weak and the stellar continuum dominates, we performed a detailed fitting of the stellar continuum. The modelling was done using pPXF (Cappellari 2017) in combination with single stellar population models (SSP) from the 2016 version of the spectral synthesis results of Bruzual & Charlot (2003). We used stellar populations with old (11 and 12.3 Gyr), intermediate (0.29, 0.64, 0.90, 1.4, 2.5 and 5 Gyr) and young (5, 25 and 100 Myr) ages. All the SSP assume a Chabrier initial mass function, solar metallicity and instantaneous starburst. We shifted the spectra to the restframe using the SDSS redshift estimate and, during the fitting procedure, we masked the spectral regions corresponding to the [O II]λλ3726,28 Å, Hδ, Hγ, Hβ, [O III]λλ4958,5007 Å, [O I]λ6300 Å, [N II]λλ6548,84 Å, Hα and [S II]λλ6717,30 Å emission lines. This procedure reproduces the stellar continuum of our sources well, including the cases where strong emission lines are present. An example of a source spectrum and stellar continuum modelling result is shown in Fig. B.2.

After subtracting the best model of the stellar continuum given by pPXF, we performed the modelling of the emission lines using Gaussian functions. As a reference model, we used the [O III]λλ4958,5007 Å doublet which is usually the strongest spectral feature associated with the ionised gas and is not contaminated by other emission lines. Each kinematical component of the model consisted of two Gaussian functions with the same width (that is σ_{[OIII] 4958} = σ_{[OIII] 5007}), fixed separation and relative amplitude ($A_{[\mathrm{OIII}]4958}=\frac{1}{3}\times A_{[\mathrm{OIII}]5007}$) according to atomic physics. We started by fitting each line of the doublet with one Gaussian component and added up to two extra Gaussian components to allow our procedure to fit complex line profiles with broad wings. Every time we added a Gaussian component in our fitting procedure, any change in χ² was evaluated with an F-Test to establish whether this provided a better fit at a confidence level > 99%. Multiple Gaussian components were required to fit the spectra of 21 sources. We considered the [O III]λ5007 Å emission line detected when the peak of the emission line model was higher than three times the standard deviation of the model residuals. If the [O III]λ5007 Å line was detected, we proceeded with fitting the other emission lines.

The best model of the [O III] line was used to fit the Hβ line, that is up to three Gaussian components per emission line. To take into account possible Hβ broad emission coming from the broad line region (BLR) of type 1 AGN, we allowed one more Gaussian component in the fit of the Hβ line if needed. In our sample, two sources required an additional broad component, and are classified as BLR. Finally, we fit the [N II]λλ6548,84 Å and the Hα lines together. Each kinematical component of this model was made by three Gaussian functions with fixed separation according to atomic physics and the same line width. We also fixed the relative line ratio of the two Gaussian functions associated with the [N II] lines ($A_{[\mathrm{N}\mathrm{I}\mathrm{I}]5648}=\frac{1}{3}\times A_{[\mathrm{N}\mathrm{I}\mathrm{I}]5684}$). Both emission lines were fit using the best model for the [O III] line. We included an extra Gaussian component to fit the Hα line when an additional component was needed to fit the Hβ line.

We extracted the fluxes of the [O III], Hα, Hβ and [N II] emission lines considering the summed flux of all the components needed to model the line. If an additional component is needed to fit the Hβ and Hα line, we extracted the flux of these lines only from the components of the [O III] model. In this way our line fluxes were not affected by the light coming from the BLR of the AGN. Overall, in the spectra of the 129 radio AGN, we have 93 [O III] detections. For the non-detections, we estimated an upper limit on the flux and used it to estimate the upper limit for the [O III] luminosity. The luminosity distribution is shown in Fig. 1. In this radio AGN sample, we have 77 Hα detections, 55 Hβ detections and 91 [N II] detections. Fifty five sources have detections for all four lines and 21 sources have detections for [O III], Hα, and [N II] but not Hβ. Fifteen sources have detections for [O III], and [N II] but not Hα and Hβ. One source has detection for [O III], and Hα but not Hβ and [N II]. One source has only a detection for [O III] and no other line.

The [O III] luminosity is also a commonly used proxy for the total radiative luminosity of an AGN. We used it as a proxy for assessing the role of radiation pressure in driving the gas kinematics. To do this, we split the sample about the median [O III] luminosity, L_[OIII] = 1.6 × 10⁴⁰ erg s⁻¹ cm⁻² Å⁻¹ (including the upper limits for non-detections). The groups are summarised in Table 1 and the results shown in Sect. 4.2.3. We did not use the extinction corrected [O III] luminosity because there were very few sources, only 55, with both Hα and Hβ detections. Thus the correction could only be performed for part of the sample. Furthermore, the relation to convert [O III] luminosity to the total radiative luminosity, discussed in Sect. 3.3, uses the observed [O III] luminosity. We note that although we have used the [O III] luminosity as a proxy for the radiative luminosity, for our sample of radio AGN, the [O III] emission is not solely due to AGN photoionisation, but can also have contribution from shock ionisation. This could affect our result for disentangling the role of radiation pressure and is a caveat of this method.

Since the [O III] line modelling requires multiple Gaussian components in some cases, a comparison between them is hard. To compare the kinematics of the individual [O III] detections in a uniform manner, following the approach of Mullaney et al. (2013), we estimated the flux-weighted average full width at half maximum (FWHM) of the [O III] λ5007 Å line (FWHM_avg), defined as

$$\begin{array}{c}\hfill {\mathit{FWHM}}_{\mathrm{avg}}={\displaystyle \sum _i}F(i)\times \mathit{FWHM}(i),\end{array}$$(1)

where F(i) is the flux ratio of the ith component to the total flux in the model, and FWHM(i) is the FWHM of the ith component. This allows for a direct comparison between the different [O III] profiles, regardless of the number of Gaussian functions needed to model the line.

Multi-component model fits can help track non-gravitational motion in the NLR, for example, a gas outflow. For instance, if they have a significantly different velocity than the systemic (core component), or broader widths than expected from gravitational motion. Throughout the paper, whenever a multi-component fit is required to fit the [O III] spectra, we used two criteria to identify the extra (“wing”) component as a gas outflow. These are – (i) the velocity difference from the core component, Δv, is at least 3σ significance, or (ii) the component’s velocity dispersion, σ_disp, is larger than the three sigma upper limit of the stellar velocity dispersion. For the stellar velocity dispersion upper limit, we used the average value for the sample, which is 348 km s⁻¹. We only used wing components whose amplitude is of at least 5σ significance, where σ_rms is the rms noise in the line free region.

3.2.2. BPT diagram

To assess the role of radiation pressure in driving outflows in our sample, we also used the ionisation state of the sources. For this we used the BPT diagnostic diagram (Baldwin et al. 1981; Kewley et al. 2006). The BPT diagram uses a combination of ionised line flux ratios to identify whether the main source of ionisation is the light from the AGN or star formation. Using the emission line fluxes obtained in Sect. 3.2.1, we can plot our sources on the BPT diagram. We used the Hα, Hβ, [N II], and [O III] line flux ratios and require that at least one of [N II] and Hα, and one of [O III] and Hβ are detected. This is needed in order to place the source on the BPT plot, even if its only an upper or lower limit. This gives us a total of 89 radio AGN, plotted in Fig. 3.

Fig. 3. BPT diagram of the sources, with sources marked according to their radio spectral shape. The solid curve is the maximum starburst line from Kewley et al. (2001). The dashed curve is the semi-empirical line from Kauffmann et al. (2003). The solid line separating high and low ionisation sources is the empirical relation from Schawinski et al. (2007).

We find that most of our sources lie in the AGN ionisation region (Seyfert or LINER), which is in agreement with the results of M17, who also found that most of the sources lie above the MIR-radio relation, suggesting that the bulk of their radio emission comes from the central AGN rather than star formation. This confirms that the contribution from star formation to the radio emission and ionisation of the ISM, is insignificant in our sample. We then grouped the sources into high and low ionisation states, using the empirical relation from Schawinski et al. (2007). Sources lying in the Seyfert region in Fig. 3 are classified as high ionisation sources and those in the LINER region are classified as low ionisation sources. We classified 23 sources as high ionisation and 58 sources as low ionisation (see Table 1). Their [O III] spectra is discussed in Sect. 4.2.3.

3.2.3. Individual outflow detections

During our modelling analysis in Sect. 3.2.1, 21 sources required a multi-component model to fit the [O III] spectra. Out of these, eighteen sources required double and three sources required triple kinematic components. As mentioned, multi-component fits hint towards the presence of kinematically disturbed and possibly outflowing gas. Using our aforementioned criteria, out of the 21 sources that had a multi-component fit, we identify 18 sources to have an outflowing component. This gives an outflow detection rate of ∼14% for our sample. We discuss these sources in Sect. 4.1.

Since we want to investigate the presence of outflows in sub-groups of sources and the number of individual outflow detections in our sample is relatively low, and we decided to use a stacking analysis for the rest of the paper. This allows us to explore, on average, the link between the gas kinematics and source properties. The stacking analysis is explained in the next section, and the results are shown in Sect. 4.2.

3.2.4. Stacking procedure

As mentioned above, we detect individual outflows in a small number (∼14%) of sources. Studies of ionised gas outflows often rely on individual detections of the outflowing gas component, and this can bias them towards the strongest outflows. But since we aim to study the presence and dependence of outflows on the properties of radio AGN, it is important to include weak outflows that might not be detected individually with the sensitivity of the SDSS spectra. To do this, we stack the stellar continuum subtracted SDSS spectra, in groups of sources based on their radio and optical properties, to follow the evolution of the ionised gas kinematics.

We calculated the stacked spectra using a weighted average of the individual source spectra, defined as:

$$\begin{array}{c}\hfill S(\lambda )=\frac{\sum _i\frac{1}{{\sigma }^2(i)}S(i,\lambda )}{\sum _i\frac{1}{{\sigma }^2(i)}},\end{array}$$(2)

where S(λ) is the stacked spectrum, S(i, λ) is the individual source spectra, and σ(i) is the noise in a line free region of the source spectra. In the stacked spectrum, the noise is expected to go down by 1/$\sqrt{N}$, where N is the number of stacked sources. In our sample, stacking all the 129 radio AGN with an average noise of 1.33 × 10⁻¹⁷ erg s⁻¹ cm⁻² Å⁻¹, gives a noise ∼0.19  ×  10⁻¹⁷ erg s⁻¹ cm⁻² Å⁻¹ in the stacked spectrum. This is in agreement with the expected theoretical decrease by about a factor of 11. The errors on the stacked spectra were estimated using the bootstrapping method. We randomly selected 5000 samples from each group with replacement and stacked them. Using the 5000 stacked spectra, we estimated the errors at each point by taking the standard deviation of the stacked values at that point. The stacked spectra for various groups of sources are shown and discussed in Sect. 4.2.

3.2.5. Stacked emission line fitting procedure

For the analysis of the stacked spectra, we have used a Bayesian inference model fitting routine to assess the quality of the fit that best describes the profile. This is a probabilistic approach, and different from our method in Sect. 3.2.1. It is better suited for assessing the average properties of a population of sources than the least squares approach. The method uses the Bayes theorem to obtain posterior probability density functions for the model parameters. In this routine, we estimate a likelihood function ℒ(θ)≡P(D|θ, M), which gives the probability of observing the data given the model parameters. We then used the Markov chain Monte Carlo methods to sample the posterior probability density functions of the model parameters such that this likelihood function was maximised. To implement this, we used the emcee package (Foreman-Mackey et al. 2013). For the [O III] doublet, our model was made up of Gaussians with fixed width, mean, and relative amplitudes of the doublet, as described in Sect. 3.2.1. To account for a non-zero baseline in the spectra, we also added a first order polynomial of the form a_px + b_p to the model.

To assess the quality of the fit we used the Bayesian information criterion (BIC; Schwarz 1978) which is estimated using a combination of the number of free parameters in the model, the number of data points and the likelihood value. This criteria is designed to give a lower value for a better fit, while punishing for adding free parameters. Therefore, the model with the least BIC is preferred. We define ΔBIC = BIC₂ − BIC₁, where BIC₂ and BIC₁ are the values for the models with higher and lower number of parameters, respectively. The general rule of thumb is this: a −2 ≤ ΔBIC < 0 means there is weak evidence that adding more parameters improves the model fit, −6 ≤ ΔBIC < −2 means there is moderate evidence that adding more parameters improves the fit, and −10 ≤ ΔBIC < −6 means there is strong evidence, and ΔBIC < −10 means there is very strong evidence that adding more parameters improves the fit (Jeffreys 1961). For our analysis, we take a conservative approach and only consider a model with more parameters to improve the fit if ΔBIC ≤ −6.

For the fitting of the stacked line presented below in Sect. 4.2, we started by fitting each line in the doublet with a single Gaussian, and performing this fit over 1000 times. The best fit parameters and the errors were then estimated from the average of the 1000 runs. These values were used to estimate the BIC. We then added another Gaussian component to the model and repeated the procedure. The best fit model was then picked depending on the change in BIC. We found that each line can be fit with one to two Gaussians. Fitting three Gaussians always worsened the fit quality, which is consistent with the fact that only three individual [O III] detections required triple components. We apply this to the analysis of the stacked spectra of the various sub-groups as described in Sect. 4.2. The results are summarised in Table 3.2.5.

3.3. Eddington ratios

Studies of ionised gas kinematics in AGN have found an increase in velocity dispersion and outflow detection rates with the Eddington ratio (for example, Mullaney et al. 2013; Woo et al. 2016; Singha et al. 2021). They conclude that this positive correlation is evidence for gas kinematics and outflows being driven by radiation pressure. To investigate this link, we calculate the Eddington ratio for our sources. The Eddington ratio can be calculated as the ratio of the total radiative output to the maximum luminosity that can be emitted by spherical accretion of gas on to the central SMBH, that is the Eddington luminosity (Best & Heckman 2012). The Eddington ratio is then defined as λ_Edd = L_rad/L_Edd.

For the radiative luminosity of our sources, we have used the relation of Heckman et al. (2004), L_rad = 3500 × L_[OIII], where L_[OIII] is the observed [O III] luminosity in the SDSS spectra. The scatter in this relation is ≈0.4 dex. In order to determine the Eddington ratios for all sources, we have used the upper limit in luminosity for [O III] non-detections. There is a caveat to using [O III] luminosity as a proxy for the radiative luminosity for our sample, that we mentioned in Sect. 3.2.1. The Eddington luminosities for AGNs can be determined using the relation L_Edd = 1.3 × 10³¹ M_BH/M_⊙ W, where M_BH is the mass of the central SMBH in units of solar mass. To obtain the SMBH mass, we used the M − σ_⋆ relation from McConnell & Ma (2013), log(M_BH/M_⊙) = 8.32 + 5.64 log(σ_⋆/200 km s⁻¹), where σ_⋆ is the stellar velocity dispersion. This relation has a scatter of ≈0.4 dex. We used the stellar velocity dispersion measurements from SDSS. Velocity dispersion values ≤100 km s⁻¹ and ≥850 km s⁻¹ were unreliable and have been removed. Finally, we have estimated the Eddington ratios for 121 sources in our sample. The distribution is shown in Fig. 4. We split the sample about the median Eddington ratio to assess its role in feedback, and show the results in Sect. 4.2.4.

Fig. 4. Eddington ratio distribution of the 121 radio AGN with reliable estimates of stellar velocity dispersion, shown in grey. The distributions for peaked and non-peaked radio AGN are shown in blue and orange respectively. Both groups have a similar range of Eddington ratios.

4. Results

The goal of this study is to connect the properties of the AGN with the one of the ionised gas. We not only consider the radio and optical properties, but also focus on the group we have identified as young radio AGN. As mentioned in the introduction, some trends with these properties have been found but, with all the parameters derived for our sample and discussed in Sect. 3, we are now able to assess their role, and use stacking a stacking analysis to disentangle the role of these properties. We start by discussing the individual outflow detections we have obtained, followed by the stacking results.

4.1. [O III] outflows

As mentioned in Sect. 3.2.3, 18 sources have an outflowing component in their [O III] spectra according to our criteria. These components are shown in Fig. 5 which plots their velocity difference from the core component vs their velocity dispersion. We find that the peaked sources have broader outflowing components, typically with a higher velocity dispersion than their non-peaked counterparts. This suggests that the kinematics of the ionised gas is linked to the radio spectral shape. Physically, this hints towards more disturbed gas at the younger stages of radio AGN life cycle. Although we see hint of a trend with the radio spectral shape in Fig. 5, the low number of individual outflow detections (∼14%) makes it hard to draw a link with the AGN life cycle. Therefore, to explore this trend for the entire sample, we stacked the optical spectra for all the sources, as mentioned in Sect. 3.2.4. Below we show the results from the stacking analysis.

Fig. 5. Velocity (with respect to systemic velocity) vs dispersion of the wing components identified as outflows. Sources are coloured according to their radio spectral shape. Horizontal lines show the average stellar velocity dispersion of the sample and its three sigma upper limit. Peaked spectrum sources show broader outflows with larger velocities, compared to non-peaked sources.

4.2. Stacking results

4.2.1. Radio spectral shape classification

The stacked [O III] spectra and model fits for the peaked and non-peaked sources are shown in Fig. 6. The peaked sources clearly show a broader spectra on the blueshifted side than the non-peaked sources. The FWHM_avg distributions of the two groups also show a higher velocity tail for the peaked sources. A two sample KS test gives a statistic of 0.34 and a p value of 0.04, which suggests that the two distributions are significantly different at 95% confidence level. This confirms our result in the previous section, that the [O III] kinematics is linked to the radio spectral shape. To check whether this difference in stacked spectra is due to the individual outflow detections, we re-performed the stacking analysis after removing the outflow detections. We found that the stacked spectra for peaked sources still had a broad component on the blueshifted side, compared to the non-peaked sources. However, the difference is not as strong as before, due to the decrease in the number of sources stacked. A two sample KS test for these groups’ distribution of [O III] average FWHM gives a statistic of 0.34 and a p value of 0.09, which suggests that the two distributions are significantly different at 90% confidence level. This confirms the fact that our result for peaked and non-peaked sources is true on average for sources, and is not affected by individual outflow detections.

Fig. 6. Stacking results for peaked and non-peaked sources. Top left: average FWHM distributions of the [O III] detections in the peaked and non-peaked groups. Bottom left: stacked [O III] spectra of peaked and non-peaked sources, including non-detections. The shaded regions show the errors in the stacked spectra. The peaked sources show a significantly broad spectra on the blueshifted side, than the non-peaked sources. Bottom right: model fit to the peaked spectra, with the core and wing components marked. Top right: residuals after subtracting the model from the peaked spectra. The shaded region marks the standard deviation of the residuals.

Our modelling analysis shows that a two component model provides the best fit for the peaked sources, whereas a single component model gives the best fit for the non-peaked sources (see Table 3.2.5 and Fig. 6). Using our criteria, it is clear that the peaked sources have a blueshifted outflow, whereas non-peaked sources do not. The outflowing component has a peak velocity of −237 ± 102 km s⁻¹, and is very turbulent, with a FWHM of 1330 ± 418 km s⁻¹ (after deconvolving velocity resolution). The outflow component also contains ≈38% of the total [O III] flux. This suggests that on average, the [O III] gas kinematics evolves with the radio jets life cycle, and that the ionised gas is most disturbed when the radio spectra is peaked, that is at the youngest stage of radio jet evolution.

As mentioned in Sect. 3.1.2, peaked sources could also be explained by the “frustration” hypothesis due to a dense ambient medium. To estimate the electron densities of the NLR region, we used the stacked [S II] doublet and the relation from Sanders et al. (2016), to obtain electron density of n_e = 237 ± 52 cm⁻³ for peaked and n_e = 243 ± 55 cm⁻³ for non-peaked sources. We note that this relation between [S II] ratio and n_e saturates at a ratio of ∼0.44 and cannot be used for gas with n_e > 10⁵ cm⁻³. However, this should not affect our results because the electron density we obtain is significantly lower than this limit. We then estimate ionised gas masses of M_ion = 4.6 × 10⁵ M_⊙ for peaked, and M_ion = 6.4 × 10⁵ M_⊙ for non-peaked sources (Revalski et al. 2022). The electron densities and ionised gas masses are similar for the two groups. This shows that peaked sources typically do not have a denser ambient medium that non-peaked sources in our sample. We also calculate that for the radio spectra to have a peak in the GHz frequencies (≳1.4 GHz), if surrounded by a medium of such electron densities, the path length would have to be ≳10 kpc. It is highly unlikely that the NLR clouds exist out to such large radial distances for majority of our sources. Therefore, we are confident that our peaked sources, on average, are not “frustrated” or absorbed due to FFA, and are genuine young sources (which nevertheless appear to be interacting with the surrounding medium, likely slowing their expansion for example, Bicknell et al. 2006, 2018; Wagner & Bicknell 2010). Some of our sources are also included in the de Vries et al. (2009) sample, who have used VLBI data to conclude that these are young radio AGN, absorbed due to SSA.

4.2.2. Radio morphology and luminosity

To assess the impact of radio morphology and luminosity on the [O III] gas kinematics, we also stacked the sources in groups of compact and extended morphology, and high and low 1.4 GHz radio luminosity (see Table 1 for details). The stacked spectra are shown in the top panel of Fig. 7. We find no significant difference in the stacked spectra based on groups of radio morphology or 1.4 GHz luminosity. The modelling analysis also shows that the best fit is provided by a single component model (see Table 2 and Fig. B.1). Although it can be seen in Fig. B.1 that there is some emission on the blue side that could be fit with another component, we found the improvement in the fit with a double component model was not significant enough to be considered the best fit model, according to our BIC criteria (ΔBIC ≤ −6). This does not mean that there are no outflows in either group. The main result is the fact that, based on the lack of a strong difference in the stacked profiles, combined with the model fit results, we do not find a strong dependence of the presence of outflows on these properties. A two sample KS test of the average FWHM distributions for compact and extended sources gives a statistic of 0.22 and a p value of 0.21. Similarly for radio luminosity, we obtain a KS statistic of 0.16 and a p value of 0.59. Thus we find no significant differences between the average FWHM of the [O III] profiles in these groups.

Fig. 7. Stacked [O III] spectra and FWHMavg distributions for groups based on (a) radio morphology, (b) central 1.4 GHz luminosity, (c) BPT diagram and (d) optical luminosity. The top panel in each sub figure shows the distribution of the [O III] detections in the groups. Bottom panel in the sub figures shows the stacked [O III] spectra (including non-detections). The fluxes have been normalised to the maximum value for easier comparison of the profile shapes. We only see a significant difference in the stacked spectra based on BPT classification (c).

To check the robustness of this result, we performed the same stacking analysis by changing our classification criteria of source morphology and using either only the criteria from Shimwell et al. (2022) or G15. This did not affect our results. The lack of significant difference in the stacked spectra of radio AGN on the basis of radio morphology and luminosity shows that these are not the dominant properties affecting [O III] gas kinematics in our sample, and the link between the gas kinematics and radio spectral shape.

4.2.3. Ionisation ratio and [O III] luminosity

Our final goal is to be able to disentangle the effects of radio properties from radiation in our sample. Therefore, to assess the role of radiation pressure in driving the [O III] gas kinematics, we stacked the spectra in groups of two classifications – the ionisation state, using the BPT diagram (Sect. 3.2.2) and the observed [O III] luminosity. The stacked spectra are shown in the bottom panel of Fig. 7. Visual inspection suggests a difference in the width of the profiles as can be seen in the figure. Our modelling analysis (Fig. B.1) also indicates a difference between high and low ionisation state groups as can be seen in Table 3.2.5. The core component of the low ionisation group has a broader FWHM of 645 ± 45 km s⁻¹ (after deconvolution of the instrument resolution) than the high ionisation group FWHM of 425 ± 33 km s⁻¹. However, the low ionisation group also has a larger average stellar velocity dispersion (σ_⋆ = 230 ± 8 km s⁻¹) than the high ionisation group (σ_⋆ = 170 ± 8 km s⁻¹). The [O III] velocity dispersion of the core component has been known to correlate with the stellar velocity dispersion for Seyferts (Greene & Ho 2005; Komossa et al. 2008), and recently, such a trend has also been found for type II AGN Woo et al. (2016, 2017). Therefore, this difference in the FWHM of the core component between the low and high ionisation groups can be related to the difference in the average stellar velocity dispersion. Interestingly, the high ionisation state profile is best fit with a double component model, whereas the low ionisation state profile is best fit with a single component model. The wing component in high ionisation group has a FWHM of 1208 ± 143 km s⁻¹ (after deconvolving velocity resolution), and is identified as an outflow according to our criteria. A two sample KS test of the average FWHM distributions between the two groups also gives a statistic of 0.45 and a p value of 0.002, suggesting a difference between the [O III] profiles at > 99% confidence level. These results suggest a link between the ionisation state and the [O III] kinematics.

However, for the groups in [O III] luminosity, we find no significant difference in the stacked spectra (Sect. 3.2.4). Both are best fit with a single component model according to our criteria, although there is some emission on the blue side in Fig. B.1 that could be fit with another component. However the improvement is not significant enough to be considered as best fit model. Therefore, on average, we find no dependence of the presence of an outflow on [O III] luminosity. Re-performing the analysis using only [O III] detections also did not affect our results. A two sample KS test of the average FWHM distributions also gave a statistic of 0.18 and a p value of 0.44, thus showing no significant difference between the two groups. These results from the ionisation state and [O III] luminosity groups appear contradictory in determining the role of radiation pressure in our sample. This could mean that radiation pressure does not play a definitive role in our sample of radio AGN, and we discuss it further in Sect. 5.2.

4.2.4. Eddington ratio

As discussed in Sect. 3.3, we split our sample into two groups of Eddington ratios to investigate their link with [O III] kinematics. The stacked spectra are shown in Fig. 8 and show a good overlap with each other. The modelling analysis (see Table 2 and Fig. B.1) shows that both spectra are best fit with a single component model. A two sample KS test of the average FWHM distributions gives a statistic of 0.34 and a p value of 0.01, suggesting a difference between the [O III] profiles at > 99% confidence level. However this difference is likely due to the larger average stellar velocity dispersion of the lower Eddington ratio group, by ≈43 km s⁻¹, as also discussed in Sect. 4.2.3.

Fig. 8. Stacked results for Eddingtion ratio. Top: average FWHM distributions of the [O III] detections and (bottom) stacked [O III] spectra for the high and low Eddington ratio groups. The fluxes have been normalised to the maximum value.

5. Discussion

Our aim has been to characterise the evolutionary stage for a sample of radio AGN, and study feedback on the ionised gas over the life cycle of these AGN. By using a stacking analysis in groups based on the optical and radio properties of the sources, we have tried to disentangle their role in driving the ionised gas outflows. In this section, we discuss our results and the link between the life cycle and feedback on the ionised gas. We then estimate the outflow energetics, and present our conclusions.

5.1. Link to radio AGN life cycle

The main result of our paper is the link between radio spectral shape and the presence of outflow, shown in Sect. 4.2.1. We find that stacking the peaked sources together shows a more prominent blueshifted component in the [O III] profile, whereas stacking non-peaked sources does not (see Fig. 6 and Table 3.2.5). We interpret this as outflowing gas in the peaked sources. This outflow is very broad, with a deconvolved FWHM = 1330 ± 418 km s⁻¹. In Sect. 4.1 and Fig. 5, we also find that individually detected outflows in peaked spectrum sources have higher velocities and widths than non-peaked sources. These results lead us to conclude that the feedback on the ionised gas is linked to the radio spectral shape of the sources. Since the spectral shape is a good indicator of the evolutionary stage, feedback is therefore linked to the life cycle of the AGN. This shows that impact of the jets evolves over the AGN life cycle. The ionised gas kinematics is most disturbed when the radio jets are young (peaked sources), and it becomes less extreme as the jets evolve (non-peaked sources).

The outflow velocity of −237 ± 102 km s⁻¹ in peaked sources is lower than those found in powerful radio galaxies. For example, for a sample of 3C radio AGN, Speranza et al. (2021) found outflowing gas with blueshifted velocities ranging from ∼300 km s⁻¹ to ∼900 km s⁻¹. Similarly, Santoro et al. (2020) found blueshifted outflow velocities typically from ∼500 km s⁻¹ to ∼1000 km s⁻¹. Although not radio AGN, but similar values have also been found in ionised gas for optically selected type 2 AGN by Bae & Woo (2016), covering the same L_[O III] range as our sample. Our lower outflow velocity for peaked sources tells us that on average, the outflowing gas in radio AGN is not as fast as in the most powerful sources. However, the broad width of the wing component, FWHM = 1330 ± 418 km s⁻¹, also suggests that the peaked sources either have a large variety in outflow velocities, or very turbulent outflows.

Our results support the picture where in the early stage of evolution, the radio jets interact with the surrounding medium, piling up energy until eventually they break out, pushing the gas to outflow velocities. As these jets grow, they likely clear a channel of gas in the ambient medium. At the later stages of evolution, there is not a lot of gas in the path of the jet to interact with in the central few kpc region, and push to high velocities, and therefore the outflows detected are not the strongest. In the last few years, hydro dynamic simulations have also suggested a similar scenario, where feedback is most prominent when the jets are newly born and they are strongly interacting with the ambient medium, and weakens as the jets grow beyond the medium (for example, Mukherjee et al. 2016, 2018; Meenakshi et al. 2022).

5.2. Relevance of other properties

Our mean stacking analysis (see Sect. 4.2.2 and Table 2) shows no dependence of the width of the stacked [O III] profile on radio morphology and luminosity. This could be due to radio morphology and luminosity not being the dominant factor in deciding the impact of the AGN on the ionised gas. This also tells us that these properties do not affect our result on the link with the radio spectral shape. As mentioned in the introduction, previous studies have found that compact sources have the most disturbed gas kinematics suggesting a link between the radio morphology and ionised gas kinematics in AGN. Since the sources are expected to be compact in their early stages of evolution, this provides a link to the AGN life cycle. However, we find that although a large fraction (∼68%) of the peaked (“young”) sources in our sample are compact, a significant fraction (∼32%) of them are also extended. In some cases, extended sources with a young core could represent multiple epochs of activity. Therefore, we find little evidence for a relation between sources being young and morphologically compact at low frequencies, in our sample. It is worth noting that the low frequency LoTSS maps with LOFAR that we use, are highly sensitive to diffuse emission around these sources, and are therefore well suited for characterising source morphology into compact and extended.

The outflow velocity for peaked sources in our sample is lower than what has been found for more powerful radio AGN, such as the 3C and 2 Jy sources (see Sect. 5.1). This would suggest a relation between the outflow velocity and radio luminosity. However, we do not detect this relation when splitting the group of peaked sources in two groups based on the 1.4 GHz luminosity. The low number of sources in each group also prevent from obtaining a conclusive result. This will require expanding the sample to higher luminosities, which will allow us to investigate this relationship in more detail. We plan to conduct such a study in the future.

The role of radiation pressure in driving feedback in our sample is not entirely clear. In Sect. 4.2.3, we find that higher ionisation state sources have an outflow component (FWHM = 1208 ± 143 km s⁻¹) in their stacked [O III] profile. The high ionisation state group also has a larger average L_[O III] by a factor of ∼10 (see Table 1). This suggests that radiation pressure also drives the outflow in the optically strongest radio AGN in our sample. However, the lack of difference in the widths of the stacked profiles in [O III] luminosity groups contradicts this role of radiation. But we cannot conclusively exclude some role of radiation pressure in our sample. In order to understand the significance of radiation pressure in depth, we need to expand the sample to a higher [O III] luminosity, Eddington ratio, for example. We plan to explore this in a future study.

Despite the presence of a broad outflow component in the high ionisation stacked spectrum, we do not expect our result on the link between feedback and radio spectral shape to be driven by radiation pressure. This is because majority of the peaked sources (∼60%) lie in the low ionisation region on the BPT diagram (see Fig. 3). Furthermore, a larger fraction of non-peaked sources (∼41%) lie in the high ionisation region than peaked sources (∼28%). Peaked sources also have a lower average [O III] luminosity than non-peaked sources, by 0.3 dex (see Table 1). Therefore, the link between [O III] profile and the radio spectral shape is due to the properties of the radio jet, and not dominated by radiation pressure.

In our analysis of stacked spectra in groups of Eddington ratio (Sect. 4.2.4) we found no significant difference in the widths of the stacked profiles. This tells us that, on average, the feedback on [O III] gas in our sample is not driven by the radiation pressure due to accretion. This is contrary to some results found before. For example, Mullaney et al. (2013) used a stacking analysis and found a positive link between the FWHM of the [O III] profile and Eddington ratio for a sample of optically selected AGN. But it is worth pointing out that their sample, comprised of type 1 AGN, typically has larger values of λ_Edd (10^−2.5–10⁰) and L_[O III] (10⁴⁰ − 10⁴⁴ erg s⁻¹). Their estimation of Eddington ratio is also different from our approach. They used the width of the Hα line to estimate L_Edd and its luminosity to estimate L_rad. More recently, Singha et al. (2021) found that presence of outflows increases with Eddington ratio for a sample of low excitation radio galaxies (LERGs). Although LERGs are comparable to the sources in our sample, they also used the parent sample of Mullaney et al. (2013), and therefore the same method for estimating the Eddington ratio. Lastly, we point out that the Eddington ratios of our sources also does not affect our result on feedback and radio spectral shape, since the distribution of this ratio is similar for peaked and non-peaked sources, as shown in Fig. 4.

Finally, studying all the above mentioned properties helps us conclude that the link between the radio spectral shape and feedback on the gas is a result of varying impact of the radio jets on the surrounding medium during their life cycle, and not driven by other radio and optical properties. The lack of an outflow in the stacked profile of non-peaked sources also suggests that extreme feedback, such as an outflow, could also be short lived (≲1 Myr).

5.3. Outflow energetics

To asses the effect of the outflow on their host galaxies in peaked sources, we need to estimate its mass and kinetic power. Since in the stacked spectrum, we only detect a blueshifted outflow, the outflow morphology is likely bi-conical on average, where the flux from the redshifted cone is attenuated by the dust in the central few kpc region. For a bi-conical outflow model, assuming constant density clouds with an electron temperature T = 10⁴ K and case B recombination, the mass outflow rate (Ṁ_out) can be estimated using the following equation from Husemann et al. (2016),

$$\begin{array}{c}\hfill \frac{{\dot{M}}_{\mathrm{out}}}{3M_{\odot }{\mathrm{yr}}^{-1}}=\left(\frac{100{\mathrm{cm}}^{-3}}{n_{\mathrm{e}}}\right)\left(\frac{L_{\mathrm{H}\beta }}{{10}^{41}{\mathrm{erg}\mathrm{s}}^{-1}}\right)\left(\frac{v_{\mathrm{out}}}{100{\mathrm{km}\mathrm{s}}^{-1}}\right)\left(\frac{\mathrm{kpc}}{D}\right),\end{array}$$(3)

where n_e is the electron density of the gas, L_Hβ is the Hβ luminosity of the wing component, v_out is the velocity of the outflowing gas, and D is the size of the outflowing region. Since we cannot spatially resolve the outflow region, we assume the radial size to be equal to half of the SDSS fibre size, which is 1.5″. At the average redshift of peaked sources (z = 0.1), this corresponds to D = 2.8 kpc. L_Hβ for the wing component is hard to determine since the signal-to-noise of the Hβ line is much lower than [O III]. Therefore, similar to Singha et al. (2021), we assume that the wing components of the Hβ and [O III] line have the same flux ratio as the core components, and estimate L_Hβ by scaling the L_[O III] of the wing component by the flux ratio of the core component. We estimate an L_Hβ, wing of 7.6 × 10³⁹ erg s⁻¹ for the wing component.

For the electron density n_e of the gas, we used the emission line ratio of the [S II]λλ6716,6731 Å doublet. Ideally, we should have only used the emission line ratios of the wing components of the [S II] doublet, because it would be representative of the outflowing gas. However, it is again extremely hard to do so due to the low signal-to-noise and blending of the doublet. In our stacked spectra for the peaked sources, we do not detect a wing component in the [S II] doublet. We therefore used the core component to estimate an electron density of 237 ± 52 cm⁻³, as discussed in Sect. 4.2.1¹.

We estimate the values of v_out using two approaches. The first approach is to assume that velocity of the outflowing gas is the centroid of the wing component, and the broadening is due to the velocity dispersion in the outflowing gas. So we set v_out = v_wing − v_core and v_turb = FWHM_wing (deconvolved). This is the most commonly used approach, but it does not account for the line-of-sight projection effects on the observed profile. Therefore, we also used a second approach from Rose et al. (2018), which assumes that the broadening of the wing component is entirely due to the different projections of the outflowing gas velocities and not the intrinsic velocity dispersion. In this approach, we set v_out = v_max, where v_max is the maximum velocity that the outflowing gas can reach, and v_turb = 0. Following Rose et al. (2018), we calculated the maximum velocity as the velocity at which the cumulative flux of the [O III] wing component (integrated from low to high velocities) is 5% of the total flux in wing component. We calculate v_max = 1184 km s⁻¹ from the stacked spectrum.

Using the first approach, we estimate a mass outflow rate of 0.09 ± 0.04 M_⊙ yr⁻¹, and with the second approach, we estimate a rate of 0.41 ± 0.13 M_⊙ yr⁻¹. We would like to note that given the uncertainties in the assumption of the outflow size and density, the true uncertainties in the mass outflow rate are likely to be larger than what we estimate. These mass outflow rates are comparable to what has been found for studies of other radio AGN. Santoro et al. (2020) studied the [O III] outflows for a sample of nine compact radio AGN, and found outflow rates varying from 0.4 to 20 M_⊙ yr⁻¹. For a powerful GPS radio source PKS 1345+12, Holt et al. (2011) found a total mass outflow rate of 7.7 M_⊙ yr⁻¹. Recently, Speranza et al. (2021) studied spatially resolved [O III] outflows for a sample of 3C radio AGN, and found outflow rates varying from ∼1−30 M_⊙ yr⁻¹. It is important to remember that all of these studies involved significantly more powerful radio AGN, and studied individual detections of outflows, which is only possible for the more massive outflows. Our stacking analysis probes also weaker outflows that are not detected individually.

Next, we estimate the kinetic power using the method from Holt et al. (2006), which includes both radial and turbulent components of the outflowing gas. Similar to their approach, we assume that the large width of the outflow component reflects the turbulent motion in the gas. We then used the formula

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{kin}}=\frac{1}{2}{\dot{M}}_{\mathrm{out}}(v_{\mathrm{out}}^2+\frac{v_{\mathrm{turb}}^2}{5.5}),\end{array}$$(4)

where v_turb is the deconvolved FWHM of the outflow component. The AGN feedback efficiency (F) can then be estimated as

$$\begin{array}{c}\hfill F=\frac{{\dot{E}}_{\mathrm{kin}}}{L_{\mathrm{mech}}},\end{array}$$(5)

where L_mech is the mechanical luminosity of the AGN. We have estimated the mechanical luminosity using the relation of Cavagnolo et al. (2010), L_mech  =  7.3  ×  10³⁶(L_1.4 GHz/ 10²⁴ W Hz⁻¹)^0.7 W. This relation was determined by comparing the radio power and the cavity power, for evolved radio jets in giant ellipticals and clusters, that made cavities in the surrounding hot, gaseous haloes as they propagated. Cavagnolo et al. (2010) estimated a scatter of ≈0.7 dex around this relation. To estimate the efficiency required by the radiative luminosity, we replaced L_mech in Eq. (5) with L_rad.

We obtain an outflow kinetic power of 1.1 ± 0.4 × 10⁴⁰ erg s⁻¹ with the first and 1.8  ±  0.6  ×  10⁴¹ erg s⁻¹ with the second approach. This gives a jet feedback efficiency of 0.01% with the first and 0.13% with the second approach. Using the radiative luminosity in Eq. (5), we obtain similar values for feedback efficiencies.

The feedback efficiencies for our sample of peaked radio AGN are low. Other studies of warm ionised gas in powerful radio AGN have found similar values of feedback efficiencies. For example, Santoro et al. (2020) estimate efficiencies of 0.001%−0.01% for the 2 Jy sample and Speranza et al. (2021) estimate efficiency of ∼0.2% for some 3C sources. These feedback efficiencies are in agreement with AGN feedback models that estimate ∼0.5% assuming 10% coupling efficiency (Hopkins & Elvis 2010), but are significantly lower than the 5−10% estimated by models that assume 100% coupling efficiency (for example Silk & Rees 1998; Di Matteo et al. 2005). It is important to remember that our estimates are only for the warm ionised gas phase of the outflow. Majority of the outflowing mass is likely to be in the molecular gas phase, where the outflowing mass can be 2−3 orders of magnitudes higher (for example Combes et al. 2013; García-Burillo et al. 2014; Morganti et al. 2015; Fiore et al. 2017; Oosterloo et al. 2019; Murthy et al. 2022). Therefore, it is crucial to study the multiphase gas to accurately characterise these outflows.

5.4. Link to episodic activity

As mentioned in Sect. 3.1.3, we find eight sources that show extended emission only in their LoTSS images (144 MHz) and not at 1.4 and 3 GHz. The lack of extended emission at higher frequencies indicates that the radio spectral index of this emission is steep, possibly due to synchrotron losses. This suggests that the extended emission is old. Out of these sources, seven have a peaked spectrum from 144 to 3000 MHz, suggesting the presence of a young radio AGN in the central region. We used the combination of these two properties to classify these seven sources as candidate restarted radio AGN.

Out of the seven candidates, six have [O III] detections, and we plot the distribution of their FWHM_avg in Fig. 9, along with the non-peaked sources. These candidates seem to have a broader [O III] profile than the non-peaked sources. A two sample KS test gives a statistic of 0.59 and a p value of 0.03. The low p value suggests that the two distributions are significantly different. This tells us that restarted sources could have more disturbed gas kinematics than evolved radio AGN. However, we note that the number of candidates in our sample are too low to conclude this with confidence. Finding concrete evidence for this trend would require detailed radio spectral studies to confirm the restarted nature of activity in the candidates. Nevertheless, this trend is in line with what we expect from our results, since restarted sources have a young radio AGN in the centre. In the picture we describe above, this would mean that after the jets switch off, the gas falls back in to the regions that were cleared, and therefore when the jets restart, there is gas again available for interaction. This suggests a link between feedback on [O III] gas and episodic activity of radio AGN.

Fig. 9. Average FWHM distributions of the [O III] profiles of candidate restarted and the non-peaked radio AGN.

6. Conclusions and future prospects

Using a combination of multi-frequency radio data to constrain the radio spectral shape, and optical spectra to model the ionised gas kinematics and AGN properties, we have linked the life cycle for a sample of 129 radio AGN to their impact on the surrounding ionised gas. The stacking analysis allowed us to disentangle the role of different AGN properties. Our main results are summarised below:

1. We find that over a radio luminosity range of ∼10²³–10²⁶ W Hz⁻¹ and up to z ∼ 0.2, the radio spectral shape is linked to the feedback on the ionised gas. Peaked sources have the most disturbed ionised gas kinematics, with a turbulent outflow of v_out ≈ −237 ± 102 km s⁻¹ and FWHM = 1330 ± 418 km s⁻¹ in their stacked spectrum.

2. Feedback on the ionised gas evolves over the radio AGN life cycle, and radio jets are most effective at disturbing the gas in the ambient medium when they are young. Lack of an outflow in the stacked spectrum of non-peaked sources suggests that these outflows are short lived.

3. The difference in the stacked profiles of peaked and non-peaked sources is not driven by radio properties such as the morphology and luminosity, and optical properties such as the ionisation ratio, radiative luminosity, and Eddington ratio.

4. On average, the mass outflow rate is 0.09 − 0.41 M_⊙ yr⁻¹ and kinetic power is 0.1 − 1.8 × 10⁴¹ erg s⁻¹ for outflows in peaked sources. These values are small, but comparable to those found for other more powerful radio AGN.

5. Although not conclusive, the presence of broader [O III] profiles in candidate restarted AGN in our sample than the non-peaked sources suggests a link between feedback and the episodic activity of radio AGN.

It is worth investigating the impact of this kind of AGN driven feedback on the host galaxy, which, as we mention above, is likely short lived. The small mass outflow rates also tell us that these outflows are likely not an efficient way to quench star formation in the host galaxy. However, jet-ISM interaction is also expected to make the ambient medium turbulent. This turbulence could keep the gas clouds from collapsing to form stars. Therefore, it would be interesting to explore a larger sample of radio AGN and their optical spectra to look for signatures of jet induced turbulence, and how that evolves with the life cycle. This could be the long term feedback on the gas, and might be more efficient at quenching star formation. We would also detect more candidate restarted radio AGN, thus investigating the link of feedback with episodic activity in more detail. We are currently working on a project to expand this sample to higher redshifts and luminosities. Although we use unresolved SDSS spectra, future studies with integral field unit observations would be crucial to study this form of feedback in detail, and evaluate its importance in the galaxy’s evolution.

Acknowledgments

We would like to thank the anonymous referee for the comments, that helped improve this paper. LOFAR data products were provided by the LOFAR Surveys Key Science project (LSKSP; https://lofar-surveys.org/) and were derived from observations with the International LOFAR Telescope (ILT). LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The efforts of the LSKSP have benefited from funding from the European Research Council, NOVA, NWO, CNRS-INSU, the SURF Co-operative, the UK Science and Technology Funding Council and the Jülich Supercomputing Centre. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. CIRADA is funded by a grant from the Canada Foundation for Innovation 2017 Innovation Fund (Project 35999), as well as by the Provinces of Ontario, British Columbia, Alberta, Manitoba and Quebec.

References

Appendix A: Continuum images

Fig. A.1. LoTSS continuum images of eight sources that were classified are compact in their FIRST images, but show clear signs of extended emission at 144 MHz. The contour levels in all images are 3${\sigma }_{\mathrm{RMS}}\times {\sqrt{2}}^n$ where n = 0, 1, 2… and σRMS is the local RMS noise in the images. The negative contours are marked by dashed lines and are at −3σRMS level. All of these sources (except J103653.01+444818.1) have a peaked source in the centre and are classified as candidate restarted sources.

Appendix B: Model fits for stacked spectra

Fig. B.1. Model fits to the stacked [O III] spectra of various groups. Bottom panel shows the data, the best fit model and the individual components. Top panel shows the residuals after subtracting the best fit model from the data and. Model parameters are given in Table 3.2.5.

Fig. B.1. continued.

Fig. B.2. An example of the SDSS spectra and stellar continuum fit for the source J135217.88+312646.4. The blue regions show the emission lines masked during the fitting. The model reproduces the stellar continuum of our data very well.

All Tables

All Figures

Fig. 1. Sample properties. Left: redshift and (middle) peak 1.4 GHz radio luminosity distributions of the 135 sources in the sample (including mergers). The dip in the number of sources in the interval 0.136 < z < 0.175 is due to the exclusion of sources affected by RFI (G14). Right: [O III] luminosity distribution for all the sources. For non-detections, we used the upper limit of the [O III] luminosity. The hatched area shows the distribution for only the [O III] detections.

In the text

Fig. 2. Colour–colour plot of spectral indices for the central region of the 129 radio AGN. x-axis shows ${\alpha }_{1400}^{144}$ derived using LoTSS and FIRST peak flux densities and y-axis shows ${\alpha }_{3000}^{1400}$ derived using FIRST and VLASS peak flux densities. Regions in the plot for peaked, non-peaked and convex spectrum sources are marked with blue, orange and grey colour respectively. The spectral shape of a source in each quadrant is drawn at the corners. The error bar shows the typical error on the spectral indices. The hatched design marks the area which we classify as flat-spectrum region. Sources lying in this region are not included in peaked, non-peaked or convex spectrum groups. [O III] detections are marked with filled circles, and non-detections are marked with hollow circles.

In the text

	Fig. 3. BPT diagram of the sources, with sources marked according to their radio spectral shape. The solid curve is the maximum starburst line from Kewley et al. (2001). The dashed curve is the semi-empirical line from Kauffmann et al. (2003). The solid line separating high and low ionisation sources is the empirical relation from Schawinski et al. (2007).
In the text

	Fig. 4. Eddington ratio distribution of the 121 radio AGN with reliable estimates of stellar velocity dispersion, shown in grey. The distributions for peaked and non-peaked radio AGN are shown in blue and orange respectively. Both groups have a similar range of Eddington ratios.
In the text

	Fig. 5. Velocity (with respect to systemic velocity) vs dispersion of the wing components identified as outflows. Sources are coloured according to their radio spectral shape. Horizontal lines show the average stellar velocity dispersion of the sample and its three sigma upper limit. Peaked spectrum sources show broader outflows with larger velocities, compared to non-peaked sources.
In the text

Fig. 6. Stacking results for peaked and non-peaked sources. Top left: average FWHM distributions of the [O III] detections in the peaked and non-peaked groups. Bottom left: stacked [O III] spectra of peaked and non-peaked sources, including non-detections. The shaded regions show the errors in the stacked spectra. The peaked sources show a significantly broad spectra on the blueshifted side, than the non-peaked sources. Bottom right: model fit to the peaked spectra, with the core and wing components marked. Top right: residuals after subtracting the model from the peaked spectra. The shaded region marks the standard deviation of the residuals.

In the text

Fig. 7. Stacked [O III] spectra and FWHMavg distributions for groups based on (a) radio morphology, (b) central 1.4 GHz luminosity, (c) BPT diagram and (d) optical luminosity. The top panel in each sub figure shows the distribution of the [O III] detections in the groups. Bottom panel in the sub figures shows the stacked [O III] spectra (including non-detections). The fluxes have been normalised to the maximum value for easier comparison of the profile shapes. We only see a significant difference in the stacked spectra based on BPT classification (c).

In the text

	Fig. 8. Stacked results for Eddingtion ratio. Top: average FWHM distributions of the [O III] detections and (bottom) stacked [O III] spectra for the high and low Eddington ratio groups. The fluxes have been normalised to the maximum value.
In the text

	Fig. 9. Average FWHM distributions of the [O III] profiles of candidate restarted and the non-peaked radio AGN.
In the text

Fig. A.1. LoTSS continuum images of eight sources that were classified are compact in their FIRST images, but show clear signs of extended emission at 144 MHz. The contour levels in all images are 3${\sigma }_{\mathrm{RMS}}\times {\sqrt{2}}^n$ where n = 0, 1, 2… and σRMS is the local RMS noise in the images. The negative contours are marked by dashed lines and are at −3σRMS level. All of these sources (except J103653.01+444818.1) have a peaked source in the centre and are classified as candidate restarted sources.

In the text

	Fig. B.1. Model fits to the stacked [O III] spectra of various groups. Bottom panel shows the data, the best fit model and the individual components. Top panel shows the residuals after subtracting the best fit model from the data and. Model parameters are given in Table 3.2.5.
In the text

	Fig. B.1. continued.
In the text

	Fig. B.2. An example of the SDSS spectra and stellar continuum fit for the source J135217.88+312646.4. The blue regions show the emission lines masked during the fitting. The model reproduces the stellar continuum of our data very well.
In the text