A TALE OF TWO FEEDBACKS: STAR FORMATION IN THE HOST GALAXIES OF RADIO AGNs

The discovery of a number of scaling relations (M_BH–σ: Ferrarese & Merritt 2000; Gebhardt et al. 2000; Merritt & Ferrarese 2001; Tremaine et al. 2002; M_BH–M_bulge: Magorrian et al. 1998; McLure & Dunlop 2002; and M_BH–L_bulge: Kormendy & Richstone 1995; Marconi & Hunt 2003) connecting the properties of the nuclear regions of galaxies to their global, host galaxy, characteristics has led to a still ongoing debate about the possible physical processes that could give rise to this connection. This link has been argued to be either causal, coincidental, or a result of stochastic processes like a series of random galactic mergers (e.g., Jahnke & Macciò 2011; Graham & Scott 2013).

In the former family of scenarios, it is believed that there is a form of "cross-talk" between the central supermassive black hole (SMBH) and its surrounding host galaxy. A robust, physical interpretation of that "cross-talk" is still widely debated. One candidate physical process that could be responsible for the interplay between the galactic nuclei and their host galaxies is generally labeled as feedback and entails a large number of potentially different mechanisms that affect and effectively regulate the growth of either the host galaxy, or the central SMBH, or both. Such regulation has been shown to be important in large cosmological simulations (e.g., Springel et al. 2005; Croton et al. 2006). The regulated growth of these two components, potentially in lockstep, has been labeled the co-evolution scenario (for a review of this, see Kormendy & Ho 2013).

One of the main ingredients of feedback models is the presence of outflows that inadvertently affect the host galaxy. The second ingredient usually common across the spectrum of feedback realizations is the presence of activity within the nucleus of the galaxy (active galactic nucleus, AGN), i.e., the ongoing accretion of matter onto the central SMBH of the galaxy (e.g., Fabian 2012). Alternatively, the aforementioned outflows, or winds, can be launched by supernova explosions (e.g., Dekel & Silk 1986; Efstathiou 2000; Scannapieco et al. 2006). Here we are interested in the role of AGNs as regulators of star formation and we are thus going to focus on AGN-related feedback mechanisms.

The study of star formation in the host galaxies of AGNs has always been of great interest due to a number of factors, among which are (a) the coincidence of the peaks of cosmic star formation and nuclear activity in terms of redshift (e.g., Hopkins & Beacom 2006, Kistler et al. 2013, and Richards et al. 2006, Aird et al. 2010, respectively), (b) the apparently necessary role of a feedback mechanism to halt the growth of halos above a given mass limit in cosmological simulations (e.g., Croton et al. 2006; Sijacki et al. 2009), and (c) the scaling relations, which imply that SMBHs and their host galaxies grow in tandem, establishing yet another link between cosmic black hole accretion and star formation. In the past, AGNs have traditionally been identified and their host galaxies studied using either optical selection (through their blue colors, e.g., Richards et al. 2002; Smith et al. 2005, or their spectral emission lines, e.g., Kauffmann et al. 2003; Richards et al. 2006) or X-ray selection (e.g., Comastri et al. 2003; Brandt & Hasinger 2005; Treister et al. 2011). To a lesser extent, active galaxies have been identified through their mid-IR excess emission (e.g., Lacy et al. 2004; Stern et al. 2005; Messias et al. 2012; Hanami et al. 2012) or their radio excess emission (e.g., Drake et al. 2003; Draper et al. 2011; Del Moro et al. 2013). Nevertheless, significant problems in studying the host galaxies of powerful AGNs arise due to the fact that emission is severely contaminated by the AGN at most wavelengths (far-IR being potentially the least affected; e.g., Hatziminaoglou et al. 2010).

In terms of outflows driven by an accreting SMBH, two main candidates exist that can affect star formation in the host galaxies of AGNs. Vigorously accreting SMBHs (also labeled as QSO-mode AGNs; e.g., Hardcastle et al. 2007) can effectively launch uncollimated flows, also called QSO winds, that originate inside the accretion disk corona or from the accretion disk itself, and can deposit a fraction of the accretion energy into the interstellar medium (ISM). Such winds have been observed in several AGNs (e.g., Rupke & Veilleux 2011, 2013) and are also widely employed in both semi-analytical and numerical simulations (e.g., Ostriker et al. 2010; Hopkins & Elvis 2010).

Alternatively, radio-loud active galaxies also provide ideal star formation quencher candidates as they are characterized by powerful, well-collimated outflows that are known to be able to deposit large quantities of mechanical energy in their surroundings (e.g., McNamara et al. 2005; Wagner & Bicknell 2011). The role of radio-loud AGNs and their jets in the evolution of galaxies, in particular with respect to star formation, has been studied intensively. Results can be divided in two broad families, one advocating for radio jets effectively inhibiting global star formation, with the other supporting a positive star formation feedback.

Results coming out of the Sloan Digital Sky Survey (SDSS) have shown that radio-selected, optically bright AGNs (the SDSS spectroscopy optical magnitude limit is 19.5 in the r-band) reside in massive galaxies and exhibit generally old stellar populations and weak signs of ongoing star formation activity (e.g., Kauffmann et al. 2003; Best et al. 2005, 2007). In a similar manner, a detailed analysis of some of the most luminous radio AGNs on the sky (from the 2 Jy and 3CRR samples) has found no strong star formation signatures in most of them, with only a handful being detected as actively star-forming galaxies by Spitzer IRS spectra (Dicken et al. 2012). Herbert et al. (2010) studied the stellar populations of a sample of intermediate redshift (z ∼ 0.5), high luminosity (>10²⁵ W Hz⁻¹ sr⁻¹) low-frequency radio AGNs, finding high-excitation, high-luminosity sources showing evidence for younger stellar populations compared to their low-excitation, low-luminosity counterparts (similar results were presented using Herschel data from Hardcastle et al. 2013). On a single source basis, individual radio AGNs have been found to exhibit suppressed star formation, or to show none at all, with evidence for their molecular gas being ionized or even blown out by radio jets (e.g., Nesvadba et al. 2010; Morganti et al. 2013). Negative feedback could be explained either through the warming up and ionization of the ISM and hence leading to less efficient star formation (e.g., Pawlik & Schaye 2009; Nesvadba et al. 2010), or through direct expulsion of the molecular gas from the galaxy, effectively removing the ingredient for stars to form (e.g., Nesvadba et al. 2006, 2011; Morganti 2010).

On the other side of the fence, there has also been a number of studies supporting jet-induced star formation, both using statistical methods and also looking at individual sources (e.g., Bicknell et al. 2000; Silk & Nusser 2010; Kalfountzou et al. 2012; Gaibler et al. 2012; Best & Heckman 2012; Zinn et al. 2013). These results can be explained by shocks driven by the radio jets in the ISM that compress the ISM and eventually lead to enhanced star formation efficiency. This mechanism has also been invoked in induced star formation from QSO winds (e.g., Ishibashi & Fabian 2012; Zubovas et al. 2013).

It is therefore apparent that although some form of feedback is needed to explain the observational results supporting co-evolution of central spheroids and their galaxies, much still remains unclear. Is there a direct, causal connection between the growth of SMBHs and their host galaxies? If so, when and how does this growth "regulation" happen? What is the net effect of radio jets in terms of star formation feedback? To that effect, we are interested in studying the star formation properties of a radio AGN sample and look for the putative link between the two. In the process, we also want to investigate whether this link is a positive or a negative one. Unlike QSO systems and their winds which usually dominate a galaxy's emission, radio AGNs offer an easier task in decoupling the emission coming from the active nucleus and that from the host galaxy.

To do this, we investigate the broadband spectral energy distributions (SEDs) of a sample of radio sources and try to decouple, by means of SED template fitting, the AGN and star formation components. In contrast to previous studies we use data from a very deep radio survey, therefore including faint radio AGNs which were missed by studies using surveys like FIRST or NVSS. As a result we are probing a broader range of radio jet powers. In addition, we employ data from the AKARI Infrared Satellite that offer an excellent spectral coverage and allow a detailed treatment of the different emission components in the infrared (namely, cold dust heated by star formation and warm/hot dust heated by the nuclear activity). In effect, and unique to our study, we can constrain both the star formation and AGN components of these radio sources simultaneously through our SED-fitting procedure. Combining the above, we are able to expand the study of AGN feedback in radio AGNs in terms of radio power, AGN luminosity, and star formation luminosity.

This paper is organized as follows: in Section 2 we describe the North Ecliptic Pole (NEP) field and the observations carried out in this field by the AKARI Infrared Satellite, in Section 3 we introduce the sample of sources we use and describe the data available, Section 4 contains the method used for the data analysis, while in Sections 5, 6, and 7 we present our results pertaining to the AGN content of the radio sources, the absolute star formation in the host galaxies of radio AGN, and finally their specific star formation and evidence of any feedback mechanism, respectively. In Section 8 we summarize our results, comparing them with other similar studies, discuss their importance, and go through possible caveats and shortcomings of our analysis. Finally in Section 9 we list our conclusions. Throughout the paper we assume the cosmological parameters H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_M = 0.27, and Ω_Λ = 0.73 (from the latest Wilkinson Microwave Anisotropy Probe release; Komatsu et al. 2011).

