© I. Delvecchio et al. 2021

1. Introduction

For nearly fifty years, astronomers have studied the observed correlation between total infrared (IR; rest-frame 8−1000 μm, i.e. L_IR) and radio (e.g., rest-frame 1.4 GHz, i.e. L_1.4 GHz) spectral luminosity arising from star formation, typically referred to as the ‘infrared-radio correlation’ (IRRC, e.g., Helou et al. 1985; de Jong et al. 1985). This tight (1σ ∼ 0.16 dex, e.g., Molnár et al. 2021) correlation is often parametrised via the IR-to-radio luminosity ratio q_IR, defined as (e.g., Helou et al. 1985; Yun et al. 2001):

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

where 3.75 × 10¹² Hz represents the central frequency over the far-infrared (FIR, rest-frame 42−122 μm) domain, usually scaled to IR in the recent literature. In the Local Universe, the IRRC (or its parametrisation q_IR) appears to hold over at least three orders of magnitude in both L_IR and L_1.4 GHz (e.g., Helou et al. 1985; Condon 1992; Yun et al. 2001). Broadly speaking, this is because the infrared emission comes from dust heated by fairly massive (≳5 M_⊙) OB stars, while radio emission arises from relativistic cosmic ray electrons (CRe) accelerated by shock waves produced when massive stars (≳8 M_⊙) explode as supernovae. Nevertheless, CRe are also subject to different cooling processes as they propagate throughout the galaxy, which are mainly caused by inverse Compton, bremsstrahlung, and ionisation losses (e.g., Murphy 2009).

Surprisingly enough, despite all such processes at play, infrared and radio emission are observed to be correlated, both in local star-forming late-type galaxies and even in merging galaxies (e.g., Condon et al. 1993, 2002; Murphy 2013). This has been a strong motivator for the use of radio-continuum emission as a dust-unbiased star formation rate (SFR) tracer also in the faint radio sky (e.g., Condon 1992; Bell 2003; Murphy et al. 2011, 2012). Moreover, measuring the offset from the IRRC has been widely used to indirectly identify radio-excess active galactic nuclei (AGN; e.g., Donley et al. 2005; Del Moro et al. 2013; Bonzini et al. 2015; Delvecchio et al. 2017).

These applications, however, significantly rely on a proper understanding of whether and how the IRRC evolves over cosmic time and across different types of galaxies. Despite its extensive application in extragalactic astronomy, the detailed physical origins of the IRRC and the nature of its cosmic evolution have long been debated (e.g., Harwit & Pacini 1975; Rickard & Harvey 1984; de Jong et al. 1985; Helou et al. 1985; Hummel et al. 1988; Condon 1992; Garrett 2002; Appleton et al. 2004; Murphy et al. 2008; Jarvis et al. 2010; Sargent et al. 2010; Ivison et al. 2010a,b; Bourne et al. 2011; Smith et al. 2014; Magnelli et al. 2015; Calistro Rivera et al. 2017; Delhaize et al. 2017; Gürkan et al. 2018; Molnár et al. 2018; Algera et al. 2020a).

For example, some studies of local star-forming galaxies (SFGs), ranging from dwarf (e.g., Wu et al. 2008) to ultra-luminous infrared galaxies (ULIRGs; L_IR > 10¹² L_⊙; e.g., Yun et al. 2001) have concluded that the IRRC remains linear across a wide range of L_IR. Conversely, other studies have argued that at low luminosities the IRRC may break down, which is consistent with a non-linear trend of the form $L_{\mathrm{IR}}\propto L_{1.4\mathrm{GHz}}^{0.75-0.90}$ (e.g., Bell 2003; Hodge et al. 2008; Davies et al. 2017; Gürkan et al. 2018), which might be partly induced by dust heating from old stellar populations (Bell 2003).

Several models have attempted to explain this non linearity. On the one hand, calorimetric models assume that galaxies are optically thick in the ultraviolet (UV), so that UV emission is fully re-emitted in the IR; likewise, CRe radiate away their total energy through synchrotron emission before escaping the galaxy (e.g., Voelk 1989). These conditions might hold in the most massive (stellar mass M_⋆ ≳ 10¹⁰ M_⊙) SFGs because of their increasing compactness (i.e. the size-mass relation $R_{\mathrm{e}}\propto M_{\star }^{0.22}$, van der Wel et al. 2014), which might enhance their ability to retain the gas ejected by stars. However, this is likely to break down towards lower M_⋆ galaxies due to smaller sizes and lower obscuration (e.g., Bourne et al. 2012). On the other hand, non-calorimetric models or the optically thin scenario (Helou & Bicay 1993; Niklas & Beck 1997; Bell 2003; Lacki et al. 2010) support the argument that several physical mechanisms cancel each other out, creating a sort of ‘conspiracy’ that keeps the IRRC unexpectedly tight and linear. Indeed, both IR and radio luminosities are expected to underestimate the total SFR in low M_⋆ and low-SFR surface density galaxies (Bell 2003), inducing a departure of the IRRC from linearity. This was not, however, observed in our study. Radio synchrotron models postulate that such small galaxies are not able to prevent CRe from escaping, causing a global deficit of radio emission at a fixed SFR. Similarly, the IR domain becomes less sensitive to SFR in low-M_⋆ galaxies (e.g., Madau & Dickinson 2014), generating an IR deficit of a similar amount, which might counterbalance the radio emission and keep the IRRC linear. Understanding the discrepancy between model predictions and observations is crucial, since the linearity (or not) of the IRRC has direct implications for using radio emission as a SFR tracer.

From an observational perspective, it is widely recognised that there is a tight relation linking SFR and M_⋆ in nearly all SFGs, namely, namely the ‘main sequence’ of star formation (MS, scatter ∼0.2−0.3 dex). This relation holds from z ∼ 5 down to the Local Universe (e.g., Brinchmann et al. 2004; Noeske et al. 2007; Elbaz et al. 2011; Whitaker et al. 2012; Speagle et al. 2014; Schreiber et al. 2015; Lee et al. 2015), showing a flattening at high M_⋆ and an evolving normalisation with redshift. Because the SFR is directly linked to L_IR, especially in massive SFGs (Kennicutt 1998), the existence of the MS gives an additional argument that studying q_IR as a function of M_⋆ could be of the utmost importance for our understanding of what drives the IRRC in galaxies.

Recent studies have corroborated the idea that the IRRC slightly, but significantly, declines with redshift (Ivison et al. 2010b; Magnelli et al. 2015; Calistro Rivera et al. 2017; Delhaize et al. 2017) in the form of q_IR ∝ (1 + z)^{[ − 0.2 : −0.1]}, although the physical explanation for such an evolution is still uncertain. Somewhat different conclusions were reached by other works (e.g., Garrett 2002; Appleton et al. 2004; Ibar et al. 2008; Jarvis et al. 2010; Sargent et al. 2010; Bourne et al. 2011) that ascribe this apparent evolution to selection effects; for instance, comparing flux-limited samples, each with a different selection function.

In this regard, we note that any selection method is sensitive to brighter, that is, more massive galaxies towards higher redshifts. By binning in redshift, only a restricted range in galaxy M_⋆ becomes detectable at each redshift for any flux limited sample, thus inducing a bias as a function of z. Therefore, it is time to simultaneously examine the evolution of the IRRC as a function of M_⋆ and redshift. We emphasise that our approach is fully empirical. However, a possible dependence of the IRRC on M_⋆ is expected from some synchrotron emission models (e.g., Lacki & Thompson 2010; Schober et al. 2017) and this might reflect some combination of the underlying physics originating the IRRC (see Sect. 5).

The main goal of the present paper is to calibrate the IRRC, for the first time, as a function of both M_⋆ and redshift over a broad range. To this end, we start from an M_⋆-selected sample of > 400 000 galaxies at 0.1 < z < 4.5 collected from deep UltraVISTA images in the Cosmic Evolution Survey (Scoville et al. 2007; centered at RA = +150.11916667; Dec = +2.20583333 (J2000)). Then we leverage the new de-blended far-IR/sub-mm data (Jin et al. 2018) recently compiled in COSMOS, which allow us to circumvent blending issues due to poor angular resolution and measure L_IR for typical MS galaxies out to z ∼ 4. In addition, we exploit the deepest radio-continuum data taken from the VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017a). Individual detections are combined with stacked flux densities of non-detections, at both IR and radio frequencies to assess the average q_IR as a function of M_⋆ and redshift.

The layout of this paper is as follows. A description of the sample selection and multi-wavelength ancillary data is given in Sect. 2. We describe the stacking analysis in Sect. 3, including measurements of L_IR (Sect. 3.1) and L_1.4 GHz (Sect. 3.4). The average q_IR as a function of M_⋆ and redshift is presented in Sect. 4, where we perform a careful subtraction of radio AGN at different M_⋆ via a recursive approach. Our main results are discussed and interpreted in Sect. 5 within the framework of previous observational studies and theoretical models. The main conclusions are summarised in Sect. 6. In addition, we test our total 3 GHz flux densities in Appendix A. A detailed comparison between radio stacking results is presented in Appendix B. In Appendix C, we discuss how the final IRRC is sensitive to our AGN subtraction method. Finally, in Appendix D we quantify how different assumptions from the literature would change our main results. Throughout this paper, magnitudes are given in the AB system (Oke 1974). We assume a Chabrier (2003) initial mass function (IMF) and a ΛCDM cosmology with Ω_m = 0.30, Ω_Λ = 0.70, and H₀ = 70 km s⁻¹ Mpc⁻¹ (Spergel et al. 2003).

2. Multi-wavelength data and sample selection

In this section, we describe the creation of a K_s prior catalogue, which we used to select our parent sample in the COSMOS field. The COSMOS field (2 deg²) boasts an exquisite photometric data set, spanning from the X-rays to the radio domain¹. The most recent collection of multiwavelength photometry comes from the COSMOS2015 catalogue (Laigle et al. 2016), which contains 1 182 108 sources extracted from a stacked YJHK_s image (blue dots in Fig. 1). In particular, this catalogue joins optical photometry from Subaru Hyper-Suprime Cam (2 deg²; Capak et al. 2007) and from the Canada-France-Hawaii Telescope Legacy Survey (CFHT-LS, central 1 deg²; McCracken et al. 2001); the near-infrared (NIR) bands Y, J, H, and K_s from UltraVISTA DR2 (down to K_s< 24.5 in the central 1.5 deg², of which 0.6 deg² are covered by ultra-deep stripes with limiting K_s < 25.2; McCracken et al. 2012), and from CFHT H and K_s observations obtained with the WIRCam (K_s < 23.9 outside the UltraVISTA area; McCracken et al. 2001). Over the full 2 deg² area, mid-infrared (MIR) photometry was obtained from the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH; Steinhardt et al. 2014; Capak et al., in prep.) using 3.6−8 μm data from the Infrared Array Camera (IRAC). We refer the reader to Laigle et al. (2016) for more details.

Fig. 1. Distribution of the full COSMOS2015 (Laigle et al. 2016) source list over the COSMOS area (blue dots). The subset of 413 678 NUVrJ–based star-forming galaxies analysed in this work (red dots) includes sources from Laigle et al. (2016) and Muzzin et al. (2013) within the UltraVISTA area, with the exception of masked regions due to saturated or contaminated photometry. See Sect. 2 for details.

In order to obtain a homogeneous galaxy selection function, we limited our study to the inner UltraVISTA DR2 area, while also excluding stars and masked regions from the COSMOS2015 catalogue with less accurate photometry, which reduces the initial sample to 45% of its size (524 061 sources). Following Jin et al. (2018), we partly fill up these blank regions by adding 22 838 unmasked K_s–selected sources from the UltraVISTA catalogue of Muzzin et al. (2013) (3σ limit of K_s < 24.35 with 2″ aperture). This ensures a more complete coverage within the UltraVISTA area, with fluctuations in the prior source density of only 2.5%. This builds on our K_s prior sample of 546 899 galaxies. Given the similar selection, we confirm that excluding the slightly shallower ∼4% subsample from Muzzin et al. (2013) leaves our results unchanged. Thus, we maintained them throughout this work.

Photometric redshifts and M_⋆ estimates were retrieved from the corresponding catalogues by fitting the optical-MIR photometry using the stellar population synthesis models of Bruzual & Charlot (2003). Both redshift and M_⋆ values represent the median of the corresponding likelihood distribution. Laigle et al. (2016) report an average photometric redshift accuracy of ⟨|Δz/(1+z)|⟩ = 0.007 at z < 3, and 0.021 at 3 < z < 6. A similar accuracy of 0.013 is reached in the catalogue from Muzzin et al. (2013) at z < 4. We further inspected a subset of 5400 sources showing a skewed redshift probability distribution function (with ≳5% chance to be offset from the median by > 0.5 × [1 + z_p]). However, we verified that removing such potential redshift interlopers does not have any impact on our results. As in Jin et al. (2018), publicly available spectroscopic redshifts were collected from the new COSMOS master spectroscopic catalogue (courtesy of M. Salvato, within the COSMOS team), and were prioritised over photometric measurements if deemed reliable (z_s quality flag > 3 ∧ |z_s − z_p|< 0.1 × (1 + z_p)). Infrared/sub-mm flux densities were de-blended and re-extracted via the prior-based fitting algorithm presented in Jin et al. (2018), which we briefly describe in Sect. 2.2.

2.1. Selecting star-forming galaxies via (NUV–r)/(r–J) colours

We aim to study the infrared-radio correlation within an M_⋆-selected sample of star-forming galaxies. To this end, we make use of the rest-frame, dust-corrected (NUV − r) and (r − J) colours available in the parent catalogues (hereafter NUVrJ). As opposed to the widely used UVJ criterion, the (NUV − r) colour is more sensitive to recent star formation (10⁶ − 10⁸ yr scales, Salim et al. 2005; Arnouts et al. 2007; Davidzon et al. 2017). Therefore, this criterion enables us to better distinguish among weakly star-forming galaxies (with specific-SFR, sSFR = SFR/M_⋆ ∼ 10⁻¹⁰ yr⁻¹) and fully passive systems (sSFR < 10⁻¹¹ yr⁻¹).

We further selected galaxies with redshift 0.1 < z < 4.5 and 10⁸ < M_⋆/M_⊙ < 10¹². This leaves us with a final sample of 413 678 star-forming galaxies (red dots in Fig. 1), out of which 22 38 (5.4%) are spectroscopically confirmed. The fraction of catastrophic failures (|z_s − z_p|> 0.15 × (1 + z_s)) is only 3.4%.

Such a sizable sample enables us to bin galaxies as a function of both M_⋆ and redshift, while maintaining good statistical power. Figure 2 shows our sample in the M_⋆–redshift diagram, highlighting the chosen grid. We note that the M_⋆ uncertainties taken from the parent catalogues incorporate the covariant errors on stellar population ages and dust reddening. These average M_⋆ uncertainties are 0.2 dex at 10⁸ < M_⋆/M_⊙ < 10⁹ and 0.1 dex above, which is far smaller than the corresponding M_⋆ bin width, thus not impacting our results. The 90% M_⋆ completeness limit (orange solid line, Laigle et al. 2016) indicates that our sample of SFGs is mostly complete down to 10¹⁰ M_⊙ out to z ∼ 4. Although we acknowledge the increasing incompleteness towards less massive galaxies in the early Universe, we believe that including them brings a valuable addition for constraining the infrared and radio properties of galaxies down to a poorly explored regime of M_⋆. This will become particularly relevant for the next generation of telescopes, such as JWST and SKA, which will routinely observe such faint sources. In addition, as we will discuss in Sect. 3.3, a very good agreement is observed between our stacked L_IR and those extrapolated from the MS relation (Schreiber et al. 2015) also at M_⋆ < 10^9.5 M_⊙, suggesting that even in this incomplete, low-M_⋆ regime our galaxies are still representative of an M_⋆-selected sample. We emphasise that the overall conclusions of this work remain unchanged when we limit ourselves to z < 3 and M_⋆ > 10^9.5 M_⊙, where our sample is largely complete. Moreover, in light of our main result, namely, that q_IR decreases with M_⋆, we anticipate that including galaxies within an incomplete M_⋆ regime would, at most, amplify the final M_⋆ dependence, thereby reinforcing our findings.

Fig. 2. Distribution of NUVrJ star-forming galaxies as a function of M⋆ and redshift. The colour scale in the central panel indicates the underlying sample size, increasing from blue to red. The grey-dashed grid encloses the 42 M⋆ − z bins into which we split our sample (413 678 objects). Galaxies within that grid are projected on the upper-left and bottom right histograms with redshift and M⋆, respectively (grey lines). The blue and red histograms represent the subsample with a signal-to-noise ratio (S/N) that is greater than 3 at 3 GHz and across all IR bands, respectively (see Sects. 2.2 and 2.3). The orange solid line marks the 90% M⋆ completeness limit of Laigle et al. (2016) for comparison.

2.2. Infrared and sub-mm de-blended data

We complemented the existing COSMOS optical-to-IRAC photometry with cutting-edge de-blended photometry from Jin et al. (2018), which is based on the de-blending algorithm developed in Liu et al. (2018) for the GOODS-North field.

The dataset includes Spitzer-MIPS 24 μm data (PI: D. Sanders; Le Floc’h et al. 2009); Herschel imaging from the PACS (100−160 μm, Poglitsch et al. 2010) and the SPIRE (250, 350, and 500 μm, Griffin et al. 2010) instruments, as part of the PEP (Lutz et al. 2011) and HerMES (Oliver et al. 2012) programmes, respectively. In addition, JCMT/SCUBA2 (850 μm) images are taken from the S2CLS programme (Cowie et al. 2017; Geach et al. 2017), the ASTE/AzTEC (1.1 mm) data are nested maps from Aretxaga et al. (2011) over a sub-area of 0.72 deg². Finally, Jin et al. (2018) also included MAMBO data (Bertoldi et al. 2007) at 1.2 mm over an area of 0.11 deg².

Briefly, Jin et al. (2018) used K_s-selected sources from the UltraVISTA survey (Sect. 2) as priors to perform a point square function (PSF) fitting of MIPS 24 μm, VLA-3 GHz (Smolčić et al. 2017a), and VLA-1.4 GHz (Schinnerer et al. 2010) images down to the 3σ level in each band. Within our final sample, this procedure identifies 67 114 MIPS 24 μm+VLA priors. Nevertheless, adopting a similar approach for extracting FIR/sub-mm flux densities of all M_⋆-selected galaxies, that is, using the full list of K_s priors, would identify up to 50 sources per beam at the resolution of the FIR/sub-mm wavelengths, leading to heavy confusion. Therefore, following Jin et al. (2018), only an M_⋆-complete subset of K_s priors was added, which ultimately prioritises IR brighter sources. This leads to a total of 136 584 K_s+MIPS 24 μm+VLA priors, that were used to de-blend and extract the Herschel, SCUBA2, and AzTEC flux densities (see Table 1). Within our final sample of 413 678 star-forming galaxies, 20 777 (5%) have a combined S/N > 3 over all FIR/sub-mm bands (10 285 at S/N > 5). These are displayed as red histograms in Fig. 2. The rest of the K_s sources are assumed to have negligible FIR/sub-mm flux densities, consistent with the background level in those bands. This is confirmed by the Gaussian-like behaviour of the noise (centered at zero) in the residual maps, after subtracting all S/N > 3 sources in each band (Jin et al. 2018). Throughout the rest of this paper, we interpret individual S/N > 3 sources as detections, while S/N < 3 sources will be stacked, as described in Sect. 3.1.

