© The Authors 2024

1. Introduction

Active galactic nuclei (AGN) are now understood to be linked to the evolution of massive galaxies in the Universe. This link manifests in the form of feedback that can quench star formation, releasing energy that can prevent cooling of the gas in hot halos, and even regulate the growth of the central supermassive black hole itself (for example Silk & Rees 1998; Fabian 2012; McNamara & Nulsen 2012). The mechanical energy released during a phase of activity can shock the host galaxy’s gas, sometimes even sweeping it up in high-velocity outflows. This mechanical energy can be propagated by accretion-driven winds, radiation pressure, collimated jets of plasma, etc. (see Fabian 2012; King & Pounds 2015; Harrison & Ramos Almeida 2024 for a review). Observational evidence of this form of feedback has been found in both radiatively efficient and inefficient AGN host galaxies.

Feedback driven by collimated jets has been widely observed in radio AGN host galaxies. These are systems where radio emission is dominated by the AGN, and cannot be entirely attributed to star formation. Although their energetic output is often dominated by mechanical energy in relativistic plasma jets, they can also release high levels of radiative energy. Over the years, many observational studies have established the presence of jet-driven feedback in radio AGN host galaxies (e.g. Morganti et al. 2021a; Murthy et al. 2022; Nesvadba et al. 2021; Ruffa et al. 2022; Schulz et al. 2021 for some recent results). A limited number of statistical studies have confirmed a positive correlation between the presence of radio emission and disturbed ionised gas kinematics, in particular Mullaney et al. (2013). An interesting aspect of this study is that their AGN sample was selected using optical emission line ratios, and these sources don’t necessarily harbour radio jets. But it highlights the broader role that jets could play in feedback. Recently, similar results have also been found for quasars and broad-line AGN by Escott et al. (2024), who obtained a larger fraction of sources with outflows in radio-detected than non-radio-detected sources. However, it is worth noting that some studies have also found no significant relation between disturbed ionised gas and radio luminosity or morphology (for example Woo et al. 2016; Ayubinia et al. 2023).

Studies of radio AGN samples have also pointed towards a positive correlation between young radio AGN and more disturbed gas kinematics (ionised and neutral phase; e.g. Holt et al. 2008; Geréb et al. 2015; Roche et al. 2016; Molyneux et al. 2019; Murthy et al. 2019; Kukreti et al. 2023), proposing a picture where the impact of jet-driven feedback evolves throughout the active phase. This picture is broadly supported by simulations of jet-ISM interaction, which find a more extreme impact on the clumpy gas medium at early stages of jet growth (Sutherland & Bicknell 2007; Mukherjee et al. 2018).

For radio AGN, observational evidence for the link between different stages of jet growth (young or evolved sources) and feedback, is still limited to powerful systems (> 10²⁶ W Hz⁻¹) or small samples. Although some trends have been found between radio emission and ionised gas kinematics for large samples (e.g. Mullaney et al. 2013; Molyneux et al. 2019), they have been focused on AGN selected using optical emission line ratios, and are not necessarily radio AGN. To understand the relation between jet growth and feedback using radio properties, selecting sources where the radio emission is dominated by jets is crucial. Our approach is to select such a sample, split it into different stages of life cycle and link it to its ionised gas kinematics.

The radio AGN life cycle has been established over the years using their spectral and morphological properties. In early stages of growth, young radio AGN are found to have a peaked radio spectrum and morphologically compact radio emission (for example Callingham et al. 2017, and see O’Dea 1998; O’Dea & Saikia 2021 for reviews). These properties have been observed in the compact steep spectrum (CSS) and gigahertz peaked spectrum (GPS) sources, and the period of this phase is found to be ∼0.001 − 1 Myr from studies based on radio spectral modelling (Murgia et al. 1999; Murgia 2003) and hot spot separation velocities (Polatidis & Conway 2002; Giroletti et al. 2003; Giroletti & Polatidis 2009). The peaked radio spectrum is either due to synchrotron self-absorption (Tingay & de Kool 2003; Snellen et al. 2000; De Vries et al. 2009a; Artyukh et al. 2008) or free-free absorption from an ionised medium (Kameno et al. 2000, 2005; Callingham et al. 2015; Mhaskey et al. 2019; Keim et al. 2019). At later stages of the life cycle (a few tens of Myr), the radio spectrum is observed to be optically thin, due to the peak shifting to lower frequencies and eventually below the observing frequencies of most radio telescopes.

A single phase of activity can last up to ∼10 − 100 Myr (see Morganti 2017 for a review). After the activity stops and there is no more injection of fresh electrons, radiative losses dominate the emission leading to steep spectra and low surface brightness diffuse emission, called the remnant stage of the life cycle (Parma et al. 2007; Saripalli et al. 2012; Brienza et al. 2017). In some cases, remnant phase emission is found along with emission from a young phase at the centre, suggesting that the new phase of activity has started even before the radio plasma from the older phase has faded below detection limits (13−15%; Jurlin et al. 2020). These are called restarted radio AGN (for example, Brienza et al. 2018; Morganti et al. 2021b,c; Kukreti et al. 2022). Therefore, the radio spectral properties change over the life cycle and can be used to characterise a sample of radio AGN into different stages of evolution.

This approach was used in Kukreti et al. (2023), hereafter K23, where we already found some evidence for more disturbed ionised gas in young radio AGN than their evolved counterparts. We characterised the evolutionary stage using their radio spectral shape, and their ionised gas kinematics using the [O III] spectra. The main result of K23 was that on average, peaked radio spectrum sources had a broader [O III] profile than sources without a peak. Since peaked sources represent a younger phase of activity, as discussed above, this suggested a link between the life cycle and feedback on ionised gas.

Although interesting, this study only used a small sample of ∼130 sources up to z ∼ 0.2 covering a radio luminosity range of ∼10²³ − 10²⁵ W Hz⁻¹. The small sample size also made it impossible to control for important AGN parameters like radio and optical luminosity, source sizes, etc. The purpose of the present paper is to investigate the link between radio AGN life cycle and feedback for a sample of 5720 radio AGN up to z ∼ 0.8 that cover a radio luminosity range from ∼10^22.5 − 10²⁸ W Hz⁻¹. The large size and range of luminosities allow us to control for source properties and disentangle their role.

The paper is structured in the following manner: Section 2 describes the approach for sample construction, Appendix A describes the modelling of the stellar continuum and emission line profiles of the sources, Appendix B describes the diagnostics used to select the radio AGN sample and Section 3 describes the radio properties of this radio AGN sample. We then present the results for different source groups in Section 4 discuss them in Section 5 and summarise our results in Section 6. Throughout the paper, we have used the ΛCDM cosmological model, with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3 and Ω_vac = 0.7.

2. Sample construction

This section outlines the approach used to construct the radio AGN sample. The sample was constructed using the Faint Images of the Radio Sky at Twenty-cm survey at 1400 MHz with ∼5.4″ resolution (FIRST; Becker et al. 1995), the LOFAR Two-metre Sky Survey at 144 MHz with 6″ resolution (LoTSS; Shimwell et al. 2017, 2022), the Very Large Array Sky Survey at 3000 MHz with ∼2.5 − 3″ resolution (VLASS; Lacy et al. 2016, 2020; Gordon et al. 2021) and the Sloan Digital Sky Survey (SDSS) which provides spectra with a 3″ size fibre. These surveys were chosen as they provide radio and optical data over similar resolutions, allowing comparison of these properties over similar physical scales.

The parent sample was the 2014 December 17 version of the FIRST survey catalogue (Helfand et al. 2015) which was cross-matched with the SDSS DR17 spectroscopic catalogue (Abdurro’uf et al. 2022) using a search radius of 3″ (approximately half the FIRST beam size). The cross-matching was done down to a peak flux density of 5 mJy in FIRST, which corresponds to L_1.4 GHz ≈ 5 × 10²³ W Hz⁻¹ at z = 0.2. This provided a sample of 29 921 sources.

We then restricted the sample to z = 0.8, after which the [O III]λλ4958, 5007 Å doublet moves outside the SDSS wavelength coverage. We removed sources below z = 0.02, since below this redshift the SDSS fibres probe too small a region in a galaxy to determine reliable parameters. To obtain accurate redshifts, we only selected sources with a fractional error in redshift less than 0.01. This narrowed the sample down to 17 405 sources.

Next, we cross-matched this sample with the LoTSS DR2 (Shimwell et al. 2022) and VLASS Quick Look epoch 2¹ (Lacy et al. 2016) catalogues. We used the recommendations provided in the VLASS catalogue user guide available online, to select sources with reliable detections. Cross-matching was done using TOPCAT (Taylor 2005), with a radius of 6″, which is comparable to the FIRST and LoTSS resolution. This gave us a sample of 6556 sources, all of which have a LoTSS cross-match and out of which 6400 have a VLASS cross-match.

One of the diagnostics used for selecting radio AGN in this paper involves the Wide-field Infrared Survey Explorer (WISE) mid-infrared colours, which trace dust heated by AGN and stars. We added WISE data using the allWISE IPAC release from November 2013 (Cutri et al. 2021). We used three bands of the WISE catalogue – W1 at 3.4 μm, W2 at 4.6 μm and W3 at 12 μm, that have an angular resolution of 6.1 − 6.5″. We applied a signal-to-noise ratio (S/N) cut of 5 for bands W1 and W2, and 3 for band W3 due to its poorer sensitivity. However, we also kept sources with S/N less than 3 in W3 and used them to determine upper limits. We selected only sources with cc_flag = 000 as suggested in the online user manual. To avoid performing k-corrections, we only cross-matched the 2310 sources up to z = 0.3, and obtained a cross-match for 2308 sources.

Radio luminosity at 1.4 GHz is commonly used in literature for radio AGN. To allow for comparison with other studies, we also estimated the total 1.4 GHz luminosity of our sources. The low sensitivity of FIRST to extended emission means that it could underestimate the total luminosity. The ∼5.4″ resolution also means that a large source can be split up into several individual components, further underestimating the total luminosity. We therefore used the combined radio catalogue from Mingo et al. (2016), who performed a detailed and accurate cross-match between FIRST and the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) sources at 1.4 GHz. The lower resolution (∼45″) of NVSS makes it more reliable for estimating total radio luminosities. We used the combined flux densities from this catalogue to estimate the k corrected 1.4 GHz luminosities.

We aim to use the [O III]λλ4958, 5007 Å doublet to trace the ionised gas kinematics, and other emission line ratios to measure the ionisation properties. For this, we used the SDSS spectra. The details of the stellar continuum and emission line modelling, [O III] profile characterisation, and stacking analysis are described in Appendix A. We remove 52 sources from this sample where either the spectra are corrupted at [O III] wavelength or an accurate characterisation of the [O III] profile is not possible. This gives a final sample of 6504 sources. A summary of the number of sources with different emission line detections is shown in Table 1.

Finally, we identify 5720 radio AGN from this final sample, using a combination of diagnostics described in Best et al. (2005), Best & Heckman (2012) and Sabater et al. (2019). These diagnostics help us construct a clean sample of radio AGN, where the AGN and not star formation dominate the radio emission. The procedure is described in Appendix B. Our radio AGN sample spans a wide range in total 1.4 GHz luminosity from ∼10^22.5 − 10²⁷ W Hz⁻¹, and has 1693 [O III] detections with an observed [O III] luminosity range of ∼10³⁹ − 10⁴⁴ erg s⁻¹. This is shown in Fig. 1a. The host galaxies have stellar masses from ∼10¹⁰ to 10^12.4 M_⊙. The host galaxies also span a wide range in g − r colour, containing both red (g − r > 0.7) and blue (g − r < 0.7) objects, with a median g − r = 1.5.

