The infrared-radio correlation of star-forming galaxies is stronglyM⋆-dependent but nearly redshift-invariant sincez∼ 4

Over the past decade, several works have used the ratio between total (rest 8−1000μm) infrared and radio (rest 1.4 GHz) luminosity in star-forming galaxies (qIR), often referred to as the infrared-radio correlation (IRRC), to calibrate the radio emission as a star formation rate (SFR) indicator. Previous studies constrained the evolution ofqIRwith redshift, finding a mild but significant decline that is yet to be understood. Here, for the first time, we calibrateqIRas a function ofbothstellar mass (M⋆) and redshift, starting from anM⋆-selected sample of > 400 000 star-forming galaxies in the COSMOS field, identified via (NUV − r)/(r − J) colours, at redshifts of 0.1 < z < 4.5. Within each (M⋆,z) bin, we stacked the deepest available infrared/sub-mm and radio images. We fit the stacked IR spectral energy distributions with typical star-forming galaxy and IR-AGN templates. We then carefully removed the radio AGN candidates via a recursive approach. We find that the IRRC evolves primarily withM⋆, with more massive galaxies displaying a systematically lowerqIR. A secondary, weaker dependence on redshift is also observed. The best-fit analytical expression is the following:qIR(M⋆, z) = (2.646 ± 0.024) × (1 + z)( − 0.023 ± 0.008)–(0.148 ± 0.013) × (log M⋆/M⊙ − 10). Adding the UV dust-uncorrected contribution to the IR as a proxy for the total SFR would further steepen theqIRdependence onM⋆. We interpret the apparent redshift decline reported in previous works as due to low-M⋆galaxies being progressively under-represented at high redshift, as a consequence of binning only in redshift and using either infrared or radio-detected samples. The lower IR/radio ratios seen in more massive galaxies are well described by their higher observed SFR surface densities. Our findings highlight the fact that using radio-synchrotron emission as a proxy for SFR requires novelM⋆-dependent recipes that will enable us to convert detections from future ultra-deep radio surveys into accurate SFR measurements down to low-M⋆galaxies with low SFR.


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): 3.75 × 10 12 [Hz] − log(L 1.4 GHz [W Hz −1 where 3.75 × 10 12 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. 1993Condon et al. , 2002Murphy 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. 2011Murphy et al. , 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). A&A 647, A123 (2021) For example, some studies of local star-forming galaxies (SFGs), ranging from dwarf (e.g., Wu et al. 2008) to ultraluminous infrared galaxies (ULIRGs; L IR > 10 12 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 IR ∝ L 0.75−0.90 1.4 GHz (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 10 M ) SFGs because of their increasing compactness (i.e. the size-mass relation R e ∝ M 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;) 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., 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 Ultra-VISTA 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 0 = 70 km s −1 Mpc −1 (Spergel et al. 2003).

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 2 ) boasts an exquisite photometric data set, spanning from the X-rays to the radio domain 1 . 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 Y JHK s image (blue dots in Fig. 1). In particular, this catalogue joins optical photometry from Subaru Hyper-Suprime Cam (2 deg 2 ; Capak et al. 2007) and from the Canada-France-Hawaii Telescope Legacy Survey (CFHT-LS, central 1 deg 2 ; 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 2 , of which 0.6 deg 2 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 2 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.
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.

Selecting star-forming galaxies via (NUV-r)/(r-J) colours
We aim to study the infrared-radio correlation within an Mselected 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 6 −10 8 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 −10 yr −1 ) and fully passive systems (sSFR < 10 −11 yr −1 ). We further selected galaxies with redshift 0.1 < z < 4.5 and 10 8 < M /M < 10 12 . This leaves us with a final sample of 413 678 star-forming galaxies (red dots in Fig. 1 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 8 < M /M < 10 9 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 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.

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.
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 ) 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 Mcomplete 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/submm 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 Gaussianlike 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.

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 −1 over 2.6 deg 2 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 2 of COSMOS, the early science release of MIGHTEE reaches a thermal noise of 2.2 µJy beam −1 , 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 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 greydashed 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. 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.

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). Notes. We note that subsets do not add up to make the final sample. . 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 /N IR > 3 over all FIR/sub-mm bands (red) and with S /N 3 GHz > 3 (blue).
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 Mselected 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).

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 2 (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 2 , 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 −1 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/submm 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.

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 S/N=4.0 S/N=22.4 8 < M* < 9 9 < M* < 9.5 9.5 < M* < 10 10 < M* < 10.5 10.5 < M* < 11 11 < M* < 12 been discussed in the literature (e.g., Bavouzet et al. 2008;Béthermin et al. 2010Béthermin et al. , 2012Béthermin et al. , 2015Kurczynski & 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 Heinis et al. 2013Heinis et al. , 2014Welikala et al. 2016). Notes. Uncertainties indicate the 1σ dispersion among all S /N > 3 stacks at a given band.
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.
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: 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.

Conversion to L IR and SFR
This section illustrates how we fit the observed FIR/submm SEDs to determine the total (8−1000 µm rest-frame) IR A123, page 6 of 29  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 χ 2 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 aposteriori the 24 µm data-point. Globally, our stacking analysis yields robust L IR estimates in 37/42 bins.
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 zdependent (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 UVuncorrected 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).
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 −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. 2012Whitaker et al. , 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). 3.5< z <4.5 This test compellingly demonstrates that our L IR can be deemed robust over the full M and redshift interval explored in this work.

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   (Schreiber et al. 2015). We observe an excellent agreement with the MS, even below M ∼ 10 9.5 M that relies upon a linear extrapolation from Schreiber et al. (2015) that is not constrained by previous data.
densities (<100 µJy), most of radio emission is thought to arise from star formation Padovani et al. 2015;Novak et al. 2017;Smolčić et al. 2017a), though some AGNrelated 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 −1 ) 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 3 . 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 A gauss /A 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: 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. 2009Ibar et al. , 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.

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 8 < M /M < 10 9 and 1.8 < z < 4.5; two at 10 9 < 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 4 ) across the SFR-M plane at M > 10 10 M . From stacking IR and radio images, they parametrised the evolution with redshift of the FIRC as: 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: 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).
When compared to the above literature, it is evident that our q IR values lie systematically below other studies at M > 10 11 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  The far-infrared luminosity used to compute q FIR was integrated between 42 and 122 µm rest-frame. This is quantified to be 1.91× smaller than the total L IR (Magnelli et al. 2015). q IR . However, it is worth showing it to quantify how much q IR changes after removing the radio AGN.

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.

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 deblended 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 nondetections (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 8 < M /M < 10 9 to 78.3% at 10 11 < M /M < 10 12 , 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.
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 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.

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 explor-  Fig. 10. q IR distribution of 3 GHz detections in the two highest M bins (left and right panels), after correcting for the internal q IR -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 10 10.5 < M /M < 10 11 and 10 11 < M /M < 10 12 , 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). Table 3. Comparison between fractions of radio-SFGs and radio-excess AGN below some threshold q thres , for the two highest M bins.
ing 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 A123, page 12 of 29 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 righthand 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 11 and 10 11 < M /M < 10 12 , 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 SFdominated 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 crosscontamination 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 3). This two-fold approach slightly improves the q IR decomposition, as highlighted by the larger cumulative fraction of radio AGN that are rejected below the q IR threshold (red open circles, 81% against the previous 70%).
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.

