© ESO 2021

1 Introduction

Bright submillimetre galaxies (SMGs; Smail et al. 1997; Hughes et al. 1998) are prime laboratories for investigating the physical processes involved in high redshift star formation. With total luminosities that imply star-formation rates (SFRs) of hundreds to thousands of solar masses per year (e.g. Barger et al. 1998; Michałowski et al. 2017) – relative to the Milky Way’s ~ 1 M_⊙ yr⁻¹ (e.g. Chomiuk & Povich 2011) – these sources are the locations of the most extreme star formation in the Universe. Understanding the nature of SMGs and the star formation mechanisms at play within them is an important step in beginning to address several open questions in the formation and evolution of massive galaxies.

In the local Universe, galaxies with very high infrared (IR) luminosities (L_IR> 10¹¹ L_⊙ and L_IR > 10¹² L_⊙ – luminous and ultraluminous infrared galaxies, LIRGs and ULIRGs, respectively) and therefore high inferred SFRs (> 10−100 M_⊙ yr⁻¹) are rare and usually the result of recent major mergers (Sanders & Mirabel 1996; Rigopoulou et al. 1999; Hopkins et al. 2013). This fits into an evolutionary picture in which mergers trigger episodes of extreme star formation, which is subsequently quenched by, for example, feedback from an active galactic nucleus (AGN), leaving a ‘red and dead’ massive elliptical galaxy (e.g. Hopkins et al. 2006). This picture can be expanded to describe more comprehensively the diversity of star-forminggalaxy populations and physical processes observed in galaxy evolution across a range of redshifts, with secularevolutionary processes also playing a role as star-formation triggers, and including galaxies with more extended, clumpy star-forming regions. While evidence from gas kinematics suggests that at high redshift (z > 1) some episodes of high SFRs in galaxies may also be triggered by mergers (e.g. Tacconi et al. 2008), the average SFR of the whole galaxy population is significantly higher than at low redshift (e.g. Madau & Dickinson 2014). Not only are galaxies with SFRs comparable to local LIRGs and ULIRGs the dominant galaxy population (Chary & Elbaz 2001; Le Floc’h et al. 2005; Magnelli et al. 2013; Bothwell et al. 2013), but a population of rare galaxies with extremely high SFRs appears at the same epoch. The fact that such intensely star-forming galaxies are common at high redshift (z > 1) suggests they may play a crucial role in galaxy evolution by rapidly assembling the high stellar masses required for the formation of the most massive galaxies we see in the nearby Universe. In this evolutionary picture, SMGs at high redshift (z > 2) are attractive candidates for the progenitors of massive elliptical galaxies at z = 0 (Lilly et al. 1999; Tacconi et al. 2008; Hickox et al. 2012; Toft et al. 2014, but see Garcia-Vergara et al. 2020). Studying how this extreme star-forming phase shapes galaxy populations at high redshift is therefore crucial in uncovering the origins and nature of massive galaxy populations in the local Universe.

To gain an understanding of the mechanisms driving star formation in these galaxies, we require observations spanning a range of wavelengths, from ultraviolet (UV) to radio, to trace physical processes occurring at different energy scales. Star-forming galaxies can emit prodigious amounts of energy from the far-infrared (FIR) to millimetre range due to the thermal heating of dust by UV photons from massive young stars (Kennicutt 1998). This re-processed emission is a robust tracer of star-formation activity in star-forming galaxies, and observations of the shape of the FIR spectral energy distribution (SED) have long been used to characterise SFRs in galaxies (e.g. Sauvage & Thuan 1992; Devereux & Hameed 1997; Chapman et al. 2010; Magnelli et al. 2012).

While much of the UV and optical emission is obscured in the most dusty sources, star formation may also be traced using radio continuum emission, which is not obscured by dust. The radio spectrum arising from star formation consists of two components: thermal, free–free emission from H II regions and non-thermal synchrotron emission. The contribution of free–free emission is only significant at high frequencies (> 10 GHz), and so in this study we focus on radio synchrotron emission, which dominates the spectrum at the low frequencies considered here. This synchrotron continuum results from supernovae exploding after a delay of several megayears following the births of populations of O and B stars, which produce cosmic rays that interact with the galaxy’s magnetic field. Since both thermal FIR emission and non-thermal radio continuum emission trace physical processes associated with star formation, one would expect these quantities to correlate. Indeed there is a well-known tight correlation between the FIR and radio luminosities of star-forming sources, the far-infrared to radio correlation (FIRC; van der Kruit 1971; de Jong et al. 1985; Helou et al. 1985; Ivison et al. 2010), which is observed consistently across over four orders of magnitude of galaxy luminosities (Yun et al. 2001).

The synchrotron radio emission of star-forming galaxies is observed to follow approximate power law behaviour at gigahertz frequencies, with a typical spectral index of α ≃−0.7 (Mauch et al. 2013)¹. The sources, acceleration, and propagation mechanisms of the cosmic ray electrons responsible for the radio continuum emission of star-forming galaxies are still not well understood, but there are several large-scale effects that can be observed to impact the shape of the radio spectrum due to conditions of the interstellar medium (ISM). Over time, a synchrotron spectrum will steepen at high frequencies due to radiative losses as high-energy electrons lose energy more rapidly than low-energy electrons (Kardashev 1962; Scheuer & Williams 1968). In the low-frequency regime, free–free absorption can flatten the spectrum below a turnover frequency at which the ISM becomes optically thick (Condon 1992; Clemens et al. 2010; Lacki 2013). Therefore, the shape of the radio synchrotron spectrum and its divergence from a simple power law hold clues as to the conditions of the ISM in and around star-forming regions.

Single-dish submillimetre observations are currently limited by the resolving power of relatively small telescopes operating at such long wavelengths, resulting in much lower-resolution observations than, for example, optical imaging. However, observations at these wavelengths also benefit from a negative K-correction (Blain & Longair 1993): as galaxies are redshifted, the observed-frame 850 μm emission traces an increasingly bright part of the galaxy’s FIR to millimetre spectrum due to the shape of the Rayleigh-Jeans tail, in effect compensating for cosmological dimming. This results in the observed submillimetre flux density remaining roughly constant for a galaxy of a given luminosity between 0.5 <z < 10 at 850 μm, enabling the detection of submillimetre sources out to very high redshifts. Both optical and radio wavelengths, however, suffer a positive K-correction across the majority of the spectrum, meaning that the flux density of sources of the same luminosity decreases with distance. Many studies have identified SMGs without counterparts at optical to near-infrared wavelengths (Simpson et al. 2014; Wang et al. 2019b; Dudzevičiūtė et al. 2020). Radio observations provide a view of star formation that is not biased by dust; however, due to the positive K-correction, observing galaxies in the radio at the high redshifts at which SMGs are most numerous (z > 2) is very challenging.