Fig. 1. Properties of the radio AGN sample. (a): L1.4 GHz vs L[O III] for the full radio AGN sample, with [O III] detections and non-detections shown separately. In the case of [O III] non-detections, the contour plot shows the upper limits for L[O III]. The green hatched region marks the luminosity space occupied by the control sample selected from Mullaney et al. (2013), described in Appendix B. The dashed and dotted lines show the radio luminosity expected from the interaction of AGN-driven winds with the ISM from the model of Nims et al. (2015), for 1% and 10% of the shock energy going into relativistic electrons. (b) Distributions of host galaxy properties for the peaked and non-peaked radio AGN in the low redshift (0.02 < z < 0.4) sample. Only the [O III] detections are shown. Vertical lines show the mean values of the distributions. The same plot for the high redshift sample is shown in Fig. C.2.

Radio emission in the so-called radio-loud AGN is attributed to jets, however, the origin is unclear for radio-quiet AGN. For these sources, radio emission can have multiple origins: star formation, innermost accretion disk coronal activity, low power jets and shocks due to interaction between AGN-driven winds and the ISM (Panessa et al. 2019; Zakamska & Greene 2014; Zakamska et al. 2016). Since we aim to study jet-driven feedback down to L_1.4 GHz ∼ 10²³ W Hz⁻¹, it is important to confirm that the radio emission in our sample is dominated by jets even at low luminosities. Using the diagnostics mentioned before, we have selected our sources to have radio emission significantly more than expected from star formation. Further, most of our sources also have radio luminosities > 10²³ W Hz⁻¹, which is more than expected from the corona (Raginski & Laor 2016). In Fig. 1a we also plot the expected radio luminosity from AGN wind shocks, using the fiducial models of Nims et al. (2015). An overwhelming majority of the sources lie much above this model, confirming that the radio emission is dominated by jets in our sample.

Our radio AGN sample spans a wide range in redshift from z = 0.02 to z = 0.8, leading to a large range of 1.4 GHz luminosities. Figure 2 shows L_1.4 GHz as a function of redshift. It can be seen that the typical radio luminosities being probed increase with redshift, as expected for a flux-limited sample. This means we are probing jets of intrinsically very different powers at the high and low redshift end. Therefore, we split the sample into two groups about z = 0.4 (see source numbers in Table 2). This allows a reasonable comparison, while also keeping enough sources in each sub-sample for a statistically significant analysis. The low redshift sample (z < 0.4) probes a wide range of luminosities (∼10^22.5 − 10²⁶ W Hz⁻¹), while still being comparable to the sample used in K23. We discuss the results for the low and high redshift samples separately in Sect. 4.

Fig. 2. Redshift vs total L1.4 GHz for the selected from the sample in Appendix B. Peaked and non-peaked sources are marked with different colours. The sample is split about the dashed line at z = 0.4 into low and high redshift sources.

To determine the role of jets in driving feedback, it is also crucial to compare the [O III] kinematics of our radio AGN with a control sample of optical AGN sources that do not host radio jets. The diagnostics we use allow us to select a clean radio AGN sample, but the 743 sources classified as SF/radio-quiet AGN could still have low luminosity radio AGN. Thus they do not provide a clean control sample. Instead, we select a control sample from the ≈24 000 type 1 and 2 AGN from Mullaney et al. (2013) up to z = 0.4. This sample overlaps with the redshift range of our low redshift sample. These sources were classified as AGN using their optical emission line widths and ratios from SDSS. To select as ‘clean’ a control sample as possible without radio jets, we only use sources with no radio detections in NVSS or FIRST. But at the sensitivity of NVSS and FIRST, a non-detection does not necessarily imply a lack of radio jets in the source. Therefore we restrict the sources to an upper limit of L_1.4 GHz < 10²³ W Hz⁻¹. Since LoTSS has higher sensitivity than FIRST and NVSS, we also removed any source with a LoTSS cross-match within a 6″ radius. This ensures that our control sample has no radio counterpart, and minimises any contamination from even low-luminosity radio jets. We further restrict the sample to 10^39.5 < L_[O III] < 10⁴² erg s⁻¹, 10¹¹ < M_⋆ < 10^11.6 M_⊙, and 180 < σ_⋆ < 300 km s⁻¹ to match it to the low redshift sample. This gives us a control sample of 514 sources. The parameter space occupied by these sources is shown in Fig. 1a.

3. Radio properties

This section discusses the radio spectra and morphology of the sample. We first estimate the spectral indices using a combination of LoTSS, FIRST and VLASS data. We then divide the sources into different spectral shape groups and link them to their evolutionary stage. Finally, we also characterise the source sizes using LoTSS and FIRST data. These properties are then linked to the [O III] kinematics further in the paper.

3.1. Characterising spectral shape

Using the high resolution of LoTSS, FIRST and VLASS we can probe the central ∼6″ region of the radio sources, over spatial scales similar to that covered by the SDSS fibre. This allows a direct comparison of the radio and optical properties. To do that, we estimate the spectral indices using the peak flux densities. For unresolved sources and sources with very weak extended emission, the peak and total flux densities are similar. However, for sources with bright extended emission, peak flux density can be significantly smaller than the total flux density. If the extended emission is optically thin and dominates the total flux density, using the total flux density to estimate the spectral index would lead to a source being classified as optically thin, even if there is a peaked spectrum source at the centre. This motivated our choice to use the peak flux density for estimating spectral indices of the radio AGN, to maximise the number of peaked spectrum sources and detect them even at the centre of large extended sources.

Estimating accurate spectral indices requires reliable flux density values and comparison over the same spatial scales. The flux densities in the VLASS epoch 1 quick look catalogue were found to be underestimated by a factor of ∼0.87 by Gordon et al. (2021). They compared the VLASS flux densities with TGSS at 150 MHz, WENSS at 330 MHz, SUMSS at 840 MHz, and FIRST at 1400 MHz to estimate the systematic flux scale offset. Assuming their corrected flux scale to be the true flux scale, we have re-scaled our VLASS epoch 2 values to match the corrected VLASS epoch 1 flux scale. For this, we first compared the flux densities of epoch 2 to epoch 1 for our sources. The comparison was done only for point sources with a deconvolved major axis size less than 2″ in both epoch data. We found a mean ratio of 0.95 for peak flux densities from epoch 1 to epoch 2, showing that the uncorrected epoch 2 flux densities were larger than the uncorrected epoch 1 values. Finally, to match the corrected scale from Gordon et al. (2021), we divided our epoch 2 values by ≈0.87/0.95 = 0.92. We scaled all our VLASS peak flux density values by this factor. Although this is a large correction factor and a more accurate flux scale correction probably requires a deeper analysis, this correction is sufficient for our purpose. Changes in this factor do not affect our results significantly as we discuss later.

The next step was to match the resolution of VLASS (∼2.5 − 3″) images to FIRST (∼5.4″). We smoothed the VLASS cutouts around our sources to match the FIRST resolution, using the same approach as in K23. We used the task IMSMOOTH in CASA to smooth all VLASS cutouts to the resolution of FIRST. We re-extracted the flux densities from these images using pyBDSF (Mohan & Rafferty 2015). Finally, the spectral indices were estimated using the formula α = log(F_1, peak/F_2, peak)/log(ν₁/ν₂). We plot all the sources on a colour-colour diagram of spectral indices from 144 to 3000 MHz, shown in Fig. 3. Spectral classification of the sources was then done using this diagram.

Fig. 3. Colour-colour plot of spectral indices for the radio AGN sample. The different quadrants are labelled with the spectral shapes of the sources lying in that region. The spectral shapes are also shown above the labels. The red square marks the region where flat-spectrum (FS) sources lie. A 1:1 line is shown in blue. While studying the link of spectral shape with gas kinematics, only HFPS, LFPS and NPS sources are considered.

As discussed in the introduction, the radio spectrum is absorbed during the young phase in the radio AGN life cycle. If the observed spectral shape is fully inverted in our frequency range (${\alpha }_{1400}^{144}>0$, ${\alpha }_{3000}^{1400}>0$), we classify the source as high frequency peaked spectrum (HFPS), except for a source in the FS region, as discussed below. They lie in the top right quadrant of the colour-colour plot in Fig. 3. These sources likely have a peak frequency larger than 3000 MHz. Some HFPS sources may have a peak frequency between our frequency range, however, we expect this group to be dominated by the inverted spectrum sources. If the observed spectral shape has a turnover within our frequency range (${\alpha }_{1400}^{144}>0$, ${\alpha }_{3000}^{1400}<0$), we label the source as low-frequency peaked spectrum (LFPS), shown in the bottom right quadrant of the colour-colour plot. These sources likely have a peak frequency within our frequency range. Again, sources lying in the FS region are not included, as discussed below. In the radio AGN life cycle, HFPS sources represent a younger phase than LFPS sources, since the peak moves to lower frequencies as the source evolves. We first combine these two groups and call them peaked spectrum (PS) sources, to study their [O III] properties. But we also study these two groups separately in more detail. Thanks to the low-frequency point provided by LoTSS, we can detect PS sources down to 144 MHz.

As radio AGN evolve, the peak in their spectrum shifts to lower frequencies, and eventually the spectrum becomes fully optically thin (${\alpha }_{1400}^{144}<0,{\alpha }_{3000}^{1400}<0$). We label these as non-peaked spectrum (NPS) sources, and they lie in the bottom left quadrant of Fig. 3, except for the sources in the FS region. The host galaxy properties of PS and NPS sources are compared in Figs. 1b, 2 and C.2. It can be clearly seen that they have similar L_1.4 GHz, L_[O III], M_⋆, D_n(4000) and σ_⋆ distributions. The σ_⋆ values used throughout the paper have been corrected for the SDSS instrumental resolution (≈70 km s⁻¹) and the resolution of the base spectra used for stellar continuum modelling (≈150 km s⁻¹). The distributions show that there are no intrinsic differences between the host galaxies of PS and NPS sources. A similar conclusion for the optical properties of GPS/CSS and Megahertz-Peaked spectrum sources was also reached by Nascimento et al. (2022).

Sources where free-free emission from the core dominates the radio emission, or where the emission is relativistically beamed due to a small jet viewing angle, can have a flat radio spectrum. These flat spectrum sources are usually defined to have a spectral index α > −0.5 (O’Dea 1998) and are marked with the red box in Fig. 3. We classify sources with $-0.55<{\alpha }_{1400}^{144}<0.05$ and $-0.65<{\alpha }_{3000}^{1400}<0.15$ as flat spectrum (FS) sources. The spectral index and radio luminosities of these sources could affected by relativistic beaming, leading to ‘false’ high luminosity sources. Therefore we do not use them in our analysis of the [O III] kinematics. Some sources also have a convex-shaped spectrum, which is optically thin at lower frequencies and optically thick at higher frequencies. Although their nature is uncertain, such spectra could be a result of multiple epochs of jet activity. We label these as convex spectrum (CS) sources and they lie in the top left quadrant of Fig. 3.

Out of the 5720 radio AGN, flux densities with a signal-to-noise ratio greater than 5 were available for 5686 sources in LoTSS and 5569 sources in VLASS. The spectral shape classifications and number of sources in each group are summarised in Table 2. We note that the flux scale correction of VLASS data could overestimate the flux densities at 3000 MHz, artificially flattening the spectral index between 1400 − 3000 MHz. This could cause some NPS sources to move into the FS and CS region of Fig. 3, and lower the number of NPS sources. However, this will not affect our results since we are still able to select a clean sample of NPS sources.