2.3. Radio data in the COSMOS field

For our analysis, we exploited data from the VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017a), one of the largest and deepest radio survey ever conducted over a medium sky area like COSMOS. With 384 h of observations, the final mosaic reaches a median root mean square (rms) of 2.3 μJy beam⁻¹ over 2.6 deg² at an angular resolution of 0.75″, the highest among radio surveys in COSMOS. A total of 10 830 sources were blindly extracted down to S/N > 5.

In addition, the COSMOS field boasts the current deepest radio-continuum data at 1.28 GHz from the MeerKAT International GHz Tiered Extragalactic Exploration (MIGHTEE, Jarvis et al. 2016) survey. With only 17 h of on-source observations over the central 1 deg² of COSMOS, the early science release of MIGHTEE reaches a thermal noise of 2.2 μJy beam⁻¹, although the effective noise is limited by confusion. This provides an excellent sensitivity to large-scale radio emission. Nevertheless, given the relatively large beam size (8.4″ × 6.8″ FWHM), MIGHTEE flux densities were de-blended using the same list of K_s + MIPS 24 μm + VLA priors applied on FIR/sub-mm images (Sect. 2.2). Similarly, VLA radio flux densities were re-extracted based on the same PSF fitting technique down to S/N > 3. However, for the sake of consistency with publicly available VLA catalogues, we take the 1.4 and 3 GHz radio flux densities of S/N > 5 sources from Schinnerer et al. (2010) and Smolčić et al. (2017a), respectively. As a sanity check, we show in Appendix A that our procedure leads to fully consistent total 3 GHz flux densities.

Given its unparalleled depth and resolution over the full COSMOS area, we primarily use the VLA 3 GHz dataset for our radio analysis, which counts 13 808 sources with S/N > 3 out of 413 678 M_⋆-selected star-forming galaxies (∼3%). These radio detections are shown as blue histograms in Fig. 2. Fainter sources will be accounted for via stacking, as described in Sect. 3.4. Nevertheless, repeating the same stacking analysis with ancillary radio datasets at 1.3 GHz (MIGHTEE) and 1.4 GHz (VLA) is essential to validating our procedure in the face of any potential variations of radio spectral index or different angular resolutions. We refer the reader to Appendix B for an extensive comparison between stacking at 3 GHz and with ancillary radio datasets, whereas throughout this paper, we use radio data only at 3 GHz.

3. Stacking analysis

The aim of this paper is to investigate how the IRRC evolves with M_⋆ and redshift simultaneously. Contrary to studies in which galaxies were individually detected at IR or radio wavelengths, leading to complex selection functions and biased samples (see discussion in Sargent et al. 2010), we start from a well-defined M_⋆-selected sample. As a consequence, our analysis makes use of IR (Sect. 3.1) and radio (Sect. 3.4) stacking. This includes a careful treatment of some common caveats concerning IR galaxy samples, such as ‘clustering bias’ (Sect. 3.2) and spectral energy distribution (SED) fitting including AGN templates (Sect. 3.3). As for stacking radio data, special care is devoted to statistically removing radio AGN from our sample (Sect. 4.2).

In addition, our notably large star-forming galaxy sample allows us to bin as a function of both M_⋆ and redshift, as shown in Fig. 3. For each bin, we also report the total number of M_⋆-selected SFGs (black), as well as the corresponding fractions having combined S/N_IR > 3 (red) and S/N_3 GHz > 3 (blue). As we can see, both fractions are a strong function of both M_⋆ and redshift. Therefore, binning along both parameters enables us to account for the fact that galaxies of distinct M_⋆ are detectable at IR and radio wavelengths over different redshift ranges. These aspects will be extensively discussed when comparing our results with previous literature (Appendix D).

Fig. 3. Number of NUVrJ–based star-forming galaxies analysed in this work, as a function of M⋆ and redshift (black). For convenience, in each bin we report the fraction of sources with combined S/NIR > 3 over all FIR/sub-mm bands (red) and with S/N3 GHz > 3 (blue).

3.1. Infrared and sub-mm stacking

In this section, we estimate the average flux densities across eight infrared and sub-mm bands, namely MIPS 24 μm, PACS 100−160 μm, SPIRE 250−350−500 μm, SCUBA 850 μm, and AzTEC 1100 μm. Similarly to other studies, we perform a median stacking on the residual maps from Jin et al. (2018), that is, after subtracting all detected sources with S/N > 3 in each band (see also Magnelli et al. 2009). Individual S/N > 3 detections will be added to stacked flux densities a posteriori through a weighted average (Eq. (3)). Median stacking strongly mitigates contamination from bright neighbors and catastrophic outliers, and thus reduces the confusion noise for the faint sources. We stress that our procedure yields very consistent results with either median or mean stacking of detections and non-detections combined (e.g., Magnelli et al. 2015; Schreiber et al. 2015), as shown in Sect. 3.3.

To produce stacked and rms images in each band, we used the publicly available IAS stacking library² (Bavouzet et al. 2008; Béthermin et al. 2010). For each band, M_⋆ bin, and redshift bin, we stack N × N pixel cutouts from the residual images, each centred on the NIR position of the M_⋆-selected priors (Sect. 2). We chose the cutout size as eight times the full-width at half maximum (FWHM) of the PSF, while for Spitzer-MIPS, we chose 13 × FWHM, since a substantial fraction of the 24 μm flux is located in the first Airy ring. Since the AzTEC map covers only a central sub-area of 0.72 deg², at 1.1 mm we only stacked within that region. We emphasise that the M_⋆, z, and SFR distribution of the SFG population within the AzTEC region is fully consistent with that derived in the rest of the COSMOS field, effectively not biasing the resulting stacked flux densities. To measure total flux densities, we followed different techniques depending on the input map. For MIPS and PACS images, we used a PSF fitting technique (e.g., Magnelli et al. 2014). A correction of 12% is further applied to account for flux losses from the high-pass-filtering processing of PACS images (e.g., Popesso et al. 2012; Magnelli et al. 2013). For SPIRE images, the photometric uncertainties are not dominated by instrumental noise but by the confusion noise caused by neighboring sources (Dole et al. 2003; Nguyen et al. 2010). Since SPIRE images, as well as those from SCUBA and AzTEC, are already brightness maps (in units of Jy beam⁻¹ or equivalent), we read the total flux from the peak brightness pixel. The total flux density was taken as the median of the input cube at the central pixel.

The uncertainties on the stacked flux densities were measured using a bootstrap technique (e.g., Béthermin et al. 2015). Within each M_⋆ − z bin, we ran our stacking procedure 100 times, in all bands. For m non-detections at a given band, in each random realisation we re-shuffle the input sample, preserving the same total m by allowing source duplication. We take the median of the resulting flux distribution as our formal stacked flux. The 1σ dispersion around this value is interpreted as the flux error. We propagate this uncertainty in quadrature with the standard deviation of the stacked map across 100 random positions within the cutout (after masking the central PSF). Although the latter component is typically sub-dominant relative to a bootstrapping dispersion, this conservative approach accounts for the strong fluctuations seen in low S/N stacked images, especially at low M_⋆.

As an example, Fig. 4 shows stacked cutouts in all IR/sub-mm bands at 0.8 < z < 1.2 (i.e. close to the median redshift of our sample) as a function of M_⋆. As expected from the tight MS relation that links M_⋆ and SFR in star-forming galaxies, stacks at low M_⋆ display lower S/N, despite having larger numbers of input sources.

Fig. 4. Stacked cutouts of NUVrJ–based SFGs at 0.8 < z < 1.2, as a function of M⋆ (left to right, expressed in log M⊙). Within each bin, we stacked only those sources with S/N < 3 at a given band. SCUBA 850 μm and AzTEC 1100 μm images are smoothed with a Gaussian kernel to ease the visualisation. Each cutout size is 8 × FWHM of the PSF, while for Spitzer-MIPS, we chose 13 × FWHM. Below each cutout, we report the corresponding S/N.

3.2. Correcting for clustering bias

The stacked flux densities calculated above can be biased high if the input sources are strongly clustered or very faint. This bias is caused by the greater probability of finding a source close to another one in the stacked sample compared to a random position. This generates an additional signal, as has extensively been discussed in the literature (e.g., Bavouzet et al. 2008; Béthermin et al. 2010, 2012, 2015; Kurczynski & Gawiser 2010; Bourne et al. 2012; Viero et al. 2013; Schreiber et al. 2015). Given the large number of stacked sources in each bin, the S/N is typically good enough to be able to correct for this effect, which becomes more prominent with increasing beam size (e.g., up to 50% for SPIRE images, see Béthermin et al. 2015). Here, we briefly describe our approach, referring the reader to Appendix A.2 of Béthermin et al. (2015) for a detailed explanation.

We model the signal from stacking as the sum of three components: a central point source with the median flux of the underlying population, a clustering component convolved with the PSF, and a residual background term (Eq. (2)). Following Béthermin et al. (2015), we attempt to separate these components via a simultaneous fit in the stacked images (Béthermin et al. 2012; Heinis et al. 2013, 2014; Welikala et al. 2016).

$$\begin{array}{c}\hfill S(x,y)=\phi \times \mathrm{PSF}(x,y)+\psi \times (\mathrm{PSF}\otimes w)(x,y)+\epsilon ,\end{array}$$(2)

where S(x, y) is the stacked image, PSF the point spread function, and w the auto-correlation function. The symbol ⊗ represents the convolution. The parameters φ, ψ, and ε are free normalisations of the source flux, clustering signal, and background term, respectively.

We parametrise the ‘clustering bias’ as bias = ψ/(φ + ψ), once we have verified that residuals (i.e. ε) are always consistent with zero within the uncertainties. We do not see any obvious M_⋆ or redshift dependence of the clustering bias. However, at a fixed wavelength, this can fluctuate significantly depending on the S/N of the input stacked image. For these reasons, we prefer to use an average clustering correction ⟨1 − bias⟩ for each band (see Table 2), drawn only from stacks with S/N > 3. For those images, we multiply the stacked flux by ⟨1 − bias⟩ at the corresponding wavelength. Only the MIPS 24 μm data are not shown, since their flux densities are not used for SED fitting in Sect. 3.3. Uncertainties on the clustering corrections were propagated quadratically with the stacked flux errors obtained in Sect. 3.1.

We stress that this method is suitable if the intrinsic source size is negligible compared to the PSF. This is especially true in SPIRE images, of which we show an example in Fig. 5. This refers to a specific bin at 0.8 < z < 1.2 and 11 < log(M_⋆/M_⊙) < 12, for which all stacks give good enough S/N. Particularly for SPIRE images, the clustering bias can make up to 40% of the total flux, and it can be recognised as a more extended and diffuse emission. However, for consistency, we extend this analysis to the full set of FIR/sub-mm data.

Fig. 5. Image decomposition of median stacks at 250, 350, and 500 μm, for a specific bin at 0.8 < z < 1.2 and 11 < log(M⋆/M⊙) < 12. From left to right: the stacked image is separated among a point source PSF, the clustering signal, and a residual background term, respectively. The colour scale is normalised to the maximum in each cutout for visual purposes. See Sect. 3.2 for details.

Lastly, the clustering-corrected source flux densities (S_stack) are combined with those of individual detections (S_i) with S/N > 3 in each band. If (m, n) is the number of stacked and detected sources, respectively, the weighted-average flux S_bin in a given bin is derived as follows:

$$\begin{array}{c}\hfill S_{\mathrm{bin}}=\frac{m\times S_{\mathrm{stack}}+{\sum }_{i=1}^n{S}_i}{n+m}·\end{array}$$(3)

Flux uncertainties were propagated in quadrature. For stacks with S/N < 3 in which we could not constrain the clustering correction, S_stack was set to the noise level of the stacked map (i.e. equal to its uncertainty). This way the weighted-average flux S_bin and its error are mainly driven by individual detections, for which the flux could be measured more accurately (Jin et al. 2018). If the combined flux S_bin has a S/N < 3, then it is set to three times the noise and interpreted as a 3σ upper limit.

3.3. Conversion to L_IR and SFR

This section illustrates how we fit the observed FIR/sub-mm SEDs to determine the total (8−1000 μm rest-frame) IR luminosity within each M_⋆ − z bin. To this end, we use the two-component SED-fitting code developed by Jin et al. (2018) (see also Liu et al. 2018). Briefly, this includes: three mid-infrared AGN torus templates from Mullaney et al. (2011); 15 dust continuum emission models by Magdis et al. (2012) that were extracted from Draine & Li (2007) to best reproduce the average SEDs of MS (14) or SB (1) galaxies at various redshifts. While the Draine & Li (2007) models were based on a number of physical parameters, the library of Magdis et al. (2012) depends exclusively on the mean radiation field ⟨U⟩ = L_IR per unit dust mass (M_d) and on whether the galaxy is on or above the MS. However, on the MS the average dust temperature strongly evolves with redshift (e.g., Magnelli et al. 2014) and directly enters M_d. Therefore, ⟨U⟩ and the SED shape both vary as a function of redshift, for which Magdis et al. (2012) empirically found as ⟨U⟩∝(1 + z)^1.15 up to z ∼ 2. More recently, Béthermin et al. (2015) revised the evolution of ⟨U⟩ with redshift out to z ∼ 4, using IR/sub-mm data in the COSMOS field, retrieving ⟨U⟩∝(1 + z)^1.8. Here, we adopt the set of 14 MS templates from Magdis et al. (2012), fit them to our data, and compare the ⟨U⟩−z trend with Béthermin et al. (2015).

The SED-fitting routine performs a simultaneous fitting using AGN and dust emission models, looking for the best-fit solution via χ² minimisation. In order to account for the typical photo-z uncertainty of the underlying galaxy population (at fixed M_⋆,z), each template is fitted to the data across a range of ±0.05 × (1 + ⟨z⟩) around the median redshift ⟨z⟩. The code keeps track of each SED solution and corresponding normalisation, generating likelihood distributions and uncertainties on for instance, L_IR, ⟨U⟩ and AGN luminosity, if any. We note that only FIR and sub-mm photometry (i.e. ignoring the MIPS 24 μm data-point) were used in the fitting procedure. This is done to avoid internal variations of the MIR dust features that cannot be captured by our limited set of templates (e.g., IR to rest-frame 8 μm ratio, IR8, Elbaz et al. 2011), which might affect the global FIR/sub-mm SED fitting. This optimisation clearly prioritises the FIR/sub-mm part of the SED, while not impacting the final L_IR estimates (e.g., Liu et al. 2018).

Figure 6 shows the best-fit star-forming galaxy template from the Magdis et al. (2012) library (green lines), as a function of M_⋆ (left to right, expressed in log M_⊙) and redshift (top to bottom). Red circles indicate the IR/sub-mm photometry, while downward arrows mark 3σ upper limits. The red dotted line is the best-fit AGN template from Mullaney et al. (2011), shown if significantly above 3σ. This is only found in the highest M_⋆ and redshift bin. Green dashed lines represent SEDs without FIR measurements and at z ≳ 1.5, for which the integrated L_IR is interpreted as 3σ upper limit (5/42 bins). Even though 24 μm has long been used as a proxy for L_IR, this is only accurate at z ≲ 1.5 (e.g., Elbaz et al. 2011; Lutz 2014). For this reason we still interpret as measurements the L_IR obtained from SEDs without FIR data, but only at z ≲ 1.5. That is the case for a few bins at the lowest M_⋆, in which the SED reproduces a-posteriori the 24 μm data-point. Globally, our stacking analysis yields robust L_IR estimates in 37/42 bins.

Fig. 6. Best-fit template obtained from a SED-fitting decomposition (green lines), as a function of M⋆ (left to right, expressed in log M⊙) and redshift (top to bottom). Red circles indicate the IR/sub-mm photometry (MIPS 24 μm, PACS 100−160 μm, SPIRE 250−350−500 μm, SCUBA 850 μm, and AzTEC 1.1 mm), while downward arrows mark the corresponding 3σ upper limits. The red dotted line is the best-fit AGN template, shown in the only bin where its significance is above 3σ. Green dashed lines represent SEDs without FIR measurements and at z ≳ 1.5, for which the integrated LIR is interpreted as 3σ upper limit (5/42 bins). MIPS 24 μm flux densities are not used in the fitting.

We find ⟨U⟩∝(1 + z)^{1.74 ± 0.18}, which is fully consistent with the revised ⟨U⟩−z trend of Béthermin et al. (2015): ⟨U⟩∝(1 + z)^{1.8 ± 0.4}. This test is reassuring since it confirms that one single z-dependent (or ⟨U⟩–dependent) MS galaxy template is fully able to reproduce the observed SED across a wide M_⋆ interval.

Given the tight correlation between L_IR and SFR, the IR data have been extensively used as proxy for SFR, assuming that most of galaxy star formation is obscured by dust (Kennicutt 1998; Kennicutt & Evans 2012). This is probably true inside the most massive star-forming galaxies (see Madau & Dickinson 2014 for a review). However, at decreasing M_⋆, galaxies become more metal-poor (e.g., Mannucci et al. 2010) and, thus, less dusty and obscured. In these systems the ultraviolet (UV) domain provides a key complementary view on the unobscured star formation (Buat et al. 2012; Cucciati et al. 2012; Burgarella et al. 2013).

On this basis, the comprehensive work by Schreiber et al. (2015) exploited IR-based SFRs (i.e. SFR_IR) and UV-uncorrected SFRs, in the deepest CANDELS fields, to calibrate the star-forming MS over an unprecedented M_⋆ (down to M_⋆ = 10^9.5 M_⊙) and redshift range (z ≲ 4). Since we carried out a similar analysis, it is worth checking whether the SFR estimates based on our IR stacking reproduce or not the MS of Schreiber et al. (2015).

For consistency, we need to collect the UV-uncorrected SFRs for our input sample. Hence, for each source we take the rest-frame NUV luminosity, L_NUV (2800 Å), given in the corresponding parent catalogue (Laigle et al. 2016 or Muzzin et al. 2013) to estimate the UV-uncorrected SFR following Kennicutt (1998): SFR_UV [M_⊙ yr⁻¹] = 3.04 × 10⁻¹⁰ L₂₈₀₀/L_⊙, already scaled to a Chabrier (2003) IMF. We verified that this conversion agrees well with more recent SFR_UV prescriptions (e.g., Hao et al. 2011; Murphy et al. 2011; see also Kennicutt & Evans 2012).

Within each M_⋆ − z bin, we simply take the median SFR_UV, and we add it to the SFR_IR corresponding to the stacked L_IR, calculated as SFR_IR = 10⁻¹⁰ L_IR/L_⊙ (Kennicutt 1998, scaled to a Chabrier 2003 IMF). Figure 7 displays our data in the SFR–M_⋆ plane, colour-coded by redshift over 0.1 < z < 4.5. At fixed M_⋆ and redshift, we show SFR_IR (circles) and the total SFR_IR + UV (open squares) for comparison. Downward arrows are 3σ upper limits scaled from L_IR. As can be seen, our data are in excellent agreement with the evolving MS relation at all redshifts (solid lines, Schreiber et al. 2015). While SFR_UV appears generally negligible compared to the total SFR, it becomes as high as SFR_IR towards low M_⋆ and low redshift (e.g., Whitaker et al. 2012, 2017). Our median values agree with Schreiber et al. (2015), even below M_⋆ ∼ 10^9.5 M_⊙, at which we extrapolate the MS relation due to lack of previous data (dashed lines). This test compellingly demonstrates that our L_IR can be deemed robust over the full M_⋆ and redshift interval explored in this work.

