© The Authors 2025

1. Introduction

Much effort over the past decade, both in observations and simulations, has been devoted to the co-evolution of galaxies and their central supermassive black holes (SMBHs). One of the first pieces of evidence leading to this investigation was the tight correlation between the masses of SMBHs and the stellar velocity dispersion of their host galaxy bulges (e.g., Magorrian et al. 1998; Ferrarese & Merritt 2000). Additional indirect evidence supporting the co-evolution scenario is the observation that the cosmic growth history of SMBHs and the galaxy population is quite similar (e.g., Silverman et al. 2008). They both peak at redshift z ∼ 2 and decline rapidly after z ∼ 1 (for a review, see Madau & Dickinson 2014). These observational results indicate that there are physical processes linking the growth of the SMBHs and their host galaxies, even though their sizes differ by orders of magnitude (parsec vs kiloparsec). However, the SMBH–galaxy co-evolution picture is still not fully understood in detail.

Accreting SMBHs are known as active galactic nuclei (AGN), and they are often very luminous with emission across the full spectral energy distribution (SED) (for a review, see Padovani et al. 2017). These features are believed to arise from energy produced by matter that is accreting onto the central SMBH. This matter may reside in an accretion disk and give rise to very strong broadened emission lines from photo-ionized gas (BLR; e.g., Baldwin 1997), a warm dusty torus contributing to near-infrared to mid-infrared continuum (e.g., Barvainis 1987), outer photoionized gas emitting high excitation narrow lines (NLR; e.g., Baldwin et al. 1981), and a corona surrounding the accretion disk that contributes more radiation in X-rays (e.g., Sun & Malkan 1989). There may also be radio jets created by relativistic charged particles (e.g., Begelman et al. 1984).

The term AGN feedback describes the impact of AGN on the evolution of their host galaxies. A negative feedback scenario is widely favored, in which the energy from the accreting material can heat or eject the gas in the host galaxies, and then prevent star formation globally (e.g., Fabian 2012, and references therein). For most numerical and semi-analytic simulations, negative AGN feedback is crucial in massive halos to prevent galaxies from growing too massive (e.g., Di Matteo et al. 2005; Croton et al. 2006). On the other hand, the observed AGN-starburst connection in Type II Seyfert galaxies and infrared AGN indicate that positive AGN feedback may also exist (e.g., Heckman et al. 2004). The outflows and shocks from AGN may also create dense gas regions and boost star formation (e.g., Salomé et al. 2015; Schutte & Reines 2022). Alternatively, the correlation between AGN activity and star formation may arise because star formation and SMBH accretion share a common gas supply (e.g., Dai et al. 2018; Li et al. 2024). Recent theoretical models suggest that the dominant feedback mode transitions from positive to negative as redshift decreases (e.g., Magliocchetti 2022; Silk et al. 2024).

In this paper, we focus on the AGN with radio jet features, known as radio AGN. For this population, kinetic (or jet, radio) mode is considered as the main feedback mechanism. These AGN mainly output mechanical energy through powerful radio jets, heat the surrounding gas, and produce weak emission lines in the interstellar medium of the host galaxy (e.g., Heckman & Best 2014). Radio AGN are more likely to be found in massive elliptical galaxies (Best et al. 2005a), and in group or cluster environments (Croston et al. 2019) than AGN that emit most of their energy in the optical or UV. They usually have low black hole accretion rates and are lack of dusty torus near the black hole (Heckman & Best 2014). In galaxy-scale simulations, various AGN jet models have been found, in principle, to be capable of quenching a galaxy and preventing gas from cooling at the centers of dark matter halos (e.g., Dubois et al. 2010; Yang & Reynolds 2016). However, the challenge remains to connect the critical free parameters in these simulations to observational data.

Radio AGN themselves have diverse observational properties. They are historically divided into Fanaroff & Riley (FR; Fanaroff & Riley 1974) types by radio morphology, where in general FR I, FR II, and FR0 (Garofalo & Singh 2019) represent edge-darkened, edge-brightened, and compact jets, respectively. Recent works suggested that the excitation properties may be more fundamentally connected to the SMBH accretion rate than the morphology (e.g., Best & Heckman 2012). Radio AGN with strong optical emission lines (e.g., [O III]λ5008) are classified as high-excitation radio galaxies (HERGs) with ongoing star formation activity in their hosts. The other population without strong emission lines are called low-excitation radio galaxies (LERGs).

The radio flux of normal star-forming galaxies with no AGN is dominated by the synchrotron continuum emitting from relativistic particles moving through magnetic fields. These particles are mainly accelerated by supernova (Biermann 1976) and SMBH accretion. The incidence of supernovae is related to the birth rate of short-lived, high-mass stars (M > 8 M_⊙), which is thought to be the origin of the tight correlation between star formation rate (SFR) and radio luminosity (hereafter radio-SFR relation). The radio-SFR relation was first observed through the constant ratios of radio and far-infrared luminosity in SFGs (e.g., de Jong et al. 1985; Condon et al. 1991), and has also been confirmed by diverse tracers of SFR (e.g., Wang et al. 2016; Brown et al. 2017; Gürkan et al. 2018; Smith et al. 2021; Heesen et al. 2022; Best et al. 2023; Das et al. 2024; Shenoy et al., in prep.). Most radio AGN classification methods are based on the excessive total radio luminosities compared to the star formation properties. Thus, the excess of radio luminosity compared to the star formation rate is used to derive the luminosity of the radio AGN jet. There are also other approaches to classification using infrared colors and spectral indices (e.g., Sabater et al. 2019, and reference therein), but these methods cannot quantitatively decompose the radio contribution from star formation and AGN jets.

Our current knowledge about the population of radio AGN is based on several all-sky radio continuum surveys at 1.4 GHz, such as the NRAO Very Large Array Sky Survey (NVSS, Condon et al. 1998) and the Faint Images of the Radio Sky at Twenty-Centimeters (FIRST, Becker et al. 1995). In the low-redshift universe (z < 1), the sensitivity of these surveys probes radio AGN with luminosities greater than L_1.4 GHz ∼ 10²³ W Hz⁻¹.

The new generation of radio telescopes and instruments, such as the LOw-Frequency ARray (LOFAR, van Haarlem et al. 2013), the Australian Square Kilometre Array Pathfinder (ASKAP, Hotan et al. 2021, MeerKAT, Jonas & MeerKAT Team 2016), are now able to discover much fainter radio AGN. The LOFAR Two Meter Sky Survey (LoTSS, Shimwell et al. 2017) is one of the deepest ongoing radio continuum large-sky surveys, which aims to observe the whole northern sky at frequencies 120−168 MHz. Compared to NVSS and FIRST at 1.4 GHz, LoTSS observes at lower frequencies, which are more sensitive to old, diffuse synchrotron emission and less affected by the dust continuum thermal emission (e.g., Penney et al. 2020). For an extragalactic source with a typical radio spectra index, it is also about one order of magnitude deeper than NVSS and FIRST, and sensitive enough to detect faint radio sources which are mainly star-forming galaxies (SFGs) at low-redshifts (Williams et al. 2019; Hardcastle et al. 2023).

To study the properties of radio AGN host galaxies, we need observations at wavelengths where emission from the stellar component and the interstellar medium (ISM) become dominant. The optical integral field unit (IFU) technology is particularly useful to study star formation and stellar populations of galaxies. Spatially resolved galaxy properties, which are not measurable in single-fiber spectroscopic surveys, are now available from IFU observations. Within the sky coverage of LoTSS, Mapping nearby Galaxies at Apache Point Observatory (MaNGA, Bundy et al. 2015) is the largest optical IFU survey that has been completed. MaNGA includes spatially resolved spectroscopic observations of about ten thousand galaxies, from which we can build a statistically significant sample of radio AGN host galaxies. The data available for MaNGA galaxies also include classifications of morphology and environment, and we can thus investigate whether these parameters affect radio AGN properties.

The combination of LoTSS and MaNGA allows us to investigate radio AGN feedback on subgalactic scales, particularly for faint radio AGN. The latest public data release, LoTSS-DR2, and the half MaNGA sample (DR16), were used by Mulcahey et al. (2022) to build a sample of 307 radio AGN and to compare the star formation, ionization maps, and age gradients of radio AGN hosts with control galaxies without radio AGN. It was found that radio AGN are generally quiescent, but are not more quenched than their controls. The authors concluded that radio AGN help host galaxies maintain quiescence but are not responsible for switching off star formation completely. Using a similar sample, Zheng et al. (2023) studied the stellar kinematics of radio AGN host galaxies. They found that the low angular momentum is important for galaxies to host luminous radio AGN.

In this paper, we update the sample using the final MaNGA data release (Abdurro’uf et al. 2022), which includes ∼6000 galaxies falling within the LoTSS-DR2 sky footprint. This roughly doubles the size of the radio AGN sample compared to Mulcahey et al. (2022). First, we aim to utilize this IFU–radio sample, which is currently the largest, to provide an improved derivation of the radio continuum–SFR relation and its intrinsic scatter. Because MaNGA provides full-coverage Hα maps with dust attenuation information from the Balmer decrement, we can overcome the limitation of single-fiber observations that only offer nuclear spectroscopic information for low-redshift galaxies. This also allows the radio continuum–SFR relation to be better constrained. Based on this relation, the excess radio luminosity compared to the prediction from its star formation rate is a physical and reliable criterion to classify radio-excess AGN (hereafter RDAGN). This relation is also used to decompose the star formation and AGN contributions to the total radio flux (i.e., to calculate the jet luminosities or upper limits).

With this RDAGN host galaxy sample in hand, our aim is to study star formation in the RDAGN hosts and to investigate how the radio jets affect the subgalactic properties of their hosts. MaNGA provides maps of the optical spectra (3600−10 300 Å) out to 1.5× the effective radius (R_e). The emission and absorption lines can be used to analyze the properties of the ionized interstellar medium (ISM) and resolved stellar populations within the host galaxy. We also use the results derived from the full spectrum fitting method (Wilkinson et al. 2017) to analyze the star formation histories (SFHs) in different regions of the host galaxies.

The structure of this paper is as follows. In Sect. 2, we introduced the LoTSS-MaNGA parent sample cross-matching and the calculation of physical parameters. In Sect. 3, we built the radio continuum–SFR relation and used it to build RDAGN and control samples. In Sect. 4, we examined the global star formation in RDAGN hosts. Additionally, we analyzed the spatially resolved mass assembly history and emission line features of these hosts in different radial bins. We discussed the main results and the implications of our study as well as the future outlook in Sect. 5. Section 6 contains an overall summary of our main results.

