A CROSS-CORRELATION ANALYSIS OF ACTIVE GALACTIC NUCLEI AND GALAXIES USING VIRTUAL OBSERVATORY: DEPENDENCE ON VIRIAL MASS OF SUPERMASSIVE BLACK HOLE

It is thought that most galaxies harbor supermassive black holes (SMBHs) in their centers, and that gas accretion onto an SMBH is the energy source of active galactic nuclei (AGNs; e.g., Richstone et al. 1998; Ferrarese & Ford 2005). There are strong correlations between the mass of an SMBH and the observational properties of its host galaxy such as velocity dispersion and stellar mass of the bulge (e.g., Magorrian et al. 1998; Ferrarese & Merritt 2000).

Growth of the BH mass is strongly coupled to the evolution of the host galaxy, but the process of "co-evolution" is still unknown. It is thought that BHs grow by accretion of gas and/or merger of BHs, but the physical mechanisms of gas inflow toward central regions of galaxies and coalescence of binary BHs are not yet revealed. In the standard hierarchical structure formation framework, major mergers of galaxies should have played important roles for the evolution of galaxies, growth of SMBHs, and AGN activity. Thus, investigating the clustering of the galaxies around AGNs is crucial to understand the evolution of SMBHs and galaxies.

Recent large-scale surveys, such as the Sloan Digital Sky Survey (SDSS), provide observational samples of over 100, 000 AGNs (Schneider et al. 2010). The auto-correlation function of AGNs was studied using the SDSS sample (Shen et al. 2007, 2009; Ross et al. 2009). The clustering of galaxies around AGNs in the areas of deep surveys was also investigated by some authors (e.g., Coil et al. 2009; Mountrichas et al. 2009, 2013). They showed that the cross-correlation function between AGNs and galaxies is similar to the auto-correlation of luminous red quiescent galaxies. Hickox et al. (2009) found that radio selected AGNs are strongly clustered, and that infrared selected AGNs are more weakly clustered than optically selected AGNs. Recently, cross-correlation between AGNs and galaxies was studied using large samples. Donoso et al. (2010) computed the cross-correlation between ∼14,000 radio-loud AGNs at z = 0.4–0.8 and a reference luminous galaxy sample, and compared the clustering amplitude of radio galaxies with that of ∼7000 SDSS quasars. They argued that radio-loud AGNs are clustered more strongly than radio-quiet ones. Krumpe et al. (2012) investigated clustering of galaxies around ∼3000 X-ray selected AGNs and ∼8000 optically selected AGNs at z < 0.5. No significant difference was found between X-ray selected and optically selected broad-line AGNs.

In order to understand the interaction between the growth of SMBHs and their surrounding environment, it is important to investigate the dependence of clustering amplitude on physical properties of BHs. Shen et al. (2009) studied the dependence of the two-point auto-correlation function of quasars on luminosity, BH mass, color, and radio loudness, and found weak or no dependence on virial BH mass. The cross-correlation function between AGNs and galaxies will achieve smaller uncertainties in clustering measurement because it has many more pairs at a given separation compared with the auto-correlation function of AGNs. Donoso et al. (2010) found a positive dependence of the cross-correlation amplitude on stellar mass M_*, but their sample is radio-loud AGNs within a narrow range of stellar mass (10¹¹ M_☉ < M_* < 10¹² M_☉). There is no previous study to investigate the dependence of the cross-correlation function between AGNs and galaxies on BH mass over a wide mass range.

In this study, we investigate the dependence of clustering amplitude of galaxies around AGNs on the BH mass to reveal a relation between the mass growth of SMBHs and the large-scale environment of their host galaxies. We have collected observational data of a large number of AGNs and galaxies by using the Japanese Virtual Observatory (JVO),¹ and computed the cross-correlation function between AGNs and galaxies. Thanks to the large sample covering a large area of the sky, we were able to perform the clustering analysis in a manner free from the cosmic variance. To investigate BH mass dependence over a wide mass range (≳ 2 dex), we also use the sample of less massive BHs by Greene & Ho (2007) in addition to the SDSS quasar catalog (Shen et al. 2011), as described below.

Throughout this paper, we use the following cosmological parameters: Ω_M = 0.3, Ω_Λ = 0.7, h = 0.7. All magnitudes are given in the AB system. All of the distances are measured in comoving coordinates.

The AGN sample was extracted from two catalogs by Shen et al. (2011) and Greene & Ho (2007), both of which contain the estimated virial mass, M_BH, of the SMBHs. Shen et al. (2011) derived the virial mass of the BHs for 105,783 quasars in the SDSS DR7 quasar catalog (Schneider et al. 2010). About half of their samples were uniformly selected by the criteria described in Richards et al. (2002), and the remaining samples were selected by a variety of earlier algorithms or serendipitous selections.

Greene & Ho (2007) derived BH masses for ∼8500 active galaxies at z < 0.35 based on SDSS DR4 spectra. They also analyzed spectra for objects classified as galaxies, not only quasars. They extracted the AGN components from spectra of galaxies and estimated BH masses. Their catalog contains more objects than the sample of Shen et al. (2011) for M_BH ≲ 10⁷ M_☉.