Although we use only the peak flux densities in the central region to calculate the spectral indices, the physical scale covered by a telescope beam increases with redshift. This leads to a larger contribution from any steep spectrum diffuse emission around the PS source, if present, to the observed peak flux density. This resolution effect means that genuine PS sources could be classified as NPS. Therefore the fraction of PS sources should decrease and the fraction of NPS sources should increase with redshift, as can be seen in Fig. 4. This contamination by PS sources makes it harder to select a clean sample of NPS sources at high redshifts. We discuss this in the context of our results in Sect. 5.1.

Fig. 4. Fraction of PS and NPS sources in redshift (z) bins for the full radio AGN sample. The fraction is the number of PS (or NPS) sources divided by the total number of sources in the bin (of all spectral shapes). The horizontal error bars show widths of the redshift bins used to estimate the fraction and the vertical error bars show uncertainties on the fraction of PS and NPS sources in the bin.

Errors in the spectral indices were estimated using a quadrature combination of the RMS noise and the flux scale errors of the surveys. We used a flux scale error of 10% for LoTSS and 5% for FIRST. Since we have also performed scaling of the VLASS flux scale, we used a conservative value of 10% for VLASS data. The median error for the spectral index from 144 − 1400 MHz was 0.05, and for 1400 − 3000 MHz was 0.15. To test the impact of these errors on the spectral shape classification, we took the following approach. We added a random value Δα_i, to each spectral index α_i, to obtain a spectral index with an offset, given by ${\alpha }_i^{\mathrm{offset}}={\alpha }_i+\mathrm{\Delta }{\alpha }_i$. This Δα_i was drawn randomly from a Gaussian distribution with unit area, zero mean and standard deviation equal to the error in α_i. The spectral indices with offsets are shown in Figs. C.3a and C.3b. We then reclassified all the sources using the colour-colour plot with the new offset spectral indices. All results presented in this paper were then tested with the classifications obtained with these spectral index offsets.

3.2. Radio source sizes

Morphologically compact radio AGN are considered to be young sources in the radio AGN life cycle. Therefore the presence of more disturbed ionised and neutral gas in these sources has been used to suggest a link between the radio AGN life cycle and feedback (e.g. Geréb et al. 2015; Maccagni et al. 2017; Molyneux et al. 2019; Murthy et al. 2019; Santoro et al. 2020). However, in this study, we aim to go further and use the radio spectral properties to characterise young radio AGN. Indeed, it can be seen in Fig. 5 that compact sizes are not always the most reliable indicator for a young source. We find that PS sources span a wide range of sizes in LoTSS, but are almost all smaller than 3″ in FIRST². This difference is due to the significantly higher sensitivity of LoTSS to extended emission than FIRST. This can also be seen in the fact that LoTSS sizes for our sources are consistently greater than FIRST sizes. NPS sources on the other hand span a wide range of sizes in both surveys. Most importantly, this shows that using only a size threshold in a single frequency continuum survey does not guarantee the selection of a clean sample of young radio AGN.

Fig. 5. Deconvolved major axis sizes from LoTSS vs. FIRST for the full radio AGN sample up to z ∼ 0.8. The red, blue and grey points mark the PS, NPS and the rest of the sources in the sample. The horizontal and vertical dashed lines mark the 3″ size. A 1:1 ratio is marked by the dashdotted line.

As mentioned before, studies have found evidence for more disturbed ionised gas in compact sources. However, there is no universal definition of a compact radio AGN. For comparison with the literature, we also explore this link in our sample. Different studies use either a ratio of major and minor axis sizes, total and peak flux densities or a combination of the two to classify a source as compact or extended. Using angular sizes for this can lead to sources of very different physical sizes being classified as compact (or extended) at low and high redshifts. To avoid this, we use the deconvolved major axis sizes from LoTSS as a proxy for the radio sizes of our sources and convert them to physical sizes using the source redshifts. Although the majority of our sample is unresolved (∼75%) in LoTSS according to the criteria of Shimwell et al. (2022), at z < 0.4 these sizes can still give reasonable upper limits. We explore the relation between [O III] kinematics and source sizes in Sect. 4.6 and discuss them in Sect. 5.4. To quantify the link between the two, we use the fraction of disturbed [O III] sources above and below 20 kpc. We choose this limit because GPS and CSS sources, which are young radio AGN, typically have sizes varying between 1 − 20 kpc (e.g. O’Dea 1998; O’Dea & Saikia 2021; Fanti et al. 1990; Dallacasa et al. 2002).

4. Results

We can now investigate the link between the radio spectral shape (a proxy for the evolutionary stage) and [O III] gas kinematics in our sample of 5720 radio AGN up to z ∼ 0.8. As we want to trace the changes in [O III] gas kinematics in different stages of the radio AGN life cycle, we restrict the sample to PS and NPS sources. This gives a final sample size of 2811 radio AGN, out of which 945 have [O III] detections (≈34%). We compare PS and NPS source properties in Figs. 1b and C.2, which shows that they have quite similar host galaxies. We remind the reader that PS sources represent a young phase of activity (∼0.1 − 1 Myr) whereas NPS sources represent more evolved sources (∼10 − 100 Myr). As mentioned in Sect. 2, we divide the sample into two groups about z = 0.4 and discuss them separately (see Table 2). We first discuss the role of radio spectral shape in driving [O III] kinematics, followed by L_1.4 GHz, L_[O III] and source sizes.

4.1. Radio spectral shape

Although our radio AGN host galaxies are selected to harbour radio jets, that does not imply a lack of radiation from the AGN that can also disturb the [O III] gas. Therefore, we first test whether the kinematical disturbances are dominantly driven by jets or radiation pressure. For this, we compare the flux-weighted average [O III] velocity dispersions (${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$) of our radio AGN sample with the control sample of optical AGN selected in Sect. 2, that have a radiatively efficient AGN. This control sample has been selected to not have any radio jets, but is matched in other host galaxy properties to the radio AGN sample. Since the control sample only goes up to z = 0.3, we use the low redshift group for comparison. Here and throughout the paper, ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ has been corrected for the SDSS instrumental resolution. Fig. 6a shows the cumulative distributions of ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ for the two samples (see Appendix A for the definition of ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$). The profiles for radio AGN are consistently broader than optical AGN with similar L_[O III] and stellar mass ranges. Selecting only type 1 or 2 AGN in the control sample, or changing the L_[O III] or M_⋆ range did not affect this result. This shows that [O III] gas is more disturbed in radio AGN host galaxies and highlights the impact of jets on disturbing the surrounding ionised gas.

Fig. 6. Cumulative distributions of the flux-weighted average [O III] velocity dispersions for [O III] detections. The inset plots show the asymmetry distributions for [O III] profiles with multiple kinematic components in every group. (a): Low redshift (0.02 < z < 0.4) radio AGN sample (PS + NPS) and the control sample of optical AGN up to z ∼ 0.4 from Mullaney et al. (2013), described in Appendix B. (b): Low redshift PS, NPS sources and the control sample. (c): High redshift (0.4 < z < 0.8) PS and NPS sources.

Next, we split the low redshift sample into PS and NPS sources, shown in Fig. 6b. There are 217 PS and 206 NPS sources with [O III] detections in the low redshift sample. We find that ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ is consistently larger for PS than NPS sources. A 2-sample KS test shows a significant difference between the two groups (KS statistic = 0.23, p-value = 2.9 × 10⁻⁵). Therefore we can reject the null hypothesis at a > 99% confidence level that the two ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ samples were drawn from the same distribution. The asymmetry parameter distributions, shown in the inset of the plot, are also significantly different (KS statistic = 0.47, p-value = 0.007). PS sources have significantly more blueshifted profiles (median asymmetry = −0.17) in comparison to NPS sources (median asymmetry = −0.05). These trends show that [O III] is more disturbed and blueshifted in PS sources than NPS sources.

To quantify the fraction of sources with ‘disturbed’ [O III] gas, we use a threshold of ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}=350$ km s⁻¹. This limit is the 3σ deviation from the average stellar velocity dispersion of the low redshift sample. Sources in the low redshift sample with ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ greater than this value are classified as kinematically disturbed, throughout the paper. We find that out of the 217 PS sources, 40 are disturbed according to this threshold, which corresponds to a fraction of $18.4_{-2.4}^{+3.0}$%³. On the other hand, out of the 206 NPS sources, only 12 are disturbed, which is a fraction of $5.8_{-1.3}^{+2.1}$%. Therefore we find that [O III] gas in PS sources is ∼3 times more likely to be kinematically disturbed than NPS sources. Using a more conservative threshold of ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}>425$ km s⁻¹ (which corresponds to an FWHM > 1000 km s⁻¹), we obtain a fraction of $8.3_{-1.5}^{+2.3}$% for PS sources and $1.5_{-0.5}^{+1.4}$% for NPS sources. Therefore our results are robust to the chosen threshold value.

We perform a similar comparison of PS and NPS sources in the high redshift sample, shown in Fig. 6c. A 2-sample KS test shows that there is no significant difference in the ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ distributions (KS statistic = 0.06, p-value = 0.86), although PS sources have slightly broader [O III] profiles than NPS sources. The asymmetry distributions too do not show any significant difference (KS statistic = 0.21, p-value = 0.21). Using the same method to define a ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ threshold as for the low redshift sample, we estimate a value of 440 km s⁻¹ for the high redshift sample, above which an [O III] profile is classified as disturbed. This value is used for the high redshift sample throughout the paper. We find that 10 out of the 148 PS sources (6.8${}_{-1.6}^{+2.7}$%) and 11 out of the 374 NPS sources (2.9${}_{-0.7}^{+1.2}$%) are kinematically disturbed. Again, the disturbed proportion is larger in PS sources, by a factor of ≈2. However, the difference between the two groups is only of ≈2σ significance. Therefore, we find only marginal evidence for more disturbed [O III] in PS sources than NPS sources at z > 0.4.

We also find that the results for low and high redshift samples are the same when we use spectral shape classifications after random offsets in spectral indices (as described in Sect. 3.1). This can be seen in the ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ distributions shown in Figs. C.3c and C.3d.

Since both the low and high redshift samples cover a wide range in L_1.4 GHz and L_[O III], there could be underlying trends that drive the differences in the [O III] profiles of PS and NPS sources. To disentangle their role, we first investigate the role of L_1.4 GHz and L_[O III] in driving the [O III] kinematics below. We then study the differences in the PS and NPS sources while controlling for L_1.4 GHz and L_[O III].

4.2. Radio and [O III] luminosities

The total energy output of an AGN is related to L_1.4 GHz and L_[O III]. Therefore, it is important to understand their role in driving the [O III] kinematics in our sample. As shown in Fig. 1, there is a correlation between L_1.4 GHz and L_[O III] of our sources, i.e. sources with high L_1.4 GHz also tend to have large L_[O III]values. Since these two properties are not independent, it is important to control for L_[O III] while studying the impact of varying L_1.4 GHz on [O III] profiles, and vice versa. To achieve this, we constructed a 2D histogram where the sources are binned in L_1.4 GHz and L_[O III] simultaneously, shown in Fig. 7. This allows us to trace the impact of varying one luminosity while controlling for another. The colour map in this plot shows the average ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ value in each bin. The number of sources and fraction of disturbed [O III] sources are also mentioned in each bin.

