The correlation of extragalactic \gamma-rays with cosmic matter density distributions from weak-gravitational lensing

# The correlation of extragalactic γ-rays with cosmic matter density distributions from weak-gravitational lensing

Masato Shirasaki National Astronomical Observatory of Japan (NAOJ), Mitaka, Tokyo 181-8588, Japan    Oscar Macias    Shunsaku Horiuchi Center for Neutrino Physics, Department of Physics, Virginia Tech, Blacksburg, Virginia 24061, USA    Naoki Yoshida Department of Physics, University of Tokyo, Tokyo 113-0033, Japan
Kavli Institute for the Physics and Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan
CREST, Japan Science and Technology Agency, 4-1-8 Honcho, Kawaguchi, Saitama, 332-0012, Japan
Chien-Hsiu Lee Subaru Telescope, NAOJ, 650 N Aohoku Pl. Hilo, HI 96720, USA    Atsushi J. Nishizawa Institute for Advanced Research, Nagoya University, Nagoya 464-8602, Aichi, Japan
###### 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, star-forming galaxies, and radio galaxies are expected to trace the underlying cosmic matter density distribution. We explore the correlation of the EGB from Fermi-LAT data with the large-scale matter density distribution from the Subaru Hyper Suprime-Cam (HSC) SSP survey. We reconstruct an unbiased surface matter density distribution at by applying weak-gravitational lensing analysis to the first-year HSC data. We then calculate the cross-correlation. Our measurement suggests overall weak correlation at angular scales of 30-60 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 cross-correlation 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 high-energy astrophysics Kraushaar et al. (1972). The extragalactic gamma-ray 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 Gamma-ray 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 gamma-ray 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Ã¡nchez-Conde (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), star-forming 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 Fermi-LAT data can provide an adequate interpretation of the origin of the EGB. Ref. Ajello et al. (2015) showed that the cumulative emission of blazars, star-forming 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, galactic-size 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 large-scale structure (LSS) tracers. For example, Ref. Allevato et al. (2014) find spatial clustering of resolved blazars detected in the 2-year all-sky survey by Fermi-LAT. 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 large-scale 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 large-scale 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 low-redshift 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 Hyper-Suprime 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, cross-correlation 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 galaxy-imaging and -ray data used, and provide details of the cross-correlation analysis. Our benchmark model of the cross correlation is discussed in Section IV. In Section V, we show the result of our cross-correlation 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 gamma-ray background

The -ray intensity along a given direction can be written as Fornengo and Regis (2014),

 Iγ(\boldmathθ)=∫dχWg(χ)g(r(χ)\boldmathθ,χ), (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,

 ⟨Iγ⟩=∫dχWg(χ), (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 three-dimensional space and we assume any source populations can be approximated as point source111This 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),

 Wg,X(χ) = ∫EmaxEmindEγe−τ(E′γ,z(χ)) ×∫ΓmaxΓmindΓXdNγ,X(E′γ,z(χ),ΓX)dEγ ×∫LmaxLmindLγLγΦX(Lγ,z(χ),ΓX),

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,

 Φb(Lγ,z=0,Γb) = Aln10Lγ[(LγL∗)γa+(LγL∗)γb]−1 ×e−[Γb−μ(Lγ)]2/2σ2,

where the parameters of , , , , and the relation of have been calibrated with Fermi-LAT resolved blazars for a given redshift-evolution scenario Ajello et al. (2015). We assume that the redshift evolution of the luminosity function is expressed in the following form,

 Φb(Lγ,z,Γb) = fe(z,Lγ)Φb(Lγ,z=0,Γb), (5) fe(z,Lγ) = [(1+z1+zc(Lγ))−p1(Lγ) +(1+z1+zc(Lγ))−p2(Lγ)]−1,

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 luminosity-density 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),

 dNγ,bdEγ=K⎡⎣(EγEb(Γb))a+(EγEb(Γb))b⎤⎦−1, (7)

where and we determine the normalization factor so that can be defined in the energy range of 0.1-100 GeV. When computing the -ray intensity from blazars, we adopt , , , and .

Note that the above model of and can reasonably explain the Fermi-LAT data of resolved blazars. The expected contribution to the mean EGB intensity above GeV will be % Ajello et al. (2015).

#### ii.1.2 Star-forming galaxies

For star-forming galaxies, we assume the power-law energy spectrum , characteristic of the LAT-detected 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 gamma-ray luminosity of the galaxies is given by

 log(Lγergs−1)=1.17log(LIR1010L⊙)+39.28, (8)