They both estimated the virial mass of BHs by means of the observed FWHM of the emission lines of Hα, Hβ, or Mg ii and the continuum luminosity at the lines. Shen et al. (2011) used Hβ estimates for z < 0.7 and Mg ii estimates for z ⩾ 0.7 as a fiducial virial mass estimate but also gave estimates based on other lines in their catalog when the lines were detected. For the samples of Shen et al. (2011), we used their fiducial mass estimate in this paper. Greene & Ho (2007) used Hα estimates. We have found that there is a systematic difference of ∼0.5 dex between the virial masses estimated in the two catalogs, mainly because they used different parameter values in the virial mass estimator. Figure 1 shows the estimated mass of 2139 AGNs which are registered in both catalogs. For the samples of Greene & Ho (2007), we have recomputed the virial mass by means of the FWHM and luminosity in Greene & Ho (2007) with the parameter values in Shen et al. (2011) for the Hα line. For the recomputed virial mass, the systematic difference is decreased as shown in the bottom panel of Figure 1. The recomputed masses by means of the data of Greene & Ho (2007), however, are still ∼0.2 dex smaller than the estimated masses in the catalog of Shen et al. (2011) on average. Half of this remaining systematic difference is because the FWHM and luminosity values were derived from the Hα line for Greene & Ho (2007) but from the Hβ line for Shen et al. (2011). The other half is because Greene & Ho (2007) measured the luminosity and FWHM of the emission lines of the extracted AGN components but Shen et al. (2011) measured the lines without eliminating the host galaxy component from the spectra. However, this difference with ∼0.2 dex does not change our main conclusions as shown later. For the AGNs registered in both catalogs, we use the data in the catalog of Shen et al. (2011).

Zoom In Zoom Out Reset image size

We used the UKIDSS DR8 Large Area Survey (LAS) catalog (Lawrence et al. 2007) for galaxy samples. For each AGN, the UKIDSS K-band data are searched around the AGN coordinates within 1°. We selected AGN samples of z > 0.1 in order to analyze the projected number density as a function of the projected distance r_p within r_p ⩽ 7 Mpc in the following. To obtain the data around the sample AGNs, we repeat to search galaxy data in the UKIDSS LAS catalog to the number of the sample AGNs. We obtained the data by accessing the UKIDSS VO service through the JVO command line tools. By means of the VO tools, we can recurrently access large data archives easily. To remove stars from the UKIDSS samples, we selected data for which the merged class flag equals 1 (galaxies) or −3 (probable galaxies). We also removed data of poor quality which have post-processing error quality bit flags larger than 255.

The limiting magnitude m_limit of the UKIDSS galaxy samples is estimated for each AGN sample, and the result is plotted in Figure 2. As can be seen from the figure, m_limit distributes around 19.7–20.3. For each AGN, we also estimated the threshold magnitude, m_th, below which detection efficiency for galaxies can be regarded as 100%. Figure 3 shows the absolute magnitude corresponding to m_limit and m_th at AGN redshift. The detailed definition of the m_limit and m_th is described in the next section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To remove the data of poor sensitivity, we adapted selection criteria of ρ₀ > 10⁻⁴ and z < 1.0, where ρ₀ is the average number density of galaxies detectable at the AGN redshift. ρ₀ was calculated by integrating the luminosity function up to the absolute magnitude corresponding to the apparent limiting magnitude m_limit as described in the next section and Shirasaki et al. (2011). Figure 4 shows the calculated ρ₀ as a function of AGN redshift. In the same figure, the number densities of the complete sample (blue) are also shown, and they are calculated by integrating the luminosity function up to the absolute magnitude corresponding to m_th. The combined AGN catalog, which is based on catalogs of Shen et al. (2011) and Greene & Ho (2007), lists 32,806 AGNs at redshift between 0.1 and 1.0. A total of 11,335 objects are distributed in the survey area of UKIDSS LAS among them. In Figures 2–4, values for areas within 1° around these 11,335 AGNs are plotted. A total of 10,482 AGNs are selected by the criterion of ρ₀ > 10⁻⁴ Mpc⁻³ for m_limit.

Zoom In Zoom Out Reset image size

To reduce the effect of foreground clusters which are accidentally located near the sample AGNs on the sky, we rejected samples with anomalous distribution of surrounding galaxies by the method described in Shirasaki et al. (2011). To reduce the effect of accidental alignment of the foreground cluster, we calculate the clustering coefficient, B_QG, around each AGN, which was defined as ξ(r) = B_QGr^−γ (Barr et al. 2003), where ξ(r) is the cross-correlation function. B_QG is calculated as B_QG = (3 − γ)/(2πC(γ))(N_total − N_bg/ρ₀)(1 Mpc)^{γ − 3}, where N_total is the total number of observed galaxies at r_p < 1 Mpc, where r_p is a perpendicular distance, and N_bg is the expected background count at r_p < 1 Mpc (see Section 3 for the definition of ξ, r_p, γ, and C(γ)). We reject AGN samples with |B_QG| > 10⁴ (30 times the clustering coefficient of the Abell class 0 objects). In addition, we select sample AGNs without the effect from the nearby cluster located in regions offset from the AGNs. We count the number density n(r_p) of UKIDSS galaxies for each circular region with Δr_p = 0.2 Mpc width around AGNs and compute their statistical error. We adopt the following criteria: the reduced χ² of the radial number density relative to the flat distribution is χ²/(n − 1) ⩽ 3; the maximum deviation, σ_max, of the n(r_p) is smaller than 5σ. In Section 4, we also show the results without these selections for comparison.

