© ESO, 2014

1. Introduction

The far-infrared (FIR) and radio luminosities of star-forming galaxies are tightly related via an empirical relation, the FIR/radio correlation (FRC; e.g. de Jong et al. 1985; Helou et al. 1985, 1988; Condon 1992; Yun et al. 2001). In the local Universe, this FRC is roughly linear across three orders of magnitude in luminosity, 10⁹ ≲ L_IR(8−1000 μm) [ L_⊙ ] ≲ 10^12.5, from dwarf galaxies to ultra-luminous infrared galaxies (ULIRGs; L_IR > 10¹² L_⊙). At lower luminosities, the FRC exhibits some signs of non-linearity (e.g. Bell et al. 2003). The FRC not only exists on galactic scales but also holds in star-forming regions within the galaxies, down to at least 0.5 kpc, albeit with some variations, as shown by resolved studies of nearby galaxies (e.g. Beck & Golla 1988; Bicay & Helou 1990; Murphy et al. 2008; Dumas et al. 2011; Tabatabaei et al. 2013a,b).

The FRC is thought to be caused by star-formation activity in galaxies. Young massive stars (≳8 M_⊙) radiate most of their energy at UV wavelengths, photons that dominate the UV continuum luminosity of galaxies and which are absorbed and re-emitted in the FIR regime by dust. After several Myrs, these young massive stars explode as supernovae (SNe), accelerating cosmic rays (CRs) into the general magnetic field of their host galaxy and resulting in diffuse synchrotron emission. Averaged over a star-formation episode, massive stars thus provide a common origin for the FIR and synchrotron emission of galaxies. This shared origin is the essence of all models constructed to explain the FRC on both global scales (e.g. the calorimeter model, Völk 1989; the conspiracy model, Lacki et al. 2010) and local scales (e.g. small-scale dynamo, Schleicher & Beck 2013; see also Niklas & Beck 1997). While many of these models are able to reproduce the basic properties of the FRC, none are consistent with all constraints provided by observations.

Despite our lack of a thorough theoretical understanding, the very tight empirical FRC (with a dispersion of s0.25 dex; e.g. Yun et al. 2001) has been used extensively for a number of astrophysical purposes. For example, the FRC observed in the local Universe has been exploited to empirically calibrate the radio luminosity as a star-formation rate (SFR) indicator, using the known L_IR-SFR relation for galaxies whose activities are not dominated by active galaxy nuclei (AGN; e.g. Condon 1992; Bell et al. 2003; Murphy et al. 2011). This allows us to estimate the level of star formation in high-redshift galaxies, taking advantage of radio observations that are often deeper than FIR surveys with better spatial resolution. The locally observed FRC has also been used to identify samples of radio-loud AGNs and to study their properties (e.g. Donley et al. 2005; Park et al. 2008; Del Moro et al. 2013). Finally, at high redshift, the local FRC has been used to estimate the distance (e.g. Carilli & Yun 1999) or the dust temperature (e.g. Chapman et al. 2005) of luminous starbursts, namely the submillimetre galaxies (SMGs; Smail et al. 1997). As these examples demonstrate, the use of the FRC has become an important tool for extragalactic astrophysics. Upcoming surveys with the Jansky Very Large Array (JVLA) will make its applications even more important.

Although the FRC is characterised well at low redshift, its form and thus its applicability at high redshift still have to be firmly demonstrated. From a theoretical point of view, we expect the FRC to break down at high redshift (i.e. z ≳ 2−3) because CR electron cooling via inverse Compton (IC) scattering off the cosmic microwave background (CMB; U_CMB ∝ (1 + z)⁴) photons is supposed to dominate over synchrotron cooling (e.g. Murphy 2009; Lacki & Thompson 2010; Schleicher & Beck 2013). However, the exact redshift and amplitude of this breakdown varies with models and with the assumed properties of high-redshift galaxies. From the observational point of view, the characteristics of the FRC at high redshift have been subject to extensive debate. Some studies have found that the FRC stays unchanged or suffers only minor variations at high redshift (e.g. Appleton et al. 2004; Ibar et al. 2008; Bourne et al. 2011) and others have found significant evolution of the FRC in the bulk of the star-forming galaxy population (e.g. Seymour et al. 2009) or for a subsample of it (e.g. SMGs; Murphy et al. 2009). Sargent et al. (2010) argue that these discrepant measurements could be explained by selection biases due to the improper treatment of flux limits from non-detections amongst radio- and FIR-selected samples. Then, applying a survival analysis to properly treat these non-detections, they conclude that the FRC remains unchanged or with little variations out to zs1.4. However, the results of Sargent et al. (2010) were still limited by relatively sparse coverage of the FIR and radio spectral energy distributions (SEDs) of high-redshift galaxies. Noticing these limitations, Ivison et al. (2010a) studied the evolution of the FRC using FIR (250, 350 and 500μm) observations from BLAST (Devlin et al. 2009) and multi-frequency radio observations (1.4 GHz and 610 MHz). Starting from a mid-infrared-selected sample and accounting for radio non-detections using a stacking analysis, they found an evolution of the FRC with redshift as ∝(1 + z)^{−0.15 ± 0.03}. Repeating a similar analysis using early observations from the Herschel Space Observatory, Ivison et al. (2010b) found support for such redshift evolution in the FRC. These findings demand modifications to any high-redshift results that adopted the local FRC. However, because results from Ivison et al. (2010a,b) were based on mid-infrared-selected samples, they might still be affected by some selection biases.

In this paper, we aim to study the FRC across 0 <z< 2.3, avoiding the biases mentioned above. To obtain a good FIR and radio spectral coverage, we use deep FIR (100, 160, 250, 350 and 500μm) observations from the Herschel Space Observatory (Pilbratt et al. 2010) and deep radio 1.4 GHz VLA and 610 MHz Giant Metre-wave Radio Telescope (GMRT) observations. To minimise selection biases and properly account for radio and FIR non-detections, we use a careful stacking analysis based on the positions of stellar-mass-selected samples of galaxies. These stellar-mass-selected samples are complete for star-forming galaxies down to 10¹⁰ M_⊙ across 0 <z< 2.3 and are built from the large wealth of multi-wavelength observations available for the blank extragalactic fields used in our analysis – GOODS-N, GOODS-S, ECDFS and COSMOS.

Besides investigating the evolution of the FRC as a function of redshift, we also aim to study its evolution in the SFR-stellar mass (M_∗) plane. Indeed, recent results have shown that the physical properties (e.g. morphology, CO-to-H₂ conversion factor, dust temperature) of star-forming galaxies vary with their positions in the SFR–M_∗ plane (e.g. Wuyts et al. 2011b; Magnelli et al. 2012b, 2014). In particular, these properties correlate with the distance of a galaxy to the so-called main sequence (MS) of the SFR–M_∗ plane, i.e. the sequence where the bulk of the star-forming galaxy population resides and which is characterised by , with 0.5 <γ< 1.0 (Brinchmann et al. 2004; Schiminovich et al. 2007; Noeske et al. 2007; Elbaz et al. 2007; Daddi et al. 2007; Pannella et al. 2009; Dunne et al. 2009; Rodighiero et al. 2010; Oliver et al. 2010; Karim et al. 2011; Mancini et al. 2011; Whitaker et al. 2012). Studying the FRC as a function of the distance of a galaxy from the MS (i.e. Δlog (SSFR)_MS = log  [ SSFR(galaxy) /SSFR_MS(M_∗,z) ]) will allow us to estimate if this correlation evolves from normal star-forming galaxies (Δlog (SSFR)_MSs0) to starbursting galaxies (Δlog (SSFR)_MSs1), as suggested by some local observations (e.g. Condon et al. 1991; but see, e.g. Yun et al. 2001).

To study the evolution of the FRC as a function of redshift and as a function of the distance of a galaxy from the MS (i.e. Δlog (SSFR)_MS), we grid the SFR–M_∗ plane in several redshift ranges and estimate for each SFR–M_∗–z bin its infrared luminosity, radio spectral index (i.e. α, where S_ν ∝ ν^−α) and radio luminosity using a FIR and radio stacking analysis. Thanks to this methodology, we are able to statistically and accurately constrain the radio spectral index and FRC of all star-forming galaxies with M_∗ > 10¹⁰ M_⊙, Δlog (SSFR)_MS > −0.3 and 0 <z< 2.3.

The paper is structured as follows. In Sect. 2 we present the Herschel, VLA and GMRT observations used in this study, as well as our stellar-mass-selected sample. Section 3 presents the FIR and radio stacking analysis used to infer the infrared luminosity, radio spectral index and radio luminosity of each of our SFR–M_∗–z bins. Evolution of the radio spectral index, , and of the FRC with Δlog (SSFR)_MS and redshift are presented in Sects. 4.1 and 4.2, respectively. Results are discussed in Sect. 5 and summarised in Sect. 6.

Throughout the paper we use a cosmology with H₀ = 71 km s^-1 Mpc^-1, Ω_Λ = 0.73 and Ω_M = 0.27.

2. Data

2.1. Herschel observations