Throughout this paper, we converted all radio luminosity and related relations to rest-frame 144 MHz assuming a radio spectral index α = −0.7 (Condon et al. 2002): L_ν1/L_ν2 = (ν₁/ν₂)^−0.7. We adopted a flat cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7. All stellar masses and star formation rates are based on the Chabrier initial mass function (Chabrier 2003). When comparing literature results based on different IMFs, we followed Madau & Dickinson (2014) and set the factor (for stellar masses and star formation rates) as Salpeter:Kroupa:Chabrier = 1:0.67:0.63.

2. Data

2.1. MaNGA

MaNGA is part of the Sloan Digital Sky Survey IV (SDSS-IV, Bundy et al. 2015), and is the largest completed IFU survey. The final data products are released in SDSS DR17 (Abdurro’uf et al. 2022), which consist of ∼10 000 unique galaxies with good quality IFU coverage out to at least 1.5 times the r-band effective radius (R_e). MaNGA galaxies are selected from the NASA-Sloan Atlas catalog (NSA, Blanton et al. 2011; Wake et al. 2017)¹, and thus have measurements of many physical parameters based on the UV-optical photometric data, such as absolute magnitude, effective radius, and position angle. We adopt the global stellar mass from the NSA catalog to avoid having to correct for the fact that the MaNGA IFU sizes vary. The stellar masses are based on the stellar population synthesis model fitted from FUV to z band photometry (Blanton & Roweis 2007). The redshift distribution of the sample is 0−0.15 and peaks at z ∼ 0.03, while the stellar mass distribution is designed to have a flat distribution between 10⁹ − 10¹¹ M_⊙ (Wake et al. 2017). The IFU cubes have an r-band full width at half maximum (FWHM) of ∼2.5 arcsec and cover 3600−10 300 Å with a spectral resolution of R ∼ 2000 (Smee et al. 2013). We note that the MaNGA sample under-samples at the faint end (r-band absolute magnitude larger than −19, Steffen et al. 2023), but this does not affect the main results in this paper since our radio AGN sample is dominated by massive galaxies (see Sect. 4.1).

All the spatially resolved optical spectra used in this paper were obtained from the MaNGA Data Analysis Pipeline (DAP, Westfall et al. 2019; Belfiore et al. 2019). Emission line, absorption line, and kinematic measurements in this paper were taken either directly from DAP-MAPS images or fitted from the DAP-LOGCUBE spectra.

The global star formation rates are calculated from the attenuation- and AGN-corrected Hα luminosity using the relation suggested by Kennicutt & Evans (2012): SFR (M_⊙ yr⁻¹) = L_HαSF (erg s⁻¹) × 10^−41.30. For dust attenuation, we use the Balmer decrement assuming the Case B recombination (Osterbrock 1989) and the reddening curve from Calzetti et al. (2000): L_Hα = L_Hα, obs × [(Hα/Hβ)_obs/2.86]^2.6. The AGN contribution to L_Hα is estimated by a re-projected optical diagnostic diagram proposed in Ji & Yan (2020), which uses the line ratios [N II]/Hα, [S II]/Hα, and [O III]/Hβ. The empirical relation can be found in Equation (6) of Jin et al. (2021). This step does not affect the star forming galaxy sample used in Sect. 3.1. For RDAGN sample, the AGN contribution correction changes about −0.2 dex (median value) on the SFRs, while their SFRs are about −2.3 dex below the star formation main sequence, thus this correction does not affect the discussions in Sect. 4.1 about the global quenching of RDAGN host galaxies.

We build the spatially resolved star formation histories based on the nonparametric, full spectra fitting results from the Firefly-VAC (Goddard et al. 2017; Neumann et al. 2022). The stellar population models (Maraston et al. 2020) used for fitting are from the MaNGA Stellar library of Milky Way stellar spectra (MaStar, Yan et al. 2019). These models are also used in the DAP and have matched spectral resolution with the galaxy spectra since they are observed by the same instrument. The age grids of MaStar models used in fitting can help us build the mass assembly history of a given spatially resolved region with a precision of at least 1 Gyr.

2.2. LoTSS

The LOFAR Two Meter Sky Survey (LoTSS) is an ongoing deep wide area radio wavelength imaging survey (Shimwell et al. 2017). It aims to observe the whole northern sky at frequencies 120−168 MHz. The expected root-mean-square (RMS) sensitivity and spatial resolution are ∼100 μJy (at favorable declinations) and 6 arcsec, respectively. The latest data release, LoTSS-DR2 (Shimwell et al. 2022), covers a 5634 square degrees sky area. We refer to Shimwell et al. (2022) for details about the data reduction. In summary, LoTSS-DR2 provides well-calibrated 6 arcsec resolution mosaic images. Additionally, it includes a catalog of approximately 5 million radio sources created by the Gaussian-fitting-based radio source extraction code PyBDSF (Mohan & Rafferty 2015).

In this paper, we use the source positions, flux density and spatial size from a value-added catalog with optical host galaxy identification information (hereafter LoTSS-VAC, Hardcastle et al. 2023). The host galaxy identification is based on a likelihood-ratio cross-match method described in Williams et al. (2019), as well as visual classifications from the citizen science project, Radio Galaxy Zoo: LOFAR². The optical and infrared information used for this cross-match is based on the DESI Legacy Imaging Surveys³ (Dey et al. 2019) and the unWISE data (Schlafly et al. 2019). About 85% LoTSS detections have an optical or infrared identification (Hardcastle et al. 2023). The total flux density and morphological parameters have been examined, and where necessary corrected, in the host identification steps. This correction is necessary for deep radio observations in the MaNGA redshift range, because many radio galaxies have large sizes or diffuse emissions, creating large radio-optical offsets which will affect the traditional cross-matching method based on source separations. Furthermore, the erroneous association of a radio galaxy with only one of its multiple components will result in its the radio flux density being significantly underestimated. For the LoTSS-MaNGA matched sample used in this paper, the LoTSS-VAC recovers 44 radio galaxies missed in the raw PyBDSF catalog and corrects the flux density of 157 galaxies, which means that the raw catalog would have lost ∼8% accuracy.

2.3. LoTSS-MaNGA parent sample

Our analysis sample includes all z > 0.01 MaNGA galaxies (to ensure a large enough IFU field of view) with LoTSS observations. Galaxies outside the LoTSS sky coverage or in the bad-quality sky areas are excluded based on the mosaic images, leaving us with 5687 MaNGA galaxies covered by LoTSS DR2. We match the MaNGA galaxy positions with the galaxy positions in LoTSS-VAC (i.e., the positions for the identified host galaxy, not the radio positions) within a radius of 9 arcsec (1.5 × beamsize and smaller than MaNGA IFU)⁴. This process results in 2655 detections (47%) with radio flux density and size measurements. We also visually check all the radio and optical images to exclude 109 (2%) galaxies with severe background source blending or having insufficient MaNGA coverage (IFU coverage is too small). For 2923 (51%) nondetections, we measure the median RMS flux within a 9 arcsec radius aperture around the MaNGA positions in LoTSS images. We note that there are several nondetections close to luminous radio sources (e.g., MaNGA 8312−3704), and they have much higher RMS levels than sources in the normal fields due to artifacts around the strong sources. Thus, we consider the 2% of nondetections with the highest RMS values as outliers and exclude them from our final sample. We adopt 8 times the RMS flux density in one beam (F_RMS) as a secure flux density upper limit for each nondetection, since the catalog source completeness at this cut level reaches more than 95% (Hale et al. 2024).

We calculate the rest-frame 144 MHz luminosity or upper limit using the formula $L_{144\mathrm{MHz}}=F_{\mathrm{LoTSS}}\times 4\pi D_{\mathrm{L}}^2\times {(1+z)}^{(-\alpha -1)}$, where D_L is the luminosity distance, z is the galaxy redshift from MaNGA observation, and α is assumed radio power-law index (−0.7, see Sect. 1). F_LoTSS is the total flux density from LoTSS-VAC for detections, or 8 × F_RMS for nondetections. This rest-frame conversion only corrects ∼2% of the value at MaNGA’s peak redshift z ∼ 0.03, which is smaller than the level of most flux density errors, thus does not affect our analysis and results.

Fig. 1 shows how the radio luminosity detection threshold varies with redshift for our sample of MaNGA galaxies with LoTSS-DR2 overlap. Purple dots represent LoTSS-detected galaxies and grey upper limits are nondetections. The rest-frame 144 MHz luminosity varies between 10²⁰ − 10²⁶ W Hz⁻¹. LoTSS DR2 reaches 90% completeness at ∼0.8 mJy (Shimwell et al. 2022), resulting in a limit of L_144 MHz ∼ 5 × 10²² W Hz⁻¹ at z = 0.15. We note that the exact limiting sensitivity also depends on the sky position, which causes the dispersion for the upper limits in Fig. 1. A few nondetections significantly above the LoTSS limit are the sources affected by strong artifacts in the LoTSS images, which are usually caused by nearby luminous sources (Shimwell et al. 2022). We try to minimize the bias from these sources by taking into account both the upper limit and its error when selecting RDAGN and control samples. This is discussed in detail in Sects. 3.2 and 3.3.

Fig. 1. Radio luminosity and upper limits of LoTSS-MaNGA sample vs their redshift. In the main panel the purple dots show radio luminosities with errors for the LoTSS detections. The gray arrows represent the upper limits of the radio luminosity for LoTSS nondetected galaxies. The solid, dash-dotted, and dashed lines represent the 90% completeness sensitivity limits of LoTSS (Shimwell et al. 2022), FIRST (Becker et al. 1995), and NVSS (Condon et al. 1998), respectively. We note that a few nondetections that are significantly above the LoTSS limit are the sources affected by nearby artifacts. The upper panel shows that the radio detection rate stays approximately constant with redshift, but decreases noticeably at z < 0.03, where the lower mass galaxies start to dominate the MaNGA sample.

The new LoTSS observations significantly improve the radio detection rate. As a comparison, about 72% and 41% of LoTSS detections will fall below the detection limits of NVSS (dashed line in Fig. 1, Condon et al. 1998) and FIRST (dash-dotted line in Fig. 1, Becker et al. 1995), respectively. The LoTSS detection rate of MaNGA galaxies (Fig. 1 upper panel) stays approximately constant with redshift, but decreases noticeably at z < 0.03. This is due to the fact that the selection of MaNGA parent sample will constrain the apparent sizes of galaxies to not over-fill the IFUs, thus at lower redshifts the sample will be dominated by physically smaller, and hence lower mass, galaxies which tend to have faint radio emission.

3. Radio-excess AGN classification

3.1. The L_144 MHz–SFR relation