We also present results of analysis for the AGN sub-sample whose deviations of surface number density of UKIDSS sources are smaller than 1.5σ at a limiting magnitude m_limit to check again the effect of foreground objects. Figure 5 shows the projected number density, n_7 Mpc, of galaxies in the area around the AGN with an angular radius corresponding to 7 Mpc (comoving) at the AGN redshit as a function of m_limit. Because the projected number density of the foreground or background galaxies is ∼30 times larger than the galaxies in the host clusters even for the highly clustered region with r₀ ∼ 10 Mpc, this criterion rarely rejects sample AGNs with the real overdensities around them, but rejects AGNs located near the foreground clusters.

Zoom In Zoom Out Reset image size

We used 9394 AGNs, which were selected by criteria described in Shirasaki et al. (2011), in the following analysis. A total of 8060 AGNs were selected by the 1.5σ criterion for n_7 Mpc. Figure 6 shows the distribution of the virial mass of the 9394 SMBHs as a function of redshift. A total of 1202 AGNs were extracted from the catalog of Greene & Ho (2007, red crosses), and the others were from the catalog of Shen et al. (2011, green circles). For the samples by Greene & Ho (2007), we plot the recomputed virial mass using the parameters of mass estimator in Shen et al. (2011). A total of 3749 objects among the samples from Shen et al. (2011) had been identified by the uniform criteria of SDSS quasars (Richards et al. 2002). We also present the results of the clustering analysis for these sub-samples. We divided our sample into four mass bins of log (M_BH/ M_☉) = 6.5–7.2, 7.2–8.0, 8.0–9.0, and 9.0–10.0, and three redshift bins of z = 0.1–0.3, 0.3–0.6, and 0.6–1.0, to see mass and redshift dependencies. The number, the averaged mass, and the averaged redshift of AGN samples for each mass range and redshift range are listed in Table 1.

Zoom In Zoom Out Reset image size

Table 1. Statistics of Fitting Parameters for Each Virial Mass and Redshift Group

Virial Mass	Redshift	nAGNa	〈log (MBH/ M☉)〉b	〈z〉c	r0d	nbge	〈ρ0〉f
log (MBH/ M☉)					( h−1 Mpc)	(Mpc−2)	(10−3Mpc−3)
All	0.1–1.0	9394	8.42	0.59	$5.8^{+0.8}_{-0.6}$	10.473 ± 0.005	1.9 ± 0.4
9.0–10.0	0.1–1.0	1331	9.23	0.72	$8.2^{+1.6}_{-1.1}$	4.542 ± 0.008	0.93 ± 0.25
	0.1–0.3	33	9.27	0.24	⋅⋅⋅  g	⋅⋅⋅	⋅⋅⋅
	0.3–0.6	347	9.24	0.48	$8.7^{+1.4}_{-1.0}$	7.29 ± 0.02	2.0 ± 0.5
	0.6–1.0	951	9.22	0.82	$7.0^{+4.9}_{-2.1}$	2.70 ± 0.01	0.38 ± 0.15
8.2–9.0	0.1–1.0	5119	8.60	0.66	$7.0^{+1.2}_{-0.8}$	6.001 ± 0.005	1.3 ± 0.3
	0.1–0.3	320	8.50	0.24	$7.3^{+0.9}_{-0.7}$	29.40 ± 0.04	5.3 ± 0.9
	0.3–0.6	1635	8.56	0.47	$6.5^{+1.0}_{-0.7}$	7.49 ± 0.01	2.1 ± 0.5
	0.6–1.0	3164	8.62	0.80	$7.6^{+2.3}_{-1.3}$	2.87 ± 0.004	0.42 ± 0.16
7.5–8.2	0.1–1.0	2278	7.91	0.44	$4.4^{+0.6}_{-0.5}$	15.69 ± 0.01	2.8 ± 0.6
	0.1–0.3	664	7.82	0.22	$4.5^{+0.6}_{-0.5}$	38.26 ± 0.03	5.7 ± 1.0
	0.3–0.6	967	7.92	0.44	$4.3^{+0.7}_{-0.6}$	8.53 ± 0.01	2.3 ± 0.5
	0.6–1.0	642	8.00	0.75	$2.5^{+1.8}_{-1.7}$	3.22 ± 0.01	0.53 ± 0.19
6.5–7.5	0.1–1.0	635	7.22	0.25	$4.6^{+0.6}_{-0.5}$	38.65 ± 0.03	5.4 ± 0.9
	0.1–0.3	513	7.18	0.20	$4.8^{+0.6}_{-0.5}$	45.69 ± 0.04	6.1 ± 1.0
	0.3–0.6	101	7.36	0.38	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
	0.6–1.0	21	7.30	0.77	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅

Notes. ^aNumber of sample AGNs. ^bAverage of logarithm of BH mass. ^cAverage redshift. ^dCorrelation length; the error contains the systematic error due to the uncertainty of ρ₀ and the 1σ statistical error. ^eAverage of the projected number density of background galaxies. ^fAverage of the averaged number density of galaxies at the AGN redshift. ^gWe do not derive parameters for the sub-sample with n_AGN < 200.

Download table as:  ASCIITypeset image

