Open Access
Issue
A&A
Volume 647, March 2021
Article Number A123
Number of page(s) 29
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202039647
Published online 18 March 2021

© I. Delvecchio et al. 2021

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

For nearly fifty years, astronomers have studied the observed correlation between total infrared (IR; rest-frame 8−1000 μm, i.e. LIR) and radio (e.g., rest-frame 1.4 GHz, i.e. L1.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 qIR, defined as (e.g., Helou et al. 1985; Yun et al. 2001):

(1)

where 3.75 × 1012 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 qIR) appears to hold over at least three orders of magnitude in both LIR and L1.4 GHz (e.g., Helou et al. 1985; Condon 1992; Yun et al. 2001). Broadly speaking, this is because the infrared emission comes from dust heated by fairly massive (≳5 M) OB stars, while radio emission arises from relativistic cosmic ray electrons (CRe) accelerated by shock waves produced when massive stars (≳8 M) explode as supernovae. Nevertheless, CRe are also subject to different cooling processes as they propagate throughout the galaxy, which are mainly caused by inverse Compton, bremsstrahlung, and ionisation losses (e.g., Murphy 2009).

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

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

For example, some studies of local star-forming galaxies (SFGs), ranging from dwarf (e.g., Wu et al. 2008) to ultra-luminous infrared galaxies (ULIRGs; LIR > 1012L; e.g., Yun et al. 2001) have concluded that the IRRC remains linear across a wide range of LIR. 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 (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 ≳ 1010M) SFGs because of their increasing compactness (i.e. the size-mass relation , van der Wel et al. 2014), which might enhance their ability to retain the gas ejected by stars. However, this is likely to break down towards lower M galaxies due to smaller sizes and lower obscuration (e.g., Bourne et al. 2012). On the other hand, non-calorimetric models or the optically thin scenario (Helou & Bicay 1993; Niklas & Beck 1997; Bell 2003; Lacki et al. 2010) support the argument that several physical mechanisms cancel each other out, creating a sort of ‘conspiracy’ that keeps the IRRC unexpectedly tight and linear. Indeed, both IR and radio luminosities are expected to underestimate the total SFR in low M and low-SFR surface density galaxies (Bell 2003), inducing a departure of the IRRC from linearity. This was not, however, observed in our study. Radio synchrotron models postulate that such small galaxies are not able to prevent CRe from escaping, causing a global deficit of radio emission at a fixed SFR. Similarly, the IR domain becomes less sensitive to SFR in low-M galaxies (e.g., Madau & Dickinson 2014), generating an IR deficit of a similar amount, which might counterbalance the radio emission and keep the IRRC linear. Understanding the discrepancy between model predictions and observations is crucial, since the linearity (or not) of the IRRC has direct implications for using radio emission as a SFR tracer.

From an observational perspective, it is widely recognised that there is a tight relation linking SFR and M in nearly all SFGs, namely, namely the ‘main sequence’ of star formation (MS, scatter ∼0.2−0.3 dex). This relation holds from z ∼ 5 down to the Local Universe (e.g., Brinchmann et al. 2004; Noeske et al. 2007; Elbaz et al. 2011; Whitaker et al. 2012; Speagle et al. 2014; Schreiber et al. 2015; Lee et al. 2015), showing a flattening at high M and an evolving normalisation with redshift. Because the SFR is directly linked to LIR, especially in massive SFGs (Kennicutt 1998), the existence of the MS gives an additional argument that studying qIR 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 qIR ∝ (1 + z)[ − 0.2 : −0.1], although the physical explanation for such an evolution is still uncertain. Somewhat different conclusions were reached by other works (e.g., Garrett 2002; Appleton et al. 2004; Ibar et al. 2008; Jarvis et al. 2010; Sargent et al. 2010; Bourne et al. 2011) that ascribe this apparent evolution to selection effects; for instance, comparing flux-limited samples, each with a different selection function.

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

The main goal of the present paper is to calibrate the IRRC, for the first time, as a function of both M and redshift over a broad range. To this end, we start from an M-selected sample of > 400 000 galaxies at 0.1 < z < 4.5 collected from deep UltraVISTA images in the Cosmic Evolution Survey (Scoville et al. 2007; centered at RA = +150.11916667; Dec = +2.20583333 (J2000)). Then we leverage the new de-blended far-IR/sub-mm data (Jin et al. 2018) recently compiled in COSMOS, which allow us to circumvent blending issues due to poor angular resolution and measure LIR 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 qIR 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 LIR (Sect. 3.1) and L1.4 GHz (Sect. 3.4). The average qIR 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 H0 = 70 km s−1 Mpc−1 (Spergel et al. 2003).

2. Multi-wavelength data and sample selection

In this section, we describe the creation of a Ks prior catalogue, which we used to select our parent sample in the COSMOS field. The COSMOS field (2 deg2) boasts an exquisite photometric data set, spanning from the X-rays to the radio domain1. 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 YJHKs image (blue dots in Fig. 1). In particular, this catalogue joins optical photometry from Subaru Hyper-Suprime Cam (2 deg2; Capak et al. 2007) and from the Canada-France-Hawaii Telescope Legacy Survey (CFHT-LS, central 1 deg2; McCracken et al. 2001); the near-infrared (NIR) bands Y, J, H, and Ks from UltraVISTA DR2 (down to Ks< 24.5 in the central 1.5 deg2, of which 0.6 deg2 are covered by ultra-deep stripes with limiting Ks < 25.2; McCracken et al. 2012), and from CFHT H and Ks observations obtained with the WIRCam (Ks < 23.9 outside the UltraVISTA area; McCracken et al. 2001). Over the full 2 deg2 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.

thumbnail Fig. 1.

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

Open with DEXTER

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 Ks–selected sources from the UltraVISTA catalogue of Muzzin et al. (2013) (3σ limit of Ks < 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 Ks 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 + zp]). 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 (zs quality flag > 3 ∧ |zs − zp|< 0.1 × (1 + zp)). Infrared/sub-mm flux densities were de-blended and re-extracted via the prior-based fitting algorithm presented in Jin et al. (2018), which we briefly describe in Sect. 2.2.

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

We aim to study the infrared-radio correlation within an M-selected sample of star-forming galaxies. To this end, we make use of the rest-frame, dust-corrected (NUV − r) and (r − J) colours available in the parent catalogues (hereafter NUVrJ). As opposed to the widely used UVJ criterion, the (NUV − r) colour is more sensitive to recent star formation (106 − 108 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 108 < M/M < 1012. This leaves us with a final sample of 413 678 star-forming galaxies (red dots in Fig. 1), out of which 22 38 (5.4%) are spectroscopically confirmed. The fraction of catastrophic failures (|zs − zp|> 0.15 × (1 + zs)) is only 3.4%.

Such a sizable sample enables us to bin galaxies as a function of both M and redshift, while maintaining good statistical power. Figure 2 shows our sample in the M–redshift diagram, highlighting the chosen grid. We note that the M uncertainties taken from the parent catalogues incorporate the covariant errors on stellar population ages and dust reddening. These average M uncertainties are 0.2 dex at 108 < M/M < 109 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 1010M 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 LIR and those extrapolated from the MS relation (Schreiber et al. 2015) also at M < 109.5M, 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 > 109.5M, where our sample is largely complete. Moreover, in light of our main result, namely, that qIR 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.

thumbnail Fig. 2.

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

Open with DEXTER

2.2. Infrared and sub-mm de-blended data

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

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

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

Table 1.

Main numbers of priors and detections that characterise our final sample of 413 678 star-forming galaxies.

2.3. Radio data in the COSMOS field

For our analysis, we exploited data from the VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017a), one of the largest and deepest radio survey ever conducted over a medium sky area like COSMOS. With 384 h of observations, the final mosaic reaches a median root mean square (rms) of 2.3 μJy beam−1 over 2.6 deg2 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 deg2 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 Ks + MIPS 24 μm + VLA priors applied on FIR/sub-mm images (Sect. 2.2). Similarly, VLA radio flux densities were re-extracted based on the same PSF fitting technique down to S/N > 3. However, for the sake of consistency with publicly available VLA catalogues, we take the 1.4 and 3 GHz radio flux densities of S/N > 5 sources from Schinnerer et al. (2010) and Smolčić et al. (2017a), respectively. As a sanity check, we show in Appendix A that our procedure leads to fully consistent total 3 GHz flux densities.

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

3. Stacking analysis

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