To obtain a more complete picture of the physical processes that shape star formation in the early Universe, very deep radio surveys overwide areas of sky are required to complement deep submillimetre surveys. Previous work has used high-resolution radio observations, typically at 1.4 GHz, as a method of pinpointing the position of submillimetre sources detected in single-dish surveys (Ivison et al. 2002; Chapman et al. 2005); however, such work has been limited by the depth of available radio sky surveys,with dedicated deep surveys over only small regions of sky, resulting in a view biased towards the brighter radio sources. Studies have largely included limited radio spectral coverage of SMGs, focusing on the nature of the FIRC and its relation to properties such as stellar mass and redshift (e.g. Yun et al. 2001; Ivison et al. 2010; Smith et al. 2014). The Low Frequency Array (LOFAR; van Haarlem et al. 2013) has opened up new ways of studying galaxies in the radio, and a number of studies have used LOFAR’s capabilities to investigate this relationship between star formation and radio luminosity in the low-frequency regime – for example Gürkan et al. (2018), Read et al. (2018), Smith et al. (2021), and Wang et al. (2019a). However, these studies generally investigate the statistical properties of large samples of galaxies, in optically selected samples at low redshift (z ≲2), rather than probing the shapes of individual radio spectra. Thomson et al. (2019) conducted an in-depth study of high-frequency (>610 GHz) spectral curvature in SMGs, finding evidence of curved spectra that they attributed to spectral ageing of the synchrotron emission from star formation; their results implied estimated starburst ages consistent with expected SMG lifetimes. Studies at low frequencies, where we may be able to observe absorption processes affecting the shape of the spectrum, have been hampered by a lack of sufficiently deep, wide-area data. More comprehensive observations of the shape of the radio spectrum, extending to lower frequencies, can provide a probe of the physical conditions that give rise to extreme star formation in SMGs. Calistro Rivera et al. (2017) exploit LOFAR’s frequency range to investigate the spectral shapes of star-forming galaxies and AGN, finding evidence of low-frequency spectral flattening in the star-forming sample. This sample is also constrained in redshift, focusing on local galaxies rather than the peak of star formation at z > 2, and so does not probe the bulk of the highly star-forming SMG population. Chyży et al. (2018) also find weak spectral flattening in local star-forming galaxies with LOFAR but largely attribute this slight effect to synchrotron losses, predicting stronger low-frequency spectral flattening due to free–free absorption at high redshift, where galaxies with high SFRs are more common.

In this study, we make use of new LOFAR deep field observations. Reaching ~ 22 μJy beam⁻¹, these observations have the potential to reveal the faint radio counterparts to high-redshift submillimetre sources at low frequencies. We select a sample of SMGs using observations from the SCUBA-2 Cosmology Legacy Survey (S2CLS; Geach et al. 2017), currently the largest area sky survey of its kind, which allows us to limit our study to the sites of the most extreme star formation. Selecting sources from the Lockman Hole field, for which we have survey coverage with both LOFAR and S2CLS, we characterise their radio spectra with additional radio data from the Jansky Very Large Array (JVLA) and the Giant Metrewave Radio Telescope (GMRT). We describe the data in Sect. 2, our sample selection in Sect. 3.1, and how we measure radio fluxes in Sect. 3.2. We briefly comment on two bright S2CLS sources that are undetected at every other wavelength in Sect. 3.3. In Sects. 3.4 onwards, we focus on the radio spectra in more detail. We investigate the diversity of radio spectral shapes and luminosities exhibited by sources in Sect. 4.

Throughout this paper we assume a flat Lambda cold dark matter (ΛCDM) cosmology with H₀ = 69.3 km s⁻¹ Mpc⁻¹ and Ω_m = 0.287 (Hinshaw et al. 2013).

2 Data

2.1 S2CLS

The S2CLS observed approximately 5 square degrees of extragalactic sky across several well-studied fields at 850 μm to a depth of ~ 1 mJy beam⁻¹, close to the SCUBA-2 confusion limit. In this work we focus on the Lockman Hole North field, centred at (α, δ) = 10^h 46^m07^s, +59°01′17″. The mapping strategy resulted in an approximately circular map of diameter 30′, with nearly uniform noise coverage over 0.28 square degrees, at an rms depth of 1.1 mJy beam⁻¹. This resultsin 126 submillimetre sources detected at a significance of > 4. Full details of the SCUBA-2 data reduction, catalogue, and source statistics are given by Geach et al. (2017). In this work we use the S2CLS source catalogue, and throughout we employ the de-boosted 850 μm flux densities.

2.2 LOFAR

We used the deep Lockman Hole image described by Tasse et al. (2021), which is based on 112 hours of LOFAR observations. The image has a central rms noise level of 22 μJy beam⁻¹ at a central frequency of 144 MHz and a resolution of 6 arcsec. Together with the images in the Boötes and European Large Area Infrared Space Observatory Survey-North 1 (ELAIS-N1) fields (not covered by S2CLS) this comprises the first data release of the LOFAR Two-metre Sky Survey (LoTSS) Deep Fields and is one of the deepest images ever made at this frequency. It offers us the best opportunity yet available to study the low-frequency properties of distant submillimetre sources.

Of all sources in the LOFAR catalogue, ~98 per cent have a candidate optical identification, selected by a combination of likelihood ratio and visual inspection using new matched-aperture multi-wavelength catalogues, as described by Kondapally et al. (2021). For the Lockman Hole, optical data are provided by the Spitzer Adaptation of the Red-Sequence Cluster Survey (SpARCS) and The Red Cluster Sequence Lensing Survey (RCSLenS) with the Canada France Hawaii Telescope (CFHT; Wilson et al. 2009; Muzzin et al. 2009; Hildebrandt et al. 2016), and there are near-infrared data from the UKIRT Infrared Deep Sky Survey - Deep Extragalactic Survey (UKIDSS DXS; Lawrence et al. 2007) as well as mid-infrared (MIR) data from the Spitzer Wide-Area Infrared Extragalactic Survey (SWIRE) and the Spitzer Extragalactic Representative Volume Survey (SERVS; Lonsdale et al. 2003; Mauduit et al. 2012). Additional FIR data come from the Spitzer Multiband Imaging Photometer (MIPS) and the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012). The Herschel catalogues use the optical, IR, or radio positions as a prior to obtain de-convolved flux densities for blended sources. (McCheyne et al., in prep.).

Spectroscopic redshifts were used where available, but the majority of redshifts were generated photometrically from the optical through MIR data in the manner described by Duncan et al. (2021). The broad spectral coverage in this range allows reasonable estimates of photometric redshift to be made (see Table 1).

2.3 Additional radio data

