The correlation of extragalactic rays
with cosmic matter density distributions
from weakgravitational lensing
Abstract
The extragalactic ray background (EGB) arises from accumulation of ray emissions from resolved and unresolved extragalactic sources as well as diffuse processes. It is important to understand the origin of the EGB in the context of cosmological structure formation. Known astrophysical ray sources such as blazars, starforming galaxies, and radio galaxies are expected to trace the underlying cosmic matter density distribution. We explore the correlation of the EGB from FermiLAT data with the largescale matter density distribution from the Subaru Hyper SuprimeCam (HSC) SSP survey. We reconstruct an unbiased surface matter density distribution at by applying weakgravitational lensing analysis to the firstyear HSC data. We then calculate the crosscorrelation. Our measurement suggests overall weak correlation at angular scales of 3060 arcmin. The result is essentially consistent with a null detection, and allows us to derive constraints on the statistical relationships between astrophysical ray sources and dark matter halos. We make forecasts for the detectability of the crosscorrelation in the final HSC data. Detection with a confidence can be achieved with the full survey covering 1,400 squared degrees planned for HSC SSP.
I Introduction
Since its discovery, the extragalactic ray has been one of the most important subjects in highenergy astrophysics Kraushaar et al. (1972). The extragalactic gammaray background (EGB) is partly explained by ray emission from resolved point sources outside the Milky Way. So far, the Large Area Telescope (LAT) onboard the Fermi Gammaray Space Telescope has detected ray sources from across the whole sky Acero et al. (2015a). Many of them are extragalactic, and hence contribute to the EGB. In addition to the resolved contributions, a diffuse and nearly isotropic emission, referred to as the isotropic gammaray background (IGRB), is necessary to provide a better fit to the EGB observations Sreekumar et al. (1998); Ackermann et al. (2015) (also see Ref. Fornasa and SÃ¡nchezConde (2015) for recent review of IGRB). Although the origin of the IGRB is still under debate, a significant fraction may originate from the emission from unresolved astrophysical sources, which are too faint to be detected on an individual basis. Possible candidates include blazars Col (2010); Harding and Abazajian (2012), starforming galaxies Thompson et al. (2006); Makiya et al. (2011), and radio galaxies Inoue (2011); Di Mauro et al. (2014).
Recently, Ref. Ackermann et al. (2015) measured the energy spectrum of the EGB in the range between 0.1 GeV to 820 GeV and showed that the EGB spectrum is well described by a power law with a photon index of and an exponential cut off of GeV. Improved modeling of a linear combination of astrophysical ray sources based on FermiLAT data can provide an adequate interpretation of the origin of the EGB. Ref. Ajello et al. (2015) showed that the cumulative emission of blazars, starforming galaxies, and radio galaxies can explain the EGB spectrum in the range 0.1–820 GeV with the cut off originating from the expected ray attenuation by the extragalactic background light Gould and SchrÃ©der (1966).
In the standard model of cosmic structure formation, galacticsize astrophysical sources are hosted by dark matter halos. Therefore, the spatial distribution of EGB emission is expected to be correlated with the underlying matter density distribution of the Universe Camera et al. (2013); Fornengo and Regis (2014). There already exist several pieces of supporting evidence that show such a correlation between the EGB and largescale structure (LSS) tracers. For example, Ref. Allevato et al. (2014) find spatial clustering of resolved blazars detected in the 2year allsky survey by FermiLAT. The cross correlations between IGRB and galaxies Xia et al. (2015) and galaxy clusters Branchini et al. (2017) have been also detected. However, these previous studies do not provide direct information on the relation between the EGB and the cosmic matter density, because the aforementioned tracers of LSS are known to be biased with respect to the underlying matter density.
Fortunately, gravitational lensing provides a physical and unbiased probe of the cosmic matter density Bartelmann and Schneider (2001). It is free from assumptions that relate the matter density to observable quantities. By contrast, e.g., for galaxies to be used as a tracer of largescale structure, some assumption needs to be made on the relation between the galaxy luminosity and mass, and even additional conditions need to be invoked such as hydrostatic equilibrium.
The correlation between the IGRB and largescale matter distribution has been detected using gravitational lensing of the cosmic microwave background (CMB), reported in Ref. Fornengo et al. (2015). In the concordance CDM cosmology, the CMB lensing effect is largely determined by the matter density distribution at redshift Lewis and Challinor (2006). Thus, it is still unclear whether there exists a strong correlation between the EGB and cosmic matter density at . Since extragalactic sources at lower redshift should appear brighter, a large fraction of the EGB is expected to originate from the structures (sources) at . This fact further motivates to use cross correlation analysis of the EGB with lowredshift structures.
In this paper, we explore the correlation between the EGB and the cosmic matter density distribution probed by an optical weak lensing survey. We achieve this by using galaxy imaging data taken from the Subaru HyperSuprime Cam (HSC) SSP survey. Unlike similar investigations in the literature Shirasaki et al. (2014, 2016); Troster et al. (2016), we work with the reconstructed surface matter density from the galaxy imaging data, i.e., lensing convergence. We then perform, for the first time, crosscorrelation analysis of the EGB (including all the resolved components) and the lensing convergence. Our measurement can provide a direct evidence of the correlation between the extragalactic rays and cosmic matter density at . The result will give important implications for the astrophysical origin of the EGB.
The rest of the paper is organized as follows. In Section II, we summarize the basics of EGB and gravitational lensing. In Section III, we describe the galaxyimaging and ray data used, and provide details of the crosscorrelation analysis. Our benchmark model of the cross correlation is discussed in Section IV. In Section V, we show the result of our crosscorrelation analysis, and discuss constraints on possible connection between astrophysical ray sources and LSS. Concluding remarks and discussions are given in Section VI. Throughout, we use the standard cosmological parameters with , the average matter density , the cosmological constant , and the amplitude of matter density fluctuations within , .
Ii Observables
ii.1 Extragalactic gammaray background
The ray intensity along a given direction can be written as Fornengo and Regis (2014),
(1) 
where denotes the radial comoving distance (function of redshift ), is the comoving angular diameter distance, is the relevant density field of the ray source, and is the window function. In this paper, we define so that the mean ray intensity can be expressed as,
(2) 
where represents an averaged quantity over the sky.
Suppose that an astrophysical source population can emit rays, we set the relevant field in Eq. (1) to be where denotes the ray intensity emitted from sources, is the dirac delta function in threedimensional space and we assume any source populations can be approximated as point source^{1}^{1}1This approximation should be valid as long as we focus on the correlation with the scale much larger than the actual size of astrophysical sources. The angular scales of interest in this paper is set to be larger than 30 arcmins.. In terms of the luminosity function, the window function for the source population can be written as Camera et al. (2013, 2015),
where is the observed ray energy, is the energy of the ray at redshift , is the ray spectrum, is the photon spectrum index, is the ray luminosity function, and the exponential factor in the integral takes into account the effect of ray attenuation during propagation owing to pair creation on diffuse extragalactic photons. For the ray optical depth , we adopt the model in Ref. Gilmore et al. (2012).
In the following, we summarize the basic quantities to compute the ray intensity using Eq. (II.1) for several astrophysical ray sources.
ii.1.1 Blazars
Following Ref. Ajello et al. (2015), we characterize the blazar population by means of a parametric description of their ray luminosity function and energy spectrum. The blazer ray luminosity function is defined as the number of sources per unit luminosity (defined in the rest frame of the source, for energies between 0.1 and 100 GeV). The luminosity function at is modeled as,
where the parameters of , , , , and the relation of have been calibrated with FermiLAT resolved blazars for a given redshiftevolution scenario Ajello et al. (2015). We assume that the redshift evolution of the luminosity function is expressed in the following form,
(5)  
where , and are assumed to depend on the ray luminosity. We use the same parameter values as in Ref. Ajello et al. (2015) in Eqs. (II.1.1) and (II.1.1). This model is called the luminositydensity dependent evolution scenario.
Blazars are known to have a characteristic, broad energy spectrum. We follow the ray energy spectrum proposed in Ref. Ajello et al. (2015),
(7) 
where and we determine the normalization factor so that can be defined in the energy range of 0.1100 GeV. When computing the ray intensity from blazars, we adopt , , , and .
Note that the above model of and can reasonably explain the FermiLAT data of resolved blazars. The expected contribution to the mean EGB intensity above GeV will be % Ajello et al. (2015).
ii.1.2 Starforming galaxies
For starforming galaxies, we assume the powerlaw energy spectrum , characteristic of the LATdetected starburst galaxies Ackermann et al. (2012). We obtain the luminosity function by scaling from the infrared luminosity function of galaxies measured in Ref. Rodighiero et al. (2010). The conversion of the infrared luminosity to the gammaray luminosity of the galaxies is given by
(8) 
where is defined in the range of 0.1100 GeV and is the infrared luminosity for 81000 m Ackermann et al. (2012). We adopt and .
ii.1.3 Radio galaxies
For radio galaxies, we follow the model of Ref. Di Mauro et al. (2014), which has established a correlation between the ray luminosity and the radiocore luminosity at 5 GHz. Using the correlation together with the radio luminosity function Willott et al. (2001), one can evaluate the contribution to the EGB from radio galaxies. We consider the bestfit relation from Ref. Di Mauro et al. (2014),
(9) 
where is defined between 0.1 and 100 GeV. We assume an average spectral index of 2.37 throughout this paper. and are adopted in Eq. (II.1).
ii.2 Gravitational lensing
ii.2.1 Basics
The weak gravitational lensing effect is commonly characterized by the distortion of image of a source object (galaxy) by the following matrix,
(10) 
where we denote the observed position of a source object as and the true position as . In the above equation, is the convergence, is the shear, and is the rotation. In the weak lensing regime (), one can relate the convergence field to the density contrast of underlying matter density field Bartelmann and Schneider (2001),
(11)  
where represents the source distribution normalized to .
ii.2.2 Reconstruction of convergence field
In optical imaging surveys, galaxies’ ellipticities are commonly used to estimate the shear component in Eq. (10). Since each component in the tensor can be expressed as the second derivative of the gravitational potential, one can reconstruct the convergence field from the observed shear field in Fourier space as
(13) 
where and are the convergence and shear in Fourier space, and is the wave vector with components and Kaiser and Squires (1993).
For a given source galaxy, one considers the relation between the observed ellipticity and the expected shear ,
(14)  
(15) 
where is the conversion factor to represent the response of the distortion of the galaxy image to a small shear Bernstein and Jarvis (2002), is the true value of cosmic shear, and and are the multiplicative and additive biases to assess possible systematic uncertainty in galaxy shape measurements. In practice, before employing the conversion in Eq. (13), one must first construct the smoothed shear field on grids Seitz and Schneider (1995),
(16) 
where locates the position of the th source galaxy, represents the inverse variance weight for the th source galaxy, and is assumed to be a Gaussian
(17) 
Using Eqs. (13) and (16), one can derive the smoothed convergence field from the observed imaging data through a Fast Fourier Transform (FFT).
The resulting field is mathematically equivalent to the filtered convergence field Schneider (1996) which is defined as
(18) 
where the smoothing filter for can be expressed as
(19) 
Iii Data
iii.1 FermiLat
We use years (from August 4, 2008 to September 4, 2015) of Pass 8 ULTRACLEANVETO photons with reconstructed energy in the 1 – 500 GeV range. Photons detected at zenith angles larger than 90 were excised in order to limit the contamination from rays generated by cosmicray interactions in the Earth limb. Moreover, data were filtered removing time periods when the instrument was not in skysurvey mode. Fermi Science Tools v10r0p5 and instrument response functions (IRFs) P8R2_ULTRACLEANVETO_V6 were used for this analysis. We included photons within a square region encompassing each of the HSC fields of view. Figure 1 shows smoothed EGB maps observed by the LAT above 1 GeV (details of the FermiLAT analysis and data processing are provided in the appendix).
iii.2 Subaru HyperSuprime Cam
Hyper SuprimeCam (HSC) is a widefield imaging camera on the prime focus of the 8.2m Subaru telescope Aihara et al. (2018); Miyazaki et al. (2018); Komiyama et al. (2018); Furusawa et al. (2018) (also see Kawanomoto, S., et al. 2018 in prep.). Among three layers in the HSC survey, the Wide layer will cover 1400 in five broad photometric bands (grizy) over 5–6 years, with excellent image quality of subarcsec seeing. In this paper, we use a catalog of galaxy shapes that has been generated for cosmological weak lensing analysis in the first year data release. We denote the data as HSCS16A. The details of galaxy shape measurements and catalog information are found in Ref. Mandelbaum et al. (2017a).
In brief, the HSCS16A shape catalog contains the shapes of 12 million galaxies selected from 137 measured with the reGaussianization method Hirata and Seljak (2003). The shape measurement is performed in the coadded band image and calibrated by simulated galaxy images similar to those used in GREAT3 Mandelbaum et al. (2015). The image simulation enables to carry out reliable shear calibration, since it takes into account realistic HSC PSFs and reproduces the observed distribution of galaxy properties with remarkable accuracy Mandelbaum et al. (2017b). We apply a conservative galaxy selection criteria for the first year science, e.g., and , giving an average raw number density of galaxies . The HSCS16A dataset consists of 6 patches; XMM, GAMA09H, GAMA15H, HECTOMAP, VVDS, and WIDE12H. While we present covergence maps for these individual patches separately, we combine our results on crosscorrelations for all 6 patches. Accurate photometric redshifts in the HSCS16A are suitable for tomographic approach as in Refs. Camera et al. (2015); Troster et al. (2016). While photometric redshifts are measured for the HSC galaxies using several different techniques, throughout this paper we use the MLZ photometric redshifts Tanaka et al. (2017). The redshift distribution of sources is commonly approximated as the sum of the probability distribution functions (PDFs) estimated from MLZ. After stacking the PDFs of the HSC galaxies, we find the mean source redshift to be . For the lensing tomography, we divide the source galaxies into three bins by their photometric redshift , , , and , corresponding to mean source redshifts of , , and , respectively.
We then reconstruct the smoothed convergence field from the HSCS16A data as described in Section II.2.2. We follow a similar approach to that of Ref. Oguri et al. (2017). Adopting a flatsky approximation, we first create a pixelized shear map for each of the 6 patches on regular grids with a grid size of 0.1 deg. We then apply the FFT and perform convolution in Fourier space to obtain the smoothed convergence field, which is referred to as the Emode convergence in the literature. The imaginary part of the reconstructed convergence corresponds to the Bmode convergence, , which can be an indicator of the existence of certain types of residual systematics in our weak lensing measurements.
In actual observations, there are missing galaxy shear data due to bright star masks and edges. Applying our method directly to such regions likely generates noisy maps. We thus determine the mask regions for each convergence map by using the smoothed number density map of the input galaxies with the same smoothing kernel as in Eq. (17). Then we mask all pixels with the smoothed galaxy number density less than 0.5 times the mean number density. In addition, we apply a conservative mask of deg about the Galactic plane in the crosscorrelation analysis presented in Section III.3. Basic characteristics of the HSCS16A data are summarized in Table 1. Figure 1 shows our reconstructed EGB (bottom panels) and convergence fields (upper panels) when we set the smoothing size of arcmin.
Name of patch  (No selection in )  ()  ()  ()  Effective area () 