Fig. 7. Colour map showing the average ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ values for the (a) low redshift (0.02 < z < 0.4) and (b) high redshift (0.4 < z < 0.8 sample in bins of L1.4 GHz and L[O III]. The colour of each bin represents the average velocity dispersion in that bin. The number of sources (the bin size) is mentioned at the top of each bin. Fraction of sources with disturbed [O III] are mentioned below the bin size, only for bins that had at least one disturbed source.

In the plot for the low redshift sample in Fig. 7a, we see a dependence of the [O III] profile widths on L_1.4 GHz and L_[O III]. We find that the highest fraction of disturbed sources and largest ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ values are found at L_1.4 GHz > 10²⁵ W Hz⁻¹ and L_[O III] > 10⁴¹ erg s⁻¹, i.e the high luminosity region. While controlling for L_[O III], the average ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ value and fraction of disturbed sources increases with L_1.4 GHz. The opposite is also true i.e. increasing L_[O III] while controlling for L_1.4 GHz also gives more disturbed [O III]. Although the differences in the fraction of disturbed sources are within 1 − 2σ, we can see a systematic trend with increasing L_1.4 GHz and L_[O III]. To avoid the effect of small number statistics, we only use bins with more than 10 sources for this result. Although not a large bin size, it is a reasonable limit given the low and high redshift sample sizes.

In the high redshift sample shown in Fig. 7b, the trend with luminosities is less significant. Although we see more disturbed profiles in the L_[O III] > 10⁴¹ erg s⁻¹ region, there seems to be no clear trend with L_1.4 GHz or L_[O III]. The fraction of disturbed sources between any two consecutive bins is within 1σ of each other. The sources with the large ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ values also do not cluster in one region of the plot, but they do lie above L_[O III] > 10⁴¹ erg s⁻¹. The fraction of disturbed sources in bins overlapping with the low redshift sample is lower for the high redshift sample, however, this is a result of the larger threshold value for defining an [O III] profile as disturbed in the high redshift sample.

4.3. Radio spectral shape and source luminosities

In the sections above, we found that the radio spectral shape, and the source luminosities (L_1.4 GHz and L_[O III]) are linked to the widths of the [O III] profiles. Given that the PS and NPS sources discussed in Sect. 4.1 cover a wide range of luminosities (see Fig. 1), it is important to disentangle the role of L_1.4 GHz and L_[O III] for these groups. We control for these luminosities and construct similar plots as Fig. 7 for the PS and NPS sources in the low and high redshift samples, shown in Fig. 8.

Fig. 8. Same plot as Fig. 7. The top row shows the plots for (a) peaked spectrum and (b) non-peaked spectrum sources in the low redshift (0.02 < z < 0.4) sample. The bottom row shows the plots for (c) peaked spectrum and (d) non-peaked spectrum sources in the high redshift (0.4 < z < 0.8) sample.

Comparing the PS and NPS sources in the low redshift sample, we see that the [O III] profiles of PS sources are consistently wider at any L_1.4 GHz and L_[O III] bin. For bins with more than 10 sources, PS sources always have a larger fraction of disturbed [O III] than NPS sources, even if we control for L_1.4 GHz and L_[O III]. This shows that the relation with radio spectral shape, discussed in Sect. 4.1, is not driven by source luminosities. The [O III] profiles of these sources are indeed linked to their radio spectral shape. In the colour maps, the difference between the [O III] profile widths of PS and NPS sources also seems to be larger in the high luminosity region. This suggests also a link between the difference in [O III] profile widths and the source luminosities. In Sect. 4.4, we assess this relationship using a stacking analysis. Similar to Sect. 4.1, these results are robust to the errors in spectral indices. We find similar trends in PS and NPS sources with spectral index offsets, as shown in Figs. C.3e and C.3f. This shows that the relation between [O III] widths and radio spectral shape is a real trend and not a result of erroneous spectral shape classification.

We construct a similar plot for the high redshift PS and NPS sources, in Fig. 8. Here we find only marginal differences in the widths of the [O III] profiles of PS and NPS sources in any bin. The fraction of disturbed sources, is almost always larger for PS than NPS sources, at any L_1.4 GHz or L_[O III] bin. This suggests that in comparison to NPS sources, more disturbed [O III] in PS sources can also be found at 0.4 < z < 0.8, however the difference is only marginal. Similar to the low redshift sample, our results are consistent with the classification after spectral index offsets, shown in Figs. C.3g and C.3h. We do not present a stacking analysis below for the high redshift sample, since we did not find any significant difference between the profiles. But discuss this result further in Sect. 5.

4.4. Stacking analysis of PS and NPS sources

In the section above, we find that the difference in the widths of [O III] profiles of low redshift PS and NPS sources is related to L_1.4 GHz and L_[O III]. But this is hard to judge from Fig. 8 alone. We also want to observe these differences in the widths of average [O III] profiles of PS and NPS sources and quantify them. For this purpose, we use a stacking analysis. We stack the [O III] profiles for low redshift PS and NPS sources while controlling for L_1.4 GHz and L_[O III] individually.

To cover most of the source luminosity range discussed in Sect. 4.2 while keeping a significant number of sources in each group, we decided to split the sample into three groups of L_1.4 GHz from 10²³ to 10²⁶ W Hz⁻¹ and L_[O III] from 10³⁹ to 10⁴² erg s⁻¹. For consistency with Sects. 4.2 and 4.3, we keep the width of each luminosity group the same, that is 1 dex on the log scale. Since we want to study the average [O III] profiles of these sources and have a small number of objects in the highest luminosity groups, we have carefully tested whether the results are reliable for the stacked profile. We have removed three sources with high S/N spectra and extremely disturbed and broad [O III] profiles (v_broad up to −600 km s⁻¹ and FWHM_broad up to 1450 km s⁻¹) that dominate the stacked profile in the 10²⁵ − 10²⁶ W Hz⁻¹ group. We do this to avoid the results being driven by a few individual extremely disturbed sources, and to get a stacked profile more representative of the true average.

The stacked [O III] profiles are now shown in Fig. 9. For the groups in L_1.4 GHz, we find that the differences in the widths of stacked profiles become more prominent with increasing L_1.4 GHz. In the 10²³ − 10²⁴ W Hz⁻¹ group, the stacked profile of PS sources is marginally broader than the NPS sources. Although in Fig. 8 we see that PS sources in this L_1.4 GHz range have wider profiles and more disturbed [O III] than NPS sources. The difference is much more significant at 10²⁴ − 10²⁵ W Hz⁻¹ and the most at 10²⁵ − 10²⁶ W Hz⁻¹. We measure the widths of the profiles by fitting models made up of Gaussian components as discussed in Appendix A. Although the figures only show the 5007 Å component of the [O III] doublet line, while fitting the models we also use the 4958 Å component for accurate characterisation of the stacked profile. The best-fit model parameters are summarised in Table C.2. Here, we focus on the widths of the low-level broad components. In the L_1.4 GHz = 10²⁴ − 10²⁵ W Hz⁻¹ group, FWHM of the broad component is 1338 km s⁻¹ for PS and 1179 km s⁻¹ for NPS sources. At L_1.4 GHz = 10²⁵ − 10²⁶ W Hz⁻¹, the FWHM of the broad component is 1639 km s⁻¹ for PS and 1091 km s⁻¹ for NPS sources. Therefore, the difference between the FWHMs of PS and NPS sources increases from ≈160 km s⁻¹ to ≈550 km s⁻¹ with L_1.4 GHz. This suggests that the ‘difference’ in the impact of young and evolved radio AGN on [O III] gas also increases with L_1.4 GHz.

Fig. 9. Stacked [O III] profiles for the peaked and non-peaked spectrum sources in the low redshift sample, while controlling for L1.4 GHz (a−c) and L[O III] (d−f). The shaded regions show the 1σ errors on the stacked profiles, estimated by stacking 1000 bootstrapped samples. The best-fit model parameters for each profile are summarised in Table C.2. (a)−(c) Stacked profiles for PS and NPS sources while controlling for L1.4 GHz: (a) 1023 − 1024 W Hz−1, (b) 1024 − 1025 W Hz−1 and (c) 1025 − 1026 W Hz−1. Stacked profiles for PS and NPS sources while controlling for L[O III]: (d) 1039 − 1040 erg s−1, (e) 1040 − 1041 erg s−1 and (f) 1041 − 1042 erg s−1.

Next, we investigate the role of L_[O III], shown in the bottom panel of Fig. 9. Similar to L_1.4 GHz, we find no difference in the stacked profiles of the lowest luminosity group from 10³⁹ − 10⁴⁰ erg s⁻¹. At higher luminosities, the difference in the broad component widths decreases marginally with increasing L_[O III]. In the stacked profiles of L_[O III] = 10⁴⁰ − 10⁴¹ erg s⁻¹ group, FWHM of broad component is 1353 km s⁻¹ for PS and 981 km s⁻¹ for NPS sources. However in the highest luminosity group of L_[O III] = 10⁴¹ − 10⁴² erg s⁻¹, FWHM of the broad component is 1434 km s⁻¹ for PS and 1202 km s⁻¹ for NPS sources. The difference between the FWHMs decreases from ≈370 km s⁻¹ to ≈230 km s⁻¹ with L_[O III]. Therefore, we find no significant change in the average profiles of PS and NPS sources with increasing L_[O III] beyond 10⁴⁰ erg s⁻¹. Including the [O III] non-detections while stacking did not affect any results.

4.5. Stacking analysis of HFPS, LFPS and NPS sources

In the sections so far, we have found significant evidence to show that [O III] gas in PS sources is more disturbed than NPS sources in the low redshift sample, even while controlling for L_1.4 GHz and L_[O III]. This supports the picture of young radio AGN having a stronger impact on the surrounding [O III] gas than their evolved counterparts. However, as discussed in Sect. 3.1, PS sources are made up of HFPS sources, which would be the youngest in the sample, and LFPS sources which would be “less” young. Although we make a sharp distinction between HFPS and LFPS sources here, it is worth noting that the peaked sources have a continuous distribution in the colour-colour plot of Fig. 3. The distinction used here is only to make bins for stacking. If the impact on [O III] kinematics in radio AGN is linked to their evolutionary stage, in principle we should be able to trace the differences in the [O III] profiles of HFPS, LFPS and NPS sources in our low redshift sample. Since the number of HFPS and LFPS sources is small, 105 and 112 respectively, we have used a stacking analysis to test whether we can recover the difference in the average [O III] profiles of these groups.

In the stacking analysis of the previous section, we find that the profiles show significant differences in the L_1.4 GHz range of 10²⁴ − 10²⁶ W Hz⁻¹. Therefore, we have restricted the sample to this luminosity range. We also excluded the four strong PS sources mentioned in the section above, out of which one is HFPS and three are LFPS sources. The stacked profiles are shown in Fig. 10 and the best-fit parameters and bin sizes are summarised in Table C.2. We observe a clear difference in the stacked profiles of the three groups. HFPS sources show the broadest profile at the base, on both blue and red sides, and a broad component FWHM of 1642 km s⁻¹. The profiles then get narrower gradually as we move to LFPS and NPS sources, which have a broad component FWHM of 1333 km s⁻¹ and 1145 km s⁻¹, respectively. This shows that our sources’ average [O III] profile gets narrower as we go from the youngest to the most evolved sources, showing the changing impact on [O III] gas. It also shows that it is possible to trace the changing impact on the average [O III] profiles, using a radio spectral shape classification. Similar to the stacking analysis in the section above, including the [O III] non-detections while stacking did not affect our results.

Fig. 10. Stacked profile for HFPS, LFPS and NPS sources in the low redshift (0.02 < z < 0.4) sample covering a L1.4 GHz range of 1024 − 1026 W Hz−1. The shaded regions show the 1σ errors estimated using the bootstrapped samples. The best-fit model parameters and group sizes are summarised in Table C.2.

4.6. Radio sizes

In this section, we investigate the relation between ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ and the total radio sizes (physical) from LoTSS, plotted in Fig. 11 for the low and high redshift samples. We do not observe any clear trend between the two quantities. In the low redshift sample, 30 out of the 256 sources (11.7${}_{-1.7}^{+2.3}$%) smaller than 20 kpc and 19 out of the 167 sources (11.4${}_{-2.0}^{+3.0}$%) larger than 20 kpc are disturbed. Similarly in the high redshift sample, 5 out of the 50 sources (10.0${}_{-2.9}^{+6.0}$%) smaller than 20 kpc and 54 out of the 472 sources (11.4${}_{-1.3}^{+1.7}$%) larger than 20 kpc are disturbed. However, it is interesting that the median ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ at different sizes is consistently larger for PS sources than NPS sources in the low redshift sample, as shown by the dashed and dashdotted lines in Fig. 11. This shows that irrespective of the radio size, a PS source at the centre typically has more disturbed [O III] gas than an NPS source. We discuss the role of radio sizes and morphology further in Sect. 5.4.

Fig. 11. LoTSS major axis physical sizes (deconvolved) vs ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ for the (a) low redshift (0.02 < z < 0.4) and (b) high redshift (0.4 < z < 0.8) sample. The different colours show the peaked and non-peaked spectrum sources. Sources unresolved in LoTSS are marked with a horizontal arrow since the sizes are upper limits in these cases. Dashed and dashdotted lines show the median values for the peaked and non-peaked sources respectively, in bins of major axis sizes. The dotted horizontal line marks the ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ threshold for classifying the [O III] in a source as disturbed.

At this stage, we remind the reader that the radio spectral shape is determined using only the peak flux densities in the central region of the sources, therefore PS or NPS does not classify the spectrum over the entire extent of the radio source. This method allows us to identify a PS at the centre of large sources, which could be candidates for restarted radio AGN. We identify 138 out of 850 PS sources up to z = 0.8 that have extended emission in their LoTSS images on scales > 50 kpc. This is done by visual inspection of the LoTSS images and also includes the [O III] non-detections. Since the young source is likely ∼0.1 − 1 Myr old, this phase of activity could not have formed the extended emission observed on > 50 kpc scales. This suggests that the PS source and the extended emission represent two distinct phases of activity, and these sources are good restarted candidates. This sample could be useful to study the cumulative impact of AGN feedback on the host galaxies over multiple epochs of activity and will be followed up in detail in a future project.

4.7. Ionised gas densities of PS and NPS sources

Although the PS and NPS host galaxies are similarly massive, the difference in the [O III] kinematics of these sources could also be due to a richer medium available for interaction with the young jets in PS sources. To test this, we estimated the electron densities for these sources using the [S II]λλ6717, 30 Å doublet method. For sources that have a ≥3σ [S II] detection, the 6717 to 6730 Å line ratio is related to the electron density as (Sanders et al. 2016):

$$\begin{array}{c}\hfill n_{\mathrm{e}}(R)=\frac{cR-ab}{a-R},\end{array}$$(1)

where n_e is the electron density, R is the flux ratio of the 6717 to 6730 Å [S II] line, a = 0.4315, b = 2107 and c = 627.1. These values are estimated for a gas with an electron temperature of 10 000 K. Since this relation saturates at electron densities below ∼10 cm⁻³ (R < 1.3) and above ∼10 000 cm⁻³ (R > 0.45), we only used sources with 0.45 < R < 1.3, to stay within the theoretical line ratio limits of this relation. We note that this method has several limitations as it assumes a fully ionised medium whereas the gas could be partially ionised. It is also not sensitive to high-density gas (see a discussion in Davies et al. 2020). Studies using auroral and transauroral lines have found significantly higher densities of the ionised gas than those estimated using the [S II]λλ6716, 6731 Å ratio (Holt et al. 2011; Santoro et al. 2020; Rose et al. 2018; Davies et al. 2020).

Out of the 566 [S II] detections in the radio AGN sample, electron densities can be estimated for 301 sources using this method, out of which 66 are PS and 90 are NPS sources up to z ≈ 0.5. In Fig. 12, we plot distributions of the estimated electron densities. We find no significant difference in the PS and NPS distributions, confirmed by a 2-sample KS test (statistic = 0.14, p value = 0.40). Although the densities could only be estimated for a subset of the PS and NPS sources, it gives some insight into the ionised gas medium. The lack of any difference between the two suggests that the ionised medium is similarly dense in both groups. We discuss this in the context of jet-driven feedback in the next section.

Fig. 12. Electron density distributions for PS and NPS sources, that have ≥3σ [S II] detections.

5. Discussion

In this paper, we have investigated the [O III] gas kinematics for a large sample of 5720 radio AGN, over a L_1.4 GHz range of 10^22.5 − 10²⁸ W Hz⁻¹, out to redshift z = 0.8. Our careful selection of radio AGN allows us to probe the feedback from jets on the [O III] gas kinematics. Using the spectral shape over a wide frequency range (144 − 3000 MHz) we have separated young and evolved radio AGN, making it possible to study their impact in detail. Below, we discuss our results from Sect. 4 in the context of jet-driven feedback in radio AGN host galaxies.

5.1. Impact of jets over the radio AGN life cycle

In our analysis, we have been able to disentangle the effects of radio jets and radiation on disturbing the ionised gas. Comparison between the widths of [O III] profiles of our radio AGN sample and the control sample in Fig. 6a shows that the presence of jets in AGN host galaxies causes significantly more disturbed [O III] gas. Evidence for a positive correlation between the presence of radio emission and disturbed [O III] kinematics has been found before, for instance, by Mullaney et al. (2013) and recently by Escott et al. (2024). However, the radio emission in their sources could have different origins: star formation, shocks, or low-power jets. On the other hand, our sample is selected to have radio emission dominated by jets since we study jet-driven feedback on [O III].

With the spectral shape classification, we were able to go one step further and study the changes in the impact on [O III] at different evolutionary stages of radio AGN. This link between the radio spectral shape and the [O III] kinematics is the most important result of this paper. In Sect. 4.1, we found that [O III] in PS sources is ∼3 times more likely to be disturbed than NPS sources in the low redshift sample, even though both groups cover similar ranges of host galaxy properties. Since the radio emission in our sample is dominated by jets, this gives insight into the link between jets and [O III] kinematics. Our results show that the impact of radio jets on surrounding [O III] gas changes as they grow, with the strongest impact on the gas kinematics occurring when they are young. At evolved stages, [O III] is less disturbed. These trends support previous evidence in this context by K23, now confirmed with a sample that is ∼44 times larger. Since the young phase of activity typically lasts for ∼0.001 − 1 Myr (see Sect. 1), this suggests that the impact of radio AGN on the warm ionised gas is most extreme in the initial stages, which are a small fraction of their total lifetime (∼10 − 100 Myr). At later stages of evolution, the impact on the ionised gas is more ‘gentle’.

In the high redshift (z > 0.4) sample, we find only marginal evidence for more disturbed [O III] in PS sources (∼2 times more likely). Although less significant (2σ), the observed trend is still in agreement with the low redshift sample. We propose that this is likely due to contamination of the NPS sources by PS sources in the high redshift sample, as mentioned in Sect. 3.1. This contamination would lead to broader [O III] in NPS sources at high redshift. This can be seen in Fig. 6, where see that the fraction of sources with ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}>350$ km s⁻¹ is larger for NPS sources in the high redshift than the low redshift sample.