The AKARI space telescope (Murakami et al. 2007) was a satellite telescope launched by ISAS/JAXA in 2006 carrying two main instruments, the InfraRed Camera (IRC; Onaka et al. 2007) and the Far-Infrared Surveyor (Kawada et al. 2007). One of the AKARI legacy fields is the NEP field with its center at [R.A., decl.] = [18:00:00,+66:33:38]. Using the IRC instrument aboard the AKARI, a two-tier survey was conducted in the NEP field, with the first, wide tier covering a total area of ∼5.4 deg² (NEP-Wide, NEPW; Kim et al. 2012), while the second tier focused on a smaller area of ∼0.67 deg² (NEP-Deep, NEPD; Takagi et al. 2012). Of particular note is the spectral coverage of the IRC. With a total of nine spectral bands, ranging from 2.4 μm (N2 band) to 24 μm (L24 band), the instrument continuously covers the whole wavelength range, including the prominent wavelength gap (9–20 μm) that characterized observations with the IRAC and MIPS instruments aboard Spitzer.

Several parallel efforts at different wavelengths have provided a rich ancillary data set to complement the observations by AKARI. Of particular note and relevance to this work for the NEPW field are deep Galaxy Evolution Explorer (GALEX) observations of a part of the NEPW field (M. Malkan 2013, private communication), deep optical observations using the Canada–France–Hawaii Telescope (CFHT; Hwang et al. 2007) and SNUCAM (Im et al. 2010) on the 1.5m telescope of Maidanak observatory (Jeon et al. 2010), near-IR observations with the FLAMINGOS instrument on the 2.1 m telescope at KPNO (Y. Jeon et al., in preparation), and radio observations at 1.4 GHz with the WSRT (White et al. 2010). In 2012, the NEP field was also observed with the SPIRE (Griffin et al. 2010) instrument aboard the Herschel Space Observatory, at wavelengths between 250 and 500 μm. Finally, several spectroscopic campaigns have also taken place, with the NEPW mainly being covered by WIYN and MMT observations (Shim et al. 2013), as well as deep observations with the DEIMOS multi-fiber spectrograph on the Keck telescope, centered on the NEPD field (T. Takagi et al., in preparation). In total more than 2000 spectroscopic redshifts are available in the NEPW field. In addition, the NEPD field, deep optical observations by Subaru telescope and deep optical and near-IR observations by the CFHT telescope are also available (Oi et al. 2014). In Figure 1 the area coverage of the main AKARI surveys is shown together with some of the ancillary data sets (CFHT, GALEX, Herschel-SPIRE).

Zoom In Zoom Out Reset image size

2.1. Herschel-SPIRE Data

As we described in Section 1, we are interested in investigating the star formation properties of radio AGNs and in particular looking for possible feedback signatures in the radio AGN host galaxies. We need to define first a sample of radio sources, identify AGN-dominated systems, decompose the contribution of star formation in their broadband energy distributions, and finally investigate possible links to their nuclear properties. To do this, we employ the multi-wavelength data available in the NEP field (see Section 2). For a summary of all the data used in this paper, see Tables 1 and 2.

3.1. Radio–IR Source Cross-matching

Our core sample is taken from the radio survey of the NEP at 1.4 GHz. To define our radio source sample we utilize the IR AKARI catalogs of Kim et al. (2012) for the NEPW, and of Oi et al. 2014, for the NEPD, which we cross-match with the radio catalog at 1.4 GHz of White et al. (2010). For the cross-matching we require a 5σ detection in one of the near-IR AKARI bands, N2 or N3, for both the NEPW and the NEPD. From the original IR band-merged AKARI NEPW catalog containing 11494 sources, we build our base near-IR NEPW catalog that includes 61165 sources, after exclusion of stellar sources using the (J − N2)–(g − i) color–color diagram (Baldry et al. 2010). The NEPD on the other hand originally contains 10313 sources at a deeper 5σ flux limit of 9.6 μJy at N2 band. Following the same process as for the NEPW, we end up with 5523 sources that form the base near-IR NEPD catalog. For the purpose of our cross-matching we therefore use in total a sample of 66688 near-IR detected sources.

The radio catalog covers an area of 1.7 deg² (compared to the total of 5.4 deg² of the total NEP-Wide AKARI survey; see Figure 1) with a beam size of 17 arcsec and a practically uniform sensitivity of 21 μJy beam⁻¹. The final catalog contains 462 radio sources at a 5σ detection limit. Two main ways to match catalogs in different wavelengths have been widely used in the literature: the Poisson probability based method of Downes et al. (1986; e.g., Ivison et al. 2007; Hodge et al. 2013) and the Bayesian based likelihood ratio method, as described in Sutherland & Saunders (1992; e.g., Rodighiero et al. 2010b; Jarvis et al. 2010). The latter requires the assumption (based usually on the data at hand) of a global magnitude distribution of a given "type" of source, according to which the probability for a candidate source to be a background source and a true counterpart is calculated. The ratio of the two is used to identify true matches. The former follows the opposite direction, in that the proximity and magnitude of each individual candidate source are used to calculate the probability that the candidate source is not a background source. Hence, although both methods use the proximity and apparent magnitude of a candidate source to identify true counterparts, the Downes et al. (1986) method does so in a more straightforward way. For very large samples, where global distributions can be more robustly constrained from the data, the likelihood ratio method potentially yields better results. For the case of the matching of the WSRT source with the AKARI catalogs we use the Poisson-probability based method of Downes et al. (1986). In short, around each radio source we calculate the Poisson probability for each near-IR source to be within a circle of radius r_c. This Poisson probability is defined as

$$\begin{equation*} P^{*}=1-e^{-\pi r^{2}N_{m}}, \end{equation*}$$

where r is the distance of the candidate counterpart from the multi-wavelength source, and N_m is the surface number density within a radius r and limiting near-IR magnitude m. r_c is defined through the positional uncertainties of the AKARI IRC instrument (assumed here σ_IRC = 0.2 arcsec; e.g., Onaka et al. 2007) and that of the WSRT array (σ_WSRT = 5 arcsec; e.g., White et al. 2010). The expected number of events (i.e., near-IR sources) with P ⩽ P* can then be approximated (for a finite search radius r_c) as

$$\begin{equation*} E=P_{c}=\pi r_{c}^{2}N_{T}, \end{equation*}$$

for P* ⩾ P_c, and

$$\begin{equation*} E=P^{*}(1+\ln {P_{c}/P^{*}}), \end{equation*}$$

for P* < P_c. P_c is a critical Poisson probability, defined by the surface number density, N_T, at the limiting magnitude of the NIR sample.

Finally, the probability of a chance cross-identification of the source can be calculated as 1 − e^−E. The near-IR candidate with the lowest such probability is chosen to be the true counterpart. The above is done iteratively, first matching sources between the NEPW and the WSRT catalogs and then a further cross-matching is done between NEPD and the WSRT catalogs. Given the overlap of the catalogs, there is a number of radio sources matched with both an NEPD and an NEPW source. In these cases we keep the NEPD match, as the depth and quality of ancillary data for the NEPD is better than for the NEPW.

Following this procedure, the final IR-radio catalog contains 321 cross-matched sources. In Figure 2 the radio 1.4 GHz flux densities and near-IR AKARI N2 band AB magnitudes of the cross-matched sample are shown. In addition we show the radio flux densities of the radio sources that were not cross-matched with an AKARI source.

Zoom In Zoom Out Reset image size

3.2. Photometric Redshifts

For the purpose of our study, the estimation of the sources' redshift is necessary. As was described in Section 2, a number of spectroscopic surveys were conducted within the NEP field. Given the generally faint optical nature of radio sources, most of the radio sources do not have spectroscopic redshifts. A small subsample of 32 cross-matched radio-IR sources have spectra. Photometric redshifts have been calculated previously for the NEPD field (e.g., Negrello et al. 2009; Takagi et al. 2010; Hanami et al. 2012), using different SED fitting methods. We extend these studies to calculate photometric redshifts using the publicly available LePhare code (Arnouts et al. 1999; Ilbert et al. 2006) and the latest catalogs available for the NEPW. For the NEPD, we use the photometric redshifts from Oi et al. (2014). While the full photometric redshift catalog for the whole NEP-Wide field will be presented elsewhere, we give here a short summary of the methodology.

For the photometric redshift estimation of NEPW sources we use the near-UV GALEX band as well as the full optical bands and near-IR bands, extending out to the W2 WISE band (4.6 μm). In particular we use the following bands: NUV, u*, g, r, i, z, B, R, I, J, H, N2, W1, and W2 (refer to Table 1 for the respective central wavelengths and sensitivities). Given the overlap between the radio and CFHT observations, 149 of the cross-matched NEPW sources have CFHT data, while the remaining 65 are covered by our Maidanak observations. The limited wavelength coverage of the Maidanak optical data (three bands in the optical) leads to very high photometric redshift uncertainties and a large fraction of catastrophic outliers. We therefore focus on the sources with CFHT data, thus narrowing down the number of available sources in the NEPW.

