© The Authors 2026

1. Introduction

Star-forming galaxies show a positive relation between total infrared (IR) luminosity (L_IR, rest-frame 8−1000 μm) and monochromatic radio luminosity (e.g., at 1.4 GHz; L_1.4), which remains remarkably tight (σ ∼ 0.26 dex; e.g., Bell 2003) and linear for over three orders of magnitude in both L_IR and L_1.4 GHz (e.g., Helou & Rowan-Robinson 1985; Condon 1992, even though a mild nonlinearity was also found in low-z galaxies, e.g., at z < 0.2 Molnár et al. 2021). This is generally known as the infrared-radio correlation (IRRC). The origin of the IRRC arises from the relationship between radio luminosity and (obscured) star formation rate (SFR). Broadly speaking, the IR emission arises from the light of massive OB stars, which is absorbed and re-emitted by dust (Madau & Dickinson 2014) and therefore directly traces the galaxy’s SFR. On the other side, radio-continuum emission at GHz frequencies originates mainly from synchrotron processes, produced by relativistic cosmic-ray electrons (CRe) accelerated by shock waves produced when massive stars (> 8 M_⊙) explode as supernovae, in the case of star-forming galaxies (SFG).

The logarithmic ratio between the two luminosities, labelled q_IR, is defined as

$$\begin{array}{c}\hfill q_{\mathrm{IR}}=log\left(\frac{L_{\mathrm{IR}}[\mathrm{W}]}{3.75\times {10}^{12}[\mathrm{Hz}]}\right)-log(L_{1.4\mathrm{GHz}}[\mathrm{W}/\mathrm{Hz}]),\end{array}$$(1)

where 3.75 × 10¹² Hz represents the central frequency over the far-infrared (e.g., 42−122 μm) regime. However, the nature of the cosmic evolution of the IRRC has long been debated (e.g., Harwit & Pacini 1975; Rickard & Harvey 1984; de Jong et al. 1985; Helou & Rowan-Robinson 1985; Hummel et al. 1988; Condon 1992; Garrett 2002; Appleton et al. 2004; Jarvis et al. 2010; Sargent et al. 2010; Ivison et al. 2010a,b; Bourne et al. 2011; Magnelli et al. 2015; Delhaize et al. 2017; Gürkan et al. 2018; Molnár et al. 2018; Algera et al. 2020; Basu et al. 2015; Rivera et al. 2017).

Several studies (Magnelli et al. 2015; Rivera et al. 2017; Delhaize et al. 2017) observed a dependence of the IRRC on the redshift, in the form of q_IR ∝ (1 + z)^δ with −0.1 < δ < −0.2. Other works (Gürkan et al. 2018; Smith et al. 2021; Delvecchio et al. 2021, D21 hereafter) have recently argued that the IRRC evolves primarily with M_★, with more massive galaxies showing a systematically lower q_IR, while the dependence on the redshift is secondary and weaker, as q_IR ∝ (1 + z)^{−0.023 ± 0.008} at a fixed stellar mass (see also Gürkan et al. 2018; Smith et al. 2021). Among the possible physical drivers of the q_IR decline observed in massive star-forming galaxies, the energy loss of cosmic-ray electrons (CRes) has been suggested to play a significant role (e.g., Schober et al. 2023; Yoon 2024). Specifically, Yoon (2024) modeled the CRe energy-loss as a function of the gas density and redshift of the host galaxy. Their findings indicate an anticorrelation between these quantities; specifically, galaxies with higher gas densities (or, equivalently, higher stellar masses; e.g., Kennicutt 1998) experience less intense energy losses, arguably because they retain their CRes more easily. On the contrary, less massive galaxies are generally most affected by CRe energy losses. As a consequence of these synchrotron energy losses, less massive galaxies have lower radio luminosities with respect to more massive ones.

One yet unexplored explanation for the M_★ evolution is that a hidden and widespread contribution from radio AGN activity might be driving the decline of the q_IR, especially at high stellar masses (log(M_★/M_⊙)≥10). One very powerful tool to assess the nature of radio emission is the brightness temperature of the emitter (T_b), defined as

$$\begin{array}{c}\hfill T_b[\mathrm{K}]=1.22\times {10}^6\times (1+z)\times \left(\frac{S_{\mathrm{obs}}}{\mathrm{\mu }\mathrm{Jy}}\right){\left(\frac{\nu }{\mathrm{GHz}}\right)}^{-2}{\left(\frac{{\theta }_x·{\theta }_y}{{\mathrm{mas}}^2}\right)}^{-1},\end{array}$$(2)

where S​_obs is the observed density flux at the frequency ν, and θ^x and θ^y are the major and minor axes of the beam (Ulvestad et al. 2005).

Traditionally, high T_b measurements (T_b ≥ 10⁵ K) require very long baseline interferometry (VLBI), which reaches milliarcsecond resolutions at GHz frequencies (e.g., Middelberg et al. 2010, 2013; Rampadarath et al. 2015; Herrera-Ruiz et al. 2017, 2018; Radcliffe et al. 2018, 2021; Njeri et al. 2023; Deane et al. 2024). Recent studies performed with the Very Long Baseline Array (VLBA) at a few milliarcsecond (mas) resolution (e.g., 10 mas ∼ 80 pc at z = 2) have identified AGNs even in galaxies whose emission is dominated by star formation at Very Large Array (VLA) scales (Herrera-Ruiz et al. 2017, 2018; Maini 2017). However, these pioneering studies have been highly biased toward radio-excess AGNs (i.e., showing excess radio emission relative to that expected from pure star formation; Smolčić et al. 2017). Our project AGN-sCAN extends this type of analysis by targeting with the VLBA a representative sample of SFGs without strong VLA radio-excess, down to a sensitivity reaching as deep as σ ≈ 5 μJy/beam at 1.4 GHz at high angular resolution (average beam size: 26 × 7 mas²). From Eq. (2), the T_b limit at > 5σ for an unresolved source at z = 0.5 is T_b ≥ 10⁵ K¹. Compact radio cores (∼mas scales) characterized by such T_b would unambiguously arise from an AGN (Condon 1992; Radcliffe et al. 2018; Morabito et al. 2022).