In addition, our notably large star-forming galaxy sample allows us to bin as a function of both M and redshift, as shown in Fig. 3. For each bin, we also report the total number of M-selected SFGs (black), as well as the corresponding fractions having combined S/NIR > 3 (red) and S/N3 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).

thumbnail Fig. 3.

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

Open with DEXTER

3.1. Infrared and sub-mm stacking

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

To produce stacked and rms images in each band, we used the publicly available IAS stacking library2 (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 deg2, 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/sub-mm bands at 0.8 < z < 1.2 (i.e. close to the median redshift of our sample) as a function of M. As expected from the tight MS relation that links M and SFR in star-forming galaxies, stacks at low M display lower S/N, despite having larger numbers of input sources.

thumbnail Fig. 4.

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

Open with DEXTER

3.2. Correcting for clustering bias

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

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

(2)

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

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

Table 2.

Average fraction of clustering signal at each FIR/sub-mm band.

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.

thumbnail Fig. 5.

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

Open with DEXTER

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

(3)

Flux uncertainties were propagated in quadrature. For stacks with S/N < 3 in which we could not constrain the clustering correction, Sstack was set to the noise level of the stacked map (i.e. equal to its uncertainty). This way the weighted-average flux Sbin 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 Sbin has a S/N < 3, then it is set to three times the noise and interpreted as a 3σ upper limit.

3.3. Conversion to LIR and SFR

This section illustrates how we fit the observed FIR/sub-mm SEDs to determine the total (8−1000 μm rest-frame) IR luminosity within each M − z bin. To this end, we use the two-component SED-fitting code developed by Jin et al. (2018) (see also Liu et al. 2018). Briefly, this includes: three mid-infrared AGN torus templates from Mullaney et al. (2011); 15 dust continuum emission models by Magdis et al. (2012) that were extracted from Draine & Li (2007) to best reproduce the average SEDs of MS (14) or SB (1) galaxies at various redshifts. While the Draine & Li (2007) models were based on a number of physical parameters, the library of Magdis et al. (2012) depends exclusively on the mean radiation field ⟨U⟩ = LIR per unit dust mass (Md) 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 Md. 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, LIR, ⟨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 LIR 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 LIR is interpreted as 3σ upper limit (5/42 bins). Even though 24 μm has long been used as a proxy for LIR, 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 LIR obtained from SEDs without FIR data, but only at z ≲ 1.5. That is the case for a few bins at the lowest M, in which the SED reproduces a-posteriori the 24 μm data-point. Globally, our stacking analysis yields robust LIR estimates in 37/42 bins.

thumbnail Fig. 6.

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

Open with DEXTER

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

Given the tight correlation between LIR 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. SFRIR) and UV-uncorrected SFRs, in the deepest CANDELS fields, to calibrate the star-forming MS over an unprecedented M (down to M = 109.5M) and redshift range (z ≲ 4). Since we carried out a similar analysis, it is worth checking whether the SFR estimates based on our IR stacking reproduce or not the MS of Schreiber et al. (2015).

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

Within each M − z bin, we simply take the median SFRUV, and we add it to the SFRIR corresponding to the stacked LIR, calculated as SFRIR = 10−10LIR/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 SFRIR (circles) and the total SFRIR + UV (open squares) for comparison. Downward arrows are 3σ upper limits scaled from LIR. 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 SFRUV appears generally negligible compared to the total SFR, it becomes as high as SFRIR towards low M and low redshift (e.g., Whitaker et al. 2012, 2017). Our median values agree with Schreiber et al. (2015), even below M ∼ 109.5M, at which we extrapolate the MS relation due to lack of previous data (dashed lines). This test compellingly demonstrates that our LIR can be deemed robust over the full M and redshift interval explored in this work.

thumbnail Fig. 7.

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

Open with DEXTER

3.4. Radio stacking at 3 GHz

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

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