As mentioned in Sect. 1, synchrotron continuum luminosity is tightly correlated with the star formation rate in SFGs. This relation can be used to select radio-excess AGN (or RDAGN for short). Using LoTSS data, the radio-SFR relation has been investigated based on different SFR tracers and over different spatial scales. For example, Gürkan et al. (2018), Smith et al. (2021), and Best et al. (2023) used SFRs from SED-fitting for SFG populations out to z ∼ 1, while Heesen et al. (2022) used SFRs converted from total infrared luminosity for resolved local group SFGs. Shenoy et al. (in prep.) used both SED-fitted SFRs as well as integrated Hα derived SFRs. Up to now, there has been no relation derived from IFU-based Hα SFRs, despite the higher accuracy of this method. For consistency, here we derive the L_144 MHz–SFR relation from the LoTSS-MaNGA sample to help us select RDAGN later.

The details of how we derive this relation can be found in Appendix A. In brief, we select a pure star-forming galaxy sample without AGN contamination of both the Hα and radio continuum luminosities. Then the L_144 MHz–SFR relation is derived by performing a linear fit to the rest-frame 144 MHz luminosities and SFRs, inspired by the efforts in Gürkan et al. (2018). The intrinsic errors σ_I of this relation are determined using the emcee (Foreman-Mackey et al. 2013) Markov chain Monte Carlo (MCMC) method.

Our derived radio-SFR relation can be expressed as:

$$\begin{array}{c}\hfill \log \left(\frac{L_{144\mathrm{MHz}}}{\mathrm{W}{\mathrm{Hz}}^{-1}}\right)=1.16\times \log \left(\frac{\mathrm{SFR}}{{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}}\right)+21.99\pm 0.23.\end{array}$$(1)

Figure 2 shows that our results generally agree with previous low-redshift relations based on different samples. For comparison, Brown et al. (2017) (blue line) used shallow TGSS sample with SFRs from single fiber Hα; Smith et al. (2021) (cyan line) and Best et al. (2023) (orange line) used a sample with deep radio observations and SFRs from multi-wavelength SED fitting; and Heesen et al. (2022) (green line) used a sample of local galaxies with SFRs measured from total infrared luminosities. Our results suggest an intrinsic scatter of σ_I ∼ 0.23 dex in this relation, consistent with other LoTSS samples (e.g., Smith et al. 2021; Cochrane et al. 2023). The slope is larger than one, which supports the prediction of calorimetry models (e.g., Chi & Wolfendale 1990).

Fig. 2. Correlation between rest-frame 144 MHz luminosity and SFR. The black dots with errors are the pure star-forming galaxies used for linear fitting. These SFGs are selected based on optical line diagnostic and radio image visual classification (see Appendix A for details). The SFRs are integrated from MaNGA Hα maps and are corrected for dust attenuation as mentioned in Sect. 2.1. The red solid line is the MCMC linear fitting result. Intrinsic errors, ±3σI, are shown as red dashed lines. The blue, cyan, green, and orange dashed lines are the linear fitting results in Brown et al. (2017), Smith et al. (2021), Heesen et al. (2022), and Best et al. (2023), respectively. Our result is in good consistency with most literature results.

This relation is found to also depend on the stellar mass (e.g., Smith et al. 2021; Heesen et al. 2022; Shenoy et al., in prep.), but the mass dependence is weak and only significant in low-mass systems (Das et al. 2024) and at higher redshift (Delvecchio et al. 2021), thus we do not include the mass parameter in our fitting. Our fitting procedure also does not include the radio nondetections, which may result in a slightly more conservative cut for radio AGN classification. Though our result is reliable enough to classify bona fide radio AGN in massive galaxies (e.g., M_* > 10^9.5 M_⊙), the mass dependence should be considered when applying it in lower mass sample.

3.2. Radio-excess AGN and their jet luminosities

We classify a MaNGA galaxy as a radio-excess AGN host if its L_144 MHz is at least 3σ level higher than the L_144 MHz predicted by its SFR and Equation (1). Considering the uncertainty of these measurements, the classification can be expressed as follows: log(L_144 MHz − 3σ_L) > 1.16 × log(SFR + 3σ_SFR)+21.99 + 3σ_I. This method results in a new sample of 616 (11%) radio AGN hosts among 5548 LoTSS-MaNGA galaxies. This way the AGN sample should all be secure AGN, but it may be incomplete at lower radio luminosities. Other detections and nondetections are considered as galaxies without secure evidence of radio excess (hereafter referred to as non-RDAGN for convenience).

We also calculate the radio jet luminosity L_AGN by subtracting the star formation contribution from the total radio luminosity: $L_{\mathrm{AGN}}=L_{144\mathrm{MHz}}^{\mathrm{observed}}-L_{144\mathrm{MHz}}^{\mathrm{SFR}}$, where $L_{144\mathrm{MHz}}^{\mathrm{SFR}}$ is the 144 MHz luminosity predicted by Equation (1). The error on the L_AGN, σ_LAGN estimate is derived from the observed uncertainty in $L_{144\mathrm{MHz}}^{\mathrm{observed}}$ and the uncertainty in $L_{144\mathrm{MHz}}^{\mathrm{SFR}}$, and thus includes both σ_SFR and σ_I.

For non-RDAGN, we need an estimate of the upper limit of L_AGN to build a reliable control sample. The SMBHs in non-RDAGN may still have weak radio activities hidden in the host star formation⁵. For non-RDAGN with radio detections, we consider that the radio continuum is dominated by star formation rather than AGN, thus we assume a maximum upper limit of L_AGN + 3σ_LAGN. For non-RDAGN not detected by LoTSS, we adopt the upper limit of radio luminosity (defined in Sect. 2.3) also as the upper limit of L_AGN. By these definitions, the non-RDAGN L_AGN upper limit is the maximum possible luminosity of the hidden AGN. A non-RDAGN with the L_AGN upper limit much lower than the L_AGN of a RDAGN can thus be confidently selected as a control galaxy. We note that the high SFR star-forming galaxies (SFGs) with radio detection will then have higher L_AGN upper limits, because they can have more radio luminosity dominated by star formation than some low SFR RDAGN. The consequences of this selection effect are discussed in Sect. 5.1.

There are two existing MaNGA radio AGN samples in the literature (Comerford et al. 2020; Mulcahey et al. 2022). The key improvements that we have implemented in this work are, (1) classifying AGN based on the self-consistent MaNGA-derived L_144 MHz–SFR relation, (2) including new measurements of L_AGN as well as upper limits in the case of nondetections. Our sample is therefore better suited to select a non-RDAGN control sample with quantified L_AGN (or upper limits) and is more accurate at distinguishing some high SFR, radio-luminous SFGs from RDAGN. Compared to Comerford et al. (2020), where the sample is matched to the NVSS/FIRST selected radio AGN catalog (Best & Heckman 2012), the new deep LoTSS observations that include many more star-forming galaxies allow us to use the more physical L_144 MHz–SFR selection rather than the L_1.4 GHz–D_n4000 classification proposed in Best et al. (2005b). In Mulcahey et al. (2022), the radio AGN sample is selected from LoTSS DR2, but the radio classification uses an empirical relation based on a non-IFU sample (Sabater et al. 2019).

Fig. 3 illustrates our radio AGN classification for the entire LoTSS-MaNGA sample. Detections above the 3σ_I limit are classified as radio-excess AGN (red), while other detections (blue) and nondetections (green) are classified as non-RDAGN with an upper limit for jet luminosity. We also show the radio AGN classified in Best & Heckman (2012) using NVSS and FIRST data (i.e., the Comerford et al. 2020 sample; black stars), demonstrating not only that our sample probes radio AGN to significantly lower luminosities, but also that some high SFR SFGs may be misclassified as AGN using their methods. A total of 77% of radio AGN in our sample are newly classified.

Fig. 3. Radio AGN classification for the whole LoTSS-MaNGA sample. The black solid and dashed lines are the radio-SFR relation (Equation (1)) and the 3σI limit plus the typical L144 MHz error, respectively. Detections with observed L144 MHz larger than the limit are classified as radio-excess AGN (red stars). Other detections (blue) and nondetections (green) are considered as radio normal galaxies (non-RDAGN). We plot the matched radio AGN classified in Best & Heckman (2012) and Comerford et al. (2020) as a comparison (black stars). Our sample includes more low-luminosity radio AGN and can distinguish some high SFR SFGs; and 77% of radio AGN in our sample are newly classified.

3.3. Control sample selection

It is well-established that galaxy age and star formation rate (SFR) are strongly correlated with stellar mass (e.g., Gallazzi et al. 2005; Guo et al. 2015). To evaluate possible different star formation histories of radio AGN host galaxies, it is essential to construct mass-controlled samples for comparison. For each radio AGN host, we create a control sample consisting of non-RDAGN (no confirmed radio excess) that meet the following criteria:

• Jet luminosity: upper-limit$(L_{\mathrm{AGN}}^{\mathrm{control}})<0.2\times L_{\mathrm{AGN}}$,

• Stellar mass: |log($M_{\ast }^{\mathrm{AGN}}/M_{\ast }^{\mathrm{control}})|<0.3$,

• Physical effective radius: |log($R_{\mathrm{e}}^{\mathrm{AGN}}/R_{\mathrm{e}}^{\mathrm{control}})|<0.2$,

• Angular diameter distance: |log($D_A^{\mathrm{AGN}}/D_A^{\mathrm{control}})|<0.2$.

The first criterion ensures that even if the control galaxy contains an unclassified radio AGN, its maximum jet luminosity is still marginal compared to that of the corresponding RDAGN (at most 20%). Other criteria are set to select as many control sample as possible while considering the uncertainty level of the measurements. For instance, the stellar mass calculation results from different methods are found to have a scatter of 0.2−0.3 dex (e.g., Małek et al. 2018; Neumann et al. 2022; Sánchez et al. 2022), the scatters among the R_e fitted from different profiles (e.g., Petrosian and Sérsic) are about 0.1−0.2 dex.

Since radio AGN tend to be more massive than non-RDAGN, there are typically more less-massive galaxies available in the control sample for a given AGN. To maintain a roughly equal median stellar mass between the AGN and control samples, for each AGN we select an equal number of more massive galaxies and less massive galaxies that meet the above criteria. Out of the 616 radio AGN, only three lack any control galaxies that meet these criteria and are discarded.

The left panel of Fig. 4 illustrates the comparison of the stellar masses of each radio AGN (M_*^AGN) and the mean stellar masses of their control samples (M_*^control). A total of 469 (77%) radio AGN have at least four control galaxies. Most radio AGN hosts exhibit matched stellar mass with their control samples. We also plot a similar comparison for the redshift and 90% light radius in Fig. 4, which shows that the RDAGN share consistent properties with their control sample. In the left panel of Fig. 4, the few dots with slightly larger M_*^AGN − M_*^control offsets are the radio AGN hosts with only one control galaxy. These galaxies constitute a small fraction (14%), and their offsets are less than 0.3 dex, thus not significantly impacting the scientific results of this study.