We supplemented the LOFAR catalogue with deep archival JVLA observations at 324 MHz and 1.4 GHz (Owen et al. 2009; Ivison et al. 2002), and at 610 MHz from the GMRT (Ibar et al. 2009). The Very Large Array (VLA) observations reach central rms noise levels of ~70 μJy beam⁻¹ and ~ 4.8 μJy beam⁻¹ at 324 MHz and 1.4 GHz, respectively, with resolutions of 6 arcsec and 1.4 arcsec. Observations from Ibar et al. (2009) at 610 MHz reach a central rms noise level of ~14 μJy beam⁻¹, with a resolution of ~6 arcsec. These survey footprints cover the S2CLS Lockman Hole coverage in its entirety. Combined, these observations are the deepest available of the Lockman Hole field across the radio spectrum, and among the deepest radio observations to date for any extragalactic survey field.

Fig. 1 Cutouts of an example S2CLS source (ID 1 in our numbering system, an 11.91 mJy 850 μm source) in each radio frequency used in this study. Each square is 50 arcsec across, with the approximate S2CLS beam size (~ 15 arcsec FWHM) marked with an orange circle. Image contrast is scaled arbitrarily for clarity.

3 Sample selection and FIR properties

3.1 Sample

Our sample consists of the 53 point sources detected at >5σ – at which significance the false detection rate falls below 1 per cent – at 850 μm in the S2CLS Lockman Hole North field, with a median flux density of S₈₅₀ = 6.45 mJy (details in Tables 1 and 2). For each submillimetre source we extracted a thumbnail image cutout in the SCUBA-2 map and at the same position in each radio map. Figure 1shows an example source with cutouts at all four radio frequencies. Due to the large beam size of SCUBA-2 (~ 15 arcsec full width at half maximum, FWHM), there is a risk of blending, where submillimetre flux from several galaxies that is unresolved within the large beam contributes to the source flux measurement (e.g. Hayward et al. 2018). We checked that each S2CLS source corresponds to a single point source in the high-resolution radio images via a visual inspection of the cutout images. By doing so,we were able to constrain any possible multiplicity to within the 1.4 arcsec resolution of the 1.4 GHz images. There are a number of sources that split up into multiple components in the 1.4 GHz images, and we discuss our treatment of them below.

Given the 6 arcsec resolution of the LOFAR image, we ran a simple positional cross-match between each submillimetre source and the LOFAR catalogue, identifying the closest LOFAR source within a generous 15 arcsec radius, equivalent to the size of the SCUBA-2 beam, and approximately corresponding to the 5σ positional error on SCUBA-2 sources at the limit of the S2CLS catalogue (Geach et al. 2017). This results in 44 matched LOFAR sources and nine sources for which there is no clear LOFAR counterpart. We calculated the corrected Poissonian probability, p, of serendipitous LOFAR matches within the search radius following Biggs et al. (2011) and using the LOFAR source surface density calculated above the detection threshold of the LOFAR catalogue sources. For identifications with p < 0.05, we assumedthat this is a robust cross-match. Three identifications were found to have p > 0.05 and were excluded, reducing the final sample size of robust S2CLS–LOFAR matches to 41 sources. Only three of these sources (7 per cent) have spectroscopic redshifts available in the LOFAR catalogue, and we used photometric redshifts where spectroscopic were not available.

Ten of the 53 objects in the original sample break up into multiple sources in the LOFAR image such that it is unclear whether the closest positional match corresponds to the correct counterpart. This gives a multiplicity fraction of ~ 20 per cent. Estimates of multiplicity based on high-resolution Atacama Large Millimeter/submillimeter Array (ALMA) follow-up of previous single-dish submillimetre surveys range from 15−40 per cent (e.g. Barger et al. 2012; Chen et al. 2013; Hodge et al. 2013; Michałowski et al. 2017), so our observed multiplicity is consistent with previous observations. It is possible that we missed multiples in the case that their angular separation is smaller than the resolution of our radio imaging; however, this is likely to affect only a small number of sources in our sample. Of those that are visibly multiples within the S2CLS beam based on the radio images, we identified cases in which the secondary source is likely a foreground contaminant by visually inspecting cutouts of the source at radio, IR, and optical wavelengths. We subsequently treated these as we did the single sources (see Sect. 3.2), ensuring that the identified counterpart is identified in source extraction. This does neglect the case in which several sources contribute submillimetre flux and only one of which has a detectable radio counterpart, a possibility that could arise if, for example, a high redshift submillimetre source is by chance aligned with a low redshift source. However, our final selected sample has a median redshift of z = 2.61 (see discussion in Sect. 4 for more details of the redshift distribution), which is consistent with SMG populations (Chapman et al. 2005; Casey et al. 2014, Sect. 4 and references therein), and so it is unlikely that low redshift interlopers contaminate our sample significantly. There is also the possibility that the submillimetre fluxes are boosted by low luminosity sources that do not have significant radio emission – in this case, for any significant contribution we would expect a detectable radio source given the FIRC, so this also seems unlikely. Without high-resolution submillimetre imaging we are unable to distinguish between a single source and the above cases, and so we continue with the assumption that in all bar one of our sources we may attribute the observed submillimetre flux to a single radio source in the LOFAR catalogue.

There are two notable sources – sources 8 and 16 – that are bright S2CLS detections but do not have significantly detected counterparts in any of the radio, optical, or IR images. We speculate that these may be very high redshift (z > 4) galaxies and discuss them in more detail in Sect. 3.3.

3.2 Radio fluxes

To measure the source flux density at the four radio frequencies, we used the Source Extraction and Photometry (SEP; Barbary 2016) Python package, an application of SExtractor (Bertin & Arnouts 1996) for Python. We made the assumption that all sources are point-like across all radio frequencies, given that inspection of the sources at the highest angular resolution (1.4 arcesc at 1.4 GHz) reveals no indication of resolved structure across the sample.

Using SEP, we measured flux densities in each of the four radio frequencies for the 41 LOFAR-detected sources described in the previoussection. These flux densities are consistent with those presented in the LOFAR catalogue. We used the LOFAR catalogue coordinates from the positional cross-match with the S2CLS source list as described above as the source coordinates, given the higher angular resolution of LOFAR. At each radio frequency, we examined the source identification resulting from SEP and ensured that the brightest source within the SEP detection ellipse is the expected counterpart to the LOFAR position so that any multiples are not incorrectly identified as the counterpart. We then took the peak pixel value from within the SEP source ellipse to bethe source flux density. Uncertainties on these measured flux densities were calculated using the off-source pixel-to-pixel rms.

Of the ten sources that break up into multiple components at higher resolution, several contain sources that are likely to be foreground galaxies based on their bright optical luminosities. Comparing luminosities across the full range ofmulti-wavelength images, we were able to determine the high redshift counterpart in nine of these ten images. In these cases, the source SMG can be identified in the LOFAR image, and thus we analysed these as describedabove, using the LOFAR coordinates for reference. In the one case where there is no clear single counterpart and there may be truly blended emission from several co-located galaxies, further analysis would require a method of partitioning fluxes that we did not attempt, and thus we excluded the affected source from the sample.

