© The Authors 2023

1. Introduction

One of the main challenges in the field of extra-galactic astronomy is understanding how the star formation rate density (SFRD – M_⊙ yr⁻¹ Mpc⁻³) evolves with redshift (see Madau & Fragos 2017; Oesch et al. 2018; Tacconi et al. 2020, for a detailed review). While we have accurate estimates of the SFRD at intermediate redshifts (1 < z < 3) using various selection methods from deep multi-wavelength data (e.g., Gruppioni et al. 2020; Novak et al. 2017; Ishigaki et al. 2018; Dunlop et al. 2017, for infrared (IR), radio, ultraviolet-UV and millimetre (mm) criteria, respectively) that agree with results from hydrodynamical simulations (e.g., Pillepich et al. 2018; Lagache 2018), uncertainties regarding higher redshifts (z > 3) still remain (e.g., Magnelli et al. 2019; Gruppioni et al. 2020). For instance, studies of UV-selected galaxy samples (e.g., McLeod et al. 2016; Ishigaki et al. 2018) estimate a steep decline in the SFRD at z > 3, while other works employing far-IR radio selection criteria (e.g., Viero et al. 2013; Rowan-Robinson et al. 2016; Novak et al. 2017; Malefahlo et al. 2022) or gamma-ray-burst-selected samples (e.g., Yüksel et al. 2008; Kistler et al. 2009) suggest a flatter behaviour.

These uncertainties may originate from our current instrumentational limitations, including low sensitivity, coarse spatial resolution, and coverage of small parts of the sky, which hinder our ability to properly characterise distant objects and obtain a statistically unbiased sample at higher redshifts. For instance, observations based on telescopes such as the Herschel Space Observatory (Pilbratt et al. 2010) or the Submillimetre Common-User Bolometer Array (SCUBA; Holland et al. 2013) preferentially select sources at intermediate redshifts (i.e. z ∼ 2 – 3; e.g. Pearson et al. 2013), and uncover only the very bright end of the far-IR luminosity function (e.g., Gruppioni et al. 2013; Wang et al. 2021). On the other hand, more sensitive telescopes, such as the Atacama Large Millimeter/submillimeter Array (ALMA; Wilson et al. 2019) or the Northern Extended Millimeter Array (NOEMA; Chenu et al. 2016) are limited to smaller areas. Thereby, these telescopes provide fewer objects (e.g., Dunlop et al. 2017; Faisst et al. 2020; Casey et al. 2021) and reveal a fainter star-forming galaxy (SFG) population (e.g., Pantoni et al. 2021). Nevertheless, recent studies using the capabilities of both single telescopes and interferometers have pushed the limits of our current understanding regarding the star formation (SF) history, exposing a mm dusty galaxy population at z > 3 (see e.g., Donevski et al. 2018; Yamaguchi et al. 2019; Gruppioni et al. 2020).

Another aspect of our current limitations regarding the accurate estimation of the SFRD at high redshifts might be hidden in the selection criteria that are currently used to detect SFGs. Various methods have been employed over the years in order to detect this galaxy population (e.g., Gruppioni et al. 2013; Bouwens et al. 2015), however, the bulk of the research is focused on the detection of Lyman-break galaxies in the rest frame UV and the application of dust-correction low-redshift empirical relations for the derivation of their star formation rates (e.g., Bouwens et al. 2015; Oesch et al. 2018). These relations might not be necessarily valid at z > 3, where the SFG population remains uncertain. Alternatively, recent studies using sensitive far-IR or mm telescopes (e.g., the ALPINE-ALMA survey; Gruppioni et al. 2020) or radio surveys (e.g., Talia et al. 2021; Enia et al. 2022) aim to reveal a an SFG population that has so far remained undetected in the optical/UV bands, which is estimated to significantly contribute to the SFRD at z > 3 (e.g., Gruppioni et al. 2020). In order to obtain accurate SFR estimates, these studies separate the SFG from the AGN population (e.g., Lanzuisi et al. 2017; Gruppioni et al. 2020) by removing the AGN contribution from the IR luminosity. However, this approach can introduce a bias in the sample. Interestingly, results from the recent surveys of the James Webb Space Telescope (JWST; e.g., Harikane et al. 2023; Nakajima et al. 2023; Pérez-González et al. 2023; Rodighiero et al. 2023) reveal a dusty-galaxy population like this at even higher redshifts (i.e., z > 6) that have not been explored before.

Previous studies, such as the work by Duncan et al. (2019), indicated that the galaxy pair and merging ratio increases with redshift (from ∼5% at z < 1 to ∼40% at z > 3 and for M_⋆ > 10^10.3 M_⊙). Therefore, in this work, we question the strategy of removing AGN-contaminated radio sources from selected samples as this approach may overlook merging systems in which the AGN emission conceals an additional SFG population that could significantly contribute to the SFRD at high redshifts. Such galaxies could be detected from the identification of far-IR bright Herschel galaxies along with a radio emission. Even though these selection criteria might favour sources with cold gas and negative k-correction, their application to fields with deep coverage (e.g., the COSMOS field; Scoville et al. 2007) along with photometric data and high spatial resolution information from telescopes such as the Hubble Space Telescope (HST) and ALMA could potentially reveal a population of SFGs that would be omitted by the current selection criteria because of the AGN contribution of a neighbour galaxy.

In an attempt to explore new techniques to assist radio surveys reaching z > 3, we present a complementary study combining the radio selection and far-IR peak spectrum criteria (e.g., Yan et al. 2022). Furthermore, we suggest that at z ∼ 3 − 4, a substantial amount of star formation occurs that is hardly revealed from optical or near-IR surveys, and which is consequently excluded by the current selection methods due to the AGN contribution in the system. In more detail, we explore a sample of z > 3 far-IR radio-emitting galaxies, detected in the COSMOS and ECDF-S fields, that might uncover interacting systems and merging galaxies with a substantial contribution of an AGN component, along with the presence of high star formation activity. If radio emission can help to pinpoint this important ‘missing’ SF component, upcoming radio surveys, such as the Evolutionary Map of the Universe (EMU) survey (e.g., Norris 2009; Norris et al. 2021) as well as novel selection criteria in the radio (e.g., Talia et al. 2021; Enia et al. 2022) could become even more important to our understanding of the SFRD at z > 3.

The paper is structured as follows: in Sect. 2 we present the method and criteria that were employed in the selection of the original sample and the subsequent reduction of the sample to the high-redshift far-IR radio candidates. Section 3 briefly highlights the photometric and spectroscopic properties of the sample, along with the results obtained from various properties of the entire sample. In Sect. 4 we compare our work with previous studies and provide a brief discussion of the results, emphasising the importance of the merging population for SFRD estimations. Finally, we provide the conclusions of this study in Sect. 5 and discuss future prospects regarding our method. Throughout this paper, we assume a Λ-CDM cosmology using H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7.

2. Sample selection

For our selection, we examined two of the most deeply explored extra-galactic fields that were observed under the projects Extended Chandra Deep Field-South (ECDF-S; Lehmer et al. 2005) and Cosmic Evolution Survey (COSMOS; Scoville et al. 2007). In past years, these two fields have been extensively observed in a wide range of electromagnetic frequencies, ranging from X-rays to radio waves. Radio data were obtained from both fields at 1.4 GHz with the VLA (VLA-COSMOS; the flux density limit is 7 μJy beam⁻¹, and VLA-ECDF-S; the flux density limit is 7.4 μJy beam⁻¹; Schinnerer et al. 2010; Smolčić et al. 2017; Miller et al. 2013) at an angular resolution of ∼2.5″.

The far-IR wavelength regime, which is necessary for our selection, is available through the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012). This survey covered a total of 380 deg² using the two main instruments on board Herschel, namely the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010; Nguyen et al. 2010; at 250, 350 and 500 μm; the confusion noise level is ∼5.8, ∼6.3, and ∼6.8 mJy, respectively) and the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010, ; at 100 and 160 μm; the detection limits are ∼7.7 and ∼14.5 mJy, respectively). The ECDFS and COSMOS fields are also included in the covered area of HerMES, comprising a total of ∼2 deg². Previous studies using HerMES data have demonstrated the efficient use of Herschel’s far-IR fluxes for detecting high-redshift (z > 3) dusty galaxies. For instance, Dowell et al. (2014) employed simple criteria for the increasing fluxes for the 250, 350, and 500 μm bands, which is indicative of a far-IR peak at longer (redshifted) > 350 μm wavelengths, revealing a dusty SF galaxy population with a mean redshift distribution of ⟨z⟩∼4.7.