Fig. 4. Left panel: Stellar mass of each radio AGN host, M*AGN, compared to the mean stellar mass of its control sample, M*control. The black dashed line is the 1:1 relation. M*control matches M*AGN for most galaxies. A few galaxies with only one control galaxy have slightly larger offsets, but our scientific results are not affected. The middle and right panels show similar comparisons for the redshift and 90 percent light radius. Both parameters are also well controlled. We note that the matched M* and redshift originate directly from the selection criteria in Sect. 3.3, while the matched R90 is the result of the matched angular distance and physical effective radius.

We do not constrain the control sample to have similar SFR with the RDAGN hosts, because our focus is on the star formation history dating back to the quenching of radio AGN host galaxies. Therefore, SFR is considered a free parameter. Nevertheless, the median difference in SFR between radio AGN hosts and their control sample is only ∼0.23 dex, which is marginal compared to the standard deviation of the radio AGN SFRs (∼1 dex).

4. Results

4.1. Radio AGN prefer massive quiescent hosts

Figure 5 illustrates the global star formation rate of our RDAGN sample. The majority of RDAGN hosts (represented by red stars) are situated below the star-forming galaxy main sequence (SFMS; black solid line). Galaxies below the 3σ_SFMS lower limit of the SFMS are defined as quenched, as they have transitioned from the SFMS to the quiescent population. Notably, the RDAGN hosts are predominantly composed of quenched galaxies (567 of 616, 92%). Furthermore, we investigated the occurrence of RDAGN in the quenched galaxy population. As shown in the upper panel of Fig. 5, the fraction of RDAGN among all quenched galaxies, f_RDAGN, exhibits an increasing trend with stellar mass, reaching at least 50% at M_* ∼ 10^11.5 M_⊙, which indicates increased recurrence of RDAGN in massive galaxies. We note here that this fraction f_RDAGN can only be considered as a lower limit at the LoTSS sensitivity, because LoTSS is still not deep enough to confirm if some weak radio excess exists in these nondetected galaxies. Deeper radio continuum surveys in the future will probe fainter radio AGN.

Fig. 5. SFR vs stellar mass diagram for the LoTSS-MaNGA sample. The background grayscale density map shows the bimodal distribution of the star formation main sequence (SFMS) and the quenched population. Radio AGN are plotted as red stars while galaxies without confirmed radio excess (non-RDAGN) are plotted as blue (detection) and green (nondetection) dots. The black solid line is the SFMS. Galaxies below the 3σSFMS lower limit of the SFMS (black dashed line) are considered as quenched. Radio AGN host galaxies are dominated by massive and quenched populations. The upper panel shows the fraction of radio AGN in all quenched galaxies. The fraction increases with stellar mass to > 50% at the high-mass end (1011.5 M⊙).

These results are consistent with our current knowledge about RDAGN in the local universe. Previous studies have indicated a preference for RDAGN to be hosted in massive, old, and elliptical galaxies (Best et al. 2005a). Recent LOFAR data also suggest that locally, all the M_* > 10¹¹ M_⊙ massive galaxies have a RDAGN inside (Sabater et al. 2019). The fraction observed in our sample is considerably smaller (41%), but is still reasonable considering that we use very different RDAGN classifications and the fraction definition compared to theirs. These results all suggested that the radio AGN hosts currently have old stellar populations and have an early and fast assembly history, which we will investigate in detail in Sect. 4.2.

In Fig. 6, we explore if there is a relationship between global quenching and the power of radio jets. We quantify the degree of quenching by measuring how far galaxies lie below the SFMS, (i.e., smaller, more negative, Δlog(SFR) indicates a more quenched galaxy): Δlog(SFR) = log(SFR_observed)−log(SFR_SFMS). Galaxies below the 3σ threshold of the SFMS (black dashed line) are classified as quenched. The upper panel of Fig. 6 displays the fraction of quenched galaxies among all RDAGN hosts (f_quenched) across different L_AGN bins.

Fig. 6. Quenching level (ΔSFR compared to the main sequence SFR at the same stellar mass) of MaNGA galaxies as a function of their radio AGN luminosities (or the upper limits). The LAGN upper limits of nondetections and are plotted as green arrows. The upper limits of non-RDAGN detections are calculated based on the star formation contribution to the total radio luminosity and are shown as blue arrows. Below the black line (3σ lower than the SFMS) are quenched galaxies. The lack of galaxies in the upper right panel indicates that almost all the hosts of luminous radio AGN are quenched. At LAGN > 1023 W Hz−1, the fraction of radio AGN hosts that are quenched galaxies, fquenched, increases with LAGN from ∼50% to 100%, as shown in the upper panel. The uncertainty shown as the gray region is due to the LAGN limits. Δlog(SFR) shows no correlation with LAGN (Pearson correlation coefficient ρ ∼ −0.11).

Due to the nature of the radio luminosity function of SFGs and radio AGN (e.g., Sabater et al. 2019), the completeness of the RDAGN sample will decrease at low L_AGN, where faint RDAGN in luminous SFGs will become indistinguishable. At L_AGN > 10²⁴ W Hz⁻¹, all the radio sources are confidently classified as RDAGN (see Figure 3), and ∼95% (140 out of 147) of them are below the −3σ limit of the SFMS, indicating that nearly all hosts of luminous RDAGN are quenched. For less luminous RDAGN, the range of f_quenched (shadow region in the upper panel) is determined by considering two extreme assumptions. First, as the lower limit that all L_AGN upper limits indicate that their hosts do not host RDAGN. Second, as the upper limit that all hosts host radio jets with L_AGN below the upper limits. Even under the first assumption, the fraction remains at least ∼50% down to AGN luminosities of ∼10²³ W Hz⁻¹. These fractions support the hypothesis that there are connections between the existence of RDAGN and the quiescence of their hosts. However, the Pearson coefficient between Δlog(SFR) and L_AGN is approximately −0.1, indicating that the quenching level is irrelevant to the radio AGN power.

We note that some radio-detected galaxies have L_AGN upper limits in Fig. 6 (blue arrows). This is because star formation dominates the total radio luminosity in those galaxies, making it hard to classify weak radio AGN based on radio luminosity. For these galaxies, the LoTSS sensitivity is no longer the limitation to confirm if there is radio AGN activity. Higher resolution star formation and radio continuum maps would be needed to decompose the AGN jets from star formation at subgalactic scales.

In summary, quiescent galaxies dominate the RDAGN hosts, with f_quenched increasing with L_AGN while Δlog(SFR) shows no correlation with L_AGN. Radio-mode AGN feedback is believed to be necessary to regulate the cooling flows in dense cluster environment and help massive galaxies maintain quiescent (e.g., Fabian 2012). Radio-mode AGN feedback models have been implemented in cosmological hydrodynamical simulations, where the existence of jet mode feedback is crucial for preventing galaxies from becoming too massive (Torrey et al. 2014). Because the jets dump most of their energy at large distances from the galaxy and because the radio jet power varies on very short timescales compared to the cooling time of the gas in dark matter halos, the luminosity of the jet may not correlate very strongly with the star formation rate in the galaxy. In addition, several observational studies have found weak correlations between the radio jet luminosity and both the kinetic power of the jet and the SMBH accretion rate (e.g., Fabian 2012; Russell et al. 2013), suggesting that the energy transfer from radio jets to surrounding gas is inefficient. These factors may explain why the radio luminosity of our RDAGN sample is irrelevant to the quenching level of the host galaxy.

4.2. Radio AGN host galaxies quenched earlier

In this subsection, we investigate the mass assembly history of the RDAGN host galaxies. It has been widely argued in simulations that negative AGN feedback effects are necessary to halt the growth of massive galaxies (e.g., Kauffmann & Haehnelt 2000; Springel et al. 2005; Schaye et al. 2015; Nelson et al. 2019). If negative feedback is universally present in real galaxies, we should observe AGN-related suppression of star formation in the star formation histories (SFHs) of galaxies.

As mentioned in Sect. 2.1, we employed the full-spectrum fitting method to obtain radially resolved nonparametric SFHs. In Sect. 4.1 we demonstrated that RDAGN hosts are mostly massive quiescent galaxies, with predominantly old stellar populations that are insensitive to recent star formation. To better illustrate the SFHs, we derive the mass accumulation history instead of the SFR history by quantifying the fraction of a galaxy’s stellar mass that was formed at different lookback times. We construct the mass assembly histories for each RDAGN host and their control galaxies in three radial bins: 0−0.5 R_e to represent the unresolved central regions, 0.5−1 R_e for the middle regions, and 1−1.5 R_e for the outskirts. The mass assembly history of individual galaxies can be somewhat uncertain, so we combine them to compare the average difference between radio galaxies and the control sample. The mass assembly histories are normalized at the observed time of t_lookback = 0, corresponding to a median redshift of z ∼ 0.07 (0.9 Gyr before today). For each RDAGN host, the controlled history is the median of its control sample’s SFHs. Then the median of all controlled histories are considered as the final results, and uncertainties are estimated using the error of the median (standard deviation divided by the square root of the sample size).

Fig. 7 shows the comparison of median mass assembly histories between RDAGN hosts and their controls. For both radio AGN hosts (red solid lines) and their controls (blue dashed lines), galaxies form the central region (0−0.5 R_e) earlier than the middle region (0.5−1 R_e) and outskirts (1−1.5 R_e), following an inside-out formation scenario. At all radii, RDAGN built their mass faster and quenched earlier, with the early quenching more pronounced in the central region. Note that this difference does not necessarily imply a direct connection between RDAGN and the evolution of their hosts, because of the very different timescales of radio jet evolution compared to the buildup of the stellar mass. It is possible that different conditions prevail during the bulk of the stellar mass build-up in radio AGN host galaxies and the control galaxies. Faster growth could be a consequence of radio AGN being located in denser environments (e.g., Peng et al. 2010). This scenario is also favored by the more significant impact in the central region; in the standard ΛCDM cosmological paradigm where structure forms hierarchically, the formation time of the inner region of the galaxy is boosted to earlier epochs in dense environments. RDAGN activity may only be triggered after their hosts exceed a certain mass threshold and may only play an important role in maintaining the quiescence at a late stage of the host galaxies’ evolution. Such a scenario is in keeping with the implementation of radio mode feedback in semi-analytical models of galaxy formation (e.g., Croton et al. 2006).

Fig. 7. Mass assembly across the cosmic lookback time for different regions of RDAGN hosts and their controls. The solid red lines are the median SFHs of radio AGN hosts and the blue dashed lines is the median of their control samples. We also plot the average mass assembly history of galaxies with halo masses of 1013 M⊙ from the Illustris simulation (Rodriguez-Gomez et al. 2016). Although the simulation relation is for the total stellar mass, it is kept fixed in all the panels for reference. The panels from left to right are the results for the regions of the whole galaxy, 0−0.5 Re, 0.5−1 Re, and 1−1.5 Re, respectively. RDAGN hosts formed faster than their controls at all radii, especially in the central region. Our spectra-fitted total mass assembly history is reasonably matched by the cosmological simulation which incorporates the radio AGN feedback model.