Section 2 describes the main physical properties of the radio targets composing the sample used in this work. Section 3 summarizes the AGN-sCAN observational strategy and data products, of which a complete and detailed description will be presented by Delvecchio et al. (in prep., D26 hereafter). Section 5 summarizes our main result: the contribution of AGN emission to the nuclear emission of SFGs, classified through the IRRC, is negligible.

We adopted a Chabrier (2003) initial mass function in the mass range 0.1−100 M_⊙. We assumed a standard ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹. We defined the radio spectral index, α, with the following convention: S​_ν ∝ ν^−α, where S​_ν is the flux density at a frequency, ν.

2. Sample selection and physical properties

This VLBA program (ID: VLBA/20B-140, PI: Delvecchio) has targeted 500 massive (M_★ ≥ 10¹⁰ M_⊙) star-forming galaxies in the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007), which is an equatorial 2 deg² field observed with all major space- and ground-based telescopes². We refer the reader to D26 for the detailed description of the survey.

The parent galaxy sample was retrieved from the COSMOS-Web survey (Casey et al. 2023), containing deep imaging in four JWST/NIRCam filters (F115W, F150W, F277W, F444W) over 0.54 deg² (Franco et al. 2025), and with JWST/MIRI 7.7 μm over 0.2 deg² (Harish et al. 2025). Given that our VLBA pointing lies within the COSMOS-Web NIRCam footprint, we updated all redshifts and counterpart properties from the latest COSMOS-Web photometric catalogue (Shuntov et al. 2025), which extends in wavelength out to Spitzer-IRAC bands. Photometric redshifts are in very good agreement with spectroscopic measurements (i.e., median absolute deviation σ_MAD ∼ 0.012 at m_F444W < 28; see Khostovan et al. 2026). We adopted the spectroscopic measurement if available and reliable (> 95% confidence level) based on the compilation by Khostovan et al. (2026).

In order to retrieve a q_IR measurement for each target, we require our galaxies to have been detected (≥5σ) at both IR and radio frequencies. The combined IR dataset comes from de-blended far-infrared/sub-mm data from Herschel (100−500 μm), SCUBA2 (450, 850 μm), AzTEC-1.1 mm, MAMBO-1.2 mm (Jin et al. 2018), and ALMA from the recent A³COSMOS survey (Liu et al. 2019; Adscheid et al. 2024), whenever available. For each source, the total IR luminosity (8−1000 μm, L_IR) is obtained via spectral energy distribution (SED) fitting, and corrected for potential IR-AGN emission (Jin et al. 2018). Radio counterparts are jointly taken from MIGHTEE 1.28 GHz (6″ resolution, Jarvis et al. 2016; Heywood et al. 2022), VLA 1.4 GHz (1.5″ resolution, Schinnerer et al. 2010), and VLA 3 GHz (0.75″ resolution, Smolčić et al. 2017) imaging. The rest-frame 1.4 GHz luminosity (L_1.4) is derived from the radio dataset with the highest available signal-to-noise (S/N), and K-corrected to the 1.4 GHz rest frame by assuming a power-law spectrum with a spectral index, α, between 3 GHz and 1.28/1.4 GHz, or assuming α = 0.75 if not jointly detected. To match the spatial resolutions between the 3 GHz and 1.4 GHz images, we downgraded the 3 GHz map by convolving it with a Gaussian kernel of 1.5″ full width at half maximum (FWHM) (similarly to, e.g., Gim et al. 2019). We refrained from using a larger FWHM (i.e., for 1.28 GHz counterparts from MeerKAT), as Delhaize et al. (2017) showed that the flux increase above 1.5″ is typically negligible. The new 3 GHz fluxes at 1.5″ resolution are higher by 6%, on average, compared to those of the native resolution (0.75″) map. This effect induced a minimal flattening in the re-calculated spectral index for the joint 1.4−3 GHz detections (for a in-depth discussion about the influence of the radio spectral index on the IRRC, see Sect. 4.4.1 in Delhaize et al. 2017).

As for the IR part, we required the combined S/N (summed in quadrature over all Herschel, SCUBA2, AzTEC, and ALMA bands) be higher than five. For the same reason, for the radio part, we required the 1.4 GHz flux to be higher than 25 μJy; that is, 5× the VLBA 1.4 GHz sensitivity at the pointing center. We restricted ourselves to galaxies at z> 0.5 in order to obtain the full galaxy size within the VLBA field of view limited by bandwidth smearing (i.e., 5.5″ diameter). We disregarded all previously detected VLBA targets from Herrera-Ruiz et al. (2017, 2018) to avoid duplications.

These selection criteria naturally identify 500 massive (M_★ ≳ 10¹⁰ M_⊙) galaxies. This sample is representative of star-forming galaxies on the IRRC, as they are broadly distributed in q_IR around the (M_★, z)-dependent IRRC from Delvecchio et al. (2021) (see also Sect. 5 for details).

3. Observation and data calibration

Our VLBA observations are summarized in the following, while further details are given in D26. In each VLBA session, the full L-band bandwidth is split into eight spectral windows of 32 MHz each, further divided into 16 spectral channels (2 MHz each) at 4 s time resolutions. This setup yields a field of view limited by a time (bandwidth) smearing of 10.82″ (2.75″), which is hence larger than the typical galaxy size at z > 0.5 (Yang et al. 2025). The final images of all targets, consisting of 4096 × 4096 pixels (at 0.001″/px scale) weighing about 67 Mbyte each³.

The pointing center was located at RA (J2000) = 10:00:36.3299 and Dec (J2000) = +02:12:58.773 and covers 30′ FWHM at 1.4 GHz. It was observed for a total of 120 h, split into 24 epochs of 5 h each, due to the limited visibility of COSMOS by all ten VLBA antennas for longer tracks. Raw data were calibrated using the VLBI PIPEline (VPIPE; Radcliffe et al. 2024, and in prep.⁴) that is designed to perform automated instrumental calibration, phase referencing, and imaging for multiphase center VLBI observations. All these steps have been successfully tested in former very long baseline surveys (EVN from Radcliffe et al. 2018, 2021; VLBA from Deane et al. 2024).

