Clustering of galaxies around AGN in the HSC Wide survey Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.

# Clustering of galaxies around AGN in the HSC Wide survey ††thanks: Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.

###### Abstract

We have measured the clustering of galaxies around active galactic nuclei (AGN) for which single-epoch virial masses of the super-massive black hole (SMBH) are available to investigate the relation between the large scale environment of AGNs and the evolution of SMBHs. The AGN samples used in this work were derived from the Sloan Digital Sky Survey (SDSS) observations and the galaxy samples were from 240 deg S15b data of the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP). The investigated redshift range is 0.6–3.0, and the masses of the SMBHs lie in the range  M. The absolute magnitude of the galaxy samples reaches to 18 at rest frame wavelength 310 nm for the low-redshift end of the samples. More than 70% of the galaxies in the analysis are blue. We found a significant dependence of the cross-correlation length on redshift, which primarily reflects the brightness dependence of the galaxy clustering. At the lowest redshifts the cross-correlation length increases from 7  Mpc around  mag to 10  Mpc beyond  mag. No significant dependence of the cross-correlation length on BH mass was found for whole galaxy samples dominated by blue galaxies, while there was an indication of BH mass dependence in the cross-correlation with red galaxies. These results provides a picture of the environment of AGNs studied in this paper being enriched with blue starforming galaxies, and a fraction of the galaxies are being evolved to red galaxies along with the evolution of SMBHs in that system.

galaxies: active — large-scale structure of universe — quasars: general

## 1 Introduction

Most galaxies have a supermassive black hole (SMBH) with mass greater than 10 M at their center (Richstone et al., 1998). Observations of the local Universe have revealed that the mass of the SMBH correlates with the several properties of the bulge component of the host galaxy (Magorrian et al., 1998; Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Ho, 2007). This observational evidence suggests that a SMBH and its host galaxy co-evolve in a coordinated way in spite of the nine orders of difference in their physical size scale (Kormendy & Ho, 2013).

SMBHs grow through the accretion of gas from their host galaxies or large-scale environment. Accretion in a secular mode, which arises through internal dynamical processes such as bar or disk instability or external processes driven by galaxy interaction, is one of the mechanisms to deliver the gas into the SMBH (Kormendy & Kennicutt, 2004). As the mass accretion rate of the secular mode cannot be as high as to maintain the activity seen in bright QSOs (Menci et al., 2014), this mode could operate in lower-luminosity AGNs.

A merger of galaxies can induce gravitational torques which drive inflows of cold gas toward the center of galaxies (e.g., Hopkins et al., 2008). Observational evidence has been obtained that shows the relation between AGN activity and a galaxy merger (e.g., Sanders et al., 1988; Treister et al., 2012). The activity of the bright QSOs could be explained by this model. Although the accretion can be the most efficient in this mode, it is also expected that AGN feedback promptly operates and stops the gas inflows when the mass of the SMBH becomes as large as 10 M (Fanidakis et al., 2013).

As an alternative process to make a SMBH evolve above 10 M, quiescent gas accretion from the hot halo (Keres et al., 2009; Fanidakis et al., 2013) and/or recycled gas from evolving stars (Ciotti & Ostriker, 2001, 2007) has been proposed. According to the predictions of semi-analytical modeling by Fanidakis et al. (2013), AGNs fueled in the hot-halo mode are located in more massive dark matter haloes than the AGNs fueled in the merger-driven starburst mode.

The mass of the host dark matter halo that AGNs reside can be inferred from the auto-correlation function of AGNs at distance scale 1 Mpc. According to the analysis of auto-correlation of SDSS QSOs by Ross et al. (2009), the halo mass is almost constant at 2 in the redshift range from 0.3 to 2.2. It is also possible to estimate the halo mass from the cross-correlation between AGNs and galaxies if the auto-correlation functions of galaxies are precisely obtained. According to the cross-correlation studies the dark matter halo mass estimated to be depending on the detection wavelength of the AGNs (e.g. Hickox et al., 2009; Krumpe et al., 2012). Studies on clustering and/or environments of AGNs have also been reported elsewhere (e.g. Croom et al., 2005; Coil et al., 2009; Silverman et al., 2009; Donoso et al., 2010; Allevato et al., 2011; Bradshaw et al., 2011; Shen et al., 2013; Mountrichas et al., 2013; Zhang et al., 2013; Georgakakis et al., 2014; Krumpe et al., 2015; Ikeda et al., 2015).

In spite of a large number of AGN clustering studies published so far, the studies focused on the relation between the mass of the central BH, which is the most fundamental property on which the history of mass accretion is imprinted, and also the properties of galaxies such as color and luminosity function around it are very limited. Since the interaction in the scale of group or cluster of galaxies can induce concurrent activities of starformation, mass accretion to the central BH, and transition to the red sequence in the constituent galaxies, the contribution of such large-scale phenomena to the evolution of SMBH can be inferred from the surrounding galaxies.

