The co-evolution of black hole growth and star formation from a cross-correlation analysis between quasars and the cosmic infrared background

We use CIB maps at 250, 350, and 500 $\micron$ (1200, 857, and 600 GHz, respectively) from the two Herschel surveys in the SDSS Stripe 82 region, i.e. the Herschel Stripe 82 Survey (HerS³³3HeRS maps and point source catalogues can be accessed from http://www.astro.caltech.edu/hers/.; Viero et al. 2013b) and the HerMES⁴⁴4The Herschel Multi-tiered Extragalactic Survey (HerMES; http://hermes.sussex.ac.uk) is a guaranteed time key project of the Herschel Space Observatory. HerMES maps and point source catalogues are publicly available at http://hedam.lam.fr/HerMES. Large-Mode Survey (HeLMS; Oliver et al. 2012). Maps were observed with the Spectral and Photometric Imaging Receiver (SPIRE) instrument (Griffin et al. 2010) aboard Herschel. HeRS covers 79 deg² along the SDSS Stripe 82 to a 5-$\sigma$ depth of 65.0, 64.5 and 74.0 mJy beam^-1 at 250, 350 and 500 $\micron$, respectively. HeLMS covers 270 deg² to a a 5-$\sigma$ depth of 64.0, 53.0 and 76.5 mJy beam^-1 at 250, 350 and 500 $\micron$. The joint HeRS and HeLMS areal coverage between -10^∘ and 37^∘ (RA) covers the full $\sim 150\,\rm deg^{2}$ subset of Stripe 82 that has the lowest level of Galactic dust emission (or cirrus) foregrounds (i.e., with $N_{\rm H}\lesssim 2\times 10^{21}\,\rm\,cm^{-2}$).

The SPIRE data, obtained from the Herschel Science Archive, were reduced using the standard ESA software and the custom-made software package, SMAP (Levenson et al. 2010, Viero et al. 2013b). Maps were made using an updated version of SMAP/SHIM, which is an iterative map-maker designed to optimally separate large-scale noise from signal. For more details, please refer to Viero et al. (2013b). The maps have a tangent plane (TAN) projection with pixel sizes of 6$\arcsec$, 8.33$\arcsec$ and 12$\arcsec$ at 250, 350, and 500 $\micron$, respectively. In comparison, the full width at half maximum (FWHM) of the SPIRE beams are 18.1$\arcsec$, 25.2$\arcsec$ and 36.6$\arcsec$ at 250, 350, and 500 $\micron$, respectively (Swinyard et al. 2010). The SPIRE maps are converted from their native unit of $\rm Jy\,beam^{-1}$ to $\rm Jy\,sr^{-1}$ by dividing them by the effective beam areas which are 0.9942, 1.765, and 3.730 $\times 10^{-8}$ steradians at 250, 350, and 500 $\micron$, respectively (SPIRE Observers’ Manual⁵⁵5http://herschel.esac.esa.int/Docs/SPIRE/html/spire_handbook.html). We note that no colour corrections from a flat-spectrum point-source calibration are applied, as they were shown in Viero et al. (2013a) to be negligible. Fig. 1 shows the three-colour SPIRE images of the HeRS and HeLMS regions. We can clearly see white foreground clouds of Galactic cirrus emission, which have been shown to correlate with HI emission from GASS 21-cm data in Viero et al. (2013b).

In this paper, we also include SPIRE maps of a subset of the HerMES level 5 and level 6 fields⁶⁶6The HerMES fields are organised, according to area and depth, into levels 1 through to 7, with level 1 fields being the smallest and deepest, and level 7 being the widest and shallowest. HeLMS is the only level 7 field., which have significant overlaps with our SDSS quasar samples. In particular, the HerMES Boötes field (11.3 deg²) and the Lockman-SWIRE field (15.2 deg²) overlap with the SDSS DR7 QSO sample. The Boötes field and the XMM-LSS field (21.6 deg²) overlap with the SDSS DR9 QSO sample. Fig. 2 shows the three-colour SPIRE images of the Boötes, Lockman-SWIRE and XMM-LSS field. Out of these three fields, the XMM-LSS field is the most contaminated by cirrus emission, but all three fields have significantly less cirrus contamination than the HeRS or HeLMS field. In addition, the level 5 and 6 fields are significantly deeper than HeRS or HeLMS. The 5-$\sigma$ noise level in the XMM-LSS and Boötes field is 25.8, 21.2 and 30.8 mJy at 250, 350 and 500 $\micron$, respectively. The Lockman-SWIRE field is even deeper and the 5-$\sigma$ noise level is 13.6, 11.2 and 16.2 mJy at 250, 350 and 500 $\micron$, respectively.

The CIB at sub-mm wavelengths is dominated by emission from DSFGs over a wide range of redshift, while other potential sources of signal, e.g., the cosmic microwave background (CMB), Sunyaev-Zel’dovich (SZ) effect, and radio galaxies, are subdominant. One exception is the emission from Galactic dust (cirrus) which cannot be ignored at FIR/sub-mm wavelengths. Therefore, the SPIRE map signal can be modelled, to first order, as the sum of the following three terms: the emission from DSFGs; Galactic cirrus emission; and instrument noise. Thus we can write

We expect that only the emission from DSFGs have intrinsic correlation with the quasars. This is because some quasars could be strong IR/sub-mm emitters from dust heated by the accreting black hole, as well as star formation in the host galaxy (e.g. Lutz et al. 2008; Shi et al. 2009; Clements et al. 2009; Serjeant et al. 2010; Bonfield et al. 2011; Dai et al. 2012). But additionally DSFGs are expected to occur as satellite galaxies in the dark matter halos where quasars are found⁷⁷7The satellite fraction of quasars (i.e. quasars hosted by satellite galaxies in the halo) is estimated to be at most a few percent (Kayo & Oguri 2012; Richardson et al. 2012; Shen et al. 2013). Therefore, most quasars are expected to be central galaxies.. A last effect is that dark matter halos correlated with the quasar-hosting halo are also expected to host DSFGs. The SPIRE instrumental noise should not correlate with quasars and therefore only contributes to the noise in the measured cross-correlation signal between quasars and the CIB. Similarly, Galactic dust emission should not have any intrinsic correlation with the quasars. However, due to the selection of the quasars, there might be (anti-)correlation between the observed quasar number density field and the intensity of the cirrus emission. There are two effects of Galactic dust on the selection of QSOs. The first effect is reddening which means QSOs with extreme reddening will not be selected as targets for spectroscopy. The second effect is extinction, which can cause the QSOs to drop out of the target list even before they reach the colour selection stage. QSO populations with steeper number counts will be more sensitive to extinction.

To explicitly check the potential correlation between the quasar samples and Galactic cirrus emission, we make use of the 100 $\micron$ maps from the Improved Reprocessing of the IRAS (Neugebauer et al. 1984) Survey (IRIS; Miville-Deschenes & Lagache 2005), a data set which corrects the original plates for calibration, zero level and striping problems. The IRIS 100 $\micron$ maps in the HerMES level 5 and level 6 fields are discussed in Viero et al. (2013) and the IRIS 100 $\micron$ maps in HeRS and HeLMS are discussed in Hajian et al. (in prep.). The 100 $\micron$ maps in all of our fields are shown in Appendix A. Here we have assumed that the dominant signal in the 100 $\micron$ maps comes from Galactic cirrus emission, not the DSFGs comprising the CIB. In other words, these 100 $\micron$ maps are used as a Galactic cirrus tracer in our fields⁸⁸8Properties of Galactic cirrus (such as column density, dust emissivity and dust temperature) can vary a lot from field to field. In some regions, the IRIS 100 $\micron$ data may not be used as a reliable Galactic cirrus tracer as a result of significant contamination from the CIB (Pénin et al. 2012).. Visual inspection of the maps in Appendix A confirms that extended diffuse Galactic cirrus emission are the dominant structures (long filamentary chains) in the 100 $\micron$ maps, especially in the HeRS and HeLMS regions. We have also computed the mean value of the 100 $\micron$ maps in the HeLMS, HeRS, XMM-LSS, Boötes and Lockman-SWIRE field to be 3.03, 2.84, 2.30, 1.42 and 1.09 MJy sr^-1, respectively. If we adopt the DIRBE measurement of the cosmic IR background at 100 $\micron$ which is $14.4\pm 6.3$ nW m^-2 sr^-1 or equivalently $0.48\pm 0.21$ MJy sr^-1 (Dole et al. 2006), then we can estimate that the fractional contribution of the CIB to the 100 $\micron$ maps in the HeLMS, HeRS, XMM-LSS, Boötes and Lockman-SWIRE field is 16%, 17%, 21%, 34% and 44%, respectively.

For our cross-correlation analysis, we use quasars from both the SDSS DR7 (Abazajian et al. 2009) and the SDSS-III DR9 (Ahn et al. 2012) spectroscopic surveys.

The cross-correlation of the two fields of interest $\rho_{1}$ (representing the quasar number density field) and $\rho_{2}$ (representing the CIB) is

where $\overline{\rho}_{1}$ is the mean number density of QSOs and $\overline{\rho}_{2}$ is the mean value of the CIB. The SPIRE maps, $\rho^{\prime}_{2}$, can be used to estimate $\rho_{2}$ but have been mean-subtracted, i.e. the maps have a mean of zero $\overline{\rho^{\prime}_{2}}=0$. The two quantities are thus related $\rho_{2}=\rho^{\prime}_{2}+C$, where $C$ is the unknown DC level of the map. To proceed we assume a value for $C$ and add this to the map⁹⁹9We add a constant $C$ to our mean-subtracted maps so that the minimum value in the map is a positive number. and will show that the choice of $C$ only affects the normalisation of our results. We define a re-normalised cross correlation $\xi^{\prime}$:

By construction, Eq. (3) gives the correlation between the QSO number density field and the FIR surface brightness field (with a normalisation constant). It is also equivalent to the stacked SPIRE map signal at the positions of the QSOs minus the mean signal of the map (which is equal to zero by design). This re-normalised cross-correlation, equivalent to stacking, is no longer dimensionless, but has the same unit as the map.

To estimate the two-point angular cross-correlation function $\xi$ between the quasar samples and the SPIRE maps, we adapt the Landy & Szalay (1993) auto-correlation estimator¹⁰¹⁰10The Landy & Szalay estimator is normally used to estimate correlation functions for point source catalogues. Here we use it to calculate the cross-correlation between a point source catalogue and an image.,

Here $D_{1}$ represents the pair-counts in the real quasar sample, $D_{2}$ the CIB at a given wavelength, $R_{1}$ the simulated quasar sample with the same angular selection effects as the real quasar sample, and $R_{2}$ the simulated CIB with the same noise properties as the real SPIRE maps. $D_{1}D_{2}(\theta)$ represents the total emission in the CIB map at a separation of $\theta$ from real quasars, $D_{1}R_{2}(\theta)$ represents the total emission in the simulated CIB map at separation $\theta$ from all real quasars, $R_{1}D_{2}(\theta)$ represents the total emission in the CIB map at separation $\theta$ from all simulated quasars, $R_{1}R_{2}(\theta)$ represents the total emission of the simulated CIB map at separation $\theta$ from all simulated quasars. Note that $D_{1}D_{2}$ is normalised by the total number of real quasars ($N_{\rm Q}$) times the total number of real SPIRE map pixels ($N_{\rm P}$), $D_{1}R_{2}$ is normalised by $N_{\rm Q}$ times the total number of simulated SPIRE map pixels (which is the same as $N_{\rm P}$), $R_{1}D_{2}$ is normalised by the total number of simulated quasars ($N_{S}$) times $N_{\rm P}$, and $R_{1}R_{2}$ is normalised by $N_{S}\times N_{\rm P}$.

The CIB map and simulated CIB map are estimated from the mean-subtracted SPIRE maps and the simulated SPIRE maps by adding the same constant number $C$ to represent the mean CIB. As in Eq. (3) we then estimate a re-normalised cross-correlation signal. So, our actual estimator is¹¹¹¹11This is similar to stacking, shown in Appendix B. We have checked that the cross-correlation signal calculated using Eq. (5) is almost identical to the stacking result.

Here $D^{\prime}_{2}=D_{2}-C$ is the real SPIRE map, and $R^{\prime}_{2}=R_{2}-C$ the simulated SPIRE map. We have checked that the cross-correlation result (calculated using Eq. (5)) as well as its error does not depend on the value of $C$ at all. Although this is a reasonable choice of estimator, we have not investigated whether it is optimal.

The simulated SPIRE are constructed by randomising the pixels of the real map. It allows us to take into account of the effect of not only instrumental noise, but also confusion noise and cirrus. The random quasar samples are generated according to the angular completeness maps in the overlapping area between the SDSS DR7 and SDSS-III DR9 quasar sample and the Herschel surveys. As mentioned before, we assume the angular completeness fraction is uniform for the DR7 quasars. So, we simply generate uniformly distribution points as the random DR7 quasar sample. Admittedly, the real angular completeness level for the DR7 quasars might vary across our fields. However, as what we are concerned with in this paper is the cross-correlation signal between the quasars and the CIB measured by Herschel-SPIRE, the variation in the DR7 QSO completeness is not expected to bias our results as long as the completeness is uncorrelated with the noise properties of the SPIRE maps. For the DR9 QSOs in Boötes and XMM-LSS, we use the angular completeness masks used in White et al. (2012) and the MANGLE software (Swanson et al. 2008) to calculate the angular completeness of the survey, which is the fraction of quasar targets in a sector for which a spectrum was obtained. For the DR9 QSOs in the HeRS and HeLMS S82 region, we simply generate uniformly distributed points as the simulated QSO sample. The total number of simulated quasars $N_{S}$ is generally over 25 times larger than the total number of real quasars $N_{\rm Q}$.

We estimate errors on the two-point angular cross-correlation function by using 100 bootstrap samples of the quasar sample in each field. The bootstrap samples which have the same size as the real QSO sample are generated by sampling with replacement the real QSO sample. The dispersion in the signal between the bootstrap quasar samples and the SPIRE maps or the IRIS maps is taken as the error on the measured cross-correlation signal (Norberg et al. 2009).

We expect that there should be no intrinsic correlation between the quasars and Galactic cirrus. However, the observed quasar number density field could potentially (anti-)correlate with the intensity of the Galactic cirrus emission as the SDSS quasars are selected based on their optical flux and colour. In Fig. 5, we plot the cross-correlation signal between the quasar samples and the Galactic cirrus emission traced by the IRIS 100 $\micron$ maps. The horizontal dotted line corresponds to no correlation. The vertical dashed line corresponds to the FWHM of the IRAS beam (Miville-Deschênes et al. 2002). Data points below the scale of the IRAS beam are strongly correlated. One quasar in Boötes happens to lie on top of a very bright point source in the IRIS map, which biases our result (the empty red circles in the top panel in Fig. 5). The filled red circles in the top panel in Fig. 5 correspond to the cross-correlation between the DR7 QSOs, excluding the aforementioned quasar, and the 100 $\micron$ map in Boötes. The reduced $\chi^{2}$ (the degrees of freedom DoF is 8), computed using data points above the scale corresponding to the FWHM of the IRAS beam, is 0.25, 0.24, 0.14 and 0.95 in the Boötes, Lockman-SWIRE, HeRS and HeLMS region, respectively. For the DR9 QSOs, the reduced $\chi^{2}$ (DoF = 8), computed using data points above the scale corresponding to the FWHM of the IRAS beam, is 0.24, 0.92, 0.29 and 2.39 in the Boötes, XMM-LSS, HeRS S82 and HeLMS S82 region, respectively. The cross-correlation signal in the HeRS S82 region is consistent with zero, but is systematically negative in the HeLMS S82 region.

The nominal $i$-band limit of the BOSS DR9 QSO sample is 21.8, which is deeper than the SDSS imaging limit of 21.3. So, in principle, it is possible that DR9 QSOs suffer more from incompleteness due to Galactic dust extinction. We calculate the cross-correlation signal between “uniform $>0$” DR9 QSOs with different $i$-band magnitude cuts (corrected for Galactic extinction), namely 21.2, 21.0 and 20.0, and IRIS 100 $\micron$ maps in the HeRS and HeLMS S82 region. The cross-correlation between different quasar sub-samples and the 100 $\micron$ maps are consistent with each other. In the HeLMS S82 region, the cross-correlation signals between different sub-samples and the 100 $\micron$ maps are still systematically below zero. This test suggests that incompleteness due to Galactic dust extinction is not the likely cause for the anti-correlation seen in the HeLMS S82 region.

In Fig. 6, we plot the cross-correlation signal between DR9 QSOs with different $i$-band magnitude cuts and IRIS 100 $\micron$ maps in the HeRS and HeLMS S82 region. We have not used the “uniform $>0$” criterion in the samples plotted in Fig. 6. The anti-correlation in the HeLMS S82 region has now disappeared. The reduced $\chi 2$ (DoF = 8), computed using data points above the scale corresponding to the FWHM of the IRAS beam, is 1.29, 0.84, 0.35 for DR9 QSOs with $i$-band cut of 21.2, 21.0 and 20.0 respectively in the HeRS region. The reduced $\chi 2$ (DoF = 8), computed using data points above the scale corresponding to the FWHM of the IRAS beam, is 0.84, 0.68, 0.78 for DR9 QSOs with $i$-band cut of 21.2, 21.0 and 20.0 respectively in the HeLMS region. This makes sense as QSOs with “uniform $>0$” correspond to the XDQSO selection which is not the main selection method used in the S82 region. We finally select the DR9 QSOs with $i$-band cut of 21.0 in the HeRS and HeLMS S82 region, but not requiring “uniform $>0$”.

In Fig. 7, we plot the combined cross-correlation signal between the quasars and the IRIS 100 $\micron$ maps, using inverse variance weighting. The filled circles represent combined cross-correlation estimates over all fields and the empty circles represent combined cross-correlation estimates using only the HeRS and HeLMS fields. We can see clearly that the two types of combined estimates are consistent with each other. The cross-correlation between either DR7 or DR9 QSOs and the 100 $\micron$ maps are consistent with zero. As we have discussed in Section 2.1, the dominant signal in the 100 $\micron$ maps are due to Galactic cirrus emission (especially in the HeRS and HeLMS regions). We conclude that the correlations between the final selected QSO samples and Galactic cirrus are consistent with zero.

In the left column of Fig. 8, we plot the cross-correlation signal between the SDSS DR7 QSOs (in the redshift range $z=[0.15,2.5]$) and the SPIRE maps at 250, 350 and 500 $\micron$ in HeRS, HeLMS, Boötes and Lockman-SWIRE. The cross-correlation signal in four different fields are consistent with each other. On small scales, the cross-correlation signal is dominated by the mean sub-mm emission of the quasars themselves, as indicated by the SPIRE point spread function (PSF) in the corresponding waveband. No fitting to the measurement is carried out at this stage, so the amplitude of the PSF is only approximate. On larger scales, significant excess signal beyond the SPIRE PSF is seen in all three wavebands in all fields, indicating that there is a population of DSFGs correlated with the quasar population. They could be mainly satellite galaxies residing in the same dark matter halo as the host galaxy of the quasar and galaxies in other halos which are correlated with the halo hosting the quasar.

In the right column Fig. 8, we plot the cross-correlation signal between the SDSS-III DR9 QSOs (in the redshift range $z=[2.2,3.5]$) and the SPIRE maps at 250, 350 and 500 $\micron$ in HeRS, HeLMS, Boötes and XMM-LSS. Again, on small scales, the cross-correlation signal is dominated by the mean sub-mm emission of the quasars themselves, as indicated by the SPIRE beams. Compared to the DR7 quasar sample, the amplitude of the SPIRE emission of the DR9 quasar sample is slightly smaller in all three bands. On larger scales, excess signal beyond the PSF is seen in all fields at all wavelengths apart from the XMM-LSS field at 350 and 500 $\micron$. The dispersion in the cross-correlation signal between different fields is also larger for the DR9 QSOs than the DR7 QSOs. It could be due to the surface density of the DR9 quasars which is around a factor of two lower than the surface density of the DR7 quasars (except in Boötes).

In Fig. 9, we plot the combined cross-correlation signal between the QSOs and the CIB at 250, 350 and 500 $\micron$, using inverse variance weighting. The filled circles represent combined cross-correlation estimates over all fields and the empty circles represent combined cross-correlation estimates using only the HeRS and HeLMS fields. The two different types of combined cross-correlation signal are consistent with each other and exhibit significant excess signal beyond the SPIRE PSF at all wavelengths. In Fig. 10, we compare the combined cross-correlation signal (using all fields) between the SDSS QSOs and the SPIRE maps at all three wavelengths. For both DR7 and DR9 quasars, the sub-mm emission of the quasars themselves decreases by almost an order magnitude from 250 to 500 $\micron$, while the change in the correlated sub-mm emission is much smaller from 250 to 500 $\micron$. This is partly to do with the broadening of the SPIRE beam from 250 to 500 $\micron$ which affects point sources (i.e. the sub-mm emission of the quasars themselves) more than extended structures (i.e. the correlated signal from DSFGs). Another possibility is that the QSOs have a higher effective dust temperature than that of the correlated DSFGs.

The halo model of the CIB anisotropies presented in this paper is identical to the model presented in Viero et al. (2013) and Hajian et al. (in prep), which is in turned developed based on Shang et al. (2013). For the sake of completeness, we outline the main elements of the CIB model here. For more details, please refer to Shang et al. (2013), Viero et al. (2013) and Hajian et al. (in prep).

The power spectra of the CIB sourced by infrared galaxies is composed of a Poisson noise (or shot noise) term and a clustered term,

The Poisson noise term in the CIB power spectrum, independent of angular scale, can be calculated from the multi-dimensional number counts of the infrared sources,

The clustered term in the angular power spectrum of the CIB is the flux averaged projection of the three-dimensional (3D) spatial power spectrum,

Here $dV_{c}$ is the comoving volume element, and $\frac{dS_{\nu}}{dz}$ is the redshift distribution of the integrated flux at a given frequency which can be calculated as

The 3D spatial power spectrum of CIB sources can be decomposed into a 1-halo and 2-halo term which describes the small-scale clustering due to pairs of galaxies in the same dark matter halo and the large-scale clustering due to pairs of galaxies in different halos respectively,

The 3D non-linear, 1-halo term of the CIB clustered term can be calculated as

Here $f_{\nu}^{\rm cen}(M,z)$ and $f_{\nu}^{\rm sat}(M,z)$ are the luminosity-weighted number of central and satellite galaxies as a function of redshift and halo mass, $\bar{j}_{\nu}(z)$ is the mean comoving specific emission coefficient as a function of redshift, $\frac{dN}{dM}(z)$ is the redshift-dependent halo mass function (taken from Tinker et al. 2008), $u_{\rm gal}(k,z|M)$ is the normalised Fourier transform of the galaxy density distribution within a halo (assumed to follow the NFW profile in this paper). Here luminosity-weighted number of central and satellite galaxies can be derived as

and

where $m$ is the subhalo mass at the time of accretion and $\frac{dn}{dm}(M,z)$ is the sub-halo mass function (taken from Tinker & Wetzel 2010). Therefore, the mean comoving specific emission coefficient is

We assume a log-normal relation between halo mass and infrared luminosity,

where $M_{\rm peak}$ and $\sigma_{L/M}$ is peak halo mass scale and $1\sigma$ range of the specific infrared luminosity per unit mass, and $L_{0}$ is the overall infrared luminosity to halo mass normalisation factor. A lower limit on the halo mass $M_{\rm min}$ is applied to the halo mass - infrared luminosity relation. To describe the spectral energy distribution (SED) of the infrared galaxy population, we use a single modified blackbody with two parameters, i.e. dust temperature $T_{\rm dust}$ and emissivity $\beta$. In addition, the infrared luminosity is assumed to evolve with redshift as $(1+z)^{\eta}$ over the range $0<z<2$ followed by a plateau at $z>2$. The 2-halo term of the CIB clustering power spectrum can be calculated as

where $P_{\rm lin}(k,z)$ is the linear dark matter power spectrum and

Here $b(M,z)$ is the linear large-scale bias. We use the prescription for the halo bias from Tinker, Wechsler & Zheng (2010) in this paper.

Due to the limited constraining power of the measured cross-correlation signal between the QSOs and the SPIRE maps (our measurements are quite noisy especially on large scales), in this paper we have used the same parameters of the CIB halo model as in Viero et al. (2013) which fits the auto- and cross-correlation power spectra of the SPIRE bands. The exact values of the parameters used in the CIB model are minimum halo mass $M_{\rm min}=10^{10.1}M_{\odot}$, dust temperature $T_{\rm dust}=24.2$ K, peak halo mass $M_{\rm peak}=10^{12.1}M_{\odot}$, $\sigma^{2}_{L/M}=0.38$, infrared luminosity to halo mass normalisation $L_{0}=10^{-1.71}L_{\odot}/M_{\odot}$, dust emissivity $\beta=1.45$ and luminosity evolution $\eta=2.19$.

Many studies have used the clustering measurements of quasars to constrain their halo occupation distributions. Generally, quasars are found to inhabit dark matter halos of masses around a few times $10^{12}M_{\odot}$ (Croom et al. 2005; Da Angela et al. 2008; Shen et al. 2009; Ross et al. 2009; White et al. 2012; Shen et al. 2012). The clustering strength of quasars seems to depend very weakly on QSO luminosity or redshift. It indicates that there is a poor correlation between halo mass and the instantaneous quasar luminosity. Richardson et al. (2012) found tentative evidence for increasing halo mass scale for quasars with increasing redshift. They found that at $z\sim 1.4$ the median halo mass hosting central quasars is around $4.1^{+0.3}_{-0.4}\times 10^{12}M_{\odot}$ while at $z\sim 3.2$ the median halo mass increases to $14.1^{+5.8}_{-6.9}\times 10^{12}M_{\odot}$. Most quasars are central galaxies according HOD analysis of quasar clustering statistics. The satellite fraction of quasars is generally found to be of the order a few percent at most, although the exact satellite fraction estimated from different studies vary significantly from each other. For example, Shen et al. (2013) using a cross-correlation study between the DR7 quasars and DR10 BOSS galaxies estimated the satellite fraction to be $0.068^{+0.034}_{-0.023}$. A similar halo model applied to the auto-correlation of the DR7 quasars combined with a binary quasar sample (Hennawi et al. 2006) in Richardson et al. (2012) inferred the satellite fraction to be $(7.4\pm 1.3)\times 10^{-4}$.

In this paper, we use a simple halo model for the quasars, widely used in the literature. The quasar mean halo occupation number is composed of central and satellite components,

A softened step function is assumed to describe the central component,

and a rolling-off power law for the satellite component,

There are five free parameters in the quasar HOD: $M_{\rm min}$ is the characteristic halo mass scale at which on average half of the dark matter halos host one quasar as the central galaxy; $\sigma_{\log M}$ is the width of the softened step function; $M_{1}$ is the approximate halo mass scale where on average halos host one satellite quasar; $\alpha$ is the power-law index; $M_{\rm cut}$ is the halo mass scale below which the number of the satellite quasars decreases exponentially.

In principle, we should also multiply a free parameter to the quasar HOD $N_{q}(M)$ to describe the duty cycle of quasars, i.e. only a certain fraction ($<1$) of the halos actually host quasars. However, if we assume that quasar duty cycle is independent of halo mass, then it will cancel the duty cycle parameter which will appear in the QSO number density $\bar{n}$ in Eq. (21) and (22). So, the quasar HOD equations (Eq. 18 – 20) should be interpreted as the number of quasars in halos of mass $M$ which host quasars (not all halos of mass $M$).

Combining the CIB halo model described in Section 4.1 and the quasar halo model described in Section 4.2, we can model the cross-correlation between the quasars and the CIB. The 1-halo term of the cross-correlation power spectrum between quasars and the CIB can be modelled as

And the 2-halo term of the cross-correlation power spectrum between quasars and the CIB can be modelled as

where

and

Under Limber’s approximation (Limber 1953), the projected angular cross-correlation power spectrum between QSOs and CIB is related to the three-dimensional spatial power spectrum,

where $dV_{c}/dz$ is the comoving volume element per unit along the redshift axis, $dS/dz$ is the redshift distribution of the cumulative CIB flux, and $dN^{\rm Q}/dz$ is the normalised redshift distribution of the QSO sample.

The angular correlation function can be expressed in terms of the angular power spectrum

where $J_{0}(x)$ is the zeroth order Bessel function. Furthermore, the angular correlation function needs to be convolved with the 2D SPIRE beam, which can be described as a 2D Gaussian function. Therefore, the final beam convolved angular cross-correlation between the QSOs and the SPIRE maps can be written as

where $\sigma$ is the standard deviation of the Gaussian beam.

As mentioned earlier, we have used the same parameters of the CIB halo model as in Viero et al. (2013) which fits the auto- and cross-correlation power spectra of the SPIRE bands. So, in this paper we are only varying the parameters related to the quasar population to fit the cross-correlation signal between the quasars and the CIB. There are in total eight free parameters in our halo model, which are the five parameters describing the halo occupation distributions of the quasars and the three parameters related to the amplitude of the mean sub-mm emissions of the quasars themselves at 250, 350 and 500 $\micron$.

In Fig. 11, we plot the combined cross-correlation signal, over all fields, between the SDSS QSOs (left column: DR7 QSOs; right column: DR9 QSOs) and the CIB (top: 250 $\micron$; middle: 350 $\micron$; bottom: 500 $\micron$) and the best-fit from the halo model. The dashed line in each panel corresponds to the best-fit mean sub-mm emission of the quasars themselves. The blue solid line corresponds to the best-fit one-halo term of the cross-correlation signal convolved with the 2D SPIRE beam, the red solid line corresponds to the best-fit two-halo term convolved with the SPIRE beam, and the black solid line corresponds to the total cross-correlation signal (i.e. the sum of the sub-mm emission of the quasars themselves and the emission from the correlated DFSGs) convolved with the SPIRE beam. The black circles are the measured cross-correlation signal. With a median redshift $\left<z\right>$ of 1.4, the mean sub-mm flux densities of the DR7 quasars is $11.12^{+1.56}_{-1.37}$, $7.08^{+1.63}_{-1.33}$ and $3.63^{+1.35}_{-0.98}$ mJy at 250, 350 and 500 $\micron$ respectively, while the mean sub-mm flux densities of the DR9 quasars ($\left<z\right>=2.5$) is $5.69^{+0.72}_{-0.64}$, $4.98^{+0.75}_{-0.65}$ and $1.76^{+0.50}_{-0.39}$ mJy at 250, 350 and 500 $\micron$ respectively. Using a modified blackbody with a dust emissivity of 1.45 combined with the mid-infrared power-law ($f_{\nu}\propto\nu^{-2}$) to describe the SED of the quasars (i.e. the same assumptions which we made regarding the SED of the DSFGs in Section 4.1), we can derive the best-fit dust temperature and total infrared luminosity for the quasars. For the DR7 quasars, the best-fit dust temperature and infrared luminosity is $39.01^{+11.00}_{-6.25}$ K and $\log_{10}L_{\rm IR}(L_{\odot})=12.38^{+0.25}_{-0.15}$ respectively. For the DR9 quasars, the best-fit dust temperature and infrared luminosity is $52.30^{+7.78}_{-5.56}$ K and $\log_{10}L_{\rm IR}(L_{\odot})=12.82^{+0.12}_{-0.09}$ respectively. Due to the degeneracy between dust temperature and dust emissivity, the best-fit temperature decreases with increasing values of dust emissivity but the best-fit infrared luminosity is almost unaffected. Therefore, the higher redshift SDSS-III DR9 quasars are on average brighter than the SDSS DR7 quasars in the infrared, possibly due to elevated star-formation activity and/or black hole accretion. Another possibility is that for the higher redshift DR9 quasars, we are starting to probe shorter wavelength range where the contribution of the accreting BH to the total infrared luminosity is larger. The dust temperature of the satellite galaxies around DR7 and DR9 quasars is fixed at 24.2 K in the halo model (see Section 4.1) as we use the same CIB halo model parameters as in Viero et al. (2013) which fits the auto- and cross-frequency channel power spectra at the SPIRE wavelengths. As we have used the same SED assumptions to describe the QSOs and the DSFGs, we can directly compare the estimated temperatures for both populations. Both SDSS DR7 and SDSS-III DR9 QSOs are found to be hotter than the correlated DSFGs, possibly due to hotter dust heated by the accreting BH.

In Fig. 12, we plot the mean halo occupation distribution of SDSS DR7 and DR9 QSOs, decomposed into its central and satellite components. The shaded region indicate the $68\%$ confidence intervals. The red, blue and green solid lines represent the best-fit value for $M_{\rm min}$ (the characteristic halo mass scale at which on average half of the halos host one central quasar), $M_{\rm cut}$ (the halo mass scale below which the HOD for the satellite quasars decays exponentially) and $M_{1}$ (the halo mass mass for hosting on average one satellite quasar) respectively. For the DR7 QSOs, the host halo mass scale for the central and satellite quasars is $M_{\rm min}=10^{12.36\pm 0.87}$ M_⊙ and $M_{1}=10^{13.60\pm 0.38}$ M_⊙, respectively. For the DR9 QSOs, the host halo mass scale for the central and satellite quasars is $M_{\rm min}=10^{12.29\pm 0.62}$ M_⊙ and $M_{1}=10^{12.82\pm 0.39}$ M_⊙, respectively. The typical halo environment of the DR7 and DR9 QSOs we find in this study is similar to previous studies using QSO auto-correlation statistics (e.g., Shen et al. 2009; White et al. 2012; Shen et al. 2012; Richarson et al. 2012).Table 2 and Table 3 list the best-fit values, uncertainties and correlations of the parameters in the halo model for the DR7 and DR9 quasars. From the HODs, we find that the satellite fraction (the ratio of satellite quasars to the total number of quasars) is $0.008^{+0.008}_{-0.005}$ and $0.065^{+0.021}_{-0.031}$ for the DR7 and DR9 quasars, respectively. Our estimates are broadly consistent with other studies in the literature, i.e. the QSO satellite fraction is found to be around a few percent at most. Thus the majority of the DR7 and DR9 QSOs are expected to be central galaxies in dark matter halos and the correlated DSFGs in the one-halo term are mostly satellite galaxies in the same halo.