In order to evaluate the significance of each VLBA source as a function of the peak S/N, we ran simulations using the Common Astronomy Software Applications (CASA, v.6.6) package. We refer the reader to D26 for more technical details on this task. Briefly, we considered the real noise map, taken from an empty sky region at the VLBA pointing centre, and injected one point-like source in the uv-plane, with known position and S/N. We iterated the simulation 100 times, keeping the input position fixed, but varying the input S/N over the 4.5 ≤ S/N ≤ 7.5 range. Only one source was injected into each realization in order to mitigate side-lobe effects and obtain reliable output fluxes. We then imaged each simulated dataset, adopting the same setup used for real data and employing BLOBCAT (Hales et al. 2012) to extract the flux and position of all blobs above a given peak S/N in each image. While the completeness at S/N > 4.5 is always observed to be higher than 90% (i.e., BLOBCAT finds the one input source in > 90% of runs), the fraction of spurious blobs is a strong function of the S/N. In particular, we defined the “spuriousness” (p_spur) as

$$\begin{array}{c}\hfill p_{\mathrm{spur}}=1-\frac{N_{\mathrm{inp}}}{N_{\mathrm{out}}},\end{array}$$(3)

where N_inp is the number of “real” (input) blobs and N_out is the number of the (output) detected blobs. We counted the output blobs within a radius of ±0.2″ from the input (real) source, similarly to what we did for real detections (see Sect. 4). We find that the spuriousness is higher than 50% at S/N < 5, while it decreases to ≈15% at S/N ∼ 5.5, and it exponentially drops below 2% at S/N ≥ 6. If counting all random blobs over the full image (4″ × 4″), the spuriousness at a fixed S/N increases proportionally to the searching area, except at S/N > 6 where it drops to below 2% as found in the central area.

4. Final sample of VLBA detections

Following Herrera-Ruiz et al. (2017) (HR2017 hereafter), we obtain the total fluxes from BLOBCAT (Hales et al. 2012) by integrating the blobs with S/N ≥ 5.5 that are found within ±0.2″; that is, ±1 VLA 3 GHz pixel from the prior VLA counterpart position (Smolčić et al. 2017), which sets our maximum VLBA-VLA separation – in HR2017 it was ±0.4″, i.e., ±1 VLA 1.4 GHz pixel (Schinnerer et al. 2010). For consistency, our chosen size of ±0.2″ is also the circular region inside which we forced the cleaning in order to maximize the S/N. Based on the above criteria, we identified four VLBA detections associated with AGN activity, with their main properties listed in Table 1. We also show them in Fig. 1, where the presence of a compact radio core is evident (e.g., purple contours). The average VLBA rms of these four detections is 6 μJy/beam. Interestingly, these sources are the faintest AGNs ever detected with VLBI techniques in SFGs. If pushing our detections threshold to a lower S/N, in particular at 5.0 < S/N < 5.5, we would identify two additional source candidates, but their high spuriousness (∼40% at their average S/N ∼ 5.2) does not allow us to claim them as secure detections.

Fig. 1. RGB images of targets VLBA ID ‘338’ (top left), ‘351’ (top right), ‘408’ (bottom left), and ‘498’ (bottom right). These were obtained by combining the F150W, F277W, and F444W JWST/NIRCam filters (Franco et al. 2025) covering a 5″ × 5″ field of view. The NIRCam pixel scale is 30 mas. The white contours are the VLA 3 GHz levels (Smolčić et al. 2017) at three, four, five, and six times the rms (2.3 μJy/beam; the VLA beam size is 0.75″ as shown by the white circle in the bottom left). In the zoomed-in images of the red panels, the purple contours are the levels at two, three, four, five, and six times the image rms detected by the VLBA, which come from a 0.2″ × 0.2″ field of view (also outlined by the red square at the centers of the galaxies). The dotted purple contour also shows the −3 rms level. The VLBA beam size for each detection is: 26 × 5 mas2 (‘338’), 26 × 5 mas2 (‘351’), 28 × 5 mas2 (‘408’), and 23 × 6 mas2 (‘498’), with a similar position angle PA ∼ 162°.

We note in Fig. 1 that the VLBA (iso-surface brightness) contours appear more complex than the elongated synthesized beam. This is because three out of four detections (all but 498) are marginally resolved (along the beam minor axis) according to BLOBCAT (Hales et al. 2012). However, given the low S/N of these detections, the limited uv coverage and the lack of a priori information on the source morphology, the imaging and deconvolution process becomes highly uncertain (e.g., Taylor et al. 1999). In such conditions, the final restored contours can be a complex blend of different spatial components and are not expected to strictly follow the shape of the clean beam; hence, some contours of faint sources would appear rounder and/or smaller than the synthesized beam. This is a known issue in deep VLBI imaging and has been observed in other VLBI surveys (e.g., Radcliffe et al. 2018). Because the three resolved detections display non-Gaussian morphologies, as recommended in Hales et al. (2012) we took total VLBA fluxes (e.g., S​_VLBA in Table 1) from the corresponding observed integrated values, but skipping the Gaussian volume correction. For the only unresolved detection (498) we instead set the total flux equal to the peak flux.

In order to make a fair comparison between our VLBA detections and non-detections, we define two control samples of VLBA non-detections (CS hereafter): (i) a M_★-matched CS (i.e., with log (M_★/M_⊙) > 10.5; 277 objects) and a combined [M_★, rms]-matched CS (i.e., a subset of the first CS, also with rms < 6.7 μJy beam⁻¹; i.e., the highest among the four detections; 42 objects). Therefore, this latter CS contains the VLBA targets that were not detected within a depth-matched area, despite them being equally as accessible as for the VLBA detections. Overall, the effective detection rate in our sample is 9% (i.e., 4/46), which is broadly consistent with the AGN duty cycle in massive star-forming galaxies (Delvecchio et al. 2022; Wang et al. 2022). We note that matching the control sample by the lowest or average rms of the detections would not have changed the effective detection rate.