Fig. 7. SFR–M⋆ relation of the NUVrJ star-forming galaxies selected in this work, colour-coded by redshift over 0.1 < z < 4.5. At fixed M⋆ and redshift, SFRIR measurements are converted from the LIR obtained from IR/sub-mm SED-fitting (circles). Downward arrows indicate 3σ upper limits. For completeness, we also show the SFRIR + UV estimates by combining SFRIR with UV-uncorrected SFRs (open squares). Solid lines mark the evolving MS relation between SFRIR + UV and M⋆ at different redshifts (Schreiber et al. 2015). We observe an excellent agreement with the MS, even below M⋆ ∼ 109.5 M⊙ that relies upon a linear extrapolation from Schreiber et al. (2015) that is not constrained by previous data.

3.4. Radio stacking at 3 GHz

In this section, we describe the equivalent stacking analysis done with radio 3 GHz (Smolčić et al. 2017a) data, in order to derive average rest-frame 1.4 GHz spectral luminosities (L_1.4 GHz) in each M_⋆ − z bin. As is done for IR stacking (Sect. 3.1), we treat detections and non-detections separately.

The total flux densities of radio sources with 3 < S/N < 5 were taken from Jin et al. (2018) (see Sect. 2.3), while for brighter sources we matched their flux densities to those of the corresponding catalogues. The purpose of this approach is twofold: using the same published flux densities for S/N > 5 detections for consistency and avoiding to deal with the effect of side-lobes from bright sources in stacked images, which might complicate total flux measurements (see Appendix A of Leslie et al. 2020 for a discussion). In addition, radio detections might contain a substantial fraction of AGN, which is expected to increase at higher M_⋆ (e.g., Heckman & Best 2014). We will carefully deal with this issue in Sect. 4.2. At relatively faint flux densities (< 100 μJy), most of radio emission is thought to arise from star formation (Bonzini et al. 2015; Padovani et al. 2015; Novak et al. 2017; Smolčić et al. 2017a), though some AGN-related radio emission might still be contributing (e.g., White et al. 2015; Jarvis et al. 2016). It is for this reason that median stacking of both detections and non-detections (e.g., Karim et al. 2011; Magnelli et al. 2015) in deep VLA-COSMOS 1.4 GHz images ought to result in minimal radio AGN contamination. This alternative approach will be tested in Appendix D.

Within the UltraVISTA area analysed in this work, the 3 GHz rms (2.3 μJy beam⁻¹) fluctuates by less than 2% (Smolčić et al. 2017b). Indeed, we anticipate that no difference between median or rms-weighted mean stacking of non-detections is observed (see Appendix A and Leslie et al. 2020), as detailed below. It is for these reasons that we chose to perform a median stacking of non-detections. Individual detections will be added a posteriori via a simple mean weighted average, as done in Eq. (3).

Our stacking routine generates cutouts with size of 8 × FWHM_3 GHz (i.e. 6″ at 3 GHz), centered on the NIR position of each input galaxy. We acknowledge that an average offset of 0.1″ was found between 3 GHz (Smolčić et al. 2017b) and UltraVISTA positions (Laigle et al. 2016), which is half the size of a pixel. To account for this systematic offset, our routine performs sub-pixel interpolation and searches for the peak flux (S_peak) within ±1 pixel from the center of the stacked image. The peak flux uncertainty is estimated via bootstrapping 100 times, as done in Sect. 3.1. We take the median of the resulting flux distribution as our formal peak brightness. The 1σ dispersion around this value is interpreted as the corresponding error. We also measured the standard deviation across 100 random positions in the stack (masking the central beam of 0.75″). This gives less conservative errors compared to a bootstrap, but it is used to derive the formal rms of the stacked map.

Total flux densities (S_tot) are calculated by fitting a 2D elliptical Gaussian function to the median stacked image, using the IDL routine MPFIT2DFUN³. Given the typically high S/N (∼10 on average) reached in the central pixel, we leave size, position angle, and normalisation of the 2D Gaussian as free parameters. We verified that adopting a circular Gaussian or forcing the normalisation to the peak flux does not significantly affect any of our stacks. The total flux was calculated by integrating over the 2D Gaussian area A_gauss. The integrated flux error was computed by multiplying the peak flux error by $\sqrt{A_{\mathrm{gauss}}/A_{\mathrm{beam}}}$, where A_beam is the known beam area, and adding a known 5% flux calibration error in quadrature (Smolčić et al. 2017a). We recall that the peak flux error already incorporates the variance of the stacked sample via bootstrapping.

In order to assess whether our sources are clearly resolved, we follow the same criterion applied to VLA 3 GHz detections (Smolčić et al. 2017a) to identify resolved sources:

$$\begin{array}{c}\hfill \frac{S_{\mathrm{tot}}}{S_{\mathrm{peak}}}>1+{(6\times {S/N}_{\mathrm{peak}})}^{-1.44},\end{array}$$(4)

where S/N_peak is simply the peak flux divided by the rms of the image. This expression was obtained empirically to define an envelope containing 95% of unresolved sources, below such threshold. We find that 31 stacks out of 42 are resolved, according to Eq. (4). For these, total flux densities are on average 1.8× higher than peak flux densities. Similarly, Bondi et al. (2018) found 77% of VLA 3 GHz detected SFGs are resolved, and this fraction does not change significantly with M_⋆ (Jiménez-Andrade et al. 2019). Of the 11 bins with unresolved emission, 3 have S/N_peak < 3; these are all among the 5 bins without L_IR estimates from IR stacking (Sect. 3.3). Analogously to our treatment of the IR measurements, we discard all those 5 bins from the rest of our analysis.

For the stacks with resolved emission, we prefer to use their integrated flux from 2D Gaussian fitting as the most accurate estimate. Instead, for unresolved stacks we use the peak flux, consistent with the treatment of 3 GHz detections (Smolčić et al. 2017a). Fitting residuals are on average 3% of the total flux and always consistent with zero within the uncertainties. We validate this approach by reproducing the total flux densities of 3 GHz detections presented in Smolčić et al. (2017a) at S/N > 5 and in Jin et al. (2018) at 3 < S/N < 5, respectively (Appendix A).

Finally, we combined the radio stacked flux densities within each M_⋆ − z bin together with individual detections, following Eq. (3). The combined 3 GHz flux densities were first scaled to 1.4 GHz assuming of S_ν ∝ ν^α, with spectral index α = −0.75 ± 0.1 (e.g., Condon 1992; Ibar et al. 2009, 2010). This assumption is discussed in Appendix B. Lastly, 1.4 GHz flux densities were converted to rest-frame 1.4 GHz spectral luminosities (L_1.4 GHz), again assuming α = −0.75. Formal L_1.4 GHz errors were calculated by propagating the uncertainties on both combined flux and spectral index.

The various checks described in Appendix B prove our L_1.4 GHz robust across the full range of M_⋆ and redshift analyzed in this work. We note that our L_1.4 GHz estimates do not necessarily trace radio emission from star formation. Indeed, radio AGN are not yet removed at this stage, and they might be potentially boosting the L_1.4 GHz. This issue is addressed in Sect. 4.2.

4. The IRRC and the contribution of radio AGN

Using the median L_IR and L_1.4 GHz obtained from stacking, we study the evolution of the IRRC as a function of M_⋆ and redshift. Logarithmic uncertainties on each luminosity were propagated quadratically to get q_IR errors. Among the 42 M_⋆ − z bins analyzed in this work, 37 yield robust estimates of L_IR and L_1.4 GHz, while the remainder are discarded from the following analysis. Unsurprisingly, these latter 5 bins (three at 10⁸ < M_⋆/M_⊙ < 10⁹ and 1.8 < z < 4.5; two at 10⁹ < M_⋆/M_⊙ < 10^9.5 and 2.5 < z < 4.5) are among the least complete in M_⋆, as highlighted in Figs. 2 and 6. Therefore their exclusion partly mitigates the M_⋆ incompleteness of the remaining sample.

4.1. q_IR before removing radio AGN

Figure 8 shows the average q_IR as a function of redshift, colour-coded in M_⋆ (stars). For comparison, other prescriptions of the evolution of the IRRC are overplotted (black lines). Bell (2003) inferred the average IRRC in local SFGs, finding q_IR = 2.64 ± 0.02 (dotted line), with a scatter of 0.26 dex. Magnelli et al. (2015) studied an M_⋆-selected sample at z ≲ 2, and constrained the evolution of the far-infrared radio correlation (FIRC, parametrised via q_FIR⁴) across the SFR–M_⋆ plane at M_⋆ > 10¹⁰ M_⊙. From stacking IR and radio images, they parametrised the evolution with redshift of the FIRC as: q_FIR = (2.35 ± 0.08) × (1 + z)^{−0.12 ± 0.04}, where the normalisation is scaled to 2.63 in the q_IR space. More recently, Delhaize et al. (2017) exploited a jointly-selected sample of IR (from Herschel PACS/SPIRE) or radio (from the VLA-COSMOS 3 GHz Large Project, Smolčić et al. 2017a) detected sources (at ≥5σ) in the COSMOS field. Through a survival analysis that accounts for non-detections in either IR or radio, they inferred the evolution of the IRRC with redshift out to z ∼ 4 as: q_IR = (2.88 ± 0.03) × (1 + z)^{−0.19 ± 0.01}. While this trend appears somewhat steeper than that of Magnelli et al. (2015), we note that Delhaize et al. (2017) did not formally remove objects with significant radio excess, while Magnelli et al. (2015) performed median radio stacking to mitigate the impact of potential outliers such as radio AGN. Nevertheless, Delhaize et al. (2017) argue that the IRRC trend with redshift would flatten if applying a 3σ-clipping: q_IR = (2.83 ± 0.02) × (1 + z)^{−0.15 ± 0.01}, which becomes fully consistent with the findings of Magnelli et al. (2015).

Fig. 8. AGN-uncorrected qIR evolution as a function of redshift (x-axis) and M⋆ (colour bar). Errors on qIR represent the 1σ scatter around the median value, estimated via bootstrapping over LIR and L1.4 GHz. For comparison, other IRRC trends with redshift are taken from the literature (black lines): Bell (2003, dotted); Magnelli et al. (2015, dashed); Delhaize et al. (2017, dot-dashed). Our qIR values still include the contribution of radio AGN. See Sect. 4.1 for details.

When compared to the above literature, it is evident that our q_IR values lie systematically below other studies at M_⋆ > 10¹¹ M_⊙, while lower M_⋆ galaxies lie closer or slightly above them. In other words, our q_IR estimates seem to display a clear M_⋆ stratification, with the most massive galaxies having typically lower q_IR than less massive counterparts. As mentioned earlier in this work, we recall that our sample, at this point, contains some fraction of radio AGN, which might be boosting the L_1.4 GHz, particularly at high M_⋆ (see e.g., Best & Heckman 2012) where radio AGN feedback is known to be prevalent. Our L_IR estimates are, instead, corrected for a potential IR-AGN contribution (Sect. 3.3). Therefore, the net effect caused by including AGN is lowering the intrinsic q_IR. Selecting typical SFGs on the MS is expected to, however, reduce the incidence of powerful radio AGN expected in massive hosts since most radio AGN at z < 1 are found to reside in quiescent galaxies (e.g., Hickox et al. 2009; Goulding et al. 2014). It is for these reasons that we caution that Fig. 8 should be taken as the AGN-uncorrected q_IR. However, it is worth showing it to quantify how much q_IR changes after removing the radio AGN.

4.2. Searching for radio AGN candidates

In this section, we describe how we carried out a detailed study aimed at identifying potential radio AGN, removing them, and ultimately deriving the intrinsic q_IR trend that is purely driven by star formation. In our radio analysis, we combined individual radio-detections (above S/N > 3) with undetected sources via a weighted average (Eq. (3)). Contrary to stacking detections and non-detections together, this formalism enables us to characterise the nature of individual radio detections, that is, whether they show excess radio emission relative to star formation.

We accept an underlying assumption that radio-undetected AGN do not significantly affect any of our radio stacks. This is supported by the excellent agreement between mean- and median-stacked L_1.4 GHz of non-detections (Fig. A.2, bottom panel). Indeed, if the contribution of radio-undetected AGN were substantial, the corresponding mean L_1.4 GHz would be significantly higher than the median L_1.4 GHz from stacking. This assumption is further supported by the fact that the fraction of identified radio AGN is a strong function of radio flux density, and the sources we stack are, by construction, faint in the radio. Algera et al. (2020c) argue that below 20 μJy (at 3 GHz), the fraction of radio-excess AGN is < 10% (see also Smolčić et al. 2017b; Novak et al. 2018). We acknowledge that our assumption does not allow us to collect a complete sample of radio AGN, especially at high redshift where the fraction of radio detections notably drops (Fig. 3). Nevertheless, we will show that any residual AGN contribution does not change our conclusions.

We briefly summarise our next steps as follows. In Sect. 4.2.1, we explore the q_IR distribution traced by individual 3 GHz detections as a function of M_⋆ and redshift. First, we identify a subset of radio detections at M_⋆ > 10^10.5 M_⊙ that is representative of an M_⋆-selected sample. Then we decompose their q_IR distribution between AGN and star formation components (Sect. 4.2.2). This enables us to subtract potential radio AGN candidates and to calibrate the intrinsic best-fit IRRC with redshift at M_⋆ > 10^10.5 M_⊙ (Sect. 4.2.3). Then we extrapolate this calibration towards lower M_⋆ bins (Sect. 4.2.4), where a similar in-depth analysis was not possible due to radio-detections being strongly incomplete in this M_⋆ regime. Finally, the intrinsic (i.e. AGN-corrected) IRRC as a function of M_⋆ and redshift is presented in Sect. 4.3.

4.2.1. The q_IR distribution of radio detections

In order to study the q_IR distribution of 3 GHz detections, we need to calculate their average L_IR as a function of M_⋆ and redshift. For convenience, we refer the reader back to Fig. 2 (blue histograms) for visualizing the distribution of 3 GHz detections in the M_⋆ − z space. Out of 13 510 radio detections among our 37 bins, 8762 (65%) have a combined S/N_IR > 3, therefore reliable L_IR measurements from SED-fitting of FIR/sub-mm de-blended photometry (Jin et al. 2018). For the remainder of the sample, we stack again their IR/sub-mm images in all bands in each M_⋆ − z bin. Stacked IR flux densities are corrected for clustering bias and converted to L_IR following the same procedure adopted for the prior M_⋆ sample (Sect. 3.1). Median stacked L_IR are retrieved for the same 37/42 bins of the full parent sample, since a stacked S/N > 3 flux was obtained in at least one FIR/sub-mm band. Then, for each source we re-scaled its median stacked L_IR to the redshift and M_⋆ of that source (assuming the MS relation), in order to reduce the variance of the underlying sample within each M_⋆ − z bin. We verified that our stacked L_IR are always systematically below the 3σL_IR upper limits inferred from FIR/sub-mm SED-fitting (Jin et al. 2018). This ensures that our stacking analysis provides more stringent constraints on the L_IR of individual non-detections.

From this analysis, we are well-placed to explore the full q_IR distribution of 3 GHz detections at different M_⋆ and redshifts. Figure 9 shows q_IR as a function of redshift, split in six M_⋆ bins. Black dots mark individual 3 GHz detections, blue stars represent the q_IR obtained by combining detections and non-detections (same as in Fig. 8), while yellow squares are the stacks of non-detections only. In each panel we report the number of 3 GHz detected sources and the fraction of them with combined S/N_IR > 3. This fraction strongly increases with M_⋆, from 3.5% at 10⁸ < M_⋆/M_⊙ < 10⁹ to 78.3% at 10¹¹ < M_⋆/M_⊙ < 10¹², which implies that at the lowest M_⋆ nearly all q_IR estimates of radio detections rely upon IR stacking. This is because the 3 GHz detection limit sets a rough threshold in SFR (if radio emission primarily arises from star formation), therefore, it is biased towards high-M_⋆ galaxies because of the MS relation. Because of these potential biases, it is essential to identify the bins in which radio detections give us access to a representative sample of M_⋆-selected galaxies.

Fig. 9. Distribution of qIR as a function of redshift, split across increasing M⋆ bins. In each panel, we compare the qIR estimates of individual radio detections (black dots) with the median stacked values of non-detections (yellow squares) and with the weighted-average qIR of detections and non-detections together (from Eq. (3), blue stars). Green upward arrows indicate the corresponding threshold qIR, lim above which radio detections become inaccessible from an M⋆-selected sample. We select relatively complete M⋆ − z bins, in which at least 70% of radio detections have qIR below the corresponding threshold. This criterion identifies 13 bins (in red square brackets). Within these bins, the peak of qIR distribution (qIR, peak) is indicated with red open circles. In the two highest M⋆ bins, the best fitting trends with redshift are shown by the red dashed lines. Each panel reports the number of individual 3 GHz sources and their fraction with S/NIR > 3. See Sect. 4.2.1 for details.

Indeed, our purpose is to use single radio detections to calibrate a threshold that best distinguishes radio AGN from radio SFGs, as a function of M_⋆ and redshift. In order to extend this calibration to our full M_⋆-selected sample, we need to make sure that our derived trends are not affected by selection biases, namely, that the radio-detected sources we rely upon are fully representative of M_⋆-selected galaxies at a given redshift. To this end, within each bin, we compare the q_IR of single radio detections against a specific threshold (q_IR, lim), corresponding to the 3 GHz survey limit at a given M_⋆ and redshift (green upward arrows in Fig. 9). This threshold is proportional to the median stacked L_IR of the full SFGs sample, divided by the 3σ luminosity limit at 1.4 GHz, scaled from 3 GHz by assuming a fixed spectral index α = −0.75 (Appendix B). Specifically, q_IR, lim indicates the limiting q_IR at which a typical MS galaxy of a given M_⋆, z and L_IR drops below the 3 GHz detection limit, which translates into a lower q_IR limit. In other words, sources with q_IR > q_IR, lim lie within an M_⋆ range that is virtually inaccessible by our 3 GHz survey. Therefore, any measurement above that threshold is not representative of an M_⋆-selected sample. Conversely, radio detections below that threshold would be seen also in an M_⋆-selected sample of SFGs.

Within this framework, we consider as complete only those M_⋆ − z bins in which at least 70% of radio detections are below q_IR, lim. This cutoff enables us to narrow down the position of the mode of the q_IR distribution, leaving us with a total of 13 complete bins (at > 70% level) across the full sample. Unsurprisingly, they are preferentially located in high-M_⋆ galaxies, or at low redshift. These are delimited by a red segment in Fig. 9.

It is quite evident that the locus populated by radio detections tends to decline with redshift, at each M_⋆. However, this behaviour is far more pronounced at low M_⋆, and likely driven by selection effects. In fact, by definition the L_1.4 GHz of radio detections increases with redshift at all M_⋆, because 3 GHz sources are drawn from a flux-limited sample. On the contrary, the L_IR of radio detections behaves differently with M_⋆: at higher M_⋆ it is mostly based on IR-detected sources, while at lower M_⋆ it comes predominantly from IR stacking. At higher M_⋆, L_IR increases with redshift similarly to L_1.4 GHz, giving rise to a nearly flat q_IR locus. At lower M_⋆, instead, L_IR stands typically below the IR detection limit, thus not bound to a monotonic redshift increase. This effect causes an apparent decrease of q_IR with redshift that is mostly driven by the radio detection limit. Indeed, a similar trend can be seen in the green arrows, which move down in redshift at each M_⋆.

