Dark matter constraints with stacked gamma rays scales with the number of galaxies

Through the DM annihilation process, $\gamma$-ray photons are generated directly or in the cascade process in which various final states (e.g. $b\bar{b}$, $W^{+}W^{-}$ and $\mu^{+}\mu^{-}$) decay into more stable particles. The $\gamma$-ray flux ${d\Phi_{\rm ann}}/{dE}$ produced by DM annihilation can be modeled as follows,

where $m_{\chi}$ is the DM mass and $\langle\sigma v\rangle$ is the velocity-averaged DM annihilation cross-section. The parameters ${\rm Br}_{i}$ and $\frac{dN_{i}}{dE}$ are the branching ratio and $\gamma$-ray energy spectrum in the $i-$th annihilation channel, respectively. In the analysis, we consider $b\bar{b}$ as a representative annihilation channel and obtain $dN_{b\bar{b}}/dE$ using the DMFIT¹¹1https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/gammamc_dif.dat [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^{\prime}$, $\Omega^{\prime}$ and $M_{\rm 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_{\rm sh}$ [42]. For clarity, we set $b_{\rm 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_{\rm 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\in dN/dz(z)$. For the particular realization of the random draw of the distance with the binding condition of $\langle d\rangle=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_{\rm st}$ galaxies. For this plot, we use the HSC-Y1 LSBG sample, where no distance is available. The $N_{\rm 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 so-called clustering redshift method. The application of this method to the data is shown in Subsection IV.2.

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_{\rm DM}(z,\theta)$ 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-correlation $\bar{w}$ with a weighting function $W(\theta)$ as

where the weight $W$ is introduced so that the signal-to-noise ratio of $\bar{w}$ can be optimized. Following [50], we empirically adopt $W=\theta^{-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; $dN_{s}^{i}/dz\simeq 1/\Delta z$ if $z_{i}<z<z_{i+1}$. After this, we rewrite Equation 5 at $z=z_{i}$ as

where $\bar{w}_{\rm DM}(z)$ can be defined similarly to Equation 5 by replacing $w(\theta)$ with $w_{\rm DM}(z,\theta)$.

Here, we assume that the biases of the LSBG and spec-z sample are constant over redshift because the redshift range that we consider is narrow. Although $\bar{w}_{\rm 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 on $\bar{w}_{\rm DM}$ is the redshift evolution. Given the narrow redshift range, in our analysis, we simply replace $\bar{w}_{\rm DM}(z)$ with the square of the linear growth factor, $D^{2}(z)$.

HSC is a wide-field camera, attached to the prime focus of the Subaru telescope, covering $\sim$$1.5^{\circ}$ 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 $\sigma$ point-source detection. In the final data release, the survey will cover 1400 $\rm deg^{2}$ sky in a depth of $i\sim$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 $\sim 200~{}\rm deg^{2}$ and produced a catalog comprising 781 LSB objects, which have a mean surface brightness in $g-$band $>24.3\rm mag/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^{\prime\prime}$, sources with the half-light radius $r_{1/2}$ satisfying $2.5^{\prime\prime}<r_{1/2}<20^{\prime\prime}$ 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 $\gamma$-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 star-forming blue galaxies. We divide the sample into red and blue LSBGs, where red is defined by $g-i\geq 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 $\times$ 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.

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 ²²2https://data.sdss.org/sas/dr13/sdss/atlas/v1/nsa_v1_0_1.fits 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.

As described in Section III.3.2, because of the very low $\gamma$-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 ${\mathcal{L}}_{\rm 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. ${\cal L}^{\rm 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.1. 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 $\sim 1.5^{\circ}$ to $\sim 0.1^{\circ}$ as $\gamma$-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 $\rho_{\rm DM}^{2}$ 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 $\Delta$ is the spherical overdensity set to be 200 and $\rho_{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 $\rm deg^{2}$. We will further discuss the parameter correlation between neighboring objects in Section V.1. In our procedure, we assume that the LSBG flux is positive definite, which implies that the data are well-described by a $\chi^{2}/2$ distribution rather than $\chi^{2}$. As such, the 95% C.L. upper limits on the cross section $\langle\sigma v\rangle_{\rm UL}$ are given when the $\Delta\log\mathcal{L}_{\rm st}(\alpha_{i,j}|\langle\sigma v\rangle,J)\sim-3.8/2$.

Following the methodology described in Subsection II.2, 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^{\circ}<\theta<1.0^{\circ}$, which is further divided in smaller logarithmic-scaled bins. We use the estimator [61],

where ${\rm D}{\rm D}$ or ${\rm D}{\rm R}$ 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,\theta)$ using the Jackknife method [62] as

where $w_{ik}^{m}$ 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 $\sim 2~{}$deg² patches. $\hat{w}$ 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 $\hat{w}^{m}_{i}$ along the angular scale, and its amplitude can be considered as a parameter ${\cal A}_{m}$. Therefore, the $dN/dz$ distribution can be stochastically determined with the likelihood function of ${\cal A}$, which is assumed to be Gaussian, possibly with a physically reasonable prior $P({\cal A}<0)=0$ or $1$ otherwise. The variance of the Gaussian likelihood function can be derived as

where $\bar{w}^{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].

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.1. We evaluate the coverage of the upper limits by using 500 Monte Carlo simulations. First, for the halo mass, we evaluate its uncertainty as $\Delta\log{M_{\rm halo}}=0.8$ at 1-$\sigma$ Gaussian error by computing the scatter of stellar-to-halo-mass conversion. In addition, we adopt a concentration parameter error of $\Delta\log{c}$ of $0.1$ at 1-$\sigma$ Gaussian error [64]. Accordingly, the total uncertainty for halo properties results in $J$-factor uncertainties of $\sim$0.9 dex at 1-$\sigma$ 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 ${\cal A}_{m}$ as $P({\cal A}_{m}|{\cal D})=P({\cal D}|{\cal A}_{m})|P({\cal 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_{\rm 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 $\langle\sigma v\rangle_{\rm 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 $\langle\sigma v\rangle_{\rm 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 ($\sim$0.4 dex at 2-$\sigma$ level) than the total halo property and $dN/dz$ measurement uncertainty.

We performed the composite analysis in Section IV.1 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 $\sim 1^{\circ}$. 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\times 10$ deg². 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,

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., $\rho_{ij}\equiv{\rm cov}(\alpha_{i},\alpha_{j})/\sigma_{i}\sigma_{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-$\sigma$ 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 $\langle\sigma v\rangle_{\rm 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 $\langle\sigma v\rangle_{\rm 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}$ $\rm GeV^{2}/cm^{5}$, and choose a specific $dN/dz$. This explains why the amplitudes of $\langle\sigma v\rangle_{\rm 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 $\langle\sigma v\rangle_{\rm UL}$.

In this subsection, we discuss a scaling relation of the statistical power on $\langle\sigma v\rangle_{\rm UL}$ as a function of the number of objects included in the stacking procedure. We first analytically show that $\langle\sigma v\rangle_{\rm UL}$ is proportional to the inverse of $N_{\rm st}$ in the limit of zero observed photon counts. We then also demonstrate that this scaling relation converges to $N_{\rm st}$ independent of the DM mass if the $N_{\rm 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; $\lambda_{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 $\gamma$-ray sources except for the LSBG flux as described in Section III.3.2. Therefore, if we fix the model parameters $\bm{\Theta}$ of the background sky, then $\lambda_{i}$ depends only on the parameter $\langle\sigma v\rangle$. We denote $\lambda_{i}=\lambda_{i}^{\rm others}(\bm{\Theta})+\lambda_{i}^{T}(\langle\sigma v\rangle)$, where $\lambda_{i}^{\rm others}$ is the total model flux of all sources without LSBG and $\lambda_{i}^{T}$ 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 $\sim$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 $\lambda_{i}^{\rm others}(\mu)=0$ and then $\log\mathcal{L}=-\sum_{i}\lambda_{i}^{T}(\langle\sigma v\rangle)$. Consequently, we obtain the composite likelihood with $N_{\rm st}$ objects using Equation 8,

where $N_{\rm st}$ is the number of objects in the composite analysis. Thus, in our approach the upper limit is proportional to $1/N_{\rm st}$ since $\lambda_{i}^{T}(\langle\sigma v\rangle)\propto\langle\sigma v\rangle$.

Further, we explore the scaling relation using the HSC-LSBG catalog and Fermi data. We calculate the constraint on $\langle\sigma v\rangle$ from randomly drawn $N_{\rm 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 $\gamma$-ray fluxes at that position. We repeat this calculation 500 times for each $N_{\rm st}$ sampling. In Figure 7, we show the median value of the ratio of $\langle\sigma v\rangle_{\rm UL}$ with $N_{\rm st}$ objects to the one with a single object for DM mass of 10 GeV, 100 GeV and 1 TeV. With $N_{\rm st}$ larger than $\sim$30, the upper limits scale with the inverse of $N_{\rm 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_{\rm 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 $\sim 10\%$ of the entire sample. Therefore, when the number of stacking objects is small, e.g. $N_{\rm st}<10$, the probability of having strongly constrained objects is less than $10\%$, and the constraining power does not scale with $1/N_{\rm 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 zero-background fluxes scales with $1/N_{\rm st}$. For DM mass of 100 GeV and 1 TeV, this scaling relation is seen, even in smaller $N_{\rm 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 photon-count statistics in high energy regimes.