5. The AGN-corrected infrared-to-1.4 GHz radio luminosity ratio

We measured q_IR exploiting the total radio VLA luminosities, where both star formation and AGN activity can, in principle, contribute. Figure 2 shows these values as a function of the stellar mass of the host galaxy. Dark stars mark the four detections, the light blue circles indicate the M_★-matched CS, while the darker blue circles represent the (M_★, rms)-matched CS. The solid black line represents the dependence of the IRRC on the stellar mass for galaxies at a ⟨z⟩∼1, as indicated by the relationship described by Delvecchio et al. (2021):

$$\begin{array}{c}\hfill q_{\mathrm{IR}}(M_{★},z)=2.646\pm 0.024\times A^{-0.023\pm 0.008}-B\times (0.148\pm 0.013),\end{array}$$(4)

Fig. 2. qIR distribution as function of host galaxy stellar mass for all 500 targets (gray points), four VLBA-detected AGNs (blue stars), and 277 non-detections part of the M★-matched sample (light blue points). The qIR of the detections account for the AGN emission detected by the VLBA (e.g., qIR, corr). The darker blue points are the 42 non-detections observed at sensitivities equal to or lower than the sensitivity of the detections, thus constituting the (M*, rms)- matched CS. The black line is the qIR versus log M★/M⊙ relation from Delvecchio et al. (2021) at ⟨z⟩∼1, and the shaded areas span the scatter of the IRRC at M★ > 1010.5 M⊙ from Delvecchio et al. (2021), of around ±0.22 dex and ±0.44 dex.

where A = (1 + z) and B = log(M_★/M_⊙).

It is important to note that the dependence on the stellar mass is stronger than that with the redshift. To measure the q_IR(M_★, z), we considered the stellar masses and redshifts presented in Sect. 3. Since the q_IR(M_★, z) distributions are similar between the M_★- and (M_★, rms) – matched CS (as expected given that the observed q_IR does not depend on the VLBA rms), in the following we consider the q_IR(M_★, z) distribution of the M_★-matched CS to increase the statistics.

If we assume that the incidence of AGNs is as high as in the deepest regions across the entire sky area mapped by the AGN-sCAN survey, we can then reason that the number of effective detections would have been seven times greater than the N_AGN = 4 detected in regions with an rms < 6.7 μJy/beam. This factor of seven was determined by noting that the subset of non-detections matched on both stellar mass and rms, (M_★, rms)-CS, is approximately seven times smaller than the subset matched solely by stellar mass (e.g., M_★-CS).

Next, we measured the difference (Δq_IR) between q_IR(M_★, z) and the observed q_IR computed using the VLA luminosities in Eq. (1). In the left panel of Fig. 3, the light purple histogram displays the distribution of Δq_IR for the (M_★, rms)-CS, while the dark purple histogram displays the distribution of Δq_IR for the detections. To reproduce the expected IRRC for pure SFGs, we built up the so-called mirror distribution (dashed orange histogram). This is obtained by mirroring the right side of the real distributions, assumed to be purely originating from SF processes, onto the left side. The errors associated with the counts of the mirror distributions (orange shaded area) were calculated with a bootstrap by extracting 1000 values from a Gaussian distribution centered on the observed number counts (N), in each bin, with a standard deviation of $\sqrt{N}$. We observe that the mirror distribution peaks around zero in both panels. This allowed us to visually assess how much the observed distribution deviates from a perfect Gaussian. In the right-hand panel of Fig. 3, we instead plot the distributions of the difference between q_IR(M_★, z) and the corrected q_IR (q_IR, corr), after subtracting the AGN emission observed with the VLBA. In this case, Δq_IR, corr equals log[(L_VLA)/(L_VLA − L_VLBA)], which is > 0 by definition.

Fig. 3. Normalized distributions of difference between the observed qIR and the one expected from Eq. (4) (Delvecchio et al. 2021) based on the stellar mass and the redshift of the host galaxy, qIR (M★, z). The light purple histogram represents the distribution of the (M★, z)-CS. The stacked dark purple histogram (i.e., dark purple histogram on top of the light purple one) represents the detections, and the hatched orange histogram is the so-called mirror distribution, interpreted as the intrinsic qIR distribution of SFGs. The bootstrapped median values of the mirror distribution were fit with a Gaussian (red curve; see Sect. 5 for further details). Left panel: qIR was computed using the VLA luminosity coming from both star formation and AGN activity. Right panel: qIR, corr is measured after subtracting the AGN radio luminosity detected by the VLBA from the total radio VLA emission.

The first thing to note is that the median values of the real distributions remain unchanged before and after the correction for the AGN contribution to the radio luminosities. Both before and after the correction, we estimate $⟨{\mathrm{\Delta }}_{\mathrm{qIR}}⟩=0.0_{-0.02}^{+0.03}$ dex. The errors on the median are calculated as the 16th and 84th percentiles of the distribution, divided by the square root of the total counts. However, it was expected that the AGN contribution would have a minimal impact on the scatter within the IRRC for these galaxies due to the low fraction of detections. Nevertheless, it is interesting to note that on an individual source basis, the AGN fluxes substantially influence q_IR, since the average difference between q_IR and q_IR, corr for the detections is about 0.14 dex (e.g., 30%), which is not negligible.