Fig. 2 Cutouts of sources 8 and 16 at radio frequencies, with S2CLS 850 μm. There are clear artefacts in the LOFAR image, which result from nearby bright sources. While low S/N detections look possible in several bands, these sources lie below the 5σ detection threshold in all bands other than 850 μm. The orange circle shows the 15 arcsec SCUBA-2 beam FWHM.

3.3 Non-detections

Two S2CLS detections – S2CLSJ104554+584754 and S2CLSJ104501+590403, sources 8 and 16, respectively, in our sample – are below the 5σ detection threshold in the LOFAR map. Inspecting these sources in all other available wavelengths, we also find them to be below the 5σ detection thresholds across the optical and IR as well as at all other radio frequencies (see Fig. 2).

While eight sources appear in the S2CLS catalogue that do not have counterparts in the LOFAR catalogue, these two sources are of particular interest due to their high signal-to-noise ratios (S/Ns). With S/Ns of 8.14 and 7.14, respectively, they are very unlikely to be false detections, the false detection rate in S2CLS being less than 0.001 per cent for a ≥ 7σ source. We focused on these two bright S2CLS sources as they are most likely to be real sources, and we did not attempt a further analysis of the fainter sources without significant multi-wavelength counterparts due to the probability of them being false detections.

A number of possible factors could lead to these bright submillimetre sources with very faint fluxes at other wavelengths: it is possible that these submillimetre detections are the result of blending of several faint sources with an angular separation smaller than the SCUBA-2 beam but larger than at other wavelengths. However, if we assume that these are indeed single sources, then these objects are either extremely reddened or they are at very high redshift (z > 4). Most likely, we are observing a combination of both redshift and reddening effects (Hughes et al. 1998; Walter et al. 2012), and thus these lone SCUBA-2 detections are candidate high redshift sources.

To estimate a limiting lower redshift at which these sources must be located to be detected at 850 μm and below the detection threshold at other wavelengths given the depths of the various surveys, we determined the detection thresholds in each of the images in which there is no detection. For the optical through to near infrared, we used 5σ detection thresholds from position-dependent depth maps (Shirley et al. 2019). For the Herschel Photodetector Array Camera and Spectrometer (PACS) and Spectral and Photometric Imaging Receive (SPIRE) images, we used the documented 5σ survey depths from Oliver et al. (2012). For Spitzer MIPS, we used the survey depths from the SWIRE second data release (DR2; Surace et al. 2005)².

As a simple illustration, we used the averaged SMG template SED from Michałowski et al. (2010) as an assumed, underlying SED for these two sources. Normalising at the observed 850 μm flux density, we redshifted the template until it lay below the detection thresholds of each wavelength (see Fig. 3). Thus we obtained a lower limit on the redshift of each source.

This results in limiting redshifts of z = 5.4 for source 8 and z = 4.7 for source 16, both of which are plausible for highly star-forming galaxies selected at 850 μm. Given the simplistic nature of our fitting procedure, and the assumptions about the underlying SED shape, we cannot give robust redshift limits on these sources; however, this is an indication of an approximate redshift that is consistent with our hypothesis that these are two submillimetre-bright, optically and radio-dim sources located at high redshift. Further investigation of these sources with submillimetre or millimetre instruments to obtain spectral line emission would be required to more accurately determine the redshifts of these objects.

Fig. 3 Calculating the lower limits on redshifts for sources 8 (left) and 16 (right), both of which are detected at 850 μm and at no other wavelength in this study. Using the 5σ detection thresholds (marked as downward pointing arrows) at optical through to FIR wavelengths, and the detection at 850 μm, we redshift the Michałowski et al. (2010) template SED until it lies below the limiting fluxes. In this way we determine limits of z = 5.4 and z = 4.7 for sources 8 and 16, respectively.

3.4 Optical to submillimetre spectral energy distributions

As described in Sect. 3.2, in the majority of cases there is a clear, single LOFAR counterpart to the submillimetre source. Therefore, to construct the SED we can simply use the LOFAR catalogue source and the corresponding flux densities and uncertainties at all optical and IR wavelengths provided in the LOFAR added-value catalogue (Kondapally et al. 2021) and the radio flux densities measured using the method described in Sect. 3.2.

Figure 4 shows an example SED for a source in our sample. We over-plot the Michałowski et al. (2010) template SED normalised at the observedframe 850 μm flux density and transformed to the photometric redshift of each source. This provides a visual comparison to a typical high redshift dusty star-forming galaxy SED, and we can see that the SEDs of sources in this study are broadly similar across the full wavelength range. As several sources in our sample are outside of Herschel survey area coverage, and therefore do not have data available with which to constrain the peak of dust spectrum, we did not conduct a full SED-fitting process but simply used the normalised star-forming galaxy templates. While this is not as robust as a fit to the observed SEDs, in the majority of cases the template traces the observed data adequately.

3.5 The far-infrared to radio correlation (FIRC)

As a quick check to confirm that none of our sources stand out as radio-loud AGN, we calculated estimates ofFIR luminosity for our sources and plotted them on the FIRC. We used a very simple method, integrating the Michałowski et al. (2010) average SMG template flux between 8 and 1000 μm in the rest frame. We then normalised this against the observed-frame 850 μm luminosity. While this is a naive approach and inevitably introduces some uncertainty, the Michałowski et al. (2010) template SED is qualitatively similar enough across the wavelength range in consideration to provide a reasonable estimate of the FIR luminosity of the sources.

The radio luminosity is calculated using the observed 1.4 GHz flux density, K-corrected according to the LOFAR catalogue redshift following . For each source, we used the value of the radio spectral index, α, resulting from a power law fit to the high-frequency radio spectrum (324 MHz–1.4 GHz) of the form S_ν ∝ ν^α (see Sect. 4). We plot our sample on the FIRC following Ivison et al. (2010), and over-plot their results in green for comparison (Fig. 5).

We notice immediately that, despite differences in the location of our sources on the FIRC compared to Ivison et al. (2010) (likely due to our FIR estimation method and selection effects, as discussed below), none of our sources appear significantly below the FIRC, and therefore we may assume that our sources do not host radio-loud AGN, or at least that their spectra are dominated by contributions from star formation. Selecting only the brightest submillimetre sources from a flux-limited sample introduces selection effects such that we bias our sample towards the FIR-bright, radio-faint population; therefore, we do not attempt to comment on the distribution of our sources with regards to the normalisation of the FIRC and its dependence on other galaxy properties, such as redshift, temperature, and stellar mass, as extensively studied by, for example, Yun et al. (2001), Ivison et al. (2010), Smith et al. (2014), and Read et al. (2018). The depth of the radio surveys used in our sample also impacts the FIRC here: many older studies are based on much shallower surveys, and the FIRC is plotted only for objects detected in both bands. With shallow survey flux limits, this biases the selection to only the brightest sources in both IR and radio. Here, selecting with S2CLS, we deliberately probed only the brightest, most star-forming IR galaxies, but with the depth of the VLA surveys we are able to detect the faint radio counterparts of these submillimetre-bright objects to lower flux limits than previous work. These faint radio sources detected in our sample havethe resulting effect that the distribution of our sources in Fig. 5 has a scatter that lies above the traditional FIRC. However, despite our sources filling a relatively small range of FIR luminosities, they span a wide range of radio luminosities, allowing us to study the variation in the radio properties of highly star-forminggalaxies. Calculating upper limits on the radio luminosities of the two sources with no radio detections, as described in Sect. 3.3, we find that these sources also lie in the same region on the FIRC as the rest of our sample, suggesting that these are likely to be normal star-forming galaxies at very high redshift, as hypothesised.