To identify potential high-redshift hybrid AGN-SFG systems, we conducted a cross-match analysis between the VLA COSMOS and ECDF-S 1.4 GHz catalogues with those of HerMES, applying a maximum separation of 18 arcsec (motivated by the beam size of Herschel at 250 μm of 18.1 arcsec). The resulting cross-matched sample contained 2292 sources for COSMOS, and for the ECDF-S field is provided 769 sources. This sample was further reduced by considering robust radio sources (signal-to-noise ratio; S/N > 5) without an indication for extended or multicomponent radio emission (according to each catalogue; for more details, see Lehmer et al. 2005; Scoville et al. 2007). To ensure the correspondence between radio and far-IR emission, we excluded from the sample sources that present multiple radio detections within one Herschel beam. Furthermore, only sources with 250, 350, and 500 μm flux errors smaller than 50% were accepted, which ensured that the shape of the far-IR emission was properly determined. This selection method resulted in 219 and 93 sources for COSMOS and ECDF-S, respectively.

Furthermore, motivated by the efficacy of studies such as Dowell et al. (2014) and Donevski et al. (2018) in selecting high-z dusty SF galaxies, we adopted a similar method to identify sources with f_350 μm/f_250 μm > 1, which would translate into a potential z > 3 redshift. This selection method resulted in a sample of 53 and 30 galaxies (see Fig. 1) for the COSMOS and ECDF-S fields, respectively. The radio 1.4 GHz and far-IR (250 μm) images of each of these sources were visually inspected in order to secure the robustness of the matching process and identify potentially uncertain or confused radio/far-IR detections. We confirmed that all 83 selected sources display compact morphologies both in the radio and in the FIR, and that the detected counterparts are robust, with average separations between the radio and far-IR detections below ∼2 arcsec.

Fig. 1. Flux density ratios of 350 and 250 μm vs. the ratios of 500 and 350 μm of HerMES in COSMOS (panel a) and ECDF-S (panel b). The sample was selected by cross-matching the HerMES catalogue with the VLA-COSMOS and VLA-ECDFS data. The dashed black lines indicate f350/f250 = 1 and f500/f350 = 1, while the points surpassing the limit of f350/f250 > 1 are presented with larger markers. The colour bar corresponds to the 350 μm flux of each source, and the legend shows the number of sources for each panel and the selected high-redshift candidates (crosses). The black lines in the subregion of each panel provide information regarding the shape of the far-IR (FIR) emission. For instance, sources with f500/f350 < 1 and f350/f250 > 1 have a peak in flux between 250 and 500 μm. The photometric properties of the sources presented in the legend are provided in Table 1.

2.1. High-redshift candidates

Application of these far-IR criteria resulted in a sample of 53 + 30 = 83 sources from both fields for which we retrieved possible photometric or spectroscopic identifications from the literature. For this exploration, we consulted Simbad (Wenger et al. 2000), the NED¹ online search tools, and the Cosmos2015 (Laigle et al. 2016) online catalogue², retrieving redshift values from sources that are within 2 arcsec of our candidates.

From this analysis, we found 18 sources in the ECDF-S and 19 sources in the COSMOS field with spectroscopic redshift determination. Five of these radio sources present a spectroscopic redshift of z_spec > 3, namely J033140.1−275631 (ECDF-S field; z_spec = 5.14; Suh et al. 2015), J100028.71+023203.7 (COSMOS field; z_spec = 4.76; Hasinger et al. 2018), [LBX2017] 214 (ECDF-S field; z_spec = 3.74; Iwasawa et al. 2012), COSMOSVLA J100256.53+021158.4 (COSMOS field; z_spec = 3.503; Pâris et al. 2014), and [DSS2017] 554 (ECDF-S field; z_spec = 3.1988; Danielson et al. 2017). The reported spectroscopic redshifts of the two highest-z sources, however, are not reliable. For instance, the candidate with the highest redshift, found in the ECDF-S, reveals contradicting spectroscopic measurements. More specifically, the z_spec = 5.137 value (derived from near-IR spectroscopy; Suh et al. 2015) contradicts Coppin et al. (2012) and Danielson et al. (2017), who estimated a lower spectroscopic value of z_spec = 1.617 (with spectra from optical and far-IR wavelengths) for this candidate. Because a high-redshift nature of this source cannot be excluded, we further explored its available spectra. This investigation indicated that the emission line responsible for the z_spec = 5.137 solution can be also interpreted as an emission line at z_spec = 1.617, which further supports the notion for the intermediate-redshift solution derived in the literature. Therefore, we considered this source as an intermediate-redshift galaxy and did not explore it further.

Likewise, for J100028.71, inspection of its DEIMOS 10k survey spectrum (Hasinger et al. 2018) revealed a very low S/N spectrum without an obvious line detection that could explain the z = 4.76 indication. We assumed an erroneous redshift estimate caused by a possible false-positive of the automated feature recognition from the DEIMOS survey. This source is however kept in this work as a high-redshift source, prompted by a photometric redshift determination of z_phot ∼ 3.2 by Miettinen et al. (2015). Inspection of the remaining three sources with a spectroscopic redshift higher than 3 revealed that the emission lines used to determine the redshifts of these sources have either a high S/N (i.e., S/N > 10) or are denoted with a secure spectroscopic quality flag. This means that these measurements are accurate.

Our search identifies 12 sources in the ECDF-S field and 34 sources in the COSMOS field based on photometric redshifts from the literature. These results are depicted in Fig. 2, where the histograms for the redshifts of each survey and the total values are plotted with different colours. This figure reveals that our selection favours radio sources around z  ∼  2 with a distribution that extends to higher redshifts. Four sources of this sample in the COSMOS field have photometric redshifts (identified in the Cosmos2015 catalogue) of z > 3 (i.e., COSMOSVLA J100233.16+020626.5 with z_phot = 3.016, COSMOSVLA J095845.94+024329.2 with z_phot = 3.172, COSMOSVLA J100028.75+023203.5 with z_phot = 3.175, and COSMOSVLA J100004.81+023045.2 with z_phot = 3.76). Similarly, in the ECDF-S sample, we identified one source with z_phot > 3, namely ALESS J033152.49−280319.1 (z_phot = 3.56; Wardlow et al. 2011). These five high-redshift candidates, along with the three sources with a spectroscopic redshift above 3, are presented along with their 1.4 GHz and far-IR HerMES flux densities in Table 1.

Fig. 2. Histograms presenting the number of sources for each redshift bin selected by the criteria for the COSMOS (empty histograms) and ECDF-S (filled histograms). The blue histograms correspond to spectroscopic values derived from the literature, and the orange histograms show the photometric values. The combined sample is depicted with an empty black histogram. The number of sources for each category is presented in the legend.

We are particularly interested in sources with the highest redshifts. We therefore decided to further investigate the five sources for which the photo-z indicates a z > 3 nature. Follow-up observations of these sources aimed to improve the redshift determination (for the photometric redshift estimates), and a more detailed analysis of their properties was initiated, using mm observatories from the northern and southern hemispheres (for the COSMOS and ECDF-S fields, respectively), the first results of which are described in the following section.

The comparison of our acquired spectroscopic/photometric redshift values for our sample with the corresponding values from alternative photometric catalogues in the literature is a valuable exercise. This comparison is presented in more detail in Appendix A and indicates that the COSMOS and the ECDF-S fields do not contain additional candidates at z > 3 that our analysis could have miss-identified at lower z, and the z > 3 nature of the five sources with photo-z considered here is confirmed. Therefore, we are more confident that the sources that are presented in the following sections are the only ones with z_phot, spec > 3 that could be selected with our selection criteria from the two fields of study.

3. Results

To characterise the morphological and physical properties of the eight high-redshift galaxies under investigation, a detailed examination of the available photometric and spectroscopic data from the literature was conducted. Simultaneously, follow-up observations were initiated, in particular, using millimetre facilities. The results are described below.

