The VLA-COSMOS 3 GHz Large Project: Average radio spectral energy distribution of highly star-forming galaxies

TL;DR

This study constructs the average radio spectral energy distribution of highly star-forming galaxies up to redshift 4, revealing a broken power-law shape that impacts radio luminosity calculations and the infrared-radio correlation trend.

Contribution

It provides the first detailed characterization of the radio SED of HSFGs across a broad redshift range, showing a broken power-law shape and its implications for star formation studies.

Findings

Radio SED of HSFGs follows a broken power-law with a break at 4.3 GHz

Spectral index steepens from 0.42 to 0.94 above 4.3 GHz

Shape of radio SED unlikely to cause the q-z trend in SFGs

Abstract

We construct the average radio spectral energy distribution (SED) of highly star-forming galaxies (HSFGs) up to z~4. Infrared and radio luminosities are bound by a tight correlation that is defined by the so-called q parameter. This infrared-radio correlation provides the basis for the use of radio luminosity as a star-formation tracer. Recent stacking and survival analysis studies find q to be decreasing with increasing redshift. It was pointed out that a possible cause of the redshift trend could be the computation of rest-frame radio luminosity via a single power-law assumption of the star-forming galaxies' (SFGs) SED.To test this, we constrained the shape of the radio SED of a sample of HSFGs. To achieve a broad rest-frame frequency range, we combined previously published VLA observations of the COSMOS field at 1.4 GHz and 3 GHz with unpublished GMRT observations at 325 MHz and 610 MHz by employing survival analysis to account for non-detections in the GMRT maps. We selected a sample of HSFGs in a broad redshift range (0.3<z<4,SFR>100M0/yr) and constructed the average radio SED. By fitting a broken power-law, we find that the spectral index changes from $\alpha_1=0.42\pm0.06$ below a rest-frame frequency of 4.3 GHz to $\alpha_2=0.94\pm0.06$ above 4.3 GHz. Our results are in line with previous low-redshift studies of HSFGs (SFR>10M0/yr) that show the SED of HSFGs to differ from the SED found for normal SFGs (SFR<10M0/yr). The difference is mainly in a steeper spectrum around 10 GHz, which could indicate a smaller fraction of thermal free-free emission. Finally, we also discuss the impact of applying this broken power-law SED in place of a simple power-law in K-corrections of HSFGs and a typical radio SED for normal SFGs drawn from the literature. We find that the shape of the radio SED is unlikely to be the root cause of the q-z trend in SFGs.

Tables

Table 1. Table 1: Rankings of different fitted models with their parameters. For simplicity, we show models without their normalization parameter and indicate the additional contribution of free-free emission with a thermal fraction at 1GHz by ‘FFE’. We show the rankings based on the Akaike and Bayesian information criteria and the log-likelihood test ( ln ⁡ L 𝐿 \ln L ).

Model	Normalized	Ranking
parameters	(normalized) flux density	AIC	BIC	$\ln L$
Broken power-lawa𝑎aa𝑎a$a$The broken power-law model was included in the table for comparison.	Eq. (3)	2	2	1
${\alpha }_1=0.53\pm 0.04$ (${\alpha }_1^{corr}=0.42\pm 0.06$)
${\alpha }_2=0.94\pm 0.06$
Model (1) - mixed modelb𝑏bb𝑏b$b$See Pacholczyk (1980) for details.	${\nu }^{-\alpha +2.1}\left(1-e^{-\frac{{\tau }_1}{{\nu }^{2.1}}}\right)$	6	5	6
$\alpha =0.74\pm 0.06$
${\tau }_1=0.2\pm 0.2$
Model (2) - foreground screen modelb𝑏bb𝑏b$b$See Pacholczyk (1980) for details.	${\nu }^{-\alpha }{e}^{-\frac{{\tau }_1}{{\nu }^{2.1}}}$	6	7	7
$\alpha =0.73\pm 0.06$
${\tau }_1=0.11\pm 0.08$
Model (3) - synchrotron aging modelb𝑏bb𝑏b$b$See Pacholczyk (1980) for details.	$\frac{{\nu }^{-\alpha }}{1+{\left(\frac{\nu }{{\nu }_b}\right)}^{\Delta \alpha }}$	4	4	5
$\alpha =0.24\pm 0.06$
$\Delta \alpha =0.55\pm 0.04$
Model (4) - synchrotron self-absorption model c𝑐cc𝑐c$c$See Condon & Ransom (2016) for details.	${\nu }^{5/2}\left(1-e^{-{\left(\frac{\nu }{{\nu }_b}\right)}^{-\alpha -5/2}}\right)$	8	8	8
$\alpha =0.70\pm 0.05$
${\nu }_b=0.76\pm 0.07\mathrm{GHz}$
Model (5) - (Lisenfeld et al. 2004) b𝑏bb𝑏b$b$See Pacholczyk (1980) for details.	${\nu }^{-\alpha }\left(1-{\left(1-\sqrt{\frac{\nu }{{\nu }_b}}\right)}^{1.1}\right)$	1	1	3
$\alpha =0.47\pm 0.03$
Model (6a) -synchrotron aging+FFEb𝑏bb𝑏b$b$See Pacholczyk (1980) for details.	$\frac{\left(1-f_{th}\left({\nu }_n\right)\right){\left(\frac{\nu }{{\nu }_n}\right)}^{-\alpha }}{1+{\left(\frac{\nu }{{\nu }_b}\right)}^{\Delta \alpha }}+f_{th}\left({\nu }_n\right){\left(\frac{\nu }{{\nu }_n}\right)}^{-0.1}$	7	5	4
$\alpha =0.24\pm 0.07$
$\Delta \alpha =0.55\pm 0.04$
$f_{th}\left(1\mathrm{GHz}\right)=0.00\pm 0.02$
Model (6b) - Lisenfeld et al. (2004)+FFEb𝑏bb𝑏b$b$See Pacholczyk (1980) for details.	$\left(1-f_{th}\left({\nu }_n\right)\right){\left(\frac{\nu }{{\nu }_n}\right)}^{-\alpha }\left(1-{\left(1-\sqrt{\frac{\nu }{{\nu }_b}}\right)}^{1.1}\right)+f_{th}\left({\nu }_n\right){\left(\frac{\nu }{{\nu }_n}\right)}^{-0.1}$	3	3	2
$\alpha =0.50\pm 0.06$
$f_{th}\left(1\mathrm{GHz}\right)=0.01\pm 0.02$

Table 2. Table 2: Summary of data reduction results of the GMRT maps. Columns b m ​ a ​ j ​ [ arcsec ] , subscript 𝑏 𝑚 𝑎 𝑗 delimited-[] arcsec b_{maj}\,\mathrm{[arcsec]}, b m ​ i ​ n ​ [ arcsec ] , and subscript 𝑏 𝑚 𝑖 𝑛 delimited-[] arcsec and b_{min}\,\mathrm{[arcsec]},\text{and} B P A [ ∘ ] BPA\,\mathrm{[^{\circ}]} show the beam major and minor axes and the beam position angle, respectively. Column Sources shows the number of sources that are detected toward the COSMOS field, column Resolved shows the number of resolved sources, column F ​ D ​ R 𝐹 𝐷 𝑅 FDR shows the false-detection rate, column Cross-matched shows the number of sources that are present in the respective GMRT catalog and in the VLA 3 ​ GHz 3 GHz 3\,\mathrm{GHz} catalog, and column HSFG Sample shows the number of sources in the subset of the GMRT catalog that is present in the highly star-forming galaxy sample described in Sect. 3 . Percentages in the table are derived relative to the column Sources .

Map	$b_{maj}\left[\mathrm{arcsec}\right]$	$b_{min}\left[\mathrm{arcsec}\right]$	$BPA{[}^{\circ }]$	Sources	Resolved	$FDR$	Cross-matched	HSFG Sample
$325\mathrm{MHz}$	$10.8$	$9.5$	$-62.8$	$634$	$177$ ($28.0\%$)	$3.0\%$	$567$ ($89.6\%$)	29
$610\mathrm{MHz}$	$5.6$	$3.9$	$65.5$	$999$	$196$ ($19.6\%$)	$1.3\%$	$870$ ($87.1\%$)	52

Table 3. Table 3: Average corrections of luminosity at 1.4 ​ GHz 1.4 GHz 1.4\,\mathrm{GHz} needed if the K-correction was performed by assuming a power-law SED with a spectral index of 0.7 0.7 0.7 . The corrections were derived for the average luminosity at 1.4 ​ GHz 1.4 GHz 1.4\,\mathrm{GHz} , as described in Sect. B , and are based on the assumption that galaxies with an SFR > 10 ​ M ⊙ / yr SFR 10 subscript M direct-product yr \,\mathrm{SFR}>10\,\mathrm{M_{\odot}}/\,\mathrm{yr} follow the broken power-law (BPL) SED and that galaxies with an SFR < 10 ​ M ⊙ / yr SFR 10 subscript M direct-product yr \,\mathrm{SFR}<10\,\mathrm{M_{\odot}}/\,\mathrm{yr} follow the SED based on Eq. ( 2 ), with the nonthermal spectral index and thermal fraction taken from Tabatabaei et al. ( 2017 ) . The percentage p < 10 subscript 𝑝 absent 10 p_{<10} is the percentage of 3 ​ GHz 3 GHz 3\,\mathrm{GHz} flux that is due to galaxies with an SFR < 10 ​ M ⊙ / yr SFR 10 subscript M direct-product yr \,\mathrm{SFR}<10\,\mathrm{M_{\odot}}/\,\mathrm{yr} . Columns of p < 10 = 0 % subscript 𝑝 absent 10 percent 0 p_{<10}=0\% and p < 10 = 100 % subscript 𝑝 absent 10 percent 100 p_{<10}=100\% are also exact corrections for single objects derived by assuming the BPL or the Tabatabaei et al. ( 2017 ) SEDs, respectively. The values in parentheses are correction C L 1.4 ​ GHz subscript 𝐶 subscript 𝐿 1.4 GHz C_{L_{1.4\,\mathrm{GHz}}} expressed in dex dex \mathrm{dex} .

$z$	$C_{L_{1.4\mathrm{GHz}}}$ ($\log C_{L_{1.4\mathrm{GHz}}}\left[\mathrm{dex}\right]$)
	0%	25%	50%	75%	100%