Radio and [O III] luminosity of an AGN are crucial parameters to determine the energy output of the system which couples with the gas in the host galaxy. Although PS and NPS sources in our sample cover a wide range in L_[O III] and L_1.4 GHz (see Fig. 1), we find that controlling for these parameters does not affect our results. At any L_[O III] or L_1.4 GHz, PS sources show wider and larger fraction of disturbed [O III] profiles than NPS sources. This is most evident in the low redshift sample. This shows that the difference we find in [O III] profiles of PS and NPS sources is not driven by source luminosities, but is indeed a result of the changing impact of jets in radio AGN along their life cycle.

Controlling for L_[O III] in the stacking analysis provides insight into the role of radiation in our sample (Fig. 9). We find that above L_[O III] = 10⁴⁰ erg s⁻¹, [O III] luminosity does not play a significant role in determining the difference between widths of stacked profiles of PS and NPS sources. PS sources are consistently broader than NPS in different L_[O III] groups. However, below 10⁴⁰ erg s⁻¹, we do not observe any difference between PS and NPS stacked profiles, likely due to the low S/N of the spectra. The strength of radiation is likely not significant for determining the difference in the impact of PS and NPS sources.

Similarly, from the stacked profiles in Fig. 9, we find that the ‘difference’ between the profile widths of PS and NPS sources, increases with L_1.4 GHz. That is, there is a positive link between the varying impact of radio AGN and the 1.4 GHz luminosity. Since L_1.4 GHz is a proxy for jet power (Cavagnolo et al. 2010; McNamara & Nulsen 2012; Hardcastle 2018), this points towards young and powerful jets being more effective at pushing the ionised gas clouds in their path to faster velocities. This is broadly in agreement with the results of some jet-ISM interaction simulations which predict that powerful jets would drill through the surrounding medium with ease while driving fast outflows (Mukherjee et al. 2016, 2018). These simulations also indicate that low-power jets, on the other hand, would struggle to break through the ISM and dissipate their energy by inducing turbulence in the ambient medium. It is not possible to make a direct comparison of these simulations with our results, and that requires spatially resolved optical spectra. However, the trends we observe for the sample appear to agree with these predictions.

The larger sample in this study compared to K23 means that we can also trace the changing impact in different regions of the colour-colour plot in Fig. 3. Our stacked profiles of HFPS, LFPS and NPS sources (see Fig. 10) show a gradual narrowing as we move from the youngest to the “less” young and the most evolved source. This can also be observed by tracking the change in the broad component FWHM of the stacked profiles. On average, within the PS source group, the impact of HFPS sources on [O III] is stronger than LFPS sources. This further supports our results on the varying impact on [O III] over the radio AGN life cycle. Such a detailed tracing of the impact at different stages of the life cycle has been possible due to the addition of low-frequency data (144 MHz) with LoTSS, which allowed the detection of a large number of PS sources. Although evidence for disturbed ionised gas has been found in powerful young AGN before (L_1.4 GHz ≳ 10²⁶ W Hz⁻¹, e.g. Gelderman & Whittle 1994; Holt et al. 2008, 2011; De Vries et al. 2009b; Santoro et al. 2020), our study expands this to a much larger sample of typical radio AGN, down to L_1.4 GHz ∼ 10²³ W Hz⁻¹. We find that the impact on ionised gas and the difference between young and evolved radio AGN is significant even down to L_1.4 GHz = 10²⁴ W Hz⁻¹ (confirming the findings for single objects like Murthy et al. 2022; Oosterloo et al. 2017; Audibert et al. 2019; Ruffa et al. 2022).

We further check that the differences we observe in [O III] profiles of PS and NPS sources are not due to any selection effects. PS sources are selected based on absorption in their radio spectra. Other than synchrotron self-absorption, this could be caused by free-free absorption from a highly dense ionised medium surrounding the radio jets. Therefore PS sources could also have more high-density ionised gas. AGN with such a dense medium in the central region would have a larger possibility of jet-ISM interactions. They would be more likely to show prominent broad wings in their spectra than NPS sources. However, as we show in Fig. 12 for a subset of sources, the ionised gas densities for PS and NPS sources are similar. Therefore, the result of [O III] kinematics with radio spectral shape is not due to a selection effect, but due to genuine differences in their interaction with the ionised gas.

5.2. Disturbed [O III] gas in NPS sources