For ∼6% of AGNs in our sample, radio counterparts were found in the Faint Images of the Radio Sky at Twenty cm source catalog (Schneider et al. 2010). The percentage of objects with radio counterparts is higher for AGNs with higher BH mass. The fraction of AGNs with a radio counterpart is 14% at M_BH > 10⁹ M_☉. For ∼12% of the sample AGNs, X-ray counterparts were found in the ROSAT catalog. The fraction of AGNs with X-ray detection is almost the same for all the BH mass ranges.

We have followed the method described in Shirasaki et al. (2011) for the cross-correlation analysis between AGNs and galaxies.

The cross-correlation function of AGNs and galaxies ξ(r) can be expressed as an excess in the number density of galaxies ρ(r) relative to the average density ρ₀ at the AGN redshift

$$\begin{equation} \xi (r) = \frac{\rho (r)}{\rho _{0}} - 1, \end{equation} \tag{ 1 }$$

where r represents the distance from an AGN.

In this analysis, the redshift of galaxies is not measured, thus a projected cross-correlation function ω(r_p) is calculated from projected number densities of galaxies n(r_p):

$$\begin{equation} \omega (r_{p}) = \frac{n(r_{p}) - n_{\rm bg}}{\rho _{0}}, \end{equation} \tag{ 2 }$$

where r_p represents the projected distance from an AGN at the redshift and n_bg represents the surface density of background/foreground galaxies.

3.1. Estimating the Average Density and the Limiting Magnitude

Estimation of ρ₀ is crucial for this analysis. However, we cannot estimate ρ₀ directly from the data themselves since the information on the redshift is lacking. In this work, we estimated ρ₀ based on the luminosity function of galaxies obtained by previous studies (Cirasuolo et al. 2007; Gabasch et al. 2004, 2006; Kochanek et al. 2001; Montero-Dorta & Prada 2009). The luminosity function ϕ(M; z, λ) is parameterized as a function of redshift z and rest-frame wavelength λ as follows.

We use the Schechter function to represent the luminosity function,

$$\begin{eqnarray} \phi (M) =& 0.4 \cdot \phi _{*} \cdot \ln (10) \cdot 10^{-0.4 (M - M_{*}) (\alpha + 1)}\nonumber\\& \cdot \exp {[-10^{-0.4 (M - M_{*})}]}. \end{eqnarray} \tag{ 3 }$$

In Figure 7, M_* and ϕ(M₀) derived from the fitting function obtained in the works by Cirasuolo et al. (2007), Gabasch et al. (2004), Gabasch et al. (2006), Kochanek et al. (2001), and Montero-Dorta & Prada (2009) are plotted as a function of redshift for each rest-frame wavelength. M₀ is a reference magnitude where the luminosity function is normalized and parameterized as a function of redshift, and it is selected at a dimmer side of the M_* magnitude. Those data points are fitted with a third-degree polynomial function of redshift as shown in Figure 7 as solid lines. The standard deviations of M_* and ϕ(M₀) from the fitted functions are 0.2 mag and 15%, respectively. For the parameter α, we used fixed values such as −1.1 for λ < 400 nm, −1.25 for 400 ⩽ λ < 1000 nm, and −1.0 for λ ⩾ 1000 nm.

Zoom In Zoom Out Reset image size

As M_* and ϕ(M₀) at an arbitrary redshift for the eight wavelength bands can be calculated using the parameterization derived above, these parameters at an arbitrary wavelength are derived by interpolating them as a function of wavelength with a cubic spline method. In this way, we parameterized the luminosity function as a function of redshift and rest-frame wavelength. We estimate ρ₀ by means of the Schechter function with the parameters at the AGN redshift and wavelength corresponding to the observed K band.

ρ₀ can be calculated by integrating ϕ(M; z, λ) to the absolute magnitude M at an AGN redshift corresponding to an apparent limiting magnitude m_limit,

$$\begin{equation} \rho _0 = \int ^{m_{\rm limit}-{\rm DM}(z)}_{m_{\rm low}-{\rm DM}(z)} \phi (M;z,\lambda) dM, \end{equation} \tag{ 4 }$$

where DM is the distance modulus and m_low is a lower boundary of the apparent magnitude.

As the limiting magnitude varies among the AGN samples, it was estimated from the measured magnitude distribution N(m) as explained below. The observed magnitude distribution N_obs(m) can be expressed as a multiplication of the true magnitude distribution N_true(m) and the detection efficiency DE(m):

$$\begin{equation} N_{\rm obs}(m) = N_{\rm true}(m) \times {\rm DE}(m). \end{equation} \tag{ 5 }$$

We model N_true(m) and DE(m) as introduced in Shirasaki et al. (2011):

$$\begin{eqnarray} N_{\rm true}(m) = \left\lbrace \begin{array}{@{}ll}c \cdot 10^{a (m-m_{\rm b})} &(m < m_{\rm b}) \\ c \cdot 10^{b (m-m_{\rm b})} &(m \ge m_{\rm b}), \\ \end{array} \right. \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} {\rm DE}(m) = \left\lbrace \begin{array}{@{}ll}1 &(m < m_{\rm th}) \\ \exp {\left(-(m-m_{\rm th})^2/\sigma _{m}^{2}\right) } &(m \ge m_{\rm th}). \\ \end{array} \right. \end{eqnarray} \tag{ 7 }$$