To study the FIR properties of galaxies in the SFR−M_∗ plane, we used deep FIR observations of the COSMOS, GOODS-N, GOODS-S and ECDFS fields provided by the Herschel Space observatory. Observations at 100 and 160μm were obtained by the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) as part of the PACS Evolutionary Probe (PEP¹; Lutz et al. 2011) guaranteed time key programme and the GOODS-Herschel (GOODS-H²; Elbaz et al. 2011) open time key programme. Observations at 250, 350 and 500μm were obtained by the Spectral and Photometric Imaging Receiver (SPIRE Griffin et al. 2010) as part of the Herschel Multi-tiered Extragalactic Survey (HerMES³Oliver et al. 2012). The PEP, PEP/GOODS-H and HerMES surveys and data reduction methods are described in Lutz et al. (2011), Magnelli et al. (2013, see also Elbaz et al. 2011 and Oliver et al. (2012), respectively. Here we only summarise the information relevant for our study.

Herschel flux densities were derived using a point-spread-function-fitting method, guided by the known position of sources detected in deep 24 μm observations from the Multiband Imaging Photometer (MIPS; Rieke et al. 2004) on board the Spitzer Space Observatory. This method provides reliable and highly complete PACS/SPIRE source catalogues (Lutz et al. 2011; Magnelli et al. 2013). We note that the small fraction of Herschel sources without a MIPS counterpart and thus missing in our source catalogues (Magdis et al. 2011) will be considered via our stacking analysis which is based on the positions of complete stellar-mass-selected samples (Sect. 3.1). The extraction of PACS sources was accomplished using the method described in Magnelli et al. (2013), while for SPIRE sources it was done using the method described in Roseboom et al. (2010), both using the same 24 μm catalogues. In GOODS-S/N, we used the GOODS 24 μm catalogue, reaching a 3σ limit of 20μJy (Magnelli et al. 2009, 2011); in the ECDFS, we used the FIDEL 24 μm catalogue, reaching a 3σ limit of 70μJy (Magnelli et al. 2009); in COSMOS, we used a 24 μm catalogue with a 3σ limit of 45μJy (Le Floc’h et al. 2009). Reliability, completeness and contamination of our PACS/SPIRE catalogues were tested using Monte Carlo simulations. Table 1 summarises the depths of all these catalogues. Note that SPIRE observations of GOODS-N/S cover much larger effective areas than those from PACS (see Table 1). However, here we restricted our study to regions with both PACS and SPIRE data.

Our MIPS-PACS-SPIRE catalogues were cross-matched with our multi-wavelength catalogues (Sect. 2.4), using their MIPS/IRAC positions and a matching radius of 1′′.

2.2. VLA observations

To study the radio properties of our galaxies, we use deep publicly available 1.4 GHz VLA observations. In the COSMOS field, we use observations from the VLA-COSMOS Project⁴ (Schinnerer et al. 2004, 2007, 2010). This project has imaged the entire COSMOS field at 1.4 GHz to a mean rms noise of s10 (15) μJy beam^-1 across its central 1(2)deg² with a resolution of (FWHM of the synthesised beam). In the GOODS-N field, we use deep 1.4 GHz VLA observations⁵ described in Morrison et al. (2010). This map has a synthesised beam size of and a mean rms noise of s3.9μJy beam^-1 near its centre, rising to s8μJy beam^-1 at a distance of 15′ from the field centre. Finally, in the GOODS-S and ECDFS fields, we use deep 1.4 GHz VLA observations⁶ presented in Miller et al. (2008, 2013). This map has a synthesised beam size of and an mean rms noise of s8μJy beam^-1. Table 2 summarises the properties of all these observations.

2.3. GMRT observations

GMRT 610 MHz observations of the GOODS-N and ECDFS fields were obtained in 2009–10 as part of the GOODS-50 project (PI: Ivison). The ECDFS observations were described by Thomson et al. (2014); those of GOODS-N covered a single pointing rather than the six used to cover ECDFS, but in most other respects were identical. The GMRT observations of ECDFS and GOODS-N reach a mean rms noise of s15μJy beam^-1. The synthesised beam measured and in the ECDFS and GOODS-N maps, respectively. Table 2 summarises the properties of all these observations.

Note that there are no GMRT observations of the COSMOS field. Therefore, constraints on the radio spectral index of star-forming galaxies (Sect. 4.1) are obtained from a somehow more limited galaxy sample, covering a sky area of s0.3 deg². These limitations could have restrained our ability to constrain the radio spectral index of faint star-forming galaxies. However, we demonstrate in Sect. 3.3 that with our current GMRT observations, we are able to estimate the radio spectral index of galaxies with M_∗ > 10¹⁰ M_⊙, Δlog (SSFR)_MS > −0.3 and 0 <z< 2.3. These constraints are sufficient for the purpose of our study. In addition, because we do not expect significant cosmic variance, we can apply these constraints obtained on a sky area of s0.3 deg² to our entire galaxy sample.

2.4. Multi-wavelength catalogues

The large wealth of multi-wavelength data used in our study is described in detail in Wuyts et al. (2011a,b; see also Magnelli et al. 2014). Here we only summarise the properties relevant for our study.

In the COSMOS field we used 36 medium and broad-band observations covering the optical to near-infrared SEDs of galaxies (Ilbert et al. 2009; Gabasch et al. 2008). We restricted these catalogues to i< 25 and to sources not flagged as problematic in the catalogue of Ilbert et al. (2009). In the GOODS-S field, we used the K_s< 24.3 (5σ) FIREWORKS catalogue, providing photometry in 16 bands from U to IRAC wavelengths. In the GOODS-N field, we used the K_s< 24.3 (3σ) catalogue created as part of the PEP survey⁷, which provides photometry in 16 bands from GALEX to IRAC wavelengths. Finally, in the ECDFS field, we used the R< 25.3 (5σ) catalogue described in Cardamone et al. (2010). This multi-wavelength catalogue provides photometry in 18 bands from U to IRAC wavelengths. In the following, ECDFS only corresponds to the outskirts of the original ECDFS region, i.e. when we refer to the ECDFS field, we implicitly exclude the central GOODS-S region.

Spectroscopic redshifts for our galaxies were taken from a combination of various studies (Cohen et al. 2000; Cristiani et al. 2000; Croom et al. 2001; Cimatti et al. 2002, 2008; Wirth et al. 2004; Cowie et al. 2004; Le Fèvre et al. 2004; Szokoly et al. 2004; van der Wel et al. 2004; Mignoli et al. 2005; Vanzella et al. 2006, 2008; Reddy et al. 2006; Barger et al. 2008; Kriek et al. 2008; Lilly et al. 2009; Treister et al. 2009; Balestra et al. 2010). For sources without a spectroscopic redshift, we used photometric redshift computed using EAZY (Brammer et al. 2008) exploiting all the available optical/near-infrared data. The quality of these photometric redshifts was assessed via comparison with spectroscopically confirmed galaxies. The median and median absolute deviation of Δz/ (1 + z) are (−0.001; 0.013) at z< 1.5, (−0.007; 0.066) at z > 1.5 and (−0.001; 0.014) at 0 <z< 2.3 (Fig. 1).

2.4.1. Stellar masses

The stellar masses of our galaxies were estimated by fitting all m data to Bruzual & Charlot (2003) templates using FAST (Fitting and Assessment of Synthetic Templates; Kriek et al. 2009). The rest-frame template error function of Brammer et al. (2008) was used to down-weight data points with m. For those stellar mass estimates we adopted a Chabrier (2003) IMF. Full details on those estimates and their limitations are given in Wuyts et al. (2011a,b).

2.4.2. Star-formation rates

To estimate the SFRs of our galaxies we used the cross-calibrated ladder of SFR indicators established in Wuyts et al. (2011a). This uses the best indicator available for each galaxy and establishes a consistent scale across all of them. For galaxies only detected in the rest-frame UV (i.e. those without a mid- or far-infrared detection), SFRs were estimated from the best fits obtained with FAST. For galaxies with detections both in the rest-frame UV and the mid-/far-infrared, SFRs were estimated by combining the unobscured and re-emitted emission from young stars. This was done following Kennicutt (1998) and adopting a Chabrier (2003) IMF: (1)where L₂₈₀₀ ≡ νL_ν(2800 Å) was computed with FAST from the best-fitting SED and the rest-frame infrared luminosity L_IR ≡ L [ 8−1000 μm ] is derived from the mid-/far-infrared observations. For galaxies with FIR detections, L_IR was inferred by fitting their FIR flux densities (i.e. those measured using PACS and SPIRE) with the SED template library of Dale & Helou (2002, DH), leaving the normalisation of each SED template as a free parameter⁸. The infrared luminosities of galaxies with only a mid-infrared detection were derived by scaling the SED template of MS galaxies (Elbaz et al. 2011) to their 24 μm flux densities. Magnelli et al. (2014) have shown that this specific SED template provides accurate 24μm-to-L_IR conversion factors for such galaxies.

2.4.3. Active galactic nuclei contamination

Active galactic nuclei can affect the observed FRC of star-forming galaxies. Such AGNs must thus be excluded from our sample. To test for the presence of AGNs, we used the deepest available Chandra and XMM-Newton X-ray observations, identifying AGNs as galaxies with L_X [ 0.5−8.0 keV ] ≥ 3 × 10⁴² erg s^-1 and L_X [ 2.0−10.0 keV ] ≥ 3 × 10⁴² erg s^-1, respectively (Bauer et al. 2004). In the GOODS-N and -S fields, X-ray observations were taken from the Chandra 2-Ms catalogues of Alexander et al. (2003) and Luo et al. (2008), respectively. In the COSMOS field, we used the XMM-Newton catalogue of Cappelluti et al. (2009). Finally, for the ECDFS, we used the 250 ks Chandra observations, which flank the 2-Ms CDFS observations (Lehmer et al. 2005). All these X-ray-selected AGNs have been removed from our sample. Note that radio-loud AGNs without X-ray detection are excluded in our radio stacking procedure, through use of median stacking (Sect. 3.2).

2.5. Final sample

Our four multi-wavelength catalogues are not homogeneously selected and are not uniform in depth, which naturally translates into different completeness limits in the SFR–M_∗ plane. These issues have been discussed and studied in Magnelli et al. (2014). They found that in the GOODS-S, GOODS-N and COSMOS fields our multi-wavelength catalogues are complete for star-forming galaxies with M_∗ > 10¹⁰ M_⊙ up to z s 2. Because the ECDFS multi-wavelength catalogue is based on deeper optical/near-infrared observations than those of the COSMOS field, we conclude that this catalogue also provides us with a complete sample of star-forming galaxies with M_∗ > 10¹⁰ M_⊙ up to z s 2. In the rest of the paper, we restrict our results and discussion to galaxies with M_∗ > 10¹⁰ M_⊙.

Our final sample contains 8846, 4753, 66070, and 254749 sources in the GOODS-N, GOODS-S⁹, ECDFS and COSMOS fields, respectively. Of these sources, 29%, 26%, 1% and 3% have a spectroscopic redshift, while the rest have photometric redshift estimates. Because we are studying the FRC, the SFR–M_∗–z bins that enter our analysis (see Sect. 3) are generally dominated by sources that have individual mid-infrared (and for part of them far-infrared) detections. In GOODS-N, GOODS-S, ECDFS and COSMOS, 19%, 28%, 9% and 12% of the galaxies have mid- or far-infrared detections, respectively. Among those sources, 60%, 45%, 5% and 11% have a spectroscopic redshift.

Figure 2 shows the number density of sources in the SFR–M_∗ plane. Over a broad range of stellar masses, star-forming galaxies (i.e. excluding massive and passive galaxies situated in the lower right part of the SFR–M_∗ plane) follow a clear SFR–M_∗ correlation. This correlation is known as the MS of star formation (Noeske et al. 2007). In the rest of the paper, we parametrise this MS using second-order polynomial functions as derived by Magnelli et al. (2014, see their Table 2). These functions are presented in Fig. 2. Comparisons between this parametrisation and those from the literature are presented and discussed in Magnelli et al. (2014). Briefly, the MS observed in our sample is consistent with the literature at M_∗ > 10¹⁰ M_⊙, i.e. within the stellar mass range of interest for our study.

3. Data analysis

The aim of this paper is to study the evolution of the FRC and radio spectral index with redshift and with respect to the position of galaxies in the SFR–M_∗ plane. Although we could base this analysis on galaxies individually detected at FIR and/or radio wavelengths, such an approach would be subject to strong limitations, mainly due to complex selection functions (see e.g. Sargent et al. 2010). Instead, we adopted a different approach based on a careful FIR and radio stacking analysis of a stellar-mass-selected sample. This allows us to probe the properties of the FRC and radio spectral index, delving well below the detection limits of current FIR and radio observations. Of course, the use of this stacking analysis has the obvious drawback that one can only study the mean properties of the FRC within the SFR–M_∗ plane, while outliers are completely missed out. This limitation has to be taken into account while discussing our results. In addition, our stacking analysis has to be performed with great care, especially for the FIR observations where large beam sizes might lead to significant flux biases if the stacked samples are strongly clustered.

3.1. Determination of far-infrared properties through stacking

To estimate the FIR properties (i.e. L_FIR and T_dust) of a given galaxy population, we stacked their Herschel observations. The stacking method adopted here is similar to that used in Magnelli et al. (2014). In the following we only summarise the key steps of this method, while for a full description we refer the reader to Magnelli et al. (2014).

For each galaxy population (i.e. for each SFR–M_∗–z bin) and for each Herschel band, we stacked the residual image (original maps from which we removed all 3σ detections) at the positions of undetected sources (i.e. sources with S_Herschel< 3σ). The stacked stamp of each galaxy was weighted with the inverse of the square of the error map. The flux densities of the final stacked images were measured by fitting with the appropriate PSF. Uncertainties on these flux densities were computed by means of a bootstrap analysis. The mean flux density (S_bin) of the corresponding galaxy population was then computed by combining the fluxes of undetected and detected sources: (2)where S_stack is the stacked flux density of the m undetected sources, and S_i is the flux density of the ith detected source (out of a total of n). We note that consistent results are obtained if we repeat the stacking analysis using the original PACS/SPIRE maps and combining all sources in a given SFR–M_∗–z bin, regardless of whether they are individually detected or not. We have verified that consistent results are obtained via a median stacking method, rather than the mean stacking described above.

To verify that the clustering properties of our stacked samples have no significant effect on our stacked FIR flux densities, we used simulations from Magnelli et al. (2014). Briefly, simulated Herschel flux densities of all sources in our final sample were estimated using the MS template of Elbaz et al. (2011), given their observed redshifts and SFRs. Simulated Herschel maps with real clustering properties were then produced using the observed positions and simulated Herschel flux densities of each sources of our fields. Whenever we stacked a given galaxy population on the real Herschel images, we also stacked at the same positions the simulated images and thus obtained a simulated stacked flux density (). Then we compared the with the expected mean flux density of this simulated population, i.e. . If ABS(()/, then the real stacked flux densities were identified as being potentially affected by clustering. This 0.5 value was empirically defined as being the threshold above which the effect of clustering would not be captured within the flux uncertainties of our typical S/Ns4 stacked flux densities. The largest clustering effects are observed at low flux densities and in the SPIRE 500 μm band, as expected. More details on these simulations can be found in Sect. 3.2.2 of Magnelli et al. (2014).

From the mean Herschel flux densities of each galaxy population we inferred their rest-frame FIR luminosities (i.e. L_FIR, where L_FIR is the integrated luminosity between 42μm and 122μm) and dust temperatures (i.e. T_dust) by fitting the available FIR photometry using the DH SED template library and a standard χ² minimisation method. From the integration of the best-fitting DH SED template, we infer L_FIR to within s0.1 dex, even in cases with only one FIR detection, because Herschel observations probe the peak of the FIR emission of galaxies (Elbaz et al. 2011; Nordon et al. 2012). Naturally, we compared these L_FIR estimates to the obscured SFR (i.e. SFR_IR) expected from our ladder of SFR indicators. We reject SFR–M_∗–z bins in which these two independent SFR_IR estimates are not consistent within 0.3 dex. Such discrepancies are only observed in few SFR–M_∗–z bins with stacked FIR flux densities with low significance, S/Ns3.

From the best-fitting DH SED template we also estimated T_dust using the pairing between dust temperature and DH templates established in Magnelli et al. (2014). The reliability of these dust temperature estimates depends on the number of FIR data points available and on whether those data points encompass the peak of the FIR emission. As in Magnelli et al. (2014), we considered as reliable only the dust temperatures inferred from at least three FIR data points encompassing the peak of the FIR emission and with .

3.2. Determination of the radio properties through stacking

To estimate the radio properties (i.e. S_{1.4 GHz} and S_{610 MHz}) of a given galaxy population, we stacked their VLA 1.4 GHz and GMRT 610 MHz observations. This stacking analysis is very similar to that employed for the Herschel observations. However, the Herschel data are very homogeneous among fields, except for the noise level, whereas the radio data differ between fields in properties such as beam shape. Hence, one cannot stack together sources regardless of their position on the sky. Instead, one can stack together sources of the same field, measure their stacked flux densities and associated errors, and finally combine information from different fields using a weighted mean. For a given galaxy population (i.e. a given SFR–M_∗–z bin), for each field and each radio band (1.4 GHz and 610 MHz), we thus proceeded as follows. We stacked all sources (regardless of whether they are detected or not) situated in a given field using the original VLA (or GMRT) images. We measured the radio stacked flux density of this galaxy population in this field by fitting their median radio stacked stamp with a 2D Gaussian function. Use of the median avoids the biasing influence of moderately radio-loud AGNs (e.g. Del Moro et al. 2013) that have not yet been removed from the source list via X-ray emission. Because the S/N of our typical radio stacked stamp was poor, in these fits we only left the normalisation of the 2D Gaussian function as a free parameter (i.e. S_Gaussian). The position, minor and major axes (i.e. [a,b]), and position angle of this 2D Gaussian function were fixed to the values found when fitting the high S/N radio stacked stamp of all galaxies within the current field, redshift bin and with M_∗ > 10¹⁰ M_⊙ and Δlog (SSFR)_MS > −0.3. The radio stacked flux density of this SFR–M_∗–z bin in this field (i.e. ) was then given by (3)where a_norm and b_norm are the minor and major axes of the radio beam in the original VLA (or GMRT) images. For the VLA observations, where the spatial resolution is relatively high and thus beam smearing by astrometric uncertainties in the stacked samples can be significant, S_radio were s1.9 times higher than S_Gaussian. Despite astrometric uncertainties, galaxy sizes marginally resolved at the s1.5″ resolution of the VLA observations might also explain part of the difference between S_radio and S_Gaussian. The appropriate uncertainty on this radio flux density was obtained from a bootstrap analysis. Finally, we combined the radio stacked flux densities of a given SFR–M_∗–z bin from different fields (i.e. GOODS-N, GOODS-S, ECDFS and COSMOS) using a weighted mean.

From their stacked 1.4 GHz flux densities, we derived the rest-frame 1.4 GHz luminosity (i.e. L_{1.4 GHz}) for our SFR–M_∗–z bins. For that, we k-corrected their stacked 1.4 GHz flux density assuming a radio spectral index of α = 0.8 (e.g. Condon 1992; Ibar et al. 2009, 2010). This value is perfectly in line with the radio spectral index observed here (i.e. ) over a broad range of redshift and Δlog (SSFR)_MS (see Sect. 4.1).

3.3. The SFR-M_∗-z parameter space

Before looking at the evolution of the FRC and radio spectral index in the SFR–M_∗–z parameter space, we need to ensure that our ability to make accurate q_FIR and measurements does not introduce significant incompleteness in any particular regions of the SFR−M_∗ plane. The left panels of Figs. 3 and 4 present the regions of the SFR–M_∗–z parameter space with accurate and q_FIR estimates from our stacking analysis, respectively. The sampling of the SFR–M_∗ plane is made with larger SFR–M_∗ bins for than for q_FIR. This is due to the fact that the GMRT observations covered a sky area of only 0.3 deg² (Sect. 2.3). This limited dataset forces us to enlarge the size of our SFR–M_∗ bins in order to increase the number of stacked sources per bin and thus improve the noise in our stacked stamps (σ_stack∝). Our estimates are considered as accurate only if S_{1.4 GHz}/N_{1.4 GHz} > 4 and S_{610 MHz}/N_{610 MHz} > 4. Our q_FIR estimates are considered as accurate only if S_{1.4 GHz}/N_{1.4 GHz} > 3 and L_FIR is reliable (see Sect. 3.1). In each of our redshift bins, our stacking analysis allows us to obtain accurate and q_FIR estimates for almost all MS and above-MS galaxies with M_∗ > 10¹⁰ M_⊙.

The right panels of Figs. 3 and 4 illustrate our ability to study the variations of and q_FIR for galaxies with M_∗ > 10¹⁰ M_⊙, respectively. In these figures we show the fraction of SFR−M_∗ bins with reliable or q_FIR estimates as a function of their Δlog(SSFR)_MS. In the rest of the paper, we consider that or q_FIR in a given Δlog (SSFR)_MS bin is fully constrained only if the completeness in this bin is ≳80%. In each redshift bin, our stacking analysis allows us to fully constrain and q_FIR in galaxies with Δlog (SSFR)_MS > −0.3.

From this analysis, we conclude that our stacking analysis provides us with a complete view on the evolution of and q_FIR up to zs2 in star-forming galaxies with M_∗ > 10¹⁰ M_⊙ and Δlog (SSFR)_MS > −0.3.

4. Results

4.1. The radio spectral index,

The radio spectral index of each SFR–M_∗–z bin was inferred using their VLA and GMRT stacked flux densities, i.e. S_{1.4 GHz} and S_{610 MHz}, respectively. Assuming that the radio spectrum follows a power law form, S_ν ∝ ν^−α, the radio spectral index is given by, (4)S_{1.4 GHz} and S_{610 MHz} are observed flux densities. Therefore, at different redshift, correspond to radio spectral indices at different rest-frame frequencies. This has to be taken into account when interpreting our results (see Sect. 5.1.2). In the following, we discuss only accurate estimates. We remind the reader that constraints on the radio spectral index of star-forming galaxies are obtained from a somewhat limited galaxy sample, covering a sky area of s0.3 deg², as GMRT observations of the COSMOS fields are not available.

The Fig. 5 presents the evolution of as a function of Δlog (SSFR)_MS up to zs2. We find no strong correlation between and Δlog (SSFR)_MS in any redshift bin (i.e. | ρ_s | ≲ 0.6) and the null hypothesis of uncorrelated data cannot be rejected with high significance (i.e. Sig.>5%). In addition, the median does not deviate significantly from the canonical value, 0.8 (e.g. Condon 1992), in any of our redshift bins. From these findings we conclude that most star-forming galaxies with M_∗ > 10¹⁰ M_⊙ and across 0 <z< 2.3 have on average a radio spectral index consistent with 0.8. However, due to the relatively low number of data points and to the large dispersion on , we cannot confidently rule out the presence of a negative, but weak, – Δlog (SSFR)_MS correlation (ρ_s< 0). Such a trend would echo (but not match) the results of Condon et al. (1991, see also Clemens et al. 2008 who found that the most extreme local starbursts have flatter radio spectra (i.e. s0.5) than local normal star-forming galaxies (i.e. s0.8). Condon et al. (1991) and Clemens et al. (2008) attribute these flat radio spectra to free-free absorption in dense nuclear starbursts.

The Fig. 6 shows that our results agree well with the range of radio spectral index observed in a large sample of sub-mJy radio galaxies (Ibar et al. 2009) and in a population of zs2 SMGs (Ibar et al. 2010, see also Thomson et al. 2014). These agreements are re-assuring because results from Ibar et al. (2009, 2010) and Thomson et al. (2014) were based on galaxies individually detected in the VLA and GMRT images.

Our results are also in line with those of Bourne et al. (2011) and Ivison et al. (2010a). Both studies are based on a stacking analysis, but Ivison et al. (2010a) applied it to a sample of 24 μm-selected (i.e. SFR-selected) galaxies while Bourne et al. (2011) applied it to a stellar-mass-selected galaxy sample (including both star-forming and quiescent galaxies). The consistencies observed between studies with different selection functions re-inforce our conclusion that does not significantly evolve across 0 <z< 2.3. We note that results from our study and those from Ivison et al. (2010a) and Bourne et al. (2011) are, however, not strictly independent since they are based on the same VLA and similar GMRT observations of ECDFS. Nevertheless, we believe that the consistencies found here are noteworthy because those studies differ in many other aspects: our stacking analysis includes observations from the GOODS-N field; our GMRT observations of ECDFS are deeper; our stellar-mass-selected sample is built using different optical-to-near-IR multi-wavelength catalogues, applying different methods to infer photometric redshifts and stellar masses.

The absence of significant evolution in with z and Δlog (SSFR)_MS is an important result for our forthcoming study of the FRC. Indeed, to k-correct our stacked 1.4 GHz flux densities into rest-frame 1.4 GHz radio luminosities, we assumed that the radio spectral index of all galaxies across 0 <z< 2.3 was equal to 0.8 (see Sect. 3.2). Significant deviations from this canonical value would have introduced artificial evolution of the FRC in the SFR–M_∗–z parameter space (Δq_FIR = −Δα × log (1 + z)).

The physical implications of the absence of significant evolution of with z and Δlog(SSFR)_MS are discussed in Sect. 5.

4.2. The FIR/radio correlation

Using the rest-frame FIR and 1.4 GHz luminosities estimated from our stacking analysis, we study the evolution of the FRC in the SFR–M_∗–z parameter space. For that, we use the parametrisation of the FRC given in Helou et al. (1988; see also Yun et al. 2001), (5)where L_FIR is the integrated FIR luminosity from rest-frame 42 to 122μm and L_{1.4 GHz} is the rest-frame 1.4 GHz radio luminosity density. Radio k-corrections are inferred assuming S_ν ∝ ν^−α and the canonical radio spectral index of α = 0.8 (see Sect. 4.1; Condon 1992). In the recent literature, alternative parametrisations of the FRC have been proposed. In particular, some studies have used the infrared luminosity from rest-frame 8 to 1000μm (L_IR) instead of L_FIR (e.g. Ivison et al. 2010a,b; Sargent et al. 2010; Bourne et al. 2011). Such different definitions have no significant impact because there is a tight relation between L_FIR and L_IR of star-forming galaxies. Over all our SFR–M_∗–z bins, we indeed found . Thus, to compare q_FIR with q_IR, we simply corrected these estimates following q_FIR = q_IR−log (1.91).

The Fig. 7 shows the evolution of q_FIR in the SFR−M_∗ plane in different redshift bins. In this figure, we do not identify any significant and systematic evolution of q_FIR in the SFR–M_∗ plane. In contrast, we notice a clear and systematic decrease in q_FIR with redshift.

In Fig. 8 we investigate the existence of more subtle evolution of q_FIR within the SFR–M_∗ plane by plotting the variation of q_FIR as a function of Δlog (SSFR)_MS. In all our redshift bins, there is a weak (0.1 ≲ | ρ_s | ≲ 0.4) correlation between q_FIR and Δlog (SSFR)_MS, though the null hypothesis of uncorrelated data can only be rejected with high significance (i.e. Sig.<5%) in two of these redshift bins (i.e. 0.5 <z< 0.8 and 0.8 <z< 1.2). In addition, at high Δlog (SSFR)_MS, the dispersion on q_FIR seems to increase.

To study further the statistical significance of a weak q_FIR– Δlog(SSFR)_MS correlation, we fit this relation with a constant and with a linear function using a Monte-Carlo approach taking into account errors both on q_FIR and Δlog(SSFR)_MS. Data points used in these fits are restricted to those situated in regions of parameter space not affected by incompleteness. For each of our 1000 Monte-Carlo realisations, we adopt new values of q_FIR and Δlog(SSFR)_MS selected into a Gaussian distribution centred on their original values and with a dispersion given by their measurement errors. To ensure that our fits are not dominated by few data points, for each Monte-Carlo realisation we resample the observed dataset, keeping its original size but randomly selecting its data points with replacement. We then fit each Monte-Carlo realisation with a constant and with a linear function. Finally, we study the mean value and dispersion of each fitting parameter across our 1000 Monte-Carlo realisations.

In all but the 0.5 <z< 0.8 redshift bin, the constant model has reduced χ² values lower than those of the linear model. In addition, in these redshift bins, the linear model has slopes (i.e. Δ [ q_FIR ]/Δ[Δlog(SSFR)_MS ]) consistent, within 1σ, with zero. This suggests that there is no significant q_FIR– Δlog(SSFR)_MS correlation in these redshift bins. In contrast, at 0.5 <z< 0.8, the linear model is statistically slightly better than the constant model (i.e. and ) and its slope is different to zero at the s3σ level (i.e. 0.22 ± 0.07). This suggests a positive but weak q_FIR– Δlog(SSFR)_MS correlation¹⁰ at 0.5 <z< 0.8. We note that it is also in this redshift bin that our measurement errors on q_FIR are the lowest. This could explain why an intrinsically weak q_FIR– Δlog(SSFR)_MS correlation can only be statistically significant at 0.5 <z< 0.8.

This Monte-Carlo approach demonstrates that the existence of a q_FIR– Δlog(SSFR)_MS correlation is yet difficult to assess from our dataset. If there exists a q_FIR– Δlog(SSFR)_MS correlation, it is intrinsically weak. Unfortunately, our observations are not good enough to firmly reveal or rule out such weak correlation in all our redshift bins. Thus, we conservatively conclude that there is no significant q_FIR– Δlog(SSFR)_MS correlation across 0 <z< 2.3, though the presence of a weak, positive trend, as observed in one of our redshift bin (0.5 <z< 0.8), cannot be firmly ruled out using our dataset. Note that debates on the existence of a weak q_FIR– Δlog(SSFR)_MS correlation also prevail in the local Universe. Condon et al. (1991) found that the most extreme local starbursts (i.e. Δlog(SSFR)_MSs1) have higher q_FIR and larger dispersions of q_FIR than normal star-forming galaxies (i.e. Δlog(SSFR)_MSs0). In contrast, Helou et al. (1985) and Yun et al. (2001) do not report any statistically significant increase in q_FIR in extreme starbursts, though they also found a larger dispersion in q_FIR for this population.

While q_FIR does not significantly evolve with Δlog(SSFR)_MS, its median value clearly decreases smoothly with redshift. To study this trend we plot in Fig. 9 the redshift evolution of the median and interquartile range of q_FIR, observed in regions of parameter space not affected by incompleteness. We find a statistically significant (Sig.<1%) redshift evolution of q_FIR. This evolution can be parametrised using, (6)where the z = 0 value of this function agrees perfectly with local observations, i.e. q_FIR(zs0)≈2.34 ± 0.26 (Yun et al. 2001). If we parametrise this redshift evolution separately for normal (Δlog(SSFR)_MS< 0.75) and starbursting galaxies (Δlog(SSFR)_MS > 0.75)¹¹, we end up with similar solutions, i.e. q_FIR(z) = (2.42 ± 0.08) × (1 + z)^{−0.16 ± 0.05} and q_FIR(z) = (2.44 ± 0.09) × (1 + z)^{−0.13 ± 0.09}, respectively. We note that while the evolution of q_FIR with redshift is statistically significant, it is moderate. Indeed, at zs2, the median value of q_FIR is still within the 1σ dispersion of local observations (Yun et al. 2001).

This redshift evolution of q_FIR could not be artificially introduced by a redshift evolution of the radio spectral index. Indeed, to create such evolution, the radio spectral index would need to change from 0.8 to 0.2 between z = 0 and z = 2, respectively (Δq_FIR = −Δα × log (1 + z)). Such extreme redshift evolution of the radio spectral index is not observed in our sample (see Fig. 6 and Sect. 4.1).

In Fig. 9 we compare our findings with those from the literature. In the past decade many papers have discussed this topic. Therefore, instead of presenting an exhaustive comparison, we compare our findings with the three papers which are, we believe, the most relevant, i.e. those based on relatively complete and well understood samples and/or valuable FIR/radio datasets. Firstly, we compare our results with those of Bourne et al. (2011). In that paper, the authors overcame the selection biases of radio- and/or infrared-selected samples by using a stellar-mass-selected sample in the ECDFS. The infrared and radio properties of their galaxies at a given stellar mass and redshift were then inferred by stacking Spitzer (24, 70 and 160μm) and VLA observations. This approach is thus very similar to that employed here. Secondly, we compare our results with those of Sargent et al. (2010). This paper used a large sample of infrared- and radio-selected galaxies in the COSMOS field. To overcome selection biases, Sargent et al. (2010) applied a careful survival analysis to the Spitzer (24 and 70μm) and VLA catalogues. Finally, we compare our results with those of Ivison et al. (2010b). In this paper, the authors used Herschel observations of the GOODS-N field to constrain the redshift evolution of q_FIR on a relative small sample of L_IR-matched galaxies.

Results from the literature are in broad agreement with our conclusions. Data from Bourne et al. (2011) exhibit a consistent q_FIR–z correlation, though with a slight overall offset of their q_FIR values (Δq_FIRs0.15). Bourne et al. (2011) also noticed that their q_FIR values were systematically higher than that expected from local observations. They attributed this offset to difference in the assumptions made to infer L_IR from Spitzer observations. Most likely, the overall offset observed here between their and our measurements has the same origin.

Results from Ivison et al. (2010b) also exhibit a clear decrease in q_FIR with redshift. However, here as well there exists some disagreement with respect to our measurements: the decrease in q_FIR starts at z > 0.7, and at 0.7 <z< 1.4, their measurements of q_FIR are significantly higher. These discrepancies could not be explained by difference in the assumptions made to derive L_FIR, because both studies used FIR observations from Herschel. Instead, these discrepancies might be explained by cosmic variance and/or by some differences in our sample selection. For example, at zs0.5, the measurement of Ivison et al. (2010b) relied on only 16 galaxies and they argued that it might need to be discounted. In addition, the sample of Ivison et al. (2010b) is L_IR-selected and not SFR−M_∗-selected, and still includes X-ray sources.

The results of Sargent et al. (2010) are perfectly in line with our findings. They found a clear decrease in q_FIR with redshift, and the amplitude and normalisation of their q_FIR–z correlation is consistent with our measurements.

From all these comparisons, we conclude that results from the literature are also broadly consistent with a decrease in q_FIR with redshift. Note, however, that these studies have mostly chosen to favour a non-evolving scenario for the FRC. This moderate redshift evolution of q_FIR, which remains within its local 1σ dispersion even at zs2, was not deemed sufficiently significant in most of these studies, which relied on relatively small samples with sparse FIR and radio spectral coverage. Our sample is sufficiently large and well controlled (i.e. complete for star-forming galaxies with M_∗ > 10¹⁰ M_⊙) with excellent FIR and radio spectral coverage to conclude that q_FIR evolves with redshift as q_FIR ∝ (1 + z)^{−0.12 ± 0.04}. Possible physical explanations and implications of this redshift evolution are discussed in Sect. 5.

In Fig. 10, we consider the reality of a correlation between q_FIR and T_dust. In all our redshift bins, there is a weak, positive q_FIR–T_dust correlation; in four bins we can reject the null hypothesis of uncorrelated data with high significance (Sig.<5%). However, using our Monte-Carlo approach to fit this relation, we find that only at 0.5 <z< 0.8 the linear model has better reduced χ² values than the constant model and has a slope different than zero at the s3σ level (i.e. Δ [ q_FIR ]/Δ [ T_dust ] = 0.023 ± 0.008). For the other redshift bins, the linear and constant models are statistically undistinguishable and the slopes of the linear model are consistent, within 1σ, with zero. Thus, we conclude that the existence of a q_FIR–T_dust correlation is statistically meaningful in only one of our redshift bin. Such q_FIR–T_dust correlation could be related to the weak, positive q_FIR– Δlog(SSFR)_MS correlation observed in the same redshift bin, because T_dust is known to be positively correlated with Δlog(SSFR)_MS (Magnelli et al. 2014). Note, however, that our estimates of T_dust might be affected/contaminated by AGN emission, complicating the interpretation on the existence of a q_FIR–T_dust correlation.

5. Discussion

5.1. Evolution of the radio spectra

Our results indicate that the radio spectral index, , does not significantly evolve with redshift, does not correlate significantly with Δlog(SSFR)_MS and is consistent everywhere with its canonical value of 0.8 (see Figs. 5 and 6). Those results are valid for all star-forming galaxies with M_∗ > 10¹⁰ M_⊙, Δlog(SSFR)_MS > −0.3 and 0 <z< 2.3.

5.1.1. AGN contamination?

Some may find the absence of significant redshift evolution in α surprising, perhaps anticipating the rapid increase in the AGN population at z > 1 (e.g. Hasinger et al. 2005; Wall et al. 2005) and their influence on the observed radio spectral index of their host galaxies. At low redshift and therefore low rest-frame frequencies, AGNs can exhibit flat radio spectra (α< 0.5; e.g. Murphy 2013) while at higher redshift and thus higher rest-frame frequencies they can exhibit steep radio spectra (α > 1.0; e.g. Huynh et al. 2007). With the increasing number of AGNs at z > 1, one might then expect variations in the radio spectral index inferred from our statistical sample. However, since our goal is to explore the properties of star-forming galaxies, and anticipating the possible impact of AGNs, we excluded the brightest X-ray AGNs from our sample and used a median radio stacking method. This should have minimised any AGN contamination. Interestingly, if we stack our bright X-ray AGNs, we find radio spectra indices in the range [0.7−1.0]¹². This suggests that the radio spectral index of the bulk of the X-rays AGN population is dominated by emission from their host galaxies. Note, however, that the AGN contribution to the radio emission of a galaxy could also have a radio spectral index of s0.75. Indeed, while the radio emission of the compact core of AGNs is supposed to have a flat spectrum at low frequencies and a steep spectrum at high frequencies, their extended radio emission associated with lobes is supposed to have a power-law spectrum with αs0.75 (Jackson & Wall 1999).

5.1.2. Nature of the radio spectra

Our VLA and GMRT observations probe different rest-frame radio frequencies at different redshifts. Consequently, even if all star-forming galaxies have the same intrinsic radio spectra, one could still expect to observe redshift evolution of in the presence of a curved spectrum. Such redshift evolution of is not observed in our sample, suggesting that in the range of rest-frame frequencies probed here, 610 MHz <ν_rest< 4.2 GHz, the radio spectra are well described by a power-law function with S_ν ∝ ν^-0.8. The radio spectra of high-redshift star-forming galaxies seem to be dominated by non-thermal optically thin synchrotron emission with the same properties as that observed in the local Universe; thermal free-free emission with relatively flat radio spectra (i.e. αs0.1) does not dominate their rest-GHz radio spectra. Unfortunately, these conclusions might be compromised if the constant value is caused by a conspiracy, involving a curved radio spectra and some intrinsic evolution with redshift.

In the local Universe, starbursts/ULIRGs have significantly flatter radio spectra (αs0.5) than normal star-forming galaxies (αs0.8; Condon et al. 1991; Clemens et al. 2008). While our sample might exhibit some flattening of the radio spectral index when moving from the MS regime (Δlog(SSFR)_MSs0) to the starburst regime (Δlog(SSFR)_MSs1), these – Δlog(SSFR)_MS correlations are weak and statistically insignificant (Sig.>5%). The observation in high-redshift starbursts of steeper radio spectra than in the local Universe echoes results for samples of SMGs, which also have αs0.8 (Ibar et al. 2010; Thomson et al. 2014). Indeed, luminous SMGs are situated well above the MS (Magnelli et al. 2012a), are believed to be strong starbursts and to be the high-redshift counterparts of local ULIRGs. The fact that high-redshift starbursts have steeper radio spectra suggests that they may not have the same ISM conditions as their local counterparts, as already hinted by numerous other studies looking at their size (e.g. Tacconi et al. 2006; Biggs & Ivison 2008; Farrah et al. 2008) and physical properties (e.g. Ivison et al. 2010c; Farrah et al. 2008). Alternatively, it could indicate that the flatter radio spectra of local ULIRGs are due to low-frequency free-free absorption, as advocated by Condon et al. (1991), a less relevant effect at the higher rest-frame frequencies probed at high redshift, though at zs2 our data are not probing above s4 GHz.

5.1.3. Expectations from theory

Using a single-zone model of CR injection, cooling and escape, Lacki et al. (2010, hereafter L10) and Lacki & Thompson (2010, hereafter LT10) explore the radio spectra of normal and starburst galaxies and their evolution with redshift. At zs0, their model predicts a relatively flat radio spectrum at rest-GHz frequencies for compact starbursts (αs0.5) and a steeper spectrum for normal star-forming galaxies (αs0.8), though this difference decreases at higher rest-frame frequencies (≫10 GHz) where both asymptote to αs1.1. The flatter radio spectrum of starbursts is attributed to more efficient Bremsstrahlung and ionisation cooling of CR electrons and positrons in their dense ISM. Thus, the L10 model predicts the observed flatter radio spectra of local ULIRGs but contradicts the hypothesis of Condon et al. (1991) about its origin via free-free absorption.

Our observations do not reveal a significant decrease in the radio spectral index as we move from normal star-forming galaxies (Δlog(SSFR)_MSs0) to starbursts (Δlog(SSFR)_MSs1). Instead, high-redshift starbursts seem to have roughly the same radio spectral index as normal star-forming galaxies. As already mentioned, this result was hinted by the observation of steeper radio spectra and more extended star-forming regions in high-redshift SMGs than in local ULIRGs. Therefore, LT10 suggested that high-redshift SMGs have larger CR scale heights (hs1 kpc), consistent with that observed in local normal star-forming galaxies. This large CR scale height leads to less efficient Bremsstrahlung and ionisation cooling of CRs and thus steeper radio spectra. Our results suggest that the bulk of the high-redshift starburst population, and not only SMGs, may have large CR scale heights, labelled puffy starbursts in LT10.

LT10 also predict evolution of the radio spectral index of galaxies with redshift. Combining intrinsic evolution of the radio spectra with k-correction considerations, they find that of normal star-forming galaxies and their puffy starbursts should evolve from s0.8 to s1.0 between zs0 and zs2. The evolution for compact starbursts is predicted to be weaker, increasing from s0.5 to s0.6 between zs0 and zs2. In both cases the steepening of the radio spectra is attributed to the increase with redshift in IC losses from the CMB. In our analysis we do not find any significant increase in with redshift. However, our large measurements errors on α () and the existence of large intrinsic dispersion as revealed by studies of individually detected high-redshift galaxies (; Ibar et al. 2009, 2010) prevents us from ruling out the possibility that high-redshift star-forming galaxies follow the expectations of LT10 for normal star-forming galaxies and puffy starbursts. Further insight into any subtle evolution of α with redshift (as well as with Δlog(SSFR)_MS) will require deeper multi-frequency radio observations.

5.2. Evolution of the FIR/radio correlation

Our results indicate that the FRC evolves with redshift as q_FIR(z) = (2.35 ± 0.08) × (1 + z)^{−0.12 ± 0.04} (see Fig. 9). They also indicate that the FRC does not significantly evolve with Δlog(SSFR)_MS, though the presence of a weak positive trend, as observed in one of our redshift bin (i.e. Δ [ q_FIR ]/Δ[Δlog(SSFR)_MS ] = 0.22 ± 0.07 at 0.5 <z< 0.8), cannot be firmly ruled out using our dataset (see Fig. 8). These results are valid for all star-forming galaxies with M_∗ > 10¹⁰ M_⊙, Δlog(SSFR)_MS > −0.3 and 0 <z< 2.3.

5.2.1. AGN contamination?

The presence of a large population of AGNs in our sample at z > 1 might be a concern for our q_FIR estimates. Anticipating this problem, we have removed X-ray AGNs from our sample and used a median radio stacking analysis to statistically exclude the relatively small population of radio-loud AGNs. However, these precautions will minimise but not completely eliminate contamination by AGNs. To further reduce this potential contamination, we repeat our stacking analysis removing from our sample AGNs selected by the IRAC colour-colour criteria of Lacy et al. (2007)¹³. This selection should exclude of our analysis obscured AGNs missed in X-ray observations (Lacy et al. 2007; Stern et al. 2005; Donley et al. 2012). We found 11 865 IRAC AGNs in our sample, of which 4347 have individual mid-infrared detections. This corresponds to s4% of our final sample and s12% of our final sample with individual mid-infrared detections. Excluding IRAC AGNs does not change qualitatively and quantitatively our results: q_FIR smoothly decreases by s0.3 across 0 <z< 2.3 and there is no significant q_FIR– Δlog(SSFR)_MS correlation. We conclude that the decrease in q_FIR is most likely not driven by AGN contamination. Note that this absence of significant contamination from AGNs echoes results from Sargent et al. (2010; see also Bonzini et al., in prep.). Indeed, using individually detected FIR and radio sources, Sargent et al. (2010) found that AGNs (X-ray selected and optically selected) and star-forming galaxies follow the same FRC (in terms of normalisation and dispersion) out to at least zs1.4.

5.2.2. Expectations from theory

Theory predicts that the FRC is driven by star-formation activity in galaxies. The UV emission of young, massive (≳8 M_⊙) stars is absorbed by dust and re-emitted in the FIR. This creates a linear correlation between SFR and L_IR, if galaxies are optically thick at UV wavelengths. After few Myrs, young massive stars explode into SNe, accelerating CRs into the general magnetic field of galaxies and resulting in diffuse synchrotron emission. Thus, when averaged over the star-formation episode (neglecting the short lag between the UV emission and explosion of the first massive young stars), we expect a clear link between SFR, FIR and radio synchrotron emission from star-forming galaxies. This constitutes the essence of the calorimeter theory first proposed in Völk (1989). In this theory, galaxies are both UV and electron calorimeters, i.e. all UV radiation from young stars is re-emitted by dust in the FIR and all CR electrons are converted within the galaxies into an observable form, mainly via synchrotron CR cooling.

The calorimeter theory has however been questioned, especially in light of the remarkably tight FRC over three orders of magnitude in luminosity. L10 proposed that – on top of calorimetry – a number of additional physical processes conspire to yield such a tight FRC, covering dwarf galaxies and ULIRGs. At low SFR surface density (Σ_SFR¹⁴), the UV calorimeter assumption fails: galaxies are not optically thick and some UV photons escape without being re-emitted in the FIR (e.g. Buat et al. 2005). At the same time, the CR electrons calorimeter assumption also fails. CR electrons escape without radiating their energy in the radio, placing the galaxies back onto the FRC. At high Σ_SFR, CR cooling via Bremsstrahlung, ionisation and IC processes become more important because of higher gas densities. CR cooling via synchrotron competes with these processes, decreasing the radio synchrotron emission of compact starbursts. However, because of the higher gas densities, compact starbursts become CR proton calorimeters. CR protons convert their energy via inelastic scattering into gamma rays, neutrinos and secondary protons and electrons. These secondary protons and electrons undergo synchrotron cooling, placing compact starbursts back onto the FRC. From this, L10 conclude that calorimetry combines with several conspiracies operating in different density regimes to produce a relatively constant (variation <0.3) FRC across the range 0.001 M_⊙ kpc^-2 yr^-1< Σ_SFR< 1000 M_⊙ kpc^-2 yr^-1.

Other theories have been proposed to explain the FRC. For example Schleicher & Beck (2013, hereafter S13; see also Niklas & Beck 1997 explained the FRC by relating star formation and magnetic field strength via turbulent magnetic field amplification, the so-called small-scale dynamo effect. This model explains the FRC of galaxies both globally and on kiloparsec scales. Contrary to the findings of L10, this model predicts that the FRC should evolve with Σ_SFR as q_FIR = −0.3 × log (Σ_SFR) + C₁.

• Is the FRC expected to evolve with Δlog (SSFR)_MS? We can relate Δlog(SSFR)_MS with Σ_SFR using estimates presented in Wuyts et al. (2011b) and based on almost the same galaxy sample used here (see Sect. 2.4). Wuyts et al. (2011b) found that Σ_SFR increases linearly with Δlog(SSFR)_MS. At zs0.1, MS galaxies (Δlog(SSFR)_MSs0) have Σ_SFRs0.02 M_⊙ kpc^-2 yr^-1 while starbursts (Δlog(SSFR)_MSs1) have Σ_SFRs0.2 M_⊙ kpc^-2 yr^-1. At zs2.0, MS galaxies have Σ_SFRs1 M_⊙ kpc^-2 yr^-1 while starbursts have Σ_SFRs10 M_⊙ kpc^-2 yr^-1 (for more details, see Fig. 4 of Wuyts et al. 2011b). Thus, in a given redshift bin, our analysis probes at least one order of magnitude in Σ_SFR. Across this range, the model of L10 expects no significant evolution¹⁵ of q_FIR, while the model of S13 expects q_FIR to decrease by 0.3. Here, we find that, if it exists, a q_FIR– Δlog(SSFR)_MS correlation is necessary weak and rather positive (Δ [ q_FIR ]/Δ[Δlog(SSFR)_MS ] = 0.22 ± 0.07). This result disfavours the model of S13.

• Is the FRC expected to evolve with redshift? LT10 and S13 also study the evolution of the FRC with redshift. At high redshift the main concern is that other CR cooling processes might start to dominate over synchrotron. In particular, IC cooling from the CMB might become dominant at high redshift (U_CMB∝(1 + z)⁴), leading to a significant increase in q_FIR (see also Murphy 2009). However, both studies concluded that the FRC should hold with no dramatic break-down (i.e. q_FIR(z) − q_FIR(z = 0) ≳ 0.5) out to relatively high redshift: up to zs2 (8) for galaxies with Σs1 M_⊙ kpc^-2 yr^-1 and up to zs3 (15) for galaxies with Σs10 M_⊙ kpc^-2 yr^-1 in S13 (LT10). Because at zs2 our galaxies have Σ_SFR > 1 M_⊙ kpc^-2 yr^-1, none of these models expect a dramatic increase in q_FIR for our sample. These expectations are in line with our results. Apart from this potential break-down of the FRC due to IC cooling from the CMB, LT10 also study the possibility of more subtle variations of the FRC with redshift. Because SMGs exhibited low q_FIR values (q_FIRs2.0; Murphy et al. 2009) and steep radio spectra (αs0.8; Ibar et al. 2010), LT10 postulate that high-redshift SMGs have larger CR scale heights than their local counterparts (hs1 kpc instead of hs0.1 kpc). Then, assuming that the magnetic field strength varies with Σ_gas and not ρ_gas, they find that such large CR scale heights decrease the CR losses via Bremsstrahlung and ionisation processes, increasing the synchrotron emission and decreasing q_FIR by s0.3. LT10 argue that such puffy starbursts are also supported by kinematic observations of SMGs, citing Tacconi et al. (2006). From this, LT10 predict two different evolution of q_FIR across the range 0.001 M_⊙ kpc^-2 yr^-1< Σ_SFR< 1000 M_⊙ kpc^-2 yr^-1. In their normal-to-compact starburst track, q_FIR remains constant as a function of Σ_SFR, while in their normal-to-puffy starburst track q_FIR significantly evolves with Σ_SFR. In this latter, normal star-forming galaxies and puffy starbursts have the same CR scale height and q_FIR smoothly decreases by 0.3 from Σ_SFR = 0.001 M_⊙ kpc^-2 yr^-1 to Σ_SFR = 1 M_⊙ kpc^-2 yr^-1 and then q_FIRs2 at Σ_SFR ≳ 1 M_⊙ kpc^-2 yr^-1. Therefore, while LT10 do not formally predict a particular redshift evolution of q_FIR, they expect that if the CR scale height of star-forming galaxies evolve with redshifts, q_FIR should evolve accordingly. At zs2, all our galaxies have Σ_SFR > 1 M_⊙ kpc^-2 yr^-1 and their q_FIR have decreased by 0.3. These observations suggest that most high-redshift star-forming galaxies may have CR scale heights of s1 kpc, not only SMGs. This is in qualitative agreement with the finding that many high-redshift MS star-forming galaxies are clumpy disks with larger disk velocity dispersions (and hence scale heights) than local spirals (e.g. Förster Schreiber et al. 2009; Newman et al. 2013). S13 also predict some subtle variations of q_FIR with redshift following, (7)where β parametrises the evolution of the typical ISM density of galaxies with redshift, ρ = ρ₀ (1 + z)^β. For MS galaxies, Σ_SFR increases by 1.7 dex from zs0 to zs2. Assuming that their CR scale heights stay the same across this redshift range (hs1 kpc) and that with Ns1−2 from the Schmidt-Kennicutt relation, one would infer that with β_MSs1.8−3.6. Using the model of S13, we therefore predict that q_FIR for MS galaxies should decrease by −1.0 (−1.6) between zs0 and zs2 for N = 2(1). Such large evolution of q_FIR is not supported by our observations. Similarly, one can predict the evolution of q_FIR for starbursts galaxies (Δlog(SSFR)_MS ≳ 0.75). Their Σ_SFR also increase by 1.7 dex from zs0 to zs2. Assuming that their CR scale heights evolve from hs0.1 kpc to hs1 kpc from zs0 to zs2, we predict – using the model of S13 – a decrease in q_FIR by −0.4 (−1.0) for N = 2(1) across this redshift range¹⁶. Again, such large evolution of q_FIR is not supported by our observations.

• Is the FRC expected to evolve with redshift in the context of the MS of star formation? There is strong observational evidence that the physical conditions reigning in the star-forming regions of high-redshift MS galaxies are similar to those of normal local star-forming galaxies, though high-redshift MS galaxies have larger Σ_SFR (Wuyts et al. 2011b; Elbaz et al. 2011; Nordon et al. 2012; Magnelli et al. 2012b, 2014). Thus, MS galaxies should follow the normal-to-puffy starburst track of LT10 across 0 <z< 2.3, with its constant CR scale height of s1 kpc. This would translate into a smooth decrease in q_FIR by 0.3 at zs2 as their Σ_SFR varies from Σ_SFRs0.02 M_⊙ kpc^-2 yr^-1 to Σ_SFRs1 M_⊙ kpc^-2 yr^-1 across this redshift range. Such predictions are supported by our observations. In contrast, the decrease in q_FIR with redshift for far-above MS galaxies (those with Δlog(SSFR)_MS ≳ 0.75, see Sect. 4.2) is more surprising. Indeed, at high redshift, far-above MS galaxies are thought to be starbursts with similar FIR properties to local ULIRG (e.g. Elbaz et al. 2011; Magnelli et al. 2012b, 2014). Consequently, one could expect these galaxies to follow the compact starburst predictions of LT10 across 0 <z< 2.3, characterised by a constant value of q_FIR over a large range of Σ_SFR. Instead, their q_FIR seems to decrease by s0.3 at zs2, suggesting that most high-redshift starbursts (and not only SMGs) have larger CR scale heights (puffy starbursts) than their local counterparts (compact starbursts). Consequently, while high-redshift starbursts seems to share some physical properties with local ULIRGs (FIR SEDs; CO-to-H₂ conversion factors), other properties seem to be significantly different (size of their star-forming regions and CR scale heights). Note, however, that there is a larger dispersion of q_FIR at high Δlog(SSFR)_MS values and that the q_FIR–z correlation is statistically less significant for far-above MS galaxies than for MS galaxies (q_FIR [ z,far-above MS] ∝ (1 + z)^{−0.13 ± 0.09} and q_FIR [ z,MS ] ∝ (1 + z)^{−0.16 ± 0.05}). This supports a more complex scenario in which a non-negligible number of high-redshift starbursts are associated with compact starbursts rather than with puffy starbursts.

5.2.3. Implications of the evolution with redshifts of the FRC

Our results have important implications for previous and future studies relying on the FRC. For example, SFRs of high-redshift galaxies determined via radio observations will be lower than expected from the local FRC. However, while the evolution of q_FIR with redshift is statistically significant, it is modest and still lies within its local 1σ dispersion at zs2. Thus, corrections of SFRs motivated by our findings should still be within the uncertainties inferred using the local 1σ dispersion of the FRC.

Our results have also important implications for our understanding of the contribution of star-forming galaxies to the radio extragalactic background. In particular, an extragalactic 3.3 GHz radio background of νI_ν ≈ 5.9 × 10^-4nWm^-2sr^-1 has been detected with ARCADE2 (Fixsen et al. 2011; Seiffert et al. 2011). The sources of this cosmic radio background (CRB) have been extensively studied in the recent years (e.g. Gervasi et al. 2008; Singal et al. 2010; Vernstrom et al. 2011, 2014; Ysard & Lagache 2012; Condon et al. 2012). All these studies report that only a small fraction of this background (s10–26%) can be accounted for by objects resolved in the deepest existing radio surveys. Thus, they conclude that if the CRB at 3.3 GHz is at the level reported by ARCADE2, it must originate from very faint radio sources still to be observed. In particular, Singal et al. (2010) speculate that these sources could have radio flux densities <10 μJy and be ordinary star-forming galaxies at z > 1 associated with an FRC that has evolved towards radio-loud (numerically lower) values. Such evolution of the FRC is qualitatively supported by our observations. However, quantitatively, the evolution of q_FIR required by Singal et al. (2010) is much larger than that observed here, i.e. Δq_FIR< −0.7 instead of the observed Δq_FIRs−0.3 at zs2. From the total infrared luminosity density across 0 <z< 4 inferred by Gruppioni et al. (2013; see also Magnelli et al. 2013) using Herschel observations¹⁷, we can estimate the contribution of star-forming galaxies to the CRB at 3.3 GHz. Assuming S_ν ∝ ν^-0.8 and q_FIR = 2.34 across 0 <z< 4, we find nWm^-2sr^-1, i.e. 5.2% of the total CRB at 3.3 GHz. Instead, using q_FIR(z) = 2.35 × (1 + z)^-0.12 at 0 <z< 2.3 and q_FIR(z > 2) = q_FIR(z = 2), we find nWm^-2sr^-1, i.e. 7.6% of the total CRB at 3.3 GHz. Finally, using q_FIR(z) = 2.35 × (1 + z)^-0.12 at 0 <z< 4, we find nWm^-2sr^-1, i.e. 7.7% of the total CRB at 3.3 GHz. Therefore, while the evolution of the FRC with redshift increases the resolved fraction of the CRB as measured by ARCADE2, it cannot entirely explain its origin. We conclude that if the CRB is at the level reported here, the contribution of star-forming galaxies that obey the FRC is at most s10% (taking into account a dispersion of the FRC of 0.3). However, note that a proper treatment of outliers of the FRC should be performed in order to estimate the entire contribution of infrared sources to the CRB.

We also estimate the contribution of X-ray AGNs to the CRB at 3.3 GHz. In each of our redshift bin, we measure the mean 1.4 GHz flux density of X-ray AGNs using our stacking analysis. Then, to obtain the 1.4 GHz radio background of X-ray AGNs, we multiply this mean radio flux by the density of X-ray AGNs detected in our GOODS-N/S samples (i.e. fields with the deepest X-ray datasets). Finally, we convert this 1.4 GHz radio background into a 3.3-GHz radio background assuming a radio spectral index of X-ray AGNs in the range [0.5−1.0]. We found that X-ray AGNs have nWm^-2sr^-1, assuming a radio spectral index of 0.5 (1.0). This corresponds to only s1.9% (1.2%) of the total CRB measured by ARCADE2 at 3.3 GHz.

6. Summary

In this paper we study the evolution of the FRC and radio spectral index () across the SFR–M_∗ plane and up to zs2. We use the deepest FIR-Herschel, 1.4 GHz (VLA) and 610 MHz (GMRT) observations available for the GOODS-N, GOODS-S, ECDFS and COSMOS fields. Infrared luminosities are inferred using the stacked PACS/SPIRE FIR photometry for each SFR–M_∗–z bin. Radio luminosities and radio spectral indices are derived using their stacked 1.4 GHz and 610 MHz flux densities. Using this methodology, we are able to overcome most of the biases affecting previous studies on the FRC. Selection biases are overcome by the use of a complete stellar-mass-selected (M_∗ ≳ 10¹⁰ M_⊙) sample of star-forming galaxies at 0 <z< 2.3. Observational biases, introduced by sparse coverage of the FIR and radio spectra, are overcome by the use of multi-wavelength FIR and radio observations. Our results, which are valid for all star-forming galaxies with M_∗ > 10¹⁰ M_⊙, Δlog(SSFR)_MS > −0.3 and 0 <z< 2.3, can be summarised as follows:

• 1.

  The radio spectral index, , does not evolve significantly with the distance of a galaxy with respect to the MS (i.e. Δlog(SSFR)_MS) nor with redshift. Instead, remains relatively constant across the SFR–M_∗–z parameter space, consistent with a canonical value of 0.8. This results suggests that the radio spectra of the bulk of the high-redshift star-forming galaxy population is dominated by non-thermal optically thin synchrotron emission, well described by a power-law function, S_ν ∝ ν^-0.8, across the range of rest-frequencies probed here, 610 MHz <ν_rest< 4.2 GHz.

• 2.

  A relatively constant radio spectral index from normal star-forming galaxies (Δlog(SSFR)_MS ≳ 0) to starbursts (Δlog(SSFR)_MS ≳ 1) is surprising in light of local observations where the radio spectra of ULIRGs are significantly flatter (αs0.5) than those of spiral galaxies (αs0.8; Condon et al. 1991; Clemens et al. 2008). However, the observation of relatively steep radio spectra in high-redshift starbursts is also supported by results from Ibar et al. (2009) and Thomson et al. (2014) who found αs0.8 for samples of individually detected high-redshift SMGs. The combination of these results suggests that most high-redshift starbursts have different ISM properties (e.g. magnetic field strength, gas densities, Σ_SFR, ...) than their local counterparts. Alternatively, it could suggest that, as advocated by Condon et al. (1991), the flatter radio spectrum of local starbursts is due to free-free absorption, less relevant at the higher rest-frequencies probed at high redshift.

• 3.

  The FRC does not evolve significantly with Δlog(SSFR)_MS, though the presence of a weak positive trend, as observed in one of our redshift bin (i.e. Δ [ q_FIR ]/Δ[Δlog(SSFR)_MS ] = 0.22 ± 0.07 at 0.5 <z< 0.8), cannot be firmly ruled out using our dataset.

• 4.

  The FRC evolves with redshift as q_FIR(z) = (2.35 ± 0.08) × (1 + z)^{−0.12 ± 0.04}. This redshift evolution of the FRC is consistent with previous findings from the literature, though most favoured a non-evolving FRC because high-redshift measurements were within its local 1σ dispersion, q_FIR(z ≈ 0) = 2.34 ± 0.26 (Yun et al. 2001).

• 5.

  The fact that the FRC still holds at high redshift, albeit with some moderate evolution, suggests that IC cooling of CR electrons and protons off photons from the CMB (U_CMB ∝ (1 + z)⁴) does not yet dominate over synchrotron cooling at zs2. The redshift evolution of the FRC suggests that the ISM properties (e.g. magnetic field strength, gas densities, Σ_SFR, ...) of star-forming galaxies evolve between zs0 and zs2.

A decrease in q_FIR with redshift was expected by some of the most up-to-date theoretical models of the FRC, e.g. L10, LT10 and S13. Such evolution is expected because Σ_SFR and ρ_gas, which control in part the normalisation of the FRC, are known to evolve significantly with redshift. While the model of S13 seems to over-estimate the decrease in q_FIR with redshift, the bulk of the high-redshift star-forming galaxy population follows the expectations of LT10 by evolving smoothly from their normal-to-compact starbursts track at zs0 to their normal-to-puffy starbursts track at zs2. This suggests that MS galaxies have a constant CR scale height of s1 kpc across 0 <z< 2.3 and that their decrease by 0.3 at zs2 because their Σ_SFR changes from Σ_SFRs0.02 M_⊙ kpc^-2 yr^-1 to Σ_SFRs1 M_⊙ kpc^-2 yr^-1 in this redshift range. A constant CR scale height for MS galaxies is consistent with the current interpretation of the MS of star formation. Indeed, despite their extreme Σ_SFR, high-redshift MS galaxies have the physical properties of local normal star-forming galaxies (e.g. Wuyts et al. 2011b; Elbaz et al. 2011; Magnelli et al. 2012b, 2014). In addition, our results suggests that most starbursts have a CR scale height of s0.1 kpc at zs0 and of s1 kpc at zs2. This translates into a smooth decrease in q_FIR by 0.3 at zs2, as starbursts move from the compact starburst to the puffy starburst predictions of LT10. Nevertheless, this vision of an homogeneous high-redshift starburst population, constituted only of puffy starbursts, is compromised by the possibility of a weak q_FIR − Δlog(SSFR)_MS correlation, combined with a larger dispersion in q_FIR at high Δlog(SSFR)_MS. These latter observations support a scenario in which a non-negligible number of high-redshift starbursts are compact.

Our results have important implications for studies relying on the local FRC. For example, SFRs of high-redshift galaxies determined via radio observations will be lower than predicted by the local FRC. However, because q_FIR is still within its local 1σ dispersion at zs2, uncertainties on previous SFR estimates should capture this evolution.

It has been postulated in the past that the CRB detected by ARCADE2 at 3.3 GHz could be dominated by z > 1 star-forming galaxies obeying an evolved, radio-loud FRC (e.g. Singal et al. 2010). However, we find that star-forming galaxies responsible for the CIB contribute at most s10% of the CRB, even after having accounted for the evolution of the FRC observed here.

---

Acknowledgments

We thank the anonymous referee for suggestions which greatly enhanced this work. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/ OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including University of Lethbridge (Canada), NAOC (China), CEA, LAM (France), IFSI, University of Padua (Italy), IAC (Spain), Stockholm Observatory (Sweden), Imperial College London, RAL, UCL-MSSL, UKATC, University of Sussex (UK), Caltech, JPL, NHSC, University of Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). Support for BM was provided by the DFG priority programme 1573 The physics of the interstellar medium. R.J.I. acknowledges support from the European Research Council in the form of Advanced Grant, cosmicism. E.I. acknowledges funding from CONICYT/FONDECYT postdoctoral project N°:3130504. F.B. and A.K. acknowledge support by the Collaborative Research Council 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG).

References

All Tables

All Figures