During this analysis, a spectral energy distribution (SED) fitting analysis was performed for each galaxy in the sample using CIGALE (Noll et al. 2009; Serra et al. 2011). The resulting fits yield two values for each modelled physical property of a galaxy (e.g., stellar mass and redshift), one corresponding to the best-fit solution and another that is derived through a Bayesian statistical analysis. In the following text, we refer to the best-fit values when presenting the sample, while the Bayesian-derived estimates along with their corresponding uncertainties are introduced in Table 3. The images of the eight sources for four different photometric bands (optical, near-IR, radio, and mm wavelengths) are presented in Fig. 4. For each wavelength regime, the observations with the highest available spatial resolution were selected.

3.1. SED fitting results

The SED fitting analysis was conducted using the CIGALE software, employing consistent input models and parameters for the entire sample. The data fluxes from the Cosmos2015 catalogue were used in the analysis, while the radio information as well as any nebular emission component were not taken into consideration. For the sources that have been identified spectroscopically in the literature (i.e., CVLA1002, CVLA958, ALESS09, and ECDFS544), the redshift input was fixed to the corresponding spectroscopic value. From the remaining galaxies, ALESS14 and CVLA1000, the redshifts were fixed to the photometric values retrieved from the literature, while for CVLA100 and CVLA028, the best-fitting SED was selected by introducing Δz = 0.1 intervals. The resulting fits, along with the model parameters used in this analysis, are presented in Fig. 3 and Table 2.

Fig. 3. CIGALE SED fitting results for the eight galaxies of our sample. The redshift values for CVLA1002, CVLA958, ALESS09, ECDFS544, ALESS14, and CVLA1000 were fixed to those derived from the literature, and the redshift for CVLA100 and CVLA028 was estimated from a range of values separated by Δz = 0.1.

To probe the significance of the AGN contribution, we applied the same fitting strategy using CIGALE without the AGN component. The updated SED fits revealed that for sources that initially required a lower AGN contribution (e.g., CVLA1000), the χ² values changed marginally by approximately ∼10%. On the other hand, for galaxies such as CVLA100, which presents a substantial AGN contribution, no accurate SED fit could be provided without an AGN component (Δχ² > 3). Consequently, this exercise demonstrates the necessity of incorporating an AGN component in the SED fit for five out of eight sources, generally leading to improved fits as compared to those obtained without an AGN component.

3.2. Sample presentation

COSMOSVLA J100256.53+021158.4. The COSMOS field deep multi-wavelength data provide the opportunity of examining the structure of this source, which appears to be composed of two galaxies: an optically bright quasar, and a fainter millimetre source (see the top left panel of Fig. 4). Spectroscopic data from the BOSS survey (Reid et al. 2016) secured the detection of the Ly-α emission line at a redshift of z = 3.503 ± 0.0007. Because the BOSS capabilities cannot resolve the two sources in study here, this archival spectrum corresponds to both galaxies. However, considering that the main source (i.e., the quasar) is the brightest galaxy at optical wavelengths, we assume that the emission line and subsequently derived redshift correspond to this source.

Fig. 4. Multi-wavelength view of the eight high-redshift galaxies of our sample, ordered according to Table 1. From left to right, we present the i- or r-band (HST, VIMOS, and VOICE data), the Spitzer 3.6 μm/VIDEO Ks-band/NICMOS 1.6 μm images, the radio VLA 1.4/3 GHz and the FIR data from ALMA and Herschel, and the red contours correspond to the VLA emission. From top to bottom, the contour levels of the radio image for each source are (7.715 × 10−6, 2.093 × 10−5), 2.827 × 10−5, (−4.710 × 10−6, 1.058 × 10−5), 5.640 × 10−6, (−4.845 × 10−6, 5.461 × 10−6), (−5.078 × 10−6, 1.407 × 10−5), (−4.633 × 10−6, 1.124 × 10−5, 2.712 × 10−5), and 2.092 × 10−5 Jy beam−1. Each panel has a size of 7.2 × 7.2 arcsec2, except for the Herschel image, which for visualisation purposes has a larger extent. Additionally, the beam size (shaded light grey area) of each instrument and a 2 arcsec scale bar are added in each panel. The ALMA data come from the projects AS2COSMOS and ALESS (Simpson et al. 2020; Hodge et al. 2013).

Furthermore, inspection of the 3 GHz data reveals that the quasar causes most of the radio flux density emitted from this system with a value of f_3 GHz-QSO = 0.052 ± 0.02 mJy (∼80% of the total), while the second source has a flux density of f_{3 GHz-2nd source} = 0.013 ± 0.009 mJy. When we assume a typical radio spectral index of 0.8 for both sources³, the corresponding radio flux density at 1.4 GHz is equal to f_{1.4 GHz-QSO} = 0.1 ± 0.04 mJy and f_{1.4 GHz-2nd source} = 0.024 ± 0.017 mJy. Therefore, if both galaxies are located at z ∼ 3.5 and the detected radio emission was generated solely from SF, the corresponding SFR (Bell 2003) would be SFR_radio-QSO = 4246 ± 1633 M_⊙ yr⁻¹ and SFR_{radio-2nd source} = 1061 ± 735 M_⊙ yr⁻¹ for the quasar and SFG, respectively (SFR_radio-total = 5307 ± 1790 M_⊙ yr⁻¹). According to Lanzuisi et al. (2017), this system of galaxies has a total stellar mass of log(M_⋆/M_⊙) = 11.69 ± 0.01 and SFR_IR-total of 1549 M_⊙ yr⁻¹. The latter value is derived from the far-IR flux integrated at the position of the quasar and is rather close to our SED-fitting, best-derived SFR estimate of 1721 M_⊙ yr⁻¹. Curiously, this value is similar to the value derived from the radio for the second source in the system (i.e., SFR_{radio-2nd source}). This is compatible with a system with a powerful quasar and a companion where most of the SF is taking place, and where no AGN is detected.

It has to be noted that for most bands, we cannot separate the brightness contribution from each of the two sources because of the sub-arcsec separation and the insufficient angular resolution of telescopes such as Spitzer. The SED fit was therefore performed using only total fluxes. A separation of the SED fitting analysis for the two sources would be non-trivial and would yield high uncertainties in the fluxes of each source. In the X-ray regime, COSMOS has been observed both with XMM-Newton (XMM-COSMOS survey; Cappelluti et al. 2009) and the Chandra X-ray Observatory (Chandra COSMOS Survey and Chandra COSMOS Legacy Survey; Elvis et al. 2009; Civano et al. 2016), reaching flux limits of f_{0.5 − 2 keV} ∼ 1.7 × 10⁻¹⁵ erg s⁻¹ cm⁻² and f_{0.5 − 2 keV} ∼ 2.7 × 10⁻¹⁶ erg s⁻¹ cm⁻², respectively. The source CVLA1002 has been observed and detected by the Chandra COSMOS survey with a flux of f_{2 − 10 keV} = 6.67 × 10⁻¹⁵ erg s⁻¹ cm⁻² (see Marchesi et al. 2016).

ALESS J033152.49−280319.1. Based on the current photometric data, this source seems to be a single galaxy, and according to the SED photometric analysis by da Cunha et al. (2015), it has a photometric redshift of z ∼ 3.38, an SFR_IR of 1202 M_⊙ yr⁻¹, and a stellar mass of 10^10.96 M_⊙. For our photometric analysis, which takes the possibility of an AGN component into account, the estimated best SFR value has a similar value of SFR_IR = 1188 M_⊙ yr⁻¹. Assuming that the galaxy is located at z = 3.38 and its radio emission is solely due to star formation, the corresponding SFR from the radio emission would be equal to SFR_radio ∼ 4213 ± 757 M_⊙ yr⁻¹. A comparison between this value and the value derived from our SED fitting analysis or by da Cunha et al. (2015) suggests that the source hosts intense SF, and an accreting SMBH explains the clear radio excess that is observed. The X-ray energies have been observed in the CDF-S field with XMM-Newton and the Chandra observatories (for more details, see Ranalli et al. 2013; Luo et al. 2017) in deep surveys, resulting in sensitivities of 6.6 × 10⁻¹⁶ erg s⁻¹ cm⁻² and 2.7 × 10⁻¹⁷ erg s⁻¹ cm⁻² at 2–10 keV, respectively. However, ALESS14 is not detected in any of these surveys, supporting the view that if an AGN exists, it contributes mostly to the radio emission, and at a much lower level at X-ray and optical/IR wavelengths. This was also suggested by the ∼20% AGN contribution to the IR emission from our CIGALE SED fit.