To calculate the photometric redshifts we use the set of CFHT galaxy SED templates from Ilbert et al. (2006) and the Polletta et al. (2007) AGN templates. The Polletta et al. (2007) library also includes two composite system templates (AGN+SB, Mrk 231 and IRAS19254-7245 South). Several stellar template libraries are used to exclude stellar contamination to our sample. We use the data set of existing spectroscopic redshifts in the NEP-Wide field to calibrate possible photometric offsets using the adaptive method described in, e.g., Ilbert et al. (2006). By concentrating on the overlap between the CFHT and NEPW coverage, we focus on 143 NEPW sources with good quality optical spectra and spectroscopically classified as normal galaxies (i.e., optical AGN are excluded), to train our photometric redshift code. The largest offset is found for the NUV band at 0.26 AB mag, with an average offset for all bands of ∼0.1 AB mag. In Figure 3 we show a comparison between the NEPW spectroscopic and photometric redshifts, for sources with good quality optical spectra and classified as galaxies (open circles) and for our sub-sample of NEPW radio sources with optical spectra (filled diamonds). The resulting photometric uncertainty, as expressed by the normalized median absolute deviation (NMAD; e.g., Ilbert et al. 2006), for the NEPW sample is σ_NMAD = 0.043, with a catastrophic outlier fraction of ∼5%. For the NEPD the uncertainty is slightly lower at σ_NMAD = 0.040 but with a higher outlier fraction at ∼8.7% (for a similar plot to Figure 3 for the NEPD, see Figure 4 in Oi et al. 2014).

Zoom In Zoom Out Reset image size

In total, 237 sources out of the 321 originally matched (∼74%) with photometric and/or spectroscopic redshifts form the basic sample which we will use in the following. In Figure 4 we show the redshift distribution for our sample, as well as a number of sub-samples. Most of the sources are found at z < 2. Moreover, spectroscopic redshifts for the NEPW and NEPD reach out to a similar redshift (∼1.4 and 1.9, respectively). Hence, our photometric redshift calibration is applicable up to that limit. For the following analysis we shall use only sources at z < 2.

Zoom In Zoom Out Reset image size

In order to first identify the AGN in our samples and then decouple the nuclear emission from a possible star formation component, we want to study the broadband SEDs of our IR-radio sources. Given the richness of the ancillary data for the NEP field and the dense and homogeneous coverage of the IR wavebands from the AKARI IRC instrument, we built detailed SEDs for our sources. We utilized all the bands mentioned in Tables 1 and 2, ranging from the far-UV to the far-IR, to that end. The radio flux density was not used in the SED fitting process as templates that cover the whole range between radio and far-UV wavelengths with consistent quality are scarce. In addition, given the limited wavelength coverage of the radio regime, little physical information can be extracted by including the radio flux density in our SED fitting. In Figure 5 we present an example of such an SED.

Zoom In Zoom Out Reset image size

For the cases of GALEX and Herschel data a two-step process was followed to cross-match the individual catalogs with the IR-radio cross-matched catalog. First a proximity search was carried out using the positions of the AKARI near-IR sources (a 5 arcsec radius was used in both cases). Positions for the Herschel and GALEX sources were taken from the band-merged catalogs, based on source extraction from the 250 μm SPIRE band and the NUV band, respectively. Next, we constructed composite thumbnail images of all of our sources using the optical (CFHT or Subaru), UV (GALEX), and far-IR (Herschel-SPIRE) data. We visually inspected each thumbnail to decide whether a match in proximity was likely a true match. Through visual inspection we assigned a "Matched" and a "Confusion" flag to each source, noting both whether this appears to be a true match and whether there is a high probability for this source to be confused (i.e., several optical and/or UV sources included within the larger beam of the SPIRE instrument). These are the different cases we considered:

• 1.  
  Good match. For sources where the NIR-radio source are found at or near (<3 arcsec) the center of the SPIRE beam, the measured SPIRE fluxes are assigned to that source. Similarly, for GALEX we consider true matches cases where the NIR-radio source is at or near (<1 arcsec) the center of the instrument's beam.