Within the UltraVISTA area analysed in this work, the 3 GHz rms (2.3 μJy beam−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 × FWHM3 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 (Speak) 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 (Stot) are calculated by fitting a 2D elliptical Gaussian function to the median stacked image, using the IDL routine MPFIT2DFUN3. 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 Agauss. The integrated flux error was computed by multiplying the peak flux error by , where Abeam 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:

(4)

where S/Npeak 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/Npeak < 3; these are all among the 5 bins without LIR estimates from IR stacking (Sect. 3.3). Analogously to our treatment of the IR measurements, we discard all those 5 bins from the rest of our analysis.

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

Finally, we combined the radio stacked flux densities within each M − z bin together with individual detections, following Eq. (3). The combined 3 GHz flux densities were first scaled to 1.4 GHz assuming of Sν ∝ να, with spectral index α = −0.75 ± 0.1 (e.g., Condon 1992; Ibar et al. 2009, 2010). This assumption is discussed in Appendix B. Lastly, 1.4 GHz flux densities were converted to rest-frame 1.4 GHz spectral luminosities (L1.4 GHz), again assuming α = −0.75. Formal L1.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 L1.4 GHz robust across the full range of M and redshift analyzed in this work. We note that our L1.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 L1.4 GHz. This issue is addressed in Sect. 4.2.

4. The IRRC and the contribution of radio AGN

Using the median LIR and L1.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 qIR errors. Among the 42 M − z bins analyzed in this work, 37 yield robust estimates of LIR and L1.4 GHz, while the remainder are discarded from the following analysis. Unsurprisingly, these latter 5 bins (three at 108 < M/M < 109 and 1.8 < z < 4.5; two at 109 < M/M < 109.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. qIR before removing radio AGN

Figure 8 shows the average qIR 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 qIR = 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 qFIR4) across the SFR–M plane at M > 1010M. From stacking IR and radio images, they parametrised the evolution with redshift of the FIRC as: qFIR = (2.35 ± 0.08) × (1 + z)−0.12 ± 0.04, where the normalisation is scaled to 2.63 in the qIR 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: qIR = (2.88 ± 0.03) × (1 + z)−0.19 ± 0.01. While this trend appears somewhat steeper than that of Magnelli et al. (2015), we note that Delhaize et al. (2017) did not formally remove objects with significant radio excess, while Magnelli et al. (2015) performed median radio stacking to mitigate the impact of potential outliers such as radio AGN. Nevertheless, Delhaize et al. (2017) argue that the IRRC trend with redshift would flatten if applying a 3σ-clipping: qIR = (2.83 ± 0.02) × (1 + z)−0.15 ± 0.01, which becomes fully consistent with the findings of Magnelli et al. (2015).

thumbnail Fig. 8.

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

Open with DEXTER

When compared to the above literature, it is evident that our qIR values lie systematically below other studies at M > 1011M, while lower M galaxies lie closer or slightly above them. In other words, our qIR estimates seem to display a clear M stratification, with the most massive galaxies having typically lower qIR 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 L1.4 GHz, particularly at high M (see e.g., Best & Heckman 2012) where radio AGN feedback is known to be prevalent. Our LIR 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 qIR. Selecting typical SFGs on the MS is expected to, however, reduce the incidence of powerful radio AGN expected in massive hosts since most radio AGN at z < 1 are found to reside in quiescent galaxies (e.g., Hickox et al. 2009; Goulding et al. 2014). It is for these reasons that we caution that Fig. 8 should be taken as the AGN-uncorrected qIR. However, it is worth showing it to quantify how much qIR changes after removing the radio AGN.

4.2. Searching for radio AGN candidates

In this section, we describe how we carried out a detailed study aimed at identifying potential radio AGN, removing them, and ultimately deriving the intrinsic qIR 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 L1.4 GHz of non-detections (Fig. A.2, bottom panel). Indeed, if the contribution of radio-undetected AGN were substantial, the corresponding mean L1.4 GHz would be significantly higher than the median L1.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 qIR distribution traced by individual 3 GHz detections as a function of M and redshift. First, we identify a subset of radio detections at M > 1010.5M that is representative of an M-selected sample. Then we decompose their qIR 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 > 1010.5M (Sect. 4.2.3). Then we extrapolate this calibration towards lower M bins (Sect. 4.2.4), where a similar in-depth analysis was not possible due to radio-detections being strongly incomplete in this M regime. Finally, the intrinsic (i.e. AGN-corrected) IRRC as a function of M and redshift is presented in Sect. 4.3.

4.2.1. The qIR distribution of radio detections

In order to study the qIR distribution of 3 GHz detections, we need to calculate their average LIR 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/NIR > 3, therefore reliable LIR measurements from SED-fitting of FIR/sub-mm de-blended photometry (Jin et al. 2018). For the remainder of the sample, we stack again their IR/sub-mm images in all bands in each M − z bin. Stacked IR flux densities are corrected for clustering bias and converted to LIR following the same procedure adopted for the prior M sample (Sect. 3.1). Median stacked LIR 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 LIR 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 LIR are always systematically below the 3σLIR 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 LIR of individual non-detections.

From this analysis, we are well-placed to explore the full qIR distribution of 3 GHz detections at different M and redshifts. Figure 9 shows qIR as a function of redshift, split in six M bins. Black dots mark individual 3 GHz detections, blue stars represent the qIR obtained by combining detections and non-detections (same as in Fig. 8), while yellow squares are the stacks of non-detections only. In each panel we report the number of 3 GHz detected sources and the fraction of them with combined S/NIR > 3. This fraction strongly increases with M, from 3.5% at 108 < M/M < 109 to 78.3% at 1011 < M/M < 1012, which implies that at the lowest M nearly all qIR 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.

thumbnail Fig. 9.

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

Open with DEXTER

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 qIR of single radio detections against a specific threshold (qIR, 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 LIR 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, qIR, lim indicates the limiting qIR at which a typical MS galaxy of a given M, z and LIR drops below the 3 GHz detection limit, which translates into a lower qIR limit. In other words, sources with qIR > qIR, 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 qIR, lim. This cutoff enables us to narrow down the position of the mode of the qIR 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 L1.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 LIR 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, LIR increases with redshift similarly to L1.4 GHz, giving rise to a nearly flat qIR locus. At lower M, instead, LIR stands typically below the IR detection limit, thus not bound to a monotonic redshift increase. This effect causes an apparent decrease of qIR 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 1010 < M/M < 1010.5 is insufficient for us to constrain a redshift trend, we only consider the remaining 12 complete bins placed at M > 1010.5. For each of them, we identify the peak of the corresponding qIR distribution of radio detections, namely, qIR, peak (see red open circles in Fig. 9). We note that qIR, 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 qIR, peak with redshift using the IDL routine MPFIT2DFUN, obtaining the best-fit expressions shown in Fig. 10.

thumbnail Fig. 10.

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

Open with DEXTER

4.2.2. Identifying radio AGN at high M

After fitting the qIR, peak trend with redshift for the two M bins, we need to account for such redshift dependence while exploring the qIR distributions of radio detections. To this end, we align the position of qIR, 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 qIR distribution is displayed in Fig. 10 for the two highest M bins (left and right panel). Each total histogram (black) includes the contribution from both galaxy- and AGN-dominated radio sources. We proceed to dissecting into the two components as follows, leaving the discussion on how radio AGN affect the redshift trend in Sect. 4.2.3.

Assuming that the peak of the distribution is populated by radio-detected SFGs, and that the intrinsic qIR distribution of SFGs is symmetric around the peak, we mirror the right-hand side of the observed qIR distribution to the left side. This symmetric function is interpreted as the intrinsic qIR 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 1010.5 < M/M < 1011 and 1011 < M/M < 1012, 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 qIR 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 qIR, 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 LIR, 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 > 1010.5M are also individually detected at IR wavelengths (see Fig. 3). Therefore, taking a single stacked LIR in each bin does not strongly impact the calibration of the SF locus.

Choosing the best dividing line between AGN and SF-dominated radio sources is a challenging, and somewhat arbitrary task. Moving the threshold to higher qIR increases the purity of SFGs to the detriment of completeness, and vice-versa for a lower threshold. Here we attempt to reach low levels of cross-contamination between the SF and AGN populations, while keeping a high completeness of the SF population. For this reason, we checked the cumulative qIR 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 (qthres) were examined: (1) qthres = qpeak–1σ; (2) qthres = qpeak–2σ; (3) qthres = qpeak–3σ; (4) qthres = qcross, AGN = SF. In this formalism, qpeak is still the peak of the SF population (blue Gaussian fit in Fig. 10), and σ its dispersion, while qcross, 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.

Table 3.

Comparison between fractions of radio-SFGs and radio-excess AGN below some threshold qthres, for the two highest M bins.

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

4.2.3. Re-calibrating the radio AGN threshold

According to the threshold defined above, we removed radio-excess AGN from our 12 complete z-bins at M > 1010.5M. Then we combined the remaining radio-detected SFGs with stacks of non-detections to compute the new L1.4 GHz in those bins, which should be free from AGN contamination. We verified that the new L1.4 GHz shifts the previously determined qIR (blue stars in Fig. 9) upward by a certain amount. In those complete bins, we fitted the AGN-corrected qIR 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.

thumbnail Fig. 11.

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

Open with DEXTER

To test the robustness of the newly derived qIR − z trend, we again shift the qIR measurements of individual detections by the offset from such a trend at each z-bin, and perform a second qIR 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 ΔqAGN = 0.43 dex). Moreover, the cumulative histograms (bottom panels) underline that this latter decomposition rejects about 81% of radio AGN below the threshold, as opposed to 70% estimated in the first step (see red open circles in Figs. 10 and 11), while missing a comparable 3−4% of SFGs. This confirms the effective improvement led by our re-calibration of the SF locus in removing radio AGN.

As shown in the updated Fig. 12 (at M > 1010.5M), subtracting radio AGN (red dots) based on this latter locus shifts all the median qIR (blue stars) exactly on the fitted qIR − 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 qIR 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σ.

thumbnail Fig. 12.

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

Open with DEXTER

While at 1011 < M/M < 1012 all z-bins (0.1 < z < 4.5) were used to constrain this trend, at 1010.5 < M/M < 1011 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 qIR estimates.

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

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

The fact that in both M bins the subtraction of radio-excess AGN leads to a flattening of the qIR − 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 qIR 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 qIR evolutionary trend could be explained by residual radio AGN activity within the SF locus. We estimate the overall fraction of ‘pure’ SFGs to be 95% at z < 1.2 and 90% at z > 1.2. Such minimal AGN contamination is probably more important at higher redshifts because SFGs are intrinsically IR brighter, so the radio-excess contrast (at fixed L1.4 GHz) is less evident. Therefore, we argue that any further correction for misclassified radio AGN would induce an even flatter qIR trend with redshift.

Finally, our approach leads to the following fractions of radio-excess AGN. At 1010.5 < M/M < 1011, radio AGN are 7.1% of all radio-detections and 2.2% of the full M sample of SFGs. At 1011 < M/M < 1012, 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).

Table 4.

Summary of the numbers and fractions of radio AGN and SFGs in different M bins after fitting the AGN-corrected qIR with redshift (Sect. 4.2.4 and Fig. 12).

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

4.2.4. Extrapolating the SF-vs.-AGN loci at low M