GAMA09H  16.8  7.93  3.84  2.72  21.53 
GAMA15H  21.5  9.32  5.29  3.85  38.04 
HECTOMAP  18.1  8.10  4.65  3.04  18.78 
VVDS  19.2  8.30  4.71  3.56  24.49 
WIDE12H  23.0  10.1  5.91  4.13  16.16 
XMM  19.9  8.61  5.12  3.57  34.52 
iii.3 Estimator of cross correlation
In order to measure the crosscorrelation, we introduce the following simple statistic for our smoothed maps,
(20) 
with
(21) 
Here, represents the smoothed convergence map reconstructed as in Section II.2.2 (or defined in Eq. [18]), and the sum in Eq. (20) is taken over all unmasked pixels for individual HSCS16A fields. Hence, the quantity is a function of the smoothing scale . It also depends on the photometricredshift selection of source galaxies in gravitational lensing analysis. We also examine the similar quantity of to study residual systematic effects in the galaxy shape measurement.
To estimate the statistical error, we utilize 200 mock shear catalogs for the HSCS16A^{2}^{2}2 We chose the number of mock catalogs so that the inverse covariance matrix in our analysis has a accuracy. Ref. Hartlap et al. (2006) examine the bias in estimating the inverse covariance matrix as a function of independent random realizations. They find that the amplitude of bias depends on the ratio of the number of bins (data vector variables) to the number of data sets. There are 16 data variables in our cross correlation analysis. With 200 random realizations, we expect that the bias in the inverse covariance matrix is as small as . . We follow the method developed in Refs. Oguri et al. (2017); Shirasaki and Yoshida (2014) to create realistic mock catalogs that incorporate the features of actual data and raytracing simulations (also see Ref. Mandelbaum et al. (2017b) and Shirasaki et al. in prep). We use the allsky raytracing simulations of gravitational lensing in Ref. Takahashi et al. (2017). Each mock data consists of 38 different source planes each separated by a comoving separation of , covering source planes up to redshift of 5.3. The angular resolution is set to be 0.43 arcmin. We use the observed photometric redshifts and angular positions of real galaxies as in Ref. Oguri et al. (2017). In short, we perform the following procedures: (1) we assign each real HSC galaxy to the nearest angular pixel in the nearest redshift source plane, (2) we randomly rotate the orientation of the galaxy to remove the real lensing effect, (3) we simulate the lensing distortion effect at the source position by adding the local lensing shear and the intrinsic shape, (4) we include the additional variance due to measurement error, and (5) repeat the procedures (1)(4) for all the source galaxies. Note that, when generating a realization, we randomly sample the source redshift from the posterior PDF of photometric redshift for each galaxy. Our mock catalogs include directly the properties of source galaxies (e.g., magnitudes, ellipticities and spatial variations in the number densities), statistical uncertainties in photometric redshifts, and also the survey geometry. We use 10 fullsky simulations to extract 200 mock HSCS16A catalogs.
For each patch of the HSCS16A data, we evaluate the crosscorrelation Eq. (20) and estimate their errors by using the standard deviation from our 200 mock catalogs. Using the measured correlations and the standard deviations for a total of 6 patches, we combine them by the inversevariance weighted method Shirasaki et al. (2016); Oguri et al. (2017). The same procedures are also applied to the 200 mock catalogs, in order to estimate the covariance of in the HSCS16A among different and source redshift selections.
Iv Analytic Model
iv.1 Formulation
Here we derive the expectation value of Eq. (20) in the standard structure formation model. Let us begin with the twopoint crosscorrelation function between and in real space. Under a flatsky approximation, it is given by
(22)  
where characters with hats represent quantities in Fourier space, is the dirac delta function in twodimensional space, , and we define the cross power spectrum as,
(23) 
Here we use the relations of Eqs. (18) and (21) in Fourier space,
(24)  
(25) 
where and are defined as Eqs. (17) and (19), respectively. For a given filter of , its Fourier counterpart is written as,
(26) 
where is the zeroth Bessel function. Hence, we can obtain the expectation value of Eq. (20) by setting in the righthand side of Eq. (22).
The cross power spectrum can be expressed as Fornengo and Regis (2014); Camera et al. (2013, 2015),
(27)  
where represents the threedimensional cross power spectrum between the matter density and the number density of astrophysical source population in the EGB. Here we ignore possible contributions from truly diffuse processes such as darkmatter annihilation and decay^{3}^{3}3 Cosmological constraints of darkmatter annihilation with our measurements are provided in the Appendix, while our method seems more suitable for studying astrophysical contributions in the EGB. . Note that while Eq. (27) does not include the smearing effect due to the ray PSF, we can properly include the effect as in Ref. Shirasaki et al. (2014).
Before going to the detailed computation, we describe our smoothing filters in Fourier space and the window function along a line of sight. Figure 2 shows the effective smoothing filter in Eq. (22). According to this figure, our filters in Eqs. (17) and (19) can be regarded as highpass and lowpass in smoothing. In the range of arcmin, we can efficiently extract the information of clustering signature between Fourier modes with narrow range of . Furthermore, Figure 3 summarizes the relevant window functions along a line of sight in our analysis. Because of the nature of gravitational lensing, we can effectively remove the contribution from very lowredshift structures. The effective window function in the computation of is broad in redshift, but we can find that the effective redshift will be of the order of corresponding to Gpc. Therefore, the main contribution in Eq. (27) will come from the clustering between matter density and ray sources at Mpc, where the linear theory in structure formation is valid with reasonable accuracy.
iv.2 A model of bias of astrophysical sources
We introduce our fiducial model of astrophysical sources. Our model is largely based on the one developed in Ref. Camera et al. (2015). This model is found to be consistent with the observed cross correlation signals between the IGRB and the gravitationallensing effect in the CMB Fornengo et al. (2015). As seen above, the expected correlation can be mostly set by the linear part of the clustering. In this case, the threedimensional cross power spectrum can be approximated as , where is the linear matter power spectrum at redshift and is the effective bias for source population . We follow the standard halomodel approach to evaluate Camera et al. (2013, 2015). Assuming correlation between the ray luminosity and host halo mass of ray source, one can find,
(28) 
where is the luminosity function of the source population (here we already perform the integral over the photon index if needed), represents the relation of host halo and the ray luminosity, and is the linear halo bias. In this paper, we adopt the model of in Ref. Sheth and Tormen (1999).
For blazars, we use the model in Ref. HÃ¼tsi et al. (2014), which has introduced a simple powerlaw scaling , to the Xray luminosity so that the model can reproduce well the abundance of Xray selected AGNs at different redshifts (see their Fig. 6 and Table 1). To relate the Xray luminosity to the ray luminosity, we introduce the relation between the bolometric blazar luminosity and disk Xray luminosity as Inoue and Totani (2009). The parameter is set to be 4.21 as follows in Ref. Harding and Abazajian (2012).
For starforming and radio galaxies, constructing the possible relation between and is more uncertain since there exists limited data. In the case of starforming galaxy, we adopt the benchmark model in Refs. Camera et al. (2013, 2015) , where it is evaluated with reasonable approximation that the mass associated to the maximum luminosity will not exceed a maximum galactic mass. Ref. Camera et al. (2015) also constructed a model for the relation for radio galaxies by using data from 12 available AGNs and possible relation in mass between darkmatter halo and supermassive black hole (SMBH) HÃ¼tsi et al. (2014); Bandara et al. (2009). The bestfit relation can be approximated as . Although one can construct more physical models, such as examined in Ref. Camera et al. (2015), we expect that our main results will be weakly sensitive to the choice of the relation for starforming and radio galaxies. This is because the dominant contribution to the EGB is expected to be from blazars, and the window function for blazars is much larger than that of other contributions at redshifts relevant in our analysis. Note that the window function in EGB is now constrained by the measurement of the EGB energy spectrum Ackermann et al. (2015); Ajello et al. (2015) and it seems difficult to increase the contribution from starforming and radio galaxies furthermore in the range of 1500 GeV.
Figure 4 shows the expected signal from our fiducial model of astrophysical ray sources. This figure clearly demonstrates that the blazar population will have a much larger contribution to than the others. The difference will be an order of in our lensing analysis. These results strongly motivate us to constrain the physical properties of blazars with our cross correlation analysis.
iv.3 Relating blazars to dark matter halos
In the following, we consider two phenomenological models to relate blazars with dark matter halos on cosmological scales.
iv.3.1 Gammaray luminosityhalo mass relation
As in Section IV.2, we shall assume that ray luminosity and host halo mass follows a relation
(29) 
which is adopted in Ref. Ando et al. (2007). Eq. (29) can be inferred from the combination of several relations. They include: (1) an empirical relation between the blazar luminosity at the peak frequency (denoted as ) and the luminosity of optical emission lines Wang et al. (2002), (2) the relation, often assumed to be in the literature Inoue and Takahara (1996); Kataoka et al. (1999), and (3) the relation with to Wang et al. (2002), where is the Eddington luminosity which is determined by the black hole mass . In addition to these relations, we use the tight correlation between black hole mass and host halo mass of Ref. Ferrarese (2002), to derive Eq. (29). Finally, we obtain
(30)  
(31) 
with
(32) 
Eqs. (30) and (31) will provide a prior information of and . Following the observational constraints in Refs. Wang et al. (2002) and Ferrarese (2002), we consider the range of and throughout this paper. Note that the parameters and in this model can be constrained by our cross correlation analysis.
iv.3.2 Linear bias model
We can directly constrain the effective bias from observed cross correlation. In this paper, we consider a simple powerlaw model,
(33) 
We find that the halo bias with a fixed mass can be approximated by Eq. (33) in the redshift range of with % level accuracy when compared to the linear model in Ref. Sheth and Tormen (1999) and Eq. (33). The bestfit relation is found to be and for , and for , and and for . Therefore, when putting constraints with our cross correlation analysis, we set the prior range of parameters and to be and , respectively.
V Result
v.1 Cross correlation
Figure 5 summarizes our crosscorrelation measurements. The error bars in this Figure are evaluated with 200 mock catalogs of galaxy shapes in HSCS16A. As seen in the colored points, we confirm that the systematic uncertainty due to imperfect knowledge of galactic ray components is smaller than the current statistical uncertainty. This is because the ROI in the HSCS16A is sufficiently far from the Galactic plane and the galactic ray contribution from decays can be well constrained. Similar results have been found in Ref. Shirasaki et al. (2016). In addition, we also examine the cross correlation with the Bmode convergence which is commonly used as an indicator of systematic uncertainties in the galaxyshape measurement. The right set of 4 panels in Figure 5 show the cross correlation with . We define the significance of our measurements, i.e., null signals, as
(34)  
where represents the cross correlation for smoothing scale with source galaxy selection , and is the covariance estimated from 200 mock catalogs. We find that the cross correlation is consistent with null detection in all cases. Table 2 summarizes the in our measurements.
Tomographic bins  