It is intriguing that NPS sources consistently show less disturbed [O III] gas than PS sources, even though the jets are still active. Since the disturbed gas is likely driven by shocks due to jet-ISM interaction (Sutherland & Bicknell 2007; Mukherjee et al. 2018), one possibility is that by the time jets have grown to larger scales, the shocked gas in the central region has cooled down, and is not as disturbed as before. In this case, the strength of the broad component could be reduced, leading to narrower observer profiles in NPS sources. Indeed, using hydrodynamical simulations of jet-ISM interaction, Meenakshi et al. (2022) found that as the jets start and propagate, they shock ionise the clouds, pushing them to high velocities, and clear the region which is then photoionised by the AGN radiation. In their scenario, once the jets have evolved to larger scales and their direct interaction with the gas clouds has stopped, the shocked gas clouds would cool down. Another possibility is that on the time scale of the NPS sources (a few to tens of Myr), the disturbed [O III] gas clouds have been pushed by jet-ISM interaction to distances beyond the coverage of the 3″ SDSS fibre. In that case, as well, the gas within the SDSS fibre would appear less disturbed.

Determining the exact nature of the disturbed gas in NPS sources is beyond the scope of this paper. Testing these possibilities would require high-resolution spatially resolved optical spectra. Indeed, spatially resolved molecular gas kinematics has provided evidence for changing effects of jet-ISM interaction as the jets expand in radio-loud (for example PKS 0023−26 by Morganti et al. 2021a) and radio-quiet (for example Girdhar et al. 2024) AGN (see Morganti et al. 2023 for an overview). We are currently performing a study with integral field unit (IFU) data for a subset of these sources to examine the size and shape of the disturbed [O III] gas emission in more detail and understand its nature.

5.3. Role of L_1.4 GHz and L_[O III]

In Sect. 4.2, we find that the width and fraction of disturbed [O III] profiles increase with both L_1.4 GHz and L_[O III]. This trend is most evident for the low redshift sample, as shown in Fig. 7a. With our approach of separating the PS and NPS sources, we can observe an even stronger trend between L_1.4 GHz and [O III] profile widths. This can be seen for the low redshift sample in Fig. 8. Therefore L_1.4 GHz, and hence the radio jet power, plays an important role in determining the impact on [O III] kinematics. The highest luminosity radio AGN have the most disturbed gas. However, this trend is not very clearly observed in the high redshift sample.

It is worth noting that some studies of AGN samples have found slightly different trends between L_1.4 GHz and [O III] profile widths. For instance, Mullaney et al. (2013) found that the average widths of [O III] profiles were largest between L_1.4 GHz = 10²³ − 10²⁵ W Hz⁻¹, and therefore [O III] was less disturbed at luminosities below and above this range. Our analysis differs in the sample selection criteria, as mentioned before. Therefore, our sample is more suited for studying the role of radio luminosity in jet-driven feedback. Our additional selection criteria of spectral shape (PS and NPS sources) means that our sample is also not contaminated by flat-spectrum sources, which could affect the results from Mullaney et al. (2013) at L_1.4 GHz > 10²⁵ W Hz⁻¹.

Previous studies have reported a correlation between the [O III] and stellar velocity dispersions (σ_O III] vs σ_⋆), suggesting that the [O III] gas kinematics are determined by the bulge gravitational potential (for example Nelson & Whittle 1996; Boroson 2002; Sexton et al. 2020). Although there is a large scatter in this correlation, it is still important to be cautious about the role of the bulge in determining [O III] kinematics (Bennert et al. 2018; Le et al. 2023). Based on this correlation, some studies have used σ_⋆ to account for the gravitational kinematics component in σ_O III]. For instance, Woo et al. (2016) have used the quantity σ_O III]/σ_⋆ to study the impact of AGN feedback for a sample of 39 000 type 2 AGN out to z ∼ 0.3. They find that the correlation between L_1.4 GHz and σ_O III] vanishes when σ_O III] is normalised using σ_⋆. Using similar techniques for a sample of ∼6000 AGN out to z ∼ 0.4, Ayubinia et al. (2023) have found the same result between L_1.4 GHz and σ_O III]/σ_⋆. They also do not find any dependence of σ_O III]/σ_⋆ on the radio-loud/radio-quiet, compact/extended groups of radio AGN. These results have suggested that radio AGN have a limited impact on [O III] kinematics, and radio AGN-driven outflows are much weaker than accretion-driven outflows (although radio AGN can provide an extra boost to the disturbed [O III] gas).

For comparison, we test our results for PS and NPS sources in the low redshift sample using σ_O III]/σ_⋆, as shown in Figs. 13a and 13b. We again find that PS sources have consistently more disturbed [O III] kinematics. In this case, we define [O III] to be disturbed if σ_O III]/σ_⋆ > 1.5, however the specific value of this limit does not affect our broad result. This shows that there is an intrinsic difference in the [O III] kinematics of PS and NPS sources. We also note that using this approach to select disturbed [O III] sources does not affect our results. The fraction of disturbed sources at any source luminosity is in agreement with the result presented in Fig. 8. Another important result is the dependence of this ratio on L_1.4 GHz and L_[O III]. We find that PS (and to some extent NPS) sources with higher luminosities have larger σ_O III]/σ_⋆. This also supports our result for the dependence of [O III] kinematics on source luminosity. However, we did not use the σ_O III]/σ_⋆ approach for the entire paper because the large scatter between σ_O III] and σ_⋆ makes it less reliable. The difference in the relation between σ_O III]/σ_⋆ and L_1.4 GHz for our sample compared to Woo et al. (2016) can be attributed to sample selection effects, as they study type 2 AGN. Although Ayubinia et al. (2023) also study radio AGN, most of their sources lie below L_1.4 GHz = 10²⁴ W Hz⁻¹ whereas most of our sources lie above it. Indeed, in Figs. 13a and 13b, we do not observe any significant dependence on radio luminosity below this limit. Our separation of PS and NPS sources also helps emphasise the role of L_1.4 GHz.

Fig. 13. Distributions of σO III]/σ⋆ and ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$. Colour map showing the average σO III]/σ⋆ values for (a) peaked and (b) non-peaked sources in the low redshift (0.02 < z < 0.4) sample, in bins of L1.4 GHz and L[O III], similar to Fig. 7a. (c) Cumulative distributions of the average [O III] velocity dispersions. Comparison for compact (< 3″) and extended (> 3″) sources according to FIRST sizes in the low redshift sample. We also plot the compact NPS sources for comparison. The inset plot shows the asymmetry distributions for the [O III] profiles with multiple kinematic components. (d) Cumulative distributions for low redshift PS and NPS sources that are compact in FIRST.

From Fig. 7a, we also conclude that increasing only L_[O III] (which is proxy for the total radiative power of an AGN; Heckman et al. 2004) also increases the [O III] profile widths and fraction of disturbed sources. Most of the radio AGN are low excitation radio galaxies (LERGs), i.e. systems with radiatively inefficient accretion and low accretion rates (Best & Heckman 2012). Our results suggest that even in such systems, radiation strength plays an important role in driving feedback on [O III] gas kinematics. For a sample of ∼800 LERGs up to z ∼ 0.3, Singha et al. (2021) also found that the majority of sources with an [O III] outflow (only ∼1.5% of their sample) had high accretion rates. They concluded that outflows in LERGs are radiative pressure driven. We note a crucial difference between our studies regarding the identification of disturbed [O III]. Singha et al. (2021) focus only on the outflow detections, which are the most extreme form of feedback, whereas we study the average velocity dispersion for the entire radio AGN sample. However, there seems to be a general agreement between the results on the role of radiation pressure.

5.4. [O III] and radio source sizes

Finally, we also studied the dependence of ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ on the radio sizes from LoTSS in Sect. 4.6, and found no significant relation between the two. We find that sources larger and smaller than 20 kpc in LoTSS are equally likely to have disturbed [O III] in both low and high redshift samples. However, it is interesting that we find larger ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ values for PS than NPS sources, at all sizes in the low redshift sample (Fig. 11). This suggests that in radio AGN host galaxies, the impact on [O III] is not determined by the extent of the radio emission, but by the emission properties in the central region.

Many studies in literature have found links between the radio morphology of AGN and the multiphase gas kinematics. In neutral gas, H I absorption studies have found broader and more asymmetric line profiles in compact radio sources compared to extended sources (for e.g. Geréb et al. 2014, 2015; Maccagni et al. 2017). Since young radio AGN usually have a compact radio morphology, this suggested the presence of more disturbed H I gas in young radio AGN. In the ionised phase, Mullaney et al. (2013) have found evidence for more disturbed [O III] in AGN with compact radio morphology. Molyneux et al. (2019) use the type 1 and 2 AGN sample up to z ∼ 0.2 from Mullaney et al. (2013), and also find more disturbed [O III] in compact radio emission sources. They conclude that the SDSS spectra of their AGN are twice more likely to have broad [O III] profiles (FWHM > 1000 km s⁻¹) if the radio emission is confined within 3″ in their FIRST images. It is worth noting that the typical radio luminosity of their sources is lower by ∼2 orders of magnitude compared to ours, as they do not select for radio AGN. The radio emission in their sources can also be due to AGN-driven winds that drive outflows and shock the ISM to produce radio emission. This could also explain the correlation between radio size and the prevalence of outflows. Another major difference in our study is the way we select young radio AGN, using a spectral shape classification down to low frequencies of 144 MHz (thanks to the addition of LoTSS data) instead of source extent. This allows us to identify young phases of activity even at the centres of extended sources.

Still, it is worth comparing their results with our low redshift sample in the context of radio morphology and spectral shape. For this, we split our low redshift sample using the FIRST deconvolved major axis sizes into compact (< 3″) and extended (> 3″) sources. It is important to note from Fig. 5 that all the PS sources in the low redshift sample are smaller than 3″ in FIRST. Therefore the compact group consists of both PS and NPS sources but the extended group consists only of NPS sources. We compare their ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ distributions in Fig. 13c and find that they are significantly different (KS statistic = 0.22, p-value = 4.8 × 10⁻³). We also find that the fraction of sources with FWHM > 1000 km s⁻¹ is 4.5${}_{-0.9}^{+1.4}$% for compact and 1.4${}_{-0.4}^{+3.1}$% for extended FIRST sources. This is in agreement with the results of Molyneux et al. (2019).

However, we propose that this result in our sample is driven by the PS sources present in the compact group. We test for this by controlling for spectral shape and only using the NPS sources for the compact group, shown in Fig. 13c. We find that the distributions of compact and extended sources are now not significantly different (KS statistic = 0.13, p-value = 0.35). The fraction of sources with FWHM > 1000 km s⁻¹ also comes down to 1.5${}_{-0.5}^{+1.9}$% for compact sources. Therefore, controlling for spectral shape removes the difference between the [O III] profile widths of compact and extended FIRST sources. This shows that the difference found before is indeed driven by PS sources. Further, in Fig. 13d we even compare the ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ distributions of only compact PS and NPS sources, and find that they are still significantly different (KS statistic = 0.20, p-value = 1.6 × 10⁻³). We conclude that controlling for FIRST sizes does not affect our results on the [O III] profiles of PS and NPS sources.

Therefore, we find that the [O III] kinematics in our radio AGN sample is not significantly dependent on the FIRST or LoTSS sizes. Instead, we consistently find that the [O III] kinematics is most strongly related to the radio spectral shape. We conclude that young radio AGN picked using their peaked radio spectrum are more likely to show disturbed [O III] profiles than evolved radio AGN with a non-peaked spectrum.

6. Summary and conclusions