Re-calibrating the radio AGN threshold
According to the threshold defined above, we removed radioexcess 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.
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  Fig. 9, except for the median q IR 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. 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σ.
While at 10 11 < M /M < 10 12 all z-bins (0.1 < z < 4.5) were used to constrain this trend, at 10 10.5 < M /M < 10 11 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 SFdominated 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 radioexcess 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 radioexcess contrast (at fixed L 1.4 GHz ) is less evident. Therefore, we Table 4. Summary of the numbers and fractions of radio AGN and SFGs in different M bins after fitting the AGN-corrected q IR with redshift (Sect. 4.2.4 and Fig. 12).

M (M ) bin
q IR −z fit # radio AGN # radio SFGs # M sample (normalisation) (% radio-det, % M sample) (% radio-det, % M sample) (1) (2) (3) (4) Notes. Columns are sorted as follows: (1)  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 11 , radio AGN are 7.1% of all radio-detections and 2.2% of the full M sample of SFGs. At 10 11 < M /M < 10 12 , 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. 2015White et al. , 2017Wong 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. 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. , 2864. Since the brightness temperature reached by VLBA observations at about 0.01 resolution exceeds 10 6 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 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.

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 radiodetected 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 COS-MOS 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.

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 powerlaw 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.
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. . 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 √ χ 2 red in each M bin, to bring  the reduced χ 2 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 A123, page 16 of 29 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 χ 2 red = 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.
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: where A = (1 + z) and B = (log M /M − 10). The corresponding χ 2 red = 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 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.

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 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).

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 9 < M /M < 10 12 , 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 IR ∝ L 0.90 1.4 GHz . 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 IR ∝ L 0.75−0.90 1.4 GHz , 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 nonlinearity 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.,  by studying the relation between q IR and SFR surface density.

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 .

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 , 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 −1 kpc −2 ) to 'starbursts' (SBs, Σ SFR 2−4 M yr −1 kpc −2 ). In particular, this model argues that q IR slightly declines with Σ SFR (up to Σ SFR ∼ 1 M yr −1 kpc −2 ) 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 −1 kpc −2 , 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 / √ A r , where A r is the axial ratio. 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).
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 9 to 10 11 M (Fig. 7), while R 2 e 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  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  model a curved radio spectrum. Second, we label as 'SBs' those galaxies that lie >4× above the MS (Sect. 5.2), while  identify them as Σ SFR 2−4 M yr −1 kpc −2 . 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 .
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 −1 kpc −2 . This behaviour is not clearly seen in our data, which instead display a smoothly declining q IR −Σ SFR trend and are nearly redshiftinvariant. 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 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 −1 kpc −2 , 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 , 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., 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.

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 MSotherwise, 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 radiodetected 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.
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,  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 , 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.

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 radiodetected 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/submm 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 counterbalances 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 AGNdominated 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 A123, page 20 of 29 (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. 2011Reines et al. , 2014Reines & 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 radiodetected 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 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.

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: 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: 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 sublinear trend L IR ∝ L 0.90 1.4 GHz 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 IR+UV ∝ L 0.81 1.4 GHz . Such a radio deficit in the dwarf-galaxy regime could be possibly linked to shorter CRe scale heights (e.g., Helou & Bicay 1993; 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 Mdependent, 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.

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 . 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 IR ∝ L 0.90 1.4 GHz , 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., 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 , a 'conspiracy' of different factors might induce a q IR flattening at Σ SFR 1 M yr −1 kpc −2 . 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. 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(Hales et al. , 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. 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  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.
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.
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  map. It covers 1.7 deg 2 with an angular resolution of 2.5 , reaching rms = 12 µJy beam −1 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 −1 at 8.4 × 6.8 resolution over 1 deg 2 in the MIGHTEE early science data, but the effective depth is limited by confusion (∼5.5 µJy beam −1 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 2 ) 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 ): S tot /S peak > 0.35 −11/(S /N 1.45 peak ) . 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.  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 11 < M /M < 10 12 . 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 11 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.

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, colourcoded 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 selfconsistent manner. Magnelli et al. (2015) 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  Fig. 10, but fitting the total q IR 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 q IR .
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 11 and 90% at at 10 11 < M /M < 10 12 .
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 radiodetected 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 11 < M /M < 10 12 the radio AGN distribution is clearly broader, but far more offset than at 10 10.5 < M /M < 10 11 , 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 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  (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. the contribution of radio-detected SFGs to our weighted average L 1.4 GHz , since they make a substantial fraction of the Mselected sample in that M bin (∼45% , Table 4). Therefore, this test proves our radio AGN subtraction broadly consistent with a totally independent approach.

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 Mselected 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 COS-MOS 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 radioexcess 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 ). 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 bestfitting 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. 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.