0.0	0.81 (-0.09)	0.89 (-0.05)	0.97 (-0.02)	1.04 ( 0.02)	1.12 ( 0.05)
0.5	0.74 (-0.13)	0.85 (-0.07)	0.95 (-0.02)	1.06 ( 0.02)	1.17 ( 0.07)
1.0	0.79 (-0.10)	0.89 (-0.05)	0.99 (-0.01)	1.08 ( 0.03)	1.18 ( 0.07)
1.5	0.83 (-0.08)	0.92 (-0.04)	1.01 ( 0.00)	1.10 ( 0.04)	1.18 ( 0.07)
2.0	0.87 (-0.06)	0.95 (-0.02)	1.02 ( 0.01)	1.10 ( 0.04)	1.18 ( 0.07)
2.5	0.90 (-0.04)	0.97 (-0.01)	1.04 ( 0.02)	1.10 ( 0.04)	1.17 ( 0.07)
3.0	0.93 (-0.03)	0.99 (-0.00)	1.04 ( 0.02)	1.10 ( 0.04)	1.15 ( 0.06)
3.5	0.96 (-0.02)	1.01 ( 0.00)	1.05 ( 0.02)	1.09 ( 0.04)	1.14 ( 0.06)
4.0	0.99 (-0.01)	1.02 ( 0.01)	1.06 ( 0.02)	1.09 ( 0.04)	1.12 ( 0.05)

Equations

$$F\left({\nu }_r|\begin{array}{c}{\alpha }_2\\ b\\ f\end{array}\right)=10^b\left[{\left(\frac{{\nu }_r}{{\nu }_n}\right)}^{-{\alpha }_2}(1-f)+f{\left(\frac{{\nu }_r}{{\nu }_n}\right)}^{-0.1}\right],$$