COSMOSVLA J100004.81+023045.2. According to the HELP catalogues (Shirley et al. 2021), this source is located at a redshift of 3.43 with an SFR_IR of 1519 M_⊙ yr⁻¹, while our SED fitting analysis predicts a lower SFR_IR of 690 M_⊙ yr⁻¹. Based on the SED fitting results, modelling of the IR emission of CVLA1000 requires an AGN contribution. According to the Cosmos2015 and Cosmos2020 (Weaver et al. 2022) catalogues, the photometric redshifts are 3.26 and 3.62, respectively, while Simpson et al. (2020) estimated a lower value of z  ∼  2.23. When we accept a redshift value of ∼ 3.3–3.6 for this source, then according to its radio flux density of f_1.4 GHz = 0.096 mJy, the SFR_radio should be ∼ 4000 M_⊙ yr⁻¹, which in comparison with the previous estimate from the IR emission strongly suggests an active SMBH at the centre of this galaxy. For the redshift solution of z = 2.23, the SED fit by Simpson et al. (2020) suggests S FRIR ∼ 135 M_⊙/yr, while the corresponding S FRradio is equal to ∼ 1545 M_⊙/yr. This value, although significantly lower than the 3723 M_⊙/yr estimate, still supports the presence of an AGN. Recent mm observations (C.-C. Chen et al. in prep) revealed a single spectral line that is compatible with the CO(3-2) transition at z = 2.48.

Due to its morphology in the near-IR, radio, and mm wavelengths, this source could be classified as a group of galaxies (see the third row of Fig. 4), however, catalogues with large samples of galaxies (i.e., Cosmos2015 and Cosmos2020) identify this galaxy as a single source due to the lack of resolution. With a radial extent of ∼1.146 arcsec, measured at the 1.6 μm – NICMOS data, and a scaling factor of 7.524 kpc ″⁻¹ (at z = 3.3), the actual size of this galaxy can be estimated as ∼8.6 kpc. This value is compatible with an ultra-compact dwarf galaxy (e.g., Rodriguez et al. 2021), however, the images reveal a source that resembles a disk rather than a spherical object. Inspection of all the available images cannot exclude the possibility of a merging system, considering that the NICMOS data depict a structure that could resemble a complex of galaxies. No detection in the X-ray bands is reported in the literature for this galaxy.

COSMOSVLA J095845.94+024329.2. Source CVLA958 was recently spectroscopically observed by our team with the Institut de Radioastronomie Millimetrique (IRAM) 30 m antenna (project code: 227-19; for more details regarding these observations, see Appendix B). Unfortunately, its mm faint nature prevented us from acquiring a secure redshift determination with the IRAM telescope. According to Cosmos2015, the photometric redshift of this galaxy is z = 3.172, while our photometric analysis suggests a minimum of χ² at 3.21. The latter redshift solution agrees with the only line that is detected from recent ALMA archival spectra (project code: 2019.1.01600.S). The detected emission line at ∼110.5 GHz has a full width at half maximum of ∼0.4 GHz (∼70.4 km s⁻¹), possibly indicating a merging system or an orderly rotating spiral galaxy (see Chen et al. 2022, for more details). The recent work by Birkin et al. (2021) also identified the emission line as due to the CO4-3 rotational transition at band 3 of ALMA, locating this source at z = 3.1755. According to our photometric SED fitting analysis, this source has a stellar mass of 10^11.07 M_⊙, an SFR of 1177 M_⊙ yr⁻¹, and a significant AGN contribution. Considering the measured flux density at 1.4 GHz of 90 μJy, we estimate that its SFR_radio (assuming that its radio emission comes purely from star formation) should be ∼ 3376 M_⊙ yr⁻¹. This value is significantly higher than the value determined from the SED fitting (see Table 3 for further details), supporting the presence of an AGN in this source as well. CVLA958 is not detected in any available dataset in the X-ray wavelengths.

COSMOSVLA J100233.16+020626.5. The source CVLA100 seems to be a single galaxy as well. Our photometric analysis places this source at a redshift of z_phot ∼ 3.3. This result agrees with the estimates of other photometric catalogues (see Laigle et al. 2016; Weaver et al. 2022). ALMA spectroscopic data acquired by our team (project code: 2017.1.01713.S; see Appendix B for more details) reveal an emission line at ∼107.3 GHz, which corresponds to the CO(4–3) emission line with the galaxy located at z = 3.29. For the derived physical properties of this source, our SED fitting analysis indicates a stellar mass of 10^10.27 M_⊙, an SFR of 566 M_⊙ yr⁻¹, and an AGN fraction of 60%. Application of Eq. (6) by Bell (2003) provides us with an estimate of the SFR_radio of 2760 M_⊙ yr⁻¹. The SED fitting models for this redshift suggest an SFR of ∼ 500 M_⊙ yr⁻¹, which could mean that the observed 1.4 GHz radio flux cannot be obtained solely from SF, but by a combination of star emission and AGN activity (see Table 3 for additional information). Similarly to CVLA958, this source is not detected at any high-energy band.

COSMOSVLA J100028.75+023203.5. As mentioned in the previous section, this source (which seems to be a single galaxy) has a spectroscopic redshift of 4.76 (see Hasinger et al. 2018), however, inspection of the DEIMOS spectrum does not reveal a secure emission line (i.e., the Ly-α line, according to the authors). Chen et al. (2022) reported the detection of an emission line at v ∼ 106.1 GHz that they interpreted as the CO(3–2) transition, placing the galaxy at z = 2.26. Alternatively, this line could be the CO(4–3) transition that places the source at z = 3.35. Our own photometric SED fitting analysis suggests a redshift of 3.27, which is consistent with the redshift provided by the HELP database and Miettinen et al. (2015). When we assume this photo-z to be the true redshift of this galaxy, our SED fit analysis proposes a stellar mass of 10^10.74 M_⊙ yr⁻¹, an SFR of 2554 M_⊙ yr⁻¹, and no AGN contribution. The derived SFR from the radio emission of this galaxy (i.e., 2166 M_⊙ yr⁻¹) further supports the argument that no radio emission is added by SMBH activity. This result establishes CVLA028 as the only galaxy in our sample without an indication of an AGN by the SED fit or radio excess. Although this galaxy has been observed by both XMM-Newton and Chandra, the images at all bands reveal no clear detection.

ALESS J033211.34−275211.9. This source is recognised in the literature as an X-ray AGN with faint radio emission. Photometric analysis from various deep surveys (e.g., Guo et al. 2020) reveals a redshift of z ∼ 3.7, while ALMA spectroscopic data (see Birkin et al. 2021) show two clear emission lines (namely the CO4-3 and [CI] lines at ∼98.22 GHz and ∼104.85 GHz, respectively). This secures the galaxy at z = 3.694. Assuming that this spectroscopic redshift and that a radio flux density of 0.043 mJy solely originate from star formation, the corresponding SFR is ∼ 1884 M_⊙ yr⁻¹. However, according to Guo et al. (2020), following the SED fitting method using CIGALE, this galaxy is a type II AGN with a stellar mass of ∼10^11.3 M_⊙, an SFR of ∼ 31 M_⊙ yr⁻¹ (our CIGALE SED fitting analysis suggests a higher value of 1194 M_⊙ yr⁻¹), and an AGN fraction of 72% (our analysis proposes a lower value of 10%). The source detected in the HST image (see Fig. 4), which is rather close to the ALMA detection (∼0.5 arcsec), seems to cause the radio and IR emission (see Wardlow et al. 2011). Therefore, this source could be a complex system of at least two galaxies possibly interacting, with the presence of a radio-quiet AGN and an SFG. Even though ALESS09 is not detected in any of the X-ray surveys, another source lies at a distance of ∼2 arcsec (named LBX2017−14; see Ranalli et al. 2013, for more details) that seems to be an X-ray bright AGN with a flux of f_{2 − 10 keV} = 6.26 × 10⁻¹⁵ erg s⁻¹ cm⁻² and no counterparts. This X-ray source was spectroscopically identified by Ranalli et al. (2013) at z = 3.74 ± 0.06, meaning that it could be part of a merging system along with the sub-millimetre and radio galaxy identified with the source ALESS09.