Our SFHs derived from stellar population synthesis fitting are consistent with some simulations including radio AGN feedback models. In Fig. 7, we plot the mass assembly history of ∼ 10^11.5 M_⊙ galaxies in a 10¹³ M_⊙ halo from the Illustris cosmology simulation (Rodriguez-Gomez et al. 2016). This simulation includes a radio-mode AGN feedback model that is found to be important in suppressing star formation in massive galaxies. These authors also found that strong and weak radio jets create similar impacts on the result (Torrey et al. 2014), which is consistent with the lack of correlation between radio luminosity and the quenching level, as discussed in Sect. 4.1. The mass assembly of the simulated galaxies (black line in Fig. 7) generally agrees with our observational results. The discrepancy between simulated and spectra-derived SFHs at t_lookback > 10 Gyr is due to that the stellar population synthesis is not sensitive to the oldest stars (i.e., the SFHs at early universe).

With the mass accumulation curves, we calculate the lookback time t_X when the galaxy accumulates X percent of the current mass. These lookback times are often used to quantify how quickly galaxies quench or how early galaxies become quenched (e.g., Pacifici et al. 2016). Since our sample predominantly consists of quiescent galaxies with no significant recent star formation, we consider t₉₀ as the lookback time when galaxies cease rapid growth. The difference can be determined through the values from RDAGN hosts and their control sample’s median values, as Δt₉₀. A positive value means that the t₉₀ of RDAGN hosts is larger (i.e., these galaxies formed earlier). We calculate the Δt₉₀ for each RDAGN host compared to its controls, and plot the distributions in Fig. 8. Panels in Fig. 8 from top to bottom are the results in radial bins 0−0.5 R_e, 0.5−1 R_e, and 1−1.5 R_e, respectively. The median values and errors on the median (standard deviation divided by square root of the sample size) are listed in the top left of each panel to indicate the significance of the differences. In Fig. 8, the distribution of Δt₉₀ is shifted towards positive values in all three radial bins, especially in the central one, consistent with the results of mass assembly histories. This indicates that the RDAGN hosts tend to become quenched ∼0.4−0.6 Gyr earlier.

Fig. 8. Distributions of the difference between the t90 from RDAGN and their controls. The three panels from top to bottom are the results for the inner, middle, and outer regions of galaxies. The median difference and its error (standard deviation divided by the square root of the sample size) are noted in the top left of each panel. Larger t90s indicate that RDAGN hosts globally become quenched earlier.

4.3. Emission line excess in the nuclei of radio AGN hosts

Previous analyses have shown relatively marginal differences in star formation properties between RDAGN hosts and their controls. In this subsection, we directly compare the observed spectra of RDAGN hosts and their controls, focusing particularly on the emission lines. Radio jets have long been known to ionize the ISM and create observable emission lines (e.g., Blandford et al. 2019). High-resolution jet simulations show that the kinetic energy in radio jets can be transferred to the surrounding gas by shocks in supersonic outflows (Weinberger et al. 2017).

Most massive quiescent galaxies do not have strong emission line features. To enhance the signal-to-noise ratio (S/N) of the emission lines (i.e., to reduce the noise), we leverage the abundant spectra from IFU observations and stack (median) them in radial bins for both RDAGN hosts and their controls. During the stacking process, all spectra are converted to the rest-frame. Spectra in each radial bin of a galaxy are first stacked by observed flux, which ensures that each galaxy contributes equally to the population stack. This approach highlights differences in emission and absorption lines between RDAGN hosts and controls. In the differential spectra, stronger emission lines in RDAGN hosts appear as emission, while weaker lines appear as absorption, and vice versa for absorption lines.

As an illustration, Fig. 9 shows the stacked differential spectra from the nuclear (within a 2.5″ diameter) regions of the radio AGN and their controls over the wavelength interval 3700 Å−8800 Å. The strong optical emission and absorption lines are indicated in green lines. In general, RDAGN hosts show redder spectra and several emission line excess or deficiency features in their nuclear regions. The most prominent features are the [N II]λ6550, 6585 and Hα emission line excess. This excess is also confirmed by other less prominent lines, such as Hβ and [S II]λ6718, 6733 doublets. For absorption lines, RDAGN hosts have clearly weaker Ca lines (CaHλ3934, CaKλ3968, and Ca IIλ8498, 8542, 8662), appearing as emission-like features. In contrast, the Mg b and Na D indices show absorption-like features, indicating RDAGN hosts have lower Na and Mg abundance. These absorption line features are not strong enough to accurately constrain the metallicity differences, so we only model the [N II] and Hα features in this work.

Fig. 9. Stacked differential spectrum between all RDAGN hosts and their controls from the nuclear regions (∼2.5″ diameter). The spectrum is shown over the rest-frame wavelength range 3700 Å−8800 Å. This wavelength range covers the most prominent emission and absorption lines from [O II]λ3727 to Ca IIλ8662, which are indicated by green dashed lines. The emission-like or absorption-like features indicate that those emission lines are stronger or weaker (respectively) in RDAGN hosts than in their control galaxy stacks. The most prominent features of emission line excess are from [N II]λ6550, 6585 and Hα. This kind of differential stacked spectra are created for different regions and different RDAGN populations, and are analyzed in Sect. 4.3.

To investigate which population of RDAGN are contributing to the stacked emission line excess and deficiency, we divide the RDAGN sample into different subsamples by jet luminosity, radio morphology, and stellar mass, respectively. In each category, we create 16 stacked different spectra from these different subsamples and from four radial bins (0−0.5 R_e, 0.5−1 R_e, 1−1.5 R_e, > 1.5 R_e). We then focus on the wavelength region between 6380 Å and 6670 Å to quantify the width and strength of the Hα + [N II] features. We subtract the continuum by fitting it polynomially and then use the single component Gaussian model to fit the equivalent width and dispersion of these three lines.

In order to ensure that our spectral stacks of radio AGN are fully complete, we use the 90% completeness cut in radio luminosity as a function of redshift to divide the radio AGN sample into four luminosity bins constrained by redshift (see Fig. 10). The median AGN luminosities for the bins are chosen to be 10²², 10²³, 10²⁴, and 10²⁵ W Hz⁻¹. In total, 388 RDAGN are left to be used for the spectral stacking analysis.

Fig. 10. Radio luminosity completeness cut for the RDAGN sample used for spectral stacking. In different radio luminosity bins, the redshift limits of the sample are defined by the 90% completeness cut. The median luminosities and redshifts of each subsample are plotted as black circles.

Figure 11 illustrates this approach when we stack the spectra from RDAGN with different L_AGN. We set luminosity bins to make the median L_AGN values equal to 10²², 10²³, 10²⁴, and 10²⁵ W Hz⁻¹. There are 70, 122, 121, and 75 RDAGN in these four luminosity bins, respectively. Panels from top to bottom represent subsamples of four L_AGN bins, from powerful RDAGN to weak RDAGN. Within each luminosity bin, the differential spectra are plotted for four radial bins (0−0.5 R_e, 0.5−1 R_e, 1−1.5 R_e, > 1.5 R_e) from left to right. We note that the outermost region, r > 1.5 R_e, is based on incomplete samples because of the field of view limit of MaNGA IFUs. The equivalent width and the signal-to-noise ratio of Hα are listed in each panel, where positive and negative values represent stronger or weaker Hα in RDAGN hosts, respectively. It can be seen that stacking can reveal equivalent width differences of less than 0.5 Å.

Fig. 11. Stacked differential spectra between RDAGN hosts and their controls in the [N II] doublets + Hα (green dashed lines) wavelength range. Stacking is done in subsamples with four LAGN bins (powerful to weak from top to bottom) and in four radial bins (center to outskirts from left to right). The red lines are the fitting results for the emission lines. During the fitting, we assume that the three emission lines have individual Gaussian profiles. The equivalent width of the Hα and its signal-to-noise ratio derived from the fitting results are listed in each panel. RDAGN hosts, especially the low-luminosity ones, have an emission line excess in their nuclear regions (leftmost panels). High LAGN RDAGN hosts show weak emission line deficiency in their outskirts (right panels).

We use a similar approach for the other two categories, radio morphology and stellar mass, and summarize the equivalent width measurements in Fig. 12. The three panels from top to bottom show the results when dividing the RDAGN host sample by L_AGN, radio morphology, and stellar mass. The calculation of equivalent width (ΔW_{Hα + [N II]}) is based on the same method as in DAP (Belfiore et al. 2019). For radio morphology, we visually classify the RDAGN into FR I, FR II, diffuse, and compact (or unresolved) morphology based on the 6″ resolution LoTSS images. The number of RDAGN in each morphology category is 81, 77, 38, and 192, respectively. In the lowest panel, we divide the RDAGN hosts into four stellar mass bins with median M_* of 10^10.5, 10^10.8, 10^11.2, and 10^11.5 M_⊙, which include 47, 50, 145, and 146 galaxies, respectively. Fig. 12 shows the radial profiles of the line excess. Among different subsamples, the compact radio sources (< 6″) show strong emission line excess, consistent with the scenario that the compact radio AGN are younger (e.g., Chilufya et al. 2024) and therefore have a greater impact on gas ionization around the central SMBHs. In addition, the lower mass RDAGN (median M_*: 10^10.5 M_⊙) have the strongest Hα + [N II] excess in their nuclear regions. This indicates that the connection between radio AGN activity and optical emission is more common when the RDAGN hosts are less massive. In the outskirts (r > R_e), we find weak emission line deficiency, mainly contributed by resolved RDAGN (FR I, FR II, and diffuse RDAGN). These are likely to be older sources where the jet impacts the ISM at larger distances.

Fig. 12. Equivalent width radial profile and the 3σ error bars of the Hα + [N II] excess or deficiency (ΔWHα + [N II]) in different RDAGN subsamples. The panels from top to bottom show the results when dividing RDAGN by jet luminosity, radio morphology, and stellar mass, respectively. The profiles of the different subsamples are shown in different colors and markers. Positive ΔWHα + [N II] represents that RDAGN have stronger emission lines and vice versa. Low-luminosity ($\overline{L_{\mathrm{AGN}}}={10}^{22}$ W Hz−1), compact (< 6″), and lower mass RDAGN ($\overline{M_{\ast }}={10}^{10.5}{\mathrm{M}}_{\odot }$) contribute most to the nuclear excess. Weak emission line deficiencies are found in the outskirts of the large-size RDAGN population.