$$F\left({\nu }_r|\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ b\\ f\\ {\nu }_b\end{array}\right)=\{\begin{array}{ll}10^b\left[{\left(\frac{{\nu }_r}{{\nu }_n}\right)}^{-{\alpha }_2}(1-f)+f{\left(\frac{{\nu }_r}{{\nu }_n}\right)}^{-0.1}\right],& {\nu }_r>{\nu }_b\\ 10^b\left[{\left(\frac{{\nu }_r}{{\nu }_n}\right)}^{-{\alpha }_1}(1-f){\left(\frac{{\nu }_b}{{\nu }_n}\right)}^{{\alpha }_1-{\alpha }_2}+f{\left(\frac{{\nu }_r}{{\nu }_n}\right)}^{-0.1}\right],& {\nu }_r<{\nu }_b\end{array},$$

$$\log F\left({\nu }_r|\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ b\\ {\nu }_b\end{array}\right)=\{\begin{array}{ll}-{\alpha }_2\log \frac{{\nu }_r}{{\nu }_n}+b,& {\nu }_r>{\nu }_b,\\ -{\alpha }_1\log \frac{{\nu }_r}{{\nu }_n}+b+({\alpha }_1-{\alpha }_2)\log \frac{{\nu }_b}{{\nu }_n},& {\nu }_r<{\nu }_b,\end{array}.$$

$$\ln \mathcal{L}(\theta )\propto -\frac{1}{2}\sum _i{\left(\frac{⟨\log F{⟩}_i-\log F({\nu }_i,\theta )}{{\sigma }_{\log F_i}}\right)}^2,$$

$$q(z)=q_0(1+z{)}^n,$$

$$L({\nu }_{e_2})=\frac{4\pi D_L^2}{1+z}S({\nu }_{o_2}),$$

$$\frac{L({\nu }_{e_2})}{L_{(}{\nu }_{e_1})}=\frac{S({\nu }_{e_2})}{S_{(}{\nu }_{e_1})}=\frac{S({\nu }_{o_2})}{S_{(}{\nu }_{o_1})}.$$

$$L({\nu }_{e_2})=\frac{4\pi D_L^2}{1+z}\frac{S({\nu }_{e_2})}{S\left({\nu }_{o_1}(1+z)\right)}S({\nu }_{o_1}).$$

$$\lg S_k=-{\alpha }_k\lg \nu +b_k,$$

$$b_k=b_0+\sum _{j=0}^{k-1}({\alpha }_{k-j}-{\alpha }_{k-(j+1)})\lg {\nu }_{k-j}=b_0+\sum _{j=0}^{k-1}\Delta {\alpha }_{k-j}\lg {\nu }_{k-j}.$$

$$\frac{S({\nu }_{e_2})}{S\left({\nu }_{o_1}(1+z)\right)}=\frac{1}{(1+z{)}^{-\alpha ({\nu }_{o_1}(1+z))}}\left(\frac{{\nu }_{e_2}^{-\alpha ({\nu }_{e_2})}}{{\nu }_{o_1}^{-\alpha ({\nu }_{o_1}(1+z))}}\right)10^{b_k({\nu }_{e_2})-b_k({\nu }_{o_1}(1+z))}.$$

$$L({\nu }_{e_2})=\frac{4\pi D_L^{2}S({\nu }_{o_1})}{(1+z{)}^{1-\alpha ({\nu }_{o_1}(1+z))}}\left(\frac{{\nu }_{e_2}^{-\alpha ({\nu }_{e_2})}}{{\nu }_{o_1}^{-\alpha ({\nu }_{o_1}(1+z))}}\right)10^{b_{k({\nu }_{e_2})}-b_{k({\nu }_{o_1}(1+z))}}.$$

$$L_{1.4\mathrm{GHz}}=4\pi D_L^2{S}_{3\mathrm{GHz}}\{\begin{array}{l}\frac{1}{(1+z{)}^{1-{\alpha }_1}}\left(\frac{1.4^{-{\alpha }_1}}{3^{-{\alpha }_1}}\right),3\mathrm{GHz}(1+z)<{\nu }_b\\ \frac{1}{(1+z{)}^{1-{\alpha }_2}}\left(\frac{1.4^{-{\alpha }_1}}{3^{-{\alpha }_2}}\right)10^{-({\alpha }_2-{\alpha }_1)\log {\nu }_b},3\mathrm{GHz}(1+z)>{\nu }_b\end{array}.$$

$$L_{1.4\mathrm{GHz}}=\frac{4\pi D_L^2}{1+z}\frac{S_{3\mathrm{GHz}}}{(1-f_{th}){\left(\frac{3(1+z)}{1.4}\right)}^{-{\alpha }_2}+f_{th}{\left(\frac{3(1+z)}{1.4}\right)}^{-0.1}}.$$

$$⟨L_{1.4\mathrm{GHz}}⟩(z)=\int f(\mathrm{SFR},z)K(\mathrm{SFR},z)S(\mathrm{SFR},z)\mathrm{d}\mathrm{SFR}.$$

$$C_{L_{1.4\mathrm{GHz}}}(p_{<10},z)=(1-p_{<10}(z))\frac{K_{BPL}(z)}{K_{\alpha =0.7}(z)}+p_{<10}(z)\frac{K_{Tab}(z)}{K_{\alpha =0.7}(z)},$$

$$p_{<10}=\frac{\underset{<10{\mathrm{M}}_{\odot }/\mathrm{yr}}{\int }f(\mathrm{SFR},z)S(\mathrm{SFR},z)\mathrm{d}\mathrm{SFR}}{\int f(\mathrm{SFR},z)S(\mathrm{SFR},z)\mathrm{d}\mathrm{SFR}}.$$

Full text

1 Introduction

Star formation rate (SFR) measurements derived from the ultraviolet and optical emission are sensitive to dust (Kennicutt 1998a), and infrared-derived SFRs are easy to understand only in the optically thick case (Kennicutt 1998a; Condon 1992). On the other hand, radio emission in star-forming galaxies (SFGs) below $\sim 30\,\mathrm{GHz}$ is unbiased by dust (Condon 1992). Furthermore, it has been found that the infrared and radio luminosities are bound by a tight correlation over many orders of magnitude (van der Kruit 1971; Helou et al. 1985; Condon 1992). This infrared-radio correlation provides the basis for radio luminosity as a star-formation tracer, as infrared luminosity is linked to the SFR (Kennicutt 1998b). It is defined by the so-called $q$ parameter, defined as $q=\log L_{TIR}/3.75\,\mathrm{THz}-\log(L_{1.4\,\mathrm{GHz}})$, where $L_{\mathrm{TIR}}$ is the total infrared luminosity ($8-1000\,\mathrm{\mu m}$) in units of $\,\mathrm{W}$ and $L_{1.4\,\mathrm{GHz}}$ is the luminosity density at $1.4\,\mathrm{GHz}$ in units of $\,\mathrm{W}/\,\mathrm{Hz}$. A frequency of $3.75\,\mathrm{THz}$ is used to obtain a dimensionless quantity. In the calorimeter model (Voelk 1989), $q$ is expected to be constant over a wide range of magnetic field strengths through electrons that lose their energy because of synchrotron and inverse Compton losses. This produces a tight correlation, even though individual galaxies may show deviations. Theoretical considerations (Murphy 2009) predict $q$ to increase with redshift, mostly because of the synchrotron electrons’ inverse Compton scattering off cosmic microwave background photons. Recent stacking and survival analysis studies find $q$ to be decreasing with increasing redshift (Ivison et al. 2010; Magnelli et al. 2015; Delhaize et al. 2017; Calistro Rivera et al. 2017). Delhaize et al. (2017) pointed out that the computation of rest-frame radio luminosity via a K-correction using only a simple, single power-law assumption of the star-forming galaxies’ spectral energy distribution (SED) might cause such a trend.

The exact shape of the radio SED of star-forming galaxies is usually assumed to be a superposition of the steep synchrotron spectrum, described by a power law, and a $\sim 10\%$ contribution at $1.4\,\mathrm{GHz}$ of a flat free-free spectrum (Condon 1992). These results are supported by observations of nearby galaxies at $3-30\,\mathrm{Mpc}$ ($z\approx 0.01$) with an $\,\mathrm{SFR}<10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ (Tabatabaei et al. 2017).

For (ultra)luminous infrared galaxies ((U)LIRGs, $\,\mathrm{SFR}>10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$), observations in the local Universe find their radio SEDs to be steep around $10\,\mathrm{GHz}$, corresponding to a spectrum that shows less flattening above $10\,\mathrm{GHz}$ than in galaxies with an $\,\mathrm{SFR}<10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ (Clemens et al. 2008; Leroy et al. 2011; Galvin et al. 2018). Both free-free self-absorption and modifications to the spectrum due to the aging of the electron population would also be present. These processes result in deviations from a simple power-law at low ($\nu_{rest}\lesssim 1\,\mathrm{GHz}$) and high ($\nu_{rest}\gtrsim\mathrm{a\,few}\,\mathrm{GHz}$) frequencies that are observationally hard to constrain.

In this paper, we constrain the shape of the average radio spectral energy distribution (SED) of a $1.4\,\mathrm{GHz}$-selected sample ($\,\mathrm{SFR}>100\,\mathrm{M_{\odot}}/\,\mathrm{yr}$) of highly star-forming galaxies (HSFGs) in the $\sim 0.4-10\,\mathrm{GHz}$ frequency range within the COSMOS field. To obtain a broad frequency range, we use the fluxes from the $1.4\,\mathrm{GHz}$ Very Large Array (VLA)-COSMOS Joint Project (Schinnerer et al. 2010) and the $3\,\mathrm{GHz}$ VLA-COSMOS Large Project (Smolčić et al. 2017a, b) and reduce previously unpublished Giant Meterwave Radio Telescope (GMRT) observations at $325\,\mathrm{MHz}$ and $610\,\mathrm{MHz}$. We construct the average radio SED by using survival analysis to account for the lower sensitivity of the GMRT maps.

In Sect. 2 we describe the VLA and GMRT datasets used in this paper, in Sect. 3 we describe the galaxy sample of HSFGs, in Sect. 4 we construct the average radio SED for HSFGs and test the applicability of our method. In Sect. 5 we report the average radio spectral indices of HSFGs, while in Sect. 6 we discuss the shape of the radio spectrum of SFGs and the impact on the infrared-radio correlation. In this paper we use the convention for the spectral index $F_{\nu}\sim\nu^{-\alpha}$, where $F_{\nu}$ is the flux density and $\alpha$ is the spectral index.

2 Data

To assess the radio spectral energy distribution, we combined catalogs of the COSMOS field at four observer-frame frequencies: $325\,\mathrm{MHz}$ and $610\,\mathrm{MHz}$, obtained with the GMRT, and at $1.4\,\mathrm{GHz}$ (Schinnerer et al. 2010) $\text{and }3\,\mathrm{GHz}$ (Smolčić et al. 2017a), obtained by the VLA. Here we briefly describe these VLA (Sect. 2.1 ) and GMRT (Sect. 2.2) data.

2.1 Available VLA data

The VLA COSMOS $3\,\mathrm{GHz}$ Large Project map (Smolčić et al. 2017a) was constructed from $384\,\mathrm{hours}$ of observations of the $2$ deg2 COSMOS field. Observations were carried out over $192$ pointings in the S-band with VLA in A and C antenna configurations. This was a wide-band observation with a total bandwidth of $2084\,\mathrm{MHz}$ derived from $16$ spectral windows, each $128\,\mathrm{MHz}$ wide. The final mosaic reached a median root-mean-square (RMS) of $2.3\,\mathrm{\mu Jy}/\mathrm{beam}$ at a resolution of $0.75\arcsec$. Considering the large bandwidth and the volume of the data, each pointing was imaged separately using the multi-scale multi-frequency synthesis algorithm (hereafter MSMF, Rau & Cornwell 2011) and then combined into a single mosaic. The MSMF algorithm had been found to be optimal for both resolution and image quality in terms of the RMS noise and sidelobe contamination (Novak et al. 2015). The $3\,\mathrm{GHz}$ catalog has been derived using blobcat (Hales et al. 2012) with a $5\sigma$ threshold. In total, 10830 sources were recovered (67 of which were multi-component sources). For more details, see Smolčić et al. (2017a).

The $1.4\,\mathrm{GHz}$ catalog is a joint catalog comprised of the VLA COSMOS $1.4\,\mathrm{GHz}$ Large (Schinnerer et al. 2007) and Deep surveys (Schinnerer et al. 2010). The joint catalog was constructed by combining the observations of the central $50\arcmin\times 50\arcmin$ subregion of the COSMOS field (VLA-COSMOS Deep Project Schinnerer et al. 2010) at a resolution of $2.5\arcsec\times 2.5\arcsec$ and previously published observations of the entire $2\deg^{2}$ COSMOS field at a resolution of $1.5\arcsec\times 1.4\arcsec$ (VLA-COSMOS Large Project Schinnerer et al. 2007). The average RMS in the resulting map was found to be $12\,\mathrm{\mu Jy}$. The large survey consisted of $240\,\mathrm{hours}$ of observations over $23$ pointings spread out over the entire COSMOS field in the L band centered at $1.4\,\mathrm{GHz}$ and with a total bandwidth of $37.5\,\mathrm{MHz}$ in VLA A configuration and $24\,\mathrm{hours}$ in the C configuration. Schinnerer et al. (2010) supplemented the $7$ central pointings with an additional $8.25\,\mathrm{hours}$ of observations per pointing using the A configuration and the same L-band configuration. These new measurements were then combined in the uv-plane with the Large Project observations. The joint catalog was constructed by using the SExtractor package (Bertin & Arnouts 1996) and the AIPS task SAD, yielding 2865 sources (Schinnerer et al. 2010).

2.2 GMRT data

The GMRT $325$ and $610\,\mathrm{MHz}$ data reduction, imaging, cataloging, and testing is described in detail in Appendix A. Here we give a brief overview of the procedure.

The overall data reduction procedure is similar for both observing frequencies. The observations were carried out with $30$ antennas with their longest baseline being $25\,\mathrm{km}$. The data were then reduced by using the source peeling and atmospheric modeling pipeline (SPAM, described in detail in Intema et al. 2017). Finally, the source extraction package blobcat was used to extract sources down to $5\sigma$, where $\sigma$ is the local RMS value computed from RMS maps produced by the AIPS task RMSD. The values of $\sigma$ are shown in the visibility functions, presented in top left panels of Figs. (11) and (12) for the $325\,\mathrm{MHz}$ and $610\,\mathrm{MHz}$ RMS maps, respectively. The fluxes were corrected for bandwidth-smearing, and resolved sources were identified.

The $325\,\mathrm{MHz}$ observations of a single pointing were carried out under the project 07SCB01 (PI: S. Croft) (at a reference frequency of $325\,\mathrm{MHz}$ and a bandwidth of $32\,\mathrm{MHz}$). The observations lasted for $45\,\mathrm{hours}$ in total, comprising four observations with a total on-source time of $\sim 40\,\mathrm{hours}$. They were reduced using the SPAM pipeline and imaged at a resolution of $10.8\times 9.5\,\mathrm{arcsec}^{2}$. A primary beam correction was applied to the pointing. We measured a median RMS of $97\,\mathrm{\mu Jy}/\mathrm{beam}$ over the $\sim 2\deg^{2}$ COSMOS field. In total, 633 sources were identified using blobcat down to $5\sigma$. By employing a peak-over-total flux criterion, we consider 177 of these sources to be resolved.

The $610\,\mathrm{MHz}$ observations were carried out at a central frequency of $608\,\mathrm{MHz}$ using a bandwidth of $32\,\mathrm{MHz}$. Observations lasting $86\,\mathrm{hours}$ (spread over eight observations with an average time on source per pointing of $\sim 4\,\mathrm{hours}$) were conducted under the project 11HRK01 (PI: H. R. Klöckner). A total of $19$ pointings were observed. The data were reduced using the SPAM pipeline, and imaged at a resolution of $5.6\times 3.9\,\mathrm{arcsec}^{2}$. Following Intema et al. (2017), primary beam correction and the average pointing error corrections were applied to each pointing prior to mosaicing. We measured a median RMS of $39\,\mathrm{\mu Jy}/\mathrm{beam}$ in the final mosaic over the $\sim 2\deg^{2}$ COSMOS field. Blobcat was used to extract 999 sources down to $5\sigma$, 196 of which we consider to be resolved.

3 Sample

To construct a complete and pure sample of star-forming galaxies over a wide range in redshift, we used the criteria described in Smolčić et al. (2017b), which we briefly summarize below. Starting from the VLA COSMOS Large Project $3\,\mathrm{GHz}$ Source Catalog, cross-correlated with multi-wavelength counterparts detected in the COSMOS field, we selected the clean star-forming galaxy sample (CSFG), as defined in Smolčić et al. (2017b), and imposed an additional cut of infrared-derived SFRs of $\,\mathrm{SFR}>100\,\mathrm{M_{\odot}}\,\mathrm{yr}^{-1}$.

Smolčić et al. (2017b) cross-correlated the $3\,\mathrm{GHz}$ Source Catalog with multi-wavelength counterparts, drawn from three catalogs: i) the COSMOS2015 catalog (Laigle et al. 2016), which contains sources detected on a $\chi^{2}$ sum of $Y$, $J$, $H$, $K_{s}$ and $z^{++}$ images, ii) the i-band catalog (Capak et al. 2007), which contains sources detected from a combination of the CFHT $i^{*}$ and Subaru $i^{+}$original point spread function (PSF) images, and iii) the IRAC catalog, which contains sources detected in the 3.6 $\mu$m band (Sanders et al. 2007). Smolčić et al. (2017b) combined counterparts from the three catalogs to maximize the number of counterparts as the i-band (IRAC) catalog is sensitive to potentially blue (red) sources undetected in the $Y$, $J$, $H$, $K_{s}$ , $z^{++}$ stack. Over the inner $1.77\,\mathrm{deg}^{2}$ subarea of the COSMOS field, which excludes all regions affected by saturated or bright sources in the optical to near-infrared bands, and thus contains sources with the best possible photometry, Smolčić et al. (2017b) found 7729 COSMOS2015, 97 i-band, and 209 IRAC counterparts to the 8696 3 GHz sources within this area. The COSMOS2015 and i-band catalogs offer precise photometric redshifts, with a joint precision of $\sigma_{\Delta z/(1+z)}=0.01$.