We extrapolate the qIR − 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 < 1010.5M are not representative of an M-selected sample. In particular, a galaxy of a given M and redshift, with infrared luminosity LIR 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 qIR offset between median measurements (blue stars) and individual radio detections (black dots). The latter lie systematically below the median qIR, 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 qIR 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 qIR is extrapolated from that calibrated at higher M, namely qIR ∝ (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 qIR. In other words, at M < 1010.5M, we assume a constant qIR − 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 qIR 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 qIR 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 qIR and search again for the best normalisation that fits the new AGN-corrected qIR measurements with redshift. We repeat this procedure twice, that is, until all median qIR are unchanged within the uncertainties, at each M. This condition sets the end of our recursion.

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

In most bins at M < 1010.5M, the AGN-corrected qIR measurements nearly coincide with those obtained from stacking non-detections alone (yellow squares). These latter values delimit the highest qIR that could be reached if removing, by definition, all radio detections. The result of similarity between the two sets of qIR 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 < 1010.5M 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 > 1010.5M (Sect. 4.2.3). Nevertheless, according to our analysis the vast majority of radio-detected dwarf galaxies (M < 109.5M, e.g., Mezcua 2017) in COSMOS are expected to be radio AGN.

Bearing this in mind, we note that the weighted average qIR (blue stars) are as yet mostly driven by non-detections (yellow squares), which outnumber individual detections (dots) by a factor of > 100 at M < 109.5M. 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 < 109.5M, the new average qIR still move upward by 0.2−0.3 dex.

4.3. The intrinsic IRRC evolves primarily with M

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

thumbnail Fig. 13.

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

Open with DEXTER

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

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

Figure 14 shows the equivalent of Fig. 13 but projected on M, with redshift bins in different colours. Error bars on each qIR were re-scaled by a factor of in each M bin, to bring the reduced χ2 of the corresponding qIR − z fit to unity. It is quite evident that an M dependence reduces substantially the scatter of the average qIR around a single trend. To better quantify this, first we bootstrapped over the uncertainties of slope and normalisation obtained from each qIR − z trend (see Table 4). Then, at each M, we interpolated the full range of bootstrapped IRRCs at z = 1, in correspondence of the 16th, 50th, and 84th percentiles. Interpolating at z = 1, besides being at roughly the median redshift of our sample, reduces the increasing divergence of each qIR − z fit at lower or higher redshifts. This leaves us with the interpolated median qIR(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 = 0.87, with an M slope close to that commonly found when fitting qIR 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.

thumbnail Fig. 14.

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

Open with DEXTER

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 qIR − M − z space. This yields the following best-fit expression:

(5)

where A = (1 + z) and B = (log M/M − 10). The corresponding = 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 Lradio–SFR) has also been highlighted in previous low-z studies (Gürkan et al. 2018; Read et al. 2018) and recently extended out to z ∼ 1 (Smith et al. 2021) using deep LOFAR-150 MHz data. A similar conclusion was independently reached in Molnár et al. (2021) when considering the dependence of the IRRC on galaxy spectral radio luminosity. To mitigate selection effects, they exploited a depth-matched sample of SFGs at z < 0.2. After performing a radio decomposition analysis in different bins of L1.4 GHz, Molnár et al. (2021) report that qIR decreases with increasing L1.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 qIR vs. z trends in the literature.

5. Discussion

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

5.1. What drives the primary M dependence?

Our main finding is that the IRRC decreases primarily with M, and only weakly with redshift. In particular, within the range 109 < M/M < 1012, the median qIR 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 L1.4 GHz − M is steeper than the dependence LIR − M (i.e. the MS). To translate this result into the corresponding IR-radio slope, we take our best qIR − M relation (Eq. (5)) at fixed redshift, and assume for simplicity a linear MS between M and SFR (i.e. LIR). This yields . 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 (, Bell 2003; Hodge et al. 2008; Davies et al. 2017; Brown et al. 2017; Gürkan et al. 2018; Molnár et al. 2021). This might challenge the idea of calibrating radio emission as a universal SFR tracer, as we discuss in Sect. 5.4.

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

5.1.1. The role of the IMF

We quantify whether a deviation from a canonical IMF slope (e.g., Chabrier 2003; n(M) ∝ M−2.35 at 0.8 < M < 100 M) could justify an M-decreasing qIR. 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 qIR.

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 qIR 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 LIR/L1.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) qIR 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 qIR with M.

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

The model proposed by Schleicher & Beck (2013) postulates that the magnetic field strength scales with SFR, boosting radio synchrotron emission during shocks or galaxy interactions (e.g., Donevski & Prodanović 2015; Tabatabaei et al. 2017). Because of the MS relation, this model implicitly predicts a net enhancement of radio emission with increasing M, as well as an increase with redshift due to higher gas density in galaxies. Related to this, the single-zone galaxy model of Lacki & Thompson (2010), which includes a detailed CR description, suggests a variation of qIR 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 qIR 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 qIR with ΣSFR and redshift.

Here, we test the above models by relating qIR and average ΣSFR measured in this work. These estimates were obtained by using the total SFRIR + 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 (Re, 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 Re, maj can be calculated as Re, maj = θFWHM/2.43 (Murphy et al. 2017). Lastly, we take the circularised radius Re = Re, maj/, where Ar is the axial ratio.

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

thumbnail Fig. 15.

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

Open with DEXTER

We calculate ΣSFR = SFRIR + UV/2 (see e.g., Jiménez-Andrade et al. 2019) and show its relation with qIR in Fig. 16, colour-coded by redshift (left panel) and M (right panel). Empty symbols highlight 7/37 bins with unresolved 3 GHz stacked emission, which translates into a lower limit in ΣSFR. By fitting only the 30 qIR 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 (qIR ∝ ( − 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).

thumbnail Fig. 16.

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

Open with DEXTER

Since more massive galaxies are characterised by more compact star formation (Elbaz et al. 2011), the decreasing qIR − Σ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 qIR − Σ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 qIR trend.

Both the slope and significance of the qIR − ΣSFR relation are consistent with the values found between qIR and M (Sect. 4.3). We argue that the declining qIR − Σ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 109 to 1011M (Fig. 7), while only increases by a factor of 1.5−2.5 in the same interval. Although this does not stand as conclusive evidence, our analysis seems to suggest that the larger SFR per unit area in more massive (and higher-z) galaxies might be driving the sub-linear behaviour of the IRRC.

The comparison with the model of Lacki & Thompson (2010) comes with a few caveats. First, we take a fixed spectral index α = −0.75 for all galaxies (which is supported by radio ancillary data in Appendix B), while Lacki & Thompson (2010) model a curved radio spectrum. Second, we label as ‘SBs’ those galaxies that lie > 4× above the MS (Sect. 5.2), while Lacki & Thompson (2010) identify them as ΣSFR ≳ 2−4 M yr−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 Lacki & Thompson (2010).

That said, their model predicts a decreasing qIR 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 qIR − ΣSFR trend and are nearly redshift-invariant. Our trend is consistent with the low-q values recently inferred by Algera et al. (2020a) in compact (Re ∼ 1 kpc) and massive (M > 1010M) sub-millimetre galaxies at 1.5 < z < 3.5. Indeed, their average qIR = 2.20 ± 0.03 lies close to the extrapolation of our best-fit qIR − ΣSFR trend at ΣSFR ∼ 100 M yr−1 kpc−2, further corroborating the relation between qIR and SFR per unit area in SFGs.

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

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

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

5.2. Does the IRRC evolve above the MS?

We investigate the behaviour of the average qIR above the MS. This is important for testing whether radio emission follows a similar enhancement as LIR when moving above the MS – otherwise, qIR 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 qIR 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 qIR 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 qFIR 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 × SFRMS (e.g., Rodighiero et al. 2011), where SFRMS 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/NIR > 3) that meet the above criterion. This is because our stacked SFRIR 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 > 1010.5M 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 LIR 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 qIR between the SB and MS populations. This leaves us with a sample of 554 SBs. As done for the full SFG sample, we performed median stacking at 3 GHz and combined the stacked signal with radio-detected SBs.

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

thumbnail Fig. 17.

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

Open with DEXTER

Although the (lack of) evolution of qIR above the MS is still debated, our decreasing qIR − ΣSFR trend (Fig. 16) would predict lower qIR 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 (Re ≲ 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 qIR 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 qIR = 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 qIR, while our results predominantly reflect the behaviour of the MS population.

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

5.3. Observing widespread AGN activity in radio-detected dwarfs

A noteworthy implication raised from our radio AGN subtraction is the possibility of widespread AGN activity within radio-detected dwarf galaxies (M < 109.5M). 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/NIR < 3) at any IR/sub-mm band (Fig. 3). This is likely a natural effect that is due to the increasing incompleteness of IR selection towards low M galaxies. From IR/sub-mm stacking, however, we obtain SFRIR > 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 L1.4 GHz values than their more massive counterparts, but still 100 times larger than those obtained from median radio stacking of non-detections. This effect fully counter-balances the high starburstiness seen in the IR, causing an overall drop of qIR in radio-detected dwarfs by over a factor of ten, with respect to the stacked population (see black dots relative to yellow squares in Figs. 9 and 12). These arguments let us suppose that radio-detected dwarfs are consistent with being AGN-dominated in the radio.

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

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

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

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

5.4. Assessing radio-synchrotron emission as a SFR tracer

In this section, we discuss the link between the IRRC and SFR in galaxies. As mentioned in Sect. 3.3, the conversion from LIR 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 LIR and L1.4 GHz is not, therefore, rigidly proportional to the SFR.

It is for this reason, we express qIR 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 (=SFRIR + UV). We thus define the parameter qSFR as:

(6)

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

We repeated the analogous qSFR decomposition analysis at M > 1010.5M to calibrate the AGN-vs.-SF locus of radio detections (Sect. 4.2). Within the two highest M bins, the best-fitting trend of qSFR with redshift has slope −0.057 ± 0.002, which is strikingly similar to that inferred for qIR (−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 qSFR − M − z plane yields the following expression:

(7)

where A = (1 + z) and B = (log M/M − 10). Similarly to the fit in the qIR 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 LIR, radio emission underestimates the total SFR by a larger factor as compared to the IR light. The sub-linear trend 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 LIR, that is, SFR. Such a radio deficit in the dwarf-galaxy regime could be possibly linked to shorter CRe scale heights (e.g., Helou & Bicay 1993; Lacki & Thompson 2010) or weaker magnetic fields (Donevski & Prodanović 2015; Tabatabaei et al. 2017) that are common in less dense SF environments.

Moreover, our adopted LIR–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 LIR. 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 qSFR − z trend or to amplify the M dependence of qSFR.

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 L1.4 GHz and SFR functions. Therefore, our results reinforce the need for M-dependent, non-linear calibrations between radio-continuum emission and SFR, in order to develop successful observing strategies for targeting low-M galaxies at radio wavelengths. These can be complemented with higher frequency observations that are more sensitive to thermal free-free emission as SFR tracer in high-redshift galaxies (e.g., Murphy et al. 2017; Penney et al. 2020; Van der Vlugt et al. 2021; Algera et al. 2020b).

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

6. Summary and conclusions

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

In each M − z bin, we performed stacking of undetected sources at both IR (Sect. 3.1) and radio (Sect. 3.4) frequencies and combined the stacked signal with individual detections a posteriori to infer the average qIR 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 > 1010.5M and extrapolated to the rest of the sample to infer the AGN-corrected IRRC (Sect. 4.3). Finally, we interpret our findings in the context of existing IRRC studies, from both models and observations. The main results of this work are listed below.

(1) The IRRC evolves primarily with M, with more massive galaxies displaying systematically lower qIR. A secondary, weaker dependence on redshift is also observed. The multi-parametric best-fitting 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). At fixed redshift, this trend translates into an IRRC of , 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 > 1010.5M 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 qIR − z trend at M > 1010.5M to a nearly flat slope. This correction nicely aligns the mode qIR of radio SFGs to the median stacked qIR of the full M sample of non-AGN galaxies. Therefore, we interpret the resulting AGN-corrected qIR 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 qIR 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 < 109.5M) 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 qIR 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 qIR in MS galaxies towards higher ΣSFR. Nevertheless, radio synchrotron models (e.g., Lacki & Thompson 2010; Schleicher & Beck 2013) predict a much stronger qIR 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 qIR 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 qIR, which is visibly at odds with our expectations. According to radio synchrotron models (Lacki & Thompson 2010), a ‘conspiracy’ of different factors might induce a qIR 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 qIR 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 qSFR − 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.