Excess emission lines in RDAGN hosts are only found in the innermost inner bins, which indicates that the excess is likely related to the central AGN. On the other hand, the deficiency is found in the outskirts at radii greater than a few kiloparsecs, where the emission lines are not likely to be affected by the central AGN, so this likely indicates a star formation deficiency in the RDAGN hosts. To examine this assumption, we examine the ionization condition of the ISM based on the velocity dispersion of Hα (σ_Hα) and the [N II]/Hα ratio.

In Fig. 13, we plot the distribution of the nuclear [N II]/Hα vs σ_Hα of all MaNGA galaxies using a logarithm scale (background gray contours). Emission lines from H II regions and AGN or shocks are marked by blue and red contours, which are used to classify galaxies according to location in the BPT diagram (e.g., Kauffmann et al. 2003). The red arrow indicates the direction where AGN or shock ionization becomes more dominant. We measure the line ratio and σ_Hα of the emission excess for the whole RDAGN sample in the nuclear bin and the deficiency in the outskirts, and plot them as a green circle and square respectively. We note that the σ_Hα measurements are shown as upper limits, because the stacking procedure will bring an unknown level of line broadening. The central line excess shows a high [N II]/Hα ratio and is located in the region predicated by AGN or shock models (e.g., D’Agostino et al. 2019). This supports the hypothesis that the central RDAGN is heating and ionizing the nuclear gas. However, the excess is weak (≤3 Å) and not correlated with jet luminosity. This suggests that the RDAGN have a weak and complex impact on the ISM properties in the central region of the galaxy. It is worth noting that the resolution of the LoTSS images and MaNGA IFUs is at the kiloparsec level, necessitating future high-resolution observations to confirm whether stronger effects are seen at smaller distances from the SMBH. In addition, high-resolution images are also key to distinguishing whether the emission is from the jet-related shock or from AGN photoionization. Compared to the central region, the deficiency in the outskirts is located much closer to the star-forming region, suggesting that the RDAGN hosts tend to be more quiescent in outer regions.

Fig. 13. Ionization condition diagnostic of the emission line excess found in the nuclear regions of RDAGN hosts (green circle) as well as the deficiency found in the outskirts (green square). The green upper limits on σHα are due to the unknown level of line broadening after stacking the spectra. The background grayscale map is the distribution of the [N II]/Hα vs σHα for the whole MaNGA sample. Star-forming-like emission and AGN or shock-dominated emission regions as classified in the BPT diagram are plotted as blue and red contours, respectively. The red arrow indicates the direction in which the ionization becomes harder due to the AGN or shock contribution. In RDAGN hosts, the line ratio suggests that the nuclear emission line excess is created by AGN or shocks, while weaker star formation may be the reason for the emission line deficiency in the outskirts.

In summary, we found evidence of emission line excess in the nuclei of RDAGN host galaxies. The excess is stronger in lower radio luminosity, more compact (smaller size of radio emission), and lower stellar mass RDAGN. The line excess is likely linked to the interaction of the jet with the gas in the central regions of the host galaxies of the RDAGN. The line ratios suggest AGN photoionization or shock ionization, and the excess is stronger towards the galaxy center. We note that the level of emission line excess in our sample is not correlated with jet power in the range 10²³ − 10²⁴ W Hz⁻¹. In the lowest luminosity bin, stronger excess is seen.

We note that the sample is not large enough to split the sample into a bin of more than one quantity. In order for the radio jet to produce an emission line excess, the host galaxy must contain gas. Our radio-selected sample, however, includes many massive galaxies that lack gas to interact with the radio jets. This may explain the stronger line excess observed in lower mass RDAGN, as lower mass galaxies usually contain more gas (e.g., Saintonge & Catinella 2022). The excess observed in the lowest radio luminosity bin may also result from the fact that these are mainly found in low-mass galaxies.

5. Discussion

In this section, we discuss in more detail some of the main unsolved problems related to the interpretation of the results in this paper and offer some ideas about how these may be solved in the future.

5.1. Selection effects

In studies focusing on the selection and analysis of RDAGN hosts, addressing the potential biases introduced by the nondetections is crucial. Simply defining the non-RDAGN sample based on the absence of radio detection or radio excess may lead to improper control samples. For instance, a higher redshift radio nondetected galaxy might still have higher luminosity jets compared to a lower redshift weak RDAGN. A SFG with an observed L_144 MHz = 10²³ W Hz⁻¹ without significant radio excess might conceal an L_AGN = 10²² W Hz⁻¹ AGN jet drowned in the emission from HII regions. Both of these two galaxies are not appropriate to be selected as part of the control sample of L_AGN = 10²² W Hz⁻¹ RDAGN, but they have been selected to be part of the control sample in most previous studies. For the sensitivity of LoTSS, this effect is only negligible at L_144 MHz > 10²⁴ W Hz⁻¹, as shown in Fig. 6.

In this work, we avoid this bias by calculating the L_AGN upper limits of all the nondetections and non-RDAGN. The upper limits can help us to select control samples that are certain to either contain no jets or, at worst, jets which are much fainter than those they are being controlled against. Nonetheless, this approach may introduce bias against galaxies with higher SFRs. For example, an L_144 MHz = 10²³ W Hz⁻¹ pure SFG without any AGN activity will not be selected as the control of an L_AGN = 10²³ W Hz⁻¹ RDAGN with the same stellar mass. Despite this bias against higher SFRs, our control sample still exhibits SFRs about 0.2 dex higher than the RDAGN hosts, as also demonstrated in Sect. 4.3. In addition, the RDAGN hosts are dominated by massive galaxies, as shown in Fig. 5, and the high SFR SFGs are rare at M_* > 10¹¹ M_⊙. These considerations ensure that our control sample selection is appropriately tailored to mitigate selection biases and provide robust comparisons with RDAGN hosts.

5.2. Disparity between the jet age and quenching time

If the radio-mode AGN feedback can indeed suppress star formation, we would expect to observe a correlation between the time when the jet activity and the quenching occurred. The timing of galaxies becoming quenched can be inferred from their mass assembly history. On the other hand, estimating the age of radio jets is more complex and typically involves two popular methods: spectral age and dynamical age estimation. The spectral age method relies on the shape of the synchrotron energy spectrum. The synchrotron radiation in a fixed magnetic field will have energy losses that scale as ∂E/∂t ∝ ν². So if there is no source driving constant particle acceleration, the spectrum will become steeper because the flux density at high frequencies will decrease faster (e.g., Harwood et al. 2015). The dynamical age is based on the size of a radio jet, because the jet speed can be derived from the environment density (e.g., Falle 1991), then the timescale can be calculated by the jet size and jet speed. We note that these two ages are often observed to have large disparity (e.g., Harwood et al. 2013).

Most of the RDAGN in this study lack sufficient frequency coverage and spatial information to accurately determine their spectral and dynamical jet ages. However, even for extremely old jets with well-calibrated dynamical and spectral ages, the jet ages are less than 1 Gyr (Turner et al. 2018), which is significantly shorter than the quiescent period observed in our RDAGN hosts. As depicted in Fig. 7, RDAGN hosts are observed to have become quenched approximately 5 Gyr ago. If both timescale measurements are accurate within factors of a few, our results suggest that the onset of radio jet activity likely occurred after the host galaxy became quenched, contradicting the scenario where AGN appear and then suppress star formation.

This then implies that the observed RDAGN activity is not responsible for the host quenching. Several recent studies have argued that different AGN populations may represent different evolution phases of AGN, and the radio AGN represent the very late stages of a cycle of AGN activity (e.g., Klindt et al. 2019; Albán et al. 2024). If so, this would explain why the radio jet activities may occur later than the quenching, since the star formation suppression happened in the AGN’s early stages when the radiation feedback is strong. Even with the short timescale of radio jet activity, we still find high RDAGN incidence in the most massive galaxies (> 50%, see Fig. 5), which indicates that radio AGN are common in massive galaxies with a high duty cycle, and the radio mode might be an important ingredient in AGN activities. However, it is still difficult to constrain the duty cycle and evolution timescale of AGN according to the observations.

5.3. The triggering of radio AGN

In Figure 5, we show that the fraction of galaxies that host radio-loud AGN is a very strong function of the stellar mass of the galaxy. In this section, we attempt to pinpoint the underlying cause for why radio AGN activity is more common in massive galaxies by including additional parameters in the analysis, namely SMBH mass, the cosmic environment, and the galaxy morphology. The SMBH mass is estimated by the M_BH − σ_* relation used in McConnell & Ma (2013). The σ_* is fitted from the MaNGA spectra using DAP (Westfall et al. 2019). We adopt the numbers of neighbors within the projected 200 kpc radius and line-of-sight Δv < 500 km s⁻¹ as an approximation for the environmental impact, which are from the MaNGA-GEMA catalog (Argudo-Fernádez et al., in prep.)⁶. We use the Hubble type visually classified by Vázquez-Mata et al. (2022) and divided galaxies into elliptical, S0, and late-type (spiral and irregular galaxies).

We calculate the incidence of RDAGN as a function of stellar mass, and plot the results for subsamples divided by the 3 additional parameters in Fig. 14. In the upper panel, the RDAGN incidence is clearly higher in high σ_* systems at fixed stellar mass. The middle and lower panels indicate that at fixed stellar mass, RDAGN prefer denser cosmic environment and early-type galaxies, which is consistent with previous findings (e.g., Best et al. 2005a). The trends as a function of σ_* at fixed stellar mass are by far the strongest ones, suggesting that black hole mass is an important parameter, after stellar mass, that controls the fraction of galaxies with radio mode AGN activity. We note that at a similar Eddington ratio, the more massive black holes will have higher luminosity and thus are more likely to be detected. Sabater et al. (2019) used a larger single-fiber observed sample and showed that the RDAGN incidence dependence on M_BH is marginal. This inconsistency is partly due to that they made a radio luminosity cut on the radio AGN sample. After applying a similar cut (L_150 MHz > 10^22.5 W Hz⁻¹), our fraction values generally agree with the results in Sabater et al. (2019) (f_RDAGN ∼ 5%−10% at M_* ∼ 10¹¹ M_⊙). However, the dependence of the RDAGN incidence on M_BH still exists. Possible reasons could be our much smaller sample size, and/or the different redshift range (z < 0.15 vs z < 0.3).

Fig. 14. Incidence of RDAGN as a function of stellar mass, shown for different galaxy populations. From top to bottom, the galaxies are divided by stellar velocity dispersion (an indicator of the black hole mass), the number of neighbors within 200 kpc (an indicator of the cosmic environment), and Hubble type (galaxy morphology). RDAGN tend to be located in host galaxies with larger stellar velocity dispersions, denser environments, and early-type morphologies. We find that most galaxies with the largest stellar dispersion (σ* > 264 km s−1) host RDAGN, which suggests that the most massive SMBHs are always accreting in the radio mode.