The CSFG sample was constructed from the full $3\,\mathrm{GHz}$ galaxy sample with COSMOS2015 or i-band counterparts as the cleanest sample of star-forming galaxies by excluding AGN through a combination of criteria, such as X-ray luminosity, mid-infrared (color-color and SED-based) criteria, $r^{+}-$NUV colors, and ($>3\sigma$) excess in the distribution of the ratio of 1.4 GHz radio luminosity and IR-based SFRs as a function of redshift (see Sect. 6.4 and Fig. 10 in Smolčić et al. 2017b). To determine the physical properties of the galaxies, a three-component SED-fitting procedure (Delvecchio et al. 2017) was applied using all of the available photometry. The fitting procedure was based on MAGPHYS (da Cunha et al. 2008) and SED3FIT (Berta et al. 2013) to account for the AGN component of the SED. SFRs were computed from the 8–1000 $\mathrm{\mu m}$ rest-frame infrared luminosity that was obtained by using the best-fit galaxy template and was then converted into SFR through the Kennicutt (1998a) conversion factor and scaled to a Chabrier (2003) initial mass function. Star formation rates are shown for the entire $3\,\mathrm{GHz}$ Large Project catalog and the subset of highly star-forming galaxies ($\,\mathrm{SFR}\geq 100\,\mathrm{M_{\odot}}/\,\mathrm{yr}$) in Fig. (1) as a function of redshift. The figure shows that a cut above $100\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ minimizes incompleteness above $z\approx 2$, which strongly influences the sample below $100\,\mathrm{M_{\odot}}/\,\mathrm{yr}$. To check whether the remaining incompleteness could influence our results, we repeated our analysis of Sect. 5 for the more complete subsample of HSFGs out to $z\sim 2$. We find no significant difference in the resulting SED and therefore conclude that the $100\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ cut sufficiently reduces incompleteness out to $z\sim 4$.

Finally, we cross-matched the $3\,\mathrm{GHz}$ catalog with the above described cuts and the $1.4\,\mathrm{GHz}$ joint catalog, to arrive at a sample of $306$ galaxies that have two reliably determined fluxes in the rest-frame frequency-flux space. We call the subset of the CSFG sample that is detected at both $1.4$ and $3\,\mathrm{GHz}$ and has $SFR>100\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ the highly star-forming galaxy (HSFG) sample. The $325\,\mathrm{MHz}$ and $610\,\mathrm{MHz}$ catalogs were then cross-matched with the HSFG sample catalog with a matching radius of 1 beam, leading to 29 and 52 sources in the HSFG sample, respectively. The summary of the number of detections and non-detections in the catalogs used to construct the radio SED of the HSFG sample is provided in Fig. (2).111The table of the cross-matched fluxes is available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/.

In summary, the selected sample is effectively a $1.4\,\mathrm{GHz}$ data-selected HSFG ($\,\mathrm{SFR}>100\,\mathrm{M_{\odot}}/\,\mathrm{yr}$). With this dataset, we associated the $3\,\mathrm{GHz}$ data, with a much lower RMS ($\sim 12\,\mathrm{\mu Jy/beam}$ and $\sim 2\,\mathrm{\mu Jy/beam}$ at $1.4\,\mathrm{GHz}$ and $3\,\mathrm{GHz}$, respectively), minimizing potential biases that could arise in the $1.4$ to $3\,\mathrm{GHz}$ spectral index distribution, given that almost all $1.4\,\mathrm{GHz}$ sources ( 90$\%$, Smolčić et al. 2017a) have counterparts at $3\,\mathrm{GHz}$. In total, the HSFG sample contains 306 galaxies, out of which $9\%$ (29) were detected at $325\,\mathrm{MHz}$ and $17\%$ (52) at $610\,\mathrm{MHz}$.

4 Average radio SED

In this section, we describe the method used to construct the radio SED and the models used to analyze it. In Sect. 4.1 we describe how we constructed the radio SED of HSFGs from the GMRT and VLA observations of the COSMOS field while accounting for lower sensitivities of the GMRT maps using survival analysis. In Sects. 4.2 and 4.3 we describe the models of the radio SED and the Monte Carlo simulations that test the applicability of the fitting procedure, respectively.

4.1 Construction of the average radio SED of HSFGs

To construct a typical SED for HSFGs in the COSMOS field, we normalized the spectra of individual galaxies at a particular rest-frame frequency. Considering that the redshifts of sources in our sample vary from $0.3$ to $4$, we obtain a well-constrained SED from rest-frame frequencies of $\sim 0.5\,\mathrm{GHz}$ to $\sim 15\,\mathrm{GHz}$.

A complication to this basic idea is that there are very few detections of HSFGs in the $325$ and $610\,\mathrm{MHz}$ maps. For this reason, we normalized the flux densities at all frequencies based on a linear fit to the $1.4$ and $3\,\mathrm{GHz}$ spectra. We chose the normalization rest-frame frequency to be the median log-frequency of sources in the HSFG sample at $1.4\,\mathrm{GHz}$ and $3\,\mathrm{GHz}$. For the HSFG sample, this frequency is $5.7\,\mathrm{GHz}$. We constrained the $325$ and $610\,\mathrm{MHz}$ flux densities of non-detections by five times the local RMS flux density of their respective maps. The problem is therefore to constrain the spectral index when a significant portion of the data is left-censored, meaning that many of the data points below a rest-frame frequency of $1\,\mathrm{GHz}$, corresponding to observer frame $325$ and $610\,\mathrm{MHz}$ frequencies, are upper limits. To correctly deal with upper limits, we decided to employ the survival analysis technique described in the following paragraph and illustrated in Fig. (3).

We combined into a single dataset the normalized fluxes and normalized upper limits of individual galaxies at different rest-frame frequencies. As we have four different observer-frame frequencies, the number of our data points is four times larger than the number of galaxies. To achieve uniform frequency binning, we used 20 bins that were equally separated in log space of the rest-frame frequencies. The number of bins was chosen empirically to achieve a satisfactory number of bins while retaining a sufficient number of data points within each bin for the statistical analysis.

For each bin, $i$, if there were no upper limits within the bin, we computed the mean of the log rest-frame frequency, $\nu_{i}=10^{\langle\log\nu_{r}\rangle_{i}}$, and its standard deviation, $\sigma_{\log\nu_{i}}$, along with the mean normalized log flux, $\langle\log F_{\nu}\rangle_{i}$, and its standard deviation, $\sigma_{\log F_{i}}$. If there were upper limits within a particular bin, we constructed the Kaplan-Meier estimate of the survival function of fluxes within the bin using the lifelines package222See https://lifelines.readthedocs.io/.. We did this by converting this left-censored dataset into a right-censored one in accordance with Feigelson & Nelson (1985) as the lifelines package handles right-censored datasets.

We then fit the survival function with the Weibull distribution (Lai et al. 2006). We have found this method to be satisfactory because the Weibull distribution survival function, $S(x|\lambda,k)=e^{-(x/\lambda)^{k}}$, where $\lambda$ and $k$ are free parameters of the model, can be easily linearized. We estimated $\langle{\log F_{\nu}}\rangle$ from $k$ and $\lambda$ and calculated $\sigma_{\log F}$ as the standard deviation of the Weibull distribution.

Survival analysis assumes that upper limits have the same underlying distribution as the data, and the Kaplan-Meier estimator further assumes that detections and upper limits are mutually independent (Kaplan & Meier 1958). The latter condition is satisfied by normalizing spectra to one frequency, as the normalization redistributes both the data and the upper limits based on the source’s normalization flux. Further tests of this procedure are described in Sect. 4.3 using Monte Carlo simulations of the properties of the HSFG sample.

4.2 SED fitting

Previous studies, listed in Sect. 1, have shown that the shape of the radio SED deviates from a simple power law. As the SED at higher frequencies ($\nu\gtrsim 10\,\mathrm{GHz}$) is expected to have a contribution not only by synchrotron, but also by free-free emission, the simplest modification to a simple power-law would be adding a thermal contribution with a flat spectral index of $0.1$ (Condon 1992) that would produce a flattening of the SED above $10\,\mathrm{GHz}$:

[TABLE]

where $\nu_{r}$ is the rest-frame frequency, and the parameters of the model are $\alpha_{2}$, the nonthermal spectral index, $b$, the normalization constant, $\nu_{n}$, the normalization frequency, and $f=f(\nu_{n})$, the thermal fraction at $\nu_{n}$. In our HSFG sample, however, we observe that the SED flattens below $3-4\,\mathrm{GHz}$ (see Sect. 5 and figures therein). To account for this, we assumed a broken power-law instead of using only the nonthermal spectral index. Therefore, the full SED would be

[TABLE]

where $\nu_{n}$ is the normalization frequency, $\nu_{b}$ the break frequency, $\alpha_{1}$ ($\alpha_{2}$) the spectral index below (above) $\nu_{b}$ , and $b$ the normalization constant.

In the HSFG sample, as described in Sect. 6.1, we do not find $f$ to be significantly greater than 0 for $\nu_{n}\sim 1\,\mathrm{GHz}$. We therefore simplified Eq. (2) to a broken power-law model for the log-normalized flux

[TABLE]

In the broken power-law model, the break frequency, $\nu_{b}$, can be both a free parameter or fixed to a particular value prior to fitting. The break frequency is derived using the Markov chain Monte Carlo (MCMC) method (implemented using the affine invariant sampler of the emcee package Foreman-Mackey et al. 2013) on the broken power-law fit with the break frequency as a free parameter. By inspecting the distribution of limits and detections, we find that the break frequency should lie above the highest frequency at which the last upper limit occurs. We thereby conclude that the break frequency would not be significantly affected by upper limits.

The MCMC code samples the log-likelihood distribution of parameters, $\theta=(\alpha_{1},\alpha_{2},b,\nu_{b}),$

[TABLE]

where the values $\langle\log F\rangle_{i}$ and $\sigma_{\log F_{i}}$ and $\nu_{i}$ are log-flux derived by using survival analysis in each log-frequency bin with the main value of log frequency in each bin being $x_{i}$. For the sake of simplicity of implementation, we considered only the errors in normalized log flux in the denominator of Eq. (4) when using the MCMC method.

The final fitting of the broken power-law model is performed using the orthogonal distance regression method (hereafter, ODR), which accounts for the errors on both axes. The $1\sigma$ confidence interval, $\sigma(\nu_{r})$, is computed as $\sigma(\nu_{r}|\theta)=\sqrt{(\nabla_{\theta}\log F)^{T}\Sigma_{\theta}\nabla_{\theta}\log F}$, with $\Sigma_{\theta}$ being the covariance matrix of the fit parameters derived using the ODR method and $\theta$ being the vector of free parameters of the broken power-law.