Shen et al. (2009) measured the clustering of QSOs for samples divided by their BH mass at redshifts 0.4–2.5 and didn’t find significant dependence except for the most massive sample for which marginally larger clustering was found in 2 level. Komiya et al. (2013) examined dependence of the AGN-galaxy clustering on the BH mass at redshift range from 0.3 to 1.0 using the UKIDSS data for the galaxy samples, and found that the cross-correlation length increases above . Krolewski & Eisenstein (2015) also examined the BH mass dependence of the AGN-galaxy clustering at redshift 0.8 using SDSS and WISE data for the galaxy samples, and found no significant relationship between clustering amplitude and BH mass. Krumpe et al. (2015) measured clustering of soft X-ray and optically selected AGNs at redshifts from 0.16 to 0.36 and detected a weak dependence on the mass at significance level of 2.7 , while they didn’t detect significant dependence on the Eddington ratio. They conclude that the mass dependence is the origin of the observed weak X-ray luminosity clustering dependence.

Extending the dataset of Komiya et al. (2013), Shirasaki et al. (2016) derived, for the first time as derived from statistically significant number of samples, color and absolute magnitude distributions of galaxies around AGNs, and found that the increase of the cross-correlation length found by Komiya et al. (2013) is due to the increase in the number density of red galaxies. Those results indicate that the most massive SMBHs are evolved in dark matter haloes more massive than the lower mass SMBHs, and the surrounding galaxies also evolve in a coordinated way with the SMBHs which are mostly located in the center of host halo. To extend further the study given in Shirasaki et al. (2016) up to redshifts 3, which covers the era of the peak of star formation and mass accretion rate to SMBHs, we require a deeper wide multi-band survey.

The Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) is a multi-band imaging survey conducted with the HSC on the 8.2m Subaru Telescope. The survey consists of three layers: Wide (1400 deg, r 26), Deep (27 deg, r27), and UltraDeep (3.5 deg, r28). The HSC-SSP Wide survey provides the first opportunity to investigate the environment of AGNs up to redshift three with unprecedented statistics. Using the unique dataset of HSC-SSP Wide survey, this paper’s purpose is to measure not only the clustering of galaxies around AGNs but also their color and luminosity distribution as a function of SMBH mass at five redshift groups from 0.6 to 3.0. The results obtained in this work would give us a unified picture of evolution of galaxies and SMBHs under the large-scale structure of the Universe at their most important stage.

Throughout this paper, we assume a cosmology with , , and . All magnitudes are given in the AB system. All the distances are measured in comoving coordinates. The correlation length is presented in unit of Mpc. The group and parameter names that are frequently refereed to in the text are summarized in Table Clustering of galaxies around AGN in the HSC Wide survey thanks: Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. of the Appendix.

## 2 Datasets

### 2.1 AGNs

The AGN samples used in this paper were drawn from the QSO properties catalog of Shen et al. (2011) (S11) and SDSS DR12 Quasar catalog (DR12Q) of Paris et al. (2017). We used S11 catalog as a reference for the black hole (BH) mass estimate to be consistent with the previous studies (Komiya et al., 2013; Shirasaki et al., 2016), and the BH masses derived from the DR12Q catalog and spectral measurements were calibrated with those derived in S11.

We divided the AGN samples according to their measured redshift and BH mass. The mass used in this analysis is based on the H, Mg II, or C IV line width as given in the S11 and DR12Q catalogs. For AGNs drawn from the DR12Q catalog, the mass was calculated according to the formula derived by Mejia-Restrepo et al. (2016) as follows:

 MBH,MgII=8.05×106 M\solar(Lλ3001044 erg s−1)0.609(FWHM(Mg II)km s−1)2, (1)
 MBH,CIV=5.71×105 M\solar(Lλ1451044 erg s−1)0.57(FWHM(C IV)km s−1)2(Lp(C III])Lp(C IV))−2.09, (2)

where and are the continuum luminosity at 300 nm and 145 nm respectively, FWHM(Mg II) and FWHM(C IV) are the full width at half-maximum of the Mg II and C IV emission lines respectively, and and are peak luminosities of the corresponding lines. We measured the continuum luminosities, and , and the luminosity ratio, directly from the SDSS spectra, while the FWHMs of the emission lines were taken from the DR12Q catalog.

To check the consistency between the BH masses derived by different relations, we compared them for the same objects. The left panel in Figure 1 shows the comparison of black hole masses given in S11, M(S11), and those calculated for the same object in DR12Q, M(DR12Q), using equation (1). We found there is a systematic offset of 0.24 dex and scatter of 0.22, so we have corrected for the offset to the mass calculated using equation (1).

The right panel in Figure 1 shows the comparison of black hole masses calculated for objects in DR12Q using equation (1), M(DR12Q), and equation (2), M(DR12Q). Systematic offset of 0.07 dex and scatter of 0.48 dex were found between them. Thus the offset of M(DR12Q) to M(S11) was estimated to be 0.17, and the corresponding correction was made to M(DR12Q). As mentioned in Mejia-Restrepo et al. (2016) and apparent from Figure 1, M has a larger uncertainty than the mass measured based on other emission lines, so we used M, whenever it was available, as a mass estimator in first priority. The AGN samples given the C IV BH mass are mostly at . The effect of the large uncertainty of the C IV BH mass is limited only to the results for .

Coatman et al. (2016) recently showed that an empirical correction to C IV BH masses based on the blueshift of C IV emission line reduces the uncertainty. To apply the correction requires an accurate measure of the AGN systemic redshift. Thus, we decided to use the relation of equation (2) which does not require precise spectroscopic measurements.

The systematic offset and scatter between and given in S11 catalog are 0.009 dex and 0.38 dex respectively as shown in Figure 10 of Shen et al. (2011). Thus, no correction to is required.

