The origin of the Far-infrared continuum of quasars
Key Words.:galaxies: evolution, (galaxies:) quasars: supermassive black holes, (galaxies:) ISM, (galaxies:) stellar content. (ISM:) dust, extinction. Radiative transfer
Context:Understanding the history of formation of quasars represents a major challenge to theoretical models. Physical insights on the connection between the central black hole and its host galaxy can be gained by means of the quasar infrared properties.
Aims:Here we investigate the origin of the far-infrared continuum of SDSS J1148+5251, using it as a prototype for the more general class of high-luminosity high-redshift quasars.
Methods:We run the radiative transfer code TRADING to follow the transfer of radiation from the central source and from stellar sources through the dusty environment of the host galaxy. We adopt simple models for the central source, including all the radiation that can travel beyond the dusty torus. The radiation from stellar sources is modeled using the code PÉGASE. The model is based on the output of the semi-analytical merger tree code, GAMETE/QSOdust, which enables to predict the evolution of the host galaxy and of its nuclear black hole, following the star formation history and chemical evolution – including dust – in all the progenitor galaxies of SDSS J1148+5251.
Results:We find that the radiation emitted by the central source, which dominates the observed spectral energy distribution from UV/optical to near and mid infrared wavelengths, can also provide an important source of heating for the dust distributed in the host galaxy, powering at least 30% and up to 70% of the observed far infrared emission at rest-frame wavelengths m. The remaining fraction is contributed by stellar sources and can only be achieved if the host galaxy is able to sustain a star formation rate of at . This points to a co-evolution scenario where, during their hierarchical assembly, the first super-massive black holes and their host galaxies first grow at the same pace until the black hole reaches a mass of and starts growing faster than its host, reaching the bright quasar phase when the black hole and stellar mass fall within the scatter of the scaling relation observed in local galaxies. This same evolutionary scenario has been recently shown to explain the properties of a larger sample of QSOs, and imply that current dynamical mass measurements may have missed an important fraction of the host galaxy stellar mass.
Conclusions:We conclude that the far-infrared luminosity of high-z quasars is a sensitive tracer of the rapidly changing physical conditions in the host galaxy. Quasars appear far-infrared bright when the host galaxy can still sustain strong starbursts, with star formation rates and progressively dim as large quasar-driven outflows deplete the host galaxy of its gas content, damping star formation and leaving the central source as the only source of dust heating.
The first comprehensive far-infrared (FIR) investigations on optically selected low- AGNs have been made using Infrared Astronomical Satellite (IRAS) and Infrared Space Observatory (ISO) data, supplemented with SCUBA and IRAM data at longer wavelengths (see e.g. Haas et al. 2000, 2003). A variety of different Spectral Energy Distribution (SED) shapes was observed, suggesting that the mid-infrared (MIR) and FIR emission originate from dust components whose heating mechanisms may be different (Haas et al. 2000). In particular, while in all sources the MIR emission could be powered by the central Active Galactic Nuclei (AGN), the relative role of starbursts and AGN in powering the FIR emission was hard to established, although Haas et al. (2003) proposed that QSOs with the largest FIR luminosities () favour a strongly dominating AGN contribution.
Progress in the field was enabled by observations with the Spitzer Space Telescope, using MIR polycyclic aromatic hydrocarbon (PAH) emission features as tracers of star formation in QSO host galaxies (Schweitzer et al. 2006; Netzer et al. 2007; Lutz et al. 2008). Indeed, the main MIR PAH emission features at 6.2, 7.7, 8.6, and 11.3 have been extensively used as extragalactic star formation indicators because of their presence in star forming environments with a wide range of physical conditions and their absence in the proximity of AGNs. While the calibration between the PAH luminosity and the star formation rate depends on the physical properties of the interstellar medium, PAH emission is sufficiently strong to be less easily outshone by the AGN emission, unlike or the 24 m continuum. Netzer et al. (2007) applied this technique to local QSOs; comparing the average SEDs of groups of QSOs that are differing in their FIR luminosity, they were able to derive a SED for the pure AGN, after subtraction of the host star formation. Lutz et al. (2008) have extended the study to a sample of mm-detected QSOs at , showing that these sources follow the same correlation between the PAH luminosity and the rest-frame FIR emission inferred for local QSOs. For mm-bright QSOs, this provides strong evidence for intense star formation rates in their hosts, sometimes exceeding 1000 . However, the mm-faint end of this QSO population may contain pure AGNs, with no strong host star formation, and the observed weak rest frame FIR emission appears consistent with pure AGN-heated dust (Lutz et al. 2008).
The launch of the Herschel Space Observatory has enabled to directly observe the galaxy rest-frame FIR emission. Using large samples of type 1 and type 2 AGN in the Herschel Multi-tiered Extragalactic Survey (HerMES) fields observed with SPIRE at m, Hatziminaoglou et al. (2010) showed that quasars out to have the same FIR colours as star-forming galaxies and that a starburst was always needed to reproduce the FIR emission (see also Elbaz et al. 2010). By matching a sample of 24 m - selected QSOs in the Spitzer Wide-area InfraRed Extragalactic Survey to HerMES, Dai et al. (2012) collected 32 QSOs with with SPIRE data. These FIR-detected QSOs show cold dust emission with in the range , dust masses of and temperatures from 18 K to 80 K, similar to ULIRGs and (sub)mm detected QSOs. They find that the FIR properties can not be predicted from shorter wavelengths, because at rest-frame frequencies m FIR-detected QSOs have similar SEDs of FIR-undetected QSOs at comparable redshifts. Assuming that star formation is powering the FIR emission, about 40% of the sample requires ; yet, for some QSOs the derived SFRs are , which is unlikely and may imply a non-negligible contribution from the AGN. Multi-components SED fits of a much larger sample of 250 m selected galaxies of HerMES with Spitzer/InfraRed Spectrograph spectra show, however, that the presence of a starburst is necessary to account for the total FIR emission, and that the presence of an AGN has no noticable effect on the FIR properties of the host galaxy (Feltre et al. 2013).
Similar studies have been recently applied to 69 QSOs at by Leipski et al. (2014), who have presented SEDs covering the rest-frame wavelengths from m to m, including Spitzer and Herschel observations. At rest-frame wavelengths , the SEDs of the 10 sources which have been detected in at least 4 of the 5 Herschel bands (with PACS at 100 and 160 m, and with SPIRE at m) are completely dominated by a luminous cold dust component with and . The origin of the FIR component is further analyzed by stacking the SEDs of Herschel-detected (10 objects detected at ), partly Herschel-detected (14 objects detected at ) and 33 Herschel-non detected sources (33 objects). At rest-frame wavelengths , the average SEDs of the Herschel-detected and partly Herschel-detected sub-samples appear very similar, but partly detected sources show a substantial drop in their average SED at longer wavelengths. The fact that some optical and MIR luminous QSOs do show FIR emission and others do not is interpreted as an indication that star formation is the dominant driver for the additional FIR component observed in Herschel-detected QSOs.
All the above studies adopt physically motivated multi-component SED fitting techniques which - however - miss the energy balance between the various components. A notable exception is the work by Dwek, Galliano & Jones (2007), who have made a simple decomposition of the observed SED of SDSS J1148 into different components. They used a simple screen model with a Galactic extinction law to calculate the spectrum of the escaping stellar radiation, choosing the magnitude of the extinction so that the total reradiated FIR emission was equal to the total energy absorbed by the dust. In this way they estimated that the inferred FIR luminosity of could be powered by star formation but requires a continuous star formation rate of 2500 for 400 Myr with a top-heavy Initial Mass Function (IMF).
Here we attempt to make a step forward by using a detailed radiative transfer (RT) model to follow the contribution of the host star formation and the central AGN to dust heating. Our aim is to to assess to what extent the observed SEDs of QSOs in the rest-frame FIR can place constraints on current scenarios for the formation and evolution of the first Super-Massive Black Holes (SMBHs) and their host galaxies.
Information on the star formation rate, the nature of the stellar populations, the mass of gas and dust of the host galaxy are obtained using the data-constrained semi-analytical merger tree model, GAMETE/QSOdust (Valiante et al. 2011, 2012, 2014). This model allows to follow the star formation history and the mass growth of the nuclear Black Hole (BH) in all the progenitor galaxies that contribute to the hierarchical build-up of QSOs. Valiante et al. (2011) show that when the model is applied to SDSS J1148+5251 (hereafter SDSS J1148), one of the best-studied QSOs at , the observed properties of the QSO (BH, gas and dust masses, SFR, stellar mass) are reproduced provided that distinct evolutionary paths are followed. This same conclusion has been recently confirmed to hold for a larger sample of QSOs at (Valiante et al. 2014) which share a common formation scenario: in all the progenitor galaxies of the final host dark matter halo, stars form according to a standard, Salpeter-like IMF via quiescent star formation and efficient merger-driven bursts. At the same time, the central BH grows via gas accretion and mergers with other BHs. When the BH mass reaches a threshold value of , a strong energy-driven wind starts to clear-up the interstellar medium of dust and gas through a large outflow, damping the star formation rate and un-obscuring the line of sight toward the QSO. Although the estimated outflow rates can be as large as a few , in agreement with current observational constraints (Maiolino et a. 2012; Cicone et al. 2015), the final star formation rates are comparably large, , due to intense merger-driven bursts. As a result, all the QSOs host galaxies are characterized by final stellar masses in the range , a factor 3 - 30 larger than the upper limits allowed by the observed dynamical and molecular gas masses.
Although large uncertainties still affect dynamical mass measurements in these high- galaxies, one possible solution to this so-called stellar mass crisis would be to assume that the cold gas is converted into stars at a smaller efficiency than adopted by the scenario described above. In this way, the BH grows faster than the host galaxy and when it reaches a final mass of , the stellar mass is still (Valiante et al. 2014). Yet, in order to reproduce the observed chemical properties of the host galaxy, in particular the dust mass, the stars must form according to a top-heavy IMF, with a characterstic mass of , that allows to maximize the stellar metal and dust yield (Valiante et al. 2011). More importantly, the final SFRs are , significantly smaller than in the standard model.
It is clear that the two scenarios discussed above result in very different predictions on the relative role of the AGN and star formation as the powering sources of the FIR emission. Hence, any attempt to gain a deeper understanding of the origin of the FIR emission in high- QSOs may provide important insights on the co-evolution of the first SMBHs and their host galaxies.
The paper is organized as follows: in section 2 we give a short description of the main features of the semi-analytical model GAMETE/QSOdust and we present the standard scenario for SDSS J1148 at , using this QSO as a prototype for the general class of high- QSOs; in section 3 we describe the numerical set-up and the radiative transfer model; in sections 4 - 6, we present the adopted model for the spectrum emitted by the central AGN and by the stars, and the assumed spatial distribution of the dust in the host galaxy; in section 7 we illustrate and critically discuss our main results and in section 8 we draw our conclusions.
For consistency with previous works (Valiante et al. 2011, 2014), we adopt a Lambda Cold Dark Matter (CDM) cosmology with , , , and km/s/Mpc. The age of the Universe at a redshift is 900 Myr and the model and observations are scaled using the luminosity distance at this redshift, 64.5 Gpc.
2 Modeling the evolution of high-redshift quasars
In this section we give a brief presentation of the best-fit models for the
formation and evolution of the QSO J1148 obtained by means of
the semi-analytical code GAMETE/QSOdust (Valiante et al. 2011, 2012, 2014).
We focus here only on those features which are particularly relevant for
the present study and address the interested reader to Valiante et al. (2011)
for a comprehensive description of the code.
Dark matter and BH evolution
High- QSOs are likely to reside in dark matter (DM) halos with masses (Fan et al. 2004; Volonteri & Rees 2006; Valiante et al. 2012; Fanidakis et al. 2013). Here we assume that J1148 is hosted by a DM halo at whose evolution, backward in time, is simulated using a binary Monte Carlo merger tree algorithm (Valiante et al. 2011). To save computational time, we only resolve dark matter halos with a virial temperature that can cool through hydrogen Lyman- emission; in the first star forming halos we plant BH seeds with a mass of and stop populating newly virialized halos once their mass falls below the mass of a halo at the same redshift. We then follow the subsequent evolution of these BH seeds as they grow through mergers and gas accretion. In our formulation, gas accretion onto the central BH is Eddington-limited and regulated by the Bondi-Hoyle accretion rate, with a constant efficiency parameter, , set to reproduce the observed BH mass of J1148 of (Barth et al. 2003; Willott, McLure and Jarvis 2003).
BH feedback and outflow rates
The gas content within each progenitor galaxy is regulated by infall and outflow rates as well as by the rate of star formation and of returned gas at the end of stellar evolution. Gas outflows can be powered by supernova explosions and BH feedback; the relative importance of these two physical processes depends on the evolutionary phase of each progenitor galaxy. At redshifts , the dominant effect is AGN feedback, since BH-driven outflows have already significantly depleted the gas content of the host galaxy, leading to a down-turn of the star formation rate (see Figs. 3 and 4 in Valiante et al. 2012). It is important to stress that, despite the very schematic implementation of AGN feedback - described as energy-driven winds - the predicted outflow rate at is in good agreement with the massive gas outflow rate of detected by Maiolino et al. (2012) and Cicone et al. (2015) through IRAM PdBI observations of the [CII] 158 m line in the host galaxy of J1148 (Valiante et al. 2012).
Star formation law
The evolution of the stellar content of the host galaxy has been investigated adopting different prescriptions for the SFR. The star formation rate can occurr through a quiescent mode, by converting a given fraction of cold gas into stars on the galaxy dynamical time, or by means of merger-induced starbursts with enhanced efficiencies, in the so-called bursting mode (see Valiante et al. 2011). The relative importance of these two modes of star formation depends on the mass ratio of the progenitor merging galaxies, . When , in equal-mass mergers, the starbursts are maximally efficient and the star formation rate can be up to 80 times greater than in the quiescent mode (see Table 2 in Valiante et al. 2011 for the numerical values of the parameters). All the models that we have explored predict a tight co-evolution of the nuclear BH and stellar content of the host galaxies in the first 600-700 Myr. However at , when the mass of the nuclear BH has grown to and AGN feedback starts to affect the gas content of the host galaxy, the rate of star formation is very sensitive to the adopted star formation efficiency. As a result, the final star formation rates and stellar masses can increase from and – for the models which adopt the lowest star formation efficiency – to and , when the largest star formation efficiency is considered. In the former case, the final stellar mass is within the limits allowed by measurements of the dynamical mass within the central of the host galaxy obtained by means of the resolved CO emission (Walter et al. 2004). Conversely, in the latter case the final stellar mass is closer to the value that would be implied by the nuclear SMBH mass of J1148 if the local BH-stellar mass relation were assumed to hold at these very high redshifts (see Figure 3 in Valiante et al. 2011). We will return to this point in Section 8.
The enrichment of the ISM is followed taking into account the release of heavy elements and dust grains by Asymptotic Giant Branch (AGB) stars and Supernovae (SN), assuming mass and metal-dependent yields (Valiante et al. 2009). We do not adopt the instantaneous recycling approximation, hence we take into account stellar evolution on the characteristic stellar lifetimes. The mass of dust is computed following the cycling of dust in the interstellar medium, taking into account grain destruction in the hot diffuse phase of the interstellar medium by interstellar shocks and grain-growth in dense molecular clouds (de Bennassuti et al. 2014; Valiante et al. 2014). Indeed, Valiante et al. (2011) show that for the adopted dust yields, stellar sources can produce up to of dust, as a result of the balance between dust production and dust destruction by the combined effect of astration and interstellar SN shocks. This is approximately a factor ten smaller than the dust mass inferred by fitting the observed rest-frame FIR-emission.
Efficient grain growth in dense interstellar clouds must be invoked to reproduce the observed dust mass
(Michałowski et al. 2010; Valiante et al. 2011).
This mechanism requires, however, that enough gas-phase metals are available to fuel grain growth (de Bennassuti et al. 2014;
Valiante et al. 2014).
The detailed analysis done by Valiante et al. (2011) convincingly showed that if the stars are formed with a Larson IMF,
the chemical properties inferred from observations of J1148 host galaxy require that either (i) the stars are formed with a characteristic stellar mass of (standard IMF) and with high-efficiency or that (ii) the stars are formed with a characteristic stellar mass of (top-heavy IMF) but with low-efficiency, (see Table 2 and Figure 7 in Valiante et al. 2011).
As we have discussed above, the two models predict different star formation rates and UV luminosities, hence FIR emission. Using the population synthesis code PÉGASE (see section 5 and references therein) and assuming that the stars form in a 50 Myr burst with solar metallicity, we estimate that (i) a SFR of leads to a FIR luminosity in the wavelength range m of for and (ii) of for . Although the star formation history that characterize the evolution of the simulated QSO is more complex and can not be approximated by a single 50 Myr burst, these conversion factors can provide a rough estimate of the FIR luminosities contributed by the stars. Given the predicted SFRs at , we find and for the standard and top-heavy IMF models, respectively. In the same wavelength range, a fit to the observed points of J1148 (see Fig.1) leads to a luminosity of (see footnote 1), close to that obtained by star formation in the host for a standard IMF model. Hence, in what follows we will adopt the standard IMF model as input to the radiative transfer calculations and discuss in section 7 the implications of our results for the other model.
3 Modeling the Spectral Energy Distribution
The observed SED of J1148 is shown in Fig. 1. Data points come from the Sloan Digital Sky Survey (Fan et al. 2003), Subaru telescope (Iwamuro et al. 2004), Spitzer (Jiang et al. 2006; Hines et al. 2006) as reported in Table 1 of Leipski et al. (2010), to which flux densities at 70, 100, 160, 250, 350 and 500 m from Herschel were added, fully sampling the peak of the rest-frame FIR emission (Leipski et al. 2013).
In the optical, the SED shows the typical power-law (close to ) of unreddened type-1 QSO, emitted by the accretion disk around the central BH. Thus the QSO can be considered as seen perpendicularly to the dusty torus surrounding the accretion disk (i.e. it is seen face-on). This dust-free view can be used to estimate the bolometric luminosity of the accretion disk, by using a template SED that extends outside of the optical wavelengths: we found a bolometric luminosity of L for the intrinsic SED used in the radiative transfer models of Stalevski et al. (2012).
The radiation travelling along directions close to edge-on passes through the torus and heats the dust to high temperatures, producing thermal radiation in the NIR. The mean QSO SED derived by Richards et al. (2006) for type-1 QSOs observed by SDSS and Spitzer indeed shows a NIR bump, but this template is unable to predict the larger output observed in this wavelength range from high-luminosity, high-redshift QSOs such as J1148 (see Fig. 7 in Leipski et al. 2013). Beyond 20 m, the emission comes from relatively colder dust in the QSO’s host (outside of the dusty torus). As discussed in Section 1, the FIR emission is generally believed to be heated by star formation in the host; however, we can not exclude that the radiation emitted by or escaping from the dusty torus of the QSO can provide an additional, non negligible contribution. Indeed, this is what we intend to assess with a detailed radiative transfer calculation.
We model the SED using the 3-D Monte Carlo Radiative Transfer code
TRADING (Bianchi 2008). The code has been developed to study dust extinction
and emission in dusty spiral galaxies (see, e.g. Bianchi & Xilouris 2011; Holwerda et al. 2012),
but it has a general structure and can be easily adapted to other astrophysical contexts.
In particular, the code uses an adaptive grid to sample the dust distribution at
different scales in environments of different densities. Despite this ability, a
full radiative transfer model from the sub-pc scale of the dusty torus up to the kpc
scales of the host galaxy is still computationally too expensive. Thus, we are
forced to make a few simplyifing assumptions:
(i) For the central source, we adopt the SED obtained by dedicated, external, models for the transfer of the radiation from the accretion disk through the dusty torus. We will present the three different models that we have explored in Section 4.
(ii) For the stellar sources in the host galaxy, we adopt the SED computed with the PÉGASE population synthesis model using as input the star formation histories, age and metallicities of the stellar populations predicted by GAMETE/QSOdust. Similarly, GAMETE/QSOdust allows to estimate the mass of gas and dust in the host. The stellar SED and the adopted spatial distributions of star, gas and dust are presented in Sections 5 and 6.
(iii) Using TRADING, we follow the transfer of radiation from the central source (accretion disk and dusty torus) and from stellar sources through the dusty environment of the host galaxy. As customary for RT models in a dusty galaxy, we consider the transfer of radiation only for Å assuming that ionizing radiation will be intercepted entirely by the neutral hydrogen in the host.
4 The central source
Several models for the radiation emitted by the dusty torus that surrounds the
accretion disk are available in the literature. These often rely on different
methods and might result in contradictory and/or non-unique solutions when fitted
to the data (see Hönig 2012 for an overview).
In this work,
our intent is not to determine the internal properties of the torus but
rather to explore central source models which provide
a reasonable agreement with the observed SED at NIR/MIR wavelengths and study their impact on the FIR emission.
To this aim we have selected three alternative models extracted from two different libraries available in the
literature: a model from Stalevski et al. (2012) that was computed using a Monte Carlo radiative transfer code similar
to TRADING (model A), and two models from the library of Hoenig & Kishimoto (2010). The first (model B) provides an
equally good agreement to the observed NIR emission as model A, but with different physical and geometrical torus
properties. The second (model C) is similar to the SED that has been recently used to fit the MIR emission of
a sample of high-z FIR-luminous QSOs, including J1148, by Leipski et al. (2013). In addition, these two classes of
models predict very different dusty torus optical depths and enable us to explore the implications of this property to the
contribution of the central source to dust heating in the host galaxy.
The first SED is taken from the simulations of Stalevski et al. (2012) that were computed using the radiative transfer code SKIRT (Baes et al. 2011)
The second model for the central source was taken from the CAT3D models of Höenig & Kishimoto (2010), where Monte Carlo simulations of individual clouds are combined with a statistical approach to compute the torus emission
The central source templates presented so far have been selected to reproduce the NIR bump at with torus emission. However, the exact nature of this emission poses a challenge to most torus models
The derived non-ionizing luminosities of the three models for the central source, that
will be used for the radiative transfer calculations, are 7.0, 6.2
and 5.4 L for model A, B, and C, respectively.
5 The stellar spectral energy distribution
As explained in section 2, the properties of the host galaxy are inferred from the output of the semi-analytical code GAMETE/QSOdust (Valiante et al. 2011, 2012, 2014). This code allows to simulate a large number of independent hierarchical histories for the DM halo that host J1148 at . Hence, the values for the star formation rate and total stellar mass for the selected model that we have quoted in Section 2 and that we report here for convenience represent the average among 50 independent merger trees, with the errors showing the associated 1 dispersion.
At the selected standard IMF model predicts a star formation rate of and a total stellar mass of (Valiante et al. 2011, 2014). The large dispersion around the mean SFR is due to sensible variations of the merger rates and the effects of AGN feedback predicted by the different hierarchical histories.
To compute the SED of the stars in the host galaxy of J1148 we use the public
code PÉGASE v2.0
Even though all the star formation histories predict a peak in the star formation rate at (see Fig. 4 in Valiante et al. 2014), the SED is largely dominated by stars formed in the last Myr of the evolution of the host, because the luminosity of the stars formed at the peak of star formation has already faded away. As expected for an optically bright QSO, Fig.1 shows that the radiation emitted by the stars at optical and NIR wavelengths is always largely subdominant with respect to that produced by the accretion disk and dusty torus.
In the radiative transfer model, stellar radiation is assumed to be emitted isotropically. The non ionizing luminosity is , for the mean and SFRs, respectively. The adopted spatial distribution for the gas, dust and stars in the host galaxy will be presented in the next section.
6 The dust distribution in the host galaxy
The radiative transfer is computed after distributing dust within a cubic volume centered
on the central source. Two dust-free regions are included in the cube:
(i) a central sphere of radius 600 pc, since the adopted SED for the central source already takes
into account radiative transfer through the dusty torus at these scales
At the selected model predicts a molecular gas mass of and a total dust mass of . Since the dominant contribution to dust enrichment comes from grain growth in dense molecular clouds, almost all of the final dust mass is associated to molecular gas (Valiante et al. 2014). Hence, we distribute the dust mass solely in spherical, homogeneous clouds simulating dust associated to molecular gas. As in Bianchi (2008), we adopt a cloud radius of about 100 pc, which is covered by about 5 cells of the adaptive grid in the standard resolution settings. With this clumpy distribution, we minimize the extinction of radiation from the central source and maximize that from the stellar component since, as we will see below, we assume the stars to be embedded in the clouds.
Interferometric observations of J1148 reveal CO emission within 2.5 kpc from the QSO’s center, which corresponds to about of molecular hydrogen (Walter et al. 2004). From these observations, we derive a molecular gas surface density of about , which we assume to represent the surface density of the individual clouds. Such high surface density clouds are not unrealistic for QSO host galaxies at high redshifts. In fact, CO excitation analyses show that the observed transitions can be modeled with a single gas component with density and temperature , suggesting that the emission comes from a very compact and central region of the host (see Carilli & Walter 2013 and references therein).
Clouds are distributed homogeneously within a radius of 3.0 kpc, so that the average surface density within 2.5 kpc from the QSO’s is similar to what observed. To be consistent with the radiative transfer models for the central source, we use the Milky Way (MW) dust composition of Draine (2003). For simplicity, we consider the mean extinction and emission properties of the models, i.e. we adopt a single grain of mean properties, rather than a full distribution of grain sizes and materials. Also, we assume thermal equilibrium emission. We will discuss the effects of these assumptions later.
Given the average molecular gas mass predicted by the model, we end up by having clouds, each one having a molecular gas mass of . The average gas to dust ratios of the clouds is , comparable to what measured for the moderately dense interstellar medium of the Milky Way (Jenkins 2009). As a result, the optical depth of the clouds along the diameter in the V-band is very large, for Draine (2003) dust composition. Thus, some edge-on lines of sight could be severely attenuated in the optical-MIR.
Finally, stellar radiation is assumed to be emitted at the center of each cloud. Given the high optical depth of clouds, the stellar radiation is entirely absorbed within each cloud. Thus, the stellar contribution to dust emission in all the models should be considered as an upper limit.
7 Results and discussion
A qualitative picture of the results of the RT calculation is presented in Fig. 3, which shows maps of the radiation emitted at 40m in the face-on direction for model A. Due to the large optical depth of the dusty clouds in the host galaxy, the central source is able to illuminate and heat only the surface of the closest clouds (left panel); the most distant clouds appear to contribute to the FIR only when the dust is heated by the stellar radiation (right panel).
More quantitatively, the full SED of J1148 in the face-on direction is shown in Fig. 4, where each panel represents the results of the RT calculation adopting one of the central source models described in Section 4. The solid lines show the predicted SED when both the radiation from the central source and from the stars contribute to heat the dust in the host galaxy, powering the FIR emission. As indicated by the shaded regions, which account for the 1 variation on the SFR among different hierarchical histories of the host galaxy (see the right panel of Fig. 1), the radiation from stellar sources affects the SED only at m; at lower wavelengths, the emission is entirely dominated by the central source. Even though it is not possible to separate the two contributions from the FIR SED (the dust emission depending in a non-linear way on the radiation absorbed by the dust from each type of sources), we can have a clue at their relative importance by comparing the above results with the RT calculation when the dust in the host galaxy is heated only by the central source (dashed line in each panel).
For model A, the observed FIR emission is marginally consistent with dust heated by the central source only.
For model B (central panel of Fig. 4) the contribution of the central source to the FIR SED is smaller, but still significant: we find (dashed line) and (shaded region). Thus, we conclude that the central source contribute to about of the FIR emission, .
Results for model C (right panel) appear to be similar, with (dashed line) and (shaded region). Hence, we find that . However, in this case, the contribution of the central source to the observed FIR emission is not only the result of dust heating in the host galaxy, but also directly through its intrinsic, unattenuated, SED. For our particular choice among the Hönig & Kishimoto (2010) models and normalization, the direct light accounts for more than half of the central source contribution to the FIR output. However, because of the rapid decrease of the SED with wavelength, alternative choices for model C (and different normalizations) might result in widely different contributions of the central source to the FIR SED. For example, in the model chosen by Leipski et al. (2013), the FIR output of direct light from the central source is about 20% of the total, similar to our choice here; in the updated model of Leipski et al. (2014), the contribution is only 5%.
The results of the RT calculation show that the contribution of
dust heating from the central source can range from moderate to significant,
depending on the adopted models and templates for the central source.
This affects the determination of the SFR from the FIR luminosity and
our ability to constrain the connection between the BH and its host
galaxy from the IR properties of QSOs.
In what follows, we discuss some critical aspects of the radiative transfer model and
how these might eventually affect the results presented above.
Limitations and uncertainties of the central source models
Among the many different models proposed in the literature, here we have explored three different templates for the radiation emitted by the dusty torus that surrounds the nuclear BH accretion disk. From a methodological point of view, the Stalevski et al. (2012) library, out of which we have selected model A, is better suited to be used as a central source of radiation in our RT model. In fact, it includes the radiation emitted by the accretion disk, so that the predicted dust-free face-on SED can be easily normalized to match the observed UV/optical data points (see section 4). Conversely, using the CAT3D library, out of which we have selected models B and C, is more challenging. This can be seen in the left panel of Fig. 5, where we have reconstructed the full face-on view of the original CAT3D models, adding to the selected dusty torus SED (models B and C) the original template used by Höenig & Kishimoto (2010) for the accretion disk. It is evident that the original accretion disk template does not show the same spectral slope as the data. Höenig & Kishimoto (2010) claim that in their models the actual shape of the spectrum is much less important than the total luminosity. Indeed, scaling their accretion disk template to the optical/UV data in the range m (dotted line), we derive a bolometric luminosity of L, which falls within the estimates reported in section 3. However, this is not the intrinsic bolometric luminosity corresponding to the dusty torus template models B and C shown in Fig. 4. In fact, scaling the CAT3D models to observations of a QSO of known luminosity distance allows to determine the physical scale of the sublimation radius (i.e. the inner edge of the dusty torus) and thus the bolometric luminosity intrinsic to each template
More critical is the problem for the dusty torus model C (dashed line in the left panel of Fig. 5).
The bolometric luminosity derived, after scaling to match
the MIR data, is L, definitely lower than what
predicted by observations. Conversely, if the model was scaled to the
optical/UV observations (i.e. if the dashed line in Fig. 5
was made to coincide with the dotted line), the MIR prediction would be
higher than the data
The results shown in the left panel of Fig. 5 suggest that choosing independent templates for the optical/UV emission of the accretion disk and for the dusty torus, as we have done for the central source models B and C, does not ensure that the heating source in the RT model is compatible with the observations. Note that this same approach has been recently followed by Leipski et al. (2013, 2014) and was applied by Barnett et al. (2015) to the SED of the most distant QSO at .
We also note that a similar problem affects the blackbdody template used by Leipski et al. (2013) to fit the observed NIR peak of high- QSO spectra (see also footnote 4 in Section 4), which is thought to arise from the radiation emitted by hot dust at the boundary between the accretion disk and the torus: in a consistent RT model, this emission should be included as an additional heating source, and would thus alter the dusty torus SED in ways that are difficult to predict without running a dedicate RT computation.
Based on the above considerations, we can reject model C
and conclude that only dusty torus SEDs peaking in the NIR (and
thus reproducing the observed bump) are acceptable, in that
they are based on intrinsic accretion disk luminosities which can
be reconciled with the observations.
Dependence on the inclination angle i
The results discussed above do not change significantly if we allow for a moderate disalignment between the polar axis of the dusty torus and the polar axis of the host galaxy. This is illustrated by the central panel of Fig. 5, where we show the predicted SEDs for models A and B using the central source only and adopting two different inclination angles. When (see also the geometric scheme in Fig. 2), the FIR luminosity contributed by the central source in model A increases only slightly, to . This is a consequence of the relatively low optical depth of the dusty torus in this model. Indeed, the left panel of Fig. 1 shows that there is a significant fraction of the radiation emitted by the accretion disk that is able to escape the dusty torus in the edge-on direction both at MIR and optical wavelengths; these two components almost equally contribute to dust heating and their effect is not sensitive to the inclination angle of the dusty-torus. In fact, when the dusty clouds in the host galaxy will be exposed to a smaller fraction of the radiation escaping the torus in the edge-on direction, , but - at the same time - will intercept a fraction of the radiation emitted by the accretion disk in the face-on direction, that is unattenuated by the dusty torus. These two components partly compensate, leading to a negligible effect of the inclination angle on the FIR emission. Given the larger optical depth of the dusty torus, in model B the net effect of the disalignment is larger (see the central panel of Fig. 5). The unattenuated UV/optical radiation from the accretion disk now intercepts the dusty clouds distributed within a solid angle . As a result, the contributions of the central source to the FIR luminosity increase to , contributing to of .
Dependence on the dust mass and composition
In the right panel of Fig.5, we show the predicted SED assuming the same heating sources as for model A, but with a factor 10 smaller dust mass, (dashed line). In this case, both the amplitude and frequency of the thermal dust emission associated to the host galaxy shift to lower values as the clouds are more transparent and the dust temperature is larger. Hence, although dust mass estimates derived by fitting the observed FIR emission are plagued by a number of uncertainties, our RT model shows that published values have not been severely over-estimated. Similarly, the mass-averaged dust temperatures computed by the RT calculations are for models A and B, respectively; these values are very close to the dust temperature estimated from the observed data points assuming optically thin dust emission (see section 2).
Throughout the analysis, we have adopted the average Draine (2003) dust model for which, with and (Bianchi 2013). Even at FIR wavelengths, some of the features predicted by the RT calculation may be affected by the adopted dust properties. Indeed, none of the three models shown in Fig. 4 can fully account for the observations at m. The excess emission could be accounted for by dust emissivity with a flatter spectral index, as predicted by SN dust models (Bianchi & Schneider 2007) and as suggested by the fits of the FIR modified blackbody in Leipski et al. (2013). However, the same excess could be the result of dust further out of the source, or passively heated clouds (without sources in them), which we have not considered. Note also that, by analyzing the optical-near infrared spectra of 33 quasars with redshifts , Gallerani et al. (2010) find that their mean extinction properties deviate from that of the SMC, which has been shown to reproduce the dust reddening of most quasars at . This may be an indication of different dust properties at high redshift.
Finally, we have checked that increasing the resolution by a factor 4 (100 pc covered by about 20 cells of the adaptive grid) does not affect our results. Similar conclusions hold for models where we have taken into account the full distribution of grain sizes and materials (rather than a single grain with mean properties), and the effects of stochastic heating (rather than thermal equilibrium emission). All the test runs show differences with respect to our reference models that are smaller than the photometric errors at m (see the dotted line in the right panel of Fig.5).
One of the main motivations of the present study was to explore to what extent the far infrared properties of high- QSOs can help constrain current scenarios for the coevolution of the BH and its host galaxy in the first Gyr of cosmic evolution.
In FIR luminous high- QSOs, such as the QSO SDSS J1148 that we have selected as a representative case, the observed FIR emission is generally attributed to an ongoing strong starburst in the host galaxy, with inferred star formation rates that can range from up to (see Table 1 in Valiante et al. 2014). Yet, these figures can only provide a loose upper limit on the star formation rate. In fact, by means of a detailed RT calculation, we show that:
the radiation emitted by the central source, which dominates the observed SED from the UV/optical to near and mid-infrared wavelengths, provides an important source of heating for the dust distributed in the host galaxy, powering at least 30 and up to 70% of the observed far infrared emission.
the FIR SED of SDSS J1148 can only be reproduced by models where the star formation rate in the host galaxy at is as large as , as predicted by the selected standard IMF model (Valiante et al. 2011). In fact, a top-heavy IMF model, with reduced star formation efficiencies, predicts a SFR at of only , too small to supply the required heating, in addition to that of the central source, and to power the observed FIR emission.
Hence, we conclude that our analysis supports a scenario for the co-evolution of the first SMBHs and their host galaxies where, during their hierarchical assembly, the evolutionary tracks followed by the systems in the plane approach the local relation from the bottom (see Fig. 5 in Valiante et al. 2014): during the first period of evolution, at redshifts , in all the progenitor galaxies stars form according to a standard, Salpeter-like IMF via quiescent star formation and merger-driven bursts. At the same time, the nuclear BHs grow through mergers with other BHs and gas accretion; hence BHs and their host galaxies grow hand-in-hand, in a symbotic way (Volonteri 2012). At , however, the BH starts growing faster than the stellar bulge and the predicted evolution progressively steepens, falling within the observed scatter of the local relation. This scenario implies that current dynamical mass measurements may have missed an important fraction of the host galaxy stellar mass.
It is important to note that the transition between the starburst-dominated stage and the QSO-dominated stage of the evolution is regulated by the effects of BH feedback: for the QSOs that we have investigated, we find that when the nuclear BH masses have grown to , strong AGN-driven outflow rates of damp the star formation rate and un-obscure the line of sight towards the QSO. At these large outflow rates, the host galaxies would be completely depleted of their gas content in less than 20 Myr, shutting down both the star formation rate and the BH activity; as a result, we predict the active and bright QSO phase to last (Valiante et al. 2014). During this relative short phase, the FIR luminosity is expected to respond to the rapidly changing physical conditions in the host galaxy. Our analysis suggest that FIR-luminous QSOs are likely tracing a phase of the evolution where the gas content in the host galaxy can still sustain strong starbursts, with . When the conditions in the ISM are such that the star formation rate fall below , the FIR emission can only be powered by dust heating from the central source, and we expect a stronger correlation between the UV/optical, NIR/MIR and the residual FIR luminosities.
It is interesting to note that this results resonate with the recent analysis by Leipski et al. (2013, 2014), where the authors find that the average SED of QSOs with Herschel detections, including SDSS J1148, are preferentially found at the high luminosity end of their sample (for and in particular for ). Moreover, as discussed in section 1, they also find that some optical and MIR luminous QSOs do show FIR emission and others do not. This is an indication that star formation is the dominant driver for the additional FIR component observed in their Herschel-detected sources. Following a similar approach, Barnett et al. (2015) has recently analyzed the SED of QSO ULAS J1120+0641 at . This object has been detected by Herschel PACS at and m but not by Herschel SPIRE at , and m, hence it would be a member of the partly (Herschel) detected objects, according to the classification by Leipski et al. (2014). On average, these sources are optically luminous AGN with strong NIR and MIR emission, but do not show a bright FIR emission, consistent with the idea that the central source is likely contributing significantly, if not dominantly, to the FIR emission. Using local scaling relations, Barnett et al. (2015) estimate a star formation rate in the range from the [CII] line luminosity and the m continuum luminosity and conclude that, at the time of observation, the BH was growing in mass 100 times faster than the stellar bulge. This is consistent with the idea that the source was likely caught in action when the QSO-driven wind has already significantly depleted the host galaxy of its gas content, suppressing star formation.
Our results are in good agreement with another study done on the same object by Li et al. (2008). By applying a 3D radiative transfer code to a hydrodynamical simulation, these authors show that the observed SED and inferred dust properties of SDSS J1148 request vigorous star formation in the merging progenitors and a substantial contribution to dust heating and the FIR luminosity from the AGN. They also suggest that the quasar host should have already built up a large stellar population by , in agreement with our best-fit evolutionary scenario described above. Our study, which is based on a semi-analytical code, complements these previous findings allowing an exploration of a larger parameter space (in terms of stellar IMF and star formation history) and a more sophisticated treatment of the chemical enrichment of the ISM by dust and metals.
Acknowledgements.We thank the anonymous Referee for her/his insightful comments, Sebastian Hönig, Christian Leipski and Dominik Riechers for the kind clarifications about their work. RS and RV acknowledge the support and hospitality of the INAF/Osservatorio Astrofisico di Arcetri during part of the work. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 306476. S.Salvadori acknowledges the support from the Netherlands Organisation for Scientific Research (NWO), VENI grant 639.041.233
- Assuming that the dust is optically thin and that it radiates as a ’grey-body’, the mass of dust can be computed as, where is the opacity coefficient per unit dust mass, is the Planck function for a dust temperature , and is the luminosity distance to the source. Using a fit to the observed points at (see section 3) and the average dust properties of Draine (2003), which we adopt for the present study, we estimate a dust mass , a dust temperature , and a total FIR luminosity .
- Stalevski et al. (2012) models are available at: https://sites.google.com/site/skirtorus
- CAT3D models are available at: http://www.sungrazer.org/CAT3D.html.
- This problem is inherent to clumpy models (see e.g. Deo et al. 2011), as smooth models reproduce this observed feature. It also depends significantly on the assumed primary source (Feltre et al. 2012).
- For a bolometric luminosity of L, Stalevski et al. (2012) adopt an external radius of the dusty torus of 15 pc. Scaling to the bolometric luminosity of J1148, L and assuming flux conservation, the radius is and we infer a radius of 600 pc.
- As customary for observations, here we neglect that the SED might be anisotropic and obtain the luminosity integrating over the whole solid angle (i.e. we multiply the integral under the SED by ).
- See Sect. 4.1 of Höenig & Kishimoto (2010) and the documentation on the dedicated webpage.
- Browsing the Stalevski et al. (2012) dataset we find similar results: if we assign the models the same bolometric luminosity as derived from the optical/UV data, then models peaking in the NIR might reproduce the observed NIR data, while models peaking in the MIR have a larger luminosity.
- Baes M., Verstappen J., De Looze I., Fritz J., Saftly W., Vidal Pérez E., Stalevski M., Valcke S., 2011, ApJS, 196, 22
- Barnett R., Warren S. J., Banerji M., McMahon R. G., Hewett P. C., Mortlock D. J. et al. 2015, A&A, in press
- Barth A.J., Martini P., Nelson C.H., Ho L.C., 2003, ApJ, 594, L95
- Bianchi S., 2008, A&A, 490, 461
- Bianchi S., Schneider R., 2007, MNRAS, 378, 973
- Bianchi S., & Xilouris E. M., 2011, A&A, 531, L11
- Bianchi S., 2013, A&A, 552, 89
- Carilli C. & Walter F. 2013, ARAA, 51, 105
- Cicone C., Maiolino R., Gallerani S., Neri R., Ferrara A., Sturm E. et al. 2015, A&A , 574, 14
- Deo R. P., Richards G. T., Nikutta R., Elitzur M., Gallagher S. C., Ivezic Z., Hines D. 2011, ApJ, 729, 108
- Dai Y. S., Bergeron J., Elvis M., Omont A., Huang J., Bock J. et al. 2012, ApJ, 753, 33
- de Bennassuti M., Schneider R., Valiante R., Salvadori S. 2014, MNRAS, 445, 3039
- Draine B. T., 2003, ARA&A, 41, 241
- Dwek E., Galliano F., Jones A.P., 2007, ApJ, 662, 927
- Elbaz D., Hwang H. S. , Magnelli B., Daddi E., Aussel H., Altieri B. et al. 2010, A&A, 518, L29
- Fan X. et al. 2003, ApJ, 125, 1649
- Fan X. et al., 2004, AJ, 128, 515
- Fanidakis N., Macció A. V., Baugh C. M., Lacey C. G., Frenk C. S. 2013, MNRAS, 436, 315
- Feltre A., Hatziminaoglou E., Fritz J., Franceschini A. 2012, MNRAS, 426, 120
- Feltre A., Hatziminaoglou E., Hernán-Caballero A., Fritz J., Franceschini A., Bock J. et al. 2013, MNRAS, 434, 2426
- Fioc M., Rocca-Volemrange B. 1997, A&A, 326, 950
- Gallerani S., Maiolino R., Juarez Y., Nagao T., Marconi A., Bianchi S., et al. 2010, A&A, 523, A85
- Haas M., Müller S. A. H., Chini R., Meisenheimer K., Klaas U., Lemke D. et al. 2000, A&A, 354, 453
- Haas M., Klaas U., Müller S. A. H., Bertoldi F., Camenzind M., Chini R. et al. 2003 A&A, 402, 87
- Hatziminaoglou E., Omont A., Stevens J. A., Amblard A., Arumugam V., Auld R. et al. 2010, A&A, 518, L33
- Hönig S., Proceedings of the Torus Workshop, 2012, held 5 December, 2012 at University of Texas, San Antonio, 2012, p. 157
- Hönig S., Kishimoto M., 2010, A&A, 523, A27
- Holwerda B. W., Bianchi S., Böker T., Radburn-Smith D., de Jong R. S., Baes M. et al. 2012, A&A, 541, L5
- Hines, D. C., Krause, O., Rieke, G. H., et al. 2006, ApJ, 641, L85
- Iwamuro, F., Kimura, M., Eto, S., et al. 2004, ApJ, 614, 69
- Jenkins E. B., 2009, ApJ, 700, 1299
- Jiang, L., Fan, X., Hines, D. C., et al. 2006, AJ, 132, 2127
- Leipski C., Meisenheimer K., Klaas U., Walter F., Nielbock M., Krause O. et al. 2010, A&A, 518, L34
- Leipski C., Meisenheimer K., Walter F., Besel M.-A., Dannerbauer H., Fan X. et al. 2013, ApJ, 772, 103L
- Leipski C., Meisenheimer K., Walter F., Klaas U., Dannerbauer H., De Rosa G. et al. 2014, ApJ, 785, 154
- Li Y., Hopkins P. F., Hernquist L., Finkbeiner D. P., Cox T. J., Springel V. et al. 2008, ApJ 678, 41
- Lutz D., Sturm E., Tacconi L. J., Valiante E., Schweitzer M., Netzer H. et al. 2008, ApJ, 684, 853
- Maiolino R., Gallerani S., Neri R., Cicone C., Ferrara A., Genzel R., 2012, MNRAS, 425, L66
- Netzer H., Lutz D., Schweitzer M., Contursi A., Sturm E., Tacconi L. J. et al. 2007, ApJ, 666, 806
- Rowan-Robinson M., 1995, MNRAS, 272, 737
- Richards G. T., Lacy M., Storrie-Lombardi L. J., Hall P. B., Gallagher S. C., Hines D. C. et al. 2006, ApJS, 166, 470
- Schweitzer M., Lutz D., Sturm E., Contursi A., Tacconi L. J., Lehnert M. D. et al. 2006, ApJ, 649, 79
- Stalevski M., Fritz J., Baes M., Nakos T., Popović L. Č., 2012, MNRAS, 420, 2756
- Valiante R., Schneider R., Bianchi S., Andersen A. C. 2009, MNRAS, 397, 1661
- Valiante R., Schneider R., Salvadori S., Bianchi S., 2011, MNRAS, 416, 1916
- Valiante R., Schneider R., Maiolino R., Salvadori S., Bianchi S., 2012, MNRAS, 427, L60
- Valiante R., Schneider R., Salvadori S., Gallerani S. 2014, MNRAS, 444, 2442
- Volonteri M. & Rees M. J. 2006, ApJ, 650, 669
- Volonteri M., 2012, Science, 337, 544
- Walter F., Carilli C., Bertoldi F., Menten K., Cox P., Lo K. Y., et al. 2004, ApJL, 615, L17
- Willott C.J., McLure R.J., Jarvis M.J., 2003, ApJ, 587, L15