4.3 Method verification through simulations

Our average SED was derived using survival analysis, where the percentage of upper limits below rest-frame $1\,\mathrm{GHz}$ is high. We tested the applicability of the survival analysis by deriving the radio spectral indices for a mock sample of galaxies. The mock samples were generated with preset spectral indices by simulating different properties of the HSFG sample and of the SED shape.

We constructed the mock sample based on the distribution of individual galaxy redshifts, $z$, a spectral index between $1.4$ and $3\,\mathrm{GHz}$, $\alpha_{2}$, and the log-flux normalization, $b$, in the real HSFG sample. We have chosen to vary the mean value of the low-frequency spectral index, $\alpha_{1}$, to achieve different degrees of flattening in the simulations while keeping the number of galaxies in the simulation equal to the number of galaxies in the real sample. In order to simulate the flattening of the spectral index below the break frequency, $\nu_{b}$, we assumed a broken power-law with a spectral index of $\alpha_{1}$ below $\nu_{b}$, with the following assumptions: 1) $\alpha_{1}$ has the same variance as $\alpha_{2}$ and 2) $\alpha_{1}$ is uncorrelated with $\alpha_{2}$, $z,$ and $b$. Before creating the full mock catalog, we removed the sources in the mock $325\,\mathrm{MHz}$ and $610\,\mathrm{MHz}$ catalogs with fluxes below $5\sigma$, where $\sigma$ is the median RMS of the respective $325\,\mathrm{MHz}$ and $610\,\mathrm{MHz}$ COSMOS maps to simulate non-detections.

Overall, by analyzing the input ($\alpha_{1}^{in}$) and output ($\alpha_{1}^{out}$) spectral index, we find that it is necessary to correct the derived spectral index by $20\%$ (downward) for $\alpha_{1}^{out}\sim 0.5$ and incrementally smaller correction for higher $\alpha_{1}^{out}$. The simulation suggests that the reason for this is the large number of upper limits occurring for lower spectral indices (i.e., flatter input slopes). The left panel of Fig. (4) shows the dependence of the input spectral index of simulations, $\alpha_{1}^{in}$, on the output spectral index of simulations, $\alpha_{1}^{out}$, while the middle panel shows the dependence of the $\alpha_{1}^{in}/\alpha_{1}^{out}$ ratio on $\alpha_{1}^{out}$ Additionally, we show the behavior of the fraction of limits below $1\,\mathrm{GHz}$ as a function of $\alpha_{1}^{in}$ in the right panel of Fig. (4). For $\alpha_{1}^{out}>1$, this correction again rises to $20\%$. Considering the results presented in the next section, this multiplicative correction is $0.8\pm 0.1$ for the real dataset, as indicated by the red points in Fig. (4). Turning the argument around, under the assumption that $\alpha_{1}$ is normally distributed, the percentage of upper limits below $1\,\mathrm{GHz}$ ($88\%$), observed in the real data, indicates $\alpha_{1}\sim 0.4$.

5 Results

Here we present the results of the radio spectral index of the HSFG sample modeled by a broken power-law. In deriving the upper limits, we used the local RMS values from the GMRT maps.

The results of this procedure are shown in the left panel of Fig. (6), yielding a break frequency of $\nu_{b}=4.3\pm 0.6\,\mathrm{GHz}$ (derived by the ODR method). Outside the $2\sigma$ confidence interval (red ellipses in the left panel of Fig. 6), the distribution starts to diverge from the simple form expected by the fit. Furthermore, the MCMC-derived parameter medians in the left panel of Fig. (6) are not aligned with the ODR-derived means and the covariance ellipses derived by ODR fitting diverge from the $\sigma$ contours of MCMC parameter estimates. This divergence probably arises because the SED shape is not equally sensitive to the change of every parameter in the model, resulting in deviations in the MCMC method that are more pronounced outside the $2\sigma$ confidence interval. Therefore, for the final result, we chose to fix the break frequency to the best-fit value provided by the ODR method and use it to fit the SED with a broken power-law with a fixed break frequency. The break frequency derived in this way is fixed to be $\nu_{b}=4.3\,\mathrm{GHz}$ for all following considerations. We fixed the value of $\nu_{b}$ to the ODR-derived best-fitting value and not to the MCMC median because of the MCMC behavior outside the $2\sigma$ confidence interval. We note that these results using break frequency values in the range from $3-6\,\mathrm{GHz}$ are not expected to affect the presented results. This is visible from the wide distribution of $\nu_{b}$ in the lower left plot in the left panel of Fig. (6). As presented in the right panel of Fig. (6), the reduced model behaves according to expectations of a normal distribution.

In Fig. (6) we show the average radio SED of HSFGs and the derived broken power-law SED with fixed $\nu_{b}$. The percentage of upper limits was $88\pm 2\%$, which is consistent with $\alpha_{1}\sim 0.5$ in Fig. (4). The spectral index changes from $\alpha_{1}=0.53\pm 0.04$ below $\nu_{b}$ to $\alpha_{2}=0.94\pm 0.06$ above $\nu_{b}.$ Considering the correction calculated using simulations in Sect. 4.3, the corrected spectral index is $\alpha_{1}^{corr}=0.42\pm 0.06$. Some of the models described in Sect 6.1 can explain this difference, for example, the difference of 0.5 is expected by model 3 assuming continuous injection of electrons (Pacholczyk 1980). Although our data clearly do not follow a simple power-law, we also report the result of the simple power-law fit of $0.67\pm 0.05$ to facilitate a quick comparison to the literature.

In the left and right panels of Fig. (8), we compute mean values of redshift and star-formation rates for different rest-frame frequency bins of our radio SED, respectively. The mean redshifts and SFRs do not vary more than one standard deviation of all the data points ($\sigma_{z}=0.3$, $\sigma_{\log\,\mathrm{SFR}}=0.6$) below $10\,\mathrm{GHz},$ and there is no clear trend between contiguous frequency bins below this frequency. Only a small fraction of the data points, corresponding to $z>2.3$, lie above $10\,\mathrm{GHz},$ and they could not influence the SED fitting procedure significantly. We therefore conclude that our average SED represents the SED of a single population of highly star-forming galaxies (with properties $z=1.7\pm 0.6$ and $\log SFR=2.4\pm 0.3$) over the whole frequency range from $300\,\mathrm{MHz}$ to $10\,\mathrm{GHz}$.

6 Discussion

In Sect. 6.1 we discuss the shape of the radio SED of star-forming galaxies and alternative, physically motivated fitting strategies. In Sect. 6.2 we discuss the impact of the shape of the radio SED on the infrared-radio correlation.

6.1 Shape of the radio SED

In the left panel of Fig. (8), we overplot the different SEDs from the literature and our broken power-law (BPL). Our results, covering a broad range of redshifts, are in overall agreement with literature SEDs derived for (U)LIRGS in the nearby Universe ($z\lessapprox 0.1,\,\,\mathrm{SFR}>10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$) (for galaxies with an $\,\mathrm{SFR}$ in the $50-150\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ range Clemens et al. 2008; Leroy et al. 2011). We find disagreement with the typical model of “normal” star-forming galaxies (NSFGs) above $\sim 10\,\mathrm{GHz}$ (Klein et al. 1988; Condon 1992; Tabatabaei et al. 2017), which may be due to the different thermal fractions in normal and highly star-forming galaxies. In the nearby Universe ($z<0.1$), the typical SED of NSFGs is described by a simple nonthermal power-law with a spectral index of $\sim 0.8$ and a thermal fraction of $8-10\%$ at $\sim 1\,\mathrm{GHz}$ (Klein et al. 1988; Condon & Yin 1990), while studies of HSFGs find a flat spectrum around $\sim 1\,\mathrm{GHz}$ and a steepening spectrum above $10\,\mathrm{GHz}$ (Clemens et al. 2008; Leroy et al. 2011).

For our HSFGs with $\,\mathrm{SFR}>100\,\mathrm{M_{\odot}}/\,\mathrm{yr}$, we find a relatively flat ($\alpha\sim 0.5$) spectral index below $\sim 4\,\mathrm{GHz}$. This result is in line with previous studies of (U)LIRGs (Clemens et al. 2008; Leroy et al. 2011; Galvin et al. 2018). There are multiple effects that could produce a flattening in the spectrum of galaxies in the radio, but only synchrotron aging and free-free absorption play a role in the $1-10\,\mathrm{GHz}$ range. Effects like synchrotron self-absorption and the Razin-Tsytovič effect (Razin, V. A. 1960; Razin, V. N. 1951), an effect that is due to the refractive index of plasma in the radio, can be ruled out because they are usually found at frequencies below $1\,\mathrm{GHz}$ (Deeg et al. 1993). These effects could be even more suppressed, as magnetic fields have been found to be stronger for higher SFRs (Tabatabaei et al. 2017). Furthermore, we do not expect to be able to successfully find evidence for any more complex shape of the SED such as multi-peak behavior of free-free absorption because of the different opacities of different star-forming regions within a galaxy, as any such behavior would be averaged out when constructing an average SED. Tabatabaei et al. (2017) instead pointed out that the nonthermal spectral index (which can be compared to our $\alpha_{1}$) flattens with increasing SFR density, which gives a possible explanation why this effect is pronounced in our HSFG sample.

Bressan et al. (2002) modeled the radio SEDs of dusty starburst galaxies that initially have $\,\mathrm{SFR}$s in the range of $100-1000\,\mathrm{M_{\odot}}/\,\mathrm{yr}$. They concluded that if the SFR changes abruptly, the core-collapse supernova rate will lag behind as a result of the finite lifetimes of stars, leading to a suppressed synchrotron emission in young starbursts. This effect can produce flat spectra below $10\,\mathrm{GHz}$, as the thermal fraction is initially very high. Prouton et al. (2004) analyzed the behavior of the thermal fraction with the age of the starburst for galaxies with an $\,\mathrm{SFR}$s from $30$ to $\sim 150\,\mathrm{M_{\odot}}/\,\mathrm{yr}$. They found that the thermal fraction decreases from the initial value of almost $100\%$ to below $10\%$ at $1.4\,\mathrm{GHz}$ within $\leavevmode\nobreak\ 10^{7.5}\,\mathrm{yr}$.

