{'fraction_non_alphanumeric': 0.03397693404172606, 'fraction_numerical': 0.013787741350265647, 'mean_word_length': 4.4755214985100045, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dark matter constraints with stacked gamma rays scales with the number of galaxies

Low-surface-brightness galaxies (LSBGs) are interesting targets for searches of dark matter emission due to their low baryonic content. However, predicting their expected dark matter emissivities is difficult because of observational challenges in their distance measurements. Here we present a stacking method that makes use of catalogs of LSBGs and maps of unresolved gamma-ray emission measured by the Fermi Gamma-Ray Space Telescope. We show that, for relatively large number of LSBGs, individual distance measurements to the LSBGs are not necessary, instead the overall distance distribution of the population is sufficient in order to impose dark matter constraints. Further, we demonstrate that the effect of the covariance between two galaxies located closely -- at an angular distance comparable to the size of the Fermi point spread function -- is negligibly small. As a case in point, we apply our pipeline to a sample of 800 faint LSBGs discovered by Hyper Suprime-Cam and find that, the 95 per cent confidence level upper limits on the dark matter annihilation cross-section scales with inverse of the number LSBGs. In light of this linear dependence with the number of objects, we argue this methodology could provide extremely powerful limits if it is applied to the more than 105 LSBGs readily available with the Legacy Survey of Space and Time.


I. INTRODUCTION
Revealing the nature of dark matter (DM) is one of the major challenges in modern cosmology and particle physics. As one of the best theoretically motivated candidates for DM, weakly interacting massive particles (WIMPs) have been considered [1] and probed via different approaches: colliders, underground experiments, and astrophysical observations [e.g. 2,3]. In the context of astrophysics, researchers often focus on probing signals from the self-annihilation or decay of WIMPs. In the early universe, WIMPs are considered to be produced in thermal equilibrium with standard model (SM) particles. As these interact with each other, WIMPs can annihilate into SM particles, particularly heavy SM particles, such as bb and W + W − [4,5]. As assumed that WIMPs are of DM abundance in the Universe, WIMP * hashimoto.daiki@f.mbox.nagoya-u.ac.jp were reported to have an annihilation cross-section of ∼ 2 × 10 −26 cm 3 /s, which is known as the thermal relic cross-section [6]. As a result of the annihilation process, γ rays can be produced directly or in secondary processes in which case more massive states decay into lighter and stable ones, particularly photons, electrons, positrons and neutrinos. Therefore, probing these particles induced by DM annihilation gives clues to specify DM properties.
The Fermi Large Area Telescope (LAT) has revealed the γ-ray sky in the energy range from about 20 MeV to approximately 1 TeV. In the fourth catalog of γ-ray point sources [7], the Fermi team presented the most detailed galactic diffuse emission model as well as a catalog containing about 5,000 objects with above 4σ significance detection. Out of all these, more than 3,700 sources have been found to be associated with pulsars, supernova remnants and/or blazars. Interestingly, Fermi -LAT observations have also been used to probe DM from different sky regions such as; the Galactic Center [8][9][10][11][12], Milky Way dwarf spheroidals (MW dSphs) [13][14][15][16][17][18], DM subhaloes [19], and the unresolved γ-ray background [20,21] (UGRB), to mention just a few. Of these targets, MW dSphs have provided the most robust DM constraint because the astrophysical background uncertainties in these objects is expected to be relatively small.
The UGRB is the residual emission obtained by subtracting the resolved point sources and galactic diffuse emission from the γ-ray observations. As aforementioned, this component has been extensively studied to probe DM annihilation signals from e.g., low-redshift galaxy catalogs, galaxy groups and galaxy clusters [22][23][24][25][26][27][28][29][30][31][32][33][34][35]. In our previous work [36], we proposed to use a lowsurface-brightness galaxy (LSBG), which has less than ∼23 mag/arcmin 2 mean surface brightness, as a new target to probe DM annihilation signals in the UGRB. In particular, LSBGs are known to be highly dominated by DM [37,38] and have more massive halos than Milky Way dSphs (∼ 10 9 M ). Furthemore, they have lower star-formation activity, γ-ray contamination from pulsars and supernova remnants, than ordinary galaxies or galaxy clusters [39]. However, owing to their faint fluxes, it is more difficult to measure precisely their redshifts than for ordinary galaxies. This explains why, in our previous study, we only used eight out of the nearly 800 LSBGs found with the Hyper Suprime-Cam (HSC) [40].
In this article, we propose a method to probe DM model parameters by using a large number of objects with unknown redshifts. Instead of individual redshifts, we focus on the redshift distribution of an entire catalog, which is obtained by the clustering redshift method. Each object's redshift is randomly assigned from the distribution in order to model the γ-ray flux induced by the DM annihilation. Using a composite likelihood analysis of the model fluxes and the UGRB data, we are able to compute constraints on the DM annihilation crosssection. In addition, we perform several validation checks of our methods using the full sample of the HSC-LSBG catalog and UGRB data from the Fermi -LAT observation. Since running a maximum-likelihood analysis with so many objects is very computationally intensive, we initially assume that the γ-ray emission from each LSBG is independent of each other. We verified the validity of this assumption by estimating the covariance of the best-fit model parameters of neighboring LSBGs. Our main analysis method is that of a stacking analysis with the population of LSBGs. We discuss the power of this technique for the problem at hand.
Because LSBGs are potentially abundant in the local universe, and the upcoming Legacy Survey of Space and Time (LSST) is expected to discover many more such systems, we anticipate that LSBGs can become competitive targets for indirect DM detection efforts. This paper is organized as follows. In Section II, we present the clustering redshift method for modeling γray emissivities of LSBGs. In Section IV, we perform a stacking analysis using co-added unresolved γ-ray maps and our sample of HSC-LSBGs. In Section V, we present a scaling relation between the 95% confidence level upper limit on the DM cross-section and the number of LSBGs.
Finally, we conclude in Section VI.

II. THEORETICAL FRAMEWORK
In this section, we present our methodology to constrain the DM annihilation cross-section using the diffuse γ-ray background and massive nearby objects having no individual distance estimates. Note that our methodology is generic and could be applied to highly DM dense regions that would appear in the Fermi data as point sources such as; dSphs, faint galaxies, and normal galaxies. As discussed in the introduction, since we do not have individual distance measurements for our targets, we require an overall probabilistic distribution of their positions along the line of sight. In addition, we assume that all objects are statistically independent, which dramatically simplifies the stacking likelihood analysis. In Subsection V A, we evaluate the validity of this assumption.
A. Predicted γ-ray flux from DM annihilations Through the DM annihilation process, γ-ray photons are generated directly or in the cascade process in which various final states (e.g. bb, W + W − and µ + µ − ) decay into more stable particles. The γ-ray flux dΦ ann /dE produced by DM annihilation can be modeled as follows, where m χ is the DM mass and σv is the velocityaveraged DM annihilation cross-section. The parameters Br i and dNi dE are the branching ratio and γ-ray energy spectrum in the i−th annihilation channel, respectively. In the analysis, we consider bb as a representative annihilation channel and obtain dN bb /dE using the DMFIT 1 [41] package included in the official fermipy analysis software. The parameter J is the so-called J-factor, which is given by the DM halo properties as follows, where s , Ω and M halo are the line-of-sight vector, solid angle and halo mass of target objects, respectively. The annihilation signal may increase if we have clumpy substructures within a halo instead of the smooth halo. This can be effectively modeled as a boost factor, b sh [42]. For clarity, we set b sh = 1 for all halos throughout this paper. The DM density profile is assumed to follow the Navarro-Frenk-White (NFW) shape [43], where c(M halo ) is the concentration parameter. The concentration parameter primarily depends on the halo mass and we compute it using the fitting formula to model the scaling relation of c(M ) with halo mass, calibrated using the high-resolution N-body simulations [44]. We convert the halo mass to the concentration by using the COLOSSUS package [45]. The full description of J-factor estimation can be found in [36]. A brief summary is as follows. We first convert the observed magnitude into the absolute magnitude. The V-band apparent magnitude can be converted from the gri system as V = g − 0.59(g − r) − 0.01 [46]. Now, the situation is that we do not have the distance to individual galaxies but an overall distribution. However, we can formally write the distance of each galaxy as a random draw from the distribution: d ∈ dN/dz(z). For the particular realization of the random draw of the distance with the binding condition of d = dN/dz, the absolute magnitude can be computed; thus, the halo mass can be derived by assuming the relations of mass-to-light ratio [47] and stellar to halo mass ratio [48].
Given that we have the overall dN/dz distribution, with measurement errors, we can simulate the effect of neglecting distance to individual galaxies. Figure 4 shows the distribution of the total J-factor after stacking N st galaxies. For this plot, we use the HSC-Y1 LSBG sample, where no distance is available. The N st LSBGs are randomly selected from the parent sample 500 times. The scatter in the figure corresponds to the 95% range due to the effect of the sample variance. In this subsection, we describe a method to estimate the dN/dz distribution from spatial clustering; the socalled clustering redshift method. The application of this method to the data is shown in Subsection IV B.
Given two galaxy samples in an overlapped region, even if the galaxies in the two galaxy samples are statistically different, they both correlate with the underlying DM distribution. In the case where one galaxy sample has known redshifts and the other does not, it is possible to take the cross-correlation between the two galaxy samples to obtain a statistical estimate of the redshift distribution of the galaxies with unknown redshifts. We denote the galaxy samples with known and unknown redshifts as 'spectroscopic (spec-z) sample' and 'photometric sample', respectively. The angular cross-correlation function can be factorized as follows [49], where subscripts s and p represent spec-z and photometric samples, respectively. dN/dz represents the redshift distribution normalized to unity, b(z) is a linear bias, and w DM (z, θ) is the angular correlation function of DM. Note that the mass of the employed DM model is high enough so as not to affect the DM clustering pattern itself. We define an integrated cross-correlationw with a weighting function W (θ) as where the weight W is introduced so that the signal-tonoise ratio ofw can be optimized. Following [50], we empirically adopt W = θ −1 . In practice, we divide the spec-z sample into narrow redshift bins so that dN s /dz can be approximated by the narrow top-hat function; wherew DM (z) can be defined similarly to Equation 5 by replacing w(θ) with w DM (z, θ).
Here, we assume that the biases of the LSBG and specz sample are constant over redshift because the redshift range that we consider is narrow. Althoughw DM (z) is fully derived from the standard cold DM theory including the non-linear matter clustering evolution, we do not need to estimate the absolute value of it. What we need to consider onw DM is the redshift evolution. Given the narrow redshift range, in our analysis, we simply replacew DM (z) with the square of the linear growth factor, D 2 (z).

III. DATA
We start by validating our methodology. For this, we apply the aforementioned pipeline to the LSBG samples [40] observed with the HSC survey. The 781 LSBG samples in ∼ 200 deg 2 of the HSC region do not have individual distance measure, but we have the spec-z samples from NASA Sloan Atlas (NSA) [51] on an almost fully overlapped sky area. Lastly, we use the Fermi -LAT unresolved extragalactic γ-ray background by applying the composite likelihood analysis to the sample of HSC-LSBGs. Below we provide details of the datasets used in our study.

A. HSC-LSBG catalog
HSC is a wide-field camera, attached to the prime focus of the Subaru telescope, covering ∼1.5 • diameter field of view with 0.17 arcsec pixel scale [52,53]. The HSC Subaru Strategic Program survey consists of three layers, wide, deep and ultradeep layers; the wide layer has five broad photometric bands g, r, i, z and y. As described in a report of the second data release of the survey by [54], the wide layer has a depth of 24.5-26.6 in the 5 filters for 5 σ point-source detection. In the final data release, the survey will cover 1400 deg 2 sky in a depth of i ∼26 mag.
Since the HSC pipeline, hscpipe, is not optimized for detecting and measuring diffuse objects, HSC images are reduced first by the hscpipe and then SExtractor is used for detection and measurements. [40] processed an HSC dataset with three broad bands (g, r and i) on a patch-by-patch basis over ∼ 200 deg 2 and produced a catalog comprising 781 LSB objects, which have a mean surface brightness in g−band > 24.3mag/arcsec 2 . Briefly, the following processes have been executed: (i) bright sources and associated diffuse lights are subtracted from images to avoid contamination to LSB objects detection; (ii) after the Gaussian smoothing with full depth at half maximum of 1 , sources with the half-light radius r 1/2 satisfying 2.5 < r 1/2 < 20 are extracted; (iii) the sources are selected by applying reasonable color cuts to remove optical artifacts and distant galaxies; (iv) by modeling the surface brightness profiles of LSBG candidates, astronomical false positive are removed; (v) by visual inspection, false candidates such as point-sources with diffuse background lights are removed and then 781 LSBGs are finally left.
To minimize a possible contamination of astrophysical γ-ray photons,it might be useful to restrict our analysis to quiescent galaxies, because such galaxies do not have ongoing star-formation activities and therefore unlikely contain high-energy astrophysical sources such as supernova remnants and AGNs, at least compared to starforming blue galaxies. We divide the sample into red and blue LSBGs, where red is defined by g − i ≥ 0.64 and blue g − i < 0.64, which include 450 and 331 objects, respectively. This color selection roughly corresponds to the galaxy age of 1 Gyr for a 0.4 × solar metallicity galaxy [40]. In Figure 1, we show the sky distribution of our LSBG sample in green dots. For a random catalog corresponding to the LSBG catalog, we employ the random catalog of the HSC photometric data, and randomly resample it such that the number density is roughly 10 times larger than the LSBG density. We also apply the bright star mask to the random catalog.

B. NSA sample
To measure the dN/dz distribution, we need reference spec-z samples in the overlapped region in both sky-coverage and redshift range. Since the LSBGs have been detected in the local universe by radio observations [38,55], optical selected LSBGs are also expected to be in the local universe; thus, we need a low-redshift spec-z sample. The NSA sample 2 is a spec-z sample obtained from the spec-z campaign of the Sloan Digital Sky Survey with the Galaxy Evolution Explorer data for the energy spectrum of the ultraviolet wavelength and includes objects up to z = 0.15 including 11,820 objects in the overlap region of the LSBG. We show the sky position and redshift distribution in Figure 1. Since the uniformity of the NSA sample is not guaranteed, we attempt to mitigate the non-uniformity as follows. First, we remove the sample from both HSC and NSA in the bright star masked regions. We checked that, after removing the masked regions, the local number density of NSA galaxies smoothed with each HSC patch has a uniform distribution in most areas of the HSC regions we are working on, except for the low-density region in the VIMOS-VLT Deep Survey (Dec.>1). We generate random catalogs including about 10 times more objects than that of the NASA galaxies. In addition, we removed the edge regions in the HSC survey footprint for safety, because the exact survey window near the boundaries is difficult to define. C. UGRB from LAT data

LAT data
We analyze the Fermi -LAT Pass 8 data obtained from 2008-08-04 to 2016-08-02, which contains several upgrades to account for a much improved knowledge of the instrument response function [56]. As recommended by the Fermi collaboration 3 , we select the photon event class P8R3 SOURCE [57], which is suitable for point-source analyses. In addition, we apply the following selections, DATA QUAL>0, LAT CONFIG==1 and P8R3 SOURCE V2, as the corresponding filter expression and the instrument response function for the event class, respectively. We use photon counts in the energy range 500 MeV to 500 GeV, where we split the data into 24 logarithmically spaced energy bins, and spatial bins of size 0.1 • × 0.1 • . We selected the width of the energy bins such that those are larger than the energy dispersion of the instrument at any given bin, but we require spatial bins to be smaller than the LAT PSF at all energies. The lower energy scale is determined to balance two opposing effects. If we decrease the low-energy limit, photons around bright sources would leak due to the broadening of the PSF, but if we increase the low-energy limit, the photon count statistics decreases.
For the composite analysis described in Subsection IV A, we select 29 patches whose centers are located in the HSC regions with the separation of at least 3 • from each other and individual patches having 10 • × 10 • sky coverage. Moreover, to avoid contaminating the γray photons produced by the Earth's atmosphere interaction with high energy cosmic rays, we exclude the photon data with zenith angles greater than 100 • . In the LAT data analysis, we use fermipy (v1.0.1) 4 [58], which is an open-source software package based on the Fermi Science Tools (v2.0.8) 5

UGRB construction and putative flux
Our procedure to reconstruct the UGRB observations and compute the likelihood profile for each LSBG object is the same one introduced in [13].
First, we perform the maximum likelihood analysis in each patch to optimize the normalization and spectral parameters of all the standard 4FGL [7] sources in the patch -including the galactic diffuse emission, isotropic emission, and extended sources.
Specifically, we adopt the standard Galactic diffuse emission model (gll iem v07.fits) and the isotropic emission template (iso P8R3 SOURCE V2 v01.txt) 6 , respectively. The isotropic template represents isotropic contributions from undetected extragalactic sources and the residual cosmic ray background.
In Figure 2, we show an example of the observed and residual counts for one of our patches. In addition, we show a histogram of the relative fluctuation of the residual counts map for the same patch. The rms of the fluctuation being is ∼ 0.04. Note that we confirm that our patches do not contain new point sources with test statistics (TS) larger than 25.
Once reconstructed the UGRB, we estimate the putative flux of each LSBG. According to the prescription described in the 2FGL catalog, the Bayesian method [59] should be applied for the likelihood analysis of very faint sources [60]. We introduce free parameters α ij representing the i-th point-source flux-amplitude in j-th energy bin. Here we define TS as T S ≡ −2∆ log L and ∆ log L is given by; where D and Θ are the LAT data and the nuisance parameters of all the γ-ray sources other than our target LSBGs. The best fit parameters α max i,j are obtained when log L(D, Θ|α i,j ) is maximized. In the above equation α i,j = 0 means no target source. We checked that the TS values of all our LSBGs are less than 1 for most energy bins, which justifies our choice of using the Bayesian method introduced in [7].

IV. RESULT
A. Composite analysis As described in Section III C 2, because of the very low γ-ray signal of our LSBGs, we compute the 95% C.L. upper limits on the DM annihilation cross-section parameter using the same approach as in [36]. Since all likelihood values obtained at each object are assumed to be independent of each other, the composite likelihood L st for the full sample of targets can be expressed as follows; where J i is the J-factor of the i-th target and the index j runs over all energy bins. L ann i,j is the log-likelihood value for the i-th LSBG, which is obtained by the difference of the log-likelihood value from the one in the case of no source.
Now we calculate the J-factor for each object to evaluate the model flux described in Section II A. For the simplicity of Equation 2, we need to consider the relationship between the PSF of the LAT instrument and the angular size of objects. The PSF (68% containment angles) decreases from ∼ 1.5 • to ∼ 0.1 • as γ-ray energy increases from 500 MeV to 500 GeV. The angular size of LSBG is smaller than a few 10 arcsecs, hence smaller than the PSF in all energy bands. Therefore, we can consider them as point-like sources in the likelihood computation. Then, the integration of ρ 2 DM over the target volume in Equation 2 is reduced to; where d A is the angular diameter distance to the object, which is given by assignment of distance randomly drawn followed the measured dN/dz distribution. Then, Equation 9 is straightforwardly calculated; where ∆ is the spherical overdensity set to be 200 and ρ c,z is the critical density at the redshift. Since our targets are regarded as point sources, our assumption is correct if the angular separations between objects are larger than the LAT PSF over all considered energy ranges; otherwise, it is incorrect because their mean surface number density is about 4 per deg 2 . We will further discuss the parameter correlation between neighboring objects in Section V A.
In our procedure, we assume that the LSBG flux is positive definite, which implies that the data are welldescribed by a χ 2 /2 distribution rather than χ 2 . As such, the 95% C.L. upper limits on the cross section σv UL are given when the ∆ log L st (α i,j | σv , J) ∼ −3.8/2. B. dN/dz measurement Following the methodology described in Subsection II B, we estimate the dN/dz distribution of the HSC-LSBG sample. For this, we divide the reference redshift sample into five equally-separated bins from z = 0 to 0.15. Then, we compute the cross correlation between the sample in each bin and the entire LSBG sample. The angular cross correlation is computed in the angular bin 0.1 • < θ < 1.0 • , which is further divided in smaller logarithmic-scaled bins. We use the estimator [61], where DD or DR represent the normalized number of pairs separated within the i−th angular bin between data and data or data and random, respectively (in our case, "data" and "random" represent the spec-z or LSBG samples in our datasets and the corresponding random samples, respectively). Subscripts p, s, and m represent the photometric sample, and reference sample in the m−th redshift bin, respectively.
The rest of the subsection is devoted to describing the covariance of the estimation, which will be used to draw the galaxy redshift randomly from the dN/dz distribution. We estimate the covariance of w(z, θ) using the Jackknife method [62] as where w m ik is the angular correlation function for the k−th Jackknife sub-sample in the i-th angular and m-th redshift bins. M is the number of the Jackknife sub-samples, and we take M = 100 which divides the entire region into ∼ 2 deg 2 patches.ŵ is the averaged correlation function over all jackknife sub-samples, Now the dN/dz at the m-th redshift bin can be obtained by integratingŵ m i along the angular scale, and its amplitude can be considered as a parameter A m . Therefore, the dN/dz distribution can be stochastically determined with the likelihood function of A, which is assumed to be Gaussian, possibly with a physically reasonable prior P (A < 0) = 0 or 1 otherwise. The variance of the Gaussian likelihood function can be derived aŝ wherew m can be computed using Equation 5. Figure 3 represents the dN/dz distribution measurement. The errors are measured using the Jackknife resampling. For the clustering analysis described in this section, we use a publicly available code, treecorr [63].

C. Evaluation of statistical uncertainty
The upper limit on the cross section in the composite analysis can be affected by the dN/dz measurement error and the halo property estimation, which is specified by the halo mass and concentration parameter as described in Section IV A. We evaluate the coverage of the upper limits by using 500 Monte Carlo simulations. First, for the halo mass, we evaluate its uncertainty as ∆ log M halo = 0.8 at 1-σ Gaussian error by computing the scatter of stellar-to-halo-mass conversion. In addition, we adopt a concentration parameter error of ∆ log c of 0.1 at 1-σ Gaussian error [64]. Accordingly, the total uncertainty for halo properties results in J-factor uncertainties of ∼0.9 dex at 1-σ error. Moreover, we randomly assign the distance to galaxies according to this distribution; therefore, the negative values of the amplitude are physically unreasonable. Therefore, as described in the previous section, we obtain the posterior function of A m as P (A m |D) = P (D|A m )|P (A m ). We finally apply linear interpolation between each center of redshift bin to the posterior distribution. We take a conservative limit of the minimum redshift of the sample corresponding to 25 Mpc, which is the minimum distance among the HSC-LSBG samples with precisely measured distance.
In Figure 4, we show the total J-factor values as a function of the number of stacked objects, N st . In the order of square, cross and circle symbols, the error-bars plot the total values including dN/dz measurement uncertainty, halo property and both, respectively. For comparison with other studies, we display σv UL at the 95% C.L. for a DM mass of 1 TeV in the right axis, and their corresponding J-factor value on the left axis. Note that when converting the J-factor to σv UL , we apply a mean flux of our UGRB sky. The upper limit is affected by both the fluctuation and target's J-factor value, however we find that, even in the lower energy regime, the scatter by the fluctuation is smaller (∼0.4 dex at 2-σ level) than the total halo property and dN/dz measurement uncertainty.

A. Correlation between neighbors
We performed the composite analysis in Section IV A assuming that the likelihood functions of individual objects are independent of each other. Such correlations are expected when the number density of the sample is high because the PSF of Fermi -LAT is ∼ 1 • . In this subsection, we demonstrate that the correlation between data at different points can be negligible.
In the HSC-Fermi sky coverage, we select 10 independent patches with the size of 10 × 10 deg 2 . From each patch, we randomly select 60 pairs of points with separations of 0.5 • to 3 • . To evaluate the correlation between the putative flux amplitudes of the paired objects, we perform the joint likelihood analysis for the pairs, which simultaneously optimizes the fluxes of the paired objects. The putative flux of an object is defined by, FIG. 5. The correlation between two amplitude parameters for each paired objects in the overall energy ranges, as a function of separation angles. The error bars represent 1-σ errors of the correlations derived from 60 independent pairs in each angular bin.
We fit the amplitude parameters to the UGRB flux and measure the covariance matrix using the fit and sed methods in fermipy. Figure 5 shows the absolute value of the correlation between two amplitude parameters in the overall energy ranges, as the function of the separation. The error-bars are computed from the 60 independent pairs that reflects the fluctuations of the residual gamma-ray flux. The cross-covariance is normalized by the diagonal terms, i.e., ρ ij ≡ cov(α i , α j )/σ i σ j . Although we expect a strong correlation on scales smaller than 1 deg, the correlation at the smallest separation, which corresponds to the HSC-LSBG mean separation, is less than 0.1 deg at 1-σ level. We note that all the cross correlation values are negative at all scales, which is the consequence of the conservation of the total flux. For further validation, we perform a composite likelihood analysis in which we obtain the likelihood profiles for the putative fluxes of all samples within a single LAT data patch simultaneously. Figure 6 compares the σv UL constraints with this simultaneous approach ('simultaneous' case) with the one obtained based on the assumption that all objects are independent of each other ('independent' case). We emphasize that the 'independent' case gives a weaker constraint on σv UL than the 'simultaneous' case. This is because the total flux conservation is imposed, which results in a larger putative flux amplitude for the 'independent' case. Note that in this calculation, we set all objects' J-factor to 10 14.5 GeV 2 /cm 5 , and choose a specific dN/dz. This explains why the amplitudes of σv UL in Figure 6 do not correspond to the ones in Figure 4.
We conclude that the correlation between neighboring points is less than 10%, on scales of order 0.5 deg. Even if we ignore this correlation, which is computationally much less expensive, we would obtain conservative constraints on the σv UL .

B. The power of stacking
In this subsection, we discuss a scaling relation of the statistical power on σv UL as a function of the number of objects included in the stacking procedure. We first analytically show that σv UL is proportional to the inverse of N st in the limit of zero observed photon counts. We then also demonstrate that this scaling relation converges to N st independent of the DM mass if the N st is sufficiently large by using HSC-LSBG catalog and Fermi observations.
Let us first revisit the Poisson likelihood for a single LSBG in the UGRB. The likelihood function for the model parameters is given by the Poisson distribution, where the index i runs over all energy bins and pixels and n i is the observed photon counts; λ i is expected photon counts for model fluxes, which is decomposed into the Galactic foreground, isotropic background and resolved point-source model fluxes as well as the flux by the DM annihilation for the single LSBG. A UGRB sky is derived from the modeling of all the γ-ray sources except for the LSBG flux as described in Section III C 2. Therefore, if we fix the model parameters Θ of the background sky, then λ i depends only on the parameter σv . We denote i is the total model flux of all sources without LSBG and λ T i is the model flux of LSBG. For simplicity, we consider that all LSBGs have the same J-factor, which means that their model fluxes are equal. In addition, we consider the fact that, in energy regimes higher than ∼30 GeV, the LAT hardly detects photons. Given that there is no photon count (n i = 0) in such energy regimes and in all pixels, we expect λ others i (µ) = 0 and then log L = − i λ T i ( σv ).
Consequently, we obtain the composite likelihood with N st objects using Equation 8, where N st is the number of objects in the composite analysis. Thus, in our approach the upper limit is proportional to 1/N st since λ T i ( σv ) ∝ σv . Further, we explore the scaling relation using the HSC-LSBG catalog and Fermi data. We calculate the constraint on σv from randomly drawn N st objects from the HSC-LSBG sample. Because here we aim to isolate the effect of number of stacking from the individual halo properties, we assume the selected sample has a same J factor, but using the actual background γ-ray fluxes at that position. We repeat this calculation 500 times for each N st sampling. In Figure 7, we show the median value of the ratio of σv UL with N st objects to the one with a single object for DM mass of 10 GeV, 100 GeV and 1 TeV. With N st larger than ∼30, the upper limits scale with the inverse of N st for all mass ranges. We have discussed this relation in our previous paper [36], which predicted the relation based only on 8 HSC LSBGs . In this paper, we confirm that prediction of scaling with N st , but furthermore we find, thanks to the large number of LSBG sample, that the relation is independent of DM mass. This is because most of the constraints are dominated by the objects having significantly low background fluxes. Such objects can be found with a certain probability, say ∼ 10% of the entire sample. Therefore, when the number of stacking objects is small, e.g. N st < 10, the probability of having strongly constrained objects is less than 10%, and the constraining power does not scale with 1/N st . However, the probability of having such objects may increase and converge to 10%. as we increase the number of stacking. Therefore, the statistical power of stacking, dominated by the objects with almost zerobackground fluxes scales with 1/N st . For DM mass of 100 GeV and 1 TeV, this scaling relation is seen, even in smaller N st . This behavior is reasonable, considering the photon-count statistics in a high energy regime. The annihilation process with more massive DM particles can produce higher energy photons, thus probing the DM annihilation for more massive DM is affected by photoncount statistics in high energy regimes.

VI. CONCLUSIONS
In this study, we proposed a new method to constrain the DM annihilation cross-section with potential γ-ray sources with unknown redshifts and the observed γ-ray background. The crucial point is that we do not need to know the individual distance to the γ-ray sources. The overall redshift distribution is sufficient to constrain the DM annihilation cross section. By applying the redshift clustering method, we obtained a redshift distribution of the whole population, instead of the redshifts of individual objects. We then randomly assigned the distance to each LSBG according to this distribution. To validate the proposed method, we measured dN/dz of full samples of the HSC-LSBG catalog including ∼800 objects and constrained the DM annihilation cross-section via a composite likelihood analysis using all LSBGs given distances and using the eight years of Fermi -LAT observations.
In the dN/dz measurement for the LSBG samples, we employed the NSA catalog as spec-z samples in overlapped regions with the HSC region. The uncertainty in the dN/dz was evaluated with Monte Carlo simulations. We found that in the composite analysis with ∼1000 objects, uncertainties of the halo property estimation and dN/dz measurement are quite small and not significantly affected to the upper limit on the cross section. Therefore, using dN/dz distribution of overall target samples instead of the individual sample's redshifts is a valid and robust approach for the DM signal search.
To reduce the computational cost, we assumed that all the spectral parameters of the different LSBGs are independent of each other. The number density of the LSBGs (∼ 0.5 • ) was comparable to the LAT PSF scale, which may lead to a correlation between spectral parameters for neighboring LSBGs, and thus may break the assumption. We validated this conjecture in two ways. First, we computed the covariance matrix between the two normalizations for each paired neighbor and checked that off-diagonals terms of the correlation matrix is smaller than 0.1 at the separation angle 0.5 • . Second, fixing the J-factor and distance of each object, we performed a composite likelihood analysis in which we obtained the likelihood profiles for putative fluxes of all samples in each LAT patch simultaneously. By comparing σv UL , we found that the constraints differ at most by a factor of 2. Moreover, if we assume that the data is indepen-dent, then the constraints become more conservative due to the relaxation of the total flux conservation condition.
Finally, we found the scaling relation of σv UL with the number of objects N st in the composite likelihood analysis. First, under an assumption of zero observed photon, we discussed the relation analytically and found that σv UL is proportional to 1/N st .Furthermore, using the real HSC-LSBG catalog and the LAT data, we computed the scaling for DM mass of 10 GeV, 100 GeV and 1 TeV and showed that for all DM masses the upper limit scales with 1/N st using N st 30. This is a significant effect which is the consequence of the Poisson statistics due to rare photons detected by the LAT and is different from the Gaussian statistics in which the scaling obeys 1/ √ N st . In future imaging surveys with next-generation telescopes such as the LSST, a large number of LS-BGs will be detected because of wider sky coverage and better sensitivity. For example, LSST has a sky coverage of ∼20,000 deg 2 and reaches the depth of ∼27.5 mag/arcsec 2 in i band and consequently has the potential to discover O(10 5 ) objects. Although, for the reference samples we need spec-z or high-precision photoz samples located in the local universe, it is expected that these will decrease the uncertainties in the estimate of σv UL due to the dN/dz measurement and the halo properties. Moreover, by increasing the statistics, we can perform a more detailed dN/dz measurement, particularly in a redshift range corresponding to a distance of less than 25 Mpc, which is the minimum distance adopted here for a conservative J-factor estimate in our likelihood procedure. We expect that in the future, this will give a more stringent constraint on the DM cross-section.