5.4. Implications for the AGN feedback scenario

Combining all the observation results found in this paper, we discuss some possible new implications for the AGN feedback scenario. We note that this discussion focuses on the moderate and faint radio AGN, where a direct connection between the radio AGN and the emission-line properties of the ISM is clearly seen. The more luminous radio jets may dump most of their energy at large distances from the galaxy, and the link to quenching effects may not be as apparent.

Because radio AGN are located in denser environments, at fixed mass, they are likely to have assembled earlier and faster. The observed higher stellar dispersion suggests that their central SMBHs are more massive, indicating possible stronger AGN activity in the past. This past AGN activity may have played an important role in quenching the host galaxies through gas heating, fast consumption of the gas, or gas outflows (e.g., Silk et al. 2024). The early-type morphology and low angular momentum (Zheng et al. 2023) indicates past galaxy-galaxy merger events, which may trigger the past AGN activity and also the fast mass growth, as well as the morphological transition from disk to spheroid (e.g., Rothberg & Joseph 2004). This is also consistent with some simulation results which suggest that episodic AGN events, rather than one time AGN event, contribute to the final quenching of the host galaxies (e.g., Piotrowska et al. 2022; Harrison & Ramos Almeida 2024).

At a later stage, the galaxy becomes quiescent due to the lack of cold gas fueling. This also prevents the central SMBH from accreting efficiently and radiatively, and the current most massive SMBHs tend to appear as radio-mode AGN. In our sample, the RDAGN hosts are dominated by quiescent galaxies and show weaker star formation activities in the outskirts, which suggests that the radio AGN may help the host galaxies to maintain quiescence. The AGN- or shock like emission excess found in the nuclear regions of RDAGN hosts indicates that radio AGN ionize and heat the surrounding gas. One remaining mystery is exactly how the weak and compact jets associated with the central excess phenomenon can maintain the quiescence of the host galaxies globally.

Finally, although a large-scale connection between radio AGN and the heating of surrounding gas has been observed in many central galaxies of massive clusters (e.g., Fabian 2012), we still need more observational evidence of the jet-gas connection in faint radio AGN on galaxy scales to obtain a full view of such physical processes in different galaxy environments and populations.

5.5. Future directions

This work is based on the radio images and optical IFUs with spatial resolutions on the order of arcseconds. In the future, new telescopes and the update of current facilities can carry out large radio sky surveys with subarcsec resolution, for example, the International LOFAR Two-metre Sky Survey⁷, the Next Generation Very Large Array (ngVLA, Di Francesco et al. 2019), and the Square Kilometre Array (SKA)⁸. At low-redshifts, they can resolve the jet features on physical scales of a few hundred of parsecs. As discussed in Sect. 4.1, the higher resolution is also essential to decompose the faint radio population in AGN and star-forming galaxies. This can help us to understand the radio AGN luminosity function in the faint end. The wide frequency coverage of different instruments can also help us build the full synchrotron SED of radio AGN, which is essential for constraining jet ages.

Meanwhile, there are more and more deep and high-resolution IFU samples. Future high-resolution IFUs can provide emission line maps around the AGN at parsec-scale resolution. Combined with high-resolution radio images, the morphology of the radio continuum and emission lines may reveal the detailed process of how the jets from AGN interact with the ISM. In addition to advancements in radio and optical observations, future large-scale surveys at other wavelengths, particularly in the X-ray and infrared, are expected to match the resolution and sensitivity achieved at optical and radio bands. This multi-wavelength approach will provide comprehensive information about the different component structures of AGN, allowing for a more complete understanding of the structure and the evolution of AGN.

6. Summary

We built the largest optical IFU – radio continuum sample to date consisting of 5548 galaxies by matching the MaNGA survey with LoTSS-DR2 observations. We revisited the tight linear SFR – radio luminosity relation using IFU data and find an intrinsic error of ∼0.23 dex. Using this relation, we identified 616 radio AGN with excessive radio emission, and calculated their jet luminosities. We analyzed their global star formation rates and their correlation with jet luminosity. We then derived the radially resolved mass assembly history by fitting the IFU spectra in radial bins and comparing the difference between radio AGN hosts and their control samples. We also stacked the spectra in different radial bins to investigate the impact of radio AGN on the emission lines of their host galaxies. Our results can be summarized as follows:

1. The linear relation in the SFR and L_144 MHz found by LoTSS-MaNGA sample is tight and can be described by Equation (1): $\log (\frac{L_{144\mathrm{MHz}}}{\mathrm{W}{\mathrm{Hz}}^{-1}})=1.16\times \log (\frac{\mathrm{SFR}}{{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}})+21.99$, which is in very good agreement with previous results based on non-IFU samples.

2. At least 41% M_* > 10¹¹ M_⊙ galaxies are found to host a radio AGN. Radio AGN host galaxies are dominated by massive quiescent galaxies. About 92% RDAGN hosts lie below the −3σ of star formation main sequence. The quenching level, ΔSFR, is not related to the jet luminosity (Pearson correlation coefficient ρ ∼ −0.1).

3. The mass assembly history of radio AGN hosts follows an inside-out scenario, similar to normal massive galaxies. Compared to a mass-matched control sample, radio AGN hosts grow faster at all radii, especially in their central regions. They grew to 90% of current mass at lookback times of 4−6 Gyr, and 0.4−0.6 Gyr earlier than their controls. The mass assembly history reconstructed from the optical spectra is consistent with that in simulations incorporating the radio-mode AGN feedback model.

4. Stacked spectra comparisons reveal emission line excess exclusively in the central region of radio AGN hosts, especially in low-mass and radio-compact RDAGN. The high ratio of the excessive [N II] and Hα supports the scenario that the emission line excess is linked to AGN or shock ionization, but it remains to be understood why the excess is not dependent upon jet luminosity.

Considering that the timescale of jet activities is usually less than 1 Gyr, the quenching of their host galaxies that happened several gigayears ago was likely not due to the current RDAGN. However, episodic cycles of radio AGN may play an important role in maintaining the global quiescence of the host galaxies and preventing the most massive galaxies from further growth. This conjecture is supported by the large fraction of quiescent galaxies in RDAGN hosts, as well as the high incidence of RDAGN in massive galaxies. The lack of a clear correlation between jet luminosity and feedback features suggests that energy transfer from radio jets to the ISM is an inefficient and complex process. Exactly how radio jets affect the star formation and interact with the gas in galaxies is yet to be fully understood. Future studies with larger samples and higher spatial resolution observations at different wavelengths are expected to shed more light on the detailed physical processes underlying radio-mode AGN feedback.

Acknowledgments

The authors would like to thank the anonymous referees for their constructive comments and suggestions. We also thank Fabrizio Arrigoni Battaia, Iker Millan Irigoyen, and Huiying Wang for their insightful discussions. PNB is grateful for support from the UK STFC via grants ST/V000594/1 and ST/Y000951/1. SS is grateful for support from the UK Science and Technology Facilities Council (STFC) via grant ST/X508408/1. 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. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss4.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. 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 that 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. This research made use of the Dutch national e-infrastructure with support of the SURF Cooperative (e-infra 180169) and the LOFAR e-infra group. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC). This research made use of the University of Hertfordshire high-performance computing facility and the LOFAR-UK computing facility located at the University of Hertfordshire and supported by STFC [ST/P000096/1], and of the Italian LOFAR IT computing infrastructure supported and operated by INAF, and by the Physics Department of Turin university (under an agreement with Consorzio Interuniversitario per la Fisica Spaziale) at the C3S Supercomputing Centre, Italy. Part of this work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number Ts17/2–1.

References

Appendix A: Fitting the L_144MHz-SFR relation

To build the L_144MHz-SFR relation, the key is to select a sample of pure star-forming galaxies with reliable SFR measurement and without radio jet contribution.