[DSS2017] 554. This source is present in the VIMOS spectroscopic survey presented by Danielson et al. (2017) with a spectroscopic redshift of 3.1988, which was used in this work as a fixed value for the subsequent SED fitting analysis. This galaxy is placed at the outskirts of the ECD-S field, meaning that it was not selected in the photometric analysis by Guo et al. (2020). Inspection of the VST Optical Imaging of the CDFS and ES1 (VOICE; Vaccari et al. 2016) data for this source in bands r, g, and i (see the bottom row of Fig. 4) suggests the possibility of multiple components. These observations might signify a merging or interacting system. The predicted SFR derived solely from the radio emission is 1520 M_⊙ yr⁻¹, while our SED fitting analysis suggests a value of 1162 M_⊙ yr⁻¹ and an AGN contribution of 20%. No clear detection has been reported regarding the possible X-ray nature of this source.

4. Discussion

We have examined a sample of eight high-redshift dusty galaxies with an AGN component (3 < z_spec, phot < 4). The selection process involved using the 1.4 GHz VLA radio catalogues from the COSMOS and ECDF-S fields in combination with the far-IR peak criterion (f_350 μm/f_250 μm > 1) using the HerMES survey. Moreover, a compact morphology was considered as an additional criterion. In the following paragraphs, we assess the efficacy of these criteria in selecting these sources. Furthermore, we provide a comparison between the derived properties of our sample with those obtained in previous studies, focusing on similar redshift ranges.

4.1. Validity of the selection criteria

While the compactness criterion merely serves as a method for avoiding local extended galaxies, the far-IR and radio criteria have been employed extensively in the past in the study of the distant Universe. For instance, Amblard et al. (2010) demonstrated that the redshift of a source tends to increase for galaxies that satisfy the criteria f_350 μm/f_250 μm > 1 and f_500 μm/f_350 μm > 1. These results are expected, considering that the far-IR emission of an SED traces the dust content of a galaxy, and more specifically, its blackbody radiation, which is shifted to higher wavelengths with increasing redshift. However, this emission is affected by the temperature of the dust as well, which consequently introduces a degeneracy that is non-trivial to break. Another point of attention is the sensitivity of the Herschel surveys, which are inherently biased towards the detection of the brightest dust-rich far-IR galaxies.

In order to evaluate this argument and selection method further, we employed the HELP catalogue, which contains far-IR observational data (i.e., 15 747 Herschel detected galaxies) along with SED fitting results from CIGALE (for a detailed description of the CIGALE parameters used in this catalogue, see Małek et al. 2018). The galaxies from this catalogue are presented in the left panel of Fig. 5, in which the blue and orange histograms correspond to the sources that follow the criterion f_350 μm/f_250 μm < 1 and f_350 μm/f_250 μm > 1, respectively. As anticipated, galaxies that adhere to the far-IR criterion are situated at higher redshifts on average. Interestingly, based on the CIGALE fits, a higher percentage (i.e., ∼45%) of the galaxies that follow the f_350 μm/f_250 μm ≥ 1 criterion exhibit SFRs of SFR > 10² M_⊙ yr⁻¹. In contrast, for the rest of the sources, a smaller percentage (∼22%) demonstrates high SFR values like this.

Fig. 5. Normalised redshift histograms of all galaxies from the HELP catalogue (panel a), separated with the criterion f350 μm/f250 μm = 1. Flux densities were estimated using the Herschel observational data as provided from the same catalogue. The vertically dotted coloured lines correspond to the average values of the same histograms, employing the flux densities from the CIGALE best model fit, and are coloured accordingly. The text in the panel denotes the percentage of the sources that present an SFR above 102 M⊙ yr−1 and are coloured according to the aforementioned criterion. Panel b illustrates the normalised histograms for the dust content of the galaxies in both the HELP and Cosmos2015 catalogues. The blue and orange histograms are derived from galaxies without and with radio VLA 1.4 GHz emission, respectively. The text of the panel demonstrates the percentage of the galaxies that present an AGN component (AGNfrac > 40%) and are placed at z > 3.

A similar exercise was repeated for the radio criterion by estimating the AGN fraction and dust mass, derived from the CIGALE fits, for each galaxy. Because the HELP catalogues do not consider the radio emission in their fits, we cross-matched this catalogue with Cosmos2015, which includes the 1.4 GHz VLA data. The results of this analysis are displayed in the right panel of Fig. 5, which shows histograms of the dust mass. The shape and mean values of these histograms indicate that the radio selection criterion does not guarantee the preferential selection of highly dusty SFGs. However, upon examining galaxies at redshifts above 3, it becomes evident that the radio selection method identifies a higher proportion of sources (i.e., ∼22 ± 9%) with a substantial AGN fraction (i.e., AGN_frac > 40%) compared to cases without detected radio emission (∼15 ± 4%).

While this percentage difference may not be statistically significant, the combination of the two far-IR and radio selection methods still provides motivation for their application in identifying high-redshift SFGs with an AGN component. The detailed analysis of the eight sources in this sample has reinforced these findings and encourages their application in large-scale surveys, such as EMU.

4.2. Comparison with previous work

Overall, our sample has a mean redshift of ⟨z⟩ = 3.37 ± 0.16, an average stellar mass of ⟨log(M_⋆/M_⊙)⟩ = 11.27 ± 0.47, ⟨L_{IR, 42 − 122 μm}⟩ = 2.1 ± 0.67 × 10³⁹ W (total IR luminosity of ⟨L_{IR, 8 − 1000 μm}⟩ = 3.9 ± 0.84 × 10³⁹ W ∼ 10¹³ L_⊙), ⟨L_1.4 GHz⟩ = 5.77 ± 1.9 × 10²⁴ W Hz⁻¹, and a mean star formation rate (derived from the 42 − 122 μm rest-frame IR luminosity of each source) of ⟨SFR_IR⟩ 1213 ± 567 M_⊙ yr⁻¹. This average SFR is consistent with previous results that reported an increasing SFR at higher redshift. For instance, Rowan-Robinson et al. (2016) and Gao et al. (2022) showed that Herschel-detected samples of hyper-luminous infrared galaxies (HLIRGs; L ∼ 10¹² − 10¹³ L_⊙) at z > 3 have SFR values ranging between ∼ 000 − 10000 M_⊙ yr⁻¹. Therefore, the inclusion of the 500 μm emission in our selection criteria naturally provides us with the high bright tail of the SFR distribution at these redshifts. Similarly, Wang et al. (2019) reported that the selection of Herschel galaxies above z > 3 revealed dusty SFR galaxies with infrared luminosities of ∼10¹² − 10¹³ L_⊙, similar to the derived values of our sample.

Furthermore, for 87.5% of the sources in our sample (seven out of eight), an AGN is indicated either by the radio excess or by the mid-IR signature obtained from the best SED fit. Although the morphology of these sources still remains uncertain, five out of eight (63%) might consist of multiple galaxies, with significant levels of both SF and AGN activity. In comparison, Gao et al. (2021) and Wang et al. (2021), who analysed a sample of HLIRGs selected using data from Herschel and the Deep Fields of the LOFAR Two-metre Sky Survey (LOTSS; Shimwell et al. 2017). described an AGN contamination of 34–64% (using six different models to fit the AGN component) for galaxies at the redshift range 3 < z < 4. These values, although model dependent, seem to be lower than the AGN contribution of our sample, which can be attributed to the fact that the LOTSS survey reaches sensitivities of 20–30 μJy beam⁻¹ at 150 MHz. When we assume a typical radio synchrotron slope, this value translates into a lower flux density value at 1.4 GHz than the COSMOS-VLA sensitivity (i.e., 7 μJy beam⁻¹). Therefore, it seems that by selecting radio sources from a deeper radio survey, the authors revealed a population of SF galaxies that receive contributions from AGN less frequently. However, it is important to emphasise that a comparison of our results with studies involving larger sample sizes should be approached with caution, considering the limited statistics that our data provides.