To further explore the potential role of radio variability, we examined the list of variable 1.5 GHz sources identified in the CHILES VERDES survey (Sarbadhicary et al. 2021). The survey provides 1−2 GHz continuum VLA data obtained between 2013 and 2019 in a 0.44 deg² area of the COSMOS field for 750 galaxies previously detected at VLA-1.4 GHz. Moderate variability (10−30% flux density variation) was observed in 18 sources, and low-variability (2−10% flux density variation) in 40 sources, for a total of 58/750 (6%) variable targets. Only one of those variable sources is in common with our VLBA 1.4 GHz targets (even though it is not a VLBA detection) and exhibits low VLA variability (±5%). Moreover, based on the comparison between VLBA-1.4 GHz (Herrera-Ruiz et al. 2017) and VLA-1.4 GHz (Schinnerer et al. 2010) fluxes taken at different epochs (2012 vs 2007, respectively), Herrera-Ruiz et al. (2017) found 13/468 (∼3%) VLBA detections with VLBA flux higher than the total VLA flux (at above > ±1σ uncertainty), which the authors interpreted as possible sign of variability. However, those 13 sources all have S​_vlba > 200 μJy, making them among the brightest of their sample. Hence, we expect a negligible fraction of radio-variable sources among our even fainter targets.

We also find that the real distribution exhibits a deviation from the Gaussianity (e.g., when compared to the mirror distribution) on the left side in both panels. However, this is minimal, as a K-S test could not reject the hypothesis that these two distributions come from the same parent sample (with a p value of 0.997).

Finally, we note that the mirror distributions are basically unchanged even after correcting the q_IR of the detections. Besides the peak, also their standard deviations are fully consistent (from σ = ±0.44 dex to σ = ±0.45 dex). Overall, these results highlight that VLBA AGN contamination, despite being non-negligible for individual detections, does not impact the peak position (and the overall shape) of the original q_IR distribution at all.

For completeness, in Appendix A, we correct the values of q_IR for the dependence on redshift. We also verified that in that case the distribution of the z-corrected q_IR distribution remains unchanged before and after the AGN correction.

6. Source counts

Following a similar approach to Deane et al. (2024) and Herrera-Ruiz et al. (2018), we computed the differential source counts, dN/dS_ν, as a function of source 1.4 GHz flux density, S​_ν, using the expression

$$\begin{array}{c}\hfill \frac{\mathit{dN}}{\mathit{dS}}=\frac{N}{\mathrm{\Omega }\mathrm{\Delta }S},\end{array}$$(5)

where N is the number of sources in a flux-density bin, Ω is the area over which a source with flux density equal to the mean flux density of the bin could be detected, and ΔS is the width of the flux-density bin.

Figure 4 displays the 1.4 GHz source count of the VLBA-detected AGN presented in this work (red star). For comparison, we also show the differential source counts derived by Herrera-Ruiz et al. (2018) using the VLBA COSMOS survey (as black dots), the derived source counts for the VLBA CANDELS GOODS-North Survey from Deane et al. (2024) (as red circles), and the VLA 3 GHz source counts split into SF and AGN subpopulations from Novak et al. (2018). The Novak et al. (2018) curves were obtained by statistically separating SFGs and radio AGNs based on the “radio-excess” criterion (Delvecchio et al. 2017) using data from the VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017). Our results are perfectly consistent with these previous works, as indeed the VLBA source counts at S​_1.4 GHz = 37 − 58 μJy follow the red curve in Novak et al. (2018), which is representative of the AGN population. At these fluxes, the AGN number counts are ∼5× lower than those of the total radio source population. This corroborates that the average AGN fraction of ∼20% among sources within this flux regime, while ∼80% of the sources are SFGs.

Fig. 4. Euclidean normalized radio-source number counts obtained by fitting with the Markov chain Monte Carlo algorithm, and the total radio luminosity function (LF) from Novak et al. (2018) using different evolving analytical LFs for AGN hosts (red line) and SF galaxies (blue line). The shaded areas encompass the 3σ errors from the χ2 fits performed on individual populations (for further details, see Novak et al. 2018). Gray circles are the number counts from Galluzzi et al. (2025) obtained from a compilation of radio observations with different instruments (e.g., VLA 1−3 GHz, ATCA, GMRT, LOFAR, MeerKAT, etc.), without any distinction between AGNs or SF. Red and black circles show the number counts from the VLBA CANDELS GOODS-North Survey (Deane et al. 2024) and VLBA 1.4 GHz COSMOS (Herrera-Ruiz et al. 2018) in the case of radio emission from AGN activity. The red star is the number count of the detections from the AGN-sCAN survey presented in this work. The gray arrow points toward the regime below our sensitivity limit, from which the number counts extrapolated from the red curve of Novak et al. (2018) of the 42 brightest AGN were extracted, in order to constrain upper limits on our results (see text in Sect. 7.1 for further details).

Our measurement in Fig. 4 (red starred symbol) is the faintest source-count measurement ever obtained from VLBI studies, and it aligns well with previous (low-resolution) extrapolations of AGN number counts (solid red line). In this context, Whittam et al. (2017) hinted at a possible flattening of the AGN number counts toward faint flux densities, meaning that compact, flat-spectrum radio cores are increasingly dominant in these regimes. Although we do not have enough statistics to robustly test this hypothesis, we remark that the excess of flat-spectrum radio cores could indicate the presence of younger radio jets, short duty cycles, weak jet collimation, or slower jet speeds (for an in-depth discussion, see Deane et al. 2024).

7. Discussion

7.1. Cumulative distribution of AGNs in the faintest regime

Our 5.5σ sensitivity limit enabled us to detect AGN activity with flux densities down to ⟨S​_1.4⟩ = 37−58 μJy. Of course, we cannot rule out a widespread contribution of AGNs at even fainter fluxes. Unveiling those ultra-faint AGNs is, however, very challenging. VLBA stacking would improve our sensitivity limit, but it is not feasible due to to the lack of prior knowledge concerning AGN positions. The VLA peak flux position was indeed not sufficient to serve as a prior, as the VLBA has over ×10 higher spatial accuracy. Similarly, HST or JWST positions are significantly less accurate than VLBA positions, with the additional caveat that they could trace star formation, which might not be spatially coincident with AGN emission. Therefore, we attempted a statistical assessment of the AGN contamination in VLBA non-detections.