By fitting the model function of N_obs(m) to the observed magnitude distribution, we obtained the model parameters a, b, c, m_b, m_th, and σ_m for each area around an AGN sample. We determine ρ₀ as

$$\begin{equation} \rho _{0} = \int ^{\infty }_{m_{\rm low} - {\rm DM}} \phi (M; z, \lambda) {\rm DE}(M+{\rm DM}) dM. \end{equation} \tag{ 8 }$$

We derive m_limit by equating the right-hand sides of Equations (4) and (8).

The uncertainty of ρ₀ determined as explained above is dominated with the uncertainties of the model parameters b and uncertainties of parameterization of the luminosity function, and they are taken into account as a systematic error in estimating the cross-correlation length. We estimate the uncertainty of b as 0.04, which comes from the standard deviation of b calculated for each AGN sample. Considering the uncertainty of b, corresponding uncertainties of m_th and σ_m are estimated by comparing these fitting parameters obtained by fixing the b parameter to b_best  ±  0.04, where b_best is the best-fitting parameter for the AGN sample. The uncertainty of ρ₀ originated from the uncertainties of m_th and σ_m (i.e., b) are calculated with error propagation. The uncertainty $\sigma _{\rho _{0},b}$ is plotted as a function of redshift in Figure 8.

Zoom In Zoom Out Reset image size

To evaluate the uncertainty of ρ₀ due to the uncertainties of parameterization of the luminosity function, we assumed uncertainties of M_* and ϕ(M₀) to be 0.2 mag and 15%, respectively, which are the standard deviation of data points in Figure 7 from the fitting functions.

$\sigma _{\rho _{0},M_{*}}$ is the uncertainty of ρ₀ due to the uncertainty of M_*, and is plotted in Figure 8. The uncertainty of ρ₀ due to the uncertainty of ϕ(M₀) is independent from redshift and is a constant value of $\sigma _{\rho _{0},\phi (M_{0})}=0.15$. Then the total uncertainty of ρ₀ is calculated as

$$\begin{equation} \sigma _{\rho _{0}}^{2} = \sigma _{\rho _{0},b}^{2} + \sigma _{\rho _{0},M_{*}}^{2} + \sigma _{\rho _{0},\phi (M_{0})}^{2}. \end{equation} \tag{ 9 }$$

3.2. Cross-correlation

We assume the power-law form for the cross-correlation function,

$$\begin{equation} \xi (r) = \left(\frac{r}{r_{0}} \right) ^{-\gamma }, \end{equation} \tag{ 10 }$$

where r₀ is a correlation length and γ is a power-law index and fixed to 1.8, which is a canonical value measured by many other works. Then the projected cross-correlation function can be expressed as

$$\begin{eqnarray} \omega (r_p) &=& 2 \int ^\infty _0 \xi (r_p,\pi) d\pi = 2 \int ^\infty _{r_p} \frac{r \xi (r)}{\sqrt{r^2-r_p^2}} dr \nonumber\\ &=& r_p \left(\frac{r_0}{r_p} \right) ^\gamma \frac{\Gamma \left(\frac{1}{2}\right)\Gamma \left(\frac{\gamma -1}{2}\right)}{\Gamma \left(\frac{\gamma }{2}\right)}, \end{eqnarray} \tag{ 11 }$$

where π and r_p are distance along and perpendicular to the line of sight, respectively, and Γ is the gamma function.

From Equations (2) and (11), the projected number density of galaxies around an AGN can be modeled as

$$\begin{equation} n(r_{p}) = C(\gamma) \times \rho _{0} \times r_{p} \left(\frac{r_{0}}{r_{p}} \right)^{\gamma } + n_{\rm bg}, \end{equation} \tag{ 12 }$$

where the term of the gamma function is represented by C(γ). By fitting this model function to the observed projected number density, we can derive the model parameters r₀ and n_bg. ρ₀ is determined by the method described in Section 3.1. Since the clustering signature for each AGN is too weak to derive the parameters individually, we applied this fitting to the averages of n(r_p) and ρ₀ for a given AGN group,

$$\begin{equation} \langle n(r_{p})\rangle = C(\gamma) \times \langle \rho _{0} \rangle \times r_{p} \left(\frac{r_{0}}{r_{p}} \right)^{\gamma } + \langle n_{\rm bg} \rangle, \end{equation} \tag{ 13 }$$

and derive r₀ and 〈n_bg〉. The uncertainty of r₀ is calculated as the square root of the square sum of the systematic error derived from the uncertainty of ρ₀ described in Equation (9) and statistical error of 1σ by fitting 〈n(r_p)〉. It should be noted that the cross-correlation function obtained by this method is not an average for the AGN group, but an average weighted with ρ₀. Thus, the result is biased to the low-z and high-sensitivity samples.