Additionally, Enia et al. (2022) reported a negligible percentage of AGN for the entire sample (0 < z < 6), which was selected using deep radio 1.4 GHz observations from the GOODS-N field (sensitivity of 2.2 μJy). This can be perceived with the parameter q_TIR = log(L_IR/3.75 × 10¹²)−log L_1.4 GHz, which is commonly used as an indicator of the presence of AGN contribution to the radio emission. Figure 6 depicts the q_TIR values for our sample and those of Enia et al. (2022), demonstrating the AGN contamination in the radio emission in most of our galaxies. In more detail, 62.5% of the sources in our sample have lower values of q_TIR than the grey region of the plot (which denotes the average redshift-dependent q_TIR). These statistics do not fully match the number of 87.5% of AGN contamination that was derived from the SED fitting and described earlier in the text. However, this discrepancy might originate from the low AGN contribution and uncertainties in the measurements of the radio and infrared fluxes of the SED of some galaxies. In contrast, the sample by Enia et al. (2022) for the redshift bin of 3 < z < 4 indicates an AGN contamination of 37.5%. This difference can be attributed to the fact that in their selection criteria, they did not require the far-IR emission (or better, far-IR peak spectrum) that we employed. As a result, only 38.4% of their total sample has a 500 μm Herschel detection.

Fig. 6. The qTIR parameter over redshift for the sample explored by Enia et al. (2022) and for our 8 galaxies (panel a). This figure was retrieved by Enia et al. (2022), published under the CC BY 4.0 international licence and was modified by adding the results from our work (square points). The dashed black line and grey shaded region correspond to the best-fit relation of their total sample and the standard deviation, respectively, while the dashed blue line shows a similar relation by Delhaize et al. (2017). All points by Enia et al. (2022) are coloured according to the contribution of AGN to the radio emission defined as $f_{\mathrm{AGN}}=1-{10}^{q_{\mathrm{tot}}-\overline{q_{\mathrm{SF}}}}$, with qtot and $\overline{q_{\mathrm{SF}}}$ corresponding to the qTIR and its evolution with redshift, respectively. The star symbols depict what is defined by Enia et al. (2022) as H-dark sources (i.e., not detected at the Hubble-WFC3 H-band). Panel b, presents the SFRs derived from the IR and radio luminosity, for the 8 sources of our sample, along with the results by Bonzini et al. (2015) (radio-quiet AGN: blue squares, radio-loud AGN: red triangles, SF galaxies: green circles) and Talia et al. (2021) (average values of their sample). All data points with redshift below 3 are coloured with transparency, while the ones with z > 3 have normal colours.

Another comparison can be conducted with the q_TIR relation derived by Delvecchio et al. (2021) and presents a mass-dependent and almost redshift-independent equation (i.e., q_TIR(M_⋆, z) = (2.646 ± 0.024)(1 + z)^{−0.023 ± 0.008} − (log M_⋆/M_⊙ − 10)(0.148 ± 0.013)). According to this relation, all of the galaxies in our sample satisfy the criterion, placing them below the mean of the SFR population. However, none of the sources in our sample follows the relation log(L_1.4 GHz/SFR_IR) = 21.984 × (1 + z)^0.013 that was derived earlier by Delvecchio et al. (2017) and corresponds to the 3σ deviation from the peak distribution of the ratio of L_1.4 GHz and SFR_IR as a function of redshift.

In the same manner, Novak et al. (2017), which selected sources from the VLA 3 GHz survey without a far-IR emission requirement, similarly presented the trend of a lower contamination of AGN (23% and 47% of the total sample and 3 < z < 4 redshift bin, respectively) than in our sample. This comparison favours the inclusion of the criterion of a far-IR rising spectrum to the radio emission when AGN at high redshifts are to be identified and studied. In contrast, an interesting comparison can be performed with studies such as that by Bonzini et al. (2015), in which the authors presented a sample with galactic properties similar to those of our galaxies. This sample was derived from the radio 1.4 GHz ECDF-S catalogues (6 μJy sensitivity) with the addition of the Herschel detection in all sources at the redshift bin of interest (i.e., z = 3 − 4). Fifteen of the 28 galaxies presented in this redshift range were categorised as radio-quiet (RQ) AGN while one as a radio-loud (RL) AGN, totalling 16 galaxies, or an AGN contamination of 73%. This value is close to the value derived from our sample, which further supports the necessity of including the far-IR criterion to detect high-z AGN. Additionally, the right panel of Fig. 6 shows that our galaxies fall in the RQ region of their work, even though our estimates indicate that at least two of our sources should have a significant AGN contribution.

In this sense, our selection method might provide an increased probability for the detection of dusty galaxies at high-z that present high SFRs with an additional AGN component, with weak radio emission, as part of a merging system. Additionally, our galaxies seem to be a more extreme group of interacting systems, having brighter IR and radio luminosities as well as stellar masses and SFRs as compared, for instance, to the work by Talia et al. (2021). This outcome might originate from the fact that their sample, once again, did not require a far-IR Herschel detection and was originally derived from the 3 GHz – COSMOS survey, which has a higher detection limit than the VLA-1.4 GHz-COSMOS survey we employed. Nevertheless, Novak et al. (2017) also explored the same field at 3 GHz and reached the same flux density sensitivity (rms of 2.3 μJy), from which they achieved more AGN detections. The main difference lies in their selection criteria, with the latter work requiring an optical/near-IR counterpart, while Talia et al. (2021) excluded these sources. The selection of radio sources with a far-IR peak spectrum and an optical counterpart might therefore reveal merging systems that at high z, present high SFRs (∼ 10³ M_⊙ yr⁻¹) and an AGN counterpart in most cases.

Nevertheless, selecting galaxies that are considered as extreme in terms of mass and luminosity will naturally limit the study to a rare sample, hence to a reduced number of selected sources. This infrequency, along with the complexity of the adopted selection method, will unavoidably provide a less accurate estimation of the SFRD at z > 3, which can be only considered as a lower limit (following the V_max method introduced by Avni & Bahcall 1980), and for our study, this is equal to SFRD > 8 × 10⁻³ M_⊙ yr⁻¹ Mpc⁻³. However, these selection criteria might be quite powerful in the aforementioned identification and in studying the interplay between AGN and SF activity at high redshift, especially with their application to the upcoming next-generation radio surveys (e.g., the EMU survey) and current/future observational data from mm facilities (e.g., Bing et al. 2023).

5. Conclusions

In this work, we present the exploration of a sample of a high-redshift (3 < z < 4) radio galaxy population that was selected with the criterion of a rising far-IR spectrum (i.e., f_350 μm > f_250 μm) along with the detection of non-extended (< 10 arcsec) radio emission at 1.4 GHz. The aim behind this selection method was to investigate the potential of combining these two criteria in identifying high-redshift dusty SFGs that are contaminated by an AGN and might be ‘hidden’ by other selection methods. A population like this, due to its extremity and low occurrence, might not contribute significantly to the SFRD of the Universe, but it can provide us with further insight into the interplay between SMBHs and SF activity at high z.

In order to explore the morphology of our sample and its physical characteristics in depth, we took advantage of the multi-wavelength deep data available in the literature from the COSMOS and ECDF-S fields. These data sets enabled us to identify that five out of eight galaxies in our sample might consist of more than one galaxy with possible interaction/merging. It has to be noted that the low angular extent of these galaxies did not permit confirmation of this statement due to the spatial resolution of the current instrumentation. Nevertheless, future observations with telescopes such as the JWST may provide further insights into this argument. Furthermore, we determined the presence of an AGN component in seven out of eight of these sources, which in most cases represents a weak radio emission, identifying these sources as radio-quiet AGN. An additional highly star-forming component seems to be present in all cases, with a typical SFR of ∼ 10³ M_⊙ yr⁻¹.

We have demonstrated that the combination of the far-IR and radio emission criteria is highly efficient in selecting merging systems at 3 < z < 4, typically with the presence of a radio-quiet AGN and a highly SFG. With upcoming wide radio surveys (e.g., the Evolutionary Map of the Universe, EMU) and next-generation telescopes such as the SKA, these results motivate further exploration of alternative selection methods. In order to advance the study of systems similar to those in our sample, particularly during the epoch of reionisation (EoR), sources with a mm rise in flux in combination with deep radio surveys might further be studied. A possible rising mm spectrum profile, which could be explored for faint individual sources using interferometric facilities (e.g., ALMA and NOEMA) and single-dish instruments (e.g., the NIKA2 camera of the IRAM 30 m antenna, or the ToLTEC of the Large Millimeter Telescope) targeting a brighter population from larger surveys could assist in the identification of the redshifted dust component of SFGs at z > 6. Similarly, the morphology of these galaxies could be further understood with infrared facilities, including the JWST.