For AGNs appearing in both catalogs, we used the redshift and black hole mass of the S11 catalog. The redshift range of the AGNs was chosen to be 0.6–3.0, so that it overlaps with the highest redshift bin (=0.6–1.0) of the previous studies (Komiya et al., 2013; Shirasaki et al., 2016) to cross check between both the results and extend up to the limit of sensitivity of this analysis. The BH mass range was set to M. We selected 6166 AGNs which are in the redshift and BH mass range and are located within the footprint of HSC-SSP Wide survey.

### 2.2 Galaxies

The galaxy sample was collected from the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) S15b Wide survey dataset. HSC-SSP is a three-layered, multi-band ( plus 4 narrow-band filters) imaging survey with the Hyper Suprime-Cam (HSC) (Miyazaki et al., 2012, 2017) on the 8.2m Subaru Telescope. The total area and the depth of observation will be 1400 deg with (Wide layer), 27 deg with (Deep layer), and 3.5 deg with (UltraDeep layer).

We used the dataset derived from Wide layer. The observed locations and effective area are summarized in Table 2.2. The typical depths of the observation are 26.8, 26.4, 26.4, 25.5, 24.7 for , , , , and bands, respectively. The detail of survey itself is described in Aihara et al. (2017a), and the content of the S15b dataset is in Aihara et al. (2017b). The S15b dataset were analyzed through the HSC pipeline (version 4.0.1) developed by the HSC software team using codes from the Large Synoptic Survey Telescope (LSST) software pipeline (Ivezic et al., 2008; Axelrod et al., 2010). The photometric and astrometric calibrations are made based on data obtained from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) 1 imaging survey (Magnier et al., 2013; Schlafly et al., 2012; Tonry et al., 2012).

The photometric magnitude used in this work is a CModel magnitude. The galactic reddening was corrected according to the dust maps derived by Schlegel et al. (1998).

The analysis performed in this paper is based on two galaxy samples; one is the -band detected sample which is drawn from all sources selected by the criteria for -band as described below and measured to be brighter than 27 magnitude in -band regardless of the detection in the other four bands; the other is the four-band detected sample which is selected by enforcing the same criteria for -bands except for the magnitude cut adapted only to -band data, regardless of the detection in -band. Since the observations in -band are shallower than the others, the detection in -band was not required to avoid the bias to redder galaxies. For the cross-correlation analysis, we used the -band detected sample. When we measure the distribution of galaxy color and luminosity around AGNs, we used the four-band detected sample. The term ”detected” used here means that the source satisfies the criteria defined below and has a non-NULL magnitude.

The criteria used to select -band detected samples are:

   iflags_pixel_edge is not True
AND iflags_pixel_interpolated_any is not True
AND iflags_pixel_saturated_any is not True
AND iflags_pixel_cr_any is not True
AND icmodel_flux_flags is not True
AND icentroid_sdss_flags is not True
AND detect_is_tract_inner is True
AND detect_is_patch_inner is True
AND deblend_nchild = 0


where iflags_pixel_edge is true if the source is near the edge of the frame; iflags_pixel_interpolated_any is true if any pixels in the source footprint have been interpolated due to saturation or cosmic rays; iflags_pixel_saturated_any is true if any pixels are saturated; iflags_pixel_cr_any is true if any pixels are masked as cosmic rays; iflags_pixel_bad is true if any pixels are masked for non-functioning or in severely vignetted areas; icmodel_flux_flags is false if the cmodel measurement fails; icentroid_sdss_flags is false if the SDSS centroiding algorithm fails; detect_is_tract_inner and detect_is_patch_inner are true if the source is in the inner region of a tract and patch, which is used to select a source of primary detection; deblend_nchild is the number of children this source was deblended into. Similar flag checks were performed for the other bands.

    flags_pixel_bright_object_center is not True


to remove the galaxy samples in bright source masks only for those located at 2 Mpc from the AGN. The bright source masks implemented in this data release are over-conservative in the choice of radius, so it has a drawback that decreases the significance of the clustering. Thus we did not adapt the bright source mask to the galaxy samples located at 2 Mpc from the AGN.

We found that there is a non-negligible number of false detections in the deblended sources, especially at larger magnitudes. To remove the false detections, we selected only the deblended sources which are brighter than 27 mag in -band and also brighter than mag, where is an -band magnitude brightest in the deblended sources which belong to the same parent. To avoid saturation in HSC photometry, no magnitude data brighter than 20 mag were used for any of the five bands.

### 2.3 AGN dataset selection

In the analysis of this paper, we treated each AGN and its surrounding galaxies as a set. Hereafter, we refer to the unit of the dataset as an AGN dataset. In this section we describe the criteria to include the AGN datasets for analysis.

For each AGN dataset, we measured the surface density of galaxies in annuli spaced by 0.2 Mpc out to 10 Mpc from the position of AGN. We kept only those AGN in which 60% of the area of all annuli at 2 Mpc around it and 80% of the area at 2Mpc were included in the survey footprint and not masked by bright sources. By this selection, among the original 6166 AGN datasets, 346 datasets were removed for the -band detected sample and 585 were removed for the four-band detected sample. Thus 5820 and 5581 AGN datasets passed this selection for -band and four-band detected samples, respectively.