• 2.  
  Obvious mismatch. For sources with obvious mismatch (sources at the edge of the instrument's beam), we assign upper limits according to the depth of the SPIRE survey of the NEP (27, 22.5, and 32.4 mJy for the 250, 350, and 500 μm bands, respectively). For GALEX we do not assign an upper limit.
• 3.  
  Heavily confused sources. For heavily confused sources we assumed sources for which more than four optical or two UV sources are enclosed within the 250 μm band beam of the SPIRE instrument. Similarly, in the case of GALEX, we consider more than two optical sources enclosed within a single GALEX beam as heavily confused. In this case, the flux of the far-IR/UV source is assigned as an upper limit.
• 4.  
  Ambiguous cases. Given the size of the SPIRE beam, there were several cases where the source cannot be easily classified in one of the above categories (light confusion and/or off-center emission). Then the mid-IR properties of the NIR-radio cross-matched source are considered. If the source is detected in at least two of the three long mid-IR bands of AKARI, then the cross-match is assumed true. In the opposite case, the flux of the initially cross-matched source is set as an upper limit. For similarly ambiguous matches with GALEX sources, the GALEX fluxes are assigned as upper limits instead of detections.

In Figure 6 we show examples of the abovementioned cases.

Zoom In Zoom Out Reset image size

In addition to detections and the Herschel-SPIRE upper limits for the confused sources, we also used AKARI and Herschel-SPIRE upper limits for all non-detected sources to further constrain the SED fitting (by definition the whole NEP field is covered by the AKARI IRC, while the Herschel-SPIRE NEP campaign has as well covered the field in its entirety). This is not the case for GALEX data, as the GALEX survey only covers part of the NEP field.

In effect, the SPIRE data provide the bottleneck of the SED fitting, defining the sub-sample of sources that have full SED coverage. In total 68 of the radio-IR sources are detected by Herschel-SPIRE. Of these, 64 have photometric and/or spectroscopic redshifts. Through the visual inspection described above, 13 sources are classified as false matches or heavily confused sources and are therefore only assigned Herschel-SPIRE upper limits. Out of the remaining 51 Herschel-SPIRE detected sources, 44 have full coverage from UV to far-IR, while the remaining 7 miss UV data. In conclusion of the 237 sources with redshifts, ∼18% have full UV-to-far-IR coverage, an additional ∼24% have full UV-to-mid-IR coverage with additional upper limits for the far-IR, ∼12% have full optical-to-mid-IR coverage and far-IR upper limits, while the remaining ∼46% is detected in a combination of optical and IR bands with upper limits for the undetected IR bands.

In Figure 7 we show the mean fractional uncertainties, η, of the different photometric bands involved in our broadband SED fitting, in an effort to estimate the error budget of our photometry. We do this separately for NEPW and NEPD sources. The highest uncertainties are seen in the AKARI mid-IR bands. However, if we consider the median values, then we see that the extreme mean values in the S11 and L24 bands are driven by a few outlier sources. As expected the next largest uncertainties are seen in the UV from GALEX and in the far-IR from Herschel. We calculate the mean fractional photometric uncertainty of our broadband SEDs to be η_avg = 6%.

Zoom In Zoom Out Reset image size

Having built our broadband SEDs following the method outlined above, we have applied a template-fitting method for the SED modeling similar to those used in the SWIRE Photometric Redshift Catalogue (Rowan-Robinson et al. 2008) or the Imperial IRAS Faint Source Catalogue (Wang & Rowan-Robinson 2009) but using the modeling code developed for Ruiz et al. (2010) and Trichas et al. (2012). This code was written in python, using the modeling and fitting software Sherpa (Freeman et al. 2001) included in CIAO 4.3 (Fruscione et al. 2006), and the Astropy package (Astropy Collaboration et al. 2013). On one hand we fitted the optical/NIR SED (wavelengths shortward of 3.2 μm) using a set of nine templates. On the other hand we fitted the IR SED using a set of four templates (to take into account the stellar contribution at short wavelengths, the best-fit optical template is included in the IR modeling). In both cases we used a χ² minimization technique. Finally, the optical and IR best-fit models were added to obtain a complete model of the SED covering the entire wavelength range.

Several physical components contribute to the emission that comprise the broadband SED of these objects, including stellar emission and star formation heated dust from the host galaxy, emission from the AGN torus, and nuclear emission from the accretion disk. In particular, the IR emission is a result of both the AGN (mainly in the near and mid-IR) and the host galaxy (far-IR being dominated by dust heated by massive young stars). As a result, the SED fitting of such objects presents a complicated task and often requires multi-component fits (e.g., Lacy et al. 2007; Seymour et al. 2008; Barthel et al. 2012). The optical to near-IR emission is mostly dominated by either a strong old stellar component, or an optical AGN, while the mid-IR to far-IR emission is mostly dominated by either a Type-2 AGN or a star formation component. Although emission from either component extends within the other half of the SED, these processes are physically distinct and can be studied quasi-separately. It is this physical motivation together with the lack of adequately good templates spanning the full wavelength range of our data that dictates the two-step SED fitting process employed here. As a result we are able to effectively separate the AGN and the galaxy component contributions for each object.

Our set of optical templates includes one elliptical galaxy, four spirals, one SB, and three QSO (see Rowan-Robinson et al. 2008 for a complete description of these templates). We included an additional extinction component to this model (using the Calzetti et al. 2000 extinction curve), with the extinction A_V as a free parameter.

We employed the set of IR templates from Rowan-Robinson et al. (2008). It includes a cirrus template (IR emission from a quiescent galaxy), an AGN hot dust torus, and two SB templates (M82 and Arp 220). Our IR model is a combination of three components: cirrus + SB + AGN. We tested this model with each SB template, and selected as the best fit the one with the lowest χ². The SED fitting takes into account the upper limits available, in addition to the detections.¹³ The use of a separate library of templates for the broadband SED fitting, compared to the fitting for the estimation of the photometric redshifts, is due to (a) the wider wavelength range of the data used in the broadband SED fitting compared to the photometric redshift estimation and (b) the need to facilitate and maximize the accuracy of the AGN and SF SED component decomposition. In particular, the Rowan-Robinson et al. (2008) templates contain only the stellar emission (or the accretion disk emission for the AGN templates) in the optical. As such we are able to decouple the fitting of the optical and IR parts of the SED, as well as the AGN and SF components. We can then combine these to create a best SED fit, while minimizing the contamination of one component on the other. This would not be possible with the Polletta et al. (2007) templates used for the photometric redshifts. In Figure 8 we show a comparison of three different sets of templates used in this study, those from Polletta et al. (2007) and Ilbert et al. (2006) for the estimation of photometric redshifts and the Rowan-Robinson et al. (2008) ones, used for the fits of the broadband SEDs. The match is very good for the wavelength range relevant for the photometric redshifts. Thus, we do not expect any systematic effects from the use of these different templates.

Zoom In Zoom Out Reset image size

In Figure 9 we show examples of template SED fittings (including the source also shown in Figure 5). After selecting the best fit (lowest χ² value) for each source, we also performed a visual check of each SED and assigned quality flags to each fit. Good fits were assigned a quality flag of 1 (the fitted SED follows all the data points closely), satisfactory fits were assigned a quality flag of 2 (one or two data points are not fitted well by the best-fit SED), while a flag value of 3 was given to bad fits (the fitted SED fails to fit more than two bands and/or large deviations are observed overall). For the following, only quality 1 (95 sources) and 2 (86 sources) fits will be considered. SED fits with flag 3 are most probably a result of either (a) wrong photometric redshift estimation or (b) a problem with the cross-matching between the different data sets. Despite the adequate range of templates employed for our fits, it is also plausible that for a subset of these sources no combination of our templates can explain their total emission leading to a bad quality fit. Finally, as was demonstrated through the process of the photometric redshift estimation, we expect that there are photometric offsets affecting each of our bands differently. Given that these can be constrained only in terms of a given template library and a training sample of spectroscopic sources, we cannot directly apply the photometric offsets calculated from LePhare to the IR-radio sources. As a result the quality of our SED fits can in some cases be further degraded.

Zoom In Zoom Out Reset image size

The following parameters are derived through the SED fitting process:

• 1.  
  optical extinction,
• 2.  
  relative contribution of each of the IR templates (cirrus, M82/Arp 220, AGN hot torus) to the 60 μm rest-frame emission,
• 3.  
  predicted flux at 0.75 (i-band), 3.6, and 60 μm rest-frame,
• 4.  
  extinction corrected optical luminosity (0.1–3 μm),
• 5.  
  IR luminosity of each IR component,
• 6.  
  total IR luminosity (L_{8–1000 μm}).

The above can be combined to derive further values such as:

• 1.  
  total AGN luminosity (L_{AGN; tot} = L_{AGN; opt} + L_{AGN; IR}),
• 2.  
  total (bolometric) source luminosity (L_{SED; tot} = L_{SED; opt} + L_{SED; IR}),
• 3.  
  total AGN fractional contribution (L_{AGN; tot}/L_{SED; tot}),
• 4.  
  IR AGN fractional contribution (L_{AGN; IR}/L_{SED; IR}),
• 5.  
  star formation luminosity (L_SF = L_M82 + L_Arp220 + L_cirrus).

Although our SED fitting code cannot explicitly calculate the uncertainty of each fit, we can estimate the uncertainties of the basic quantities derived through the SED fitting. There are at least three sources of uncertainty in our SED fits, (1) photometric uncertainty, (2) redshift uncertainty, and (3) SED fitting uncertainty. In Figures 3 and 7 we have explored the first two. The SED fitting uncertainty is less straightforward to quantify, since here we are using template SEDs, rather than physically derived model SEDs. Given sufficient points in the infrared, the total IR luminosity can be well constrained. Even with just two measurements in the mid-IR, the L_{IR; tot} can be constrained within a factor of two (e.g., Rowan-Robinson et al. 2005; Siebenmorgen & Krügel 2007). Given the dense sampling of the mid-IR spectrum from the AKARI IRC, we expect therefore the IR luminosities to be well constrained. However, we need also to consider the propagated photometric redshift error to the luminosities calculated through the SED fitting. This is found to be the dominant source of uncertainty, with a value of η_L = 0.15. Taking into account these sources of uncertainty, a quadratic sum gives a total uncertainty of $\sqrt{{{\eta _{z}^{2}+\eta _{L}^{2}+\eta _{{\rm avg}}^{2}}}}=0.167$, or 16.7%.¹⁴ We shall use this fractional uncertainty in the following to estimate typical uncertainties for derived values such as luminosities and star formation rates (SFRs).

4.1. Spectroscopic Sub-sample as a Benchmark

4.2. Importance of far-IR Data for Constraining AGN and SFRs

In the previous section we investigated the AGN recovery success rate using the spectroscopic sub-sample as a benchmark. Given the rich multi-wavelength data set at our disposal, we study how the availability of certain wavelengths affects the outcome of our SED fitting, especially in terms of the AGN component. In particular we are interested in the usefulness of far-IR data from Herschel. As we already mentioned, the contribution of the AGN component to the far-IR part of the SED is expected to be minimal. On the other hand, far-IR emission is unaffected from obscuration and is thus crucial for constraining the long-term star formation activity of a galaxy. To that end, for sources with available SPIRE data, we compare SED fits with and without the Herschel-SPIRE data points. In Figure 10 we show the comparison between the two cases for the total IR luminosity (L_{8–1000 μm};left) and the AGN component IR luminosity (right). As the SED fitting process involves separate fitting of the optical/NIR and mid-IR/far-IR parts of the broadband SED, the determination of an optical AGN is excluded from this comparison. As expected, the total IR luminosity is significantly underestimated in the absence of far-IR points, despite the presence of mid-IR measurement up to 24 μm. The difference is around half an order of magnitude, apparently independent of the total IR luminosity itself. Given the way that our SED fitting code treats upper limits, we can infer a maximum factor of two under-estimation for the far-IR luminosity (and hence SFRs) for sources with only upper limits in the SPIRE bands.

Zoom In Zoom Out Reset image size

The dusty AGN hot torus component is more significantly affected by the inclusion of far-IR measurements. The most obvious difference concerns a sizeable fraction of sources that, while in the absence of far-IR points are fitted with an IR AGN component, in actuality appear not to have an AGN IR component. These sources are shown as leftward facing arrows in the right panel of Figure 10. These sources are exclusively found at the lower AGN luminosity range. Nevertheless, those sources that do appear to be genuinely hosting an AGN component, have their luminosities accurately measured even in the absence of far-IR data points. From this we infer that the far-IR emission coming from the AGN is negligible and therefore SEDs covering up to mid-IR wavelength are sufficient to reliably constrain the AGN luminosity. In particular, especially at AGN luminosities higher than 10¹⁰L_☉, the availability of far-IR data does not seem to affect the determination of AGN properties, in terms of neither its presence nor its luminosity.

As a final remark here, it should be noted that given that we are using empirical templates to estimate AGN luminosities and that AGN, especially at the lower end of the AGN luminosity function, are poorly studied in the far-IR, we cannot draw any robust conclusion concerning the contribution of an AGN component to the far-IR. In a subsequent paper we shall study this question in more depth and assess the importance of any AGN contribution at these very long IR wavelengths.

4.3. Identifying AGN and Constraining SF

We are interested in investigating the AGN content of our radio sources and in particular in identifying those systems that show a composite behavior, their SEDs revealing signatures of both nuclear and star-forming activity. This is done purely and explicitly in terms of our SED fitting results, unlike previous studies where individual colors where employed to the same end (e.g., Takagi et al. 2007; Hanami et al. 2012). Out of the 237 sources fitted, 47 are fitted with one of the optical AGN templates, 111 are fitted with the IR AGN template, while 106 sources do not require any AGN contribution for their SED fits. There are 27 sources for which both an optical and IR AGN template is used for their SED fit. The majority of the sources has its optical emission fitted by one of the spiral galaxy templates (135), while only 33 sources are fitted with an elliptical template galaxy.

5.1. Radio Loudness

In Figure 11 we look at a measure of radio loudness, that of the ratio between rest-frame radio and near-IR fluxes (e.g., Norris et al. 2011). 3.6 μm fluxes are extracted from our best-fit SED. The rest-frame radio flux density is calculated by assuming an average radio spectral index α = 0.8 (typical of star-forming galaxies; Condon 1992). Most of our sources are consistent with star-forming galaxies (green and blue loci in Figure 11), while only a handful of the sources are found above the radio-loud QSO locus. Four sources with very high radio-to-NIR ratios can be classified as IFRGs (e.g., Norris et al. 2006). A few of the AKARI-WSRT sources are consistent with the submillimeter galaxy population (shown with filled red stars in Figure 11). At the radio flux densities of the WSRT NEP survey, it appears our sample is not dominated by radio-loud galaxies, but rather by star-forming and composite systems. The fact that we are using a uniform radio spectral index of 0.8 may be biasing the position of true radio-loud sources on these plots.

Zoom In Zoom Out Reset image size

5.2. AGN Contribution in Radio Sources

By using the SED fitting results we can for the first time see how the AGN luminosity fractional contribution behaves as a function of radio luminosity. The results are shown in Figure 12. For the highest radio luminosities there is a moderate excess of sources with a significant contribution to their bolometric luminosity by an AGN component. These are the radio AGN systems we are interested in terms of their star formation properties. We also see that the rise in the AGN contribution in terms of bolometric luminosity agrees well with previous studies for the local universe and it starts roughly at around 10⁴⁰ erg s⁻¹ (e.g., Mauch & Sadler 2007). However, we also see significant scatter, with AGNs identified at the lowest radio luminosities (presumably radio-quiet AGNs) as well as non-AGNs found in the highest radio luminosities (presumably vigorously star-forming galaxies). This drives the point that at the faintest radio flux densities the fraction of AGNs is in comparison smaller than that found in the bright radio surveys (like the FIRST and the NVSS), while the radio population is mainly dominated by star-forming (e.g., Seymour et al. 2004; Mainieri et al. 2008) or potentially composite galaxies (e.g., Strazzullo et al. 2010).

Zoom In Zoom Out Reset image size

Even at the lowest radio luminosities, where star formation should dominate, we find on average 10% contribution to the energy output from an AGN component. It should be noted however that, despite the constraints coming from the Herschel-SPIRE upper flux limits, a lack or scarcity of far-IR detections might be leading to an overestimation of the AGN component for some of these sources (as was shown in Section 4.2). Conversely, even for the highest radio luminosities for most of these sources the bolometric luminosity is actually dominated by their stellar component. For both bright and faint ends of the radio luminosity, we suffer from low-number statistics. We note in addition that the luminosity limit defined by, e.g., Mauch & Sadler (2007; shown in Figure 12) applies to the local universe. Given the evolution of both the SF and AGN luminosity functions with redshift, a shift of the dividing line toward higher luminosities at higher redshift is expected. However, since most of our sources are found at z < 1, we expect this not to be a dominant effect.

5.3. Mid-IR Colors and AGN Content

We are interested in comparing our sample with other AGNs selected in the IR. To this end, we use the near- to mid-IR colors of our sources from WISE. In Figure 13 we plot a color–color diagram utilizing the WISE detected sources in our sample. Different regions of the (W1–W2) versus (W2–W3) color–color plot have been found to be populated by different classes of objects (e.g., Wright et al. 2010; shown with differently colored contours in Figure 13). We see that our most IR-luminous sources are found within the QSO and ULIRG regions, although avoiding the ULIRG colors. This implies a significant contribution from an AGN component for these sources, resulting in a relatively flatter near-IR spectrum compared to star formation dominated systems with strong PAH features in the W1 band. Apart from these brightest sources, the rest of the sources cover the whole range of WISE colors for spiral galaxies, with a significant fraction falling within the starburst region. Our sources do not exhibit strong red W1–W2 colors, indicating the absence of heavily extincted objects. It should be noted here that the position of a source on the color–color space is not only related to its SED shape but also to its redshift.

Zoom In Zoom Out Reset image size

Assef et al. (2010) and Stern et al. (2012) have shown that AGN-dominated systems exhibit W1–W2 colors of the order of unity (flat spectrum in the near-IR). Depending on the fractional contribution of the host galaxy and the dust extinction suffered by the AGN, the W1–W2 are modulated 0.8, up to 1.8 for a range of redshifts between zero and two (Figures 1 and 2 in Stern et al. 2012). A W1–W2 cutoff of 0.8 gives an AGN sample of ∼95% reliability and higher than ∼80% completeness (Stern et al. 2012). As most of our sources lie at z < 1.4, where the distinction between the AGNs and non-active galaxies is relatively clear, the effect of the redshift range of the sample should not affect the distribution of sources on the WISE color–color plot significantly. As such, we conclude, from the mid-IR as well, that our sample lacks many purely AGN-dominated systems (W1 − W2 > 0.8).

We now turn to the star formation in the host galaxies of our IR-radio sources. To calculate the star formation rate we employ the infrared luminosity-to-SFR relation, as described in Kennicutt (1998):

$$\begin{equation*} {\rm SFR}_{{\rm IR}}[M_{\odot }\,{\rm yr^{-1}}]\,=\,4.5\,\times \,10^{-44}L_{{\rm IR}}\,[{\rm erg\,s^{-1}}], \end{equation*}$$

where L_IR is usually defined as the integrated luminosity in the wavelength range 8–1000 μm. It is known that this relation refers to starburst galaxies, for which relatively young stellar populations created continuously within the last 10–100 Myr dominate the emission between 10 and 120 μm. Although this is a theoretically derived relation, assuming a Salpeter initial mass function, empirical calibrations of the relation between SFRs and L_{8–1000 μm} exist in the literature. These are usually within 30% of the Kennicutt (1998) relation and therefore broadly agree with each other. This will be taken into consideration for the following.

Although the actual fraction is still under debate, it is known that the AGN can significantly contribute to both the mid-IR and, to lesser degree, to the far-IR luminosity of a source. It is therefore crucial to distinguish between the star formation and the AGN components in our SEDs. Therefore, in order to calculate the SFR for these sources, we first calculate the IR luminosity due to star formation. This is calculated as the sum of the infrared luminosities of the M82, Arp220, and cirrus templates used for the SED fit of each source.¹⁵ Having calculated the SFR for each object we plot these as a function of radio luminosity in Figure 14. We can observe that radio-faint sources follow the expected relation between radio-emission and SFR closely. This relation arises due to the dominating non-thermal synchrotron emission from supernovae within a star-forming galaxy, that peaks at radio wavelengths (e.g., Condon 1992). The data show a systematic offset from the Condon (1992) line, especially at low radio luminosities. If we consider the broken power-law relation of Bell (2003) this discrepancy is eliminated.

Zoom In Zoom Out Reset image size

However, as we cross the luminosity boundary above which radio AGNs should dominate, a drop in the SFR is seen with respect to the expected relation. This is due to the increasingly important contribution of AGN radio jets to the radio luminosity at 1.4 GHz. Sources that fall significantly off the expected relation should be dominated by their radio jets. It is interesting to note that although the SFR appears to taper off at the highest radio luminosities (with a decreasing tendency), it is still in absolute terms higher than the SFR of pure star-forming galaxies at lower radio luminosities (and implied lower redshifts).

We are now interested in looking at the same relation but with respect to the bolometric AGN luminosity, as it is derived from the SED fitting. This is plotted in Figure 15, where the x- and y-axis are the same as in Figure 14 but with the additional color scaling reflecting the bolometric luminosity of the AGN (here defined as L_{AGN; opt} + L_{AGN; IR}). We observe a disconnect between the most radio-luminous sources and those with the most luminous AGN component in their SEDs. Interestingly, sources with the most luminous AGN components also show the highest SFR compared to all other sources. In addition there is a trend that AGN-fitted sources show higher SFR than non-AGN ones. This becomes especially evident at radio luminosities above 10³⁹ erg s⁻¹. Focusing on the high radio luminosity end, it becomes apparent that there are sources with considerable radio excess (i.e., their emission is dominated by a radio jet) that are not picked up as AGNs from the SED fitting (these are the sources included in the dark gray diamond point at L_1.4 GHz = 10⁴¹ erg s⁻¹). These sources therefore exhibit optical and IR emission dominated by their stellar components rather than from their active nucleus. We can use the SFR–radio relation from Bell (2003) to pick out sources that show radio-excess (and should consequently harbor a radio AGN), independent of whether they are classified as AGN through the SED fitting. In Figure 15 these sources are found below the dotted line that represents the limit of a radio luminosity five times higher than expected for a given SFR. We will return to this point later on.

Zoom In Zoom Out Reset image size

We elaborate further on the apparent very high SFR for the most luminous AGNs in our sample in Figure 16, where the IR luminosity from star formation is plotted versus total AGN luminosity. We see a clear correlation between the two quantities. Sources with higher AGN luminosity also show stronger star formation. Fitting a line through all the sources we get a slope of 0.33 ± 0.05. If we fit a line only to sources with observed u − r colors typical for late-type galaxies (u − r < 2.2; e.g., Strateva et al. 2001) we get a slope of 0.38 ± 0.08. The Pearson correlation coefficient for all sources and for late-type galaxies are 0.61 and 0.63, respectively. Although not expressed in terms of positive AGN feedback, a positive relation between the two luminosities is put forward by Netzer (2009) through the study of a large sample of type 2 AGNs. This kind of relation appears to grossly underestimate the star formation of our sources especially at the lower end of AGN luminosities. A similar, slightly steeper, relation is predicted from Zubovas et al. (2013), who model the positive feedback of an AGN interacting with and inducing SF in a dense molecular gas disk within which the AGN is embedded. Finally, we also compare our data points with the positive AGN feedback model of Ishibashi & Fabian (2012), which again assumes the compression of the ISM from an AGN outflow. In that model the SFR depends on the central AGN luminosity (which drives the outflow), the SF efficiency, and the gas fraction of the galaxy. In Figure 16 we plot the solutions of this model for a range of SF efficiencies (ranging from 0.01 to 0.1, typical of normal SF galaxies and starbursts, respectively). Although at the low-AGN-luminosity end the relation appears flatter than those of both Netzer (2009) and Zubovas et al. (2013), it still falls significantly below the SFRs inferred from our data.

Zoom In Zoom Out Reset image size

The selection of this sample in the radio may be biasing our results toward higher SFRs. As we showed previously, these systems are largely dominated by their SF components and therefore are preferentially selected in the radio, over less vigorously star-forming galaxies. However, especially at the lower AGN luminosity regime the discrepancy between the models/previous observations and our results is around two orders of magnitude. This will be discussed further later on. Finally, given that this correlation is recovered for all AGNs in our sample (not only the radio-loud ones), we cannot interpret Figure 16 as evidence for radio jet positive feedback.

On the other hand, we see a good agreement with the simple theoretical model from Hickox et al. (2014; shown within the gray shaded locus for a range of redshifts from the local universe to z ∼ 2) that assumes a one-to-one positive relation between SF and AGN accretion. Indeed most of our sources, which extend to z ∼ 2, appear to be included within the two lines, with a significant fraction of sources (and namely those at intermediate radio luminosities) lying above the z ∼ 2 line of the model. It should be noted that the normalization of this model is done assuming a constant SFR to accretion rate ratio, derived for host galaxies hosting X-ray AGNs (e.g., Rafferty et al. 2011; Mullaney et al. 2012). Whether such a normalization is valid for our sample of radio-selected sources is not obvious.

It is important to differentiate between the different host galaxies of AGNs. In particular it is known that radio galaxies tend to live in galaxies more massive (e.g., Kauffmann et al. 2003; Best et al. 2005) than where other AGNs live. Moreover, there is a known dependence between SFR and stellar mass of a galaxy (also known as the main sequence of star formation), as well as the evolution of the SFR per unit mass with cosmic time (e.g., Elbaz et al. 2011). We thus need to constrain the stellar masses of our sample. To do this we use the observed SEDs to interpolate the rest-frame luminosity of each source at 2.2 μm. Old stellar populations emit the bulk of their light at NIR wavelengths and thus 2.2 μm rest-frame luminosity should be a good proxy for stellar mass. To perform the conversion between the 2.2 μm luminosity and the stellar mass we use a constant light-to-mass ratio value of 0.85 (for a compilation of references of observationally constrained mass-to-light ratios see, e.g., Portinari et al. 2004). In Figure 17 we plot the distribution of stellar masses in our sample as a function of L_1.4 GHz. For comparison we also show the stellar masses from the SDSS-FIRST sample from Best et al. (2005; these are derived from spectral features, e.g., the D4000 break, combined with stellar population synthesis models). The stellar masses of our sample are on average similar to those of the SDSS, despite the much brighter optical selection of the SDSS galaxies compared to the AKARI-WSRT sample. This can in part be attributed to the near-IR selection of our basic sample, which preferentially picks up older, more massive, galaxies. We do recover the positive correlation between L_1.4 GHz and M_stel, in that more radio-luminous sources reside in heavier host galaxies (Figure 17). It is interesting to note that the luminosity distributions of the two samples are almost the same. The AKARI-WSRT sample extends out to high redshifts while probing a relatively faint radio population. On the other hand the SDSS-FIRST sample probes a significant fraction of the entire sky in a very shallow manner, recovering very bright radio sources at relatively low redshifts.

Zoom In Zoom Out Reset image size

The specific SFR (sSFR) describes the SFR per unit of stellar mass and is a mass-normalized measure of star formation activity in galaxies. We can calculate this measure of the sources in our sample. As a result we plot in Figure 18 the sSFR of the AKARI-WSRT sources as a function of their redshift. Here the uncertainties for the redshift are taken from the σ_NMAD calculated in Section 3 while for the sSFR we take the propagated uncertainty of the SFR and the stellar mass, the former being equal to the mean SED fitting uncertainty (η_SED = 0.07) and the latter assigned as equal to that of the photometric uncertainty (η_phot = 0.06). In effect we are placing our sample on the main sequence of star-forming galaxies as defined by Elbaz et al. (2011). We see that for the local universe (z < 0.4) there is considerable scatter about the main sequence with both vigorous star-forming galaxies and "dead" galaxies present. This trend continues at higher redshifts, with however the observational effect of our flux-limited surveys leading to sources with the lowest sSFR being missed at redshifts above 0.5. Finally, there is evidence for more radio-luminous sources lying closer to and below the passive-galaxy regime than less radio-luminous ones.

Zoom In Zoom Out Reset image size

To explore the effects of our radio selection on the parameter space defined by the specific SFR and the redshift of a source, we calculate for a given redshift the minimum SFR that our radio survey can detect (with a 5σ sensitivity of around 0.1 mJy). We use the Bell (2003) relation to calculate this. As this relation is a result of a best fit to the data, scatter of the order of a factor of few is expected. In the next step we use a range of stellar masses to calculate, for a given stellar mass, the minimum specific SFR that a pure star-forming galaxy needs in order to be detected in our radio survey. In Figure 18 we plot these selection loci for two stellar mass limits that effectively contain most of our sources. We see that as expected from the stellar mass distribution of the sample (Figure 17) most of the sources fall above the locus of M_* = 10¹²M_☉. Those that fall below this line tend to have high radio luminosities, implying the presence of a jet.

In Figure 19 we again plot the sSFR as a function of redshift, this time averaging over radio luminosity bins. Let us focus on the left panel first. Here, for each radio luminosity bin we separate sources with and without an AGN component in their best-fit SEDs. Thus we want to investigate both the effect on star formation of the presence of both an AGN and a radio jet. Two observations can be directly made from this figure. For a given radio luminosity bin, sources with an AGN component appear to have higher sSFR than those with no AGN component. This trend appears reversed only in the first radio luminosity bin, which also has the lower average redshift. This is probably related to the 4 points in Figure 16 with relatively lower SFRs than the implied trend, given their AGN luminosities. However, as the mean redshifts for AGN sources are higher than non-AGN ones, an increased sSFR due to redshift evolution of the "Main Sequence" is expected. Even taking into account the uncertainties in sSFR, redshift evolution cannot account fully for the difference between the AGN and non-AGN samples.

Zoom In Zoom Out Reset image size

For the first four bins of radio luminosity (black, blue, cyan, and green stars) we observe a mild rise of sSFR with radio luminosity (on par with the expected redshift evolution and within the uncertainties of the points). This is not surprising since, at these radio luminosities (L_1.4 GHz < 10⁴⁰ erg s⁻¹), the radio-selection leads to preferential inclusion of galaxies with strong SF. The effect appears opposite for the next two bins of radio luminosity (orange and red stars). Namely, at radio luminosities higher than 10⁴⁰ erg s⁻¹ (shown in orange color in Figure 19) we see that for sources with an AGN component in their best-fit SED, increasing radio luminosity leads to a decrease of the sSFR. As can be seen when comparing the green, orange, and red stars, even though these sources have on average similar total AGN luminosities, at increasing radio luminosities their sSFR decreases.

In Figure 15 we showed that there is a sizeable fraction of sources that exhibit radio excess but are not classified as AGNs through their broadband SED fitting. We can define a sample of radio-excess sources comprising the points found below the dotted line in Figure 15. The excess radio emission of these sources implies the presence of a jet and therefore the presence of an AGN. Points above the line are classified as sources with no radio excess. A fraction of these are AGNs, whose radio emission is consistent with coming from SF in their host galaxies. However, the presence of a weak nuclear radio emission cannot be excluded. Using this new classification, we now turn to the right panel of Figure 19. Similar to the left panel, we plot the average sSFRs within radio luminosity bins but for radio-excess, all no radio-excess, and no radio-excess sources with an AGN component.

On face value, we see that when comparing between radio-excess and no radio-excess sources, the former lie exclusively at the low end of the "Main Sequence," with the latter being predominantly above it. Given our radio selection, perhaps this is not surprising. We do not see any significant trend with radio luminosity for either group of sources, beyond the expected evolution with increasing redshift. Finally, it is interesting to note that for radio-excess sources the scatter of sSFRs within each bin is very small, while the scatter for the non-excess ones appears significantly higher, as reflected by their respective uncertainties.

The comparison between the two panels of Figure 19 is not a straightforward one. A clean separation between different types of sources is not possible, as for example sources classified as AGNs in the left panel (filled circles) contain both radio-excess and non-radio-excess sources. In addition in the left panel, sources without an AGN component (from their SED fit) include sources with radio excess and hence a radio AGN. Conversely, in the right panel, radio-excess sources (filled stars) comprise both AGNs identified from the SED fitting, as well as sources that according to their SED do not have an AGN contribution. Clearly these are different flavors of AGNs and an independent third measure of activity is needed to break the ambiguity.

We showed in Figure 17 that these sources cover a wide range of stellar masses, as well as that the stellar mass correlates with the radio luminosity. Therefore, when comparing between different radio luminosity bins, we are potentially comparing between sources not only at different redshifts but also at different stellar masses. In an effort to address this issue, we narrow down our selection both in terms of mass and radio luminosity, selecting a sub-sample which spans one order of magnitude in stellar mass (10¹¹ M_☉ < M_* < 10¹² M_☉) and three orders of magnitude in radio luminosity (10^38.5 < L_1.4 GHz < 10^41.5 erg s⁻¹). We select these ranges, as within the M_*–L_1.4 GHz space, the sources are distributed relatively homogeneously, in effect alleviating any bias induced by the dependence of the stellar mass on radio luminosity. We now repeat the exercise of Figure 19 for this sub-sample, calculating the average specific SFRs for the differently selected AGN and non-AGN sources. These are shown as a function of radio luminosity in Figures 20 and 21.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 20 sources without an AGN SED component show sSFRs consistent with passive galaxies. On the contrary, for all radio luminosity bins, sources with an AGN component show higher sSFRs. In addition, we observe a drop in the sSFR when going from the 10⁴⁰ erg s⁻¹ radio luminosity bin to the 10⁴¹ erg s⁻¹ one. Moreover we note that the actual bolometric AGN luminosity (shown in numbers next to each point) is approximately the same (within errors) for all three bins. Again we observe, that for the highest radio luminosity bin, the sSFR is consistent with that of an average star-forming galaxy at that redshift.

In Figure 21 we turn to the radio-excess sources and compare them, within our narrow stellar mass bin, to the no radio-excess ones. As in Figure 19 (right panel), radio-excess sources show low sSFRs, with the first radio luminosity bin having a sSFR well below the "Main Sequence." When compared to sources without a radio excess (and hence absence of a radio jet), we see that on average they show higher sSFRs. We see therefore, that even when considering the trends between redshift, radio luminosity, and stellar mass, the behaviors seen in Figure 19 persist.

Let us summarize briefly the main findings before discussing further.

• 1.  
  Our SED fitting can constrain the presence of nuclear activity, showing that there is an excess of AGNs with increasing radio luminosity (Figure 12).
• 2.  
  Most of the sources in our sample exhibit mid-IR colors typical of late-type galaxies, with some pure AGN and early-type galaxies thrown in the mix. Given our selection criteria, we conclude that most of our sources are composite systems which are not dominated by their nuclear emission (Figure 13).
• 3.  
  We have shown that there is a strong positive correlation between star formation IR luminosity and AGN luminosity, with the most luminous AGNs in our sample exhibiting the highest SFR among the entire sample (Figures 15 and 16).
• 4.  
  On average, when looking at sources of similar bolometric AGN luminosity and radio luminosities around and above ∼10⁴⁰ erg s⁻¹, those with stronger radio emission show significantly lower specific star formation. In particular, radio-moderate AGNs appear to inhabit starburst-like galaxies at redshifts of around z ∼ 1, while at the same redshift, radio-luminous AGNs show specific star formation typical of normal star-forming galaxies at that redshift (Figure 19, left panel).
• 5.  
  When using radio-excess as a selection for radio AGNs, we find that at all radio luminosities these systems exhibit low sSFR, although still in agreement with normal star-forming galaxies at their respective redshifts (Figure 19, right panel).
• 6.  
  Even if we consider a narrower stellar mass and radio luminosity range, we recover the same trends of lower sSFR in the host galaxies of AGNs at the highest radio luminosities (Figures 20 and 21).

8.1. Comparison with Other Studies

8.2. Causality versus Coincidence

8.3. Drawbacks and Caveats

In this section we want to briefly outline potential issues that might deduct from the robustness of our results and discuss the degree they may affect our work. These can be split into three wide categories.

8.3.1. Selection Effects

Our basic selection is done in the near-IR and the radio. We start from a complete radio sample (down to a flux density of ∼100 μJy), doing an additional near-IR selection by cross-matching with the AKARI catalog. By definition radio selection picks up a combination of star-forming galaxies, starbursts, and radio AGNs. Given the high sensitivity and relatively small area coverage, we do not expect a large number of very radio-luminous AGN-dominated systems to be included in our sample. At 1.4 GHz radio emission can be produced both from star-forming galaxies and AGN radio jets. This leads to a radio-luminosity-dependent selection effect that may be affecting our analysis of quantities as a function of radio luminosity. Therefore one of the main issues that we are facing is the accurate identification and separation of star formation and nuclear activity.

We address this issue by a twofold approach, selecting AGNs either through their optical/IR SED shape or through their radio excess. These two methods are not mutually exclusive, but should instead be viewed as complimentary and in practice should select different "flavors" of AGNs. As a result, both of them individually suffer some degree of contamination (i.e., AGN sources classified as SF) and incompleteness. This is exhibited in Figure 15 where there appear to exist "non-AGN" sources, according to the SED fitting, with radio luminosities characteristic of radio AGNs. As a result, e.g., in the left panel of Figure 19 the "non-AGN" points (shown with open circles), especially at the two highest radio luminosity bins, would be contaminated by AGN sources which are missed from our SED fitting. This is addressed in the right panel of the same figure, where radio excess is used as the means of picking out radio jet-dominated systems. Nevertheless, even for those systems, a fraction of the radio luminosity might be coming from a star formation component.

A second selection bias which needs to be considered is that induced by our cross-matching with the AKARI near-IR catalogs. By definition, sources below the survey limit of the AKARI-NEPW and -NEPD in the N2 and N3 bands will not be included in these catalogs and their radio counterparts will not be cross-matched and therefore excluded from the analysis. This population of very faint near-IR sources likely represents local less massive galaxies and/or high redshift sources. The former are likely to be star-forming galaxies (since local universe radio AGNs reside in the massive end of the galaxy mass function) and their exclusion therefore should not impact our results. The latter would be either starbursting galaxies or intermediate power radio AGNs, which at intermediate radio luminosities and faint optical/NIR luminosities would be classified as radio-loud systems. Therefore, we may be missing a fraction of radio-loud (or radio-excess) sources at intermediate redshifts.

A third selection bias concerns the exclusion of radio-loud AGNs of type Fanaroff–Riley II (e.g., Urry & Padovani 1995). These radio-loud AGNs are famous for their extended radio-lobes, mainly seen at lower frequencies and at scales reaching out to Mpc. The radio power of these AGNs is mainly emitted by the radio lobes, with the central source, the core, being relatively faint or undetected. In such a case our cross-matching with the AKARI catalog would lead to a non-match. Alternatively, the radio lobes might themselves be falsely cross-matched with an AKARI source. In the first scenario, the source is completely missed by our study. In the latter case, this would introduce false radio sources in our sample, which however, given the radio spectrum of radio lobes, would not show high radio luminosity at 1.4 GHz. Additionally, these sources are prevalent at flux densities above 1 mJy (e.g., Windhorst et al. 1993), while the bulk of our sample is found at sub-mJy levels. A visual inspection of the radio map does not reveal such spurious associations. Moreover, with only 16 such sources reported by White et al. (2010) and with a cross-matched fraction of ∼70%, any effect on our results is negligible.

Finally, our radio selection may be inducing some bias toward higher star-forming galaxies, with low-SFR and dead galaxies being excluded since they fall below our radio detection limit. This was partly addressed in Figure 18 through the loci of lower limit sSFR for given stellar masses. We revisit this effect in Figure 22, where we again plot the "Main Sequence" of star formation. In addition, we use sources from our sample that show no radio-excess (above the dotted line in Figure 15) to calculate the minimum sSFR range that our sample should exhibit in order to be picked up by our radio selection. In effect, for each source in the no radio-excess sample we assign a radio flux equal to our detection limit (0.1 mJy), while retaining its redshift and stellar mass. Hence for each source we can calculate a minimum sSFR value. The range (within redshift bins) is shown with the green shaded area in Figure 22. As it becomes obvious, the "Main Sequence" is well included within the minimum range that our radio selection can probe. We therefore conclude that we are not missing a significant fraction of low sSFR sources. This is due to the exceptional depth of our radio data.

Zoom In Zoom Out Reset image size

8.3.2. Observational Effects

Although the NEP field has a wealth of multi-wavelength data, the edges of the NEPW suffer from incomplete and/or shallow optical observations. As a result we were unable to calculate reliable photometric redshifts for a significant fraction of the originally cross-matched radio sources. Although this definitely impacts the completeness and statistical robustness of this study, it does not induce any directional bias to our results. The sources missed through this observational caveat cover a wide range of radio flux densities and optical magnitudes.

A second effect pertains to the availability of far-IR data. The entire NEP field was surveyed by the Herschel-SPIRE instrument, though in a shallow manner. As a result, the bulk of the AKARI-WSRT sources remained undetected. As we showed in Figure 10 the inclusion of far-IR points is crucial both for the correct identification of AGN components and for the accurate determination of the total IR luminosity of a source. Both these quantities are crucial to this analysis. For the SED fitting of sources undetected by SPIRE we used upper limits for each SPIRE band derived according to the average sensitivity of the survey. In this way we can constrain the possible far-IR emission of all sources. Given however the statistical nature of such an upper limit, we expect a certain level of contamination from areas of the sky where data might have been problematic or where confusion did not allow proper extraction of sources. For such cases we expect that we are underestimating the upper-limits in the far-IR and as a result we are underestimating the SFR for the sources affected by this. The fraction of non-detection in the far-IR bands increases with increasing radio luminosity. Therefore we expect that sSFRs for the two higher radio luminosity bins might be somewhat underestimated.

A third effect concerns the interdependence of many of the quantities studied here with redshift. We expect a redshift and therefore a luminosity segregation of the selected sample, in that at lower redshifts (i.e., lower radio luminosities) our selection preferentially picks up more star-forming galaxies, while at intermediate and higher redshifts our selection is dominated by starburst galaxies and radio AGNs. As such, we can in part attribute the modest increase of AGN contribution with radio luminosity (Figure 12) to their common dependence on redshift. In addition, we note that the correlation between the SF and AGN luminosities shown in Figure 16 is also affected by this effect.

8.3.3. Analysis Effects

By definition, SED fitting requires a number of assumptions which by and in themselves affect our results. An often debated aspect of template SED fitting is the choice of templates. It has been shown that the use of excessively many templates leads to degeneracies and affects the results negatively. On the other hand, the exclusion of a certain "type" of template could lead to strong directional biases in the analysis. The templates used here are drawn from Rowan-Robinson et al. (2008) and have been sequentially developed and improved over several iterations (Rowan-Robinson 2003; Rowan-Robinson et al. 2004, 2005; Babbedge et al. 2004, 2006). Although based on empirical templates, they have been reproduced by means of stellar population synthesis modeling and as such have both a physical meaning and match observations well. These templates have been used to derive photometric redshifts for the SWIRE fields resulting in a redshift accuracy better than 3.5% and a very low outlier fraction (∼1%), exhibiting their high quality. Given our sample selection in the radio, the inclusion of a single elliptical, "dead," galaxy might be decreasing our ability to effectively model the SEDs of radio AGNs which predominantly should live in early-type galaxies (e.g., Best & Heckman 2012). However, the WISE colors of our sample, combined with the relatively low to intermediate radio luminosities, imply that the majority of our radio sources do not actually reside in early-type galaxies.

A second point of interest concerns the way we derive stellar masses for our sources. As was already described, we interpolated the observed SEDs to calculate luminosities at 2.2 μm and then used a constant mass-to-light ratio to transform that luminosity into stellar mass. Several points of concern might arise:

• 1.  
  at near-IR wavelengths AGN contamination is important for sources with a strong AGN component,
• 2.  
  the presence of a strong AGN component might affect the mass-to-light ratio,
• 3.  
  the mass-to-light ratio has been shown to not be constant and depending on several factors, including star formation histories and the initial mass function assumed,
• 4.  
  SED interpolation is subject to observational uncertainties and provides a rough estimate of the K-correction term needed for accurate rest-frame emission determination.

As was mentioned previously, most of our sources do not show a dominant power law in their optical/near-IR SEDs. Instead, for most cases the 1.6 μm bump can be identified, leading us to conclude that for most sources old stars dominate the near-IR part of the SED. In addition, given the lack of very powerful, QSO-like, AGNs in our sample, we do not expect the mass-to-light ratio to be affected by the AGN component. On the other hand, the light-to-mass ratio conversion is indeed highly dependent on the assumptions made and has been shown to vary among morphologically different galaxies. However, given the scope of this study, unless there is systematic difference between AGN and non-AGN and radio-faint and radio-bright sources, this would not lead to a directional bias. Although it is conceivable that there is an early-type late-type difference between bright and faint radio-galaxies, this is not supported by our data and therefore we do not expect this effect to be important. Finally, given our error budget in estimating stellar masses (light-to-mass ratio variation, AGN contamination, photometric uncertainties) we believe uncertainties due to the interpolation to not be dominant. A more physical modeling of the SEDs (using for example stellar synthesis modeling) would alleviate these issues. However, currently there is no satisfactory physical SED modeling tool to account for both star formation and AGN components in a consistent and accurate way. The fact that we recover the relation between stellar masses and radio luminosity gives us confidence about the accuracy of our method.

We have used a rich multi-wavelength data set to construct and model the broadband SEDs of a sample of radio sources in the NEP field. With the advantage of the excellent wavelength coverage of the AKARI IRC instrument and by fitting each SED with a combination of templates, both star-forming and nuclear active galaxies, we constrained these two components in each of the radio sources. We studied the relation between star formation and nuclear activity and in particular looked for evidence supporting either positive or negative feedback in radio AGNs. On face value we have found both. On closer inspection, we concluded that the positive trend between the AGN and star formation luminosity does not appear to satisfy either the assumptions or the expected relation of models of radio jet, or general AGN, induced star formation. This leads us to believe that such a correlation is either due to the coincidence of these two activities due to a third parameter which is not investigated here, due to fueling of AGNs through stellar feedback, or potentially a combination of selection and redshift effects. Visual inspection of all the optical images of our sample did not reveal, to first order, an excess of disturbed or interacting systems and therefore mergers as the third parameter which induces both SF and AGN activity is unlikely for our sources. Given the selection, luminosity ranges, and colors of our sources, we speculate that the stellar feedback scenario is plausible. However, with the current data set, we cannot exclude redshift and/or selection effects to be strongly affecting this result.

On the other hand, we find support for negative feedback that can be attributed to the presence of a radio jet. While on average radio-excess sources show low specific star formation, the sequence of optical/IR-identified AGNs with increasing radio luminosity showing progressively lower sSFRs points toward the idea of suppressed star formation in the host galaxies of radio-luminous AGNs. However, the luminosity-dependent selection function (i.e., picking progressively fewer SF galaxies at higher radio luminosities) in the radio may be biasing our results, producing in part the gradual decrease of specific SF with increasing radio luminosity. Ultimately, higher resolution and/or higher-frequency observations in the radio are needed to clearly differentiate between the different radio emission components and constrain the power of the radio jet at radio luminosities <10⁴⁰ erg s⁻¹.

Even so, it should be underlined that star formation is not quenched in these radio-luminous systems. On the contrary, in absolute terms, the SFR exhibited by these radio-loud AGNs is comparable to, if not in excess of, the lower radio luminosity (lower redshift) of non-AGN sources. As such the host galaxies of these radio AGNs are consistent with normally star-forming galaxies in their respective cosmic epochs. This is in perfect agreement with theoretical expectation of the "maintenance" nature of "radio-mode" feedback as opposed to the more energetic "QSO-mode" responsible for the putative transformation of sources into the red sequence. The finding that radio-excess-selected radio AGNs show low sSFRs (although still on the "Main Sequence" of star formation), independent of their radio luminosity, might imply that the effect of radio jets is different in different AGN activity regimes.

Taking this a step further and combining our results with the literature we can see an emerging picture of galaxy evolution where nuclear activity and star formation, at least in the circum-nuclear region, is intimately connected. Although the triggering process of AGNs is under debate, it is accepted that some gas reservoir is needed for both stars to be formed and SMBHs to be fed. As a result of this fundamental underlying connection, we see correlations between AGN and star formation luminosities. It appears that although for the brightest AGN sources mergers are predominantly important, lower power radio sources and composite systems like the ones we are studying here are potentially associated with secular processes instead. We have shown results pointing toward the suppression of star formation in the host galaxies of radio-luminous AGN systems, with SFRs however consistent with the "Main Sequence" of star formation at their respective redshifts. At higher bolometric and/or AGN luminosities (perhaps at later evolutionary stages of the sources we study here) efficient and energetic AGN winds and/or radiation pressure from an Eddington-limited accreting SMBH may totally quench star formation, transforming these galaxies into "red and dead" quiescent galaxies in the local universe.

Future observatories like SPICA and the James Webb Space Telescope will be able to cover a much wider parameter space both in terms of sensitivity and wavelength range and will therefore allow the expansion of this or similar studies to both the fainter end of the optical luminosity function as well as the less massive end of the mass function of galaxies. In the nearer term, the acquisition of spectroscopic redshifts for this sample of sources and a full optical coverage of the NEPW area would lead to a clear improvement of the results. In addition to new means of AGN classification, optical spectra will allow the estimation of the accretion rates of these radio AGNs and therefore place constraints on the characteristic properties of the "radio-mode" feedback. In a subsequent paper we will present the stacking analysis of the Herschel-SPIRE and PACS data of this radio sample, placing further stringent constraints on the SFR and specific SFR of our sample.

We thank the anonymous referee whose comments have significantly improved this manuscript. M.K. thanks Mar Mezcua for interesting discussions that have improved this manuscript. This work was supported by the National Research Foundation of Korea (NRF) grant No. 2008-0060544, funded by the Korean government (MSIP). This research is based on observations with AKARI, a JAXA project with the participation of ESA. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This research has made use of NASA's Astrophysics Data System Bibliographic Services. This research has also made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.