Acknowledgments

This work was supported by Fundação para a Ciência e a Tecnologia (FCT) through the research grants PTDC/FIS-AST/29245/2017, 1220UID/FIS/04434/2019, UIDB/04434/2020, and UIDP/04434/2020. C.P. acknowledges support from DL 57/2016 (P2460) from the ‘Departamento de Física, Faculdade de Ciências da Universidade de Lisboa’. This paper makes use of ALMA data. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan); together with NRC (Canada), MOST, and ASIAA (Taiwan); and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. This research made extensive use of core PYTHON packages for Astronomy as well as TOPCAT. We thank the anonymous referee for all the valuable comments, which improved the quality of the paper.

References

Appendix A: Comparison of the retrieved z with photometric catalogues

In this appendix, we present a comparison between the redshift values estimated by various photometric catalogues and our analysis, as described above. In more detail, the catalogues consist of the Cosmos2015 (Laigle et al. 2016), the Cosmos2020 (Weaver et al. 2022, employing all possible redshift estimations, namely the classic and Farmer methods using the SED fitting codes Lephare and EAZY) and the HELP catalogues (Shirley et al. 2021). The redshift values derived from these catalogues are compared with the values that our team acquired from the literature and from the IRAM/ALMA observations of the aforementioned two targets (i.e. CVLA958 and CVLA100; see Sect. 3). For the cross match between catalogues, a radius of 1 arcsec was employed in order to ensure that each source from the catalogues corresponded to the correct source from our sample. The results of these comparisons are presented in Fig. A.1.

Fig. A.1. Comparison betthe orange points depict the spectroscopic values from the literature. Notably, the Cosmos2020 catalogue yields two distinct sets of results beween the Cosmos2015 (a), Cosmos2020 (b and c), and HELP (d) catalogues and our redshift values for our sample. The blue points in each plot correspond to the photometric z values from each catalogue, and cause two different SED fitting codes were used: Lephare and EAZY, denoted by circles and triangles, respectively.

In our findings, certain sources exhibit lower redshift values in the catalogues than in the present study. These differences can be attributed to the low S/N of the photometric points that were used in the SED fitting. Furthermore, the authors of these catalogues may not have individually examined each source due to the extensive volume of data. For example, the source with coordinates (Ra, Dec) = (150.02005417, 2.51254722) is depicted in the Cosmos2020 comparison plots (see Fig. A.1) with a redshift of z = 2.48, derived from the EAZY SED fit (χ² of ∼15). However, our work adopts a redshift value of 3.46 obtained from the Lephare SED fit solution for the same source. This choice was made due to the lower χ² values obtained from both the literature Lephare analysis (χ² = 3.76) and our own CIGALE fits (χ² = 2.24).

Additionally, this figure shows that according to Cosmos2015, two galaxies that were identified by our search as z < 3 sources have a higher photometric redshift value. In more detail, the sources are identified as UVISTADR1 J100131.33+020325.2 and UVISTADR1 J100114.75+023516.3 and have a z_{phot, cosmos2015} of 3.152 and 3.007, respectively. However, the more detailed photometric analysis by Muzzin et al. (2013) proposed lower values for both sources of ∼1.8 and ∼2.7, respectively. Similarly, for the cosmos2020-Farmer method, a source named UVISTADR1 J095834.49+015504.9 was estimated to be at z_{phot, cosmos2020, EAZY} ∼ 3.1, however, Muzzin et al. (2013) reported a lower value of ∼2.5. The latter redshift results can be further supported by the redshift estimate by Lephare cosmos2020 of z_{phot, cosmos2020, Lephare} ∼ 2.2. These discrepancies between various photometric catalogues can mainly be attributed to the faint nature of these sources and to the subsequent low fit quality of their SED. Finally, according to the HELP catalogue, no further source could be added to our z > 3 sample.

Appendix B: Spectroscopic follow-up of CVLA958 and CVLA100

Because of the uncertainty of the high-redshift nature of CVLA958 and CVLA100, these sources were originally included in a follow-up spectroscopic study by our team with the aim to understand their nature and potentially employ this knowledge to detect galaxies at higher redshifts. In this appendix, we present the retrieved spectroscopic data for these two galaxies, bearing in mind the redshift solution of z ∼ 3 that the current photometric analysis indicated. The original idea of the spectra acquired with the Institute Radio Astronomie Millimétrique (IRAM) 30m antenna and the ALMA telescopes was the detection of two possible CO emission lines that would constrain the redshift of each galaxy.

B.1. IRAM observations

The IRAM⁴ institute is composed of two facilities, namely the interferometer NOrthern Extended Millimeter Array (NOEMA; Chenu et al. 2016) and a single-dish telescope with a diameter of 30m, located in the Sierra Nevada, Spain. The adopted strategy to retrieve the millimetre spectrum of our galaxies was to scan the 3 mm band with the spectral line receiver Eight MIxer Receiver (EMIR; Carter et al. 2012) of the 30m antenna, along with the Fast Fourier Transform Spectrometers (FTS200; Klein et al. 2012) and Wideband Line Multiple Autocorrelator (WILMA⁵) backends (in parallel) at the wavelength range of 73-117 GHz (band E0/E090). Given the possible AGN nature of the targeted sources, the CO spectral line energy distribution will likely peak at mid-J transitions (e.g. for the frequency range of band E0, the CO J = 4-3 line could be at z ∼ 2.7 − 5 and the CO J = 5-4 line at z ∼ 3.7 − 6.5). The expected redshift span for this candidate isz ≥ 3, and therefore, the 3 mm spectral range provides an excellent trade-off of weather versus line detection to search for ideally two emission lines required to determine the source redshift. Nevertheless, in the following analysis, only the data from the FTS units were investigated since they provide a higher spectral resolution than WILMA, which is necessary to resolve faint emission lines.

The project code is 227-19, and the observations were initially scheduled for 19-24 December 2019. However, considering that the weather conditions were not appropriate to observe for the full length of our original allocated time, additional observing time was provided during 21-23 March 2020, 12-14 May 2020, and finally, 14-17 July 2020. Out of a total of 103 hours provided for this source, only ∼8 hours of observations generated spectra with acceptable quality to be considered in the assembly of the final spectrum. The rest of the observing time was lost due to instrumentation malfunctions and weather conditions (strong wind and high water vapour).

The final spectrum of CVLA958 was constructed using only data from the last three days of observations (i.e. 21, 22, and 23 December 2019). We further inspected the histograms of the opacity of the clouds (τ), rms random noise temperature (T_rms), and receiver and system temperature (T_rec and T_sys, respectively) from these observations in order to exclude possible outliers in the scans. After removing these outliers and correcting for the platforming issue (e.g. Klein et al. 2008, 2009), we generated the final spectrum for the 3 mm band, depicted in the left panel of Fig. B.1. In this plot, the total spectrum (blue line) and its filtered version for a Gaussian filter of 50 km/s (orange line) are presented, and the colours are denoted in the legend. The full applied FTS configuration is also shown, as is the location of the CO lines, accepting the recent spectroscopic redshift value of 3.1755 by Birkin et al. (2021) (dashed red lines). The figure clearly shows that no secure line can be attributed to any feature because the S/N does not permit this determination.

Fig. B.1. Reduced millimetre spectra of CVLA958 and CVLA100 observed by the IRAM 30m and ALMA telescopes. The IRAM spectrum for CVLA958 is shown, with its unbinned reduced spectrum presented in blue, and the orange line corresponds to the Gaussian filtered spectrum using a width of 50 km/s (left panel). The CO lines for the redshift solution of z = 3.1755 provided by Birkin et al. (2021) are depicted with the dashed red lines. The applied FTS setups are presented at the top of the plot and are coloured according to the adopted configuration and observing date. The ALMA spectrum for CVLA100 at band 3 is shown in the right panel. The vertical black line corresponds to the possible CO(4-3) emission line at z = 3.29, located at ν = 107.3 GHz. The inset cutout illustrates the velocity-integrated map of the line emission, and the Spitzer 3.6 μm contour (white isocontour) is overlaid.