We derived a quantitative threshold of how faint the required flux limit would have to be, by conservatively assuming that our VLBA non-detections from the (M_★, rms)-matched CS are the first 42 brightest AGNs expected just below the faintest flux of our VLBA detections; i.e., at S​_1.4 < 37 μJy. To do so, we converted the Euclidean-normalized 1.4 GHz source counts into a dimensionless quantity, re-scaled to the effective sky area of the (M_★, rms)-matched CS (i.e., 9.7 × 10⁻⁶ sr). This yields a number of expected AGNs in the same area, as a function of 1.4 GHz flux. We made sure that the relation between differential N_AGN and S​_1.4 yielded an integrated N_AGN = 4 in the range 36−58 μJy (i.e., matching our VLBA detections; red star in Fig. 4). Then, we extrapolated the same trend down and below our nominal VLBA flux-density limit. Finally, we computed the cumulative distribution of N_AGN at S​_1.4 < 37 μJy, and we found out that 42 sources would be cumulated at a flux density of 5 μJy, which is the minimum flux density reached by the ≈42 brightest AGNs. We stress again that this is an over-conservative assumption, since it is very unlikely that our non-detections are all clustered just below the flux limit. However, at these low flux densities, the brightness temperature drops to T_b ≪ 10⁵ K (see Eq. (2)). This suggests that, even if we were able to detect such faint sources with the VLBA, the origin of their emission would be ambiguous between AGNs or star formation. Consequently, we conclude that pushing VLBA observation to sub-μJy sensitivity would not be an effective strategy for revealing a residual hidden population of radio AGNs. For instance, according to the ngVLA exposure time calculator⁵, we note that Band 1 (2.3 GHz observations) would reach brightness temperatures of T_b ≥ 10⁵ K (at 5σ; i.e., ≈5μJy/beam) in about 16 h on source. Hence, fainter fluxes, as expected for the bulk of our VLBA non-detections, would be detected with a lower T_B < 10⁵ K, thus carrying an ambiguous origin despite their easy detection.

Furthermore, even in this conservative scenario in which our non-detections are the brightest AGNs below our sensitivity limit, the ratio S​_VLBA/S​_VLA is 17% at S​_VLA ≈ 70 μJy. S​_VLBA is the mean flux of the 42 brightest AGNs obtained from the cumulative distribution of the AGN number counts in the 5 μJy ≤ S​_1.4 ≤ 36 μJy range, while S​_VLA is the mean VLA 1.4 GHz flux density of the 42 non-detections. This ratio would translate into a global upper q_IR shift of < 0.1 dex. We note that, based on the M_★-dependent IRRC, q_IR decreases by 0.22 dex within the same M_★ range of the analysed sample (log M_★/M_⊙ > 10.5, Delvecchio et al. 2021). Therefore, given that radio AGN contamination (including VLBA detections and non-detections) counteracts this effect by < 0.1 dex, the decreasing q_IR with M_★ is not primarily explained by AGN contamination.

Among other possible mechanisms that have been proposed to explain the decline of q_IR at stellar masses log(M_★/M_⊙)≥10.5 in SFGs, there is an increasing energy-loss rate of high-energy CR electrons, which can be due to particularly efficient supernova-driven turbulence (Schober et al. 2023). The CR energy loss could also depend on the gas surface density and the redshift of the host galaxy, as modeled by Yoon (2024). Finally, Ponnada et al. (2025) found that in many cases “radio-excess” objects can be better understood as “IR-dimmed” objects with longer-lived radio contributions at low-z from Type-Ia supernovae and intermittent black-hole accretion in weak SFGs.

7.2. Diffuse AGN emission as a possible caveat

A potential caveat of this analysis is that the AGN flux could be underestimated because of extended emission (beyond the VLBA beam; e.g., 200 pc × 50 pc at z = 1), that is resolved out by the VLBA (e.g., Morabito et al. 2025). We first tested that the integrated fluxes remain consistent across various tapering and Briggs schemes in our images. We also note that additional AGN flux over more extended physical scales would correspond to a lower T_b; hence, T_b < 10⁵ K (see Eq. 2).

A statistical assessment of this AGN emission at intermediate scales comes from Herrera-Ruiz et al. (2017). The vast majority of their VLBA detections (> 90%) are classified as AGNs based on the radio excess observed from the VLA fluxes. Briefly, the radio excess is determined by the difference between the radio luminosity predicted by the IRRC at the redshift and stellar mass of the host galaxy (which traces star formation only) and the observed radio luminosity measured by the VLA, which includes contributions from both star formation and AGNs. Assuming that this radio-excess traces the full AGN-related flux over the entire galaxy size, the comparison between total VLBA and AGN-driven VLA flux indicates how much flux is lost on larger scales. We note (see Fig. 11 in Herrera-Ruiz et al. 2017) that at fluxes of S​_1.4 < 50 μJy, the VLBA flux is already > 75% of the total VLA AGN flux; hence, it is highly compact. By extrapolating this trend toward fainter flux regimes (as for our VLBA targets), we would expect that the VLBA flux is capturing the vast majority (if not all) of the total AGN flux. Finally, the agreement between the predicted and observed number count in this work (Fig. 4) further supports the hypothesis that the VLBA effectively captures the full AGN emission from our targets.

8. Conclusions

We present our 1.4 GHz VLBA observations carried out in the context of the AGN-sCAN survey (PI: Delvecchio), targeting a sample of 500 SFGs in COSMOS at z > 0.5 located around the IRRC (see Fig. A.1). By exploiting this sample down to a limiting T_B ≳ 10⁵ K, we were able to test whether the decrease of q_IR with galaxy stellar mass, as observed in Delvecchio et al. (2021), could be explained by residual radio contamination from radio-faint AGN buried in the total radio emission of massive star-forming galaxies. To do so, we investigated the difference in the q_IR distribution before and after correcting for a VLBA AGN contribution.

To perform a fair comparison between detections and non-detections, we built up a (M_★, rms)-matched CS of 42 pure SFGs, which did not show AGN emission as detected by the VLBA at 1.4 GHz and have a stellar mass of log(M_*/M_⊙)≥10.5.

Our findings can be summarized as follows:

• (a)

  We revealed a radio core in four targets placed in the deepest area. The effective AGN detection rate (rms-matched) in SFGs around the IRRC is ∼9%. By subtracting the VLBA AGN contribution at 1.4 GHz, we find that the q_IR distribution remains essentially unchanged. Indeed, the median and dispersion values of the distributions are entirely consistent with each other before and after applying this correction (as shown in the left and right panels of Fig. 3). Given that such an AGN correction is negligible for log(M_★/M_⊙)≥10.5, we expect that it will be even less significant at lower stellar masses.

• (b)

  As shown in Fig. 4, the AGN number counts at 37−58 μJy lie well on the extrapolation of the AGN number counts from Novak et al. (2018), and the AGN contribution amounts to 20% of the total number count of the entire radio population at these fluxes. Moreover, even in the extreme case that our VLBA non-detections were the brightest AGNs just below the detection threshold, we still predict a maximum q_IR shift of < +0.1 dex. This upper limit is significantly smaller than the 0.22 dex decline of q_IR within the M_★ range of the VLBA detections, thus excluding the possibility that radio-AGN contamination is the primary driver of the M_★-dependent q_IR.

Finally, we note that pursuing deeper VLBA observations would not be an effective strategy to unveil a residual hidden population of even fainter AGNs. Other potential physical mechanisms that could explain the q_IR decline in massive galaxies include CRe transport (Yoon 2024) and a longer lived radio contribution from Type Ia supernovae (Ponnada et al. 2025). Nevertheless, our results highlight that the SFR-radio relation is robust against contamination from AGN activity, once the so-called radio-excess AGNs are removed.

Acknowledgments

The authors thank the referee for a constructive report that helped us clarify the results and implications of this work. G. P., I. D., M. G., F. U., and I. P. acknowledge funding by the European Union – NextGenerationEU, RRF M4C2 1.1, Project 2022JZJBHM: “AGN-sCAN: zooming-in on the AGN-galaxy connection since the cosmic noon” – CUP C53D23001120006. C. S. acknowledges financial support from INAF – Ricerca Fondamentale 2024 (Ob.Fu. 1.05.24.07.04), from the Italian Ministry of Foreign Affairs and International Cooperation (grant PGRZA23GR03), and by the Italian Ministry of University and Research (grant FIS-2023-01611, CUP C53C25000300001). IP also acknowledges support from INAF under the Large Grant 2022 funding scheme (project “MeerKAT and LOFAR Team up: a Unique Radio Window on Galaxy/AGN co-Evolution”). W.W. acknowledges the grant support through JWST programs. Support for programs JWST-GO-03045 and JWST-GO-03950 were provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. This research has also been supported by the European Regional Development Fund under grant agreement PK.1.1.10.0007 (DATACROSS). This work is based, in part, on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST.

References

Appendix A: Evolution of q_IR with redshift

In Fig. A.1 we show that q_IR evolves with the redshift following a power law (1 + z)^A (e.g. black line). This is because our targets were drawn from a flux-limited sample of joint IR and radio detections. From a non-linear squares (NLS) method, we find A = −0.13 ± 0.03, in agreement with e.g. Delhaize et al. (2017). The same figure shows that our 4 VLBA detections (blue stars) are placed at different redshifts. The statistic is too low to explore the observed q_IR distribution in multiple redshift bins, and for this reason, we decided to re-scale the q_IR values to the same redshift, choosing a reference value z₀. We compute the z-corrected q_IR using the expression:

$$\begin{array}{c}\hfill q_{\mathrm{IR},{\mathrm{z}}_0}=q_{\mathrm{IR}}(z){\left(\frac{1+z_0}{1+z}\right)}^A\end{array}$$(A.1)

where z₀ is the reference redshift, that we set to an arbitrary value of z₀ = 0.5. Finally, we show the 42 SFG (green dots in Fig. A.1) of the (M_★, rms) - CS, by considering those VLBA observations in regions deep as or deeper than 6.7 μJy beam⁻¹ (i.e., the poorest sensitivity reached among the detections). We find that the distributions of the mass-matched CS and the mass- and rms-matched CS spread randomly around the IRRC, indicating that systematic effects are absent.

Fig. A.1. qIR as a function of redshift for all the 500 targets (grey points), the 4 VLBA-detected AGN (dark blue stars) and the 277 non-detections that constitute the M★-matched sample (red points). The green dots are the 42 non-detections which make the (M★, rms)-matched control sample. The black line is the NLS fit of all the targets described by the power law (1 + z)−0.13.

Fig. A.2. Normalized qIR distributions re-scaled to z0 = 0.5 (qIR, z0) for the control sample of VLBA non-detections (yellow histogram), while the 4 detections (light blue histogram) are shown both before (bottom panel) and after the VLBA AGN correction (top panel). The red Gaussian is the fit of the yellow+blue (or yellow) histograms, while the black curve is the Gaussian fit to the bootstrapped mirror distribution.

Then, we compute the corrected q_IR, z0 (q_{IR, z0, corr}) of the 4 detections by subtracting the VLBA fluxes arising from the AGN, which is log(L_{1.4GHz, VLA} − L_{1.4GHz, VLBA}). Figure A.2 shows the distribution of the q_{IR, z0, corr} of the 4 detections (blue histogram) and the q_IR, z0 of the control sample (yellow histogram). The bin size is 0.2 dex, that is about twice larger than the average q_IR uncertainty.

The red curve is the Gaussian fit of the real distribution. The gray curve is obtained by extracting 1000 values from the so-called “mirror distribution” (hatched histogram) using a bootstrap method, which would represent the distribution observed if all galaxies were star-forming. The errors used in the bootstrap are the 16^th/84^th percentiles of the 1000 values. The mirror distribution is obtained by unfolding the right-hand side of the distribution (above the peak, which is assumed to originate from pure SF) on the left side. A K-S test was not able to exclude that the mirror and the real distribution are drawn from different parent samples (with a p-value of 0.99). Accordingly, the median value of the real distribution, which is μ_R = 2.67 ${}_{-0.02}^{0.02}$, and the median of the bootstrapped-mirror distribution (${\mu }_{\overline{M}}$ = 2.69 ${}_{-0.04}^{0.04}$), are consistent between the errors on the medians, computed as the 16^th/84^th percentiles of the respective distributions divided by the root squared number of the total counts. Thus, as expected, the real distribution matches the mirrored one after subtracting the AGN contribution from the q_IR, z0.