By the analysis described above, we have estimated the scale length of the AGN–galaxy cross-correlation for the whole sample to be $r_0= 5.8^{+0.8}_{-0.6} \,h^{-1}\,{\rm Mpc}$. This is comparable with or slightly smaller than the results of the previous studies of AGN–galaxy cross-correlation (r₀ = 5.95 ± 0.90 h⁻¹ Mpc for X-ray AGNs at z = 0.7–1.4, Coil et al. 2009; 6.98 ± 0.6 h⁻¹ Mpc for optical AGNs at z < 1, Mountrichas et al. 2009; 6.0 ± 0.5 h⁻¹ Mpc for optical AGNs at 0.2 < z < 0.6, Padmanabhan et al. 2009). Donoso et al. (2010) found r₀ = 8.35 ± 0.09 for radio AGNs and r₀ = 5.02  ±  0.24 h⁻¹ Mpc for optical AGNs at z = 0.4–0.8. Krumpe et al. (2012) derived that $r_0=6.91^{+0.17}_{-0.18}\,h^{-1}\,{\rm Mpc}$ at z = 0.16–0.36 between SDSS AGNs and the SDSS main-galaxy sample and $r_0=7.21^{+0.21}_{-0.22}\,h^{-1}\,{\rm Mpc}$ at z = 0.36–0.5 for SDSS AGNs and luminous red galaxies. Mountrichas et al. (2013) derived that $r_0=4.4^{+0.1}_{-0.2}$ for X-ray selected AGNs at z ∼ 0.1.

Now, we present the results of the cross-correlation analysis adopted for the four mass ranges described in Figure 6 and Table 1, to see the dependence of the clustering amplitude on BH mass. Figures 9 and 10 show the measured projected number density n(r_p) and projected cross-correlation function ω(r_p), respectively, for each mass range. Four panels represent results for the mass ranges of log (M_BH/ M_☉) = 9.0–10.0 (top left), 8.2–9.0 (top right), 7.5–8.2 (bottom left), and 6.5–7.5 (bottom right). The projected number density of the areas of each circular ring is plotted with the Poisson error bars. The solid lines denote least-χ² fitting by the power law (Equation (13)). The fitting parameters are summarized in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 11 shows the scale length r₀ of the cross-correlation as a function of the virial mass of SMBHs. Error bars of r₀ denote the square root of the square sum of the systematic error calculated from $\sigma _{\rho _0}$ described in Section 3.1 and the statistical error of 1σ. We can see a trend that r₀ increases as mass increases for M_BH > 10⁸ M_☉. When we consider only the statistical error, the significance of difference of r₀ is 9.6σ between log (M_BH/ M_☉) = 7.5–8.2 bin and 8.2–9.0 bin, and 3.4σ between 8.2–9.0 bin and 9.0–10.0 bin. For M_BH ≲ 10⁸ M_☉, we cannot see the significant mass dependence. These trends are also seen for the data set grouped with finer mass ranges.

Zoom In Zoom Out Reset image size

Figure 12 shows r₀ for AGN samples from the catalog of Shen et al. (2011, green circles), which are identified by the uniform criteria of SDSS (Richards et al. 2002) and for samples of Greene & Ho (2007, blue triangles). As seen in the figure, the mass dependencies of both sub-samples are quite similar. We can see, however, a little offset between the results for the two sub-samples. This can be due to the systematic difference of the BH mass estimate with ∼0.2 dex between the two catalogs as mentioned in Section 2.

Zoom In Zoom Out Reset image size

We rejected sample AGNs with anomalous galaxy distribution around the AGNs, as described in Section 2 to remove the effect from foreground objects. To see whether we are removing real clustering by this rejection, we present the results of the analysis done both with and without the selection, in Figure 13. The blue triangles show the results without any sample rejection. We cannot find a significant clustering signal at log (M_BH/ M_☉) < 7.5 for this analysis. The red squares are the fiducial results with the selection criteria described in Shirasaki et al. (2011). For the green circles, we also reject AGN samples with an extremely high or low projected number density (n_7 Mpc) of galaxies around them.

Zoom In Zoom Out Reset image size

The analysis without selection gives slightly smaller r₀. This indicates that the clustering feature is weakened by foreground contamination for the analysis without sample rejection. The analysis with the n_7 Mpc criterion gives $r_0=6.3^{+0.9}_{-0.6}\,h^{-1}\,{\rm Mpc}$ for the whole mass range. This is ∼0.5 h⁻¹ Mpc larger than the result without the n_7 Mpc criterion. On the other hand, we see almost the same relative mass dependence for all three cases.

4.1. AGN Redshift

As seen in Figure 6, distribution of the BH mass in our sample depends on the redshift. This is because SMBHs with larger mass tend to be luminous and can be observed even at higher redshift. To see the dependence on the redshift, we divided the samples into three groups with redshift ranges of z = 0.1–0.3, z = 0.3–0.6, and z = 0.6–1.0. Figure 14 shows the estimated r₀ for the three redshift ranges and four mass ranges. As seen in the top panel of this figure, r₀ is not dependent on the redshift. Therefore, the increasing trend of r₀ seen in Figure 11 would not be due to the redshift bias. We could also see the mass dependence for the sub-samples of low redshift (z = 0.1–0.3, red squares) and intermediate redshift (z = 0.3–0.6, green circles). For the higher redshift sample (z = 0.6–1.0), the estimated error is too large to see the dependence on BH mass. We summarize the estimated r₀, 〈ρ₀〉, and 〈n_bg〉 for each mass range and redshift range in Table 1.

Zoom In Zoom Out Reset image size