The high-redshift candidate CVLA100 was also observed by our team with the IRAM 30m antenna as part of the project 198-16 and during the time period of 23-27 March 2017. In these observations, we also scanned the E0 band of EMIR (frequency range of 73–117 GHz) using the FTS200 backend. Additional time was allocated for our source, and observations were conducted remotely on 26 May 2107, aiming at the 2 mm spectral window. Similarly, with the acquired spectrum of CVLA958, no secure emission line could be detected for CVLA100. In order to further test possible lines that are not apparent with visual inspection, we developed a simple convolution algorithm to search for emission from CO lines at a given redshift. Essentially, the approach of this code is to sum the observed spectra at the positions where a CO detection was expected at a given redshift, and subsequently vary the redshift in steps of 0.01. Application of this analysis, presents no emission line detection, considering that the resulting values for specific redshifts ranging from z = 0 to z = 6 do not deviate significantly from zero.

B.2. ALMA observations

No emission was detected in the IRAM spectrum for CVLA100, our team decided to secure additional observing time for the same candidate using the ALMA telescope. More specifically, the source was observed with multiple tunings in bands 3 (84-116 GHz) and 4 (125-163 GHz) in various configurations (longest baselines from 300 m to 1400 m), under the project code 2017.1.01713.S. With acceptable weather conditions (water vapour of ≤3 mm at all times; wind at ∼3 m/s), our observations achieved a mean angular resolution of ∼1 arcsec and a mean continuum sensitivity of ∼0.2 mJy after 4.1 hours at band 3 and 1.4 hours of integration time at band 4. By analysing the spectrum, extracted using the CASA software (CASA Team et al. 2022), we detected a feature at 107.3 GHz which, considering the numerous CO emission lines, assisted in constraining the redshift of the source. One plausible solution indicates that this feature corresponds to the redshifted CO(4-3) emission line at z = 3.3. This proposed redshift possibility closely agrees with the suggested solution of z = 3.15 obtained from our photometric analysis and the referenced literature catalogues (i.e. HELP, Cosmos2015, and cosmos2020).

All Tables

All Figures

Fig. 1. Flux density ratios of 350 and 250 μm vs. the ratios of 500 and 350 μm of HerMES in COSMOS (panel a) and ECDF-S (panel b). The sample was selected by cross-matching the HerMES catalogue with the VLA-COSMOS and VLA-ECDFS data. The dashed black lines indicate f350/f250 = 1 and f500/f350 = 1, while the points surpassing the limit of f350/f250 > 1 are presented with larger markers. The colour bar corresponds to the 350 μm flux of each source, and the legend shows the number of sources for each panel and the selected high-redshift candidates (crosses). The black lines in the subregion of each panel provide information regarding the shape of the far-IR (FIR) emission. For instance, sources with f500/f350 < 1 and f350/f250 > 1 have a peak in flux between 250 and 500 μm. The photometric properties of the sources presented in the legend are provided in Table 1.

In the text

Fig. 2. Histograms presenting the number of sources for each redshift bin selected by the criteria for the COSMOS (empty histograms) and ECDF-S (filled histograms). The blue histograms correspond to spectroscopic values derived from the literature, and the orange histograms show the photometric values. The combined sample is depicted with an empty black histogram. The number of sources for each category is presented in the legend.

In the text

	Fig. 3. CIGALE SED fitting results for the eight galaxies of our sample. The redshift values for CVLA1002, CVLA958, ALESS09, ECDFS544, ALESS14, and CVLA1000 were fixed to those derived from the literature, and the redshift for CVLA100 and CVLA028 was estimated from a range of values separated by Δz = 0.1.
In the text

Fig. 4. Multi-wavelength view of the eight high-redshift galaxies of our sample, ordered according to Table 1. From left to right, we present the i- or r-band (HST, VIMOS, and VOICE data), the Spitzer 3.6 μm/VIDEO Ks-band/NICMOS 1.6 μm images, the radio VLA 1.4/3 GHz and the FIR data from ALMA and Herschel, and the red contours correspond to the VLA emission. From top to bottom, the contour levels of the radio image for each source are (7.715 × 10−6, 2.093 × 10−5), 2.827 × 10−5, (−4.710 × 10−6, 1.058 × 10−5), 5.640 × 10−6, (−4.845 × 10−6, 5.461 × 10−6), (−5.078 × 10−6, 1.407 × 10−5), (−4.633 × 10−6, 1.124 × 10−5, 2.712 × 10−5), and 2.092 × 10−5 Jy beam−1. Each panel has a size of 7.2 × 7.2 arcsec2, except for the Herschel image, which for visualisation purposes has a larger extent. Additionally, the beam size (shaded light grey area) of each instrument and a 2 arcsec scale bar are added in each panel. The ALMA data come from the projects AS2COSMOS and ALESS (Simpson et al. 2020; Hodge et al. 2013).

In the text

Fig. 5. Normalised redshift histograms of all galaxies from the HELP catalogue (panel a), separated with the criterion f350 μm/f250 μm = 1. Flux densities were estimated using the Herschel observational data as provided from the same catalogue. The vertically dotted coloured lines correspond to the average values of the same histograms, employing the flux densities from the CIGALE best model fit, and are coloured accordingly. The text in the panel denotes the percentage of the sources that present an SFR above 102 M⊙ yr−1 and are coloured according to the aforementioned criterion. Panel b illustrates the normalised histograms for the dust content of the galaxies in both the HELP and Cosmos2015 catalogues. The blue and orange histograms are derived from galaxies without and with radio VLA 1.4 GHz emission, respectively. The text of the panel demonstrates the percentage of the galaxies that present an AGN component (AGNfrac > 40%) and are placed at z > 3.

In the text

Fig. 6. The qTIR parameter over redshift for the sample explored by Enia et al. (2022) and for our 8 galaxies (panel a). This figure was retrieved by Enia et al. (2022), published under the CC BY 4.0 international licence and was modified by adding the results from our work (square points). The dashed black line and grey shaded region correspond to the best-fit relation of their total sample and the standard deviation, respectively, while the dashed blue line shows a similar relation by Delhaize et al. (2017). All points by Enia et al. (2022) are coloured according to the contribution of AGN to the radio emission defined as $f_{\mathrm{AGN}}=1-{10}^{q_{\mathrm{tot}}-\overline{q_{\mathrm{SF}}}}$, with qtot and $\overline{q_{\mathrm{SF}}}$ corresponding to the qTIR and its evolution with redshift, respectively. The star symbols depict what is defined by Enia et al. (2022) as H-dark sources (i.e., not detected at the Hubble-WFC3 H-band). Panel b, presents the SFRs derived from the IR and radio luminosity, for the 8 sources of our sample, along with the results by Bonzini et al. (2015) (radio-quiet AGN: blue squares, radio-loud AGN: red triangles, SF galaxies: green circles) and Talia et al. (2021) (average values of their sample). All data points with redshift below 3 are coloured with transparency, while the ones with z > 3 have normal colours.

In the text

Fig. A.1. Comparison betthe orange points depict the spectroscopic values from the literature. Notably, the Cosmos2020 catalogue yields two distinct sets of results beween the Cosmos2015 (a), Cosmos2020 (b and c), and HELP (d) catalogues and our redshift values for our sample. The blue points in each plot correspond to the photometric z values from each catalogue, and cause two different SED fitting codes were used: Lephare and EAZY, denoted by circles and triangles, respectively.

In the text

Fig. B.1. Reduced millimetre spectra of CVLA958 and CVLA100 observed by the IRAM 30m and ALMA telescopes. The IRAM spectrum for CVLA958 is shown, with its unbinned reduced spectrum presented in blue, and the orange line corresponds to the Gaussian filtered spectrum using a width of 50 km/s (left panel). The CO lines for the redshift solution of z = 3.1755 provided by Birkin et al. (2021) are depicted with the dashed red lines. The applied FTS setups are presented at the top of the plot and are coloured according to the adopted configuration and observing date. The ALMA spectrum for CVLA100 at band 3 is shown in the right panel. The vertical black line corresponds to the possible CO(4-3) emission line at z = 3.29, located at ν = 107.3 GHz. The inset cutout illustrates the velocity-integrated map of the line emission, and the Spitzer 3.6 μm contour (white isocontour) is overlaid.

In the text