We calculated the FIRC parameter, q_IR, which parameterises the FIRC, as follows: (1)

We obtained a mean value of q_IR = 2.24 ± 0.29 for the objects in our sample. This is consistent with the comparable value of 2.41 ± 0.2 obtained by Ivison et al. (2010), despite our very simple method for estimating FIR luminosities and the selection effects discussed above.

To explore the range of the radio properties of these sources, we look in more detail at the radio spectra in the following section.

Fig. 4 Observed-frame SEDs of four example sources from our sample. Blue points show the flux densities across the spectrum, with radio flux densities measured as described in Sect. 3.2 and optical–FIR from the LOFAR catalogue. Uncertainties are plotted, though in many cases they are smaller than the size of the plotted points. The average SMG template from Michałowski et al. (2010) (orange) is over-plotted, normalised at the observed-frame 850 μm SCUBA-2 data point. We show several examples to demonstrate that, in most cases, the template traces the IR dust spectrum well (though in others, less well) but that due to the lack of Herschel data covering the peak of the dust spectrum in several sources, we use this normalised template to estimate FIR luminosities instead of template fitting.

Fig. 5 Far-infrared to radio correlation. Sources in our sample with extremely flat low-frequency radio spectral indices are plotted with light-coloured squares, and the positions of the ‘normal’ sources are plotted with dark circles. The FIRC relation calculated from our data following Eq. (1) is plotted as a dark dashed line, while the shaded region shows the 1σ and 3σ variances, respectively. Data from Ivison et al. (2010) are plotted in green crosses for comparison, with the FIRC from that work plotted as a dotted green line.

4 Radio spectra

The simplest approach to modelling the radio SED is to fit a power law, S_ν ∝ ν^α, and allow the spectral index, α, to vary. Fitting a power law across all available radio data for each source, we measured a median³ spectral index across the whole sample of − 0.86 ± 0.06, consistent with values of α in previous studies (Condon 1992; Sirothia et al. 2009; Mauch et al. 2013). If we divide the SED into two sections, low frequency (150–324 MHz) and high frequency (324 MHz–1.4 GHz), we find median spectral indices of − 0.56 ± 0.16 at low frequency and − 0.97 ± 0.15 at high frequency. This choice of dividing frequency (324 MHz) is still firmly in the traditional low-frequency regime; however, if we were to split the sample at a higher frequency (e.g. 610 MHz), we would risk averaging over the part of the spectrum where we might see the effects of spectral flattening most prominently. In addition, our radio data span the range 150 MHz–1.4 GHz, and we did not include any higher frequency radio data in our analysis, so our use of the terms ‘low frequency’ and ‘high frequency’ are purely relative in this case. The lower average spectral index at low frequency reveals an average spectral flattening, but across the sample we see a range of different spectral shapes in the radio regime (Fig. 6). Extremely star-forming galaxies have been found to have steeper spectral indices in our high-frequency range than in previous studies (Galvin et al. 2018), and models by Lacki & Thompson (2010) also predict steeper synchrotron spectra in high redshift ‘puffy’ SMGs than in compact starbursts.

This variety suggests that there is not one single set of physical conditions responsible for the shape of the radio spectrum, but that it is more likely that a variety of physical conditions contribute to the observed diversity of radio spectral shapes. With only four observed frequencies, we did not attempt to fit a more complex model but instead examined how the slope of the spectrum changes with frequency. We constructed a radio colour–colour plot using the observed-frame radio flux densities at 150 MHz, 324 MHz, and 1.4 GHz. Defining a colour, equivalent to the local spectral index, as the relationship between the two flux densities, (2)

we plot the low-frequency colour (α_low, ν₁ = 150 MHz, ν₂ = 324 MHz) against high-frequency colour (α_high, ν₁ = 324 MHz, ν₂ = 1.4 GHz) to obtain a diagnostic plot from which we may read a measure of spectral curvature (Fig. 7). A non-negligible number of sources have significantly flatter low-frequency spectral slopes than expected from typical estimates of power-law-type radio spectra with α ≃−0.8.

We selected the nine sources with α_low > −0.25 (marked in Fig. 7 with a dashed red line) to investigate further and performed a number of comparisons to see whether this extremely flat α_low sub-sample differs from the parent sample in any observational parameters. We note that conventionally a more conservative cut is made in spectral index to define flat-spectrum sources (typically α =−0.5); however, we divided our sample as such so as to investigate the most extreme spectral flattening. Taking a more traditional division between flat and steep spectrum sources at α = −0.5 does not, however, significantly affect our conclusions, and in fact a flat-spectrum sample with α >−0.5 is comparable to the extremely flat-spectrum α > −0.25 sample in terms of luminosity and redshift distribution (as in Fig. 8) as well as the other properties explored in this study.

As an initial check, we confirmed that the distribution of these sources is similar to that of the parent sample in both luminosity and redshift (Fig. 8). The extremely flat-spectrum sample is statistically indistinguishably distributed from the parent sample in both, with a Kolmogarov-Smirnov two-sample test (KS test) giving p-values of p = 0.98 and p = 0.61 for IR luminosity and redshift, respectively.

Next we investigated the MIR colours of these sources. If the difference in radio spectral shape is correlated with a difference in dust temperatures, we might expect to see this sub-sample in a distinct region of an MIR colour–colour diagram. Mid-infrared colour photometry can be used as a diagnostic to identify the presence of AGN in star-forming galaxies (Laurent et al. 2000). The MIR spectrum of a star-forming galaxy is typically driven by emission from warm dust (T ~ 25–50 K) heated by H II regions associated with recent star formation, as well as photodissociation regions giving rise to sharply peaked emission features from polycyclic aromatic hydrocarbons (PAHs). A luminous AGN will contribute significantly to the short-wavelength emission of the MIR spectrum via dust in the AGN torus being heated to much higher temperatures of up to T ~1500 K. Colour–colour plots constructed from Spitzer Infrared Array Camera (IRAC) photometry can be used to distinguish whether the dominant contribution is from AGN, in which case the MIR colours will be very red (Lacy et al. 2004). We plot the positions of our sources in MIR colour–colour plots using the 3.6, 4.5, 5.8, and 8 μm IRAC flux densities (Fig. 9).