For the HSFGs in the COSMOS field, we find a steep ($\alpha\sim 1$) spectrum above $4\,\mathrm{GHz}$ that does not exhibit signs of free-free emission. A steep spectral index has also been found by Galvin et al. (2018) ($\alpha=1.06$) for galaxies with an $\,\mathrm{SFR}$ in the range of $\sim 50-250\,\mathrm{M_{\odot}}/\,\mathrm{yr}$, which they associated with a steeper injection index of the CR distribution. Clemens et al. (2008) found a similar trend (for galaxies with an $\,\mathrm{SFR}\sim 50-150\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ ) that they tried to explain using the continuous injection model, which can explain differences in spectral indices up to $\Delta\alpha\approx 0.5$.

Since the SED of HSFGs has been found to differ from a simple power-law, we fit the following models often used in literature (e.g., Calzetti et al. 1994; Condon 1992; Clemens et al. 2008; Condon & Ransom 2016; Lisenfeld et al. 2004; Klein et al. 2018) that might explain this behavior:

Mixed model

• •

This model assumes that the power-law synchrotron spectrum with a spectral index $\alpha$ is attenuated by free-free absorption with opacity $\tau_{\nu}=\frac{\tau_{1}}{\nu^{2.1}}$ arising from the same region. 2. 2.

Foreground screen model

• •

This model assumes that the power-law synchrotron spectrum with a spectral index $\alpha$ is attenuated by free-free absorption with opacity $\tau_{\nu}=\frac{\tau_{1}}{\nu^{2.1}}$ arising from a region that is in our line of sight, but is not cospatial with the emitting region. 3. 3.

Synchrotron aging model

• •

In contrast to the previous two models, this model does not take into account free-free absorption of synchrotron radiation initially described by a simple power-law. This model instead explains the shape of the SED by taking into account that the cosmic ray electron population loses energy with time. The model predicts a spectrum that was initially described by an injection index that produces a power-law SED with a spectral index $\alpha$ that gradually steepens above a frequency $\nu_{b}$ to a spectral index $\alpha+\Delta\alpha$ due to aging. 4. 4.

Synchrotron self-absorption

• •

The synchrotron self-absorption model has the same form as model 1, with the opacity of the form $\tau_{\nu}=(\nu/\nu_{b})^{-(\alpha+5/2)}$, where $\alpha$ is the synchrotron spectral index and $\nu_{b}$ the break frequency. 5. 5.

Lisenfeld et al. (2004) model

• •

In this model, the electrons are propagated through and eventually out of the galactic halo by means of convective winds. 6. 6.

Free-free emission

• •

Models 1-5 can be extended by including an additional contribution of thermal free-free emission, $f_{th}$, as described by Eq. (2). We consider the following:

• •

1. (a)

Synchrotron aging model with free-free emission 2. (b)

Lisenfeld et al. (2004) model with free-free emission

Among the models listed above, free-free absorption and synchrotron aging models would reproduce flattening at lower frequencies, while models (6a)-(6b) can be used to estimate the contribution of free-free emission to the spectrum. The parameterization of the fitted SEDs and their parameters can be found in Table (1). We find that models (1) and (2), which include free-free absorption, fit the low-frequency data well. They overestimate the measured flux at higher frequencies (by $\sim 0.1\,\mathrm{dex}$ at $10\,\mathrm{GHz}$), however. The synchrotron aging model (3) agrees with the broken power-law at both low and high frequencies. The low-frequency spectral index corresponds to the value typically found in supernova remnants (SNRs) and the Galactic plane in the Milky Way (Green 2014; Klein et al. 2018). We tried to constrain the thermal fraction with the broken power-law model by using Eq. (2), but found a strong degeneracy between parameters $f_{th}$ and $\alpha_{2}$ and therefore cannot draw meaningful results from this model. This degeneracy is expected for this class of models (Transtrum et al. 2010) and can only be broken by more observations at even higher frequency ($>10\,\mathrm{GHz}$).

To compare how well these models explain our dataset, we used the log-likelihood test and the Akaike and Bayesian information criteria (AIC and BIC, respectively). The greatest log likelihood means that the broken power-law is optimal in terms of the fitting. In contrast to the log-likelihood test, AIC and BIC criteria favor models with a smaller number of free parameters, penalizing the models with a higher number of free parameters by higher AIC and BIC values. In the right panel of Fig. (8) we show the ranking in our statistical analyses (log likelihood, AIC, and BIC) for the models in Table (1). The results reported in Table (1) were computed using the ODR method restricted to a $2\sigma$ subset of the respective model’s parameter space, derived using the MCMC procedure.

We found that our broken power-law model has the greatest log likelihood and the second lowest AIC and BIC values. The BPL model shares the same average rank of the AIC, BIC, and log-likelihood ranks with model (5) because model (5) has fewer free parameters and is therefore favored by the AIC and BIC values, while it ranks third on the log-likelihood test. As shown in Fig. 15, the broken power-law model and models (3), (5), (6a), and (6b) describe the average radio SED points sufficiently well. We therefore conclude that we can favor models ranking the lowest on all three tests. In addition to the broken power-law, the highest ranking models are model (5) and its extension that includes free-free emission (model 6b). Additionally, literature models of HSFGs (Clemens et al. 2008; Leroy et al. 2011) score better in all three tests than the standard NSFG model (Condon 1992). This comparison means that HSFGs do not have the same SED shape as the NSFGs. The main difference to the NSFGs is the behavior at higher frequencies, which might be due to a smaller thermal fraction at higher SFR. Since our empirical broken power-law describes our dataset better than simple theoretical models (models 1-3), and different goodness-of-fit criteria sort their preference differently, we conclude that neither simple free-free absorption over the whole sample nor synchrotron aging are the only drivers of the radio SED of HSFGs.

6.2 Infrared-radio correlation

Recent stacking and survival analysis studies (Ivison et al. 2010; Magnelli et al. 2015; Delhaize et al. 2017; Calistro Rivera et al. 2017) have found the infrared-radio correlation (described by the $q$ parameter) to decrease with increasing redshift. Delhaize et al. (2017) mentioned the possibility that the computation of rest-frame radio luminosity via a K-correction using a single power-law assumption of the star-forming galaxies’ SED could produce a redshift-dependent trend.

Here we calculate the infrared-radio correlation by using a K-correction based on our best-fitting broken power-law SED for galaxies with an $SFR>10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ and the combination of thermal and nonthermal emission, following Eq. (2), with a nonthermal spectral index of $\alpha=0.97$ and a $10\%$ thermal fraction (Tabatabaei et al. 2017) for $SFR<10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ (described in detail in Appendix B). We adopted a threshold between these two types of SFGs at $10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$, considering that our results are consistent with studies of (U)LIRGs ($SFR>10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ Clemens et al. 2008; Leroy et al. 2011). We find that the broken power-law derived luminosity is $20\%$ lower than that derived by assuming a $0.7$ spectral index for redshifts below $z\sim 1$. We also find that the difference is negligible for sources above $z\gtrsim 2$. We constrained the infrared-radio correlation using the same procedure as applied by Delhaize et al. (2017). In the infrared, we used Herschel-detected $\geq 5\sigma$ objects, with the COSMOS2015 (Laigle et al. 2016) $24\mu m$ catalog as priors to minimize blending, while in the radio, we used the Smolčić et al. (2017b) $3\,\mathrm{GHz}$ Large Project Counterpart catalog. If an object is undetected in the radio, we estimated its upper $q$ limit as five times the local RMS value in the $3\,\mathrm{GHz}$ map convolved to a resolution at which its peak flux ceases to change significantly with increased convolution, while we used lower $q$ limits based on upper limits to the infrared luminosity, derived through integration of the best-fit SED. For details, see Delhaize et al. (2017). In short, we used the joint radio- and infrared-detected samples of star-forming galaxies in the COSMOS field, and then used the survival analysis method to find the mean $q$ value for different redshift bins and fit a power law of the form

[TABLE]

where $q_{0}$ is the value of $q$ derived for $z=0$ and $n$ is the exponent. Contrary to what would be expected from theory (Murphy 2009), but in line with previous studies (Delhaize et al. 2017; Calistro Rivera et al. 2017), the $q$ value, shown in Fig. (9) and computed taking into account the K-corrections described above, decreases with redshift with $n=-0.17\pm 0.02$, and $q_{0}=2.88\pm 0.03$. In Fig. (9) we show the $q-z$ trend derived under various SED assumptions, as described above. In all the cases, we find that $q$ decreases with increasing redshift. We therefore conclude that the shape of the radio SED, as assumed here, probably does not explain the $q-z$ trend. For comparison, in Fig. (9), we also show the resulting $q-z$ trend based on a simple power-law fit with $\alpha=0.79$ for galaxies with an $\,\mathrm{SFR}<10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ (Tabatabaei et al. 2017). This fit is less complex, but has a lower parameter uncertainty. We find that this simplification does not significantly alter the resulting $q-z$ trend ($n=-0.18\pm 0.01$ and $q_{0}=2.89\pm 0.03$.).

Molnár et al. (2018) have studied the impact of morphology on the evolution of the $q$ parameter (for star-forming galaxies covering a broad $\,\mathrm{SFR}$ range of $1-1000\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ and in the redshift range of $0.2-1.5$) and found that spheroid-dominated star-forming galaxies show a steep redshift behavior ($n=-0.19\pm 0.02$), while disk-dominated galaxies do not ($n=-0.04\pm 0.01$). They argued that in order for the radio spectral index of spheroid-dominated galaxies to explain the $q-z$ behavior, it would have to decrease with redshift from $\alpha=0.7$ at $z=0.8$ to $\alpha=0.45$ at $z=1.5$. To investigate this, we have derived an average radio SED for complete samples of star-forming galaxies based on morphology. To match the selection performed by Molnár et al. (2018), we selected SFGs with an $SFR>20\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ to be fairly complete out to a redshift of $z=1.5$, with the requirement that we included only galaxies marked as either spheroid- or disk-dominated, based on the classification of Molnár et al. (2018). We show the position in the $\,\mathrm{SFR}-z$ plane of their respective complete subsets in the left panel of Fig. (10) and their respective broken power-law fits in the right panel. Overall, we find that the median MCMC-derived break frequency is $3\,\mathrm{GHz}$ and that the differences in the parameters of the broken power-law fit between spheroid ($\alpha_{2}=0.9\pm 0.3$, $\alpha_{1}=0.5\pm 0.1$) and disk-dominated ($\alpha_{2}=1\pm 0.2$, $\alpha_{1}=0.51\pm 0.09$) galaxies are smaller than $1\sigma$. This suggests that a difference in radio-SEDs of these morphologically different galaxies, and thus a potentially different K-correction, is not the cause of the differing trend found by Molnár et al. (2018).