Since the single complete z-bin found at 10¹⁰ < M_⋆/M_⊙ < 10^10.5 is insufficient for us to constrain a redshift trend, we only consider the remaining 12 complete bins placed at M_⋆ > 10^10.5. For each of them, we identify the peak of the corresponding q_IR distribution of radio detections, namely, q_IR, peak (see red open circles in Fig. 9). We note that q_IR, peak represents the mode of radio detections, rather than the average that is, instead, potentially affected by underlying radio AGN (Sect. 4.2.2). Then we fitted a power law trend of q_IR, peak with redshift using the IDL routine MPFIT2DFUN, obtaining the best-fit expressions shown in Fig. 10.

Fig. 10. qIR distribution of 3 GHz detections in the two highest M⋆ bins (left and right panels), after correcting for the internal qIR–redshift dependence. Only sources within complete z-bins were considered at each M⋆. The total distribution (black histogram) was dissected among SF-dominated (blue) and AGN-dominated (red) radio sources, and fitted with two separate Gaussian functions (dot-dashed lines). While the SF population was fitted first, the AGN part was fitted in a second step from the total–SF residuals. The 1σ dispersion attributed to SF equals 0.20 and 0.23 dex at 1010.5 < M⋆/M⊙ < 1011 and 1011 < M⋆/M⊙ < 1012, respectively. Bottom panels: corresponding cumulative Gaussian fits, both normalised to unity. The vertical green dotted line marks the threshold that best separates between the SF and AGN populations (see Table 3), which rejects about 70% AGN and only 3−4% SFGs (open circles).

4.2.2. Identifying radio AGN at high M_⋆

After fitting the q_IR, peak trend with redshift for the two M_⋆ bins, we need to account for such redshift dependence while exploring the q_IR distributions of radio detections. To this end, we align the position of q_IR, peak in each redshift bin to match the best-fit redshift trend. This allows us to marginalise over the internal redshift trend, and merge all radio detection homogeneously within the same M_⋆ bin. The resulting redshift-corrected q_IR distribution is displayed in Fig. 10 for the two highest M_⋆ bins (left and right panel). Each total histogram (black) includes the contribution from both galaxy- and AGN-dominated radio sources. We proceed to dissecting into the two components as follows, leaving the discussion on how radio AGN affect the redshift trend in Sect. 4.2.3.

Assuming that the peak of the distribution is populated by radio-detected SFGs, and that the intrinsic q_IR distribution of SFGs is symmetric around the peak, we mirror the right-hand side of the observed q_IR distribution to the left side. This symmetric function is interpreted as the intrinsic q_IR distribution of SFGs (blue histogram). We fitted it with a Gaussian function, leaving normalisation and dispersion free to vary. The Gaussian fit yields a dispersion of 0.20 and 0.23 dex at 10^10.5 < M_⋆/M_⊙ < 10¹¹ and 10¹¹ < M_⋆/M_⊙ < 10¹², respectively (blue dot-dashed lines). The residual histogram (total-SF) is then fitted with a second Gaussian function (red dot-dashed lines), which parametrises the additional radio-excess population ascribed to AGN. We attempted to fit the AGN population with other non-Gaussian functions since the lowest q_IR tail is not perfectly reproduced with a Gaussian shape. However, we stress that our purpose is separating the two populations statistically and prioritise a clean identification of SFGs, while a proper characterisation of the shape of the AGN population is beyond the scope of this paper.

We note that our fitting approach relies on the assumption that q_IR, peak is entirely attributed to SF. Therefore, by mirroring and fitting the SF Gaussian first, it is possible that we might be underestimating the intrinsic relative fraction of radio AGN. We discuss this potential issue in Appendix C, though we anticipate that our main findings could only be reinforced if addressing this effect.

Another possible caveat of our approach lies in the assumption that IR-undetected sources are represented by a single stacked L_IR, rescaled, however, to the M_⋆ and redshift of each object based on the MS relation. However, we checked that the distribution of radio detections that are also IR-detected displays an average scatter of 0.22 dex, as for the full radio-detected sample shown in Fig. 10. This is because the vast majority of radio sources at M_⋆ > 10^10.5 M_⊙ are also individually detected at IR wavelengths (see Fig. 3). Therefore, taking a single stacked L_IR in each bin does not strongly impact the calibration of the SF locus.

Choosing the best dividing line between AGN and SF-dominated radio sources is a challenging, and somewhat arbitrary task. Moving the threshold to higher q_IR increases the purity of SFGs to the detriment of completeness, and vice-versa for a lower threshold. Here we attempt to reach low levels of cross-contamination between the SF and AGN populations, while keeping a high completeness of the SF population. For this reason, we checked the cumulative q_IR distribution drawn by the two Gaussian fits (AGN in red, SF in blue), as shown at the bottom of Fig. 10, each normalised to unity.

Four different thresholds (q_thres) were examined: (1) q_thres = q_peak–1σ; (2) q_thres = q_peak–2σ; (3) q_thres = q_peak–3σ; (4) q_thres = q_{cross, AGN = SF}. In this formalism, q_peak is still the peak of the SF population (blue Gaussian fit in Fig. 10), and σ its dispersion, while q_{cross, AGN = SF} represents the cross-over value at which the numbers of radio AGN and radio SFGs match each other. For each threshold, in Table 3 we report the cumulative fractions of SF and AGN populations lying below it. Qualitatively speaking, an ideal compromise consists of a low fraction of SF galaxies and a high fraction of AGN below the threshold.

This comparison highlights that the best trade-off between cross-contamination and completeness is given by the threshold q_thres = q_peak–2σ, in both M_⋆ bins. This method rejects about 70% of potential radio-excess AGN, and only 3−4% of SFGs, which we believe is quite acceptable. The offset from the corresponding q_peak value is on average 0.43 dex (that we deem robust in Sect. 4.2.3), which implies that our radio-excess AGN should have statistically at least 63% of their total radio emission arising from AGN activity.

4.2.3. Re-calibrating the radio AGN threshold

According to the threshold defined above, we removed radio-excess AGN from our 12 complete z-bins at M_⋆ > 10^10.5 M_⋆. Then we combined the remaining radio-detected SFGs with stacks of non-detections to compute the new L_1.4 GHz in those bins, which should be free from AGN contamination. We verified that the new L_1.4 GHz shifts the previously determined q_IR (blue stars in Fig. 9) upward by a certain amount. In those complete bins, we fitted the AGN-corrected q_IR with redshift, obtaining a significantly flatter relation than before, as shown in Fig. 11. This suggests that the steeper redshift trend seen before (Sect. 4.2.1) might be driven by radio AGN contamination, while the intrinsic redshift trend is significantly flatter, and possibly M_⋆ invariant.

Fig. 11. Same as Fig. 10, but normalising the peak to the flatter qIR − z trend calibrated after removing AGN (Sect. 4.2.3). This two-fold approach slightly improves the qIR decomposition, as highlighted by the larger cumulative fraction of radio AGN that are rejected below the qIR threshold (red open circles, 81% against the previous 70%).

To test the robustness of the newly derived q_IR − z trend, we again shift the q_IR measurements of individual detections by the offset from such a trend at each z-bin, and perform a second q_IR decomposition, as shown in Fig. 11. The Gaussian fit that parametrises star formation is nearly unchanged, with a dispersion of 0.21−0.22 dex in the two highest M_⋆ bins. The 2σ threshold below the peak is also very similar: 0.42 and 0.44 dex in the two bins (therefore we use an average Δq_AGN = 0.43 dex). Moreover, the cumulative histograms (bottom panels) underline that this latter decomposition rejects about 81% of radio AGN below the threshold, as opposed to 70% estimated in the first step (see red open circles in Figs. 10 and 11), while missing a comparable 3−4% of SFGs. This confirms the effective improvement led by our re-calibration of the SF locus in removing radio AGN.

As shown in the updated Fig. 12 (at M_⋆ > 10^10.5 M_⊙), subtracting radio AGN (red dots) based on this latter locus shifts all the median q_IR (blue stars) exactly on the fitted q_IR − z trend (blue solid lines). This agreement suggests that no further AGN subtraction is needed in those complete bins. Therefore, we can confidently assume that the new median q_IR coincide with the instrinsic peak of the SF population. Given the robustness of our analysis, we compute a weighted-average redshift slope among the two highest M_⋆ bins, by simply weighing each slope by the inverse square of its uncertainty. This way, we obtain an average slope of −0.055 ± 0.018, that is, not flat at a significance of 3σ.

Fig. 12. Distribution of qIR as a function of redshift and M⋆, after removing radio AGN (red dots). Symbols are the same as in Fig. 9, except for the median qIR estimates (blue stars), which are here re-calculated after removing radio-excess AGN (red dots). Fractions of radio-detected AGN and SFGs are reported in each M⋆-bin, as well as the fraction of AGN relative to the full M⋆ sample analysed in this work (in brackets). The blue and red solid lines denote the intrinsic IRRC of SFGs and the locus below which we classify radio sources as AGN (0.43 dex below the IRRC), respectively.

While at 10¹¹ < M_⋆/M_⊙ < 10¹² all z-bins (0.1 < z < 4.5) were used to constrain this trend, at 10^10.5 < M_⋆/M_⊙ < 10¹¹ we only used the first 5/7 z-bins (0.1 < z < 2.5). We now extrapolate the same relation also at 2.5 < z < 4.5, finding a good agreement with the median q_IR estimates.

The resulting fractions of radio AGN identified in the two highest M_⋆ bins should be quite representative of the overall incidence of radio AGN in these galaxies. This is suggested by the tightness of the SF Gaussian fit (σ ∼ 0.21−0.22 dex), that we interpret as the instrinsic scatter of the IRRC in these galaxies. Therefore, radio-undetected AGN that are not captured in our analysis, if any, are expected to be mostly composite (AGN+SF) radio sources whose total emission is predominantly arising from star formation processes.

It is worth noting that about 20% radio AGN still lie within our clean sample of SFGs, as shown in Fig. 11. As highlighted in Molnár et al. (2021), while this high-q tail of AGN is SF-dominated in the radio, it could add to the intrinsic scatter of the underlying pure SFG sample. Therefore, our inferred scatter of 0.21−0.22 dex could be slightly overestimated (see e.g., 0.16 dex in Molnár et al. 2021 for local SFGs), which would also be due to larger uncertainties on L_1.4 GHz and L_IR than in the Local Universe.

The fact that in both M_⋆ bins the subtraction of radio-excess AGN leads to a flattening of the q_IR − z trend might be also induced by a larger relative fraction of radio AGN with increasing redshift. As a sanity check, in both M_⋆ bins we split and decomposed the q_IR distribution of Fig. 11 separately at z < 1.2 and z > 1.2, examining the evolution of the relative fraction of radio AGN. Although we do find that radio AGN are slightly more prevalent at higher redshifts (i.e. on average from 12% at z < 1.2 to 18% at z > 1.2), we confirm that the dispersion of the SF population is redshift-invariant (∼0.20 dex), both before and after removing radio AGN. This implies that the relative offset between the AGN and SF loci is unchanged, therefore, our sample of > 2σ radio-excess AGN is globally preserved.

After removing those AGN, the flatter, yet declining q_IR evolutionary trend could be explained by residual radio AGN activity within the SF locus. We estimate the overall fraction of ‘pure’ SFGs to be 95% at z < 1.2 and 90% at z > 1.2. Such minimal AGN contamination is probably more important at higher redshifts because SFGs are intrinsically IR brighter, so the radio-excess contrast (at fixed L_1.4 GHz) is less evident. Therefore, we argue that any further correction for misclassified radio AGN would induce an even flatter q_IR trend with redshift.

Finally, our approach leads to the following fractions of radio-excess AGN. At 10^10.5 < M_⋆/M_⊙ < 10¹¹, radio AGN are 7.1% of all radio-detections and 2.2% of the full M_⋆ sample of SFGs. At 10¹¹ < M_⋆/M_⊙ < 10¹², radio AGN are 11.7% of all radio-detections and 6.0% of the full M_⋆ sample of SFGs (see Table 4). These numbers are consistent with the known prevalence of radio AGN in the most massive galaxies (e.g., Heckman & Best 2014; Hardcastle & Croston 2020). An increasing incidence of (X-ray) AGN activity with M_⋆ has also been reported in recent studies (Aird et al. 2019; Delvecchio et al. 2020; Carraro et al. 2020), possibly driven by the ability of the dark matter halo mass to regulate the amount of cold gas that trickles to the central black hole (Delvecchio et al. 2019).

Our empirical threshold identifies as radio-excess AGN sources with at least 63% of the total radio emission arising from AGN activity. Therefore, radio sources with lower, yet substantial AGN contribution could still be mis-classified as radio-SFGs (e.g., White et al. 2015, 2017; Wong et al. 2016). We attempt at quantifying this fraction by comparing our classification against ancillary VLBA data in the COSMOS field (Herrera Ruiz et al. 2017, 2018). This excellent dataset contains 468 VLBA sources detected at > 5σ, targeted from a pre-selected sample of VLA-COSMOS 1.4 GHz sources at S/N_1.4 > 5.5 (Schinnerer et al. 2010, 2864 sources). Since the brightness temperature reached by VLBA observations at about 0.01″ resolution exceeds 10⁶ K, detections are most likely to be radio AGN (Herrera Ruiz et al. 2017). Therefore, this sample provides an unambiguous method to test our source classification, though for a very tiny fraction of our sample with 1.4 GHz flux S_1.4 > 55 μJy, typically hosted in massive galaxies (M_⋆ > 10¹⁰ M_⊙). Out of 13 510 3 GHz radio detections among our 37 bins, we found only 189 VLBA counterparts within 0.5″ search radius. A fraction as high as 90% (170/189) were identified as ‘radio-excess AGN’ based on our recursive approach. The remaining 10% AGN mis-classified as SFGs from our approach are all IR-detected sources with typically high SFRs, which clearly reduces the apparent contrast between AGN- and SF-driven radio emission at arcsec scales. Although limited to a relatively bright and highly incomplete subsample, the comparison with VLBA data further demonstrates the reliability of our radio AGN identification method.

4.2.4. Extrapolating the SF-vs.-AGN loci at low M_⋆

We extrapolate the q_IR − z trend of non-AGN galaxies calibrated in the previous section towards less massive counterparts. As mentioned in Sect. 4.2.1, 3 GHz detections placed at M_⋆ < 10^10.5 M_⊙ are not representative of an M_⋆-selected sample. In particular, a galaxy of a given M_⋆ and redshift, with infrared luminosity L_IR of a typical MS galaxy would likely fall below the 3 GHz detection limit, as indicated by the green arrows in Fig. 9. Radio detections at these masses are therefore quite peculiar relative to the overall galaxy population.

This is further suggested in Fig. 9 by the q_IR offset between median measurements (blue stars) and individual radio detections (black dots). The latter lie systematically below the median q_IR, deviating more and more at lower M_⋆. For these reasons, we refrain from calibrating the IRRC directly on those radio detections. We prefer to use the median q_IR values as benchmark since they should be sensitive to a more representative sample of galaxies of that M_⋆.

We proceed as follows. Within each M_⋆ bin, the redshift trend of q_IR is extrapolated from that calibrated at higher M_⋆, namely q_IR ∝ (1 + z)^{−0.055 ± 0.018} (Sect. 4.2.3). It is only the normalisation that is left free to vary, in order to best fit the median q_IR. In other words, at M_⋆ < 10^10.5 M_⊙, we assume a constant q_IR − z slope. This approach is preferable to leaving also the slope as a free parameter since the small number of bins is insufficient for us to constrain the redshift trend as accurately as previously done with single detections. However, we stress that if we leave the slope free when fitting the q_IR in each M_⋆ bin, we always obtain slopes that are consistent between zero and −0.055, within 1σ uncertainties.

Following the iterative approach already tested at higher M_⋆, the best-fitting trend of q_IR with redshift enables us to identify radio AGN as sources lying > 0.43 dex below the best-fit SF locus. After subtracting those radio AGN, we re-calculate the weighted-average q_IR and search again for the best normalisation that fits the new AGN-corrected q_IR measurements with redshift. We repeat this procedure twice, that is, until all median q_IR are unchanged within the uncertainties, at each M_⋆. This condition sets the end of our recursion.

The final, AGN-corrected q_IR are shown in Fig. 12 for all M_⋆ bins (blue stars). This plot highlights the sample of radio-detected AGN that was removed (red dots) and the final SF locus (blue solid lines) that we eventually inferred after subtracting those AGN. The numbers of radio-detected AGN and SFGs are reported in each panel for convenience.

In most bins at M_⋆ < 10^10.5 M_⊙, the AGN-corrected q_IR measurements nearly coincide with those obtained from stacking non-detections alone (yellow squares). These latter values delimit the highest q_IR that could be reached if removing, by definition, all radio detections. The result of similarity between the two sets of q_IR measurements is due to a heavy subtraction of radio AGN from the sample of radio detections. Within the sample of radio detections, the fraction of radio AGN identified at M_⋆ < 10^10.5 M_⊙ increases with decreasing M_⋆. From the first to the fourth M_⋆ bin, these fractions are: 94.5%, 80.8%, 51.0%, and 14.6%, respectively. However, when compared to the size of our full M_⋆ sample in each bin, they drop to (in the same order): 0.4%, 0.6%, 1.1%, and 1.7%, respectively (see Table 4). These latter numbers are consistent with a decreasing incidence of radio AGN towards lower M_⋆ systems, following the trend constrained at M_⋆ > 10^10.5 M_⊙ (Sect. 4.2.3). Nevertheless, according to our analysis the vast majority of radio-detected dwarf galaxies (M_⋆ < 10^9.5 M_⊙, e.g., Mezcua 2017) in COSMOS are expected to be radio AGN.

Bearing this in mind, we note that the weighted average q_IR (blue stars) are as yet mostly driven by non-detections (yellow squares), which outnumber individual detections (dots) by a factor of > 100 at M_⋆ < 10^9.5 M_⊙. However, those few radio detections (mostly radio AGN, red dots) stand out from the stacks of non-detections (yellow squares) typically by over a factor of ten, up to one-hundred. As a consequence, after removing radio AGN at M_⋆ < 10^9.5 M_⊙, the new average q_IR still move upward by 0.2−0.3 dex.

4.3. The intrinsic IRRC evolves primarily with M_⋆

After correcting our combined L_1.4 GHz measurements for radio AGN contamination, we are finally able to examine the evolution of the intrinsic IRRC as a function of M_⋆ and redshift, as presented in Fig. 13. For each M_⋆ bin, we show the best-fit power-law trend, whose slope was directly inferred in Sect. 4.2.3 in the two highest M_⋆ bins (i.e. −0.055 ± 0.018). We verified that our median L_IR estimates are, instead, totally unchanged after removing radio-excess AGN, as expected given their minimal fraction relative to the parent M_⋆-selected SFG sample.

Fig. 13. Intrinsic (i.e. AGN-corrected) qIR evolution as a function of redshift (x-axis) and M⋆ (colour bar). The LIR estimates are the same reported in Fig. 8, while L1.4 GHz measurements have been re-calculated after excluding radio-detected AGN (Sect. 4.2). For comparison, other IRRC trends with redshift are taken from the literature (black lines): Bell (2003, dotted); Magnelli et al. (2015, dashed); Delhaize et al. (2017, dot-dashed) and their AGN-corrected version after removing 2σ outliers (triple dot-dashed lines).