To explore the effects of the correction from the AGN contribution, we study q_IR, z0 in two different scenarios: (1) without the VLBA correction for both detections and non-detections; and (2) only for the control sample of non-detections. Interestingly, the real and mirrored distributions are still statistically indistinguishable in the (1) case. The bottom panel of Fig. A.2, indeed, shows that μ_R = 2.65 ${}_{-0.01}^{0.02}$ and ${\mu }_{\overline{M}}$ = 2.68 ${}_{-0.02}^{+0.02}$, which are (again) consistent within the error bars between each other and with the values in Fig. A.2. However, the difference between μ_R and μ_M is ∼ 0.03 dex and larger than when considering the AGN-corrected q_IR, z0. We also observe a slight shift on the left side between the red and black curves, reflecting the fact that all the detections have q_IR, z0 < μ_R. The shift, indeed, disappears in the central panel of Fig. A.2 when we consider SFG-only from the control sample. Finally, the values of μ_R are identical among the (1) and (2) cases. In conclusion, we can robustly claim that VLBA AGN contamination, despite being non-negligible for individual detections, does not impact at all the peak position (and the overall shape) of the original q_IR distribution.

All Tables

All Figures

Fig. 1. RGB images of targets VLBA ID ‘338’ (top left), ‘351’ (top right), ‘408’ (bottom left), and ‘498’ (bottom right). These were obtained by combining the F150W, F277W, and F444W JWST/NIRCam filters (Franco et al. 2025) covering a 5″ × 5″ field of view. The NIRCam pixel scale is 30 mas. The white contours are the VLA 3 GHz levels (Smolčić et al. 2017) at three, four, five, and six times the rms (2.3 μJy/beam; the VLA beam size is 0.75″ as shown by the white circle in the bottom left). In the zoomed-in images of the red panels, the purple contours are the levels at two, three, four, five, and six times the image rms detected by the VLBA, which come from a 0.2″ × 0.2″ field of view (also outlined by the red square at the centers of the galaxies). The dotted purple contour also shows the −3 rms level. The VLBA beam size for each detection is: 26 × 5 mas2 (‘338’), 26 × 5 mas2 (‘351’), 28 × 5 mas2 (‘408’), and 23 × 6 mas2 (‘498’), with a similar position angle PA ∼ 162°.

In the text

Fig. 2. qIR distribution as function of host galaxy stellar mass for all 500 targets (gray points), four VLBA-detected AGNs (blue stars), and 277 non-detections part of the M★-matched sample (light blue points). The qIR of the detections account for the AGN emission detected by the VLBA (e.g., qIR, corr). The darker blue points are the 42 non-detections observed at sensitivities equal to or lower than the sensitivity of the detections, thus constituting the (M*, rms)- matched CS. The black line is the qIR versus log M★/M⊙ relation from Delvecchio et al. (2021) at ⟨z⟩∼1, and the shaded areas span the scatter of the IRRC at M★ > 1010.5 M⊙ from Delvecchio et al. (2021), of around ±0.22 dex and ±0.44 dex.

In the text

Fig. 3. Normalized distributions of difference between the observed qIR and the one expected from Eq. (4) (Delvecchio et al. 2021) based on the stellar mass and the redshift of the host galaxy, qIR (M★, z). The light purple histogram represents the distribution of the (M★, z)-CS. The stacked dark purple histogram (i.e., dark purple histogram on top of the light purple one) represents the detections, and the hatched orange histogram is the so-called mirror distribution, interpreted as the intrinsic qIR distribution of SFGs. The bootstrapped median values of the mirror distribution were fit with a Gaussian (red curve; see Sect. 5 for further details). Left panel: qIR was computed using the VLA luminosity coming from both star formation and AGN activity. Right panel: qIR, corr is measured after subtracting the AGN radio luminosity detected by the VLBA from the total radio VLA emission.

In the text

Fig. 4. Euclidean normalized radio-source number counts obtained by fitting with the Markov chain Monte Carlo algorithm, and the total radio luminosity function (LF) from Novak et al. (2018) using different evolving analytical LFs for AGN hosts (red line) and SF galaxies (blue line). The shaded areas encompass the 3σ errors from the χ2 fits performed on individual populations (for further details, see Novak et al. 2018). Gray circles are the number counts from Galluzzi et al. (2025) obtained from a compilation of radio observations with different instruments (e.g., VLA 1−3 GHz, ATCA, GMRT, LOFAR, MeerKAT, etc.), without any distinction between AGNs or SF. Red and black circles show the number counts from the VLBA CANDELS GOODS-North Survey (Deane et al. 2024) and VLBA 1.4 GHz COSMOS (Herrera-Ruiz et al. 2018) in the case of radio emission from AGN activity. The red star is the number count of the detections from the AGN-sCAN survey presented in this work. The gray arrow points toward the regime below our sensitivity limit, from which the number counts extrapolated from the red curve of Novak et al. (2018) of the 42 brightest AGN were extracted, in order to constrain upper limits on our results (see text in Sect. 7.1 for further details).

In the text

	Fig. A.1. qIR as a function of redshift for all the 500 targets (grey points), the 4 VLBA-detected AGN (dark blue stars) and the 277 non-detections that constitute the M★-matched sample (red points). The green dots are the 42 non-detections which make the (M★, rms)-matched control sample. The black line is the NLS fit of all the targets described by the power law (1 + z)−0.13.
In the text

Fig. A.2. Normalized qIR distributions re-scaled to z0 = 0.5 (qIR, z0) for the control sample of VLBA non-detections (yellow histogram), while the 4 detections (light blue histogram) are shown both before (bottom panel) and after the VLBA AGN correction (top panel). The red Gaussian is the fit of the yellow+blue (or yellow) histograms, while the black curve is the Gaussian fit to the bootstrapped mirror distribution.

In the text