1

An exhaustive overview of the COSMOS field is available at: http://cosmos.astro.caltech.edu/

4

The far-infrared luminosity used to compute qFIR was integrated between 42 and 122 μm rest-frame. This is quantified to be 1.91× smaller than the total LIR (Magnelli et al. 2015).

Acknowledgments

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

References

  1. Aird, J., Coil, A. L., & Georgakakis, A. 2019, MNRAS, 484, 4360 [Google Scholar]
  2. Algera, H. S. B., Smail, I., Dudzevičiūtė, U., et al. 2020a, ApJ, 903, 138 [Google Scholar]
  3. Algera, H. S. B., Hodge, J. A., Riechers, D., et al. 2020b, ApJ, submitted [arXiv:2012.08499] [Google Scholar]
  4. Algera, H. S. B., van der Vlugt, D., Hodge, J. A., et al. 2020c, ApJ, 903, 139 [Google Scholar]
  5. Appleton, P. N., Fadda, D. T., Marleau, F. R., et al. 2004, ApJS, 154, 147 [Google Scholar]
  6. Aretxaga, I., Wilson, G. W., Aguilar, E., et al. 2011, MNRAS, 415, 3831 [Google Scholar]
  7. Arnouts, S., Walcher, C. J., Le Fèvre, O., et al. 2007, A&A, 476, 137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Baugh, C. M., Lacey, C. G., Frenk, C. S., et al. 2005, MNRAS, 356, 1191 [Google Scholar]
  9. Bavouzet, N., Dole, H., Le Floc’h, E., et al. 2008, A&A, 479, 83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bell, E. F. 2003, ApJ, 586, 794 [Google Scholar]
  11. Bertoldi, F., Carilli, C., Aravena, M., et al. 2007, ApJS, 172, 132 [Google Scholar]
  12. Best, P. N., & Heckman, T. M. 2012, MNRAS, 421, 1569 [Google Scholar]
  13. Béthermin, M., Dole, H., Beelen, A., & Aussel, H. 2010, A&A, 512, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Béthermin, M., Le Floc’h, E., Ilbert, O., et al. 2012, A&A, 542, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Bonaldi, A., Bonato, M., Galluzzi, V., et al. 2019, MNRAS, 482, 2 [Google Scholar]
  17. Bondi, M., Zamorani, G., Ciliegi, P., et al. 2018, A&A, 618, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Bonzini, M., Mainieri, V., Padovani, P., et al. 2015, MNRAS, 453, 1079 [Google Scholar]
  19. Bourne, N., Dunne, L., Ivison, R. J., et al. 2011, MNRAS, 410, 1155 [Google Scholar]
  20. Bourne, N., Maddox, S. J., Dunne, L., et al. 2012, MNRAS, 421, 3027 [Google Scholar]
  21. Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151 [Google Scholar]
  22. Brown, M. J. I., Moustakas, J., Kennicutt, R. C., et al. 2017, ApJ, 847, 136 [Google Scholar]
  23. Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [Google Scholar]
  24. Buat, V., Noll, S., Burgarella, D., et al. 2012, A&A, 545, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Burgarella, D., Buat, V., Gruppioni, C., et al. 2013, A&A, 554, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Calistro Rivera, G., Williams, W. L., Hardcastle, M. J., et al. 2017, MNRAS, 469, 3468 [Google Scholar]
  27. Capak, P., Aussel, H., Ajiki, M., et al. 2007, ApJS, 172, 99 [Google Scholar]
  28. Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2012, Nature, 484, 485 [Google Scholar]
  29. Carraro, R., Rodighiero, G., Cassata, P., et al. 2020, A&A, 642, A65 [EDP Sciences] [Google Scholar]
  30. Ceraj, L., Smolčić, V., Delvecchio, I., et al. 2018, A&A, 620, A192 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  32. Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62 [Google Scholar]
  33. Condon, J. J. 1992, ARA&A, 30, 575 [Google Scholar]
  34. Condon, J. J., Huang, Z. P., Yin, Q. F., & Thuan, T. X. 1991, ApJ, 378, 65 [Google Scholar]
  35. Condon, J. J., Helou, G., Sanders, D. B., & Soifer, B. T. 1993, AJ, 105, 1730 [Google Scholar]
  36. Condon, J. J., Cotton, W. D., & Broderick, J. J. 2002, AJ, 124, 675 [Google Scholar]
  37. Cowie, L. L., Barger, A. J., Hsu, L. Y., et al. 2017, ApJ, 837, 139 [Google Scholar]
  38. Cucciati, O., Tresse, L., Ilbert, O., et al. 2012, A&A, 539, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Dabringhausen, J., Kroupa, P., & Baumgardt, H. 2009, MNRAS, 394, 1529 [Google Scholar]
  40. Davé, R. 2008, MNRAS, 385, 147 [Google Scholar]
  41. Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Davies, L. J. M., Huynh, M. T., Hopkins, A. M., et al. 2017, MNRAS, 466, 2312 [Google Scholar]
  43. de Jong, T., Klein, U., Wielebinski, R., & Wunderlich, E. 1985, A&A, 147, L6 [NASA ADS] [Google Scholar]
  44. Deason, A., Wetzel, A., & Garrison-Kimmel, S. 2014, ApJ, 794, 115 [Google Scholar]
  45. Del Moro, A., Alexander, D. M., Mullaney, J. R., et al. 2013, A&A, 549, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Delhaize, J., Heywood, I., Prescott, M., et al. 2021, MNRAS, 501, 3833 [Google Scholar]
  48. Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2017, A&A, 602, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2018, MNRAS, 481, 4971 [Google Scholar]
  50. Delvecchio, I., Daddi, E., Shankar, F., et al. 2019, ApJ, 885, L36 [Google Scholar]
  51. Delvecchio, I., Daddi, E., Aird, J., et al. 2020, ApJ, 892, 17 [Google Scholar]
  52. Dole, H., Lagache, G., & Puget, J. L. 2003, ApJ, 585, 617 [Google Scholar]
  53. Donevski, D., & Prodanović, T. 2015, MNRAS, 453, 638 [Google Scholar]
  54. Donley, J. L., Rieke, G. H., Rigby, J. R., & Pérez-González, P. G. 2005, ApJ, 634, 169 [Google Scholar]
  55. Draine, B. T., & Li, A. 2007, ApJ, 657, 810 [Google Scholar]
  56. Dubois, Y., Volonteri, M., Silk, J., et al. 2015, MNRAS, 452, 1502 [Google Scholar]
  57. Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Fakhouri, O., Ma, C.-P., & Boylan-Kolchin, M. 2010, MNRAS, 406, 2267 [Google Scholar]
  59. Franco, M., Elbaz, D., Zhou, L., et al. 2020, A&A, 643, A53 [EDP Sciences] [Google Scholar]
  60. Garrett, M. A. 2002, A&A, 384, L19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Garrison-Kimmel, S., Rocha, M., Boylan-Kolchin, M., Bullock, J. S., & Lally, J. 2013, MNRAS, 433, 3539 [Google Scholar]
  62. Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789 [Google Scholar]
  63. Genzel, R., Burkert, A., Bouché, N., et al. 2008, ApJ, 687, 59 [Google Scholar]
  64. Goulding, A. D., Forman, W. R., Hickox, R. C., et al. 2014, ApJ, 783, 40 [Google Scholar]
  65. Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3 [EDP Sciences] [Google Scholar]
  66. Gürkan, G., Hardcastle, M. J., Smith, D. J. B., et al. 2018, MNRAS, 475, 3010 [Google Scholar]
  67. Hales, C. A., Murphy, T., Curran, J. R., et al. 2012, MNRAS, 425, 979 [Google Scholar]
  68. Hales, C. A., Norris, R. P., Gaensler, B. M., et al. 2014, MNRAS, 441, 2555 [Google Scholar]
  69. Hao, C.-N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124 [Google Scholar]
  70. Hardcastle, M. J., & Croston, J. H. 2020, New Astron. Rev., 88, 101539 [Google Scholar]
  71. Harwit, M., & Pacini, F. 1975, ApJ, 200, L127 [Google Scholar]
  72. Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589 [Google Scholar]
  73. Heinis, S., Buat, V., Béthermin, M., et al. 2013, MNRAS, 429, 1113 [Google Scholar]
  74. Heinis, S., Buat, V., Béthermin, M., et al. 2014, MNRAS, 437, 1268 [Google Scholar]
  75. Helou, G., & Bicay, M. D. 1993, ApJ, 415, 93 [Google Scholar]
  76. Helou, G., Soifer, B. T., & Rowan-Robinson, M. 1985, ApJ, 298, L7 [Google Scholar]
  77. Herrera Ruiz, N., Middelberg, E., Deller, A., et al. 2017, A&A, 607, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Herrera Ruiz, N., Middelberg, E., Deller, A., et al. 2018, A&A, 616, A128 [EDP Sciences] [Google Scholar]
  79. Hickox, R. C., Jones, C., Forman, W. R., et al. 2009, ApJ, 696, 891 [Google Scholar]
  80. Hodge, J. A., Becker, R. H., White, R. L., & de Vries, W. H. 2008, AJ, 136, 1097 [Google Scholar]
  81. Hopkins, A. M., & Beacom, J. F. 2006, ApJ, 651, 142 [Google Scholar]
  82. Hummel, E., Davies, R. D., Wolstencroft, R. D., van der Hulst, J. M., & Pedlar, A. 1988, A&A, 199, 91 [NASA ADS] [Google Scholar]
  83. Ibar, E., Cirasuolo, M., Ivison, R., et al. 2008, MNRAS, 386, 953 [NASA ADS] [CrossRef] [Google Scholar]
  84. Ibar, E., Ivison, R. J., Biggs, A. D., et al. 2009, MNRAS, 397, 281 [Google Scholar]
  85. Ibar, E., Ivison, R. J., Best, P. N., et al. 2010, MNRAS, 401, L53 [NASA ADS] [CrossRef] [Google Scholar]
  86. Ivison, R. J., Magnelli, B., Ibar, E., et al. 2010a, A&A, 518, L31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Ivison, R. J., Swinbank, A. M., Swinyard, B., et al. 2010b, A&A, 518, L35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Jarvis, M. J., Smith, D. J. B., Bonfield, D. G., et al. 2010, MNRAS, 409, 92 [Google Scholar]
  89. Jarvis, M., Taylor, R., Agudo, I., et al. 2016, MeerKAT Science: On the Pathway to the SKA, 6 [Google Scholar]
  90. Jiménez-Andrade, E. F., Magnelli, B., Karim, A., et al. 2019, A&A, 625, A114 [EDP Sciences] [Google Scholar]
  91. Jin, S., Daddi, E., Liu, D., et al. 2018, ApJ, 864, 56 [Google Scholar]
  92. Karim, A., Schinnerer, E., Martínez-Sansigre, A., et al. 2011, ApJ, 730, 61 [Google Scholar]
  93. Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739 [NASA ADS] [Google Scholar]
  94. Kaviraj, S., Martin, G., & Silk, J. 2019, MNRAS, 489, L12 [Google Scholar]
  95. Keller, B. W., Wadsley, J., & Couchman, H. M. P. 2016, MNRAS, 463, 1431 [Google Scholar]
  96. Kennicutt, R. C., Jr. 1998, ApJ, 498, 541 [Google Scholar]
  97. Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531 [Google Scholar]
  98. Koudmani, S., Henden, N. A., & Sijacki, D. 2020, MNRAS, submitted [arXiv:2007.10342] [Google Scholar]
  99. Kurczynski, P., & Gawiser, E. 2010, AJ, 139, 1592 [Google Scholar]
  100. Lacki, B. C., & Thompson, T. A. 2010, ApJ, 717, 196 [Google Scholar]
  101. Lacki, B. C., Thompson, T. A., & Quataert, E. 2010, ApJ, 717, 1 [Google Scholar]
  102. Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24 [Google Scholar]
  103. Le Floc’h, E., Aussel, H., Ilbert, O., et al. 2009, ApJ, 703, 222 [Google Scholar]
  104. Lee, N., Sanders, D. B., Casey, C. M., et al. 2015, ApJ, 801, 80 [Google Scholar]
  105. Lehmer, B. D., Basu-Zych, A. R., Mineo, S., et al. 2016, ApJ, 825, 7 [Google Scholar]
  106. Leslie, S. K., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58 [Google Scholar]
  107. Liu, D., Daddi, E., Dickinson, M., et al. 2018, ApJ, 853, 172 [Google Scholar]
  108. Lutz, D. 2014, ARA&A, 52, 373 [Google Scholar]
  109. Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, A&A, 532, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
  111. Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6 [Google Scholar]
  112. Magnelli, B., Elbaz, D., Chary, R. R., et al. 2009, A&A, 496, 57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  114. Magnelli, B., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Magnelli, B., Ivison, R. J., Lutz, D., et al. 2015, A&A, 573, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Mancuso, C., Lapi, A., Cai, Z.-Y., et al. 2015, ApJ, 810, 72 [Google Scholar]
  117. Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115 [Google Scholar]
  118. Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34 [Google Scholar]
  119. Marleau, F. R., Clancy, D., Habas, R., & Bianconi, M. 2017, A&A, 602, A28 [EDP Sciences] [Google Scholar]
  120. Massardi, M., Bonaldi, A., Negrello, M., et al. 2010, MNRAS, 404, 532 [NASA ADS] [Google Scholar]
  121. Mauch, T., Cotton, W. D., Condon, J. J., et al. 2020, ApJ, 888, 61 [Google Scholar]
  122. McCracken, H. J., Le Fèvre, O., Brodwin, M., et al. 2001, A&A, 376, 756 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  123. McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  124. Mezcua, M. 2017, Int. J. Mod. Phys. D, 26, 1730021 [Google Scholar]
  125. Mezcua, M., Civano, F., Fabbiano, G., Miyaji, T., & Marchesi, S. 2016, ApJ, 817, 20 [Google Scholar]
  126. Mezcua, M., Suh, H., & Civano, F. 2019, MNRAS, 488, 685 [Google Scholar]
  127. Mezcua, M., & Domínguez Sánchez, H. 2020, ApJ, 898, L30 [Google Scholar]
  128. Molnár, D. C., Sargent, M. T., Delhaize, J., et al. 2018, MNRAS, 475, 827 [Google Scholar]
  129. Molnár, D. C., Sargent, M. T., Leslie, S., et al. 2021, MNRAS, submitted [arXiv:2103.04803] [Google Scholar]
  130. Mullaney, J. R., Alexander, D. M., Goulding, A. D., & Hickox, R. C. 2011, MNRAS, 414, 1082 [Google Scholar]
  131. Murphy, E. J. 2009, ApJ, 706, 482 [Google Scholar]
  132. Murphy, E. J. 2013, ApJ, 777, 58 [Google Scholar]
  133. Murphy, E. J., Helou, G., Kenney, J. D. P., Armus, L., & Braun, R. 2008, ApJ, 678, 828 [Google Scholar]
  134. Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67 [Google Scholar]
  135. Murphy, E. J., Bremseth, J., Mason, B. S., et al. 2012, ApJ, 761, 97 [Google Scholar]
  136. Murphy, E. J., Stierwalt, S., Armus, L., Condon, J. J., & Evans, A. S. 2013, ApJ, 768, 2 [Google Scholar]
  137. Murphy, E. J., Momjian, E., Condon, J. J., et al. 2017, ApJ, 839, 35 [Google Scholar]
  138. Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJS, 206, 8 [Google Scholar]
  139. Nelson, E. J., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2016, ApJ, 828, 27 [Google Scholar]
  140. Nguyen, H. T., Schulz, B., Levenson, L., et al. 2010, A&A, 518, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  141. Niklas, S., & Beck, R. 1997, A&A, 320, 54 [NASA ADS] [Google Scholar]
  142. Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43 [Google Scholar]
  143. Novak, M., Smolčić, V., Delhaize, J., et al. 2017, A&A, 602, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  144. Novak, M., Smolčić, V., Schinnerer, E., et al. 2018, A&A, 614, A47 [EDP Sciences] [Google Scholar]
  145. Oke, J. B. 1974, ApJS, 27, 21 [Google Scholar]
  146. Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614 [Google Scholar]
  147. Padovani, P., Bonzini, M., Kellermann, K. I., et al. 2015, MNRAS, 452, 1263 [Google Scholar]
  148. Pannella, M., Elbaz, D., Daddi, E., et al. 2015, ApJ, 807, 141 [Google Scholar]
  149. Penney, J. I., Blain, A. W., Assef, R. J., et al. 2020, MNRAS, 496, 1565 [Google Scholar]
  150. Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  151. Popesso, P., Magnelli, B., Buttiglione, S., et al. 2012, ArXiv e-prints [arXiv:1211.4257] [Google Scholar]
  152. Read, S. C., Smith, D. J. B., Gürkan, G., et al. 2018, MNRAS, 480, 5625 [Google Scholar]
  153. Reines, A. E., & Deller, A. T. 2012, ApJ, 750, L24 [Google Scholar]
  154. Reines, A. E., Sivakoff, G. R., Johnson, K. E., & Brogan, C. L. 2011, Nature, 470, 66 [Google Scholar]
  155. Reines, A. E., Greene, J. E., & Geha, M. 2013, ApJ, 775, 116 [Google Scholar]
  156. Reines, A. E., Plotkin, R. M., Russell, T. D., et al. 2014, ApJ, 787, L30 [Google Scholar]
  157. Rickard, L. J., & Harvey, P. M. 1984, AJ, 89, 1520 [Google Scholar]
  158. Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40 [Google Scholar]
  159. Salim, S., Charlot, S., Rich, R. M., et al. 2005, ApJ, 619, L39 [Google Scholar]
  160. Sargent, M. T., Schinnerer, E., Murphy, E., et al. 2010, ApJ, 714, L190 [Google Scholar]
  161. Sargent, M. T., Béthermin, M., Daddi, E., & Elbaz, D. 2012, ApJ, 747, L31 [Google Scholar]
  162. Schinnerer, E., Sargent, M. T., Bondi, M., et al. 2010, ApJS, 188, 384 [Google Scholar]
  163. Schleicher, D. R. G., & Beck, R. 2013, A&A, 556, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  164. Schober, J., Schleicher, D. R. G., & Klessen, R. S. 2016, ApJ, 827, 109 [Google Scholar]
  165. Schober, J., Schleicher, D. R. G., & Klessen, R. S. 2017, MNRAS, 468, 946 [Google Scholar]
  166. Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  167. Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1 [Google Scholar]
  168. Silk, J. 2017, ApJ, 839, L13 [Google Scholar]
  169. Smith, D. J. B., Jarvis, M. J., Hardcastle, M. J., et al. 2014, MNRAS, 445, 2232 [Google Scholar]
  170. Smith, D. J. B., Haskell, P., Gürkan, G., et al. 2021, A&A, 648, A6 [EDP Sciences] [Google Scholar]
  171. Smolčić, V., Novak, M., Bondi, M., et al. 2017a, A&A, 602, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  172. Smolčić, V., Delvecchio, I., Zamorani, G., et al. 2017b, A&A, 602, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  173. Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15 [Google Scholar]
  174. Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175 [Google Scholar]
  175. Steinhardt, C. L., Speagle, J. S., Capak, P., et al. 2014, ApJ, 791, L25 [Google Scholar]
  176. Tabatabaei, F. S., Schinnerer, E., Krause, M., et al. 2017, ApJ, 836, 185 [Google Scholar]
  177. Van der Vlugt, D., Algera, H. S. B., Hodge, J. A., et al. 2021, ApJ, 907, 5 [Google Scholar]
  178. van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28 [Google Scholar]
  179. van Dokkum, P. G. 2008, ApJ, 674, 29 [Google Scholar]
  180. Viero, M. P., Wang, L., Zemcov, M., et al. 2013, ApJ, 772, 77 [Google Scholar]
  181. Voelk, H. J. 1989, A&A, 218, 67 [NASA ADS] [Google Scholar]
  182. Welikala, N., Béthermin, M., Guery, D., et al. 2016, MNRAS, 455, 1629 [Google Scholar]
  183. Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29 [Google Scholar]
  184. Whitaker, K. E., Pope, A., Cybulski, R., et al. 2017, ApJ, 850, 208 [Google Scholar]
  185. White, R. L., Helfand, D. J., Becker, R. H., Glikman, E., & de Vries, W. 2007, ApJ, 654, 99 [Google Scholar]
  186. White, S. V., Jarvis, M. J., Häußler, B., & Maddox, N. 2015, MNRAS, 448, 2665 [Google Scholar]
  187. White, S. V., Jarvis, M. J., Kalfountzou, E., et al. 2017, MNRAS, 468, 217 [Google Scholar]
  188. Wong, O. I., Koss, M. J., Schawinski, K., et al. 2016, MNRAS, 460, 1588 [Google Scholar]
  189. Wu, Y., Charmandaris, V., Houck, J. R., et al. 2008, ApJ, 676, 970 [Google Scholar]
  190. Yun, M. S., Reddy, N. A., & Condon, J. J. 2001, ApJ, 554, 803 [Google Scholar]