The colour bar highlights a clear stratification of q_IR with M_⋆, with more massive galaxies showing systematically lower q_IR values. This behaviour was already seen in Fig. 8 before removing radio AGN, but here it suggests that some additional mechanisms unrelated to AGN activity might be boosting (reducing) radio emission in more (less) massive systems, relative to the IR.

For comparison, other IRRC trends with redshift are reported from Bell (2003, dotted line), Magnelli et al. (2015, dashed line), and Delhaize et al. (2017, dot-dashed line). Since Delhaize et al. (2017) did not remove radio-excess AGN, we also show their AGN-corrected relation by removing 2σ outliers (as reported in Delvecchio et al. 2018): q_IR ∝ (1 + z)^{−0.12 ± 0.01} (triple dot-dashed line). This trend is flatter than the previous one, more consistent with that of Magnelli et al. (2015) and more appropriate for a comparison with our approach. In the following, we examine the significance of the M_⋆ dependence at fixed redshift and we provide a multi-parametric fit as a function of both parameters.

Figure 14 shows the equivalent of Fig. 13 but projected on M_⋆, with redshift bins in different colours. Error bars on each q_IR were re-scaled by a factor of ${\sqrt{\chi }}_{\mathrm{red}}^2$ in each M_⋆ bin, to bring the reduced χ² of the corresponding q_IR − z fit to unity. It is quite evident that an M_⋆ dependence reduces substantially the scatter of the average q_IR around a single trend. To better quantify this, first we bootstrapped over the uncertainties of slope and normalisation obtained from each q_IR − z trend (see Table 4). Then, at each M_⋆, we interpolated the full range of bootstrapped IRRCs at z = 1, in correspondence of the 16th, 50th, and 84th percentiles. Interpolating at z = 1, besides being at roughly the median redshift of our sample, reduces the increasing divergence of each q_IR − z fit at lower or higher redshifts. This leaves us with the interpolated median q_IR(z = 1) as a function of M_⋆ (black open squares). Error bars indicate the uncertainty on the median value. The black dashed line marks the corresponding linear best-fit trend: q(M_⋆)_z = 1 ∝ ( − 0.124 ± 0.015) × log(M_⋆/M_⊙ − 10). This function yields a ${\chi }_{\mathrm{red}}^2$ = 0.87, with an M_⋆ slope close to that commonly found when fitting q_IR as a function of redshift (e.g., Magnelli et al. 2015), and significant at over 8σ. Even though the interpolated fit at z = 1 is purely indicative, this check suggests that M_⋆ might be the primary driver of the evolution of the IRRC across redshift.

Fig. 14. Distribution of AGN-corrected qIR as a function M⋆, colour-coded by redshift (stars). At each M⋆, open squares indicate the median qIR values at z = 1, obtained after propagating the uncertainties of slope and normalisation of the corresponding qIR − z fit and interpolating each at z = 1. These values were fitted with a linear function in log − log space (black dashed line).

Moreover, in order to incorporate the dependence of the IRRC on both M_⋆ and redshift simultaneously, we performed a multi-parametric fit in the three-dimensional q_IR − M_⋆ − z space. This yields the following best-fit expression:

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

where A = (1 + z) and B = (log M_⋆/M_⊙ − 10). The corresponding ${\chi }_{\mathrm{red}}^2$ = 0.90. The best-fit slopes with redshift and M_⋆ are significant at 2.9σ and 11σ levels, respectively. This further strengthens the need for a primary M_⋆ dependence, followed by a weaker and less significant redshift dependence. These numbers and confidence levels refer to the median trend. However, we acknowledge that when assuming a constant IRRC scatter of 0.21−0.22 dex across all M_⋆ galaxies, the weak co-dependence on redshift could be easily washed out. This dilution might also hide a mildly increasing redshift trend, which could be expected by Inverse Compton cooling of cosmic ray electrons (Murphy 2009). Nevertheless, the main argument of our analysis is to demonstrate how previously reported best-fitting IRRC trends with redshift are likely a red herring, whereas the M_⋆ (or a related proxy) is a better predictor of the average IRRC in SFGs.

The need for an additional M_⋆-dependence of the IRRC (in the form of L_radio–SFR) has also been highlighted in previous low-z studies (Gürkan et al. 2018; Read et al. 2018) and recently extended out to z ∼ 1 (Smith et al. 2021) using deep LOFAR-150 MHz data. A similar conclusion was independently reached in Molnár et al. (2021) when considering the dependence of the IRRC on galaxy spectral radio luminosity. To mitigate selection effects, they exploited a depth-matched sample of SFGs at z < 0.2. After performing a radio decomposition analysis in different bins of L_1.4 GHz, Molnár et al. (2021) report that q_IR decreases with increasing L_1.4 GHz. Assuming that radio emission comes predominantly from star formation, this is in line with our inferred M_⋆ dependence, since more massive SFGs are also brighter in radio (Leslie et al. 2020). This further corroborates the idea that the IRRC varies across different types of galaxies, at fixed redshift (but see e.g., Pannella et al. 2015 for an alternative interpretation). Therefore, we conclude that our results are in qualitative agreement with Molnár et al. (2021), who also demonstrate the implications of such a non-linearity for decreasing q_IR vs. z trends in the literature.

5. Discussion

The main result of this work is the finding that the IRRC primarily evolves with M_⋆, and only weakly with redshift (Eq. (5)). While the M_⋆ dependence has not been explored in detail so far, except in the Local Universe (e.g., Gürkan et al. 2018, see Sect. 5.1), for several years, much effort has been devoted to the understanding the mild, but significant decline of the IRRC with redshift from both an observational (e.g., Murphy 2009; Ivison et al. 2010a; Sargent et al. 2010; Magnelli et al. 2015; Delhaize et al. 2017; Calistro Rivera et al. 2017; Molnár et al. 2018) and a theoretical (Lacki & Thompson 2010; Schleicher & Beck 2013; Schober et al. 2016; Bonaldi et al. 2019) perspective. In Appendix D, we expand on the role played by various assumptions in deriving different IRRC trends presented in the literature. In this section, instead, we interpret our results and discuss the many implications of our findings in the context of the origin and evolution of the IRRC. In particular, we split the discussion in several sections, each focusing on a specific issue. First, we explore some physical interpretations of the origin of an M_⋆ and redshift-dependent IRRC (Sect. 5.1). We further investigate the possible evolution of the IRRC above the MS (Sect. 5.2). In addition, we present a discussion on the incidence of AGN activity (Sect. 5.3). Finally, we comment on the use of radio emission as a SFR tracer in light of our results (Sect. 5.4).

5.1. What drives the primary M_⋆ dependence?

Our main finding is that the IRRC decreases primarily with M_⋆, and only weakly with redshift. In particular, within the range 10⁹ < M_⋆/M_⊙ < 10¹², the median q_IR decreases by 0.25 dex (a factor of 1.8), at fixed redshift and with high significance (∼10σ, see Eq. (5)). This suggests that the dependence L_1.4 GHz − M_⋆ is steeper than the dependence L_IR − M_⋆ (i.e. the MS). To translate this result into the corresponding IR-radio slope, we take our best q_IR − M_⋆ relation (Eq. (5)) at fixed redshift, and assume for simplicity a linear MS between M_⋆ and SFR (i.e. L_IR). This yields $L_{\mathrm{IR}}\propto L_{1.4\mathrm{GHz}}^{0.90}$. In recent years, the deviation from a linear trend has been gaining increasing momentum thanks to several studies finding a similar sub-linear behaviour in the Local Universe ($L_{\mathrm{IR}}\propto L_{1.4\mathrm{GHz}}^{0.75-0.90}$, Bell 2003; Hodge et al. 2008; Davies et al. 2017; Brown et al. 2017; Gürkan et al. 2018; Molnár et al. 2021). This might challenge the idea of calibrating radio emission as a universal SFR tracer, as we discuss in Sect. 5.4.

Here, we explore some physical parameters behind this non-linearity that might induce an M_⋆-evolving q_IR similar to our findings. First, we discuss the possible role of a top-heavy IMF. Later, we test some radio synchrotron models (e.g., Lacki & Thompson 2010) by studying the relation between q_IR and SFR surface density.

5.1.1. The role of the IMF

We quantify whether a deviation from a canonical IMF slope (e.g., Chabrier 2003; n(M) ∝ M^−2.35 at 0.8 < M < 100 M_⊙) could justify an M_⋆-decreasing q_IR. In particular, we note that reprocessed IR light comes predominantly from stars with M > 5 M_⊙, while radio synchrotron emission comes from more massive stars with M > 8 M_⊙. We check whether a systematically flatter IMF in more massive galaxies could explain the observed decreasing q_IR.

A top-heavy IMF has been directly constrained only in massive early-type galaxies at z ∼ 0 (Cappellari et al. 2012) from the comparison between dynamical masses and optical light, but only proposed or indirectly inferred otherwise (e.g., Baugh et al. 2005; Hopkins & Beacom 2006; Davé 2008; van Dokkum 2008; Dabringhausen et al. 2009). To quantify the change of q_IR as a function of IMF slope, we integrate the IMF over the ranges 5−100 M_⊙ and 8−100 M_⊙, with varying IMF slope. The ratio between the two integrals is somewhat proportional to L_IR/L_1.4 GHz. However, we find only 8% variation of the integral ratio across the full range of slopes [−2.35, 0], as compared to 80% (i.e. 0.25 dex) q_IR variation across all M_⋆. In line with the conclusions of Murphy (2009), we argue that a top-heavy IMF in the most massive galaxies proves insufficient to explain the evolving q_IR with M_⋆.

5.1.2. The role of the SFR surface density (Σ_SFR)

The model proposed by Schleicher & Beck (2013) postulates that the magnetic field strength scales with SFR, boosting radio synchrotron emission during shocks or galaxy interactions (e.g., Donevski & Prodanović 2015; Tabatabaei et al. 2017). Because of the MS relation, this model implicitly predicts a net enhancement of radio emission with increasing M_⋆, as well as an increase with redshift due to higher gas density in galaxies. Related to this, the single-zone galaxy model of Lacki & Thompson (2010), which includes a detailed CR description, suggests a variation of q_IR as a function of SFR surface density (Σ_SFR), from ‘normal galaxies’ (Σ_SFR ≲ 0.06 M_⊙ yr⁻¹ kpc⁻²) to ‘starbursts’ (SBs, Σ_SFR ≳ 2−4 M_⊙ yr⁻¹ kpc⁻²). In particular, this model argues that q_IR slightly declines with Σ_SFR (up to Σ_SFR ∼ 1 M_⊙ yr⁻¹ kpc⁻²) due to the escape of CRe, generating a radio dimming especially in lower M_⋆ (or smaller size) galaxies. This effect is also expected to become more pronounced with redshift, due to IC scattering off the CMB that is expected to dominate over synchrotron cooling (Murphy 2009). Conversely, at Σ_SFR ≳ 1 M_⊙ yr⁻¹ kpc⁻², Lacki & Thompson (2010) invoke a ‘conspiracy’ of ionisation losses to balance spectral ageing, and additional synchrotron emission from secondary CRs, that together flatten q_IR with Σ_SFR and redshift.

Here, we test the above models by relating q_IR and average Σ_SFR measured in this work. These estimates were obtained by using the total SFR_IR + UV calculated from IR stacking and adding the dust-uncorrected UV contribution (Sect. 3.3). Galaxy sizes are drawn from median radio stacking of non-detections, carried out in Sect. 3.4 at each M_⋆ − z bin via 2D elliptical Gaussian fitting. Though these measurements do not include the contribution of single 3 GHz detections, they come from about 97% of all M_⋆-selected galaxies in our sample, hence, they should be statistically representative of their average radio properties. This approach implicitly assumes that radio emission encloses the total star formation of the host, which is quite plausible, especially in high-M_⋆ galaxies, where the dominant obscured SF traced by IR is also seen in the radio (e.g., Jiménez-Andrade et al. 2019). To scale angular sizes of θ_FWHM into effective radius (R_e, enclosing half of the total flux density), we assume that our galaxies follow a disc-like surface brightness profile (Sérsic index n = 1), as found for MS galaxies (e.g., Nelson et al. 2016). Under this assumption, the major-axis R_e, maj can be calculated as R_e, maj = θ_FWHM/2.43 (Murphy et al. 2017). Lastly, we take the circularised radius R_e = R_e, maj/$\sqrt{A_{\mathrm{r}}}$, where A_r is the axial ratio.

Figure 15 displays our median stacked 3 GHz size measurements (or upper limits) as a function of redshift and M_⋆. Error bars are obtained from the IDL routine MPFIT2DFUN. Upper limits are shown for unresolved stacks and correspond to the angular 3 GHz beam size (0.75″, grey dashed line), except for the highest M_⋆ bin at z < 0.5 that was convolved with a Gaussian kernel of 3″ FWHM (see Appendix B). Our R_e measurements are fully consistent with VLA 3 GHz sizes that were independently derived in a recent study of Jiménez-Andrade et al. (2019). The authors used the same VLA 3 GHz COSMOS images to construct a M_⋆-complete sample of 3184 radio-detected SFGs with M_⋆ > 10^10.5 M_⊙, most of which lie around the MS relation (Schreiber et al. 2015). The best-fitting R_e trend with redshift reported by Jiménez-Andrade et al. (2019) for MS galaxies (blue solid line, R_e ∝ (1 + z)^{−0.26 ± 0.08}) is broadly consistent with our evolutionary trend based on median 3 GHz stacks (orange solid line, R_e ∝ (1 + z)^{−0.18 ± 0.07}). Our slightly larger size measurements are likely due to radio-detected SFGs (Jiménez-Andrade et al. 2019) having a more centrally peaked surface brightness compared to our stacks (Bondi et al. 2018).

Fig. 15. Distribution of 3 GHz effective radius (in kpc) as a function of redshift and colour-coded by M⋆ (stars). Size measurements are taken from median stacked 3 GHz images of non-detections. Upper limits are given for unresolved stacks and correspond to the angular 3 GHz beam size (0.75″, grey dashed line). We observe a clear increase of Re with galaxy M⋆. The bins at M⋆ > 1010.5 M⊙ with resolved emission are fitted with a power-law redshift trend, which yields Re ∝ (1 + z)−0.18 ± 0.07 (orange solid line). A comparison study by Jiménez-Andrade et al. (2019) is shown (blue solid line) for 3 GHz detected SFGs at similar M⋆ in COSMOS, obtaining Re ∝ (1 + z)−0.26 ± 0.08.

We calculate Σ_SFR = SFR_IR + UV/2$\pi R_{\text{e}}^2$ (see e.g., Jiménez-Andrade et al. 2019) and show its relation with q_IR in Fig. 16, colour-coded by redshift (left panel) and M_⋆ (right panel). Empty symbols highlight 7/37 bins with unresolved 3 GHz stacked emission, which translates into a lower limit in Σ_SFR. By fitting only the 30 q_IR and Σ_SFR measurements, we obtain a significant anti-correlation similar in slope to that observed with M_⋆ (Sect. 4.3), marked by the black dashed line (q_IR ∝ ( − 0.13 ± 0.02) × logΣ_SFR). From the M_⋆ and redshift of each bin, we obtain a surrogate trend with rest-frame optical (5000 Å) sizes estimated indirectly from the van der Wel et al. (2014) scaling relation for SFGs (grey dotted line). We note that the mild difference between the two latter trends originates from the 2× smaller radio sizes compared to rest-frame optical sizes (Bondi et al. 2018).

Fig. 16. Evolution of qIR as a function of SFR surface density (ΣSFR = SFRIR + UV/2$\pi R_{\text{e}}^2$), colour coded by redshift (left panel) and by M⋆ (right panel). The SFRIR + UV estimates are taken from Sect. 3.3, while the effective radius Re is measured from stacked 3 GHz images via 2D Gaussian fitting. This plot shows a significant anti-correlation similar in slope to that observed with M⋆ (Sect. 4.3), marked by the black dashed line (qIR ∝ ( − 0.13 ± 0.02) × logΣSFR). For comparison, the best-fit trend with rest-frame optical (5000 Å) sizes estimated from van der Wel et al. (2014) scaling relation is shown (grey dotted line).

Since more massive galaxies are characterised by more compact star formation (Elbaz et al. 2011), the decreasing q_IR − Σ_SFR trend is linked to that with M_⋆. Nevertheless, unlike the trend with optical sizes, our Σ_SFR measurements are not bound to M_⋆ ‘by construction’, but rather measured from independent tracers (IR+UV and 3 GHz data). We thus stress that our proposed q_IR − Σ_SFR dependence is simply meant to be a more physical proxy for the observed M_⋆ dependence. At fixed M_⋆, the SFR surface density increases with redshift (left panel of Fig. 16), in qualitative agreement with our (weakly) decreasing q_IR trend.

Both the slope and significance of the q_IR − Σ_SFR relation are consistent with the values found between q_IR and M_⋆ (Sect. 4.3). We argue that the declining q_IR − Σ_SFR slope is primarily driven by the SFR, and only weakly by radio sizes. Indeed, at fixed redshift, the SFR(IR+UV) increases along the MS by a factor > 30 from 10⁹ to 10¹¹ M_⊙ (Fig. 7), while $R_{\text{e}}^2$ only increases by a factor of 1.5−2.5 in the same interval. Although this does not stand as conclusive evidence, our analysis seems to suggest that the larger SFR per unit area in more massive (and higher-z) galaxies might be driving the sub-linear behaviour of the IRRC.

The comparison with the model of Lacki & Thompson (2010) comes with a few caveats. First, we take a fixed spectral index α = −0.75 for all galaxies (which is supported by radio ancillary data in Appendix B), while Lacki & Thompson (2010) model a curved radio spectrum. Second, we label as ‘SBs’ those galaxies that lie > 4× above the MS (Sect. 5.2), while Lacki & Thompson (2010) identify them as Σ_SFR ≳ 2−4 M_⊙ yr⁻¹ kpc⁻². As we show in Fig. 16, the vast majority of our MS galaxies exhibit a value for Σ_SFR that is below the ‘SB’ threshold of Lacki & Thompson (2010).

That said, their model predicts a decreasing q_IR with Σ_SFR (see their Fig. 1) that steepens with redshift, then followed by a flattening (or a reversal) at Σ_SFR ≳ 1 M_⊙ yr⁻¹ kpc⁻². This behaviour is not clearly seen in our data, which instead display a smoothly declining q_IR − Σ_SFR trend and are nearly redshift-invariant. Our trend is consistent with the low-q values recently inferred by Algera et al. (2020a) in compact (R_e ∼ 1 kpc) and massive (M_⋆ > 10¹⁰ M_⊙) sub-millimetre galaxies at 1.5 < z < 3.5. Indeed, their average q_IR = 2.20 ± 0.03 lies close to the extrapolation of our best-fit q_IR − Σ_SFR trend at Σ_SFR ∼ 100 M_⊙ yr⁻¹ kpc⁻², further corroborating the relation between q_IR and SFR per unit area in SFGs.