7 Summary

We constructed the average radio spectral energy distribution for a $1.4\,\mathrm{GHz}$-selected sample of highly star-forming galaxies with an $\,\mathrm{SFR}>100\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ in the COSMOS field. To achieve a broad rest-frame frequency range, we combined previously published VLA observations at $1.4\,\mathrm{GHz}$ and $3\,\mathrm{GHz}$ with unpublished GMRT observations at $325\,\mathrm{MHz}$ and $610\,\mathrm{MHz}$ by employing survival analysis to account for non-detections in the GMRT maps caused by their higher RMS values. By fitting a broken power-law to the SED, we find the spectral index to change from $\alpha_{1}=0.42\pm 0.06$ below $4.3\,\mathrm{GHz}$ to $\alpha_{2}=0.94\pm 0.06$ above $4.3\,\mathrm{GHz}$. Our results are in line with previous low-redshift studies of star-forming galaxies with an $\,\mathrm{SFR}>10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ that show that their SED differs from the one found in normal star-forming galaxies by having a steeper spectral index around $10\,\mathrm{GHz}$, which could imply a smaller thermal fraction than in normal star-forming galaxies.

We further constructed the IR-radio correlation by using our broken power-law SED for galaxies with an $\,\mathrm{SFR}>10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ and an SED based on a steep nonthermal synchrotron spectral index ($\alpha=0.97$) and a $10\%$ thermal fraction at $1.4\,\mathrm{GHz}$ (Tabatabaei et al. 2017) for galaxies with $\,\mathrm{SFR}<10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$. We find that the shape of the radio-SED is probably not the root cause of the redshift trend found by previous studies, as we see a clear trend of the IR-radio correlation decreasing with increasing redshift also when a more sophisticated K-correction is applied to the data.

Acknowledgements

This research was funded by the European Union’s Seventh Frame-work program under grant agreement 337595 (ERC Starting Grant, ‘CoSMass’).

EV acknowledges funding from the DFG grant BE1837/13-1. EV also acknowledges support of the Collaborative Research Center 956, subproject A1 and C4, funded by the Deutsche Forschungsgemeinschaft (DFG). JD acknowledges financial assistance from the South African SKA Project (SKA SA; www.ska.ac.za).

Appendix A GMRT data reduction

Here we describe the calibration, imaging, source extraction, and validation (including all of the corrections applied) of the GMRT 325 and 610 MHz data.

Observations and imaging

The GMRT observations of the $2\,\mathrm{deg^{2}}$ COSMOS field were conducted using 30 antennas, their longest baseline being $25\,\mathrm{km}$. The channel width of observations was $125\,\mathrm{kHz},$ with a total bandwidth of $32\,\mathrm{MHz}$. The observations were reduced with the source peeling and atmospheric modeling pipeline (SPAM, described in detail in Intema et al. 2017). It accounts for radio interference at low radio frequencies and for direction-dependent effects, such as ionospheric dispersive delay.

The $325\,\mathrm{MHz}$ observations were carried out in a single pointing (see the top middle panel of Fig. 11) under the project 07SCB01 (PI: S. Croft). The observations were allocated $45\,\mathrm{hours}$ in total and comprised four observations with an average time per observation of $\sim 120\min$ and a total time on target of $\sim 40\,\mathrm{hours}$. Data reduction and imaging were performed using the spam pipeline, with the final map being imaged at a resolution of $10.8\times 9.5\,\mathrm{arcsec}^{2}$. A primary beam correction was then performed using the parameters from the GMRT manual. The RMS as a function of fractional area is shown in the top left panel of Fig. (11). We find a median RMS of $97\,\mathrm{\mu Jy}/\mathrm{beam}$ over the $\sim 2\deg^{2}$ COSMOS field, uniformly rising toward the edges of the map.

The $610\,\mathrm{MHz}$ observations were carried out over 19 pointings (see the top middle panel of Fig. 12) under the project 11HRK01 (PI: H. R. Klöckner). The observations were allocated $86\,\mathrm{hours}$ in total and were spread over eight observations with an average time on source of $\sim 4\,\mathrm{hours}$ per pointing. We found that four pointings (labeled 14, 16, 22, and 23) deviated by more than $2\sigma$ from the mean beam axes and beam position angles, and these were therefore excluded. The beam properties of the remaining pointings are shown in the middle row of Fig. (12). After excluding the four pointings, we used one common restoring beam for the remaining pointings. The pointing images were imaged at a common resolution of $5.6\times 3.9\,\mathrm{arcsec}^{2}$. A primary beam correction was then performed on each pointing image before we combined them into one mosaic using the parameters from the GMRT manual, following the procedure outlined in Smolcic et al. (2018). The RMS as a function of fractional area is shown in the top left panel of Fig. (12), yielding a median RMS of $39\,\mathrm{\mu Jy}/\mathrm{beam}$ over the $2\deg^{2}$ COSMOS field.

A summary of the GMRT data reduction statistics can be found in Table (2)333The final catalogs will be publicly available on IRSA.. Below we describe the source extraction method, reliability tests, and corrections we applied to the catalogs.

Source extraction

Blobcat (Hales et al. 2012) was used to extract sources down to five times the local RMS value ($5\sigma$) in each map. The local RMS values were taken from the RMS map produced by the aips task rmsd. We used default blobcat input parameters except for the minimum blob size, which was set to $3$ pixels in right ascension and declination, as described in Smolčić et al. (2017a). Even though the imaged pointing extends beyond the COSMOS field, we restricted the catalog and further analysis to the $2\,\mathrm{deg^{2}}$ COSMOS field.

We find 633 source components in the $325\,\mathrm{MHz}$ map and 999 source components in the $610\,\mathrm{MHz}$ map. After accounting for multicomponent sources and blending, the number of sources in the COSMOS field is 633 (see below for details) in the $325\,\mathrm{MHz}$ map, while the number of sources remaining in the $610\,\mathrm{MHz}$ map is 986. All of the sources detected by blobcat within the $2\,\mathrm{deg^{2}}$ COSMOS field are shown in the top middle panels of Figs. (11) and (12). Below we describe the various tests and corrections performed to generate the final source catalogs.

Bandwidth smearing

Given the finite bandwidth of the map, a certain amount of bandwidth smearing is expected. To assess the effect of bandwidth smearing on the extracted fluxes, we computed the peak-over-total flux ratio, $S_{P}/S_{T}$.

For the $325\,\mathrm{MHz}$ map, which consists of only one observed pointing, in the middle left panel of Fig. (11) we show the $S_{P}/S_{T}$ ratio as a function of the distance to the pointing center. Bandwidth smearing reduces the value of peak fluxes compared to the total fluxes with increasing distance to the pointing center (Bridle & Schwab 1999). We binned the data, calculated the median $S_{P}/S_{T}$ , and used the expected functional dependence of the flux ratio and distance to the pointing center in the presence of bandwidth smearing (as explained in detail in Bridle & Schwab 1999) to account for bandwidth smearing of the peak fluxes. We used this relation to correct the peak flux for each source, and we verified that the $S_{P}/S_{T}$ ratio does not depend on distance to the pointing center after the correction. The final total fluxes are consistent with the shallower Smolčić et al. (2014) $324\,\mathrm{MHz}$ Catalog (top right panel of Fig. 11).

In the top right panel of Fig. (12), we show the peak-over-total flux ratio, $S_{P}/S_{T}$, as a function of the S/N of sources extracted from the $610\,\mathrm{MHz}$ mosaic, combined from 19 pointings. We corrected for bandwidth smearing in the $610\,\mathrm{MHz}$ mosaic by dividing the peak flux of each source with the median value of the peak-over-total flux distribution.

Resolved and unresolved sources

In radio maps, ideal unresolved sources with infinitely high S/Ns have peak fluxes in $\mathrm{Jy/beam}$ equal to their total fluxes because their flux is contained within the area that is covered by the beam. However, because of the noise properties of the map, we expect the $S_{T}/S_{P}$ ratio for unresolved sources to increasingly scatter around unity with decreasing S/N. We thus discerned resolved sources in the $325\,\mathrm{MHz}$ map and the $610\,\mathrm{MHz}$ mosaic through the following procedure. We first fit the fifth percentile in the total-over-peak ratios that are lower than 1 as a function of the S/N, as the total flux can only be lower than 1 as a result of statistical effects. We did this by fitting all detections with a total flux lower than the peak flux with a relation of the form $S_{T}/S_{P}=1-a(S/N)^{b}$, as is commonly used in the literature (see, e.g., Smolčić et al. 2017a; Bondi et al. 2003). Second, we assumed that the noise is symmetric around unity. We therefore mirrored this envelope around $S_{T}/S_{P}=1$ and considered all sources above it resolved.

The best-fit parameters are $a=3.58$ and $b=-1.13$ and $a=1.26$ and $b=-0.76$ for the $325\,\mathrm{MHz}$and $610\,\mathrm{MHz}$ catalogs, respectively. Therefore, we consider 177 (196) sources with $S_{T}/S_{P}>1+3.58(S/N)^{-1.13}$ ($S_{T}/S_{P}>1+1.26(S/N)^{-0.76}$) to be resolved at $325\,\mathrm{MHz}$ ($610\,\mathrm{MHz}$), as shown in the middle right (left) panel of Fig. (11) (Fig. 12). We set their total fluxes of unresolved sources to the values of their respective peak fluxes. We note that the fraction of resolved sources (see Table 2) is higher in the $610\,\mathrm{MHz}$ catalog because the $610\,\mathrm{MHz}$ map is $1.3-1.8$ times more sensitive than the $325\,\mathrm{MHz}$ map (derived by assuming spectral indices in the range from 1 to 0.5).

Multicomponent sources and deblending