No selection in  1.84 (4)  1.07 (4) 
0.82 (4)  1.14 (4)  
1.96 (4)  3.92 (4)  
4.31 (4)  0.61 (4)  
Combined  7.41 (16)  8.95 (16) 
v.2 Implications for blazarhalo connection
We then use the null detection of the cross correlation to place constraints on the blazarhalo connection. Assuming the observed data follow multivariate Gaussian with covariance matrix , we introduce the statistics,
(35)  
where represents theoretical template of interest for observed . Our theoretical prediction is computed as described in Section IV. We set the parameters of interest as in Eqs. (29) or (33). The data consists of four bins in with four different sourcegalaxy selections in the lensing analysis. Hence, we have degrees of freedom, whereas there are two parameters in our model, in Eq. (29) or in Eq. (33). The constraints on the parameters can be obtained by , corresponding to a 68% confidence level.
Let us first consider the model of Eq. (29). In this case, our cross correlation measurements yield
(36) 
Considering the galactic constraints in Ref. Ferrarese (2002), we may assume (the slope in relation is set to be to ). If , Eq. (36) reads . From this constraint, we can infer the typical halo mass of the blazar population contributing to the EGB. The average luminosity of blazars can be obtained as,
(37) 
where we find at relevant redshift of . Using our constraint, we expect that the average host halo mass of blazars at will be larger than , which is consistent with recent measurements of blazar spatial clustering Allevato et al. (2014).
Furthermore, we can constrain the effective blazar bias in Eq. (33). The constraint can be well approximated as,
(38) 
Figure 6 highlights our constraint on the blazar bias. We compare it with other measurements of bias of active galactic nuclei (AGN). In this figure, we consider two representative scenarios. One assumes that blazars tend to live in denser environments in LSS with typical host halo mass of and the redshift evolution in host halo mass will be negligible. In this scenario, the parameter will be close to 1. The second is that there are no redshift evolutions in blazar bias over , i.e., , corresponding to rapid redshift evolution in host halo mass. We find that our constraints for both scenarios are consistent with other measurements of AGN bias at different wavelengths. This implies that the unified hypothesis of AGN Urry and Padovani (1995) will still be supported, although the cross correlation analysis with future lensing surveys will allow us to put a tighter constraint on blazar bias (or detect it), leading to cosmological tests of the unified AGN scheme.
Figure 7 represents the expected signal of cross correlation between the EGB and the smoothed lensing convergence in the HSC final data release. Here we assume that the covariance will scale with survey area and that the survey coverage will become 10 times as large as HSCS16A. We find that the HSC final data can detect the cross correlation signal due to clustering of the blazar population in the EGB. In Figure 7, we consider the fiducial model of blazars in Section IV.2. When using source galaxies in the HSC final data with the redshift of , we can directly probe the cosmic matter density at , and the significance of the EGB correlation will reach a level. This indicates that we can measure the blazar bias at relevant redshifts with accuracy.
An independent bias measurement of AGN populations is needed to understand the nature of AGN power sources (e.g. Martini and Weinberg (2001)). Black hole accretion is among the most popular mechanism for powering AGN activities (e.g., Salpeter (1964)). In the Black hole accretion scenario, the lifetime of AGNs can be inferred from the balance between gravitational force and radiative pressure, while the expected lifetime contains unknown factors including radiative efficiency and mass accretion rate. Hence, observational determination of the AGN lifetime is a key task to improve our understanding of the nature of AGNs.
For a given redshift , the comoving number density of active AGNs can be expressed as where is the lifetime of AGN, is the age of Universe, and is the comoving number density of halos hosting AGNs. Since rare and massive objects are highly biased tracers of the underlying mass distribution Kaiser (1984); Bardeen et al. (1986), the bias measurement of AGNs can be used to constrain the typical halo mass of active AGNs and evaluate . Assuming the blazars in EGB will reside in halos with mass greater than , one can easily derive the relation between observed bias and lifetime as follows:
(39)  
(40) 
where is the halo mass function at redshift . In this paper, we adopt the model of in Ref. Sheth et al. (2001). Figure 8 shows the relation between and for different redshifts. Suppose that the blazar in EGB has the bias factor of 2, the cross correlation analysis between HSC final data and Fermi can determine the bias with uncertainty of . The gray filled regime in Figure 8 shows this expected constraint in blazar bias by our proposed statistics. Since the relevant redshift in our analysis should be , we can put a constraint of the lifetime of blazars to be , leading to years at . It is worth noting that our measurement includes all blazars contributing to EGB and thus our approach is complementary to other approaches using resolved objects (e.g. Allevato et al. (2014)).
Vi Conclusion and Discussion
In this paper, we have calculated the crosscorrelation of the EGB (i.e., the whole extragalactic ray radiation including resolved sources) and the cosmic matter density at for the first time. One may naively expect strong correlation between the two if the EGB is mostly contributed by known astrophysical sources Ajello et al. (2015). The crosscorrelation at degree angular scales is found to be consistent with a null detection given the current statistical uncertainties. The result still allows us to constrain possible relations between the dominant source population of the EGB and the underlying mass distribution. We have tested and confirmed that our measurement is robust in terms of residual systematic effects in both the ray and lensing analyses.
Provided that the astrophysical models in Section II.1 can reproduce the observed EGB spectrum well, we argue that the correlation of blazars and the largescale matter density dominates the measured crosscorrelation at degree scales. This allows deriving an estimate of the blazar bias with respect to the underlying matter density at . The clustering bias can be used to infer the host environment of the blazar population in the hierarchical structure formation model, and also to examine the unified AGN scenario Urry and Padovani (1995). Interestingly, our result is consistent with other observations of AGN bias at different wavelengths (FIG. 6).
Assuming that the statistical uncertainty of the crosscorrelation is reduced with survey area, we conclude that the HSC final data enables us to detect the correlation between blazars and the matter density at a level. This will offer an invaluable opportunity to identify the dominant contribution to the EGB at , and provide a stringent statistical test of the unified AGN scenario at ray energy scales. If the blazar bias is determined accurately, one can develop realistic models of blazar population in the context of structure formation in a similar manner to that applied to SDSS galaxies van den Bosch et al. (2007).
Identifying the EGB sources remains a major challenge. We need other information to reveal possible contributions from starforming galaxies, radio galaxies, and diffuse processes. Joint measurement of multiple crosscorrelation functions will be the key to this end. It is suggested that the crosscorrelation of local galaxies and the cosmic infrared background helps to infer the starforming activity in galaxies over a wider range of redshifts Xia et al. (2015); Feng et al. (2017); Cuoco et al. (2017), while observations of galaxy clusters, radio galaxies, and the CMB lensing effect will tighten the constraints on the contribution from misaligned AGN Fornengo et al. (2015); Branchini et al. (2017). Furthermore, a tomographic approach using ray energies and redshifts of LSS tracers will be viable with upcoming galaxy surveys Ando (2014); Camera et al. (2015). Ultimately, by putting all these information together, we will be able to place stringent constraints on the particle nature of dark matter. Our cosmological probe is complementary to ray measurements of local galaxies Hooper and Linden (2011); Ackermann et al. (2014).
Acknowledgements.
The authors thank Yutaka Komiyama for providing useful comments on the manuscript. SH is supported by the U.S. Department of Energy under award number desc0018327. NY acknowledges financial support from JST CREST (JPMJCR1414). Numerical computations presented in this paper were in part carried out on the generalpurpose PC farm at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. The Hyper SuprimeCam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. The PanSTARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the PanSTARRS Project Office, the MaxPlanck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the HarvardSmithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST1238877, the University of Maryland, and Eotvos Lorand University (ELTE). This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software athttp://dm.lsst.org
.
Based [in part] on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by the Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan.
Appendix A FermiLat Analysis Methods
The analysis of the full FermiLAT data sets encompassing each of the HSC regions was divided into four smaller regions, each of which was modeled by a diffuse ray background and ray point source model. Namely, we divided each region into four contiguous patches with overlapping boundaries for which we started the fit with a sky model that includes all pointlike and extended LAT sources listed in the 3FGL Acero et al. (2015b) catalog as well as with models for the Galactic diffuse and isotropic emission. In our baseline model the Galactic diffuse emission was modeled by the standard LAT diffuse emission model gll_iem_v06.fits, and as a proxy for the residual background and extragalactic ray radiation spectra we used a single isotropic component with the spectral shape given by model iso_P8R2_ULTRACLEANVETO_V6_v06.txt^{4}^{4}4http://fermi.gsfc.nasa.gov/ssc/.
We used the Fermipy Python package^{5}^{5}5 http://fermipy.readthedocs.io/en/latest/. in conjunction with standard Fermi Science Tools to fit and characterize the sources included in the sky models corresponding to each region of interest (ROI). In each ROI we bin the LAT data in 8 energy bins per decade and with a spatial pixel size of . Each ROI is analyzed separately; however, in order to perform our cross correlation analysis in the full HSC regions we merge the patches from the analyses of the different ROIs.
Using the initial baseline model mentioned above, the first step of our method is to find the best spectral parameters for all free sources within our ROIs. To obtain convergence, all the fits were performed hierarchically; freeing first the normalization of the sources with the highest intensities followed by the lower ones within the ROIs. The fitting consecutively restarts from the updated bestfit models and repeats the same procedure this time for the spectral shape parameters. Due to the years of LAT data used in our analysis compared to the 4 years of the 3FGL Acero et al. (2015b) catalog, we expected to find new ray point source candidates in our ROIs. The significance of each new point source candidate was evaluated using the test statistic , where and are the likelihoods of the background (null hypothesis) and the hypothesis being tested (alternative hypothesis: background plus source). All point source candidates that were found to have a TS were included to our sky model of the respective ROI. The new Fermipy package includes special routines to perform automated point source searches, and in this study we followed the standard method recommended in the Fermipy documentation. We then used this package to refine the positions and the spectral parameters of the new point source candidates. Since the set of four patches making up every HSC ROI have some overlapping regions at the boundaries, we remove any duplicate point source candidates found in more than one patch. In this study we detected 48 new point source candidates with TS in the combined the HSC ROIs. The positions and spatial overlaps with the RomaBZCAT Massaro et al. (2009) Blazar catalog are shown in Table 3. All the characteristics of the new point source candidates are provided in FITS files with this article.
The EGB photons in the ROIs were obtained by subtracting the bestfit Galactic diffuse emission model from the photon counts maps. We note that the EGB images obtained in this way could still contain some isotropic detector backgrounds. However, our analysis is able to reproduce well the EGB rays derived by the FermiLAT team (see, e.g., Ref. Shirasaki et al. (2016) for an example of this method). This shows that most of the detector cosmicray induced backgrounds are safely removed by our conservative photon selection filters.
We estimated the systematic uncertainties in the Galactic diffuse emission model in a similar fashion to that explored in the Fermi collaboration paper on the EGB Ackermann et al. (2015). Namely, we reran our pipeline using the alternative foreground Models A, B and C as described in the appendix of Ref. Ackermann et al. (2015). Such foreground models encompass a very generous range of the systematics associated with this kind of analysis, and provides a test in rigor comparable to that performed by the Fermi team.
Name  Association  TS  