We followed Donley et al. (2012) and over-plotted the region defined therein to identify galaxies hosting AGN using their IR colours. Ourmain interest in investigating the MIR colours is to determine whether the difference in radio spectral slope could result from AGN contributions. We have already established that there is no clear signature of radio-loud AGN in anysources; however, by inspecting the MIR colours we can probe their dust temperature distributions, which could reveal AGN activity with low radio luminosities that could affect the shape of the radio spectrum. A number of our sources lie within this region, suggesting that these sources may contain AGN. There are also a number of other mechanisms that could affect the MIR properties in addition to AGN activity – it should be noted that we cannot disentangle any differences in MIR emission due to the relative timescales of star-formation and AGN activity, nor that of the merger stage. The range of redshifts across the sample will also result in a spread of MIR colours due to the location of PAH features in the observed frame. Nevertheless, there is no clear distributional difference between the extremely flat-spectrum sample and the rest of the sample – a KS test giving p = 0.4 in the S_5.8 ∕S_3.6 colour and p = 0.95 in the S_8.0∕S_4.5 colour – with neither sample clearly occupying a distinct region in the colour–colour plot, suggesting that the cause of this radio spectral flattening must be something other than simply AGN contributions to the spectrum. As the sample is of only a small number of sources, it is difficult to make a statistically robust statement on their distribution; however, there is nothing to suggest that the extremely flat-spectrum sources are special in any way as regards their IR properties.

As there are no discernible statistically significant differences between the extremely flat α_low sources and the rest of the sample in our observable properties, we propose that this difference in radio spectral shape may be related to properties of the galaxies that we are unable to detect with our unresolved, galaxy-averaged observations.

Fig. 6 Examples of radio SEDs from our sample, with data at 150 MHz from LOFAR, 324 MHz and 1.4 GHz from the JVLA, and 610 MHz from the GMRT (see Sect. 2.3). There is a range of luminosities and spectral shapes. Orange lines show a simple power law with spectral index α = −0.7 normalised to the 610 MHz point as a reference and to demonstrate the spectral curvature present in many of these radio spectra. Error bars are shown but are largely within the size of the plotted points.

Fig. 7 Radio colour–colour plot showing the low-frequency (150 MHz –324 MHz) spectral slope on the x-axis and the high-frequency (324 MHz–1.4 GHz) spectral slope on the y-axis. The dashed line indicates αlow =−0.25, the cutoff above which we define sources to have extremely flat low-frequency spectra. Sources above this threshold are marked with squares in the plots throughout this paper. Sources are numbered with the IDs assigned in this study, as detailed in Table 1.

Fig. 8 Normalised histograms showing the distribution of sources in IR luminosity and redshift, with the full sample shown in dark and the extremely flat-spectrum sources shown in a light outline. The median of each sample is also shown, with a dashed line showing the median of the full sample and a dotted line showing the median of the extremely flat-spectrum sample.

Fig. 9 IRAC colour–colour plot, with the AGN ‘wedge’ from Donley et al. (2012) shown with a dashed line. Light squares show the extremely flat-spectrum sources.

Low-frequency flattening of radio spectra

There are several possible scenarios in which the low-frequency radio spectral slope is flattened relative to an α =−0.7 power law spectrum resulting from the conditions of the ISM, as observed in the nine sources described in Sect. 4. Synchrotron self-absorption can play a role in flattening the low-frequency spectrum; however, as the galaxies in our sample do not have AGN-dominated spectra, this is unlikely to make a significant contribution. In the case of low-luminosity AGN contribution, the level of flattening observed in our extremely flat-spectrum sources could only be obtained with a combination of the maximum possible AGN fractional contribution and self-absorption tuned to occur exactly at our breaking frequency of 324 MHz. This would be an unlikely coincidence, and as such we conclude that it is unlikely for AGN contributions alone to cause the observed flattening. A more likely possibility is that, if the source of synchrotron emission is embedded in an ISM that is sufficiently dense and clumpy, free–free absorption begins to have a significant effect as we observe the lower-frequency spectrum. Nearby star-forming galaxies have been observed to have significant free–free absorbed spectra and clumpy star-forming regions (e.g. Lenc & Tingay 2009; Rampadarath et al. 2014) but it is only now with the deep low-frequency radio data that have become available with LOFAR that we are able to investigate this at high redshift.

To model free–free absorbed spectra, we assumed typical properties of the ISM of SMGs and computed the effect this has on an α = −0.7 power law spectrum at the median (estimated) redshift z = 2.6 of the sample. We calculated the optical depth τ = ∫ κdl at the rest-frame frequencies that correspond to our observed frequencies for a galaxy at the median redshift of the sample (z = 2.6), where κ is the free–free absorption coefficient as defined by Condon (1992): (3)

where n_e is electron density and T_e electron temperature. Estimates of electron densities in star-forming galaxies from emission line diagnostics are poorly constrained and range from 10 to 400cm⁻³ (Parkin et al. 2013; Farrah et al. 2013; Masters et al. 2014; Kaasinen et al. 2017). We assumed an electron density of n_e = 50 cm⁻³, which is on the low end of star-forming electron density estimates – however, as the length scale of relevance is inversely proportional to for an equivalent level of absorption, if we were to assume a higher value of n_e we would find length scales that are smaller by a factor of . This choice then gives us an approximate upper limit for the length scales required for free–free absorption to flatten the spectrum at low frequency, as observed. We assumed an electron temperature of T_e = 10⁴K (Downes et al. 1980). Of course, the properties of the ISM in an individual galaxy will vary in density and temperature across the galaxy’s extent as well as across our sample, which covers a range of redshifts; we expect to see variation in the nature of the ISM. This simple model of the effects of free–free absorption assumes a single slab of absorbing material; due to the nature of our observations, it would be unrealistic to attempt to model any more realistic distributions of dense gas in the ISM. However, more complex distributions of gas mixed with the radio-emitting plasma in simulations of free–free absorption in radio spectra (e.g. Bicknell et al. 2018; Varenius et al. 2014) produce very similar results overall.

Figure 10 shows the effect of free–free absorption on such a spectrum, over absorbing columns of between zero and 300 pc. In this model, the level of spectral flattening that we would be able to observe occurs at rest-frame frequencies of < 1 GHz, which corresponds to observed-frame frequencies of ≲280 MHz for the median redshift, z = 2.6, of our sample. This in part motivated our choice of cutoff frequency to measure the spectral slope as measuring at 610 MHz would average out the flattening we see most strongly at these lower frequencies. We return to the radio colour–colour plot in Fig. 7 and plot the colours of this model spectrum over these increasing absorbing columns. This simple model demonstrates that, with our assumed electron densities and temperature, an absorbing column of several hundred parsecs in size would be sufficient to flatten the low-frequency radio spectral index to the degree that we observe in several sources. As noted above, if we increase the electron density in our model, this would be consistent with smaller-scale absorbing columns (approximately tens of parsecs).