To remove the possible redshift bias, we also present the mass dependence for sub-samples with the normalized redshift distributions. We constructed sub-samples as follows. For redshift bins with Δz = 0.1, the selection probability is determined to give the same redshift distribution for the four mass ranges. We randomly selected AGN samples following the selection probability and constructed a set of sub-samples. Figure 15 shows one example of the sub-sample. We constructed 10 sets of sub-samples and measured r₀ for them. For these sub-samples, there is no correlation between redshift and BH mass.

Zoom In Zoom Out Reset image size

Figure 16 shows the mass dependence of r₀ for the resampled AGNs with normalized redshift distribution. We plot medians of the 10 sets of sub-samples as green circles. The error bars show maximum and minimum values for the 10 sets. The resampled AGNs also show similar mass dependence as Figure 11. Large-mass SMBHs show strong clustering amplitude at M_BH > 10⁸ M_☉. The clustering amplitude for AGNs with log (M_BH/ M_☉) = 6.5–7.5 and 7.5–8.2 is almost the same although the error bars are very large. The relative mass dependence for the sub-samples is free from redshift bias since the sub-samples for the four mass ranges have the same redshift distribution. We note that there is also no correlation between BH mass and luminosity for these sub-samples. Therefore, the BH mass dependence is thought to be due neither to redshift bias nor to luminosity bias.

Zoom In Zoom Out Reset image size

4.2. AGN Luminosity

We divided our sample into three luminosity ranges of L₅₁₀₀ < 10^44.5 erg s⁻¹, 10^44.5 erg s⁻¹ ⩽L₅₁₀₀ < 10^44.8 erg s⁻¹, and 10^44.8 erg s⁻¹ ⩽L₅₁₀₀ to see luminosity dependence, where L₅₁₀₀ is the monochromatic continuum luminosity at rest-frame 5100 Å. The top panel of Figure 17 shows the dependence of r₀ on luminosity. We cannot find significant luminosity dependence for all four mass ranges. On the other hand, we could see the mass dependence at M_BH > 10⁸ M_☉ for all three luminosity ranges, as seen in the bottom panel. Therefore, we can conclude that the BH mass dependence seen in Figure 11 is not due to the dependence on luminosity.

Zoom In Zoom Out Reset image size

Shen et al. (2009) argued that the amplitude of AGN–AGN auto-correlation depends weakly on optical luminosity. Donoso et al. (2010) found that the clustering amplitude varies with radio luminosity on scales less than ∼1 Mpc but is almost independent of luminosity for the larger scale.

These results indicate that the clustering amplitude on a large scale depends on the mass of SMBHs but weakly depends on luminosity.

4.3. Completeness and Luminosity of the Galaxy Sample

As described in Section 3.1, we correct the incompleteness of the faint end of the galaxy sample by estimating the detection efficiency DE(m) based on the magnitude distribution of the UKIDSS sample for each area around an AGN. There may be a criticism that the correction of the incompleteness of the galaxy sample is somehow biased to the luminosity of the galaxy sample. Figure 18 shows the cross-correlation length calculated for a complete galaxy sample that consists of bright galaxies with m < m_th, where m_th is a threshold magnitude, below which DE(m) = 1 (see Equation (7)). For this analysis, ρ₀ is also recomputed by integrating the luminosity function up to m_th. This result also shows that the cross-correlation increases above M_BH ∼ 10⁸ M_☉. Therefore, the increasing trend is not due to the ambiguity that comes from using an incomplete galaxy sample.

Zoom In Zoom Out Reset image size

It is known that brighter galaxies tend to cluster more strongly than dimmer galaxies. Thus, it might be possible to explain the larger cross-correlation length for more massive SMBHs by the bias due to the galaxy brightness. Figure 19 shows the cross-correlation length calculated for luminosity-limited samples. We selected the luminosity-limited galaxy samples which are defined as M ≡ m − DM(z_AGN) < −22 and M < −23.5 for each AGN, where DM(z_AGN) represents distance modulus for the AGN redshift. M is not the absolute magnitude for foreground and background galaxies, but they should make no contribution to the clustering signal and not affect r₀. For these analyses, we selected AGN samples with M_th ≡ m_th − DM > −22.0 and M_th > −23.5, respectively (see also Figure 3). Therefore, these absolute-magnitude-limited galaxy samples are also "complete"(m < m_th). Although the error bar is relatively large, we can see the similar trend of the cross-correlation length against virial mass of SMBHs as Figure 11. The significance of the difference of r₀ by considering statistical error is 1.9σ for mass ranges between log (M_BH/ M_☉) = 7.5–8.2 and 8.2–9.0, and 1.7σ for between 8.2–9.0 and 9.0–10.0, for the analysis with the M < −23.5 sample, and 3.8σ for between log (M_BH/ M_☉) = 7.5–8.2 and 8.2–9.0 for the M < −22.0 sample. Therefore, we can conclude that the increase of the cross-correlation length seen in Figure 11 is not only due to the bias related with the galaxy brightness. The estimated scale length is larger than Figure 11 since brighter galaxies are more strongly clustered.

Zoom In Zoom Out Reset image size