If the resolution of the map is good enough, source extraction algorithms start to separate jets, lobes, and star-forming regions of radio sources into separate sources. We therefore need to associate these components with a single source, termed multicomponent. To quantify this effect in our GMRT maps and find multicomponent sources, we made use of the $3\,\mathrm{GHz}$ source catalog. This catalog was used because it reaches a much higher depth (it is $4-13$ times more sensitive than the $325\,\mathrm{MHz}$ map and $3-7$ times more sensitive than the $610\,\mathrm{MHz}$ map for spectral indices of 1-0.5) and a higher resolution ($0.75\arcsec$ at $3\,\mathrm{GHz}$ vs $\sim 10\arcsec$ at $325\,\mathrm{MHz}$ and $\sim 5.6\arcsec$ at $610\,\mathrm{MHz}$). If the source has multiple components, that is, components that are not connected by detectable emission, in this higher resolution $3\,\mathrm{GHz}$ map (see Smolčić et al. 2017a, for details), it is a good candidate for a multicomponent source in the map of lower resolution. Conversely, we do not expect single-component sources at $3\,\mathrm{GHz}$ to appear as multicomponent at lower resolution.

By visually inspecting cutouts of the maps extracted at positions of multicomponent sources in the $3\,\mathrm{GHz}$ source catalog, we found components that needed to be combined into one multicomponent source. The method used to detect multicomponent sources was independently verified using the algorithm described in detail in Vardoulaki (2018 in review). This algorithm generated $3\,\mathrm{GHz}$ cutouts that were convolved to a $10\arcsec$ resolution, and we then used them to create masks for the component-finding algorithm. The algorithm searches for positions within the mask and within a radius that is defined by the linear projected size of the source (linear projected sizes are measured by a semi-automatic machine-learning technique and have been verified by eye, Vardoulaki 2018 in prep.). To be consistent with Smolčić et al. (2017a), the total flux was measured after the multicomponent source was identified using the AIPS task TVSTAT after clipping all the pixels below $2\sigma$. This is consistent with the method used for all other radio source catalogs in the COSMOS field (Schinnerer et al. 2007, 2010; Smolčić et al. 2014, 2017a). In the $325\,\mathrm{MHz}$ map, we find only 2 components that needed to be combined into 1 multicomponent source, and in the $610\,\mathrm{MHz}$ mosaic we find 25 components that needed to be combined into 12 multicomponent sources. As our analysis is based on the $3\,\mathrm{GHz}$ source catalog, in the final catalog, we set the position of the multicomponent source in the GMRT catalogs to its respective position in the $3\,\mathrm{GHz}$ catalog, its peak flux to -99, and we replaced its total flux by the value derived above and set the multicomponent flag to 1.

We further identified two sources that were blended into one source at $325\,\mathrm{MHz}$. This number of blended sources was expected, based on the analysis of the $3\,\mathrm{GHz}$ source catalog (middle right panel of Fig. 11). We deblend the identified blended source in the $325\,\mathrm{MHz}$ map by using blobcat to be consistent with other flux estimates of the catalog. To ensure that blobcat separated the blended source into two sources, we set the $fS/N$ parameter to 5.78. We tested these measured fluxes using the CASA task IMFIT with visually identified masks that separated the two sources. The sum of total fluxes of the two sources derived using blobcat equals the total flux of the blended source. In the final catalog, we replaced the blended source with two entries in the catalog and the blobcat -derived properties of the respective sources.

To summarize, we associated two source components with one source and deblended one source into two separate sources in the $325\,\mathrm{MHz}$ catalog. We also associated 25 source components with 12 multicomponent sources in the $610\,\mathrm{MHz}$ catalog. The total number of sources in the GMRT catalogs of the COSMOS field is therefore 633 and 986, respectively.

False-detection rate

As noise in radio-interferometric images is non-Gaussian and source detection is based on a $5\sigma$-base peak detection, we statistically expect that a small percentage of detections might be random noise fluctuations. To assess the false-detection rate, we assumed that the background noise is symmetric around zero and ran blobcat on the inverted maps (i.e., map multiplied by -1), with the same input parameters as in the original maps. Given that sources with negative fluxes are not expected in the radio map, any source extracted from the inverted map can be reliably identified as a false detection. We find 19 “detections” in the inverted $325\,\mathrm{MHz}$ map, yielding a false-detection rate of $3\%$ and 13 “detections” in the inverted $610\,\mathrm{MHz}$ map, yielding a false-detection rate of $1.3\%$.

Astrometry

To test the astrometric accuracy of the data, we show in the bottom row panels of Figs. (11) and (12) the offsets of the source positions in the $325\,\mathrm{MHz}$ and $610\,\mathrm{MHz}$ catalogs, respectively, from those present in the $3\,\mathrm{GHz}$ catalog. The mean offsets in right ascension are within one tenth of the beam size of the GMRT maps, which means that the astrometric accuracy is good.

Final catalogs

In the final catalog we report the positions of sources as well as their peak (in $\mathrm{Jy/beam}$), and total fluxes (in $\mathrm{Jy}$). For the multicomponent sources, the total flux of both components is indicated in one catalog entry. The ’Multicomponent’ column indicates the presence of more than one component to the source. For resolved sources, the column ’Resolved’ is set to 1.

Appendix B K-corrections

Generally, the luminosity at a rest-frame frequency $\nu_{e_{2}}$, $L(\nu_{e_{2}})$, is connected to the flux measured at the observed frequency $\nu_{o_{2}}$, $S(\nu_{o_{2}})$, as

[TABLE]

where $z$ is the object redshift and the observed frequency is connected to the emitted frequency by $\nu_{e}=\nu_{o}(1+z).$ This relation does not depend explicitly on the shape of the SED, so we have the relation connecting ratios of luminosities at two different frequencies

[TABLE]

We can model $\frac{S(\nu_{e_{2}})}{S_{(}\nu_{e_{1}})}$ by an SED. Inserting Eq. (7) into the first equation, we obtain

[TABLE]

We can now proceed to find the K-correction for a broken power-law.

B.1 Broken power-law

We assumed that the spectrum is that of a broken power-law with spectral indices $\alpha_{0}$, $\alpha_{1}$, and so on, with each new spectral index, $\alpha_{k}$, beginning above a fixed frequency $\nu_{k}$. In this case, the flux for each frequency range $k$ can be expressed as

[TABLE]

where $b_{k}$ is the $k$-th bin normalization constant. Using the fact that the flux needs to be a continuous function of $\nu$, we can roll back $b_{k}$ to the first bin, $b_{0}$ as

[TABLE]

The bin number, $k$, is easily determined from $\nu$ and the list of break frequencies $\nu_{i}$, giving functions $k(\nu)$ and $\alpha(\nu)$. We can repeat this process for each frequency $\nu_{e_{1}}$ and $\nu_{e_{2}}$, yielding

[TABLE]

Finally, we obtain the luminosity as

[TABLE]

For the broken power-law with a single break frequency, $\nu_{e_{2}}=1.4\,\mathrm{GHz}$ and $\nu_{o_{1}}=3\,\mathrm{GHz}$, we can express Eq. (12) as

[TABLE]

The left panel of Fig. (13) shows the luminosity at $1.4\,\mathrm{GHz}$ rest-frame derived using various K-corrections. From the right panel of the same figure, we see that luminosities are by up to $26\%$ overestimated compared to the broken power-law for $z\leq 1$ if we assume a single spectral index of $0.7$.

B.2 K-corrections and thermal fraction

We used the results of Tabatabaei et al. (2017) to derive luminosities of normal star-forming galaxies. They have found that the SED of normal star-forming galaxies has a steep nonthermal spectral index ($\alpha_{2}=0.97$) and a $\sim 10\%$ thermal fraction at $1.4\,\mathrm{GHz}$. The K-correction in this case is

[TABLE]

As discussed in Sect. 6.2, we compared the $q-z$ trends derived using this K-correction and also the K-correction based on the Tabatabaei et al. (2017) simple power-law spectral index of $0.79$. In the next section, we derive the average correction of radio luminosity needed for a sample of galaxies consisting of both NSFGs and HSFGs.

B.3 Average K-corrections

To better understand the impact of the different K-corrections on the $q-z$ trend, we also derived corrections of $L_{1.4\,\mathrm{GHz}}$ for an average $1.4\,\mathrm{GHz}$ luminosity by assuming a distribution of fluxes, $f(\,\mathrm{SFR},z)$, that is dependent on both $\,\mathrm{SFR}$ and $z$. We started by rewriting Eq. (12) for SFGs as $L_{1.4\,\mathrm{GHz}}=K(\,\mathrm{SFR},z)S_{3\,\mathrm{GHz}}$, where $S_{3\,\mathrm{GHz}}$ is the observed flux density. As discussed in Sect. 6.2, we used a K-correction, $K_{BPL}$, based on the broken power-law SED for galaxies with an $\,\mathrm{SFR}>10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ and the Tabatabaei et al. (2017) SED-based K-correction, $K_{Tab}$, for $\,\mathrm{SFR}<10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$. The average $1.4\,\mathrm{GHz}$ luminosity at redshift $z$ is

[TABLE]

We divided this average luminosity by the average luminosity, $\langle L_{1.4\,\mathrm{GHz}}\rangle(z)_{\alpha}$, derived by assuming a K-correction, $K_{\alpha}$, based on a simple power-law SED with a single spectral index $\alpha=0.7$, yielding the correction of $L_{1.4\,\mathrm{GHz}}$,

[TABLE]

where we have introduced the fraction of sources with an SFR below $10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$, $p_{<10}$, defined as

[TABLE]

Given the finite number of sources in our sample, this fraction is estimated as the sum of fluxes of all sources within a redshift bin that have an $SFR<10\,\mathrm{M_{\odot}}/\,\mathrm{yr}$ divided by the sum of fluxes of all sources within the redshift bin. We see that $C_{L_{1.4\,\mathrm{GHz}}}(p_{<10},z)$ is a linear interpolation between K-corrections derived by assuming a broken power-law SED and the Tabatabaei et al. (2017) SED and can be useful for quickly estimating the impact of using a simple K-correction based on an $\alpha=0.7$ power-law SED instead of a more elaborate SED. Average corrections for selected redshifts and $p_{<10}$ values can be found in Table (3). In Fig. (14) we calculated average corrections for a subset of galaxies from the (Delhaize et al. 2017) sample that were detected in the VLA-COSMOS $3\,\mathrm{GHz}$ Counterpart catalog (Smolčić et al. 2017b). The corrections are significant only below $z\sim 1$. The fact that the corrections are not needed for higher redshift bins indicates that the shape of the radio SED cannot explain the $q-z$ relationship.

Appendix C Radio SED of HSFGs with models from Table 1