The spatial uniformity of the galaxy samples in the AGN field was also examined to identify the AGN datasets that are significantly contaminated by nearby galaxies and stellar groups or showing non-uniformity for any other reason. For this purpose, we calculated two parameters for the radial number density distribution of galaxies, and , where is a square sum of the deviation from the number density distribution fitted to the observed data using equation (6), which will be derived in section 3.1, is a maximum deviation from the density distribution. Adapted criteria for those parameters are: and . By this selection, 268 (236) AGN datasets were removed from the 5820 (5581) AGN datasets for the -band (four-band) detected sample. Thus 5552 (5345) AGN datasets passed all the above selections.

Figure 2 shows the mass vs redshift distribution of the AGNs that were selected according to the conditions described above and used in the cross-correlation analysis in this paper. Figures 3 and 4 show histograms of redshift and black hole mass, respectively. We divided the redshift range into five groups as shown in Figure 2, which we call z0, z1, z2, z3 and z4 for 0.6–1.0, 1.0–1.5, 1.5–2.0, 2.0–2.5, and 2.5–3.0, respectively. For each redshift group, except for z4, the BH mass range was divided into two groups M8 and M9 so that each group has a similar number of AGNs rather than being divided at the same mass. This is because the effect of sample variance becomes dominant when the number of AGN samples is small and it is difficult to homogenize the samples among different mass groups by the division at constant mass for all the redshifts. The mass group was divided at = 8.4, 8.8, 9.0, and 8.9 for z0, z1, z2 and z4 redshift groups, respectively, which allows us to make the statistical uncertainties even for both mass groups. The mass dependence for the z4 redshift group was not examined, since the number of samples in the z4 group is too small to do so.

Table 2.3 shows the number of AGN datasets which passed all the selection criteria described above for each dataset. The -band detected galaxy samples were used for cross-correlation analysis, and the four-band detected samples were used for deriving color distribution and luminosity function of clustering galaxies, and they were also used to measure the dependence of cross-correlation length on galaxy luminosity. The redshift-matched samples were used for examining the BH mass dependence at each redshift.

## 3 Analysis method

### 3.1 Cross-correlation length

The cross-correlation function between AGNs and galaxies was calculated using the method described in our previous papers (Shirasaki et al., 2011; Komiya et al., 2013; Shirasaki et al., 2016). The analysis method is briefly described here.

We assumed the power-law form of the cross-correlation function,

 ξ(r)=(rr0)−γ, (3)

where is a cross-correlation length and is a power-law index fixed to 1.8, which is a typical value found in previous studies on AGN-galaxy cross-correlation studies (e.g., Hickox et al., 2009; Coil et al., 2009; Krumpe et al., 2015). The projected cross-correlation function is expressed in an analytical form as:

 ω(rp)=rp(r0rp)γΓ(12)Γ(γ−12)Γ(γ2), (4)

where is a projected distance from AGN and is the Gamma function. is related to the observed average surface number density of galaxies as:

 ω(rp)=n(rp)−nbg¯ρ0, (5)

where is an average surface number density of background galaxies integrated along the line of sight, and is an average of the number density of galaxies at AGN redshifts. From equations (4) and (5), the observed average surface number density of galaxies around AGNs can be expressed as:

 n(rp)=C(γ)¯ρ0rp(r0rp)γ+nbg, (6)

where . We fitted the model function of equation (6) to the surface number density derived by averaging over a group of AGN datasets, then obtained best fit parameters for the cross-correlation length and the background density .

In calculating , we masked the central 10” around each AGN to avoid blending problems with the AGN. Around a bright source there is a region where the number density of galaxy selected by the criteria described in section 2.2 is significantly reduced due to blending with the bright source and/or increase in the background noise level. That region needs to be counted as a dead region in calculating an effective area around AGNs.

To estimate the dead region, we adapted a bright source flag to remove the galaxies near bright sources, and the area of the masked region were calculated using the random catalog which was created avoiding the masked region with number density 100 arcmin. The random catalog we used is the one included in the S15b database. The bright source mask was applied only for those which are located 2 Mpc from the AGN as explained in section 2.2. To estimate the dead region at area of 2 Mpc, we calculated the correction factor from the average of the ratio of number density derived using masked () and un-masked () galaxy samples as follows:

To estimate in equation (6), for each AGN dataset was estimated from the luminosity function which was derived by parametrizing the luminosity functions in the literature. The luminosity function is expressed with the Schechter function (Schechter, 1976):

 ϕ(M)=0.4×ln10 ϕ∗10−0.4(M−M∗)(α+1)exp{−10−0.4(M−M∗)}. (9)

The parameters , , and were parametrized as a function of redshift at rest-frame wavelengths 150 nm, 280 nm, SDSS , , and band so that they approximated the parameters derived in literature (Gabasch et al., 2004, 2006; Dahlen et al., 2005, 2007; Parsa et al., 2016) for redshift = 0.5–3.0. In total 40 sets of luminosity function parameters from the literature were used to determine the parametrization. Then they were interpolated as a function of wavelength.

The redshift parametrization was performed with the following functions:

 ϕ∗=ϕ−180.4×ln10×100.4(18+M∗)(α+1)exp{−100.4(18+M∗)}, (10)
 ϕ−18=10(a0+a1z), (11)
 M∗=b0+b1z+b2z2, (12)
 α=c0+c1z, (13)

where represents luminosity density at . We used for the parametrization instead of using , because correlates with parameter and is strongly affected by the uncertainty of . , on the other hand, correlates with more weakly than , and its dependence on redshift and wavelength is rather small. The coefficients of equations (11)–(13) used in this analysis are summarized in Table 3.1. The comparison of these parametrization with the parameters in literature is shown in Figures 57.