This spatial scale of dense, star-forming clumps is consistent with the highest-resolution observations of SMGs (Swinbank et al. 2010; Ikarashi et al. 2015; Dessauges-Zavadsky et al. 2019). A number of studies imaging SMGs at the highest resolutions have found evidence of clumpy substructure at spatial scales of several hundred parsecs (e.g. Iono et al. 2016; Tadaki et al. 2018); however, for the most part it is only possible to resolve down to scales of approximately tens of kiloparsecs in the handful of submillimetre sources that are strongly gravitationally lensed and observed with interferometric instruments. One such strongly lensed submillimetre source that has been well studied, SDP.81, displays a spatially non-uniform dust continuum, with several bright 100−500 pc scale clumps (Rybak et al. 2015; Hatsukade et al. 2015; Swinbank et al. 2015) and ‘knots’ (Tamura et al. 2015) on < 100 pc scales. The ‘Cosmic Snake’ is another – with magnification affording resolved scales as small as 30 pc; Dessauges-Zavadsky et al. (2019) find several molecular clouds on scales between 30–210 pc.

These scales correspond approximately to the Jeans length for the gas densities typical of star-forming galaxies (Shimakawa et al. 2015; Hatsukade et al. 2015), with giant molecular gas clouds collapsing as a result of Toomre instabilities (Toomre 1964). However, we caution the reader that there has been debate as to whether observed ‘clumpy’ substructure is biased by effects relating to interferometric observations. The case for ~100 pc scale substructure in the ‘Cosmic Eyelash’ (Swinbank et al. 2010) has recently been challenged, with work showing that the inferred structure may be due to filtering and resolution effects amplifying spurious features in low S/N interferometric images (Cava et al. 2018; Gullberg et al. 2018; Ivison et al. 2020).

Our extreme low-frequency spectral flattened sample with α > −0.25 consists of only nine of the 42 sources in the whole sample, implying that if our assumption that this spectral flattening is caused by free–free absorption is correct, the majority of the population of submillimetre-bright sources might not be expected to show strong evidence of dense, clumpy star-forming regions on these scales. The effect of absorption on our unresolved galaxy-averaged observations must depend on the fraction of low-frequency radio emission that is embedded within and absorbed by high density gas; if a sufficient fraction is able to escape, we will see this flattening to a lesser degree. Thus there are several explanations for the lack of free–free absorption signatures in the majority of our galaxy sample. It could be due to an intrinsically smoother, more diffuse distribution of star-forming material (as suggested in e.g. Hodge et al. 2016; Ivison et al. 2020). There is likely also a dependence on age – in younger star-forming regions, the radio emission is more likely to still be contained within a dense gaseous structure – and geometry may also play a role. Effects of environment and the merger stage may affect the distribution of gas in our galaxies, and we spanned a large range in redshift over which we would expect the nature of the ISM to evolve. Additionally, our estimated ~100 pc length scale is based on a conservatively low estimate of electron density, as previously discussed; higher assumed electron densities would lead to clumps on scales below the resolution limits of observations of even the most highly magnified lensed sources.

Due to the serendipitous nature of locating gravitationally lensed galaxies, there are very few sources in which observations of sufficiently high resolution can be made to detect substructure on the ~100 pc scales implied by our calculation. Since only ~20 per cent of sources in our sample display this low-frequency radio spectral flattening, a much larger sample of hundreds of galaxies observed at sub-kiloparsec scales would be required to detect many sources with this substructure if our assumptions are correct that this observed radio spectral flattening is due to free–free absorption, and assuming that structure does exist on scales large enough to be detected with current instrumental capabilities.

Fig. 10 Our simple model of free–free absorption and its effect on spectral indices. Left: effect of free–free absorption on a simple power law spectrum with radio spectral index α = −0.7, with increasing lengths of absorbing column, for our observing frequencies in the observed frame. We assume an electron density ne = 50cm−3 and temperature T = 104K. Right: as Fig. 7, but with the effect of free–free absorption on a simple power law spectrum plotted in orange, with increasing lengths of absorbing column from zero (left) to 300 pc (right).

5 Conclusions

Taking advantage of new, deep LOFAR images, we investigated the low-frequency radio spectra of a sample of highly star-forming galaxies selected at 850 μm from the S2CLS. Our conclusions are as follows:

• 1.

  We find that this sample of SMGs displays a range of radio luminosities and spectral shapes despite being selected at a very narrow range of submillimetre fluxes, implying that the radio properties of this sample do not follow a tight correlation with the star-formation properties we might infer from submillimetre observations alone.

• 2.

  We find evidence of radio spectral flattening at low frequencies (α_low = −0.47 ± 0.16 on average, and nine sources with α_low < −0.25). These flat-spectrum sources are indistinguishable from the full sample in their distributions of redshift and IR luminosity, as well as in their IR colours. In the absence of any clear observational differences between sources with flat low-frequency spectra and the rest of the sample, we infer that this must be due to underlying properties of the galaxies that we cannot observe in our unresolved imaging. We suggest that this radio spectral flattening may be due to free–free absorption arising from non-uniform, clumpy distributions of ionised gas in which star formation is embedded. Taking typical values of electron density and temperature, we estimate that clumps must be of the order of a few hundred parsecs in size (comparable to substructure that has been observed in strongly lensed submillimetre sources) to account for the observed radio spectral curvature. This presents an additional piece of evidence in favour of clumpy star formation in high redshift galaxies that does not depend on the angular resolution of morphological imaging but can be detected in the galaxy-averaged properties in the radio SED. Due to the serendipitous nature of imaging lensed sources, there are few observations reaching high enough resolution to detect structure on this scale, and so larger samples of high-resolution images (e.g. using submillimetre interferometric instruments such as ALMA) would be required to test whether the proportion of galaxies in our sample that exhibit this radio spectral flattening (~20 per cent) is typical of the submillimetre population in general.

• 3.

  In addition to the range of radio spectra observed in this sample, we also find two bright submillimetre sources (>7σ detections) that are undetected at all other wavelengths, from optical through to radio. We speculate that, due to samplingthe peak of the thermal dust emission spectrum at 850 μm, these objects are located at high redshift (z > 4) and are too faint due to cosmological dimming at all other observed wavelengths. We propose them as interesting candidates for (sub)millimetre interferometric follow-up to determine their redshifts with certainty.

Finding this variety of spectral shapes and evidence of free–free absorption in spectral flattening at low frequencies is consistent with high-resolution observations of clumpy star-forming regions in submillimetre galaxies. A larger sample of submillimetre galaxies with low-frequency radio observations and high-resolution interferometric follow-up would be beneficial to further investigating the nature of this spectral flattening.

Acknowledgements