Using the integral field spectra from MaNGA, the SFR can be derived from the Hα line, and corrected for dust attenuation on subgalactic scales. There are galaxies in which the Hα is not dominated by star formation but instead is possibly contaminated by narrow-line AGN or diffuse ionized gas. We exclude these galaxies using the widely used emission line diagnostic (BPT diagram, Baldwin et al. 1981 based on the [N II]λ6585/Hα and [O III]λ5008/Hβ line ratios measured from the central 2.5″ MaNGA spaxels. Specifically, our pure-SFG sample only includes galaxies locating in the star-forming region of the BPT diagram, using the separation diagnostic derived by (Kauffmann et al. 2003), and the galaxies should also have nuclear Hα equivalent widths larger than 3 Å, which ensures the SFR calculation is not affected by the diffused gas from evolved stars (Cid Fernandes et al. 2011).

Then we visually check the LoTSS image of these galaxies and exclude those showing jet-like features and require the signal-to-noise ratio of radio detection larger than 5. We note that this visual classification cannot exclude unresolved radio AGN. However, they are mostly quiescent galaxies in low-redshift (e.g., Fig. 5) and thus will be excluded from the SFG sample by the Hα equivalent width cut. Considering that there are only very few unresolved radio AGN in the final SFG sample (e.g., some outliers in Fig. 5), the radio AGN contamination for the fitting should be marginal. These unresolved radio AGN could be one source of the intrinsic scatter of the relation and could make the relation slightly overestimated. This way we get a sample of 906 secure SFGs with reliable SFR and radio luminosity measurements.

Fig. A.1 shows the emission line selection of this SFG sample. ΔBPT is the distance to the division line in the BPT diagram (Equation 1 in Kauffmann et al. 2003), and the negative values represent that galaxies locate in the star-forming region. As a result of all the constrains, this SFG sample has a SFR limit about 0.1 M_⊙ yr⁻¹ (see Fig. 2), thus our result about L_144MHz-SFR relation may not be applicable for SFGs with SFRs lower than this value.

Fig. A.1. SFG sample selected by the emission line diagnostic and radio morphology. The y-axis is the equivalent width of Hα in the nuclear region measured from MaNGA spectra. The x-axis is the distance to the division line in the BPT diagram (Equation 1 in Kauffmann et al. 2003), based on the [N II]λ6585/Hα and [O III]λ5008/Hβ ratios. Negative and positive ΔBPT mean that galaxies locate in the star-forming region and AGN region in the BPT diagram, respectively. This Hα-ΔBPT diagram can classify the galaxies into SFG, AGN, and quiescent galaxies. Our pure SFG sample is plotted as stars and color-coded by their radio luminosities. The background grayscale contour is the distribution of the whole MaNGA sample.

We set three parameters in the linear fitting of the L_144MHz-SFR relation. In addition to the slope and constant, we introduce an intrinsic scatter parameter, which should be the scatter level of the L_144MHz-SFR relation after excluding the effect from measurement errors. Our fitting follows the method in Hogg et al. (2010). The Markov chain Monte Carlo (MCMC) sampling result is shown in Fig. A.2. The intrinsic error is found to be small, about 0.23 dex, which indicates that the radio emission is a good tracer for galaxy global SFR. The relation is listed in Equation 1. At the 3σ level, our result is consistent with previous results based on different SFR tracers and radio data from different frequencies (Wang et al. 2016; Brown et al. 2017; Gürkan et al. 2018; Smith et al. 2021; Heesen et al. 2022; Best et al. 2023).

Fig. A.2. Markov chain Monte Carlo sampling of our fitting for the L144MHz-SFR relation. The three parameters are the slope, the constant value, and the intrinsic scatter in logarithm scale.

All Figures

Fig. 1. Radio luminosity and upper limits of LoTSS-MaNGA sample vs their redshift. In the main panel the purple dots show radio luminosities with errors for the LoTSS detections. The gray arrows represent the upper limits of the radio luminosity for LoTSS nondetected galaxies. The solid, dash-dotted, and dashed lines represent the 90% completeness sensitivity limits of LoTSS (Shimwell et al. 2022), FIRST (Becker et al. 1995), and NVSS (Condon et al. 1998), respectively. We note that a few nondetections that are significantly above the LoTSS limit are the sources affected by nearby artifacts. The upper panel shows that the radio detection rate stays approximately constant with redshift, but decreases noticeably at z < 0.03, where the lower mass galaxies start to dominate the MaNGA sample.

In the text

Fig. 2. Correlation between rest-frame 144 MHz luminosity and SFR. The black dots with errors are the pure star-forming galaxies used for linear fitting. These SFGs are selected based on optical line diagnostic and radio image visual classification (see Appendix A for details). The SFRs are integrated from MaNGA Hα maps and are corrected for dust attenuation as mentioned in Sect. 2.1. The red solid line is the MCMC linear fitting result. Intrinsic errors, ±3σI, are shown as red dashed lines. The blue, cyan, green, and orange dashed lines are the linear fitting results in Brown et al. (2017), Smith et al. (2021), Heesen et al. (2022), and Best et al. (2023), respectively. Our result is in good consistency with most literature results.

In the text

Fig. 3. Radio AGN classification for the whole LoTSS-MaNGA sample. The black solid and dashed lines are the radio-SFR relation (Equation (1)) and the 3σI limit plus the typical L144 MHz error, respectively. Detections with observed L144 MHz larger than the limit are classified as radio-excess AGN (red stars). Other detections (blue) and nondetections (green) are considered as radio normal galaxies (non-RDAGN). We plot the matched radio AGN classified in Best & Heckman (2012) and Comerford et al. (2020) as a comparison (black stars). Our sample includes more low-luminosity radio AGN and can distinguish some high SFR SFGs; and 77% of radio AGN in our sample are newly classified.

In the text

Fig. 4. Left panel: Stellar mass of each radio AGN host, M*AGN, compared to the mean stellar mass of its control sample, M*control. The black dashed line is the 1:1 relation. M*control matches M*AGN for most galaxies. A few galaxies with only one control galaxy have slightly larger offsets, but our scientific results are not affected. The middle and right panels show similar comparisons for the redshift and 90 percent light radius. Both parameters are also well controlled. We note that the matched M* and redshift originate directly from the selection criteria in Sect. 3.3, while the matched R90 is the result of the matched angular distance and physical effective radius.

In the text

Fig. 5. SFR vs stellar mass diagram for the LoTSS-MaNGA sample. The background grayscale density map shows the bimodal distribution of the star formation main sequence (SFMS) and the quenched population. Radio AGN are plotted as red stars while galaxies without confirmed radio excess (non-RDAGN) are plotted as blue (detection) and green (nondetection) dots. The black solid line is the SFMS. Galaxies below the 3σSFMS lower limit of the SFMS (black dashed line) are considered as quenched. Radio AGN host galaxies are dominated by massive and quenched populations. The upper panel shows the fraction of radio AGN in all quenched galaxies. The fraction increases with stellar mass to > 50% at the high-mass end (1011.5 M⊙).

In the text

Fig. 6. Quenching level (ΔSFR compared to the main sequence SFR at the same stellar mass) of MaNGA galaxies as a function of their radio AGN luminosities (or the upper limits). The LAGN upper limits of nondetections and are plotted as green arrows. The upper limits of non-RDAGN detections are calculated based on the star formation contribution to the total radio luminosity and are shown as blue arrows. Below the black line (3σ lower than the SFMS) are quenched galaxies. The lack of galaxies in the upper right panel indicates that almost all the hosts of luminous radio AGN are quenched. At LAGN > 1023 W Hz−1, the fraction of radio AGN hosts that are quenched galaxies, fquenched, increases with LAGN from ∼50% to 100%, as shown in the upper panel. The uncertainty shown as the gray region is due to the LAGN limits. Δlog(SFR) shows no correlation with LAGN (Pearson correlation coefficient ρ ∼ −0.11).

In the text

Fig. 7. Mass assembly across the cosmic lookback time for different regions of RDAGN hosts and their controls. The solid red lines are the median SFHs of radio AGN hosts and the blue dashed lines is the median of their control samples. We also plot the average mass assembly history of galaxies with halo masses of 1013 M⊙ from the Illustris simulation (Rodriguez-Gomez et al. 2016). Although the simulation relation is for the total stellar mass, it is kept fixed in all the panels for reference. The panels from left to right are the results for the regions of the whole galaxy, 0−0.5 Re, 0.5−1 Re, and 1−1.5 Re, respectively. RDAGN hosts formed faster than their controls at all radii, especially in the central region. Our spectra-fitted total mass assembly history is reasonably matched by the cosmological simulation which incorporates the radio AGN feedback model.

In the text

Fig. 8. Distributions of the difference between the t90 from RDAGN and their controls. The three panels from top to bottom are the results for the inner, middle, and outer regions of galaxies. The median difference and its error (standard deviation divided by the square root of the sample size) are noted in the top left of each panel. Larger t90s indicate that RDAGN hosts globally become quenched earlier.

In the text

Fig. 9. Stacked differential spectrum between all RDAGN hosts and their controls from the nuclear regions (∼2.5″ diameter). The spectrum is shown over the rest-frame wavelength range 3700 Å−8800 Å. This wavelength range covers the most prominent emission and absorption lines from [O II]λ3727 to Ca IIλ8662, which are indicated by green dashed lines. The emission-like or absorption-like features indicate that those emission lines are stronger or weaker (respectively) in RDAGN hosts than in their control galaxy stacks. The most prominent features of emission line excess are from [N II]λ6550, 6585 and Hα. This kind of differential stacked spectra are created for different regions and different RDAGN populations, and are analyzed in Sect. 4.3.

In the text

	Fig. 10. Radio luminosity completeness cut for the RDAGN sample used for spectral stacking. In different radio luminosity bins, the redshift limits of the sample are defined by the 90% completeness cut. The median luminosities and redshifts of each subsample are plotted as black circles.
In the text

Fig. 11. Stacked differential spectra between RDAGN hosts and their controls in the [N II] doublets + Hα (green dashed lines) wavelength range. Stacking is done in subsamples with four LAGN bins (powerful to weak from top to bottom) and in four radial bins (center to outskirts from left to right). The red lines are the fitting results for the emission lines. During the fitting, we assume that the three emission lines have individual Gaussian profiles. The equivalent width of the Hα and its signal-to-noise ratio derived from the fitting results are listed in each panel. RDAGN hosts, especially the low-luminosity ones, have an emission line excess in their nuclear regions (leftmost panels). High LAGN RDAGN hosts show weak emission line deficiency in their outskirts (right panels).

In the text

Fig. 12. Equivalent width radial profile and the 3σ error bars of the Hα + [N II] excess or deficiency (ΔWHα + [N II]) in different RDAGN subsamples. The panels from top to bottom show the results when dividing RDAGN by jet luminosity, radio morphology, and stellar mass, respectively. The profiles of the different subsamples are shown in different colors and markers. Positive ΔWHα + [N II] represents that RDAGN have stronger emission lines and vice versa. Low-luminosity ($\overline{L_{\mathrm{AGN}}}={10}^{22}$ W Hz−1), compact (< 6″), and lower mass RDAGN ($\overline{M_{\ast }}={10}^{10.5}{\mathrm{M}}_{\odot }$) contribute most to the nuclear excess. Weak emission line deficiencies are found in the outskirts of the large-size RDAGN population.

In the text

Fig. 13. Ionization condition diagnostic of the emission line excess found in the nuclear regions of RDAGN hosts (green circle) as well as the deficiency found in the outskirts (green square). The green upper limits on σHα are due to the unknown level of line broadening after stacking the spectra. The background grayscale map is the distribution of the [N II]/Hα vs σHα for the whole MaNGA sample. Star-forming-like emission and AGN or shock-dominated emission regions as classified in the BPT diagram are plotted as blue and red contours, respectively. The red arrow indicates the direction in which the ionization becomes harder due to the AGN or shock contribution. In RDAGN hosts, the line ratio suggests that the nuclear emission line excess is created by AGN or shocks, while weaker star formation may be the reason for the emission line deficiency in the outskirts.

In the text

Fig. 14. Incidence of RDAGN as a function of stellar mass, shown for different galaxy populations. From top to bottom, the galaxies are divided by stellar velocity dispersion (an indicator of the black hole mass), the number of neighbors within 200 kpc (an indicator of the cosmic environment), and Hubble type (galaxy morphology). RDAGN tend to be located in host galaxies with larger stellar velocity dispersions, denser environments, and early-type morphologies. We find that most galaxies with the largest stellar dispersion (σ* > 264 km s−1) host RDAGN, which suggests that the most massive SMBHs are always accreting in the radio mode.

In the text

Fig. A.1. SFG sample selected by the emission line diagnostic and radio morphology. The y-axis is the equivalent width of Hα in the nuclear region measured from MaNGA spectra. The x-axis is the distance to the division line in the BPT diagram (Equation 1 in Kauffmann et al. 2003), based on the [N II]λ6585/Hα and [O III]λ5008/Hβ ratios. Negative and positive ΔBPT mean that galaxies locate in the star-forming region and AGN region in the BPT diagram, respectively. This Hα-ΔBPT diagram can classify the galaxies into SFG, AGN, and quiescent galaxies. Our pure SFG sample is plotted as stars and color-coded by their radio luminosities. The background grayscale contour is the distribution of the whole MaNGA sample.

In the text

	Fig. A.2. Markov chain Monte Carlo sampling of our fitting for the L144MHz-SFR relation. The three parameters are the slope, the constant value, and the intrinsic scatter in logarithm scale.
In the text