The root mean square (RMS) of the difference between parameters from the literature and those calculated by parametrization of equations (10)–(13) are: 0.15, 0.28 mag, and 0.14 for , , and , respectively. The error of by 0.15 corresponds to a systematic error of by 20% for .

The error of significantly affects the luminosity densities at the bright end of , as the slope of the luminosity function becomes steeper there. Thus, comparison was made between the number densities at calculated using the parameters from the literature and those calculated with the above parametrization . According to the comparison, RMS of the relative error is 0.36 for 40 sets of luminosity functions of the literature, which corresponds to an upper deviation of cross-correlation length of , that is 28% of the cross correlation length, and a lower deviation of , i.e. 16%, for . When the comparison is made by extending to a lower luminosity side as , the RMS of the relative error is 0.16 for 20 sets of luminosity functions of , which corresponds to upper and lower deviation of 10% and 8% for , respectively.

Then is calculated as an integral of the product of the luminosity function and a detection efficiency DE() for apparent magnitude . The detection efficiency was modeled with the following function:

 DE(m)=⎧⎪ ⎪⎨⎪ ⎪⎩1(m

where and represent the threshold magnitudes for a decrease in the detection efficiency, and represent the attenuation widths. Two attenuation functions which have different threshold magnitudes and attenuation widths were required to fit the data accurately.

These parameters were obtained for each AGN dataset by fitting a model function to the observed magnitude distribution, at 21–27. We assumed a power law form for the model function of the magnitude distribution which is expected for an ideal observation of 100% detection efficiency at any magnitude:

 Norg(m)=c×10a(m−23), (15)

Then the function is fitted to the observed magnitude distribution . The residuals between them are within the statistical error.

Using the obtained from the fitting, is calculated as:

 ρ0=∫mupp−DMmlow−DMϕ(M)DE(M+DM)dM, (16)

where represents the distance modulus for the AGN redshift, and and represent the magnitude range of the galaxy samples.

### 3.2 Color and absolute magnitude distributions for galaxies

The color () and absolute magnitude () distributions for galaxies around AGNs were derived by subtracting the distribution in an offset region () from that at an central region of the AGN field () as follows (Shirasaki et al., 2016):

 n(D)=ncen(D)−noff(D) (17)
 n(M)=ncen(M)−noff(M) (18)

The offset region in this work is defined as an annular region at a projected distance from 7 to 9.8 Mpc from the AGN, and the central region is from 0.2 to 2 Mpc.

To calculate the color and absolute magnitude for the same rest frame bandpass at different redshifts, we performed SED fitting using the EAZY software developed by Brammer et al. (2008) and calculated color and magnitude at fixed bandpasses. For the redshift z0, z1, and z2 groups (0.6–2.0), the color is defined as and the absolute magnitude is , where , and represent absolute magnitudes at wavelengths 270, 380, and 310 nm, respectively. For the redshift z2, z3, and z4 groups (1.5–3.0), the color is defined as and the absolute magnitude is , where , , and represent absolute magnitudes at wavelengths 165, 270, and 210 nm, respectively. For redshift z2 group, we calculated the distributions for both definitions for comparison between them.

Although EAZY can be used to estimate photometric redshift (photo-z), it was used here to interpolate the observed SEDs of all the sources assuming its redshift to be the same as the AGN redshift. This is because the photo-z is not well constrained at redshifts explored in this work, as the structure around 400 nm, which is primarily used to determine the photo-z, is out of the observed passband. For this reason, photo-z wasn’t used in our analysis and the clustering of galaxies was measured by stacking the galaxy distribution around AGNs and subtracting the DC offset, which is a contribution from foreground/background galaxies, from the distribution.

## 4 Results

### 4.1 Cross-correlation function

First we present the cross-correlation functions derived by using the -band detected galaxy samples for each redshift group. To test also for the dependence on galaxy luminosity, two galaxy samples were constructed with different threshold magnitudes. The threshold magnitudes were chosen to be and , where () represents the absolute magnitude where averaged detection efficiency estimated as described in section 3.1 is 90% (50%). The absolute magnitude is measured in -band observer frame as calculated with , where is -band apparent magnitude and is distance modulus for the AGN redshift. The values of and are summarized in the second and third columns of Table 4.1.

Figure 8 shows the distributions of average number density of galaxies with absolute magnitude brighter than for five redshift groups. The model functions of equation (6) fitted to the observations are also plotted in the same figure. The corresponding projected cross-correlation functions calculated using equation (5) are shown in Figure 9. The estimated fitting parameters are summarized in upper parts of Table 4.1 for each threshold magnitude.

The statistical error of was estimated by the jackknife method. Jackknife resamplings were made by omitting, in turn, each of the AGN datasets. Then was calculated for each jackknife resampling and its error was estimated from their variance:

 σ2r0=N−1NN∑i=1(r0,i−¯r0)2, (19)

where is a cross-correlation length obtained for the i-th jackknife sample, is the average of , and is the number of jackknife samples.

The cross-correlation lengths are plotted as a function of redshift in Figure 10 for each threshold sample. The result shows that the cross-correlation length increases as the redshift increases, and also is larger for more luminous galaxy samples. Therefore, it is important for examining the BH mass dependence after matching the distribution of redshift and galaxy luminosity for different BH mass groups.

To reduce the effect of redshift dependence in the comparisons between the different BH mass groups, we constructed redshift-matched samples for each redshift group. The redshift-matched samples were constructed by selecting the same number of AGN datasets for each 0.02 bin.

The projected cross-correlation functions calculated for the redshift-matched samples are shown in Figure 11 for four redshift groups z0, z1, z2 and z3. The mass dependence for the z4 group was not examined due to poor statistics. In this figure, the circles represent lower BH mass groups (M8) and the squares represent higher BH mass groups (M9). Galaxy samples brighter than were used. The estimated fitting parameters are summarized in the bottom part of Table 4.1.

The cross-correlation lengths are plotted as a function of BH mass in Figure 12. From these results, no significant difference is seen between the two mass groups.

For comparison with the previous results (Shirasaki et al., 2016) derived using UKIDSS and SDSS data for galaxy samples, the result obtained for red galaxy samples () as described in the next subsection are also plotted. In the work of Shirasaki et al. (2016), they calculated cross-correlation functions between AGNs and galaxies, which were mostly red galaxies, at redshifts 1.0, and obtained the result that the cross-correlation length depends on BH mass but does not depend on redshift. Our current result at higher redshifts indicates its strong dependence on redshift and not on BH mass. The different behavior partly comes from the difference in the properties of galaxy samples; that is, the galaxy samples in the previous work are dominated by red-type galaxies and mostly dimmer than , while those in the current work are dominated by blue-type galaxies and mostly consist of galaxies brighter than at higher redshifts. As shown in the following subsection the clustering of galaxies significantly increases above , which introduces redshift dependence to the cross-correlation length.

### 4.2 Galaxy luminosity dependence of cross-correlation length

As larger clustering was indicated for more luminous galaxies as shown in Figure 10, we calculated cross-correlation lengths for flux-limited samples derived from the four-band detected samples to examine the dependence on the luminosity of galaxies more closely. To calculate the absolute magnitude for the same rest frame bandpass at different redshifts, we performed SED fitting using the EAZY software developed by Brammer et al. (2008) and calculated the magnitude at a fixed bandpass. The observed magnitudes at , , , , and bands are used in the SED fitting. The flux-limited galaxy samples were constructed based on absolute magnitude for z0, z1 and z2 and for z2, z3 and z4, where and represents the absolute magnitude at wavelength 310 nm and 220 nm, respectively.

Figure 13 shows the cross-correlation length as a function of average absolute magnitude of the galaxy sample brighter than the threshold magnitude in or , and Figure 14 shows it as a function of difference from the characteristic magnitude of the parametrized luminosity function derived in section 3.1. In these figures, dotted lines represent an empirical formula expressed as:

 r0=r0,min+exp(−M+20σM), (20)

where represents an asymptotic cross-correlation length at low luminosities, represents or , and corresponds to a slope of the curve at larger luminosities.

From these figures, it is apparent that the cross-correlation length rapidly increases at magnitudes brighter than , and it also increases with increasing redshift when compared for the same . corresponds to the cross-correlation length calculated for the whole galaxy samples of each redshift group. Due to the difference in the detection threshold of galaxy luminosity, the average luminosity for all the galaxy samples is larger for the higher redshift groups than for the lower ones. Thus the difference of is due to the difference in galaxy luminosity as well as the difference in redshift.

### 4.3 Color distributions of galaxies around AGNs

Color distributions of galaxies around AGNs were derived by the subtraction method described in section 3.2, and the results are shown in Figure 15 for color parameter at z0 to z2 and in Figure 16 for at z2 to z4. Redshift-matched AGN samples were used for this analysis. The galaxy samples used are the four-band detected samples, and the threshold magnitude for the galaxy was chosen to be , which is an absolute magnitude where average detection efficiency is 50%. The measured values of are summarized in the fourth and sixth columns of Table 4.1.

The panels at the left end column in the figures show the excess color density of or distributions for combined M8+M9 groups, and they are placed in increasing order of their redshifts from top to bottom. The two panels on the right hand side of each row are those for M8 and M9 groups at the corresponding redshift, respectively. In the case of the z4 group, only the distribution of the whole mass group is shown, as the number of samples is too small to split the sample into two mass groups for the analysis of mass dependence. is defined as a difference between magnitudes at rest frame wavelengths of 270 nm and 380 nm for redshift z0, z1 and z2 groups, and is a color between magnitudes at 165 nm and 270 nm for redshift groups of z2, z3 and z4. The magnitudes were calculated by performing SED fitting using the EAZY software Brammer et al. (2008) as described in section 3.2.

The color distributions were fitted with a double Gaussian model in which the two Gaussian functions correspond to red (smaller ) and blue (larger ) galaxy types. The fitted model functions are plotted in the same figures. As can be seen from the figure, blue galaxies dominate the clustered component in our galaxy samples.

In the fit to the distribution for z2-(M8+M9) group shown in Figure 15, the Gaussian parameters for each component were not well constrained. This is partly because there is a tendency for the distribution of the blue component to shift to the redder side as the redshift increases, and the components overlap in higher fractions for that group. Thus we fixed the Gaussian parameters of the red component to those obtained for z1-(M8+M9) group. In the case of z3 and z4 shown in Figure 16, the red component is barely visible, thus the Gaussian parameters for the red component were fixed to those obtained for distribution of z2-(M8+M9).

The fitting for the M8 and M9 samples were performed by fixing the Gaussian parameters to those obtained for the combined M8+M9 samples except for the z0 group. In the case of z0, the peak positions of the blue component are slightly different between the M8 and M9 groups. Thus, the Gaussian parameters for the blue component were derived independently. Those fitting parameters are summarized in Table 4.3.

The fractions of the blue galaxy for the combined mass groups are around 0.7–1.0, and there is a marginal indication that the blue fraction increases at higher redshifts. Care needs to be taken in comparison of the different redshifts by noting the difference in the galaxy samples. The galaxy samples at higher redshift are biased to more luminous galaxies, and it is also expected that the fraction of red galaxies increases along with the increase in luminosity as observed in the color-magnitude diagrams of the local galaxies. Thus the observed redshift dependence of the blue fraction is biased to the red component at higher redshift, and it is expected that the increase in the blue fraction is more or less enhanced if the bias is corrected.

As already mentioned above, the peak position of the blue galaxy is found to shift to the redder side as redshift increases. We also found that the peak position shifts to the redder side for brighter galaxies of the same redshift, so the shift of the peak position to the redder side at higher redshift is partly due to the dependence of the color on the absolute magnitude of the galaxy.

The obtained blue galaxy fractions are plotted in Figure 17 as a function of BH mass for each redshift group. A decreasing trend of blue fraction as BH mass is found for all redshifts except for z1. The result of z1 might have suffered from the sample variance. This can be tested by adding more samples from the upcoming HSC-SSP survey.

In the previous work of Shirasaki et al. (2016) that used SDSS and UKIDSS for galaxy samples, the blue galaxy fraction was less than  0.2 at redshift group z0 of this work, as the galaxy sample was strongly biased to red-type galaxies. In that work they also obtained the result that red galaxies more luminous than were strongly clustered around AGNs of  M. To test the consistency with the previous result on the cross-correlation with red galaxies, we calculated the cross-correlation lengths for red galaxies in the z0-M8 and z0-M9 groups separately.

We extracted red galaxies defined as , and restricted their brightness to . To calculate the cross-correlation function, we should know the average number density of these galaxy types at the corresponding redshift, . It was estimated as using the blue fraction determined for z0-M8 group. Then we obtained the cross-correlation lengths for red galaxy samples as for z0-M8, and for z0-M9. The values of are plotted in Figure 12 and compared with the previous work. They are consistent with each other.

### 4.4 Luminosity functions of galaxies around AGNs

Next we examined the absolute magnitude distributions of the galaxy clustered around AGNs. The left (right) panel of Figure 18 shows the distributions of () for z0, z1 and z2 (z2, z3 and z4) redshift groups. The excess densities were calculated as in equation (18) and corrected for their detection efficiencies calculated as described in section 3.1. They are plotted for .

In the same figure luminosity functions scaled by a factor corresponding to the estimated cross-correlation length are also plotted with solid lines. The luminosity function was calculated by using the parameterization derived in section 3.1. The scaling factor was calculated as:

 nIN−nOUT¯¯¯ρ0=2C(γ)rγ03−γ⎛⎝r3−γmax,IN−r3−γmin,INr2max,IN−r2min,IN−r3−γmax,OUT−r3−γmin,OUTr2max,OUT−r2min,OUT,⎞⎠ (21)

where () represents the surface number density at an inner (outer) region of the AGN field defined by annulus from to (from to ), and the right-hand side of the equation is calculated by integrating the equation (6). The boundary radii are , , , and . The arrows at the top of the figure indicate 90% detection limits. These results show overdensity against expectations from the luminosity function at magnitudes brighter than , which results in larger clustering of bright galaxies as obtained in section 4.2.

To test if the overdensity can be explained by the uncertainty of parameter used in the parametrization of luminosity function model, we measured the parameters for the observed magnitude distributions by fitting the Schechter function to them and compared them with the values given by the parametrization. In the fitting was fixed to . The results are summarized in Table 4.4 and compared in Figure 19 with the parametrization (thick solid lines) and literature values (open markers). The obtained are typically smaller (i.e. brighter) than given by the parameterization by 1 mag, well exceeding the uncertainty of , 0.28 mag, expected for the parametrization. Thus it is unlikely that the overdensity is attributable only to the uncertainty of the assumed luminosity function.

To investigate the galaxy type that contributes to the overdensity at the bright end, we derived the magnitude distributions for blue and red galaxies separately. They were classified at or . Figure 20 shows the comparisons between the magnitude distributions for blue and red galaxies. Both galaxy types contribute to the overdensity.

parameters measured for the distributions are summarized in Table 4.4 and plotted in Figure 21. The obtained parameters for red galaxies are systematically brighter than those for blue galaxies by more than two sigma at redshift z0, z1, and z2 for measurements at 310 nm. No statistically significant difference is found for the values for z2, z3, and z4 measured at 220 nm, which is mostly due to poor statistics and lower resolution of parameters on the deconvolution of the two components.

We also investigated whether there is a difference in the luminosity functions of clustering galaxies for lower and higher BH mass groups. Figure 22 shows the comparisons between the magnitude distributions for M8 and M9 mass groups, and parameters measured for them are summarized in Table 4.4 and plotted in Figure 23.

The obtained parameters for M9 group are systematically brighter than those for M8 group by more than two sigma at redshifts z0 and z1. No statistically significant difference is found for the values at redshift z2 (both for 310 nm and 220 nm) and z3. The results may indicate the enhancement of bright galaxies around AGNs with higher BH mass () at redshift 1.5, while at redshift 1.5 the luminosity function of galaxy is similar around AGNs with BH mass of .

### 4.5 AGN linear bias and dark matter halo mass hosting AGNs

The AGN-galaxy cross-correlation functions can be used to measure the linear bias of AGN distribution then to estimate the average mass of dark matter haloes hosting the AGNs. Using the relation between the halo mass and its bias obtained by numerical simulation, we can estimate the halo mass by assuming that the AGNs are located in haloes with a similar bias.

To derive the linear bias of AGNs from the AGN-galaxy cross-correlation function, we should know also about the auto-correlation function of galaxies. It is not, however, determined only from the dataset used in this work. Thus, we assume that auto-correlation functions measured by Zehavi et al. (2011), which were derived by using the SDSS galaxies at redshifts 0.25, doesn’t change up to redshift 3. Since the galaxy auto-correlation functions measured at higher redshifts (e.g., Ishikawa et al., 2015) are almost the same as those measured by Zehavi et al. (2011), it can be a good approximation.

The auto-correlation function of galaxies was assumed to be expressed as:

 ξGG=(1−fblue)ξGG,red+fblueξGG,blue (22)

where and are the auto-correlation function of red and blue galaxies. We adapted and from Zehavi et al. (2011). represents the fraction of blue galaxies and the values described in Table 4.3 were adapted. Then auto-correlation function of AGNs is calculated as:

 ξAA=ξ2AG/ξGG. (23)

In the previous subsections we showed that the galaxies around AGNs show larger clustering at higher luminosities (i.e. ) than ordinary galaxies (i.e. ). As the galaxy auto-correlation function used here is the one for , we also need to derive the AGN-galaxy cross-correlation function for galaxies.

The primary cause of the larger clustering of our galaxy samples is due to the evolution of, i.e. decrease in, the parameter in the luminosity function at neighbors of AGNs as shown in the section 4.4. To reduce the effect of the evolution and derive for galaxies, we calculated the cross-correlation lengths from the magnitude distributions of the clustering galaxies which were obtained by fitting the Schecher function to the observation in estimating for them. The parameter of the Schecher function was fixed to 1.2. The magnitude distributions were shown in Figure 18 as dashed lines.

is expressed as:

 nfit(M)=f(r′0)ϕ(M;M∗=M∗(clust)), (24)

where is the multiplying factor calculated as in equation (21) for a given correlation length , ) is the luminosity function as given by , which is measured for galaxies clustring around AGNs. The values of and are given by the parametrization used in this work. Thus , which is the cross-correlation length corrected for the evolution and thus the one for galaxies, is determined by solving the equation (24) for . The obtained values of are summarized in Table 4.5.