In this study, we have selected a sample of 5720 radio AGN up to z ∼ 0.8 where the radio emission is dominated by jets. We have shown that the presence of radio jets in AGN disturbs the [O III] gas more compared to AGN without jets. Using a wide-frequency radio spectral analysis, we have been able to identify sources in different stages of their life cycle. We have then established a link between the radio spectral shape (proxy for AGN evolutionary stage), and the [O III] kinematics. Our main findings are summarised below:

1. [O III] in PS sources is ∼3 times more likely to be disturbed than in NPS sources at z < 0.4. This shows that jets in young radio AGN have a stronger impact on the surrounding ionised gas. At z > 0.4, we find marginal evidence for PS sources having ∼2 times more disturbed [O III].

2. The fraction of disturbed sources and the average [O III] velocity dispersion is consistently higher for PS sources, even when controlling for 1.4 GHz and [O III] luminosity. This shows that the feedback is indeed linked to the evolutionary stage of radio AGN. The most extreme impact of radio AGN is likely limited to the initial stages of evolution.

3. The difference between the widths of stacked [O III] profiles of PS and NPS sources increases with 1.4 GHz luminosity. This shows that young and powerful jets have the strongest impact on [O III] gas.

4. A lower fraction of disturbed [O III] sources in the NPS group suggests that the shocked gas has cooled down by the time jets evolve to large scales. Although the shocked gas may have been pushed out to distances beyond the coverage of the SDSS fibre. Testing this hypothesis would require studying the spatially resolved gas kinematics.

5. Using a stacking analysis, we find that the average [O III] profile becomes narrower as the source spectral shape changes from completely inverted, to peaked and finally to non-peaked. We propose that this traces, in detail, the changes in the impact of jets along the radio AGN life cycle.

6. We find that the fraction of disturbed [O III] sources and average velocity dispersion increases with both 1.4 GHz and [O III] luminosity, even while controlling for one of the two. This highlights the role of the total energetic output in driving feedback and shows that in radio AGN, radiation pressure can also add to the effect of jets on ionised gas.

7. We do not find any strong dependence of [O III] kinematics on the source sizes. We find that any dependence is driven by PS and NPS sources, and controlling for spectral shape removes it. We suggest that the trends observed between disturbed [O III] and compact radio AGN are driven by PS sources.

Our study shows that the impact of jets on the host galaxy’s gas could change as the jets evolve. Understanding the exact changes in the ionised gas kinematics and morphology over galaxy scales, and the eventual impact on star formation, requires detailed studies of the spatially resolved gas and stellar properties. We plan to perform such studies using IFU data from MaNGA (Bundy et al. 2015) and the GMOS-N spectrograph on the GEMINI telescope, for a subset of sources from this sample.

Acknowledgments

We would like to thank the anonymous referee for their comments, that helped improve this paper. P.K. acknowledges support through an Emmy Noether Grant (German Research Foundation) of Dominika Wylezalek (WY 179/1-1). 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 are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefited from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland; The Istituto Nazionale di Astrofisica (INAF), Italy. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under a 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: Modelling the stellar continuum and emission line characterisation

To obtain the emission line spectra, we first modelled the stellar continuum and removed it from the observed SDSS spectra. Modelling was done using the STARLIGHT spectral synthesis code⁴ (Cid Fernandes et al. 2005, 2009). Single stellar population models of Bruzual & Charlot (2003) with young (1, 3.2, 5, 10, 25, 40 and 100 Myr), intermediate (0.29, 0.64, 0.9, 1.43, 2.5 and 5 Gyr) and old (11 and 13 Gyr) age populations were used. We assumed a Chabrier initial mass function, and three different metallicities: 0.004, 0.02 and 0.05 in mass-fraction notation. Before modelling the stellar continuum, the spectra were shifted to rest-frame wavelength, and the following emission lines were masked: [O II]λλ3726,28Å, Hδ, Hγ, Hβ, [O III]λλ4958,5007Å, [O I]λ6300 Å, [N II] λλ6548,84Å, Hα and [S II]λλ6717,30 Å. A visual inspection revealed that this procedure could reproduce the stellar continuum emission of the sources well, and can be used for subtraction to obtain the emission line spectra.

After obtaining the continuum subtracted spectra, modelling of the line profiles was done using Gaussian functions. For reference, we used the [O III] doublet which is usually the strongest emission line. We used kinematic components made up of two Gaussian functions, with a fixed amplitude ratio ($A_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]4958}=\frac{1}{3}\times A_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]5007}$) and equal widths (${\sigma }_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]4958}={\sigma }_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]5007}$) according to atomic physics. We started with a single kinematic component and added up to two components depending on whether adding an extra component significantly improved the fit. We judged this using an F-test, and considered an extra component to significantly improve the fit, if it improved the χ² at > 99% confidence level. We considered a source as an [O III] detection if its amplitude was ≥3σ, where σ is the standard deviation of the residuals after subtracting the best-fit model. We obtained 2 363 [O III] detections, out of which 410 required double and 120 required triple components.

For sources with [O III] detections, we then fit the Hβ line. We used the [O III] best-fit model and allowed one more component to account for any possible broad line region (BLR) emission in the Hβ line from type 1 AGN. Out of all the Hβ detections, 173 sources required an extra BLR component in the Hβ line. Similarly using the best-fit [O III] model, we then fit the [N II]λλ6548,84 Å and the Hα lines. The separation between the [N II] and Hα lines, and the relative amplitudes of the two components of the [N II] line ($A_{[\mathrm{N}\mathrm{I}\mathrm{I}]5648}=\frac{1}{3}\times A_{[\mathrm{N}\mathrm{I}\mathrm{I}]5684}$) were fixed. We again allowed for an extra broad component in the Hα line, whenever one was needed for the Hβ line. Lastly, we fit the [S II]λλ6717,30 Å doublet, with a fixed velocity separation and the same velocity dispersion and number of components as the best fit [O III] model. The number of different emission line detections in the full sample and the radio AGN sample selected in the next section, are summarised in Table 1.

Next, we estimated the emission line fluxes by adding up fluxes of the Gaussian components in the best-fit model. In the case of Hβ and Hα line fits that required an extra BLR component, we neglected the flux of the BLR component as we did not want the light from the BLR region to affect our fluxes. For non-detections, we used a 3σ value to estimate an upper limit on the fluxes.

We removed 22 sources from the sample where the spectra were corrupted at the [O III] wavelengths. We further removed 30 sources where the Hβ profile required an extra broad component, and the emission was strong enough to leak into the [O III] profile. This contaminates the profile width and flux estimates. Finally, we were left with a sample of 6 504 sources.

[O III] profile characterisation and stacking: The multi-component fitting routine makes it challenging to compare [O III] profiles with different number of kinematic components in their best-fit model. To overcome this, we estimated a flux weighted average velocity dispersion ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$, similar to Mullaney et al. (2013), given as

$$\begin{array}{c}\hfill {\overline{\sigma }}_{[\mathrm{O}iii]}={\left({\displaystyle \sum _i}{(f(i)\times {\sigma }_{[\mathrm{O}iii]}(i))}^2\right)}^{\frac{1}{2}},\end{array}$$(A.1)

where f(i) is the flux ratio of the i th component to the total flux of the line and σ_O III](i) is the velocity dispersion of the i th component in the best-fit model. Using this one parameter we can compare profiles with different numbers of components. In addition, it is also useful to assess the shape of the line profiles and determine any asymmetry, if present. The asymmetry helps trace kinematically disturbed gas in the sources. We measured the asymmetry of the lines using the definition from Liu et al. (2013), given as

$$\begin{array}{c}\hfill A=\frac{(v_{90}-v_{50})-(v_{50}-v_{10})}{v_{90}-v_{10}},\end{array}$$(A.2)

where v₉₀, v₅₀, and v₁₀ are the velocities at 90, 50 and 10 per cent of the cumulative flux. In this convention, a heavily blueshifted profile has A < 0 whereas a heavily redshifted profile has A > 0 and a symmetric profile has A = 0.

In this paper, we also use a stacking analysis to measure the average differences between [O III] profiles of different source groups. Using a stacking approach, we can control for various source properties and disentangle their role in driving the [O III] kinematics. It also suppresses the noise (by 1/$\sqrt{N}$ for N samples) thus allowing the detection of low-level broad components, which would be otherwise missed in individual fits. Stacking requires accurate spectroscopic redshift measurements to align the rest-frame spectra, which we have for the entire sample. We estimated the stacked spectrum using the equation

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

where S(i, λ) is the individual source spectrum, σ(i) is the noise in a line-free region of the source spectrum and S(λ) is the stacked spectrum. We estimate errors on the average stacked spectrum using a bootstrapping approach. We randomly select sources in each bin to create a bootstrapped sample and stack their spectra, repeating this process 1000 times. We use the standard deviation of the distribution at each wavelength as the error for that wavelength point.

We also estimate the impact of errors in redshift on the stacked [O III] profiles. To do this, we add a random value to each redshift, to obtain a redshift with an offset. This value is drawn randomly from a Gaussian distribution with unit area, zero mean and standard deviation equal to the error in reported in SDSS for that redshift. We then re-perform the stacking analysis to test the impact of redshift errors. We find that it has a negligible impact on the final stacked profiles.

Appendix B: Selecting radio AGN

To study the relation between radio AGN life cycle and their impact on [O III] gas, we first need to select a clean sample of radio AGN. Radio emission from a galaxy can originate either from the AGN itself or star formation (SF). In the local universe, star-forming galaxies dominate a radio-selected sample up to L_1.4GHz ≈ 10²³ W Hz⁻¹ (e.g. Sadler et al. 2002; Best et al. 2005). This makes it important to separate them from radio AGN. Even AGN identified using optical emission line ratios can have radio emission dominated by star formation in the host galaxy. In this section, we discuss the approach used for selecting a clean sample of radio AGN from the 6 504 radio selected sources in Sect. 2. We followed the method outlined in Best et al. (2005), Best & Heckman (2012) and Sabater et al. (2019), which combines four diagnostic methods. The plot for each diagnostic is shown in Fig. C.1.

First is the ‘D_n(4000) versus L_1.4GHz/M_⋆’ method, shown in Fig. C.1 and developed by Best et al. (2005). D_n(4000) is the 4000 Å break in the spectrum of the galaxy (sensitive to the presence of a young stellar population), and L_1.4GHz/M_⋆ is the ratio of the radio luminosity to the stellar mass. This method relies on the idea that star-forming galaxies with various star formation histories, would lie in a similar region on a plot of these two quantities. This is because both of these quantities are dependent on the specific star formation rate. However, radio AGN have jets that add to the radio emission (L_1.4GHz), allowing us to separate the two populations. Best et al. (2005) obtained the track of a galaxy with an exponentially decreasing star formation rate on this diagram, which can be used to separate the two sources of radio emission. They then added 0.225 to the D_n(4000) values of this track to construct a conservative division. This division was later modified by Kauffmann et al. (2008), who used a straight line division at L_1.4GHz/M_⋆ < 12.2. For our sample, we estimated D_n(4000) as the flux ratio between 4000-4100 Å and 3850-3950 Å. We used the stellar masses from the ‘Portsmouth SED-fit’ catalogue by Maraston et al. (2013), which were available for 5 535 sources. The stellar masses were estimated from broad-band spectral energy distribution (SED) fitting of stellar population models using u, g, r, i and z band photometry and spectroscopic redshifts from SDSS. We classified the sources as radio AGN if they lie above the criterion from Kauffmann et al. (2008)⁵, and as star-forming (i.e. galaxies where the radio emission is dominated by star formation) if they lie below it.