In the previous studies, Shen et al. (2009) have shown that most massive SMBHs are more strongly clustered than the remainders from auto-correlation analysis of quasars. It has also been shown that radio selected AGNs are more strongly clustered than the cases for the optically selected AGNs (Hickox et al. 2009; Donoso et al. 2010), and the characteristic BH mass of radio AGNs is higher than optically selected ones. For 2300 quasars at z = 0.6–1.2, Zhang et al. (2013) found that the clustering amplitude is larger for quasars with more massive BHs at a significance of ∼1σ, using photometric data from SDSS Stripe 82. Our results are consistent with these previous studies. These may indicate that the environments of galaxies have played an important role for the growth of high-mass SMBHs. The clustering amplitude is relevant to the mass of the host dark-matter halo and the frequency of major merger. If the mass growth of SMBHs is mainly driven by the major mergers of galaxies, massive SMBHs are expected to be in massive halos.

In contrast, we did not find a significant luminosity dependence. Luminosity is thought to represent gas accretion activity at this time. The activity of SMBHs is thought to be a transient event and not strongly correlate with large-scale structure. On the other hand, BH mass is thought to represent cumulative accretion history and merger history of BHs, and can be related with large-scale environment.

For less massive BHs with M_BH ≲ 10⁸ M_☉, the significant correlation between r₀ and M_BH is not seen in our study. BH mass has been thought to be correlated with the mass of dark-matter halos (e.g., Ferrarese 2002). Recently, however, Kormendy & Bender (2011) found that the mass of SMBHs does not directly correlate with dark-matter halo mass, at least for low-mass SMBHs in disk galaxies, based on the observations of nearby SMBHs for which the mass of the host dark halos is derived by the stellar kinematics. Graham (2012) and Scott et al. (2013) examined the relation between M_BH and mass of the host spheroid, M_sph, and found that the slope of the log M_BH − log M_sph relation changes at M_BH = 2 × 10⁸ M_☉.

One possible scenario to explain the absence of the positive mass dependence for M_BH ≲ 10⁸ M_☉ would be that the less massive BHs could be formed in the isolated galaxies by secular processes. If they have grown by secular processes, then the mass of a seed BH should be much larger than a typical stellar mass BH (for a review, see, e.g., Volonteri 2010). Some authors (Lodato & Natarajan 2006; Begelman et al. 2006) argue that SMBHs with 10⁴–10⁶ M_☉ are formed through the direct collapse of pre-galactic gas at z > 10. Such a heavy seed BH can grow to ∼10⁷–10⁸ M_☉ by a few Gyr without a BH merger under the assumption that the mass accretion rate is ∼0.1 times the Eddington rate. Another scenario would be that they are in a growing phase by the major mergers of galaxies.

For the clustering amplitude of less massive SMBHs, the contribution of the AGNs hosted in satellite galaxies in the massive dark-matter halos can also be important, while most massive SMBHs are in the central regions of dark-matter halos. Padmanabhan et al. (2009) argue that the satellite fraction of quasars is more than 25% in their sample, which is selected from the SDSS catalog. The percentage should be larger for the less massive BHs. Further investigation of less massive BHs is crucial to understand the formation and evolution mechanisms of an SMBH.

We have investigated the clustering of galaxies around 9394 AGNs for z = 0.1–1. We obtained the galaxy data of UKIDSS LAS by means of VO tools. Our results are free from the effect of cosmic variance owing to the large sample covering the large area of the sky. The estimated correlation length ranges between 4 and 10 h⁻¹ Mpc depending on BH mass, and depends neither on redshift for z = 0.1–1 nor on luminosity. The results may indicate that higher mass BHs reside in a more clustered environment for M_BH > 10⁸ M_☉.

While our results show a positive mass dependence for M_BH > 10⁸ M_☉, our results for M_BH ≲ 10⁸ M_☉ show no significant mass dependence. Although our sample of less massive BHs is small, this would give a critical mass scale for the emergence of an environment effect for BH growth.

In this study, the redshift range where the BH mass dependence of the cross-correlation is measured with good accuracy is limited to below z ∼ 0.6. This is because the number density of UKIDSS LAS galaxies is too low at higher redshift, and also the systematic error due to the uncertainty of the M_* parameter of the luminosity function is too large. To understand the relation of the BH mass accretion history and its environment, it is crucial to observe its evolution up to at least redshift two where the number density of QSOs is maximal and the mass accretion is expected to be the most prosperous. To reduce the uncertainty, it is crucial to perform deeper observations so that the limiting magnitude reaches well beyond the characteristic luminosity of galaxies at AGN redshifts. At redshifts larger than 0.5, the dominant factor to the uncertainty of ρ₀ is an uncertainty of M_* parameterization and the ρ₀ uncertainty becomes larger than 20%. Observations two magnitudes deeper will extend the redshift range where the uncertainty of ρ₀ is less than 20% up to 1.0. Future instruments such as Hyper Suprime-Cam (HSC) can measure the AGN environment more accurately with good statistics. When the survey is performed with 26 mag in the r band with HSC, we can estimate ρ₀ with accuracy less than 20% up to redshift 2, and as a result can estimate r₀ with an accuracy less than 10%. Such a deep and wide survey would reveal the mechanism of AGN evolution at an important epoch at which its activity was the most prosperous.

We appreciate S. Eguchi and M. Enoki for their useful discussions. Results are based on data obtained from the Japanese Virtual Observatory, which is operated by the Astronomy Data Center, National Astronomical Observatory of Japan. This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France. This work is based on data obtained as part of the UKIRT Infrared Deep Survey. Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web site is http://www.sdss.org/.