Explaining low-q_IR SBs, a further fine-tuning in the model of Lacki & Thompson (2010), involves invoking ‘puffy SBs’ with larger disc scale height (h = 1 kpc, i.e. SMG-like) than canonical ‘compact SBs’ (h = 100 pc, i.e. ULIRG-like). Indeed, in puffy SBs, CRe can travel longer distances before escaping the galaxy, creating secondary hadrons that induce an extra boost of radio emission. However, this process should globally steepen the observed radio spectra (α ∼ −0.9:−1.0), due to bremsstrahlung and ionisation losses being weak with respect to synchrotron and IC losses. This prediction, again, has not been confirmed by our data (Fig. B.2).

Another speculative hypothesis could be linked to the amplification of the magnetic field strength at higher SFRs, which boosts radio emission in more massive galaxies along the MS (e.g., Tabatabaei et al. 2017); however, a fine-tuned balance between concomitant CRe losses and secondary CRe production is also required (Algera et al. 2020a).

In summary, our checks cannot firmly elucidate the main physical driver of the IRRC with M_⋆, but they seem to support an empirical link between q_IR and SFR surface density. Our findings do not seem to follow the q_IR flattening or spectral index variations with Σ_SFR predicted by models (e.g., Lacki & Thompson 2010; Schleicher & Beck 2013). Of course, our data do not have enough statistical power to discern all the underlying physical mechanisms and spectral variations that the model obviously addresses. We postulate that a more detailed data-vs.-model comparison would require depth-matched observations at multiple radio frequencies of massive compact galaxies.

5.2. Does the IRRC evolve above the MS?

We investigate the behaviour of the average q_IR above the MS. This is important for testing whether radio emission follows a similar enhancement as L_IR when moving above the MS – otherwise, q_IR is not shown to be a good tracer of starburstiness (i.e. offset from the MS). This issue is still highly debated. For instance, Condon et al. (1991) found that the most extreme ULIRGs at z ∼ 0 have higher q_IR and larger scatter compared to the MS population, which can be explained by flatter radio spectra due to free-free absorption (see also Murphy et al. 2013). On a different note, Helou et al. (1985) and Yun et al. (2001) do not report any significant deviation of q_IR in local SB galaxies, although they also observed a larger scatter for this population. More recently, Magnelli et al. (2015) found a mild (+0.2 dex) enhancement of q_FIR above the MS, though not significant. Such apparent tension is probably also due to different definitions of ‘starburst’ galaxies and different sample selections.

For sake of consistency with Magnelli et al. (2015), in this section, we define ‘SBs’ as galaxies with SFR > 4 × SFR_MS (e.g., Rodighiero et al. 2011), where SFR_MS corresponds to the SFR predicted by the MS (Schreiber et al. 2015), at each M_⋆ and redshift. Our measured SFR estimates come from IR+UV, as described in Sect. 3.3. However, following Carraro et al. (2020), we select as SBs only individually IR-detected galaxies (S/N_IR > 3) that meet the above criterion. This is because our stacked SFR_IR estimates are mostly dominated by MS galaxies, while the SB subsample is likely washed out in all median stacks. Especially at low M_⋆ and high redshift, this approach yields an incomplete SB sample due to galaxies being IR-fainter. In order to mitigate possible selection biases, we only focus on SB galaxies with M_⋆ > 10^10.5 M_⊙ and z ≲ 2.5. This interval is set to ensure that all SB galaxies (i.e. lying > 4× above the MS) lie above the limiting L_IR of Herschel PACS+SPIRE data in COSMOS (Béthermin et al. 2015) and are, thus, IR-detected. We further remove radio-excess AGN (pre-identified in Sect. 4.2) from the SB subsample of radio detections in order to consider only bona fide SFGs and fairly compare the AGN-corrected q_IR between the SB and MS populations. This leaves us with a sample of 554 SBs. As done for the full SFG sample, we performed median stacking at 3 GHz and combined the stacked signal with radio-detected SBs.

Figure 17 shows the resulting q_IR of SBs (circles) relative to the full SFG sample (MS+SB, stars) out to z ≲ 2.5, at M_⋆ > 10^10.5 M_⊙. For comparison, some previous IRRC trends are reported (black lines), as in Fig. 13. While some possible hints of (∼0.05 dex) higher q_IR in SBs could be present, these are consistent with MS analogues within 1σ in all bins. Therefore, this test suggests that q_IR evolves primarily with M_⋆, irrespective of whether a galaxy is on or above the MS.

Fig. 17. Comparison of qIR between our full SFG sample (MS+SB, stars) and the SB subsample (circles), as a function of redshift. To mitigate the incompleteness of an IR-based selection of SBs, we only show bins with M⋆ > 1010.5 M⊙. Black lines highlight best-fit IRRC trends from the literature for comparison. This test suggests that qIR evolves irrespective of whether a galaxy is on or above the MS.

Although the (lack of) evolution of q_IR above the MS is still debated, our decreasing q_IR − Σ_SFR trend (Fig. 16) would predict lower q_IR in SB than in MS galaxies, due to SBs being more compact. However, we note that our IR-detected SBs are both > 4× more star-forming and smaller in size (R_e ≲ 1 kpc at z < 2, see Jiménez-Andrade et al. 2019) than MS analogues. Therefore, both parameters add to boost Σ_SFR. Since our SBs appear mostly unresolved in 3 GHz stacks, we can only place lower Σ_SFR limits which prevent us from ruling out a possible flattening of q_IR at the largest SFR surface densities. Higher resolution radio observations of these objects would be crucial for testing such a behaviour.

On a side note, the sample of sub-millimetre galaxies for which Algera et al. (2020a) obtained an average q_IR = 2.20 includes SFGs within a factor of three from the MS relation, thus, they are not formally SBs. It might be possible that SB galaxies follow a different regime of q_IR, while our results predominantly reflect the behaviour of the MS population.

As mentioned in Sect. 5.1.2, Lacki & Thompson (2010) distinguish between ‘compact SBs’ (h = 100 pc, ULIRG-like) that are similar to local merging galaxies, and ‘puffy SBs’ (h = 1 kpc, SMG-like) with lower q_IR values, that are more common at high-z (Genzel et al. 2008). Nevertheless, a compact to puffy SB transition above the MS should be reflected to a steepening of their radio spectra indices (Lacki & Thompson 2010), although we are unable to discern it from our data. In this respect, Magnelli et al. (2015) did not report any significant spectral index variation above the MS. Therefore, we caution that a simple dependence of q_IR on the SF compactness might not be suitable for unveiling the physics behind the IRRC in SBs, which might be also be connected to the geometry of the SF regions or multiple mechanisms at play.

5.3. Observing widespread AGN activity in radio-detected dwarfs

A noteworthy implication raised from our radio AGN subtraction is the possibility of widespread AGN activity within radio-detected dwarf galaxies (M_⋆ < 10^9.5 M_⊙). As highlighted in Sect. 4.2.4 and Table 4, about 90% of radio-detected dwarfs are classified as radio AGN. This fraction drops down to only ∼0.5% relative to the full M_⋆ sample of dwarfs. Such huge difference suggests that radio-detected dwarfs are a quite peculiar and not representative subsample of these low-M_⋆ galaxies.

From an IR perspective, nearly all radio-detected dwarfs (> 99%) are completely undetected (S/N_IR < 3) at any IR/sub-mm band (Fig. 3). This is likely a natural effect that is due to the increasing incompleteness of IR selection towards low M_⋆ galaxies. From IR/sub-mm stacking, however, we obtain SFR_IR > 4× higher than the MS relation, placing these sources in the SB region (e.g., Rodighiero et al. 2011; Sargent et al. 2012). This might apparently support the scenario of a SF-driven origin of radio emission in dwarfs.

Nevertheless, on the radio side, these sources display on average lower L_1.4 GHz values than their more massive counterparts, but still 100 times larger than those obtained from median radio stacking of non-detections. This effect fully counter-balances the high starburstiness seen in the IR, causing an overall drop of q_IR in radio-detected dwarfs by over a factor of ten, with respect to the stacked population (see black dots relative to yellow squares in Figs. 9 and 12). These arguments let us suppose that radio-detected dwarfs are consistent with being AGN-dominated in the radio.

While there is broad consensus on the prevalence of radio AGN within massive galaxies (e.g., Heckman & Best 2014), in which AGN-driven feedback could hamper star formation, little is known about its incidence and impact in dwarfs. These systems are thought to host the pristine relics of the first black hole seeds, whose growth has been long believed to be disfavoured by SNa-driven feedback (e.g., Reines et al. 2013; Dubois et al. 2015; Mezcua et al. 2016; Marleau et al. 2017). However, there is mounting evidence that AGN feedback may also play a role at the low-mass end of the galaxy population.

From a theoretical perspective, cosmological simulations typically find that starbursting dwarf galaxies triggered by major mergers can be very frequent (Fakhouri et al. 2010; Deason et al. 2014). These events can induce widespread AGN feedback at low-M_⋆ regimes, that could help solve the so-called ‘too big to fail’ problem, whereby simulated dwarfs outnumber by several factors their observed counterparts (Garrison-Kimmel et al. 2013; Kaviraj et al. 2017). This excess number cannot be suppressed via SNa feedback alone, but through additional AGN feedback (Keller et al. 2016; Silk 2017; Koudmani et al. 2020).

To search for observational AGN signatures in dwarf galaxies, spatially-resolved emission line diagnostics (Mezcua & Domínguez Sánchez 2020), along with deep X-ray and high angular resolution radio observations have been used (e.g., Reines et al. 2011, 2014; Reines & Deller 2012; Mezcua et al. 2019). In the Local Universe, these campaigns led to the confirmation of on-going AGN activity in starbursting dwarf galaxies (Reines & Deller 2012). At higher redshifts, Mezcua et al. (2019) performed a statistical study of radio-detected dwarf galaxies at z < 3.4 using deep VLA-COSMOS 3 GHz data (Smolčić et al. 2017a). They isolated a sample of 35 bona-fide dwarf galaxies, which displayed radio jets powers and efficiencies as high as those of more massive galaxies. These studies argue that AGN feedback may be more common than previously thought, and potentially impactful for regulating galaxy star formation (Kaviraj et al. 2019). Our findings that most radio-detected dwarfs stand above the MS and display excess radio emission are, therefore, not surprising, and remain in broad agreement with the above literature.

To investigate the possible AGN nature of our radio-detected dwarfs, we stacked them using deep Chandra images (Civano et al. 2016; Marchesi et al. 2016) at different redshifts, finding no X-ray detection in any of them. However, the 3σ upper limits are about 5× higher than the level of X-ray emission predicted by star formation (Lehmer et al. 2016), which does not rule out their AGN nature.

5.4. Assessing radio-synchrotron emission as a SFR tracer

In this section, we discuss the link between the IRRC and SFR in galaxies. As mentioned in Sect. 3.3, the conversion from L_IR to SFR is quite accurate in massive galaxies, while towards less massive and less obscured systems, the UV may contribute as much as the IR to the global SFR. The observed correlation between L_IR and L_1.4 GHz is not, therefore, rigidly proportional to the SFR.

It is for this reason, we express q_IR through a slightly different formalism that accounts for the addition of dust-uncorrected UV emission, in order to study the connection between radio emission and total SFR (=SFR_IR + UV). We thus define the parameter q_SFR as:

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

where L_SFR is simply the SFR_IR + UV scaled back to spectral luminosity units. This formalism enables us to keep similar units as for q_IR, while switching from luminosity to total SFR.

We repeated the analogous q_SFR decomposition analysis at M_⋆ > 10^10.5 M_⊙ to calibrate the AGN-vs.-SF locus of radio detections (Sect. 4.2). Within the two highest M_⋆ bins, the best-fitting trend of q_SFR with redshift has slope −0.057 ± 0.002, which is strikingly similar to that inferred for q_IR (−0.055 ± 0.018, Sect. 4.2.3). Then we extrapolated such trend at lower M_⋆ bins (Sect. 4.2.4) and recursively removed radio AGN to derive the AGN-corrected IRRC.

Using the same approach as for Eq. (5), the multi-parametric fitting in the q_SFR − M_⋆ − z plane yields the following expression:

$$\begin{array}{c}\hfill q_{\mathrm{SFR}}(M_{\star },z)=(2.743\pm 0.034)\times A^{(-0.025\pm 0.012)}-B\times (0.234\pm 0.017),\end{array}$$(7)

where A = (1 + z) and B = (log M_⋆/M_⊙ − 10). Similarly to the fit in the q_IR space, the redshift dependence is weaker, and less significant than the M_⋆ dependence, which is unsurprisingly steeper than before. This suggests that radio emission drops considerably more than SFR in low-M_⋆ non-AGN galaxies. Reversing the argument, at fixed L_IR, radio emission underestimates the total SFR by a larger factor as compared to the IR light. The sub-linear trend $L_{\mathrm{IR}}\propto L_{1.4\mathrm{GHz}}^{0.90}$ that we inferred in our analysis (see also Bell 2003; Hodge et al. 2008; Davies et al. 2017; Brown et al. 2017; Gürkan et al. 2018) becomes even steeper when adding the UV contribution to L_IR, that is, SFR${}_{\mathrm{IR}+\mathrm{UV}}\propto L_{1.4\mathrm{GHz}}^{0.81}$. Such a radio deficit in the dwarf-galaxy regime could be possibly linked to shorter CRe scale heights (e.g., Helou & Bicay 1993; Lacki & Thompson 2010) or weaker magnetic fields (Donevski & Prodanović 2015; Tabatabaei et al. 2017) that are common in less dense SF environments.

Moreover, our adopted L_IR–SFR conversion does not account for the ‘cirrus’ emission associated with cold dust heated by old (> A-type) stellar populations, which might lower the intrinsic SFR at fixed L_IR. However, this effect is expected to contribute in low-sSFR galaxies, that is, at high M_⋆ and low redshift (e.g., Yun et al. 2001). Hence, we expect it (if any) to further flatten the q_SFR − z trend or to amplify the M_⋆ dependence of q_SFR.

In addition, we note that the lower efficiency in producing synchrotron emission in low-SFR, low-M_⋆ galaxies is already factored in recent synchrotron emission models of SFGs (e.g., Massardi et al. 2010; Mancuso et al. 2015; Bonaldi et al. 2019) based on empirical matching between local L_1.4 GHz and SFR functions. Therefore, our results reinforce the need for M_⋆-dependent, non-linear calibrations between radio-continuum emission and SFR, in order to develop successful observing strategies for targeting low-M_⋆ galaxies at radio wavelengths. These can be complemented with higher frequency observations that are more sensitive to thermal free-free emission as SFR tracer in high-redshift galaxies (e.g., Murphy et al. 2017; Penney et al. 2020; Van der Vlugt et al. 2021; Algera et al. 2020b).

These considerations are relevant in the context of the forthcoming SKA. In particular, the SKA mid-frequency receivers will be equipped with five bands, of which the SKA Band2 (0.95−1.76 GHz) will be the workhorse for radio-continuum based SFR measurements. Even the faintest and least massive galaxies in our sample will be routinely observed by SKA, probing diverse populations of SFGs (and composite AGN+SF objects). Our findings highlight that a detailed understanding of the physics behind the relation between radio synchrotron emission and SFR is fundamental for fully exploiting the unique SKA capabilities in terms of depth and angular resolution.

6. Summary and conclusions

In this work, we calibrate the IRRC of SFGs as a function of both M_⋆ and redshift, out to z ∼ 4. Starting from an M_⋆-selected sample of 413 78 galaxies SFGs selected via (NUV − r)/(r − J) colours in the COSMOS field, we leverage new de-blended IR/sub-mm data (Jin et al. 2018), as well as deep radio images from the VLA COSMOS 3 GHz Large Project (Smolčić et al. 2017a).

In each M_⋆ − z bin, we performed stacking of undetected sources at both IR (Sect. 3.1) and radio (Sect. 3.4) frequencies and combined the stacked signal with individual detections a posteriori to infer the average q_IR as a function of M_⋆ and redshift (Sect. 4.1). We develop a recursive approach for identifying and then subtracting radio-excess AGN in different M_⋆ and redshift bins (Sect. 4.2). This technique is calibrated on a (> 70%) M_⋆-complete subsample of 3 GHz detections at M_⋆ > 10^10.5 M_⊙ and extrapolated to the rest of the sample to infer the AGN-corrected IRRC (Sect. 4.3). Finally, we interpret our findings in the context of existing IRRC studies, from both models and observations. The main results of this work are listed below.

(1) The IRRC evolves primarily with M_⋆, with more massive galaxies displaying systematically lower q_IR. A secondary, weaker dependence on redshift is also observed. The multi-parametric best-fitting expression is the following: q_IR(M_⋆,z) = (2.646 ± 0.024) × (1 + z)^{( − 0.023 ± 0.008)} − (0.148 ± 0.013) × (log M_⋆/M_⊙ − 10). At fixed redshift, this trend translates into an IRRC of $L_{\mathrm{IR}}\propto L_{1.4\mathrm{GHz}}^{0.90}$, which corroborates the similar sub-linear behaviour reported in the literature (e.g., Bell 2003; Hodge et al. 2008; Gürkan et al. 2018). The typical scatter of the IRRC at M_⋆ > 10^10.5 M_⊙ is around 0.21−0.22 dex (a factor of 1.7), which is consistent with other studies (Yun et al. 2001; Bell 2003; Molnár et al. 2021) and roughly constant with M_⋆ and z.

(2) Our recursive approach for removing radio AGN enables us to statistically decompose radio-detected SFGs and AGN (Figs. 10 and 11) as a function of M_⋆ and redshift. Removing radio AGN substantially flattens the observed q_IR − z trend at M_⋆ > 10^10.5 M_⊙ to a nearly flat slope. This correction nicely aligns the mode q_IR of radio SFGs to the median stacked q_IR of the full M_⋆ sample of non-AGN galaxies. Therefore, we interpret the resulting AGN-corrected q_IR measurements as robust against further AGN removal. We acknowledge that residual radio AGN activity within radio-detected SFGs (10−20%) could be possible. Nevertheless, we expect this effect to, at most, further flatten out the evolution of q_IR with redshift and to induce an even steeper M_⋆ dependence, thus reinforcing our main findings.

(3) The fraction of radio AGN identified within the full M_⋆ sample strongly increases with M_⋆, spanning from 0.4% to 6% across the full range (Table 4), in agreement with previous studies (e.g., Heckman & Best 2014). However, when limited to 3 GHz detected sources, about 90% of radio-detected dwarfs (M_⋆ < 10^9.5 M_⊙) are radio-excess AGN. We argue this is likely a selection effect induced by our 3 GHz-limited being biased towards the brightest radio sources in such low-M_⋆ systems. We test the reliability of our radio AGN identification owing to available VLBA data of radio AGN (Herrera Ruiz et al. 2017), confirming the AGN nature for 90% of them.

(4) We examined the evolution of q_IR as a function of SFR surface density (Σ_SFR), as a proxy for M_⋆, finding a very similar trend both in slope and statistical significance. In agreement with recent observations of high-redshift dusty SFGs (Algera et al. 2020a), our results support a decreasing q_IR in MS galaxies towards higher Σ_SFR. Nevertheless, radio synchrotron models (e.g., Lacki & Thompson 2010; Schleicher & Beck 2013) predict a much stronger q_IR evolution with redshift, and Σ_SFR- (i.e. M_⋆-) dependent radio spectral indices, neither of which are seen in our data. Another possibility is related to magnetic field amplification in massive highly SFGs (Tabatabaei et al. 2017).