The second diagnostic uses L_Hα versus L_1.4GHz to classify sources, shown in Fig. C.1. It is based on the idea that for star-forming galaxies, both Hα and radio luminosities are related to the star formation rate of massive stars. Therefore sources with excess radio emission (coming from jets) can be separated from sources where both Hα and radio luminosities come from star formation. We used the division from Best & Heckman (2012), which is a conservative division designed to produce clean radio AGN samples, and given as log(L_Hα/L_⊙) = 1.12× (log(L_1.4GHz/W Hz⁻¹)–17.5). We classify sources with larger radio luminosity than predicted from this relation as radio AGN and the rest as star-forming. We use upper limits on Hα luminosities for non-detections. Overall, we were able to classify 4 027 sources using this diagnostic.

The third diagnostic is the BPT diagram (Baldwin et al. 1981) shown in Fig. C.1, which uses the [N II]/Hα and [O III]/Hβ emission line ratios to separate AGN and star-forming sources of ionisation. We used the relation from Kauffmann et al. (2003) to classify sources as AGN if they lie above the relation, and star forming if they lie below it. Although we do not have emission line detections for the majority of the sample, it can still provide useful information for sources with upper limits. We classified 1 503 sources using the BPT diagram.

The fourth diagnostic is the WISE colour-colour diagram, shown in Fig. C.1. Star-forming galaxies can be separated from radio AGN host galaxies (typically ellipticals with low levels of star formation) in their W2 − W3 colour (Yan et al. 2013). Herpich et al. (2016) used a division of W2 − W3 (AB) = 0.7 to separate galaxies with and without star formation, however, the division is not very clean. Sabater et al. (2019) modified it to W2 − W3 (AB) = 0.8 for LoTSS sources following a comparison with the radio AGN in the Herschel-ATLAS North Galactic Plane field selected by Gürkan et al. (2018). For this work, we used the same division as Sabater et al. (2019). Since this is a crude division, it has less weight in the overall classification and we use it only to resolve cases where other diagnostics give contradictory results. Since we did not perform k-corrections to the WISE magnitudes, we use them only for z < 0.3 sources. We classified 2 019 sources using this method.

In each diagnostic, a source can be labelled as an AGN (or radio AGN), SF/radio-quiet AGN or unclassified, leading to 81 different possible combinations. A summary of the source classifications by each diagnostic is given in Table B.1. We use these combinations to give a final classification to each source, following the approach of Sabater et al. (2019) and Best & Heckman (2012). We refer the reader to these papers for further details on the classification. Sources that could not be classified with any diagnostic, or were classified as SF using only one diagnostic, were assigned a final classification of radio AGN if their 1.4 GHz luminosity was greater than 10²⁵ W Hz⁻¹. The different combinations, number of sources in each, and their final classification are summarised in Table C.1. Overall, we classified 5 720 sources as radio AGN, 726 as SF/radio-quiet AGN and 58 were unclassified.

Appendix C: Supplementary tables and figures

Fig. C.1. Diagnostic plots for selecting radio AGN. Red points mark the sources classified as radio AGN after combining all four diagnostics, and black points show the star-forming/radio-quiet AGN (SF/RQAGN) sources. (a)Dn(4000) vs L1.4GHz/M⋆ plot. The dashed curve marks the division used for radio AGN and SF/RQAGN from Best & Heckman (2012). (b)LHα versus L1.4GHz plot, with the separation line from Best & Heckman (2012). (c) BPT diagram for the sources. The dashed curve shows the semi-empirical relation from Kauffmann et al. (2003). (d) WISE colour-colour plot for sources with the division line from Sabater et al. (2019).

Fig. C.2. Distributions of host galaxy properties for the peaked and non-peaked radio AGN in the high redshift sample. Only the [O III] detections are shown. Vertical lines show the mean values of the distributions.

Fig. C.3. Results for the radio AGN sample after adding offsets to spectral indices. (a)−(b) Comparison of spectral indices before and after adding the offsets described in Sect. 3.1. Cumulative distributions of ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ after reclassifying sources with the spectral index offset for the (c) low redshift and (d) high redshift sample. Colour maps showing the average ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ values and the fraction of disturbed sources in bins of L1.4GHz and L[O III], again after reclassifying sources with the spectral index offset. (e) PS and (f) NPS sources in the low-redshift (0.02 < z < 0.4) sample. (g) PS and (h) NPS sources in the high-redshift (0.4 < z < 0.8) sample.

All Tables

All Figures

Fig. 1. Properties of the radio AGN sample. (a): L1.4 GHz vs L[O III] for the full radio AGN sample, with [O III] detections and non-detections shown separately. In the case of [O III] non-detections, the contour plot shows the upper limits for L[O III]. The green hatched region marks the luminosity space occupied by the control sample selected from Mullaney et al. (2013), described in Appendix B. The dashed and dotted lines show the radio luminosity expected from the interaction of AGN-driven winds with the ISM from the model of Nims et al. (2015), for 1% and 10% of the shock energy going into relativistic electrons. (b) Distributions of host galaxy properties for the peaked and non-peaked radio AGN in the low redshift (0.02 < z < 0.4) sample. Only the [O III] detections are shown. Vertical lines show the mean values of the distributions. The same plot for the high redshift sample is shown in Fig. C.2.

In the text

	Fig. 2. Redshift vs total L1.4 GHz for the selected from the sample in Appendix B. Peaked and non-peaked sources are marked with different colours. The sample is split about the dashed line at z = 0.4 into low and high redshift sources.
In the text

Fig. 3. Colour-colour plot of spectral indices for the radio AGN sample. The different quadrants are labelled with the spectral shapes of the sources lying in that region. The spectral shapes are also shown above the labels. The red square marks the region where flat-spectrum (FS) sources lie. A 1:1 line is shown in blue. While studying the link of spectral shape with gas kinematics, only HFPS, LFPS and NPS sources are considered.

In the text

Fig. 4. Fraction of PS and NPS sources in redshift (z) bins for the full radio AGN sample. The fraction is the number of PS (or NPS) sources divided by the total number of sources in the bin (of all spectral shapes). The horizontal error bars show widths of the redshift bins used to estimate the fraction and the vertical error bars show uncertainties on the fraction of PS and NPS sources in the bin.

In the text

	Fig. 5. Deconvolved major axis sizes from LoTSS vs. FIRST for the full radio AGN sample up to z ∼ 0.8. The red, blue and grey points mark the PS, NPS and the rest of the sources in the sample. The horizontal and vertical dashed lines mark the 3″ size. A 1:1 ratio is marked by the dashdotted line.
In the text

Fig. 6. Cumulative distributions of the flux-weighted average [O III] velocity dispersions for [O III] detections. The inset plots show the asymmetry distributions for [O III] profiles with multiple kinematic components in every group. (a): Low redshift (0.02 < z < 0.4) radio AGN sample (PS + NPS) and the control sample of optical AGN up to z ∼ 0.4 from Mullaney et al. (2013), described in Appendix B. (b): Low redshift PS, NPS sources and the control sample. (c): High redshift (0.4 < z < 0.8) PS and NPS sources.

In the text

Fig. 7. Colour map showing the average ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ values for the (a) low redshift (0.02 < z < 0.4) and (b) high redshift (0.4 < z < 0.8 sample in bins of L1.4 GHz and L[O III]. The colour of each bin represents the average velocity dispersion in that bin. The number of sources (the bin size) is mentioned at the top of each bin. Fraction of sources with disturbed [O III] are mentioned below the bin size, only for bins that had at least one disturbed source.

In the text

	Fig. 8. Same plot as Fig. 7. The top row shows the plots for (a) peaked spectrum and (b) non-peaked spectrum sources in the low redshift (0.02 < z < 0.4) sample. The bottom row shows the plots for (c) peaked spectrum and (d) non-peaked spectrum sources in the high redshift (0.4 < z < 0.8) sample.
In the text

Fig. 9. Stacked [O III] profiles for the peaked and non-peaked spectrum sources in the low redshift sample, while controlling for L1.4 GHz (a−c) and L[O III] (d−f). The shaded regions show the 1σ errors on the stacked profiles, estimated by stacking 1000 bootstrapped samples. The best-fit model parameters for each profile are summarised in Table C.2. (a)−(c) Stacked profiles for PS and NPS sources while controlling for L1.4 GHz: (a) 1023 − 1024 W Hz−1, (b) 1024 − 1025 W Hz−1 and (c) 1025 − 1026 W Hz−1. Stacked profiles for PS and NPS sources while controlling for L[O III]: (d) 1039 − 1040 erg s−1, (e) 1040 − 1041 erg s−1 and (f) 1041 − 1042 erg s−1.

In the text

	Fig. 10. Stacked profile for HFPS, LFPS and NPS sources in the low redshift (0.02 < z < 0.4) sample covering a L1.4 GHz range of 1024 − 1026 W Hz−1. The shaded regions show the 1σ errors estimated using the bootstrapped samples. The best-fit model parameters and group sizes are summarised in Table C.2.
In the text

Fig. 11. LoTSS major axis physical sizes (deconvolved) vs ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ for the (a) low redshift (0.02 < z < 0.4) and (b) high redshift (0.4 < z < 0.8) sample. The different colours show the peaked and non-peaked spectrum sources. Sources unresolved in LoTSS are marked with a horizontal arrow since the sizes are upper limits in these cases. Dashed and dashdotted lines show the median values for the peaked and non-peaked sources respectively, in bins of major axis sizes. The dotted horizontal line marks the ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ threshold for classifying the [O III] in a source as disturbed.

In the text

	Fig. 12. Electron density distributions for PS and NPS sources, that have ≥3σ [S II] detections.
In the text

Fig. 13. Distributions of σO III]/σ⋆ and ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$. Colour map showing the average σO III]/σ⋆ values for (a) peaked and (b) non-peaked sources in the low redshift (0.02 < z < 0.4) sample, in bins of L1.4 GHz and L[O III], similar to Fig. 7a. (c) Cumulative distributions of the average [O III] velocity dispersions. Comparison for compact (< 3″) and extended (> 3″) sources according to FIRST sizes in the low redshift sample. We also plot the compact NPS sources for comparison. The inset plot shows the asymmetry distributions for the [O III] profiles with multiple kinematic components. (d) Cumulative distributions for low redshift PS and NPS sources that are compact in FIRST.

In the text

Fig. C.1. Diagnostic plots for selecting radio AGN. Red points mark the sources classified as radio AGN after combining all four diagnostics, and black points show the star-forming/radio-quiet AGN (SF/RQAGN) sources. (a)Dn(4000) vs L1.4GHz/M⋆ plot. The dashed curve marks the division used for radio AGN and SF/RQAGN from Best & Heckman (2012). (b)LHα versus L1.4GHz plot, with the separation line from Best & Heckman (2012). (c) BPT diagram for the sources. The dashed curve shows the semi-empirical relation from Kauffmann et al. (2003). (d) WISE colour-colour plot for sources with the division line from Sabater et al. (2019).

In the text

	Fig. C.2. Distributions of host galaxy properties for the peaked and non-peaked radio AGN in the high redshift sample. Only the [O III] detections are shown. Vertical lines show the mean values of the distributions.
In the text

Fig. C.3. Results for the radio AGN sample after adding offsets to spectral indices. (a)−(b) Comparison of spectral indices before and after adding the offsets described in Sect. 3.1. Cumulative distributions of ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ after reclassifying sources with the spectral index offset for the (c) low redshift and (d) high redshift sample. Colour maps showing the average ${\overline{\sigma }}_{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}$ values and the fraction of disturbed sources in bins of L1.4GHz and L[O III], again after reclassifying sources with the spectral index offset. (e) PS and (f) NPS sources in the low-redshift (0.02 < z < 0.4) sample. (g) PS and (h) NPS sources in the high-redshift (0.4 < z < 0.8) sample.

In the text