where is defined in the range of 0.1-100 GeV and is the infrared luminosity for 8-1000 m Ackermann et al. (2012). We adopt and .

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 radio-core 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 best-fit relation from Ref. Di Mauro et al. (2014),

 log(Lγergs−1)=1.008log(Lr,coreergs−1)+2.00, (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,

 Aij=∂xitrue∂xjobs≡(1−κ−γ1−γ2−ω−γ2+ω1−κ+γ1), (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),

 κ(\boldmathθ) = ∫∞0dχWκ(χ)δm(r(χ)\boldmathθ,χ), (11) Wκ(χ) = 32(H0c)2Ωm0(1+z(χ))r(χ) ×∫∞χdχ′p(χ′)r(χ′−χ)r(χ′),

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

 ^κ(\boldmathℓ)=ℓ21−ℓ22ℓ21+ℓ22^γ1(\boldmathℓ)+2ℓ1ℓ2ℓ21+ℓ22^γ2(\boldmathℓ), (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 ,

 ~γα = eobs,α2R, (14) ~γα = (1+m)γtrue,α+cα, (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),

 γsm,α(\boldmathθ)=∑iwi(~γi,α−ci,α)W(\boldmathθ−\boldmathθi)∑iwi(1+mi)W(% \boldmathθ−\boldmathθi), (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

 κsm(\boldmathθ)=∫d2ϕκ(\boldmathϕ)Uκ(\boldmathθ−% \boldmathϕ), (18)

where the smoothing filter for can be expressed as

 Uκ(θ)=2∫∞θdθ′W(θ′)θ′−W(θ). (19)

## Iii Data

### iii.1 Fermi-Lat

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 cosmic-ray interactions in the Earth limb. Moreover, data were filtered removing time periods when the instrument was not in sky-survey 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 Fermi-LAT analysis and data processing are provided in the appendix).

### iii.2 Subaru Hyper-Suprime Cam

Hyper Suprime-Cam (HSC) is a wide-field 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 sub-arcsec 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 re-Gaussianization 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 cross-correlations 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 flat-sky 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 E-mode convergence in the literature. The imaginary part of the reconstructed convergence corresponds to the B-mode 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 cross-correlation 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.

### iii.3 Estimator of cross correlation

In order to measure the cross-correlation, we introduce the following simple statistic for our smoothed maps,

 ⟨κEIγ⟩(θG)≡∑iκsm(\boldmathθi)Iγ,sm(% \boldmathθi)∑i1, (20)

with

 Iγ,sm(\boldmathθ)=∫d2ϕIγ(\boldmathϕ)W(\boldmathθ−% \boldmathϕ). (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 photometric-redshift 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 HSCS16A222 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 ray-tracing simulations (also see Ref. Mandelbaum et al. (2017b) and Shirasaki et al. in prep). We use the all-sky ray-tracing 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 full-sky simulations to extract 200 mock HSCS16A catalogs.

For each patch of the HSCS16A data, we evaluate the cross-correlation 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 inverse-variance 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 two-point cross-correlation function between and in real space. Under a flat-sky approximation, it is given by

 ⟨κsm(\boldmathθ1)Iγ,sm(\boldmathθ2)⟩ = ∫d2ℓ1(2π)2∫d2ℓ2(2π)2ei\boldmathℓ1⋅\boldmathθ1−i\boldmathℓ2⋅\boldmathθ2 (22) ×⟨^κsm(\boldmathℓ1)^I∗γ,sm(% \boldmathℓ2)⟩, = ∫d2ℓ1(2π)2ei\boldmathℓ1⋅\boldmathθ12CκIγ(ℓ1) ×^W(ℓ1,θG)^Uκ(ℓ1,θG),

where characters with hats represent quantities in Fourier space, is the dirac delta function in two-dimensional space, , and we define the cross power spectrum as,

 ⟨^κ(\boldmathℓ1)^Iγ(\boldmathℓ2)⟩ ≡ CκIγ(ℓ1)(2π)2δ(2)(% \boldmathℓ1−\boldmathℓ2). (23)

Here we use the relations of Eqs. (18) and (21) in Fourier space,

 ^κsm(\boldmathℓ)=^Uκ(ℓ,θG)^κ(\boldmathℓ), (24) ^Iγ,sm(\boldmathℓ)=^W(ℓ,θG)^Iγ(\boldmathℓ), (25)

where and are defined as Eqs. (17) and (19), respectively. For a given filter of , its Fourier counterpart is written as,

 ^f(ℓ)=∫∞0ℓdℓ2πf(θ)J0(ℓθ), (26)

where is the zero-th Bessel function. Hence, we can obtain the expectation value of Eq. (20) by setting in the right-hand side of Eq. (22).

The cross power spectrum can be expressed as Fornengo and Regis (2014); Camera et al. (2013, 2015),

 CκIγ(ℓ) = ∑X∫dχr2(χ)Wg,X(χ)Wκ(χ) (27) ×Pmg,X(ℓ+1/2r(χ),z(χ)),

where represents the three-dimensional 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 dark-matter annihilation and decay333 Cosmological constraints of dark-matter 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 high-pass and low-pass 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 low-redshift 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 gravitational-lensing 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 three-dimensional 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 halo-model approach to evaluate Camera et al. (2013, 2015). Assuming correlation between the -ray luminosity and host halo mass of -ray source, one can find,

 beff,X(z)=∫dLγbh(Mh,X(Lγ,z),z)LγΦX(Lγ,z)∫dLγLγΦX(Lγ,z), (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 power-law scaling , to the X-ray luminosity so that the model can reproduce well the abundance of X-ray selected AGNs at different redshifts (see their Fig. 6 and Table 1). To relate the X-ray luminosity to the -ray luminosity, we introduce the relation between the bolometric blazar luminosity and disk X-ray luminosity as Inoue and Totani (2009). The parameter is set to be 4.21 as follows in Ref. Harding and Abazajian (2012).

For star-forming and radio galaxies, constructing the possible relation between and is more uncertain since there exists limited data. In the case of star-forming 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 dark-matter halo and super-massive black hole (SMBH) HÃ¼tsi et al. (2014); Bandara et al. (2009). The best-fit 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 star-forming 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 star-forming and radio galaxies furthermore in the range of 1-500 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 Gamma-ray luminosity-halo mass relation

As in Section IV.2, we shall assume that -ray luminosity and host halo mass follows a relation

 log(Lγ1044.4ergs−1)=A0+A1log(Mh1012M⊙), (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

 A0 = 0.64[log(Llines10−3LEdd)+log(α00.1)] (30) −log(LpkLγ), A1 = 1.088(α11.7), (31)

with

 log(MBH108M⊙)=α0+α1log(Mh1012M⊙). (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 power-law model,

 logbeff,b(z)=b0+b1log(1+z). (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 best-fit 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 cross-correlation 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 B-mode convergence which is commonly used as an indicator of systematic uncertainties in the galaxy-shape 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

 χ20 = ∑s,s′∑i,j\boldmathC−1ij(s,s′) (34) ×⟨κEIγ⟩(θG,i;s)⟨κEIγ⟩(θG,j;s′),

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.

### v.2 Implications for blazar-halo connection

We then use the null detection of the cross correlation to place constraints on the blazar-halo connection. Assuming the observed data follow multivariate Gaussian with covariance matrix , we introduce the statistics,

 χ2 = ∑s,s′∑i,j\boldmathC−1ij(s,s′){⟨κEIγ⟩(θG,i;s)−μi,s(\boldmathp)} (35) ×{⟨κEIγ⟩(θG,j;s′)−μj,s′(\boldmathp% )},

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 source-galaxy 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

 A1<−0.18A0+1.80(1σlevel). (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,

 ⟨Lγ⟩(z)=∫dLγLγΦb(Lγ,z)∫dLγΦb(Lγ,z), (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,

 b1<−11.77b0.0970+14.62(1σlevel). (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:

 tb(Mmin,z)tU(z) = ∫dLγLγΦb(Lγ,z)∫MmindMnh(M,z), (39) beff,b(Mmin,z) = ∫MmindMbh(M,z)nh(M,z)∫MmindMnh(M,z), (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 cross-correlation 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 cross-correlation 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 large-scale matter density dominates the measured cross-correlation 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 cross-correlation 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 star-forming galaxies, radio galaxies, and diffuse processes. Joint measurement of multiple cross-correlation functions will be the key to this end. It is suggested that the cross-correlation of local galaxies and the cosmic infrared background helps to infer the star-forming 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 de-sc0018327. NY acknowledges financial support from JST CREST (JPMJCR1414). Numerical computations presented in this paper were in part carried out on the general-purpose PC farm at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. The Hyper Suprime-Cam (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 Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck 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 Harvard-Smithsonian 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. AST-1238877, 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 at http://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 Fermi-Lat Analysis Methods

The analysis of the full Fermi-LAT 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 point-like 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.

We used the Fermipy Python package 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 best-fit 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 Roma-BZCAT 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 best-fit 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 Fermi-LAT team (see, e.g., Ref. Shirasaki et al. (2016) for an example of this method). This shows that most of the detector cosmic-ray 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.

## Appendix B Cross-Correlation With Isotropic Gamma-Ray Components and Constraints on Dark Matter Annihilation

In this appendix, we present the cross-correlation analysis of the lensing convergence and the isotropic -ray background (IGRB) including the unresolved components in Fermi-LAT data.

The lensing convergence is equivalent to the cosmic surface mass density with a weight function along line-of-sight direction . The details of the reconstruction algorithm of from the HSCS16A galaxy-imaging data are given in Section II.2.2. We use the cross-correlation estimator in Eq. (20) to the IGRB measured by Fermi-LAT 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 inverse-weighed method to combine the cross-correlation signals for individual HSCS16A patches together. In addition to the normal masks as described in Section III.3, we also place circular masks around Fermi-LAT -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.

The left panels in Figure 9 summarize our cross correlation measurements of the IGRB and . We work with four different source-galaxy selections in the lensing analysis. Regardless of the choice of smoothing scales and galaxy selections, we find that the IGRB- cross-correlation is consistent with a null detection. We also confirm that the cross-correlation measurements of the IRGB and the B-mode 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 power-spectrum 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ánchez-Conde 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 cross-section , 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 Milky-Way 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:

1. Our cross-correlation analysis uses linear-scale clustering alone (see Section IV.1) that contains less information on the internal mass distribution of DM halos.

2. 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. Cross-correlation 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 large-scale structure for indirect detection of DM annihilation).

Overall, the cross-correlation measurement presented in this paper is a more promising tool to explore expected correlations of astrophysical -ray sources with the cosmic matter density.