(5) We compare the average q_IR between MS galaxies and an M_⋆-complete subsample of SBs with SFRs > 4× above the MS (Sect. 5.2). Despite SBs being more compact than MS analogues (Jiménez-Andrade et al. 2019), we do not observe a significant difference in q_IR, which is visibly at odds with our expectations. According to radio synchrotron models (Lacki & Thompson 2010), a ‘conspiracy’ of different factors might induce a q_IR flattening at Σ_SFR ≳ 1 M_⊙ yr⁻¹ kpc⁻². Our findings do not seem to support this prediction. However, our current data do not allow us to discriminate between various model scenarios in this Σ_SFR regime. Alternatively, we postulate that SB galaxies might follow a different q_IR relation with Σ_SFR than MS analogues, in which multiple mechanisms could play a role.

(6) We verified that adding the UV dust-uncorrected contribution to the IR as a proxy for the total SFR would further steepen the q_SFR − M_⋆ trend, leaving the evolution with redshift unchanged. These findings imply that using radio-synchrotron emission as a SFR tracer requires M_⋆-dependent conversion factors. Finally, our results can be useful to make accurate calibrations for future radio-continuum surveys as SFR machines down to dwarf galaxy regimes, especially in the upcoming SKA era.

Acknowledgments

The authors are grateful to the referee for a detailed and constructive report that greatly helped us clarify the results and implications of this work. ID is supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 788679. ID thanks R. Gobat for useful discussions. MJJ acknowledges support from the UK Science and Technology Facilities Council [ST/N000919/1], the Oxford Hintze Centre for Astrophysical Surveys which is funded through generous support from the Hintze Family Charitable Foundation and a visiting Professorship from SARAO. SJ acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) under grant AYA2017-84061-P, co-financed by FEDER (European Regional Development Funds). DL acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 694343). IHW acknowledges support from the Oxford Hintze Centre for Astrophysical Surveys which is funded through generous support from the Hintze Family Charitable Foundation. RC acknowledges financial support from CONICYT Doctorado Nacional N° 21161487 and the Max-Planck Society through a Partner Group grant with MPA. JD acknowledges the financial assistance of SARAO. MN acknowledges support from the ERC Advanced Grant 740246 (Cosmic Gas). IP acknowledges financial support from the Italian Ministry of Foreign Affairs and International Cooperation (MAECI Grant Number ZA18GR02) and the South African Department of Science and Technology’s National Research Foundation (DST-NRF Grant Number 113121) as part of the ISARP RADIOSKY2020 Joint Research Scheme. SMR hereby acknowledged the financial assistance of the National Research Foundation (NRF) towards this research. JS acknowledges the funding from the Swiss National Science Foundation under Grant No. 185863. Opinions expressed and conclusions arrived at, are those of the author and are not necessarily to be attributed to the NRF. The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. We acknowledge use of the Inter-University Institute for Data Intensive Astronomy (IDIA) data intensive research cloud for data processing. IDIA is a South African university partnership involving the University of Cape Town, the University of Pretoria and the University of the Western Cape. The authors acknowledge the Centre for High Performance Computing (CHPC), South Africa, for providing computational resources to this research project.

References

Appendix A: Testing total radio flux densities

We validate our total flux estimation against individual detections taken from published VLA catalogues at 3 GHz. At S/N > 5 we used the catalogue of Smolčić et al. (2017a), while total flux densities at 3 < S/N < 5 were taken from Jin et al. (2018). After excluding the 67/10830 multi-component sources identified in Smolčić et al. (2017a), we calculate peak and total flux densities of each source, following the approach described in Sect. 3.4. In Fig. A.1, we display the comparison between total flux densities (dots), highlighting the corresponding median ratio at various intervals (squares). It is worth noting that Smolčić et al. (2017a, red) used the software BLOBCAT (Hales et al. 2012, 2014) to sum over all blobs identified in the 3 GHz image above a certain S/N cut. Therefore, it is best suited for non-Gaussian shapes. On the other hand, our approach described in Sect. 3.4 assumes a 2D elliptical Gaussian, with size, angle, and normalisation being free to vary. Despite the different techniques, we find a good agreement at S/N > 5, with a logarithmic offset of −0.028 dex and dispersion of 0.12 dex. At 3 < S/N < 5, total flux densities from Jin et al. (2018, blue) were computed via Gaussian PSF fitting, using a circular beam size of 0.75″. Despite the low S/N regime, we also observe a fair agreement, with an offset of −0.05 dex and dispersion of 0.24 dex. This check proves our total flux densities fully consistent with the published values of Smolčić et al. (2017a) and Jin et al. (2018) for individual 3 GHz detections down to S/N ∼ 3.

Fig. A.1. Comparison of total flux densities of 3 GHz detections between our procedure and catalogue flux densities, both at S/N > 5 (Smolčić et al. 2017a, red dots) and at 3 < S/N < 5 (Jin et al. 2018, blue dots). Squares highlight the median ratio at various intervals. The global offset and dispersion suggest a good agreement within the uncertainties down to S/N ∼ 3.

We further demonstrate that our choice of performing median stacking at 3 GHz, rather than rms-weighted mean stacking, does not impact our final L_1.4 GHz estimates. A comparison between median and mean L_1.4 GHz is presented in Fig. A.2. The top panel displays L_1.4 GHz from the combined flux of detections and non-detections (see Eq. (3)), while the bottom panel refers to the case of purely undetected sources. Colours indicate different M_⋆ bins. Only stacks in which peak flux densities have S/N > 3 are shown. No systematics is observed, at any M_⋆, between mean and median stacked L_1.4 GHz. This is consistent with White et al. (2007), who showed that, in the noise-dominated regime, the stacked median traces the population mean. Moreover, such excellent agreement confirms that the uniform 3 GHz sensitivity across the full map ensures that either stacking method can reliably recover the average flux of the underlying galaxy population.

Fig. A.2. Top panel: comparison between median L1.4 GHz (x-axis) and rms-weighted mean L1.4 GHz (y-axis) for combined 3 GHz detections and non-detections, following Eq. (3). Colours indicate various M⋆ bins. Bottom panel: same comparison, but referring to 3 GHz undetected sources only.

The fact that non-detections (bottom panel) display consistent L_1.4 GHz between mean and median stacking suggests that, if any, radio AGN do not dominate the total radio emission in our stacks. The same argument cannot be implicitly extended to the combined flux densities since these mean weighted-average flux densities could be biased towards fewer and brighter radio detections, which reduces the statistical weight of non-detections. Indeed, L_1.4 GHz of radio detections (top panel) are always > 3× larger than L_1.4 GHz of non-detections (bottom panel), despite the smaller numbers. This partly smooths over the initial fluctuations between mean and median stacking, thus delivering an even tighter agreement, as we observe.

Appendix B: Stacking ancillary VLA and MIGHTEE data

Here we perform a radio stacking analysis, as for 3 GHz data (Sect. 3.4), in order to check whether our 3 GHz based L_1.4 GHz are stable against different resolutions or spectral frequencies. We exploit VLA data from the 1.4 GHz Deep Project (Schinnerer et al. 2010) map. It covers 1.7 deg² with an angular resolution of 2.5″, reaching rms = 12 μJy beam⁻¹ in the central 50″ × 50″. A total of 2864 sources were blindly extracted down to S/N > 5. In addition, we make use of 1.3 GHz MIGHTEE (Jarvis et al. 2016; Heywood et al., in prep.) data. The MIGHTEE images formally reach 2.2 μJy beam⁻¹ at 8.4″ × 6.8″ resolution over 1 deg² in the MIGHTEE early science data, but the effective depth is limited by confusion (∼5.5 μJy beam⁻¹ in the central part).

Source flux densities in VLA 1.4 GHz and MIGHTEE 1.3 GHz maps were re-extracted, using K_s + MIPS 24 positional priors. While the angular resolution at VLA 1.4 GHz is high enough to yield a negligible fraction of overlapping priors within the beam, MIGHTEE data suffer from blending issues. To this end, MIGHTEE flux densities were de-blended as in Jin et al. (2018) down to the 3σ level. Then, individual S/N > 3 detections were removed from the original image, and we used the residual map for stacking 1.3 GHz non-detections. Of course, only sources within the MIGHTEE area (central 1 deg²) were stacked, containing roughly half of the sample size used for VLA stacking.

The stacking analysis follows the same reasoning and assumptions presented in Sect. 3.4. Stacked MIGHTEE flux densities are measured in the central pixel, that is assumed to trace the total flux. VLA 1.4 GHz peak flux densities were, instead, scaled to total flux densities as done at 3 GHz. Nonetheless, a different, yet empirical relation was adopted to identify resolved sources at VLA 1.4 GHz, calibrated on 1.4 GHz detections (Schinnerer et al. 2010): S_tot/$S_{\mathrm{peak}}>0.{35}^{-11/({S/N}_{\mathrm{peak}}^{1.45})}$. Because of the larger beam size compared to 3 GHz, we find fewer resolved stacks (17/25). Total flux densities were converted to rest-frame L_1.4 GHz assuming α = −0.75 ± 0.1, that was propagated along with flux errors to deliver reliable L_1.4 GHz uncertainties. Upper limits at 3σ were assigned for S/N < 3 stacks.

Figure B.1 shows stacked cutouts at 0.8 < z < 1.2 at VLA 3 GHz (top), VLA 1.4 GHz (middle) and MIGHTEE 1.3 GHz (bottom) data, as a function of M_⋆ (increasing from left to right). While stacking at 3 GHz delivers S/N > 3 in 39/42 bins, only 25/42 and 29/42 have S/N > 3 in VLA 1.4 GHz and MIGHTEE 1.3 GHz stacked images, respectively. For VLA 1.4 GHz, the small number of bins is attributed to shallower than 3 GHz observations. For MIGHTEE, instead, this is probably induced by the confusion-limited signal in the stacks due to the larger MeerKAT primary beam at 1.3 GHz (e.g., Mauch et al. 2020). Nevertheless, because VLA 3 GHz data are much less sensitive than MeerKAT to large-scale radio emission, total radio flux densities might be underestimated at 3 GHz. This issue can be, however, especially relevant at low redshift (z < 0.5) and for resolved or multi-component radio sources (e.g., Delhaize et al. 2021). In fact, a visual inspection of the median 3 GHz stacks of non-detections does not clearly reveal any missing flux in the residual images at the scales of the MIGHTEE beam, except in the bin at z < 0.5 and 10¹¹ < M_⋆/M_⊙ < 10¹². To quantify this effect, we convolved all the original 3 GHz stacked cutouts with a Gaussian kernel of 3″ FWHM, re-calculating the total flux densities and comparing them with the previous measurements. The reason why this specific beam width was chosen is that beyond 3″ it has already been shown that no significant missing flux is recorded at 3 GHz (see Table 2 in Delhaize et al. 2017). Of course, this convolution drastically reduces the global S/N of the final stacks, leaving us with S/N > 3 in only 16/42 bins (as opposed to 39/42 before). However, only the bin at the lowest z and highest M_⋆ displays, on average, a total flux density that is larger by 0.3 dex, while the other bins show consistent estimates within the uncertainties. Since no extra flux is visible in the new residual image, we replaced the total flux density of that single bin with the 3″ convolved value and used this value in the rest of our analysis. In any case, we stress that the final L_1.4 GHz obtained by combining both detections and non-detections is unchanged since the fraction of radio detections is about 56% at z < 0.5 and M_⋆ > 10¹¹ M_⊙ (see Fig. 2), thus, washing out the difference in the stacked flux. As a consequence, this effect has no impact on the rest of our analysis. In addition, we emphasise that any extra missing 3 GHz flux at low redshift would further strengthen our final redshift-invariant IRRC (Sect. 4.3). This motivates our choice of using primarily VLA 3 GHz images for our analysis.

Fig. B.1. Stacked cutouts of our sample at 0.8 < z < 1.2, as a function of M⋆ (left to right, expressed in log M⊙). Only individual undetected sources (S/N < 3) are stacked. Top, middle and bottom rows: VLA 3 GHz, VLA 1.4 GHz and MIGHTEE 1.3 GHz data, respectively. Each cutout size corresponds to 8 × FWHM of the beam. Below each cutout we report the corresponding S/N of the total stacked flux.

Considerations on the radio spectral index

We briefly discuss and test our assumption of taking a single spectral index α = −0.75 by comparing L_1.4 GHz estimates independently inferred from stacking the three above datasets. In Fig. B.2, we compare L_1.4 GHz obtained from 3 GHz stacks (x-axis) and ancillary radio stacks (y-axis), using either VLA (1.4 GHz, circles) or MIGHTEE (1.3 GHz, squares) data, colour-coded by M_⋆. Downward arrows with open symbols mark 3σ upper limits where the stacked S/N < 3. We find a good agreement between all the various datasets, suggesting that using a single α = −0.75 is a reasonable assumption across the full M_⋆ range explored in this work. As a sanity check, the median spectral index traced by VLA 3 GHz and MIGHTEE 1.3 GHz individual detections is −0.77, in agreement with our assumption. However, we prefer to adopt a fixed α = −0.75 in our calculation to treat both radio detections and non-detections in a self-consistent manner.

Fig. B.2. Comparison between rest-frame 1.4 GHz spectral luminosity L1.4 GHz obtained from 3 GHz stacks (x-axis) and ancillary radio stacks (y-axis) using VLA (1.4 GHz, circles) and MIGHTEE (1.3 GHz, squares) data. We assumed a single spectral index α = −0.75 to scale flux densities from 3 GHz to 1.4 GHz. Colours indicate different M⋆ ranges. Downward arrows with open symbols mark 3σ upper limits if S/N < 3. The broad agreement between the various datasets suggests that using a single α = −0.75 is a reasonable assumption across the full M⋆ range explored in this work.

Magnelli et al. (2015) measured the average spectral index exploiting VLA 1.4 GHz and GMRT 610 MHz data for an M_⋆-selected galaxy sample. They found that the observed 610 MHz−1.4 GHz slope, which probes closer to the rest-frame of 1.4 GHz than our 3 GHz data, does not seem to change with M_⋆ or SFR, at least out to z ∼ 2. More recently, Calistro Rivera et al. (2017) exploited Low Frequency Array (LOFAR) data at 150 MHz in the Boötes field, out to z ∼ 2.5. Interestingly, they observed a spectral flattening of the radio SED of SFGs in the observed range [150 MHz−1.4 GHz] (see also Read et al. 2018; Gürkan et al. 2018). However, they argue that this feature should not affect the k-correction for the rest-frame 1.4 GHz luminosities L_1.4 GHz. Therefore, these studies provide mounting evidence that using a single power-law spectral index α = −0.75 at our frequency is a reasonable assumption.

Appendix C: Impact of a different radio AGN-vs.-SFG fitting approach

We discuss a potential caveat related to our AGN-vs.-SF decomposition presented in Sect. 4.2.2. Specifically, our procedure relies on the assumption that the mode of the observed q_IR distribution (q_IR, peak) of radio detections is entirely attributed to SF. Though this is supported by a number of previous studies arguing that radio AGN are a sub-dominant population in the sub-mJy regime (e.g., Padovani et al. 2015; Smolčić et al. 2017a; Novak et al. 2018; Ceraj et al. 2018; Algera et al. 2020c), the contribution of radio-faint AGN to q_IR, peak might not be negligible. If this is the case, by mirroring and fitting the SF Gaussian first, it is possible that we are underestimating the true fraction of radio AGN relative to SFGs. To quantify this potential issue and test how much it would affect our final M_⋆-dependence of q_IR, here, we follow a different approach.

The observed q_IR distribution is fitted with two Gaussian functions simultaneously, which parametrise the contribution of SFGs and radio-excess AGN. Contrary to Sect. 4.2.2, we do not set the SF peak to q_IR, peak, but we leave it free to vary along with the dispersion and normalisation for both functions. In this simultaneous fitting we give equal input weights to all bins, regardless of the number of sources in each. This approach is thus expected to return a rather conservative AGN contribution relative to SFGs.

The results are shown in Fig. C.1 for the two highest M_⋆ bins. As in Fig. 10, the best-fit SF (blue) and AGN (red) Gaussians head up to reproduce the total distribution (black). However, we clearly notice two main differences compared to the previous approach. Firstly, the AGN distribution is far broader than the SF distribution in both M_⋆ bins. Secondly, the relative fraction of radio AGN that we mis-classify as SFGs (red tail at > q_peak − 2σ) is as high as 40−70% – hence, it is much higher than the 30% obtained in Sect. 4.2.2 when fitting and mirroring the SF part first. This is clearly displayed by the cumulative AGN fraction in the bottom panels. Instead, the relative fractions of ‘pure’ SFGs above the 2σ threshold are about 80% at 10^10.5 < M_⋆/M_⊙ < 10¹¹ and 90% at at 10¹¹ < M_⋆/M_⊙ < 10¹².

Fig. C.1. Same as Fig. 10, but fitting the total qIR distribution of 3 GHz detections (black histogram) simultaneously with a SF Gaussian (blue) and an AGN Gaussian (red dashed). The 1σ dispersion has remained unchanged to about 0.21 dex. Bottom panels: corresponding cumulative Gaussian fits, both normalised to unity. The vertical green dotted line marks the 2σ threshold of 0.42 dex below which we consider a source as radio-excess AGN. As opposed to Fig. 10, this cutoff removes only 30−50% of the total radio AGN population. However, we estimate this effect to be more prevalent in lower M⋆ galaxies. This implies that accounting for such mis-classified radio AGN would likely strengthen our final M⋆-dependent qIR.

Despite the lower level of purity of the SFG population, we emphasise that the main results of this paper are quite robust against the AGN versus SF fitting procedure. Indeed, both peak and dispersion (∼0.21 dex) of the SF population are essentially unchanged, the peak being identical to q_IR, peak and the dispersion reaching ∼0.21 dex. Therefore, it is reasonable to assume that the mode of the observed q_IR distribution is attributed to radio-detected SFGs. Related to this, the threshold q_peak − 2σ is still equal to 0.42 dex, implying that roughly the same exact sources as in Sect. 4.2.2 would be identified as radio-excess AGN. This agreement demonstrates that our recursive radio AGN removal would lead to the same final IRRC, regardless of the assumed shape of the AGN distribution.

If we were able to statistically remove the underlying radio AGN contribution within the SF population (though impossible with the present data), this would systematically increase q_IR by a larger amount towards lower M_⋆ galaxies. Indeed, at 10¹¹ < M_⋆/M_⊙ < 10¹² the radio AGN distribution is clearly broader, but far more offset than at 10^10.5 < M_⋆/M_⊙ < 10¹¹, thus at higher M_⋆ the two populations are more distinguishable. As a consequence, we argue that a proper correction for such an effect would further exacerbate the M_⋆ stratification of q_IR reported in this work.

Appendix D: Differences compared to the literature

Our best-fit relation of q_IR as a function of M_⋆ and redshift (Eq. (5) in Sect. 4.3) is fully consistent with the average q_IR value measured in local SFGs (i.e. 2.64 in Bell 2003) for a typical galaxy with M_⋆ ∼ 10¹⁰ M_⊙. At higher redshifts, instead, our average q_IR measurements follow flatter evolutionary trends compared to previous studies (Fig. 13), while the best-fit normalisation appears broadly consistent with the literature only at M_⋆ > 10^10.5 M_⊙. In order to interpret these differences in a quantitative fashion, we identify three key points that combined differentiate our approach from that adopted in the previous literature: (i) removing radio AGN via a recursive approach in each M_⋆ and redshift bin; (ii) exploiting an M_⋆-selected sample of SFGs; (iii) binning the derived q_IR as a function of both M_⋆ and redshift. To test our results against different techniques, we expand on each of these aspects below.