[deg]  [deg]  or spatial overlap  
PSJ0922.8+0432  227.74  35.48  44  
PSJ0932.7+1040  222.40  40.57  5BZBJ0932+1042  28 
PSJ0930.7+0031  233.24  35.02  5BZQJ0930+0034  76 
PSJ0833.60455  229.77  20.07  35  
PSJ0848.00704  233.77  21.97  5BZBJ08470703  101 
PSJ0819.50755  230.58  15.53  5BZBJ08190756  36 
PSJ1420.2+0612  352.20  60.28  5BZBJ1420+0614  37 
PSJ1410.1+0202  343.20  58.62  5BZBJ1410+0203  139 
PSJ1356.6+0237  338.22  60.92  48  
PSJ1443.50959  342.98  43.99  34  
PSJ1405.10642  333.41  51.74  30  
PSJ1359.10658  331.03  52.17  29  
PSJ1401.00914  330.23  49.90  5BZQJ14010916  61 
PSJ1425.30119  345.15  53.66  5BZBJ14250118  33 
PSJ1644.2+4544  71.33  40.84  5BZGJ1644+4546  42 
PSJ1607.8+4947  77.85  46.40  36  
PSJ1604.2+5008  78.59  46.91  5BZBJ1603+5009  54 
PSJ1607.8+4950  77.90  46.39  42  
PSJ1536.0+3744  60.58  54.02  5BZUJ1536+3742  253 
PSJ1529.5+3813  61.70  55.23  5BZBJ1529+3812  29 
PSJ1602.1+3328  53.77  48.71  5BZUJ1602+3326  33 
PSJ2249.9+0451  75.76  46.59  31  
PSJ2226.5+0207  67.11  44.48  28  
PSJ2235.20631  59.21  51.68  38  
PSJ2301.00157  71.81  53.49  5BZQJ23010158  136 
PSJ2225.70803  55.01  50.63  45  
PSJ2156.10036  57.71  40.32  5BZUJ21560037  127 
PSJ2211.00004  61.32  42.96  5BZBJ22110003  35 
PSJ2148.00734  48.46  42.42  5BZBJ21480733  80 
PSJ1219.7+0446  282.89  66.39  5BZBJ1219+0446  47 
PSJ1219.7+0547  282.00  67.39  37  
PSJ1215.1+0734  277.57  68.63  5BZGJ1215+0732  53 
PSJ1216.1+0931  276.01  70.52  5BZGJ1216+0929  58 
PSJ1216.00243  285.59  58.96  5BZBJ12160243  36 
PSJ1207.60104  280.65  59.89  5BZQJ12070106  69 
PSJ1136.00426  270.09  53.55  5BZQJ11350428  45 
PSJ1129.30530  268.53  51.82  59  
PSJ1131.00945  272.35  48.29  45  
PSJ0237.3+0206  168.27  51.19  30  
PSJ0239.8+0417  166.90  49.10  5BZQJ0239+0416  95 
PSJ0245.00257  176.18  53.65  55  
PSJ0208.60045  161.24  57.79  5BZBJ02080047  38 
PSJ0220.90841  176.01  61.94  5BZBJ02200842  38 
PSJ0226.60553  174.07  58.96  63  
PSJ0241.00506  177.65  55.86  5BZBJ02400504  51 
PSJ0205.80955  172.03  65.42  33  
PSJ0140.70758  156.59  67.57  5BZBJ01400758  42 
PSJ0142.60543  154.87  65.36  5BZBJ01420544  35 
New point source candidates in the HSC ROIs found with a in years of FermiLAT data. This table is also provided as a FITS file in the online material. The horizontal dividers separate the new point source candidates in the GAMA09H, GAMA15H, HECTOMAP, VVDS, WIDE12H and XMM patches, from top to bottom. The fourth column displays associations or spatial overlaps with tentative multiwavelength counterparts. For this purpose we use the RomaBZCAT blazar catalog Massaro et al. (2009). To identify possible multiwavelength counterparts to the new point source candidates we searched in the seed locations within the 68% containment of the point spread function for one of our highenergy bands, which comes to . 
Appendix B CrossCorrelation With Isotropic GammaRay Components and Constraints on Dark Matter Annihilation
In this appendix, we present the crosscorrelation analysis of the lensing convergence and the isotropic ray background (IGRB) including the unresolved components in FermiLAT data.
The lensing convergence is equivalent to the cosmic surface mass density with a weight function along lineofsight direction . The details of the reconstruction algorithm of from the HSCS16A galaxyimaging data are given in Section II.2.2. We use the crosscorrelation estimator in Eq. (20) to the IGRB measured by FermiLAT and the lensing convergence from HSCS16A. As in Section III.3, we evaluate the statistical uncertainty of the cross correlation with 200 HSCS16A mock catalogs and use the inverseweighed method to combine the crosscorrelation signals for individual HSCS16A patches together. In addition to the normal masks as described in Section III.3, we also place circular masks around FermiLAT ray point sources, with a radius of 1 deg. The circular masks are introduced in order to minimize possible contaminants from residual  rays from resolved bright sources.
Tomographic bins  

