The mean star formation rates of unobscured QSOs: searching for evidence of suppressed or enhanced star formation
We investigate the mean star formation rates (SFRs) in the host galaxies of 3000 optically selected QSOs from the SDSS survey within the Herschel–ATLAS fields, and a radio-luminous sub-sample, covering the redshift range of 0.2–2.5. Using WISE & Herschel photometry (12 – 500) we construct composite SEDs in bins of redshift and AGN luminosity. We perform SED fitting to measure the mean infrared luminosity due to star formation, removing the contamination from AGN emission. We find that the mean SFRs show a weak positive trend with increasing AGN luminosity. However, we demonstrate that the observed trend could be due to an increase in black hole (BH) mass (and a consequent increase of inferred stellar mass) with increasing AGN luminosity. We compare to a sample of X-ray selected AGN and find that the two populations have consistent mean SFRs when matched in AGN luminosity and redshift. On the basis of the available virial BH masses, and the evolving BH mass to stellar mass relationship, we find that the mean SFRs of our QSO sample are consistent with those of main sequence star-forming galaxies. Similarly the radio-luminous QSOs have mean SFRs that are consistent with both the overall QSO sample and with star-forming galaxies on the main sequence. In conclusion, on average QSOs reside on the main sequence of star-forming galaxies, and the observed positive trend between the mean SFRs and AGN luminosity can be attributed to BH mass and redshift dependencies.
keywords:(galaxies:) quasars: supermassive black holes – galaxies: star formation – galaxies: active – galaxies: evolution
The co-evolution of a galaxy and its central supermassive black hole (BH) is a case argued by both empirical observations (e.g., the correlation of the mass of the BH and the galaxy spheroid) and results from cosmological models of galaxy evolution (see Alexander & Hickox 2012; Fabian 2012; Kormendy & Ho 2013 for reviews). This co-evolution of the galaxy and the central BH could be a result of a connection between the processes of star-formation, and BH growth. The former is commonly quantified using the star formation rate (SFR), and the latter by the luminosity of the active galactic nucleus (AGN; visible during episodes of BH growth). Since both processes are primarily fuelled by the cold gas supply within the galaxy, we may expect a first order connection between the two processes. However, models of galaxy evolution require a more interactive connection, with the AGN having a regulating role over the amount of available cold gas, and hence the SFR of the galaxy (e.g., Di Matteo, Springel & Hernquist 2005; Bower et al. 2006; Genel et al. 2014; Schaye et al. 2015).
To investigate if the AGN has indeed a regulatory role on the SFR of a galaxy there have been many studies on the star-forming properties of galaxies hosting AGN (see Harrison 2017 review). With observations from the Herschel space observatory (Herschel; Pilbratt et al. 2010) we can place strong constraints on the far-infrared emission of galaxies (FIR; ), which traces the reprocessed emission from the dusty star-forming regions (see Lutz 2014; Casey, Narayanan & Cooray 2014). Combining Herschel FIR observations with deep X-ray or optical observations, it is possible to independently constrain the AGN power in the X-ray and optical, while placing strong constraints on the SFR of the host in the FIR. However, since it is also possible for the AGN to contribute to the FIR luminosity due to the thermal re-radiation of obscuring dust from the surrounding torus (e.g. Antonucci 1993), it is important to decompose the AGN and star-formation emission at infrared wavelengths (e.g. Netzer et al. 2007; Mullaney et al. 2011; Del Moro et al. 2013; Delvecchio et al. 2014).
The majority of FIR studies of X-ray selected AGN that reach moderate to high AGN luminosities () find that the mean SFRs as a function of AGN luminosity show flat trends independently of redshift, up to 3 (e.g. Mullaney et al. 2012; Harrison et al. 2012; Rosario et al. 2012; Azadi et al. 2015; Stanley et al. 2015; Lanzuisi et al. 2017). Although this is in discrepancy with some earlier studies reporting negative trends between the mean SFRs and AGN luminosity (e.g., Page et al. 2012), an analysis by Harrison et al. (2012) demonstrated how these results are driven by small number statistics. Indeed, following studies (e.g., Azadi et al. 2015, Stanley et al. 2015, Lanzuisi et al. 2017), that used large samples of X-ray selected AGN all converge to the same results of a flat trend between the mean SFRs and AGN luminosity. In Stanley et al. (2015) we demonstrated how the flat trends can be reproduced by empirical “toy-models” that assume AGN live in star-forming galaxies (Aird et al. 2013; Hickox et al. 2014), but with AGN activity as a stochastic process, with the probability of an AGN at a given luminosity defined by the observed Eddington ratio distribution (e.g., Aird et al. 2012).
Recently hydrodynamical simulations of both isolated mergers and of full cosmological volumes have also been able to reproduce the observed flat trend between the average SFR and AGN luminosity for populations of galaxies hosting low to moderate AGN luminosities (i.e., ; e.g., Volonteri et al. 2015; McAlpine et al. 2017). In agreement with the simple “toy-models”, these simulations find that AGN luminosities can vary over several orders of magnitude for a fixed SFR (or stellar mass). However, in the simulations the underlying connection between these two processes is non universal and can be sensitive to different feeding and feedback prescriptions invoked by the simulations (e.g., Thacker et al. 2014). Crucial tests of these simulations will be to correctly reproduce the SFRs for the galaxies that host the most luminous AGN, such as Quasi-Stellar-Objects (QSOs), luminous in the optical (with ), and/or very luminous in the radio (roughly ). Such AGN have the most energetic outputs, and may be the most likely to impact directly upon the star formation of their host galaxies (e.g., Bower et al. 2017).Â
FIR studies of optically selected QSOs at are finding that they tend to live in galaxies with ongoing star formation (e.g., Kalfountzou et al. 2014; Netzer et al. 2016; Gürkan et al. 2015; Harris et al. 2016) at levels consistent with those of the star-forming population (e.g., Rosario et al. 2013). When looking at the mean SFR as a function of the bolometric AGN luminosity some studies argue for a positive correlation (e.g. Bonfield et al. 2011; Rosario et al. 2013; Kalfountzou et al. 2014; Gürkan et al. 2015; Harris et al. 2016). However, when the QSOs are selected to be FIR luminous, the mean SFR shows a flat trend with the bolometric AGN luminosity (e.g, Pitchford et al. 2016).
The most powerful AGN can sometimes also be traced by their radio emission. Powerful radio AGN can be selected in multiple ways such as a simple radio luminosity cut (e.g., McAlpine, Jarvis & Bonfield 2013; Magliocchetti et al. 2014), based on their radio loudness (i.e., ratio of radio to optical luminosity; ; Kellermann et al. 1989), which is used to split between radio-loud () and radio-quiet AGN, or based on their excitation level (or radiative efficiency), between low-excitation (radiatively inefficient) and high-excitation (radiatively efficient) radio galaxies (LERGs and HERGs respectively; Best & Heckman 2012 and references therein). FIR studies of radio AGN, with samples of HERG type AGN, find that at their hosts have ongoing star formation, independent of selection methods (e.g., Seymour et al. 2011; Karouzos et al. 2014; Kalfountzou et al. 2014; Magliocchetti et al. 2014; Drouart et al. 2014; Gürkan et al. 2015; Drouart et al. 2016; Podigachoski et al. 2016). Studies taking a luminosity cut where only the most luminous radio AGN are selected find evidence of intense FIR emission and star formation, at similar levels to the radio selected star forming galaxies, at redshifts of (e.g., Magliocchetti et al. 2014; Magliocchetti et al. 2016). Studies selecting radio-loud AGN are showing evidence of a positive trend of mean SFRs with both radio AGN luminosity (e.g., Karouzos et al. 2014), optically derived AGN bolometric luminosity (e.g., Kalfountzou et al. 2014; Gürkan et al. 2015), and AGN torus luminosity (e.g, Podigachoski et al. 2016). However, it is worth noting that LERG type AGN tend to show lower SFRs than HERG type AGN (e.g., Hardcastle et al. 2013; Gürkan et al. 2015).
A key limitation in the majority of previous studies is that they have not simultaneously taken into account the observed stellar mass and redshift dependencies of SFR observed for the global galaxy population. The average global SFR of galaxies increases with increasing redshift up to 2–3 where we observe the peak of cosmic star formation. The increase of the typical SFR with redshift has also been established for QSO samples, through studying the SFR volume density (e.g., Serjeant et al. 2010) and through the use of maximum likelihood estimators to establish a correlation (e.g., Bonfield et al. 2011). Furthermore, there is a well known stellar mass dependency of the SFR, for actively star-forming systems, which is called the main sequence of star-forming galaxies (e.g., Noeske et al. 2007; Elbaz et al. 2007; Whitaker et al. 2012; Schreiber et al. 2015). Indeed, some studies found that the BH mass (and the inferred stellar mass) is an important factor when studying the SFRs of QSOs (e.g., Rosario et al. 2013; Harris et al. 2016). These effects could be driving the observed correlations of the SFR with AGN luminosity, and need to be simultaneously taken into account when investigating such trends. An additional source of uncertainty in some studies on the SFRs of galaxies hosting AGN, is the fact that observed powerful AGN could be contributing significantly to the FIR luminosities (e.g. Drouart et al. 2014; Symeonidis et al. 2016). Not removing the potential AGN contamination to the FIR photometry used to derive SFRs can cause an artificial boost in the SFR values.
In this work, we aim to overcome the limitations outlined above. We define the mean SFRs of more than 3000 optical QSOs, selected based on their broad optical emission lines, at , and a sub-sample of 258 radio-luminous QSOs of 10, over the redshift range of 0.2 2.5. Although not selected based on the excitation level criteria, our sample consists of HERG type AGN. We compare our results to the normal star-forming galaxies of the same epoch, expanding the work of Rosario et al. (2013) to higher and lower redshifts. Furthermore, we expand the – plane of Stanley et al. (2015) to higher AGN luminosities. In our analysis we will simultaneously take into account of both redshift, and stellar mass dependencies, and remove AGN contamination from the IR luminosity. The paper is organised as follows: In section 2 we define the sample and photometry used in our work. In section 3 we present the methods followed, and in section 4 we present our initial results. Finally, in section 5 we discuss our methods and the results of our analysis, and in section 6 we present the conclusions of this work. Throughout this paper we assume , , , and a Chabrier (2003) initial mass function (IMF), unless otherwise specified.
2 Sample & Data used
The aim of this work is to constrain the mean SFRs as a function of AGN bolometric luminosity, out to very high AGN luminosities (; see Figure 1), in addition to investigating dependencies of the mean SFRs on the presence of a radio-luminous AGN.
FIR photometry provides one of the best measures of the star formation rate, as it traces the peak of the dust-reprocessed emission from star-forming regions (e.g., Kennicutt 1998; Calzetti et al. 2010; Domínguez Sánchez et al. 2014; Rosario et al. 2016). Furthermore, when studying QSO samples the optical-UV is no longer an option for determining the star formation as the QSO light dominates at these wavelengths. We use FIR data from the Herschel-ATLAS observational program (H-ATLAS; Eales et al. 2010; section 2.2) that covered the fields of GAMA09, GAMA12, and GAMA15 in its Phase 1, and the north and south galactic poles (NGP, and SGP respectively) in its Phase 2 observations. The Herschel-ATLAS fields benefit from multi-wavelength coverage (see Bourne et al. 2016 for a detailed description of all accompanying data), with excellent optical (SDSS; section 2.1), MIR and FIR photometry (WISE and Herschel; section 2.2), and radio observations (FIRST; section 2.3). We use the available data to draw a sample of optically selected QSOs from the SDSS survey with WISE and Herschel coverage (see section 2.1), determine a radio-luminous sub-sample of QSOs using the FIRST survey (see section 2.3), and measure their SFRs using WISE and Herschel observations. As we only study the fields that have overlap with the SDSS survey area, we exclude the SGP field.
2.1 Optical/SDSS QSOs
To define our QSO sample we use the publicly available SDSS data release 7 (DR7) QSO catalogue as presented in Shen et al. (2011) (see also Schneider et al. 2010 for original selection of QSOs). We chose this release as it includes the spectral analysis and virial BH mass estimates.
To provide a measurement of the power of the QSOs we use the AGN bolometric luminosity () as given in Shen et al. (2011), which has been derived from , , and , for sources at redshifts of 0.7, 0.71.9, and 1.9 respectively, using the spectral fits and bolometric corrections from the composite SED in Richards et al. (2006) (BC 9.26, BC 5.15 and BC 3.81; see Shen et al. 2011). All the QSOs of our sample have bolometric luminosities of (see Figure 1). We constrain the sample of QSOs within the regions covered by H-ATLAS.
We also make use of the virial BH mass () estimates from Shen et al. (2011), from which we estimate the stellar masses (see section 5.2.1 and Eq. 4). The have been calculated using the FWHM and continuum luminosities of the H, Mg II, and C IV lines (see section 3 of Shen et al. 2011). Specifically, the is estimated from H for sources with redshifts of 0.7, from Mg II for sources with 0.71.9, and from C IV for sources with 1.9. In Figure 2 we show the values of our sample as a function of . For comparison we also indicate three different levels of the Eddington ratio, , the ratio of over the Eddington luminosity. The of our sample covers a dynamic range of 3 orders of magnitude, with a mean and median value of 0.34 and 0.24, respectively.
Overall, this study looks at sources with redshifts of 0.2–2.5, and includes a total of 3026 QSOs, with BH masses and AGN bolometric luminosities of predominantly , and , respectively.
2.2 Mid-infrared and Far-infrared photometry
For our analysis we stack the matched-filter-smoothed PACS and SPIRE image products provided by the H-ATLAS team (see Valiante et al. 2016) for the four fields of GAMA09 (54 deg), GAMA12 (54 deg), GAMA15 (54 deg), and NGP (150 deg) that overlap with the SDSS survey. Detailed information on the construction of the images is presented in Valiante et al. (2016). The images used in our analysis have had the large scale background subtracted (i.e., the cirrus emission within our galaxy), and each pixel contains the best estimate of the flux density of a point source at that position, making them ideal for stacking analyses. In addition to the images there are also noise maps available that provide the instrumental noise at each pixel.
To define the MIR properties of our sample we use the WISE all-sky survey (Wright et al. 2010).111The WISE all-sky catalogue is available at: http://irsa.ipac.caltech.edu/Missions/wise.html Using a radius of 1 we match to the optical positions of our QSO sample described in section 2.1, with a spurious match fraction of 0.1%. 222To chose the matching radius and estimate the spurious match fraction, we follow the procedure outlined below. First we take all matches between the two catalogues that are within 20, and produce the distribution of the number of matches in bins of increasing separation. The shape of the distribution has a characteristic shape, with a peak around 0 separation, followed by a fairly steep decrease until it reaches a minimum in the number of matches. Once the separation passes the point of minimum matches, there is a steady increase in the number of matches as the separation increases. We chose the matching radius to be the separation where the minimum in the distribution occurs, and use the slope of increasing number of matches at the large separation end of the distribution to extrapolate to the smaller separations and estimate the number of spurious matches within the chosen matching radius. We find that 94.2% of our sources have a WISE counterpart. Sources in the catalogue with less than a 2 significance at a given band, have been attributed an upper limit defined by the integrated flux density measurement plus two times the measurement uncertainty. In the cases were the flux density is negative then the upper limit is defined as two times the measurement uncertainty (see the explanatory supplement to the WISE All-Sky data release, accessible through the link given in footnote 1). For our SED fitting analysis (section 3.3) we use the W3 and W4 bands at 12 and 22, respectively.
2.3 Radio data and classification
To determine the radio luminosities of our QSO sample we use the FIRST radio catalogue (Becker, White & Helfand 1995), which covers the full sky area observed by SDSS, to a sensitivity of 1 mJy. To identify the radio detected QSOs we matched the SDSS QSO catalogue to the FIRST catalogue using a 2 radius, to minimise the number of spurious matches, with a resulting spurious match fraction of 1.4% (see footnote 2). We calculate the 1.4GHz luminosity () from the catalogued flux densities, using the following equation:
where D is the luminosity distance, is the catalogued flux density, and assuming with a spectral index of 0.8. In Figure 3 we plot the radio luminosity of the radio detected sources as a function of redshift.
We classify sources as radio-luminous AGN, using a luminosity lower limit cut of 10 for 0.8, and 10 for 0.8 (see Figure 3). Based on work from McAlpine, Jarvis & Bonfield (2013), Magliocchetti et al. (2016) argue that the radio luminosity beyond which the contribution by star-forming galaxies to the total radio luminosity function becomes negligible, increases from 10 in the local Universe up to at redshift of 1.8, after which it remains constant. Our luminosity cut is always higher than these thresholds, indicating that we are selecting sources where the AGN is dominating the radio emission, and do not expect star forming galaxies to be contaminating our selection. Furthermore, in section 4.3 we demonstrate how the radio luminosities of this sample are 1–3 orders of magnitude higher than the radio luminosities predicted from the IR luminosities due to star-formation. Therefore we are selecting only AGN dominated radio sources. Within the redshift range studied here ( 0.2–2.5), there are 258 QSOs classified as radio-luminous.
For this study we measure the average SFRs of 3026 optical QSOs as a function of their bolometric luminosity and redshift. We use multi-wavelength photometry covering the MIR–FIR wavelengths (12–500) to perform SED fitting. With the sample of QSOs explored in this study we can extend the SFR – plane of Stanley et al. (2015) by an order of magnitude in AGN luminosity, with 3026 sources covering the luminosities of 10–10. Following Stanley et al. (2015), we have divided our sample in four redshift ranges, 0.2–0.5, 0.5–0.8, 0.8–1.5, and 1.5–2.5, which then are split in bins of roughly equal number of sources (80–100 sources; see Table 1). For each – bin we performed stacking analysis in the Herschel PACS and SPIRE bands to estimate the mean 100, 160, 250, 350, and 500 fluxes (section 3.1). We also calculate the mean 12 and 22 WISE fluxes (section 3.2), and mean bolometric AGN luminosities from the optical data (see section 2.2). We then used the mean fluxes of each – bin to perform composite SED fitting to decompose the IR luminosity into the AGN and star formation contributions (section 3.3). The combination of the multi-wavelength stacking and SED fitting, provides constraints on the mean IR luminosity due to star formation free from the possible AGN contamination, and the uncertainties of monochromatic estimations.
3.1 Stacking Herschel photometry
In this section we describe the methods followed to calculate the mean stacked flux density for each – bin in our analysis. For each bin we perform a weighted mean stack of the H-ATLAS PACS-100, 160, and SPIRE-250, 350, and 500 images at the optical positions of the SDSS QSOs. In all cases we regrid the images to pixels of 1, so as to have more accurate central positioning. We used the noise maps to define the weighting on the mean, by taking the inverse of the noise as the weight, to take into account that instrumental noise changes within the maps. The equation for the weighted mean of each pixel in the stacked image is:
where is the flux density of a pixel in the stacked image, is the flux density of the equivalent pixel at all images used in the stack, and is the inverted flux density at the equivalent pixel of the noise map. We note that the results do not change if we take to be the inverse variance, with a difference of 2%.
From the mean stacked image (see Figure 4) we measure the mean flux density of each – bin. For the PACS stacks, which are in units of Jy/pixel, we integrate the flux density within an aperture of 3 radius and use the recommended aperture corrections of 2.63, and 3.57, for 100 and 160, respectively (Valiante et al. 2016). For the SPIRE stacks, which are in units of Jy/beam, we take the flux density of the central pixel.
To ensure that a stacked flux density measurement is significant, and above the noise, we perform random stacks within the image. Random stacks are stacks that are calculated for a number of random positions on the map. Because each bin includes a different number of sources from each field, we perform random stacks for each bin individually, and require that the number of random positions to be taken from each field is the same as that used to produce the stack image for the sources in the bin. We perform 10000 random stacks of the maps following the same procedure as for the stack images of the sources, to create a distribution of randomly stacked values. Examples of the resulting random stack distributions for all the bands are shown in Figure 4. The resulting random stack distributions for the SPIRE bands are not centred on zero, but are positively offset by typical values of 1.3 mJy, 2 mJy, and 0.5 mJy, for the 250, 350, and 500, respectively. The offset is caused by the fact that random stacks will include positive flux density from the confused background (i.e., blending of faint sources). These are taken into account for the science stacks below. We fit a Gaussian to each random stack distribution, and from the fit we calculate the of the distribution. We use the 3 of the random-stack distribution plus the non-zero offset as our detection limit. If the stacked flux density measurement is above the defined limit then it is a detection and we use its absolute value, if it is below the limit we take an upper limit equal to the 3 value of the random stack distribution.
The offset of the random-stack distribution described above reflects a boosting in flux density from the confusion background that causes a boost in flux density of the individual bins. We therefore remove this offset from the stacked flux density in all bands for all z– bins. QSOs are well known for their clustering (e.g., White et al. 2012 and references therein), which may cause an additional boost to the stacked flux densities. In Wang et al. (2015) it was found that due to clustering of other dusty star-forming galaxies around optical QSOs there is a 8–13% contamination to the 250–500 flux density, respectively. To place an estimate on the possible contamination due to neighbouring sources, we fit the radial light profile of the stacked images using a combination of the PSF model and a constant contamination factor constrained at longer radii (see Appendix B). We find that the contamination derived using our simple method is equivalent, to the offsets found within the random stack distributions of our bins. Consequently, the contamination measured here is still only constraining the confusion background of our stacks. It is possible that there is additional contamination due to clustering that our method is not able to constrain. However, an additional contamination of 10% in the SPIRE bands will not affect our final results on the IR luminosity by more than their 1 uncertainties.
The uncertainties on the mean fluxes are estimated using the bootstrap technique. We perform 1000 re-samplings for each bin, by randomly selecting 80% of the sources in each bin, and calculate the mean flux density of each. From the resulting distribution of mean flux densities we can define the 1 uncertainties by taking the 16th and 84th percentiles (see examples in Figure 4).
3.2 Mean flux densities of the Wise counterparts
For each – bin of our sample we took the mean flux densities at 12 and 22 for the sources with a WISE counterpart. The fraction of sources with upper limits in the 12 and 22 bands ranges between the – bins, with a median of 1.3% and 32% respectively. When present the limits show a random enough distribution amongst the measured flux densities to allow us to use the non-parametric Kaplan-Meier estimator for the calculation of the mean of each bin, including both upper limits and measured flux densities (K-M method; e.g., Stanley et al. 2015; see Feigelson & Nelson 1985 for more details). We use this method for the estimation of the mean WISE fluxes in each bin of our sample. We chose this method over stacking the WISE photometry, as the source extraction that has been performed by the WISE team has taken into account of instrumental effects (Wright et al. 2010), providing good quality photometry. To test how the modest fraction of sources with WISE upper limits could affect the uncertainties on our estimations, we take two extreme cases, where all the upper limit sources are given a value of 0, and where all upper limit sources are assumed detections at that limit. We calculate the mean flux density for both approaches in all – bins, and find that the range between the two approaches is less than 0.15 mJy in the 12 band, and less than 2 mJy in the 22 band, and the mean calculated with the K-M method always lies within the range of these values. Based on this we trust that the K-M method is giving realistic results. We use bootstrap re-sampling to estimate the 1 uncertainties on the means. We note that the bootstrap uncertainties on the mean fluxes are always smaller than the range estimated for the extreme cases above.
3.3 Composite SED fitting
In Figure 5 we show how the Herschel bands cover the peak of the star-forming templates at the redshifts of interest, making them essential for the estimation of the SFRs. However, the AGN could also be contributing to the FIR fluxes of each bin, especially at higher redshifts (see Figure 5). For this reason we perform SED fitting to the WISE-12 and 22, PACS-100, 160, and SPIRE-250, 350, and 500 mean flux densities of each – bin, and decompose the AGN and star formation contributions to the IR luminosity.
We follow the methods described in Stanley et al. (2015), which we briefly outline here. We simultaneously fit an AGN template and a set of star-forming templates, and leave the normalisation of the star-forming and AGN template as free parameters in the fit. The set of star-forming templates includes the five originally defined in Mullaney et al. (2011), extended by Del Moro et al. (2013) to cover a wider wavelength range (i.e., 3–10 ; however for the purposes of our SED fitting we are only fitting within the 3–1000 wavelength region), as well as the Arp220 galaxy template from Silva et al. (1998) (see Figure 5). The AGN template used in our fitting analysis was defined in Mullaney et al. (2011) based on a sample of X-ray AGN, and is shown in Figure 5. For each – bin we perform two sets of SED fitting, one using only the six different star-forming templates, and the other using the combination of the AGN and the star-forming templates. Using the BIC parameter (Bayesian Information Criteria; Schwarz 1978) to compare the two sets of fits, we determine if a fit requires the AGN component, and find that all of our bins require the presence of the AGN counterpart in their IR SEDs. The fit with the minimum BIC value is taken to be the best-fitting result. Examples of best-fit SEDs for bins at the four different redshift ranges are given in Figure 6; the resulting best-fit SEDs for all the – bins are shown in Appendix A of this paper.
From the resulting best-fit SEDs we calculate the mean IR luminosity due to star formation of each bin, , by integrating the SF component over 8–1000. The same is also done to estimate the mean IR (8–1000) luminosity of the AGN () of each bin. To determine the uncertainty on the , and , we propagate the error on the fit, and the range of luminosities of the fits within that can be argued to be equally good fits to the best fit (e.g., Liddle 2004). From the calculated values we estimate the corresponding mean SFR values by using the Kennicutt (1998) relation corrected for a Chabrier (2003) IMF. For both the SED fitting analyses and the calculation of the IR luminosities, we use the mean redshift of the sources in each – bin.
We can see from Figure 6 that as we move towards higher redshifts the strong AGN component, present in all our fits, becomes increasingly dominant in the FIR bands. Indeed, as we show in section 4.1 the AGN can contribute up to 60% to the total IR flux density at redshifts of 2. However, as our SED fits show a strong AGN component, the results of this analysis will be dependent on the AGN template of choice.
As an initial test on the suitability of the AGN template of choice, we compare the bolometric AGN luminosity derived from our fitted AGN components to that derived from the optical. To do this we use the 6 rest-frame luminosity of the fitted AGN components of our bins, and convert to an AGN bolometric luminosity with a bolometric correction factor of 8 (following Richards et al. 2006). We find that the IR derived bolometric AGN luminosity is consistent to the optical-derived bolometric AGN luminosity within a scatter of a factor of 1.5 around the “1–1” line. Consequently, we trust that the AGN template that we use is reliable for subtracting the AGN contribution to the total IR emission, and therefore for calculating the SFRs of this sample. To further verify our approach, in section 5.1.2 we perform some additional tests using different AGN templates (from Mor & Netzer 2012, and Symeonidis et al. 2016; see dot-dashed and dashed curves in Figure 5) to evaluate the effect the choice of AGN template has on our results. We find that our choice of AGN template is fair and our results will not change significantly for different AGN templates.
To examine if different selection methods could affect the shape of our resulting composite SEDs we have split each – bin based on two different selections, and repeat the stacking and SED fitting analysis described in this section. First we have used the WISE colour classification of Stern et al. (2012) for MIR AGN, and find that the majority of sources in the bins (ranging between 49â-98% for the different bins) are selected as MIR AGN. The resulting composite SEDs of the MIR AGN selected sub-sample, as well as the resulting values are consistent to those of the overall sample within a factor of 1.2 for 30/34 of the bins, while the rest lie within a factor of 2–3.3. A second selection was based on wether or not a source is detected at 250 in the 5 point source catalogue of Valiante et al. (2016),333We have matched the optical positions of the QSOs to the 5 point source catalogue of Valiante et al. (2016) using a matching radius of 4. The matching radius was chosen based on the method described in footnote 2. thus selecting FIR luminous sources. Unsurprisingly the majority of our sample is undetected in the FIR, with only 8–19 FIR detected sources in each bin. The resulting composite SEDs of the FIR undetected sub-samples, as well as the resulting values, are consistent to those of the overall sample, supporting the idea that the few FIR luminous sources in each bin are not driving the means significantly (something also demonstrated by the bootstrap distributions of our stacks; see Figure 4). Overall, we find that our mean composite SEDs are representative of the full sample and their shape, i.e. the combination of strong AGN and SF components, is still seen when splitting the sample on the MIR or FIR properties, and are not driven by biases caused by MIR or FIR bright sources.
The method followed in this study is significantly different to the one favoured by a number of previous studies performed with H-ATLAS (e.g. Hardcastle et al. 2013; Kalfountzou et al. 2014; Gürkan et al. 2015). In those studies the authors have been using monochromatic 250 luminosities to derive SFRs, where the other FIR bands are used to derive the temperature of the modified black-body SED using FIR colours. As this method does not take into account the AGN emission at FIR wavelengths, there is a level of uncertainty on the IR luminosity due to star formation, as there is possible AGN contamination; see section 4.1. However, as discussed in sections 4.2 & 4.3 our results are in general agreement with the SFR values reported in previous work, including those mentioned above.
Here we present the main results of our study on the mean SFRs of QSOs, following the analysis presented in section 3. Initially, we compare our results of mean SFRs from our composite SEDs, to those from a monochromatic derivation at 250 (section 4.1). We then investigate the SFR properties of our full QSO sample (section 4.2), and the radio-luminous sub-sample (section 4.3).
4.1 Multi-band SED fitting versus single band derivation for the calculation of star formation luminosities
A common method of previous studies in estimating the SFRs of AGN and QSOs, is using stacking at observed frame 250 from which the IR luminosity and SFRs are then inferred (e.g., Harrison et al. 2012; Kalfountzou et al. 2014; Gürkan et al. 2015). In this section we compare our results from the multi-wavelength composite SED fitting to the single band 250 derivation, where we do not take into account the contribution from the AGN. To derive the average IR luminosities (integrated over rest-frame 8–1000) from the 250 stacked fluxes, we normalise from the 6 SF galaxy templates that we used in our SED fitting method (see Figure 5, and section 3.3), to the mean flux density at 250, and take the mean of the resulting IR luminosities of the 6 SF galaxy templates (referred to as ).
In Figure 7 we compare the results of the two methods described above: (1) the mean IR luminosity derived from the observed frame 250 photometry and (2) the multi-wavelength SED fitting and decomposition method followed in our analyses. We find that for redshifts of a single-band derivation from the 250 band is not affected significantly by the AGN, with a median offset of a factor of 1.2. At redshifts of we see a more luminosity-dependent effect, with the being affected by the AGN by an increasing factor with AGN luminosity, reaching up to a factor of 2 overestimation at the highest luminosities (10). At the highest redshifts of the is consistently overestimated by a factor of 2–2.5. Similar results on the contribution of the AGN to the total IR luminosity have also been found for higher redshift QSOs (; see Schneider et al. 2015).
4.2 The mean SFRs of optical QSOs as a function of the bolometric AGN luminosity
As described in section 3, we split our sample in bins of redshift and , for which we then estimate the mean () through multi-wavelength stacking and SED fitting that decomposes the AGN and star-forming components.
In Figure 8(a) we present our results on as a function of and redshift. We see a positive trend of the as a function of of more than an order of magnitude, something also observed in previous studies (e.g., Bonfield et al. 2011 Rosario et al. 2013; Kalfountzou et al. 2014; Karouzos et al. 2014; Gürkan et al. 2015). However, when splitting in redshift ranges, we find that the observed trend is largely due to the redshift evolution of typical SFR values. Within each redshift range we still see a slight positive trend of with , with the factor of increase ranging from 1.6–6.3 (0.2–0.8dex), with the highest redshift range of showing the largest increase with .
In Figure 8(b) we show the as a function of the mean BH mass () of each bin, and see a positive trend of the with at all redshifts. This is in agreement with the results of Harris et al. (2016) on QSOs at higher redshifts (). Since the BH masses and stellar masses of the galaxies correlate (e.g., Kormendy & Ho 2013), an increase in likely reflects an increase in stellar mass. Consequently, the observed positive trend of the with (Figure 8(a)) could also be a result of increasing stellar masses (i.e., as seen for the star-forming galaxy population; e.g., Schreiber et al. 2015) rather than AGN luminosity. We explore this further in section 5.2.1.
In our previous work (Stanley et al. 2015) we constrained the for a sample of X-ray AGN in bins of redshift and . The sample of X-ray AGN covers 3 orders of magnitude in of both moderate and high luminosity AGN (). The sample of high luminosity optical QSOs in this work is ideal to extend the – plane as defined in Stanley et al. (2015) to the highest region. It also allows us to search for systematic differences between the two populations of AGN. In Figure 9 we plot the as a function of for both the X-ray AGN and optical QSOs extending the – plane to 4 orders of magnitude. Where there is overlap in AGN luminosity between the X-ray selected AGN sample of Stanley et al. (2015) and our current sample of optical QSOs, we see a good agreement in values.444We note that there is a relative uncertainty between AGN bolometric luminosities when calculated from different photometry. To estimate the possible uncertainty between the estimates of the bolometric luminosity of our optical QSOs and the X-ray AGN sample of Stanley et al. (2015), we use 2XMM to SDSS DR7 cross-correlated catalogue from Pineau et al. (2011). We take the X-ray hard band flux density and calculate a bolometric luminosity, and compare to the bolometric luminosity from the optical measurements. We take the ratio of the two, and find that there is a median offset of 3.6 (or 0.56 dex). However, despite the uncertainty on comparing these samples, the observed trends will not be significantly affected. As this is a different sample to those we compare here, and there is no definitively correct bolometric luminosity correction, we do not apply this correction, but we do indicate it in Figure 9. At the highest redshift range of the values of 4 bins at of our QSO sample seem in disagreement to those of the X-ray AGN, although they are still consistent within the scatter of the X-ray AGN sample in that redshift range. However, in this comparison the two samples have not been matched in stellar mass, and this may drive some of the differences between the values (see section 5.2.1 for a discussion on the effect of mass on the expected values). Overall, this comparison shows that the two populations of AGN have consistent mean SFRs at fixed , and that the values as a function of of our sample complement and extend the trends observed for the X-ray AGN sample.
4.3 The mean SFRs of Radio-luminous QSOs
In parallel to our analysis of the full sample of QSOs, we also analysed a sub-sample of radio-luminous QSOs selected based on a radio luminosity () cut (see Figure 2 and section 2.3). As we show below, the radio luminosities of our sample are at least an order of magnitude above those corresponding to the of our bins, and so we are confident that these radio luminosities are dominated by emission associated with the AGN and not the star formation. For each redshift range we split the sample in bins of roughly equal numbers (15–54; see Table 2). Due to the limited number of sources we can only have two bins in each redshift range. For each bin we follow the procedures described in section 3 to estimate .
In Figure 10 we plot as a function of in each bin. We also plot the IR-radio correlation for star-forming galaxies (from Magnelli et al. 2014; Pannella et al. 2015) multiplied by factors of 50, 500, and 5000, to demonstrate how the radio luminosities of our sample are a factor of 10–5000 above those corresponding to their values. In agreement with previous results on radio selected AGN (e.g., Seymour et al. 2011; Karouzos et al. 2014; Kalfountzou et al. 2014; Magliocchetti et al. 2014; Drouart et al. 2014; Gürkan et al. 2015; Drouart et al. 2016; Podigachoski et al. 2016), we find that the radio-luminous QSOs of our sample live in galaxies with significant on-going star formation. Even though we only have two luminosity bins in each redshift range, the values as a function of are suggestive of a flat trend, further implying that the radio luminosity does not originate from the star formation in these systems and also indicating the lack of a direct relationship between the star formation emission of the galaxy and the radio-emission of the QSOs. This is also found in previous studies with different sample selections to ours, such as Seymour et al. (2011), and Drouart et al. (2016).
When comparing the radio-luminous QSOs to the overall QSO sample (dominated by radio-quiet QSOs), the values are consistent within scatter, and show similar trends with redshift. This result is also in agreement with previous work by Kalfountzou et al. (2014), comparing radio-loud to radio-quiet QSOs, at similar redshifts and .
However, in our SED fitting analyses we do not take into account of the synchrotron component that can be present for radio-luminous QSOs. Consequently, our results on the could still be contaminated by synchrotron emission due to the AGN. Assuming a conservative spectral index of , we take the 1.4GHz flux density of the sources in each bin and integrate over 8–1000 to calculate the total IR luminosity due to synchrotron emission for each source. We compare the mean for each – bin to the corresponding values and find that the lower bins of each redshift range are contaminated by %, but the higher bins of each redshift range can be contaminated by 30–100% making them highly uncertain, with the most uncertain bins having . More detailed analyses using multi-wavelength radio photometry to constrain the spectral index of the sources, is required to best constrain these results.
In this section we explore the caveats in our method and the possible implications on our results (section 5.1). Following this, we discuss the possible drivers of the weak positive trends of the mean SFR with AGN luminosity seen in our results (section 5.2).
5.1 Verification of our methods
5.1.1 SED broadening
In our SED fitting approach we assume that the observed-frame wavelengths correspond to the rest-frame wavelength of the mean redshift of a given – bin, for all of the sources within the bin. That is, we do not take into account modest k-corrections due to the different redshifts of the sources within the stack. This may result in some broadening of the average SED that we did not take into account. However, as our – bins have fairly narrow redshift ranges (see Tables 1 & 2) and there is a fairly even scatter around the mean redshift of the bins, we expect that overall there should not be significant broadening effects. To test this, we shift each of our AGN and star-forming templates to the redshift of each source in our – bins. For each – bin we then take the mean of all the redshifted SED templates to get a mean SED shape, and compare to the original SED template shifted at the mean redshift of the – bin. We find that the shape of the mean redshifted SED templates is the same to the original template when shifted to the mean redshift, apart from some smoothing of the SF template PAH features. Consequently, our results on are not affected by SED broadening effects.
5.1.2 The choice of AGN template
Since the resulting composite SEDs of our sample show such a strong AGN component, our results may be sensitive to the AGN template that we assume. For this reason, we repeat our analysis using two different AGN templates, that of Mor & Netzer (2012) and that of Symeonidis et al. (2016). The template of Mor & Netzer (2012), derived from a QSO sample with similar methods to Mullaney et al. (2011), has a steeper drop-off at longer wavelengths compared to our default template (see dot-dashed curve in Figure 5). The Symeonidis et al. (2016) template, also derived from a QSO sample, has a more gradual drop-off at longer wavelengths compared to our default template (see dashed curve in Figure 5). The varying contribution of the AGN template at longer wavelengths may affect the values estimated by our SED fitting approach.
In the first case, fitting with the Mor & Netzer (2012) AGN template we find that over all – bins the results on the do not change significantly, with a maximum increase in of a factor of , and all bins remaining consistent within the 1 to the original results. In the second case where we fit using the Symeonidis et al. (2016) AGN template, we find that up to redshifts of 1.5, the results from using the two templates are consistent within our estimated 1 errors. However, at the highest 1.5–2.5 redshift range, the results using the Symeonidis et al. (2016) template show a much larger scatter to that of our original results, and show no sign of the correlation observed in our original results. Overall, the values resulting from fitting with the Symeonidis et al. (2016) template are within a factor of 3 for 31/34 of the –, with the remaining 3/34 bins, that are in the highest redshift range showing a difference of a factor of 8. This result highlights that our results for the highest redshift range are the most sensitive to the choice of AGN template used in the SED fitting analysis. However, in the very recent studies of Lani, Netzer & Lutz (2017) and Lyu & Rieke (2017), that looked at the IR AGN SED of PG quasars following a similar approach to that of Symeonidis et al. (2016), it was found that the IR AGN SED has a steeper drop-off at long wavelengths than that argued for in Symeonidis et al. (2016). Indeed, the shape is more similar to that of the Mullaney et al. (2011) and Mor & Netzer (2012) AGN templates. Consequently, even though our highest redshift range is the most sensitive in the AGN template of choice, it is not likely that our results are as strongly affected as the use of a Symeonidis et al. (2016) type template would suggest.
5.2 Understanding the observed trends between the mean SFR and AGN properties.
5.2.1 Comparing to the main sequence of star-forming galaxies
The main sequence of star-forming galaxies is defined from the observed correlation between SFR and stellar mass, and has been found to evolve with redshift (e.g., Noeske et al. 2007; Elbaz et al. 2011; Schreiber et al. 2015). Consequently, the SFR of a normal star-forming galaxy will be dependent on its stellar mass and redshift. In this subsection we test the simple hypothesis that on average QSOs lie on the main sequence of star-forming galaxies. This follows from Stanley et al. (2015), where we showed that when taking into account of the stellar masses and redshifts of the X-ray AGN sample, their mean SFRs are consistent with the main sequence of star-forming galaxies. By comparing our results to the mean SFRs of main sequence galaxies with the same redshift and stellar masses, we can test if the QSO sample shows systematic differences to the overall star-forming population. Furthermore, we can determine if the trends we observe are simply driven by the galaxy properties.
For each – bin of our sample we use the BH masses and redshifts of the individual sources to estimate the IR luminosity of main sequence galaxies () corresponding to the properties of each. We use Eq. 9 of Schreiber et al. (2015) to calculate the :
where , , and (we note that Schreiber et al. 2015 assume a Salpeter IMF for Eq. 3 which we take into account here). The 1 scatter in the relation is +/- 0.3dex and remains out to at least a redshift of 4 (Schreiber et al. 2015). As can be seen in the above equation, to estimate the we need a measurement of the stellar masses of our sample. As our sample consists of QSOs, where the QSO emission overpowers that of the host galaxy in the optical, it is very unreliable to use SED fitting methods to the optical photometry to calculate stellar masses. However, BH masses are available for all of the QSOs in our sample from Shen et al. (2011) (see section 2.1), and can be used to infer stellar masses. To convert the BH masses to stellar masses we make use of the equation defined in Bennert et al. (2011), which includes an empirically derived redshift evolution term for redshifts of . 555We remind readers that in this paper we assume a Chabrier IMF, the same is assumed in Bennert et al. (2011). However, the equation defining the main sequence of star-forming galaxies is defined for a Salpeter IMF. For this reason we multiply the stellar masses calculated using Eq. 4 by a factor of 1.8 to correct to a Salpeter IMF (Bruzual & Charlot 2003) in order to use in Eq. 3.
To establish if our optical QSOs are consistent with being a randomly selected sample from the main sequence of star-forming galaxies, we follow a similar approach to Rosario et al. (2013). For each – bin, we perform a Monte-Carlo (MC) calculation of the corresponding to the redshifts and masses of the sources in the bin. Using Eq. 4, we define a distribution of possible stellar masses for each QSO based on their BH mass, and pick a random value from the distribution. The width of the stellar mass distribution includes both the scatter in Eq. 4 and the error on the BH mass (provided by Shen et al. 2011; see section 2.1). Based on the chosen stellar mass, and the known redshift of the source we define a log-normal distribution of values centred at the luminosity from Eq. 3, with a of 0.3dex (Schreiber et al. 2015). We pick a random value from the distribution of values for the source. We repeat this approach for all sources in a – bin and then calculate the of the bin. The above process is repeated 10,000 times for each bin, and results in a distribution of from which we can define the mean and 1 range of the possible values for a given z– bin.
In Figure 11 & 12 we plot the results for in comparison to the of the QSO sample. In Figure 11 we plot as a function of in comparison to the results for main sequence star-forming galaxies for each redshift range. With the coloured lines we show the , and the coloured shaded regions correspond to the 1 uncertainty on distribution for each bin from our MC calculation. Additionally, we take the ratio of the from our analysis over that of the main sequence (). We show the / ratio as a function of in Figure 12(a), propagating the errors of the two variables. With the line we show the expected ratio for the main sequence, while the dashed lines indicate the range covered by the scatter of the main sequence relation as defined by Schreiber et al. (2015). From these two figures we can see an apparent trend in the values of QSOs relative to those of the main sequence star-forming galaxies, as a function of redshift. At the highest redshift range of 1.5 2.5, the values are systematically below the main sequence by an average factor of 0.75 (or 1.33 if taking the inverse ratio). Moving to intermediate redshifts of 0.8 1.5 the values become consistent with those of the main sequence, while at redshifts of 0.8 the values move above those of the main sequence by a factor of 1.5. However, even though some of the means are not consistent within their errors, they are still consistent within the factor of 2 scatter of the main sequence (see Figure 12(a)).
Following the same approach as for the full QSO sample, we estimate the expected IR luminosity of main sequence star-forming galaxies () at the same redshift and stellar mass (estimated from the available ) for our radio-luminous QSO sample, and compare to their . In Figure 12(b) we show the / ratio as a function of . We find that the radio-luminous QSOs have values consistent with those of the main sequence within the factor of 2 scatter of the main sequence relation, that show the same redshift dependence as the full sample. Similar results were shown by Drouart et al. (2014) at 2.5, following a similar SED fitting approach, for a smaller sample of 70 powerful radio-galaxies from the Herschel radio galaxy evolution sample (HERGÉ). Additionally a number of studies have argued for radio-AGN/QSOs living in star-forming galaxies up to redshifts of 5 (e.g., Drouart et al. 2014, Rees et al. 2016, Magliocchetti et al. 2016) following a variety of approaches.
The relative offset between our results and the main sequence of star-forming galaxies will be dependent on the M– relation and redshift evolution assumed. We have performed the same procedure for two other cases of M– relations with and without redshift evolution. In one case we used the Kormendy & Ho (2013) relation, with no redshift evolution and assuming M. In the second case we used the Häring & Rix (2004) relation, following the redshift evolution of Merloni et al. (2010). In both cases the trends of with are similar to those seen when using the Bennert et al. (2011) relation. The / values remain within a factor of 1.25 of those estimated with the Bennert et al. (2011) relation, and remain within the factor of 2 scatter of the main sequence. Furthermore, it is also possible for QSOs in our sample to be caught in a phase of having a larger than that expected by the local M–