Using , can be expressed as:

 ξAG=(r/r′AG)−γ. (25)

calculated by equation (23) was fitted to the power law function of the form . The obtained and values are summarized in Table 4.5.

The linear bias of AGNs is calculated as:

 bAGN=σ8,AGNσ8,DM, (26)

where and are the rms fluctuations of the AGN and the dark matter density within spheres of comoving radius of 8 Mpc, respectively.

 σ8,AGN=J2(γ)1/2(rAA8)γ/2, (27)
 J2(γ)=72(3−γ)(4−γ)(6−γ)2γ, (28)

and is

 σ8,DM=σ8D(z)D(0), (29)

where is the linear growth factor given as

 D(z)=5ΩmE(z)2∫∞z1+yE3(y)dy, (30)
 E(z)2=Ωm(1+z)3+ΩΛ. (31)

Then is calculated from as:

 bAGN=(rAA8)γ/2J2(γ)1/2(σ8D(z)D(0))−1. (32)

The values of are summarized in Table 4.5.

The relation between dark matter halo mass and bias of dark matter halo is derived by Sheth et al. (2001) as:

 b=1+1√aδc[√a(aν2)+√ab(aν2)1−c−(aν2)c(aν2)c+b(1−c)(1−c/2)] (33)

where, = 0.707, = 0.5, = 0.6, = 1.686, which is a critical over density required for collapse in the spherical model. is defined as:

 ν=δcσ(Mh)D(0)D(z) (34)

is rms density fluctuation of dark matter halo with mass , and is parametrized by (van den Bosch, 2002):

 σ(Mh)=σ8f(u)f(u8) (35)

where with ,

 u=3.804×10−4Γ(MhΩm)1/3, (36)
 f(u)=64.087(1+1.074u0.3−1.581u0.4+0.954u0.5−0.185u0.6)−10 (37)

Thus the mass of the dark matter halo can be estimated using Equation (33) from the AGN bias parameter. The mass estimates are summarized in Table 4.5.

The obtained AGN biases are plotted as a function of redshift with solid circles in Figure 24. The AGN biases obtained in the literature (Hickox et al., 2009; Krumpe et al., 2012; Ross et al., 2009; Croom et al., 2005) are also shown for the comparison. In the figure, the bias expected for three halo masses, = 12.0, 12.5, and 13.0, is plotted with dashed lines. The obtained bias at redshift 2.0 is consistent with those obtained in other studies, and hosted by dark matter haloes with average mass in the range of M. The large scatters of our results at 2.0 is mostly due to the dominance of galaxies in our samples, which introduce large error in the correction of to .