Appendix A: Testing total radio flux densities

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

thumbnail Fig. A.1.

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

Open with DEXTER

We further demonstrate that our choice of performing median stacking at 3 GHz, rather than rms-weighted mean stacking, does not impact our final L1.4 GHz estimates. A comparison between median and mean L1.4 GHz is presented in Fig. A.2. The top panel displays L1.4 GHz from the combined flux of detections and non-detections (see Eq. (3)), while the bottom panel refers to the case of purely undetected sources. Colours indicate different M bins. Only stacks in which peak flux densities have S/N > 3 are shown. No systematics is observed, at any M, between mean and median stacked L1.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.

thumbnail Fig. A.2.

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

Open with DEXTER

The fact that non-detections (bottom panel) display consistent L1.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, L1.4 GHz of radio detections (top panel) are always > 3× larger than L1.4 GHz of non-detections (bottom panel), despite the smaller numbers. This partly smooths over the initial fluctuations between mean and median stacking, thus delivering an even tighter agreement, as we observe.

Appendix B: Stacking ancillary VLA and MIGHTEE data

Here we perform a radio stacking analysis, as for 3 GHz data (Sect. 3.4), in order to check whether our 3 GHz based L1.4 GHz are stable against different resolutions or spectral frequencies. We exploit VLA data from the 1.4 GHz Deep Project (Schinnerer et al. 2010) map. It covers 1.7 deg2 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 deg2 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 Ks + 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 deg2) were stacked, containing roughly half of the sample size used for VLA stacking.