No selection in  1.24 (4)  5.50 (4) 
0.72 (4)  6.57 (4)  
1.13 (4)  5.87 (4)  
4.24 (4)  0.99 (4)  
Combined  7.06 (16)  16.6 (16) 
The left panels in Figure 9 summarize our cross correlation measurements of the IGRB and . We work with four different sourcegalaxy selections in the lensing analysis. Regardless of the choice of smoothing scales and galaxy selections, we find that the IGRB crosscorrelation is consistent with a null detection. We also confirm that the crosscorrelation measurements of the IRGB and the Bmode convergence is consistent with a null detection, suggesting that our lensing analysis is not compromised by systematic errors. Table 4 summarizes the significance of our results using the definition in Eq. (34).
We can derived constraints on dark matter (DM) annihilation directly using the cosmic matter density distributions. Using the formulation as described in Section IV.1, we compute the expected correlation originating from DM annihilation. We adopt the model of cross powerspectrum between the matter density and rays from DM annihilation as in Ref. Shirasaki et al. (2016) (see their Eq. [14]). We then consider a conservative model of substructures in DM halos SánchezConde and Prada (2014). Following a similar likelihood analysis as in Section V.2, we can constrain the two parameters in DM annihilation, DM particle mass and the annihilation crosssection , for a given annihilation channel. The right panel in Figure 9 shows our cosmological constraints on DM annihilation. For comparison, we also present the constraints from similar analyses with shallower but wider galaxy survey data Troster et al. (2016), and from analyses of the MilkyWay satellite galaxies Ackermann et al. (2014). Since our measurements in the HSCS16A cover a smaller sky area than Ref. Troster et al. (2016), our constraints are found to be weaker. Also it is unlikely for us to be able to place constraint on DM annihilation at a similar level as the local measurements Ackermann et al. (2014), even if we use the HSC final data release, because the increased sky coverage can improve the statistical uncertainty only by a factor of .
There are two reasons why our analysis does not provide tight constraints on DM annihilation:

Our crosscorrelation analysis uses linearscale clustering alone (see Section IV.1) that contains less information on the internal mass distribution of DM halos.

HSC is suited for measuring the matter distributions in the distant universe. The median source redshift in HSC is larger than those in other data in the literature. Crosscorrelation of the IGRB with objects or the matter density at lower redshifts can be more optimal for indirect detection of DM annihilation (also see Ref. Shirasaki et al. (2015) for optimal targets in largescale structure for indirect detection of DM annihilation).
Overall, the crosscorrelation measurement presented in this paper is a more promising tool to explore expected correlations of astrophysical ray sources with the cosmic matter density.
References
 Kraushaar et al. (1972) W. L. Kraushaar, G. W. Clark, G. P. Garmire, R. Borken, P. Higbie, V. Leong, and T. Thorsos, Astrophys.J. 177, 341 (1972).
 Acero et al. (2015a) F. Acero et al. (FermiLAT), Astrophys. J. Suppl. 218, 23 (2015a), arXiv:1501.02003 [astroph.HE] .
 Sreekumar et al. (1998) P. Sreekumar et al. (EGRET Collaboration), Astrophys.J. 494, 523 (1998), arXiv:astroph/9709257 [astroph] .
 Ackermann et al. (2015) M. Ackermann et al. (FermiLAT), Astrophys. J. 799, 86 (2015), arXiv:1410.3696 [astroph.HE] .
 Fornasa and SÃ¡nchezConde (2015) M. Fornasa and M. A. SÃ¡nchezConde, Phys. Rept. 598, 1 (2015), arXiv:1502.02866 [astroph.CO] .
 Col (2010) Astrophys. J. 720, 435 (2010), arXiv:1003.0895 [astroph.CO] .
 Harding and Abazajian (2012) J. P. Harding and K. N. Abazajian, JCAP 1211, 026 (2012), arXiv:1206.4734 [astroph.HE] .
 Thompson et al. (2006) T. A. Thompson, E. Quataert, and E. Waxman, Astrophys. J. 654, 219 (2006), arXiv:astroph/0606665 [astroph] .
 Makiya et al. (2011) R. Makiya, T. Totani, and M. A. R. Kobayashi, Astrophys. J. 728, 158 (2011), arXiv:1005.1390 [astroph.HE] .
 Inoue (2011) Y. Inoue, Astrophys. J. 733, 66 (2011), arXiv:1103.3946 [astroph.HE] .
 Di Mauro et al. (2014) M. Di Mauro, F. Calore, F. Donato, M. Ajello, and L. Latronico, Astrophys. J. 780, 161 (2014), arXiv:1304.0908 [astroph.HE] .
 Ajello et al. (2015) M. Ajello et al., Astrophys. J. 800, L27 (2015), arXiv:1501.05301 [astroph.HE] .