D.1. Radio AGN subtraction

In Sect. 4.2, we performed a recursive subtraction of radio AGN as a function of M_⋆ and redshift, carefully calibrated on high-M_⋆ galaxies, and then extrapolated to lower M_⋆ analogues. However, other studies followed alternative approaches to discard radio AGN when deriving the intrinsic IRRC. For instance, Magnelli et al. (2015) performed median stacking of both radio detections and non-detections out to z ∼ 2. This method strongly reduces the contribution of a few bright outliers, assuming that the bulk radio population is made of SFGs. This assumption is quite reasonable, since Magnelli et al. (2015) started from an M_⋆-selected sample, of which radio detections make a negligible fraction.

We compare our q_IR with mock measurements obtained by following the stacking method of Magnelli et al. (2015), but applied to the sample used in our work. Figure D.1 displays the final L_1.4 GHz estimates that we obtained after removing radio AGN (x-axis) against those derived from median radio stacking (Magnelli et al. 2015, y-axis). We note that our L_IR estimates and Magnelli et al. (2015)’s were instead calculated through a fully consistent approach, therefore only a difference in L_1.4 GHz might lead to systematics in the final q_IR trends. The colour bar highlights the average M_⋆ of each bin. Out of 37 bins analyzed in this work, 35 yield a S/N > 3 from median 3 GHz stacking (circles), while 3σ upper limits are shown for the remaining bins (downward arrows). This comparison clearly reveals a very good agreement between final 1.4 GHz luminosities, with all measurements being consistent within the uncertainties. Despite the different approaches, the agreement extends down to dwarf galaxies, supporting the AGN nature of most radio-detected sources (Sect. 4.2.4). A possible (though not significant) deviation of ∼0.1 dex might be present at the highest M_⋆, with our measurements returning slightly higher L_1.4 GHz measurements than those of Magnelli et al. (2015). This might be ascribed to the contribution of radio-detected SFGs to our weighted average L_1.4 GHz, since they make a substantial fraction of the M_⋆-selected sample in that M_⋆ bin (∼45%, Table 4). Therefore, this test proves our radio AGN subtraction broadly consistent with a totally independent approach.

Fig. D.1. Comparison between AGN-corrected L1.4 GHz from this work (x-axis) and median L1.4 GHz obtained from stacking detections and non-detections together (Magnelli et al. 2015, y-axis). Different M⋆ ranges are colour-coded, while downward arrows mark 3σ upper limits for 2/37 bins. Despite these different approaches, we notice a very good agreement in all bins, that strengthens the reliability of our recursive AGN subtraction.

D.2. Sample selection and binning

An additional aspect worth testing is whether different sample selections lead to distinct IRRC trends. We started from an M_⋆-selected sample of SFGs based on K_s-band priors, that typically reaches much deeper than any infrared or radio survey, compared to an average galaxy SED. A rare exception is represented by very high-redshift (z > 4) or heavily dust-obscured systems, which are visible only in IRAC (e.g., Davidzon et al. 2017) or deep ALMA imaging (e.g., Franco et al. 2020). For this reason, studies that derived the IRRC based on exclusive or joint samples of radio/IR detections, are partly biased against low-M_⋆ galaxies. For instance, the work of Delhaize et al. (2017) was based on a jointly-selected infrared (from Herschel, with S/N ≥ 5 in at least one PACS or SPIRE band) and radio (VLA 3 GHz with S/N ≥ 5; Smolčić et al. 2017a) sample of SFGs in the COSMOS field, out to z ∼ 5. By performing double-censored survival analysis to account for sources undetected at either radio or FIR wavelengths, they found an evolving q_IR ∝ (1 + z)^{−0.19 ± 0.01}, which flattens to q_IR ∝ (1 + z)^{−0.12 ± 0.01} after removing 2σ outliers (as reported in Delvecchio et al. 2018), particularly radio-excess AGN. We repeat our IR and radio stacking analysis using the same sample of SFGs from Delhaize et al. (2017) (9575 sources), to demonstrate that our analysis leads to consistent results when matching the input sample.

We split the sample of Delhaize et al. (2017) among the same seven redshift bins analyzed in this work. For each, we perform median stacking of 3 GHz and IR images in all bands, combining both detections and non-detections. This approach should be comparable to the search for the median value carried out via survival analysis (Delhaize et al. 2017). Although we do not formally remove radio AGN in this check, we showed in Appendix D.1 that median radio stacking yields broadly consistent results (see Magnelli et al. 2015). Figure D.2 displays the median q_IR obtained by stacking the SFG sample of Delhaize et al. (2017) in different redshift bins (stars). This yields a best-fitting q_IR ∝ (1 + z)^{−0.11 ± 0.05}, that is fully consistent with the flatter trend of Delhaize et al. (2017) after removing 2σ outliers (triple dot-dashed line). This check proves our technique solid against different sample selections from the literature.

Fig. D.2. Median qIR as a function of redshift obtained by analysing the SFG sample of Delhaize et al. (2017, stars). Black lines indicate the median qIR − z trend of Delhaize et al. (2017) before (dot-dashed) and after (triple dot-dashed) removing 2σ outliers. The grey solid line marks the resulting best-fit qIR trend with redshift, that is highly consistent with that of Delhaize et al. (2017) after removing radio AGN. Numbers below each star denote the median M⋆ of the underlying sample.

Figure D.2 also highlights the important role played by the binning grid in driving a declining IRRC with redshift. In particular, the colour-coded M_⋆ clearly indicates how a joint IR and radio selection is sensitive to increasing galaxy M_⋆ with redshift. Moreover, the scatter of the IRRC reported by Delhaize et al. (2017) is around 0.35 dex, while the dispersion that we measured at M_⋆ > 10^10.5 M_⊙ (Sect. 4.2.2) is only 0.21−0.22 dex. This is similar to the value reported by Bell (2003) (i.e. 0.26 dex) for nearby galaxies, recently narrowed down to 0.16 dex in Molnár et al. (2021). A possible reason for the smaller than 0.35 dex dispersion in our study might be that we are splitting SFGs among different M_⋆, each carrying an intrinsically smaller dispersion compared to the full SFG sample. Because of the decreasing q_IR with M_⋆, binning only as a function of redshift leads to a mixture of different galaxy M_⋆ that results into a larger global dispersion.

All Tables

All Figures

Fig. 1. Distribution of the full COSMOS2015 (Laigle et al. 2016) source list over the COSMOS area (blue dots). The subset of 413 678 NUVrJ–based star-forming galaxies analysed in this work (red dots) includes sources from Laigle et al. (2016) and Muzzin et al. (2013) within the UltraVISTA area, with the exception of masked regions due to saturated or contaminated photometry. See Sect. 2 for details.

In the text

Fig. 2. Distribution of NUVrJ star-forming galaxies as a function of M⋆ and redshift. The colour scale in the central panel indicates the underlying sample size, increasing from blue to red. The grey-dashed grid encloses the 42 M⋆ − z bins into which we split our sample (413 678 objects). Galaxies within that grid are projected on the upper-left and bottom right histograms with redshift and M⋆, respectively (grey lines). The blue and red histograms represent the subsample with a signal-to-noise ratio (S/N) that is greater than 3 at 3 GHz and across all IR bands, respectively (see Sects. 2.2 and 2.3). The orange solid line marks the 90% M⋆ completeness limit of Laigle et al. (2016) for comparison.

In the text

	Fig. 3. Number of NUVrJ–based star-forming galaxies analysed in this work, as a function of M⋆ and redshift (black). For convenience, in each bin we report the fraction of sources with combined S/NIR > 3 over all FIR/sub-mm bands (red) and with S/N3 GHz > 3 (blue).
In the text

Fig. 4. Stacked cutouts of NUVrJ–based SFGs at 0.8 < z < 1.2, as a function of M⋆ (left to right, expressed in log M⊙). Within each bin, we stacked only those sources with S/N < 3 at a given band. SCUBA 850 μm and AzTEC 1100 μm images are smoothed with a Gaussian kernel to ease the visualisation. Each cutout size is 8 × FWHM of the PSF, while for Spitzer-MIPS, we chose 13 × FWHM. Below each cutout, we report the corresponding S/N.

In the text

	Fig. 5. Image decomposition of median stacks at 250, 350, and 500 μm, for a specific bin at 0.8 < z < 1.2 and 11 < log(M⋆/M⊙) < 12. From left to right: the stacked image is separated among a point source PSF, the clustering signal, and a residual background term, respectively. The colour scale is normalised to the maximum in each cutout for visual purposes. See Sect. 3.2 for details.
In the text

Fig. 6. Best-fit template obtained from a SED-fitting decomposition (green lines), as a function of M⋆ (left to right, expressed in log M⊙) and redshift (top to bottom). Red circles indicate the IR/sub-mm photometry (MIPS 24 μm, PACS 100−160 μm, SPIRE 250−350−500 μm, SCUBA 850 μm, and AzTEC 1.1 mm), while downward arrows mark the corresponding 3σ upper limits. The red dotted line is the best-fit AGN template, shown in the only bin where its significance is above 3σ. Green dashed lines represent SEDs without FIR measurements and at z ≳ 1.5, for which the integrated LIR is interpreted as 3σ upper limit (5/42 bins). MIPS 24 μm flux densities are not used in the fitting.

In the text

Fig. 7. SFR–M⋆ relation of the NUVrJ star-forming galaxies selected in this work, colour-coded by redshift over 0.1 < z < 4.5. At fixed M⋆ and redshift, SFRIR measurements are converted from the LIR obtained from IR/sub-mm SED-fitting (circles). Downward arrows indicate 3σ upper limits. For completeness, we also show the SFRIR + UV estimates by combining SFRIR with UV-uncorrected SFRs (open squares). Solid lines mark the evolving MS relation between SFRIR + UV and M⋆ at different redshifts (Schreiber et al. 2015). We observe an excellent agreement with the MS, even below M⋆ ∼ 109.5 M⊙ that relies upon a linear extrapolation from Schreiber et al. (2015) that is not constrained by previous data.

In the text

Fig. 8. AGN-uncorrected qIR evolution as a function of redshift (x-axis) and M⋆ (colour bar). Errors on qIR represent the 1σ scatter around the median value, estimated via bootstrapping over LIR and L1.4 GHz. For comparison, other IRRC trends with redshift are taken from the literature (black lines): Bell (2003, dotted); Magnelli et al. (2015, dashed); Delhaize et al. (2017, dot-dashed). Our qIR values still include the contribution of radio AGN. See Sect. 4.1 for details.

In the text

Fig. 9. Distribution of qIR as a function of redshift, split across increasing M⋆ bins. In each panel, we compare the qIR estimates of individual radio detections (black dots) with the median stacked values of non-detections (yellow squares) and with the weighted-average qIR of detections and non-detections together (from Eq. (3), blue stars). Green upward arrows indicate the corresponding threshold qIR, lim above which radio detections become inaccessible from an M⋆-selected sample. We select relatively complete M⋆ − z bins, in which at least 70% of radio detections have qIR below the corresponding threshold. This criterion identifies 13 bins (in red square brackets). Within these bins, the peak of qIR distribution (qIR, peak) is indicated with red open circles. In the two highest M⋆ bins, the best fitting trends with redshift are shown by the red dashed lines. Each panel reports the number of individual 3 GHz sources and their fraction with S/NIR > 3. See Sect. 4.2.1 for details.

In the text

Fig. 10. qIR distribution of 3 GHz detections in the two highest M⋆ bins (left and right panels), after correcting for the internal qIR–redshift dependence. Only sources within complete z-bins were considered at each M⋆. The total distribution (black histogram) was dissected among SF-dominated (blue) and AGN-dominated (red) radio sources, and fitted with two separate Gaussian functions (dot-dashed lines). While the SF population was fitted first, the AGN part was fitted in a second step from the total–SF residuals. The 1σ dispersion attributed to SF equals 0.20 and 0.23 dex at 1010.5 < M⋆/M⊙ < 1011 and 1011 < M⋆/M⊙ < 1012, respectively. Bottom panels: corresponding cumulative Gaussian fits, both normalised to unity. The vertical green dotted line marks the threshold that best separates between the SF and AGN populations (see Table 3), which rejects about 70% AGN and only 3−4% SFGs (open circles).

In the text

	Fig. 11. Same as Fig. 10, but normalising the peak to the flatter qIR − z trend calibrated after removing AGN (Sect. 4.2.3). This two-fold approach slightly improves the qIR decomposition, as highlighted by the larger cumulative fraction of radio AGN that are rejected below the qIR threshold (red open circles, 81% against the previous 70%).
In the text

Fig. 12. Distribution of qIR as a function of redshift and M⋆, after removing radio AGN (red dots). Symbols are the same as in Fig. 9, except for the median qIR estimates (blue stars), which are here re-calculated after removing radio-excess AGN (red dots). Fractions of radio-detected AGN and SFGs are reported in each M⋆-bin, as well as the fraction of AGN relative to the full M⋆ sample analysed in this work (in brackets). The blue and red solid lines denote the intrinsic IRRC of SFGs and the locus below which we classify radio sources as AGN (0.43 dex below the IRRC), respectively.

In the text

Fig. 13. Intrinsic (i.e. AGN-corrected) qIR evolution as a function of redshift (x-axis) and M⋆ (colour bar). The LIR estimates are the same reported in Fig. 8, while L1.4 GHz measurements have been re-calculated after excluding radio-detected AGN (Sect. 4.2). For comparison, other IRRC trends with redshift are taken from the literature (black lines): Bell (2003, dotted); Magnelli et al. (2015, dashed); Delhaize et al. (2017, dot-dashed) and their AGN-corrected version after removing 2σ outliers (triple dot-dashed lines).

In the text

	Fig. 14. Distribution of AGN-corrected qIR as a function M⋆, colour-coded by redshift (stars). At each M⋆, open squares indicate the median qIR values at z = 1, obtained after propagating the uncertainties of slope and normalisation of the corresponding qIR − z fit and interpolating each at z = 1. These values were fitted with a linear function in log − log space (black dashed line).
In the text

Fig. 15. Distribution of 3 GHz effective radius (in kpc) as a function of redshift and colour-coded by M⋆ (stars). Size measurements are taken from median stacked 3 GHz images of non-detections. Upper limits are given for unresolved stacks and correspond to the angular 3 GHz beam size (0.75″, grey dashed line). We observe a clear increase of Re with galaxy M⋆. The bins at M⋆ > 1010.5 M⊙ with resolved emission are fitted with a power-law redshift trend, which yields Re ∝ (1 + z)−0.18 ± 0.07 (orange solid line). A comparison study by Jiménez-Andrade et al. (2019) is shown (blue solid line) for 3 GHz detected SFGs at similar M⋆ in COSMOS, obtaining Re ∝ (1 + z)−0.26 ± 0.08.

In the text

Fig. 16. Evolution of qIR as a function of SFR surface density (ΣSFR = SFRIR + UV/2$\pi R_{\text{e}}^2$), colour coded by redshift (left panel) and by M⋆ (right panel). The SFRIR + UV estimates are taken from Sect. 3.3, while the effective radius Re is measured from stacked 3 GHz images via 2D Gaussian fitting. This plot shows a significant anti-correlation similar in slope to that observed with M⋆ (Sect. 4.3), marked by the black dashed line (qIR ∝ ( − 0.13 ± 0.02) × logΣSFR). For comparison, the best-fit trend with rest-frame optical (5000 Å) sizes estimated from van der Wel et al. (2014) scaling relation is shown (grey dotted line).

In the text

Fig. 17. Comparison of qIR between our full SFG sample (MS+SB, stars) and the SB subsample (circles), as a function of redshift. To mitigate the incompleteness of an IR-based selection of SBs, we only show bins with M⋆ > 1010.5 M⊙. Black lines highlight best-fit IRRC trends from the literature for comparison. This test suggests that qIR evolves irrespective of whether a galaxy is on or above the MS.

In the text

	Fig. A.1. Comparison of total flux densities of 3 GHz detections between our procedure and catalogue flux densities, both at S/N > 5 (Smolčić et al. 2017a, red dots) and at 3 < S/N < 5 (Jin et al. 2018, blue dots). Squares highlight the median ratio at various intervals. The global offset and dispersion suggest a good agreement within the uncertainties down to S/N ∼ 3.
In the text

	Fig. A.2. Top panel: comparison between median L1.4 GHz (x-axis) and rms-weighted mean L1.4 GHz (y-axis) for combined 3 GHz detections and non-detections, following Eq. (3). Colours indicate various M⋆ bins. Bottom panel: same comparison, but referring to 3 GHz undetected sources only.
In the text

	Fig. B.1. Stacked cutouts of our sample at 0.8 < z < 1.2, as a function of M⋆ (left to right, expressed in log M⊙). Only individual undetected sources (S/N < 3) are stacked. Top, middle and bottom rows: VLA 3 GHz, VLA 1.4 GHz and MIGHTEE 1.3 GHz data, respectively. Each cutout size corresponds to 8 × FWHM of the beam. Below each cutout we report the corresponding S/N of the total stacked flux.
In the text

Fig. B.2. Comparison between rest-frame 1.4 GHz spectral luminosity L1.4 GHz obtained from 3 GHz stacks (x-axis) and ancillary radio stacks (y-axis) using VLA (1.4 GHz, circles) and MIGHTEE (1.3 GHz, squares) data. We assumed a single spectral index α = −0.75 to scale flux densities from 3 GHz to 1.4 GHz. Colours indicate different M⋆ ranges. Downward arrows with open symbols mark 3σ upper limits if S/N < 3. The broad agreement between the various datasets suggests that using a single α = −0.75 is a reasonable assumption across the full M⋆ range explored in this work.

In the text

Fig. C.1. Same as Fig. 10, but fitting the total qIR distribution of 3 GHz detections (black histogram) simultaneously with a SF Gaussian (blue) and an AGN Gaussian (red dashed). The 1σ dispersion has remained unchanged to about 0.21 dex. Bottom panels: corresponding cumulative Gaussian fits, both normalised to unity. The vertical green dotted line marks the 2σ threshold of 0.42 dex below which we consider a source as radio-excess AGN. As opposed to Fig. 10, this cutoff removes only 30−50% of the total radio AGN population. However, we estimate this effect to be more prevalent in lower M⋆ galaxies. This implies that accounting for such mis-classified radio AGN would likely strengthen our final M⋆-dependent qIR.

In the text

Fig. D.1. Comparison between AGN-corrected L1.4 GHz from this work (x-axis) and median L1.4 GHz obtained from stacking detections and non-detections together (Magnelli et al. 2015, y-axis). Different M⋆ ranges are colour-coded, while downward arrows mark 3σ upper limits for 2/37 bins. Despite these different approaches, we notice a very good agreement in all bins, that strengthens the reliability of our recursive AGN subtraction.

In the text

Fig. D.2. Median qIR as a function of redshift obtained by analysing the SFG sample of Delhaize et al. (2017, stars). Black lines indicate the median qIR − z trend of Delhaize et al. (2017) before (dot-dashed) and after (triple dot-dashed) removing 2σ outliers. The grey solid line marks the resulting best-fit qIR trend with redshift, that is highly consistent with that of Delhaize et al. (2017) after removing radio AGN. Numbers below each star denote the median M⋆ of the underlying sample.

In the text