The stacking analysis follows the same reasoning and assumptions presented in Sect. 3.4. Stacked MIGHTEE flux densities are measured in the central pixel, that is assumed to trace the total flux. VLA 1.4 GHz peak flux densities were, instead, scaled to total flux densities as done at 3 GHz. Nonetheless, a different, yet empirical relation was adopted to identify resolved sources at VLA 1.4 GHz, calibrated on 1.4 GHz detections (Schinnerer et al. 2010): Stot/. 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 L1.4 GHz assuming α = −0.75 ± 0.1, that was propagated along with flux errors to deliver reliable L1.4 GHz uncertainties. Upper limits at 3σ were assigned for S/N < 3 stacks.

Figure B.1 shows stacked cutouts at 0.8 < z < 1.2 at VLA 3 GHz (top), VLA 1.4 GHz (middle) and MIGHTEE 1.3 GHz (bottom) data, as a function of M (increasing from left to right). While stacking at 3 GHz delivers S/N > 3 in 39/42 bins, only 25/42 and 29/42 have S/N > 3 in VLA 1.4 GHz and MIGHTEE 1.3 GHz stacked images, respectively. For VLA 1.4 GHz, the small number of bins is attributed to shallower than 3 GHz observations. For MIGHTEE, instead, this is probably induced by the confusion-limited signal in the stacks due to the larger MeerKAT primary beam at 1.3 GHz (e.g., Mauch et al. 2020). Nevertheless, because VLA 3 GHz data are much less sensitive than MeerKAT to large-scale radio emission, total radio flux densities might be underestimated at 3 GHz. This issue can be, however, especially relevant at low redshift (z < 0.5) and for resolved or multi-component radio sources (e.g., Delhaize et al. 2021). In fact, a visual inspection of the median 3 GHz stacks of non-detections does not clearly reveal any missing flux in the residual images at the scales of the MIGHTEE beam, except in the bin at z < 0.5 and 1011 < M/M < 1012. 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 L1.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 > 1011M (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.

thumbnail Fig. B.1.

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

Open with DEXTER

Considerations on the radio spectral index

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

thumbnail Fig. B.2.

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

Open with DEXTER

Magnelli et al. (2015) measured the average spectral index exploiting VLA 1.4 GHz and GMRT 610 MHz data for an M-selected galaxy sample. They found that the observed 610 MHz−1.4 GHz slope, which probes closer to the rest-frame of 1.4 GHz than our 3 GHz data, does not seem to change with M or SFR, at least out to z ∼ 2. More recently, Calistro Rivera et al. (2017) exploited Low Frequency Array (LOFAR) data at 150 MHz in the Boötes field, out to z ∼ 2.5. Interestingly, they observed a spectral flattening of the radio SED of SFGs in the observed range [150 MHz−1.4 GHz] (see also Read et al. 2018; Gürkan et al. 2018). However, they argue that this feature should not affect the k-correction for the rest-frame 1.4 GHz luminosities L1.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 qIR distribution (qIR, 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 qIR, 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 qIR, here, we follow a different approach.

The observed qIR 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 qIR, peak, but we leave it free to vary along with the dispersion and normalisation for both functions. In this simultaneous fitting we give equal input weights to all bins, regardless of the number of sources in each. This approach is thus expected to return a rather conservative AGN contribution relative to SFGs.

The results are shown in Fig. C.1 for the two highest M bins. As in Fig. 10, the best-fit SF (blue) and AGN (red) Gaussians head up to reproduce the total distribution (black). However, we clearly notice two main differences compared to the previous approach. Firstly, the AGN distribution is far broader than the SF distribution in both M bins. Secondly, the relative fraction of radio AGN that we mis-classify as SFGs (red tail at > qpeak − 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 1010.5 < M/M < 1011 and 90% at at 1011 < M/M < 1012.

thumbnail Fig. C.1.

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

Open with DEXTER

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 qIR, peak and the dispersion reaching ∼0.21 dex. Therefore, it is reasonable to assume that the mode of the observed qIR distribution is attributed to radio-detected SFGs. Related to this, the threshold qpeak − 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 qIR by a larger amount towards lower M galaxies. Indeed, at 1011 < M/M < 1012 the radio AGN distribution is clearly broader, but far more offset than at 1010.5 < M/M < 1011, 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 qIR reported in this work.

Appendix D: Differences compared to the literature

Our best-fit relation of qIR as a function of M and redshift (Eq. (5) in Sect. 4.3) is fully consistent with the average qIR value measured in local SFGs (i.e. 2.64 in Bell 2003) for a typical galaxy with M ∼ 1010M. At higher redshifts, instead, our average qIR 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 > 1010.5M. 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 qIR 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 qIR 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 L1.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 LIR estimates and Magnelli et al. (2015)’s were instead calculated through a fully consistent approach, therefore only a difference in L1.4 GHz might lead to systematics in the final qIR 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 L1.4 GHz measurements than those of Magnelli et al. (2015). This might be ascribed to the contribution of radio-detected SFGs to our weighted average L1.4 GHz, since they make a substantial fraction of the M-selected sample in that M bin (∼45%, Table 4). Therefore, this test proves our radio AGN subtraction broadly consistent with a totally independent approach.

thumbnail Fig. D.1.

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

Open with DEXTER

D.2. Sample selection and binning

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

We split the sample of Delhaize et al. (2017) among the same seven redshift bins analyzed in this work. For each, we perform median stacking of 3 GHz and IR images in all bands, combining both detections and non-detections. This approach should be comparable to the search for the median value carried out via survival analysis (Delhaize et al. 2017). Although we do not formally remove radio AGN in this check, we showed in Appendix D.1 that median radio stacking yields broadly consistent results (see Magnelli et al. 2015). Figure D.2 displays the median qIR obtained by stacking the SFG sample of Delhaize et al. (2017) in different redshift bins (stars). This yields a best-fitting qIR ∝ (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.

thumbnail Fig. D.2.

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

Open with DEXTER

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 > 1010.5M (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 qIR with M, binning only as a function of redshift leads to a mixture of different galaxy M that results into a larger global dispersion.

All Tables

Table 1.

Main numbers of priors and detections that characterise our final sample of 413 678 star-forming galaxies.

Table 2.

Average fraction of clustering signal at each FIR/sub-mm band.

Table 3.

Comparison between fractions of radio-SFGs and radio-excess AGN below some threshold qthres, for the two highest M bins.

Table 4.

Summary of the numbers and fractions of radio AGN and SFGs in different M bins after fitting the AGN-corrected qIR with redshift (Sect. 4.2.4 and Fig. 12).

All Figures

thumbnail Fig. 1.

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

Open with DEXTER
In the text
thumbnail Fig. 2.

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

Open with DEXTER
In the text
thumbnail Fig. 3.

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

Open with DEXTER
In the text
thumbnail Fig. 4.

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

Open with DEXTER
In the text
thumbnail Fig. 5.

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

Open with DEXTER
In the text
thumbnail Fig. 6.

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

Open with DEXTER
In the text
thumbnail Fig. 7.

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

Open with DEXTER
In the text
thumbnail Fig. 8.

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

Open with DEXTER
In the text
thumbnail Fig. 9.

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

Open with DEXTER
In the text
thumbnail Fig. 10.

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

Open with DEXTER
In the text
thumbnail Fig. 11.

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

Open with DEXTER
In the text
thumbnail Fig. 12.

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

Open with DEXTER
In the text
thumbnail Fig. 13.

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

Open with DEXTER
In the text
thumbnail Fig. 14.

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

Open with DEXTER
In the text
thumbnail Fig. 15.

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

Open with DEXTER
In the text
thumbnail Fig. 16.

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

Open with DEXTER
In the text
thumbnail Fig. 17.

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

Open with DEXTER
In the text
thumbnail Fig. A.1.

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

Open with DEXTER
In the text
thumbnail Fig. A.2.

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

Open with DEXTER
In the text
thumbnail Fig. B.1.

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

Open with DEXTER
In the text
thumbnail Fig. B.2.

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

Open with DEXTER
In the text
thumbnail Fig. C.1.

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

Open with DEXTER
In the text
thumbnail Fig. D.1.

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

Open with DEXTER
In the text
thumbnail Fig. D.2.

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

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.