J.R., M.J.H., P.N.B., I.M., K.C. and R.K. acknowledge support from the UK Science and Technology Facilities Council [STFC: ST/N504105/1, ST/R000905/1, ST/R000972/1, ST/R505146/1 and ST/R504737/1]. J.E.G. is supported by the Royal Society through a University Research Fellowship. K.J.D. and H.R. acknowledge support from the ERC Advanced Investigator programme NewClusters 321271. M.J.J. acknowledges support from the UK Science and Technology Facilities Council [ST/N000919/1] and the Oxford Hintze Centre for Astrophysical Surveys which is funded through generous support from the Hintze Family Charitable Foundation. M.B. acknowledges support from INAF under PRIN SKA/CTA FORECaST and from the Ministero degli Affari Esteri della Cooperazione Internazionale - Direzione Generale per la Promozione del Sistema Paese Progetto di Grande Rilevanza ZA18GR02. IP acknowledges support from INAF under the SKA/CTA PRIN “FORECaST” and the PRIN MAIN STREAM “SAuROS” projects. K.C. is supported by a Royal Society Leverhulme Senior Research Fellowship (SRF/R1/191013). LOFAR, the Low Frequency Array designed and constructed by ASTRON, has facilities in several countries, which are owned by various parties (each with their own funding sources), and are collectively operated by the International LOFAR Telescope (ILT) foundation under a joint scientific policy. The ILT resources have benefited from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; the Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland. Part of this work was carried out on the Dutch national e-infrastructure with the support of the SURF Cooperative through grant e-infra 160022 & 160152. The LOFAR software and dedicated reduction packages on https://github.com/apmechev/GRID_LRT were deployed on the e-infrastructure by the LOFAR e-infragroup, consisting of J. B. R. Oonk (ASTRON and Leiden Observatory), A. P. Mechev (Leiden Observatory) and T. Shimwell (ASTRON) with support from N. Danezi (SURFsara) and C. Schrijvers (SURFsara). This research has made use of the University of Hertfordshire high-performance computing facility (https://uhhpc.herts.ac.uk/) and the LOFAR-UK compute facility, located at the University of Hertfordshire and supported by STFC [ST/P000096/1]. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC). The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; the Operation, Maintenance and Upgrading Fund for Astronomical Telescopes and Facility Instruments, budgeted from the Ministry of Finance (MOF) of China and administrated by the Chinese Academy of Sciences (CAS), as well as the National Key R&D Program of China (No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. The authors wish to recognize the very significant role that the summit of Maunakea has within the Indigenous Hawaiian (Kānaka Maoli) community, and acknowledge the continued use of indigenous land for telescope facilities. We made use of SciPy (Jones et al. 2001), NumPy (Van Der Walt et al. 2011) and matplotlib, a Python library for publication-quality graphics (Hunter 2007), as well as Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013).

References

All Tables

All Figures

	Fig. 1 Cutouts of an example S2CLS source (ID 1 in our numbering system, an 11.91 mJy 850 μm source) in each radio frequency used in this study. Each square is 50 arcsec across, with the approximate S2CLS beam size (~ 15 arcsec FWHM) marked with an orange circle. Image contrast is scaled arbitrarily for clarity.
In the text

	Fig. 2 Cutouts of sources 8 and 16 at radio frequencies, with S2CLS 850 μm. There are clear artefacts in the LOFAR image, which result from nearby bright sources. While low S/N detections look possible in several bands, these sources lie below the 5σ detection threshold in all bands other than 850 μm. The orange circle shows the 15 arcsec SCUBA-2 beam FWHM.
In the text

Fig. 3 Calculating the lower limits on redshifts for sources 8 (left) and 16 (right), both of which are detected at 850 μm and at no other wavelength in this study. Using the 5σ detection thresholds (marked as downward pointing arrows) at optical through to FIR wavelengths, and the detection at 850 μm, we redshift the Michałowski et al. (2010) template SED until it lies below the limiting fluxes. In this way we determine limits of z = 5.4 and z = 4.7 for sources 8 and 16, respectively.

In the text

Fig. 4 Observed-frame SEDs of four example sources from our sample. Blue points show the flux densities across the spectrum, with radio flux densities measured as described in Sect. 3.2 and optical–FIR from the LOFAR catalogue. Uncertainties are plotted, though in many cases they are smaller than the size of the plotted points. The average SMG template from Michałowski et al. (2010) (orange) is over-plotted, normalised at the observed-frame 850 μm SCUBA-2 data point. We show several examples to demonstrate that, in most cases, the template traces the IR dust spectrum well (though in others, less well) but that due to the lack of Herschel data covering the peak of the dust spectrum in several sources, we use this normalised template to estimate FIR luminosities instead of template fitting.

In the text

Fig. 5 Far-infrared to radio correlation. Sources in our sample with extremely flat low-frequency radio spectral indices are plotted with light-coloured squares, and the positions of the ‘normal’ sources are plotted with dark circles. The FIRC relation calculated from our data following Eq. (1) is plotted as a dark dashed line, while the shaded region shows the 1σ and 3σ variances, respectively. Data from Ivison et al. (2010) are plotted in green crosses for comparison, with the FIRC from that work plotted as a dotted green line.

In the text

Fig. 6 Examples of radio SEDs from our sample, with data at 150 MHz from LOFAR, 324 MHz and 1.4 GHz from the JVLA, and 610 MHz from the GMRT (see Sect. 2.3). There is a range of luminosities and spectral shapes. Orange lines show a simple power law with spectral index α = −0.7 normalised to the 610 MHz point as a reference and to demonstrate the spectral curvature present in many of these radio spectra. Error bars are shown but are largely within the size of the plotted points.

In the text

Fig. 7 Radio colour–colour plot showing the low-frequency (150 MHz –324 MHz) spectral slope on the x-axis and the high-frequency (324 MHz–1.4 GHz) spectral slope on the y-axis. The dashed line indicates αlow =−0.25, the cutoff above which we define sources to have extremely flat low-frequency spectra. Sources above this threshold are marked with squares in the plots throughout this paper. Sources are numbered with the IDs assigned in this study, as detailed in Table 1.

In the text

	Fig. 8 Normalised histograms showing the distribution of sources in IR luminosity and redshift, with the full sample shown in dark and the extremely flat-spectrum sources shown in a light outline. The median of each sample is also shown, with a dashed line showing the median of the full sample and a dotted line showing the median of the extremely flat-spectrum sample.
In the text

	Fig. 9 IRAC colour–colour plot, with the AGN ‘wedge’ from Donley et al. (2012) shown with a dashed line. Light squares show the extremely flat-spectrum sources.
In the text

Fig. 10 Our simple model of free–free absorption and its effect on spectral indices. Left: effect of free–free absorption on a simple power law spectrum with radio spectral index α = −0.7, with increasing lengths of absorbing column, for our observing frequencies in the observed frame. We assume an electron density ne = 50cm−3 and temperature T = 104K. Right: as Fig. 7, but with the effect of free–free absorption on a simple power law spectrum plotted in orange, with increasing lengths of absorbing column from zero (left) to 300 pc (right).

In the text