Constraining dark matter annihilation with HSC Low Surface Brightness Galaxies

TL;DR

This study explores low surface brightness galaxies as new targets for dark matter detection, analyzing gamma-ray data to set constraints on dark matter annihilation, and forecasts improved sensitivity with future detections.

Contribution

It introduces LSBGs as novel targets for dark matter searches and provides initial gamma-ray constraints using Subaru and Fermi data.

Findings

No gamma-ray excess detected in 8 LSBGs.

Constraints set on dark matter annihilation cross-section at 10^-23 cm^3/s.

Future LSBG detections could improve sensitivity to 3*10^-25 cm^3/s.

Abstract

Searches for dark matter annihilation signals have been carried out in a number of target regions such as the Galactic Center and Milky Way dwarf spheroidal galaxies (dSphs), among a few others. Here we propose low surface brightness galaxies (LSBGs) asnovel targets for the indirect detection of dark matter emission. In particular, LSBGs are known to have very large dark matter contents and be less contaminated by extragalactic gamma-ray sources (e.g., blazars) compared to star forming galaxies. We report on an analysis that uses eight LSBGs (detected by Subaru Hyper Suprime-Cam survey data) with known redshifts to conduct a search for gamma-ray emission at the positions of these new objects in Fermi Large Area Telescope data. We found no excesses of gamma-ray emission and set constraints on the dark matter annihilation cross-section. We exclude (at the 95% C.L.) dark matter scenarios predicting a cross-section higher than 10^-23[cm^3/s] for dark matter particles of mass 10 GeV self-annihilating in the b_b channel. Although this constraint is weaker than the ones reported in recent studies using other targets, we note that in the near future, the number of detections of new LSBGs will increase by a few orders of magnitude. We forecast that with the use of the full catalog of soon-to-be-detected LSBGs the constraint will reach cross-section sensitivities of ~ 3*10^-25 [cm^3/s] for dark matter particles with masses less than 10 GeV.

Figures

Tables

Table 1. Table 1 : Parameters of LSBGs in our analysis. Each parameter represents object ID, Galactic longitude, Galactic latitude, redshifts, i 𝑖 i -band magnitudes, errors of i 𝑖 i -band magnitudes, color diagram for g − r 𝑔 𝑟 g-r , g − i 𝑔 𝑖 g-i and comoving distances of LSBGs, respectively.

LSBG ID	($l$ [deg], $b$ [deg])	redshift	$m_i$	$m_{i,\mathrm{err}}$	$g-r$	$g-i$	Distance[Mpc]
171	(62.628, -45.915)	0.0439	17.55	0.24	0.49	0.77	186
285	(178.250, -57.202)	0.00581	17.47	0.24	0.44	0.6	24.9
456	(351.210, 54.493)	0.0286	17.18	0.24	0.4	0.57	122
464	(348.724, 55.429)	0.0257	16.95	0.24	0.37	0.58	110
575	(224.099, 24.123)	0.00695	18.45	0.24	0.29	0.3	29.7
613	(339.731, 57.465)	0.0244	19.17	0.24	0.19	0.29	104
729	(336.533, 56.860)	0.0251	17.57	0.24	0.36	0.55	107
750	(276.818, 59.451)	0.00862	18.38	0.24	0.24	0.33	36.9

Table 2. Table 2 : Stellar mass, halo mass and J 𝐽 J -factor value for eight LSBGs. Each parameter is a median value in the Monte Carlo simulation. Errors of stellar mass, halo mass and J 𝐽 J -factor are ∼ 0.2 , ∼ 0.4 \sim 0.2,\sim 0.4 and ∼ 0.7 similar-to absent 0.7 \sim 0.7 in logarithmic, respectively.

LSBG ID	${\log }_{10}\left(M_{\ast }\left[M_{\odot }\right]\right)$	${\log }_{10}(M_{\mathrm{halo}}$$\left[M_{\odot }\right]$)	${\log }_{10}(J\left[{\mathrm{GeV}}^2{\mathrm{cm}}^{-5}\right]$)
171	9.3	11.1	14.4
285	7.6	10.4	15.7
456	9.1	11.1	14.7
464	9.1	11.0	14.8
575	7.4	10.3	15.5
613	9.3	10.7	14.7
729	9.3	10.7	14.8
750	7.6	10.4	15.4

Equations

$$\frac{d\Phi_{\gamma}}{dE_{\gamma}}=J\times\frac{\langle\sigma v\rangle}{8\pi m_{\chi}^{2}}\sum_{i}Br_{i}\frac{dN_{i}}{dE_{\gamma}^{\prime}}\Bigg{|}_{E_{\gamma}^{\prime}=(1+z)E_{\gamma}},$$

$$V=g-0.59(g-r)-\mathrm{0.01.}$$

$$M_{\ast }=\frac{M_{\odot }}{L_{\odot }}{L}_{\mathrm{LSBG}}.$$

$$J=[1+b_{\mathrm{sh}}(M_{\mathrm{halo}})]{\int }_{s}d{s}^{\prime }{\int }_{\Omega }d{\Omega }^{\prime }{\rho }_{\mathrm{DM}}^2(s^{\prime },{\Omega }^{\prime })$$

$${\rho }_{\mathrm{DM}}(r)=\frac{{\rho }_s}{cr/r_{\mathrm{vir}}{\left[(cr/r_{\mathrm{vir}})+1\right]}^2},$$

$${\int }_{V}dV{\rho }_{\mathrm{DM}}^2(r)=\frac{M_{\mathrm{halo}}}{9}\Delta {\rho }_c(z)c^3{\left(\log (1+c)-\frac{c}{1+c}\right)}^{-2}\left(1-\frac{1}{(1+c{)}^3}\right),$$

$$\int ds\int d\Omega {\rho }_{\mathrm{DM}}^2(s,\Omega )\to \int dV{\rho }_{\mathrm{DM}}^2(r)/d_A^2.$$

$$J$$

$$\log {\mathcal{L}}_{\mathrm{stack}}(\mathcal{D}|\{⟨\sigma v⟩,J_1,J_2,...\})=\sum _{i,j}\log {\mathcal{L}}_{i,j,}^{\mathrm{ann}}({\mathcal{D}}_j|\{⟨\sigma v⟩,J_j\}),$$

$$f=\frac{n_{\mathrm{sig}}}{b_{\mathrm{eff}}},$$

$$b_{\mathrm{eff}}=\frac{N}{\left({\sum }_k\frac{P_{\mathrm{sig},k}^2(\mu )}{P_{\mathrm{bkg},k}(\theta )}\right)-1},$$

Full text

1 Introduction

It has been discovered that dark matter (DM) accounts for $\sim 25\%$ of the energy density of the Universe [1]. Despite great efforts over several decades DM has not been identified yet. Unravelling the nature of DM is one of the most important subjects in astrophysics. Weakly Interacting Massive Particles (WIMPs) have been considered one of the best theoretically motivated candidates for a dominant fraction of DM (e.g., [2, 3]). They can be produced in thermal equilibrium through interactions with standard model particles in the early Universe, and therefore some possible processes in which WIMPs self-annihilate into standard model particles can exist [4, 5, 6]. To explain the present abundance of DM, a total cross-section for such annihilation processes is estimated to be $\sim 3\times 10^{-26}~{}[\rm cm^{3}/s]$ for DM masses of order a few GeV (for a more accurate calculation, see [7]). In these processes, $\gamma$ rays are produced through primary or secondary processes (cascades into unstable states that lead to $\gamma$-ray production). Although $\gamma$ rays induced by those processes are expected to be rare, we can search for DM self-annihilation signals in regions of high DM concentration in the Universe.

Such $\gamma$-rays can be observed by detectors like the Fermi Large Area Telescope (Fermi-LAT) on board the Fermi Gamma-Ray Space Telescope. Fermi-LAT is the most sensitive $\gamma$ rays telescope in the energy range of 20 MeV to 1 TeV ([8, 9, 10, 11]). Thanks to its unprecedented angular resolution and sensitivity, Fermi-LAT has detected thousands of $\gamma$-ray sources [12, 13]. Close to the Galactic plane the Galactic diffuse emission produced by the interaction of energetic cosmic rays with interstellar gas, ambient photons and magnetic fields dominates. While at high latitudes, the isotropic $\gamma$-ray emission$-$thought to be of extragalactic origin$-$is the most dominant component and have been investigated about its origin (e.g., [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]).

In recent studies, nearby objects like the center of the Milky Way, dwarf spheroidal galaxies (dSphs), the Local group of galaxies and nearby galaxy clusters have been used as a targets for indirect searches. For example, several different teams [28, 29, 30, 31, 32, 33, 34, 35], have found an excess $\gamma$-ray emission in the Galactic Center which has been shown to roughly match the expected characteristics of a DM annihilation signal.

Although the origin of the Galactic Center excess has not yet been singled out, recent studies have presented strong evidence that the origin of this signal is related to the stars (e.g., [36, 37, 38]). It is true that the Galactic center is where the most luminous annihilating DM emission is expected, however, it is difficult to disentangle a DM signal from this region due to large background and foreground uncertainties.

The authors of [39, 40, 41, 42, 43, 44] have searched for DM emission in Milky Way dSphs, which are perhaps the most desirable targets to date because they are relatively close to us and have less astronomical $\gamma$-ray contamination. In fact, in the case of the Milky Way dSphs, the strongest constraints on the velocity averaged DM cross-section have been provided.

In this work, we propose a new target for probing DM annihilation; low surface brightness galaxies (LSBG), which are known to have surface brightness of order  $>23~{}{\rm mag/arcmin^{2}}$ in most cases. Due to their very low brightness, LSBGs are more difficult to be identified than ordinary galaxies. Currently, surveys with high sensitivity like Subaru Hyper Suprime-Cam (HSC), has made it possible to discover nearby LSBGs in large amounts. The mass-to-light ratio of LSBG is likely to be large (50-100), thus, those systems are highly dominated by DM [45, 46]. LSBGs have more massive halos than Milky Way dSphs and present less star-formation, pulsar or supernovae emission than those of ordinary galaxies or galaxy clusters [47]. Therefore, it is expected that they are mostly quiescent in $\gamma$-rays.

Since LSBGs have a smaller angular size than the point spread function (PSF) of the Fermi-LAT, it is not possible to decompose LSBGs into their substructures. So in practice, LSBGs can be treated as point sources in our $\gamma$-ray analysis. Importantly, since our sample of LSBG are located at high latitudes, the impact of uncertainties in the Galactic diffuse emission is expected to be small. This will naturally allow us to obtain DM limits that can be regarded as being robust to systematic uncertainties in the Galactic diffuse emission model.

In this analysis, we use an LSBG catalog [48] including $\sim 800$ objects by Subaru HSC observations of $\sim 200~{}{\rm deg^{2}}$ of the sky. Unfortunately, we currently have only eight LSBGs with known accurate redshifts, hence here we only employ this small sample of LSBGs. Using $\gamma$ ray observations taken with the Fermi-LAT, we constrain the DM velocity averaged cross-section for DM masses of order a few GeV. Next generation telescopes like the Large Synoptic Survey Telescope (LSST) is expected to discover large amounts of LSBGs due to the large survey areas it can cover and sufficient depths, which will improve our DM constraints significantly. In the future, the importance of searching for DM annihilation using LSBGs should increase and it is worthwhile to probe the DM annihilation using LSBGs at present.

This paper is organized as follows. In Section 2, we introduce a sample of HSC LSBGs and $\gamma$-ray observation by Fermi-LAT used in our analysis, In Section 3, we derive formulations of the $\gamma$-ray flux induced by the DM annihilation within LSBGs and describe the method of analysis to probe the annihilation cross-section of DM. Our results are shown in Section 4 and finally, summary and discussion are described in Section 5.

2 Data

2.1 HSC LSBG sample

HSC is a wide field imaging camera installed on the prime focus of the Subaru telescope [49, 50].With this imaging camera and superb site conditions, seeing $\sim 0^{\prime\prime}.6$, HSC survey enables us to measure the shape and photometry for large number of galaxies out to redshift $z\sim 1.5$. The HSC survey has three survey layers of Wide, Deep and Ultra deep, and in the Wide layer, we observe in five visible broad bands, $g,r,i,z$ and $y$. The observation data set of the Wide layer for the first internal data release S16A [51] has survey areas of $\sim 200~{}{\rm deg^{2}}$ of the sky, which is eventually going to reach to $1400~{}{\rm deg^{2}}$ in the final data release. Limiting magnitudes for $5\sigma$ point-source detection of $g,r,i,z,$ and $y$-band are 26.5, 26.1, 25.9, 25.1, and 24.4, respectively. The observation data is processed by an open-source software called hscpipe (see [52] for the detail), which has been developed as a pipeline for the LSST data [53, 54].

In [48], the LSBG catalog has been produced using the observation data set of the Wide layer with S16A data. In order to avoid the limitation of survey areas by requiring all bands, they have used only three bands ($g,r,i$-band), because the survey progress is different by photometric bands. Using hscpipe, they produce sky-subtracted co-add images and divide these images into equivalent rectangular regions called tracts (each tract has an area of $1.7~{}{\rm deg^{2}}$). Further, tracts are divided into $9\times 9$ grids of patches and each patch is pixelized into $4200\times 4200$ pixels. To identify LSBG objects, the validity of the sky-subtraction is important because of mean surface brightness of target objects fainter than the brightness of the night sky. As described in [52], an algorithm for the sky-subtraction used in hscpipe causes over-subtracting the background around extended sources ($>1^{\prime}$). Therefore, LSBGs found in [48] can be biased against detection of LSBGs around such extended sources.

hscpipe is optimized to identify faint and small objects like distant galaxies. If we use the pipeline for extended objects, the single extended object are decomposed into multiple child objects, so called “shredding”. Typically, LSBG in the HSC LSBG catalog is divided into more than 10 child objects. Moreover, LSBG brightness is comparable to the sky background noise and thus the shape and surface brightness of LSBGs are likely misestimated. Therefore in procedures to identify LSBGs, they have developed a pipeline based on the LSST codebase mainly, instead of use of hscpipe. Also they have used SExtractor [55] to compose an initial catalog of LSBG candidates and then select them based on size and color measurements. Finally imfit [56] is used to refine their parameters. Above procedures are carried out on a patch-by-patch basis. LSBGs are finally defined such that the mean surface brightness is larger than 24.3 mag arcmin*-2*. Finally, 781 LSBGs are detected within a HSC wide S16A footprint. The LSBGs identification procedure is summarized as follows,

First of all, bright sources and associated diffuse lights are subtracted from images because they might mimic the LSBGs. Those subtracted region are replaced to the randomized background noise. 2. 2.

After smoothing with the Gaussian kernel (FWHM = $1^{\prime\prime}$), sources with $2.5^{\prime\prime}<r_{1/2}<20^{\prime\prime}$, where $r_{1/2}$ is the half-light radius, are extracted. Furthermore, by applying reasonable color cuts, optical artifacts and high-redshift galaxies are removed. 3. 3.

By modeling the surface brightness profiles of LSBG candidates, astronomical false positive are removed. Finally, remaining false candidates, which is typically point-like sources with diffuse background lights, are removed by visual inspection. Final LSBG samples are 781 objects.

Basically, these procedures are applied to the $i$-bands images; however, in order to reduce the false detection including any artificial effects, all LSBG candidates are required to be detected in $g$-band images as well. In this work, we use only eight LSBGs out of the 781 HSC LSBGs, where accurate distances to the LSBGs are measured. [48, 57]. In Table 1, parameters of these eight LSBGs are summarized and in Figure 1, each LSBG image used in this paper is shown.

2.2 LAT data reduction

In this work, we use 8 years (2008-08-04 to 2016-08-02) of Fermi-LAT observations. We have selected the photon event class P8R3 SOURCE, which is the recommended class of events for point source analysis. Only events within the energy range 500 MeV to 500 GeV are selected, and these are further binned in 26 logarithmically spaced energy bins in our analysis. The LAT possesses a low angular resolution for energies $\lesssim 500$ MeV which could introduce a bias in the analysis due to possible point source confusion, while at high energies the photon count statistics decreases. The chosen low energy limit is a compromise between sensitivity and statistics. For each LSBG shown in Tab. 1 we select a patch of the sky of size $10^{\circ}\times 10^{\circ}$ with centroid in each target object. In the construction of the counts maps we have employed spatial bins of size $0.1^{\circ}$ and the corresponding instrument response functions (IRFs) for our photon class event  [9], P8R3_SOURCE_V2 222https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/IRF_overview.html. As recommended in the cicerone333https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone, we applied the quality cuts DATA_QUAL>0 && LAT_CONFIG==1. Furthermore, we reject events with zenith angles larger than $100^{\circ}$ to avoid contamination of photons produced by interaction of cosmic rays with the Earth’s atmosphere. The analysis of the LAT data was done with fermipy444https://fermipy.readthedocs.io/en/latest/ [58], which is an open-source software package based on the Fermi Science Tools555https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ (v11r5p3), including some additional high-level tools for common analysis tasks.

In our analysis, we first perform maximum likelihood runs to estimate the parameters for all $\gamma$-ray sources in our ROIs. The Galactic foreground and background $\gamma$-ray emission is modeled with the standard diffuse Galactic foreground model, gll_iem_v07.fits, and the isotropic emission model, iso_P8R3_SOURCE_V2_v01.txt, respectively. As for the resolved $\gamma$-ray sources, we derive them from the 4FGL catalog [13]. The amplitude and spectral shape of all the $\gamma$-ray sources lying in our ROIs were varied in our maximum likelihood runs while the sources which reside outside our ROIs and within region of $15^{\circ}\times 15^{\circ}$ centered on LSBG were fixed to the their catalog values. Figure 2 displays the photon count maps around each LSBG as observed by the Fermi-LAT. The corresponding residual maps obtained after subtracting the best-fit background/foreground emission model are also shown in the same figure. These images are constructed using the full energy range (500 MeV to 500 GeV). We note that photons counts above a few of tens of GeVs hardly contribute to the overall observed emission.

Following the procedure of [40], we search for excess emission at the position of each LSBG. This is done separately in each energy bin in order to derive flux constraints that are not sensitive to their assumed spectrum. To obtain convergence in our maximum likelihood fits, we utilize the following parameter relaxation method: First, we optimize the amplitudes of all $\gamma$-ray sources in order of higher intensity sources down to fainter ones. Next, we repeat our maximum likelihood run starting from the updated source parameters and freeing all the spectral shape parameters. The last step of our procedure consists in running the fermipy pipeline for the flux upper limits calculations.

In the computation of the flux upper limits we follow the prescription given in the 2FGL catalog [59]. Namely, we employ a Bayesian [60] method recommended for analyses of very faint $\gamma$-ray point sources. For a rigorous application of this method, we first checked that our LSBGs had test statistics (TS) values of less than one in most energy bins. In addition, since the point-source fluxes are assumed to be positive, the likelihood function can be approximated to a $\chi^{2}/2$ distribution. In this sense, the 95% confidence level (C.L.) flux upper limits are obtained by varying the spectral normalization parameter of the putative LSBGs until $\log\mathcal{L}_{0}-\log\mathcal{L}\sim 2$ (where $\mathcal{L}_{0}$ is the likelihood for zero-flux). (see Section 3.5 of [59]).

Figure 3 shows the likelihood profiles for all our target LSBGs. We note that some of our ROIs overlap with each other, however, in the computation of the combined limits with a joint likelihood procedure we assume that all of our ROIs are independent. Details on the methods utilized for converting flux upper limits into constraints for the DM annihilation cross-section are given in the next section.

3 Methods

3.1 Modelling of the Dark Matter annihilation signal

In this section, we provide details on our modelling methods for the expected $\gamma$-ray emission from DM annihilation as well as on the pipeline used to constrain the DM model parameters.

The predicted $\gamma$-ray annihilation flux from each LSBG can be written as

[TABLE]

where $z$ is the redshift of the LSBGs, $\langle\sigma v\rangle$ is the ensemble average of the annihilation cross-section and relative velocity of DM particles, $E_{\gamma},E_{\gamma}^{\prime}$ are observed and emitted $\gamma$-ray energy, $m_{\chi}$ is DM mass, $dN/dE_{\gamma}$ is the $\gamma$-ray energy spectrum and $Br_{i}$ is the branching ratio in the $i$-th annihilation channel, respectively. The J-factor, $J$ is fully specified by astronomical parameters of target halos, such as the halo mass or distance to the halo. We note that any astronomical factors to be considered for the $\gamma$-ray flux are absorbed in the J-factor. Other factors in the right-hand side of Eq. (3.1) are fully given by the particle physics nature of the DM candidate. In general, DM particles can self-annihilate in various final state particles such as the $b\bar{b},\tau^{+}\tau^{-},\mu^{+}\mu^{-}$ and $W^{+}W^{-}$ as well as $\gamma\gamma$ and $\gamma Z$, which in turn can produce secondary $\gamma$-ray emission through their subsequent cascade decay. In our analysis, we consider the single annihilation channel $b\bar{b}$, which is a representative channel and benchmark to assess the sensitivity of the LAT to a putative DM signal from our search targets. The DM spectra is obtained with the DMFIT [61] package provided within fermipy. The DMFIT666https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/gammamc_dif.dat itself provides an interpolation to DM $\gamma$-ray spectra tables extracted from DarkSUSY [62].

Amongst the various effects expected to impact the $J$-factors, the mass estimate should include the largest uncertainty. Here we describe the way to estimate the halo mass for the eight LSBGs with HSC optical data. In particular, we convert the observed $g,r$ and $i$ band magnitudes into $V$ band magnitude using [63],

[TABLE]

Given the luminosity distance $d$, we can convert apparent magnitude to absolute magnitude $M_{V}$, by using the well known relation $M_{V}=V+5-5\log_{10}d$. Also, by assuming the mass-to-light ratio to be unity [64] we obtain

[TABLE]

Finally, we apply the stellar to halo mass ratio [65] for the crude estimate of the total halo mass $M_{\rm halo}$ of our LSBG sample. In Table 2, stellar mass and halo mass of the eight LSBGs are summarized. Below we describe how the $M_{\rm halo}$ estimates are implemented in the J-factor calculations.

The J-factor is defined as the line of sight integral of the squared DM density profile,

[TABLE]

where $\Omega$, $\rho_{\rm DM}$ and $M_{\rm halo}$ are the solid angle of the target object, the DM density profile and the halo mass, respectively. The the boost factor $b_{\rm sh}(M_{\rm halo})$ takes into account the excess of the annihilation rate due to the substructure present in the halo. In our analysis, we employ the model proposed in Ref. [66] for the calculation of the boost factor and assume $b_{\rm sh}(M_{\rm halo})=1$ for all 8 LSBGs. As for the dark matter profile, we assume a smooth NFW density profile [67],

[TABLE]

where $r_{\rm vir}$ is the virial radius, $\rho_{s}$ and $c$ are the normalization and concentration. We adopt the mass-concentration relation obtained in [68]. The parameters $\rho_{s}$ and $c$ can be uniquely determined once the DM halo mass is specified. It can be shown that the volume integral of the $\rho^{2}_{\rm DM}$ can be reduced to

[TABLE]

where $\Delta$ is the overdensity of spherical collapse to be assumed as $200$ and and $\rho_{c}$ is the critical density of the Universe. Since the angular size the LSBGs considered in this study is much smaller than the PSF of the LAT instrument in all energy bands, we can in practice treat all the LSBGs as $\gamma$-ray point sources in our analysis pipeline. It follows that the line of sight integral can be replaced as,

[TABLE]

Thus, eq.(3.4) can be simplified as

[TABLE]

where $d_{A}$ is the comoving angular diameter distance to the LSBG. To evaluate the uncertainty on the $J$-factor, we perform Monte-Carlo error estimation by assuming 1-$\sigma$ Gaussian error, $\Delta\log_{10}c=0.1$ [65]. Halo mass uncertainties can be converted from $i$-band magnitude uncertainty from Table 1 which can read as 1-$\sigma$ Gaussian errors with $\Delta\log_{10}M_{\rm halo}\sim 0.4$. We ignore the error on the distance as it is sufficiently accurately measured. Furthermore, astrophysical uncertainties in the matter profile of LSBGs result in J-factor uncertainties of at most 0.2 dex [69].

3.2 Composite Analysis

As it can be inferred from Eq. 3.1, the predicted DM energy annihilation spectra is independent of the search target, only the J-factors present differences in our sample of LSBGs. This characteristic allows us to combine data from all individual LSBGs in order to set stronger constraints on the DM model parameters. After first analyzing every individual LSBG, a composite analysis is subsequently performed following the methods employed by the Fermi team in dSphs [40]. Specifically, our pipeline adds together the photon counts from the signal and background regions for each LSBG and then computes upper limits on the signal region following the same prescription as for individual sources. We note that this composite likelihood technique takes into account that every LSBG has a different J-factor. All the LSBGs are found to have very low statistical significance detection in the Fermi data. Therefore, in the computation of the 95% C.L. upper limits of the DM annihilation cross-section, we employ the Bayesian [60] procedure recommended in the 2FGL catalog [59] for analyses of dim point sources.

Using our bin-by-bin method we perform likelihood scans as a function of putative DM fluxes associated to each LSBG in every individual energy bin. Since all the likelihood values obtained at each ROI are assumed to be independent of each other, the joint likelihood ${\mathcal{L}}_{\rm stack}$ for the full sample of targets can be written as

[TABLE]

where ${\cal L}^{\rm ann}_{i,j}$ is the log-likelihood value obtained for the $j$-th LSBG, $i$-th energy bin and $J_{j}$ is the $J$-factor of one particular LSBG. In this formulation, we assume that our all LSBGs are statistically independent from each other. Since our sample of LSBGs are treated as point sources, this assumption is correct if the angular separations between LSBGs are larger than the LAT PSF at the energy range in question. We will further discuss in Sec. 4.1.

As described in Section 2.2, in our procedure we assume that the flux of the LSBGs is positive definite. This implies that in the limit of large number counts, the data is well described by a $\chi^{2}/2$ distribution. In this sense the 95% C.L. upper limits on the velocity averaged DM cross section can be obtained when the total likelihood $\Delta\log{\mathcal{L}}_{\rm stack}({\cal D}|\langle\sigma v\rangle)\sim 3.8/2$ (for 1 degree of freedom). We estimate the astrophysical uncertainties associated to the upper limits computed this way by performing Monte Carlo simulations of the J-factor distributions for the eight LSBGs considered in this work.

4 Result

4.1 Dark Matter Constraints

As shown in Figure 3, the sample of LSBGs considered in this work were not significantly detected in $\gamma$-rays. Therefore, we use them to impose constraints on the DM annihilation cross section. The computation of the upper limits is done using the methods explained in the previous section. Figure 4 displays the 95% C.L. upper limits on $\langle\sigma v\rangle$ as obtained for each individual LSBG (dashed lines) using their respective median J-factor. The joint upper limit for the full sample (black solid line) was obtained with the joint likelihood method described in Eq. 3.9. The green band displays the impact on our limits due to astrophysical uncertainties in the DM model parameters. As discussed earlier, these come mainly from uncertainties in the halo mass and the matter concentration of the LSBGs. As it can be seen, the objects LSBG-285 and LSBG-575 place the strongest constraints on $\langle\sigma v\rangle$ and also provide dominant contributions to the joint constraint. Even after stacking over the full sample of eight LSBGs, our joint constraints are weaker than the ones obtained using more traditional targets like dSphs or nearby galaxy groups and clusters of galaxies [40, 69]. We note that in our stacking procedure we assume that the log-likelihood values obtained for each LSBG are statistically independent. However we found that there are two pairs of LSBGs – (LSBG-456, LSBG-464) and (LSBG-613, LSBG-729) – that are separated by only $\sim 2^{\circ}$. This separation is comparable to the size of LAT PSF in our low energy bins and it is possible that these objects might bias our results. In order to estimate the impact that this assumption had on our results we computed the corresponding cross-covariance for both pairs of LSBGs. In particular, we found that the impact on our upper limits was at the $\sim 10\%$ level. This demonstrates that the statistical independence assumption is appropriate for our analysis.

In addition, we have estimated the systematic uncertainties introduced by our Galactic diffuse emission model. We followed the rigorous approach recommended in Ref. [23]. In that reference three different Galactic diffuse emission model were constructed using the Cosmic Rays (CR) propagation code GALPROP777http://galprop.stanford.edu. The three different model named Model A, B and C encapsulate a wide range of uncertainties in the interstellar gas column density distribution, CRs source distribution and energetics as well as the diffusion coefficient and Galactic magnetic fields. More details about those models are given in [23]. Using the alternative Galactic diffuse emission models we repeated our $\langle\sigma v\rangle$ upper limits calculation and found that those are affected at the few percent level for DM mass values smaller than 100 GeV, while no difference was apparent for larger DM masses. This is shown by the black dotted lines in Figure 4.

4.2 Forecast for future surveys

Even though we obtain a relatively weak constraint on the DM cross section with our sample of 8 LSBGs, we expect fairly stringent constraints with a much larger sample. As already mentioned in Sec.2, 781 LSBGs are found in the S16A HSC-wide footprint which roughly covers $200~{}{\rm deg^{2}}$. We estimate that ${\cal O}(10^{4})$ LSBGs will eventually be found within the $1400~{}{\rm deg^{2}}$ area of the HSC wide survey. In addition, future optical imaging surveys like LSST are expected to have better sensitivities which may make it possible to identify more nearby LSBGs. For example, the LSST will cover $\sim 2\times 10^{4}~{}{\rm deg^{2}}$ of the sky with the limiting magnitude $\sim 27$ mag in $g,r$ and $i$-band. If the LSBGs are discovered by the LSST survey with the same rate of the HSC, we might be able to disentangle of order ${\cal O}(10^{5})$ LSBGs in the near future. However, incorporating all those objects in an analysis of this kind is unrealistic because it is difficult to obtain accurate redshift measurements for all of them. Conservatively, we assume that only a fraction of the soon-to-be-discovered LSBGs will have reliable distance measurements. In particular, we assume that future LSBG catalogs will have the same redshift measurement rate as that of the current HSC-LSBG catalog (8 objects within $\sim 200\;{\rm deg}^{2}$).

We further assume that future observations of LSBGs will have a similar redshift and mass distribution to those in our current HSC catalog. Under this considerations, we estimate that the number of LSBGs to be detected is 60 and 400 for the ultimate HSC catalog and LSST, respectively. In the case of the LSST forecast, we consider only half of the entire LSST footprint to avoid the Galactic plane region, as this is a region that is highly contaminated by Galactic diffuse $\gamma$-ray emission.

Figure 5 shows our predicted constraints on the DM cross section using the ultimate HSC catalog and forthcoming LSST objects.

In our forecast analysis we perform a joint analysis assuming that objects in each survey area have the same likelihood as our eight LSBGs. In Figure 5 we show the expected constraints for the ultimate HSC and LSST surveys along with the limits obtained with a simple scaling of the number of LSBGs. This scaling is done using the results of each survey area from the analysis of the current catalog of 8 LSBGs with measured redshifts.

To understand the benefits of additional LSBG discoveries for DM searches it is useful to consider the signal fraction $f$ defined [70] as

[TABLE]

which approximately denotes the size of a DM signal ($n_{\rm sig}$ in units of photon counts) relative to the effective background photons ($b_{\rm eff}$). The quantity $f$ defined this way can be seen as a generalization of the signal-to-noise-ratio. In addition, the effective background can be estimated from the maximum likelihood fit covariance matrix (e.g., Eq.F.2 in Ref. [70]) and is given by

[TABLE]

where the index $k$ runs over all the energy bins and pixels in each ROI and $N$ is the total number of observed photons in the corresponding ROI. The normalized signal and background model components can be straightforwardly computed as $P_{{\rm sig,}i}({\bf\mu})=\lambda_{{\rm sig,}i}({\bf\mu})/\sum_{k}\lambda_{{\rm sig,}k}({\bf\mu})$ and $P_{{\rm bkg,}i}({\bf\theta})=\lambda_{{\rm bkg,}i}({\bf\theta})/\sum_{k}\lambda_{{\rm bkg,}k}({\bf\theta})$, where $\lambda_{\rm sig}$ and $\lambda_{\rm bkg}$ represent the expected photon flux from a DM signal and the background flux, respectively. The signal parameters $\mu$ are relevant DM parameters like $\langle\sigma v\rangle$ or $J$-factor. Equivalently, the nuisance parameter $\theta$ can be either the spectral normalization or index of all the $\gamma$-ray background/foreground sources in the ROIs (see more details in the Appendix). One noteworthy property of $b_{\rm eff}$ is that if the expected signal model and background model differ significantly such that $\sum_{k}P_{{\rm sig,}k}^{2}({\bf\mu})/P_{{\rm bkg,}k}({\bf\theta})\gg 1$ then $b_{\rm eff}$ will be proportionally less than $N$. In our case, this condition can be satisfied away from the Galactic plane due to the rapidly falling energy spectra of the background and smaller overall intensity. We also note that in the high energy regime the number of observed photons $N$ is drastically reduced thus making the search for DM in LSBGs signal limited targets.

As explained in a previous section, for a certain LSBG, the number of signal events $n_{\rm sig}$ is proportional to the product of the J-factor and $\langle\sigma v\rangle$. By making the simple assumption that future LSBGs will have the same J-factors as those currently known, is easy to get that the total $n_{\rm sig}$ considered in the joint likelihood analysis is proportional to the number of LSBGs $N_{\rm LSBG}$ times $\langle\sigma v\rangle$. Upper limits at the 95% CL on the total number of signal events $n_{\rm sig}$ are then obtained by increasing $n_{\rm sig}$ (that is, increasing $\langle\sigma v\rangle$), while varying all other parameters, until $-\Delta\log{\cal L}_{\rm stack}=2.0$ from the best-fit point.

We emphasise that it is essential to increase the number of spectroscopically confirmed LSBGs with future observations in order to obtain better and more efficient constraints on the DM annihilation cross sections using our proposed method.

We also note that the forecast shown in Figure 5 does not consider the actual distribution of the $J$-factor, which may affect the power of the predicted constraints. As can be seen in Eq. (3.4), the $J$-factor depends on the inverse square of the distance to the object but is only proportional to the DM halo mass. This means that in case that future surveys find LSBGs closer to us than our current catalog of 8 LSBGs, then we might be able to obtain DM constraints that are stronger than those shown in the conservative predictions of Figure 5.

5 Summary and Discussion

For the first time, we have searched for dark matter emission in LSBGs using 8-years of observations of $\gamma$-ray data collected by Fermi-LAT. No statistically significant signal was detected in these objects and thus we imposed new constraints on the DM annihilation cross-section. For this we used 8 LSBGs recently discovered by the HSC survey with reliable distance measurements. The advantages of using LSBGs for constraining DM properties over other possible targets are: First, they are expected to have a low level of unresolved point point source contamination within. Secondly, they are typically 10 times more massive than e.g. Milky Way dSphs. Lastly, the expected number of LSBGs will be very large (in fact much larger than that of dSphs). The main disadvantage of this novel dark matter target, compared to dSphs, is that they are much farther away and so they typically have smaller J-factors. However, we note that LSBGs are still located in the local Universe and the sensitivities can be incremented by increasing the detection of nearby LSBGs. Since the angular extensions of our LSBGs is sufficiently smaller than the Fermi PSF, in practice we were able to treat them as $\gamma$-ray point sources. This made our analysis less sensitive to the model systematics and potential biases.

We derived 95 % C.L. upper limit on $\langle\sigma v\rangle$ for each LSBG and also with a joint likelihood analysis using the full set of 8 objects. We found that $\sim 10^{-23}~{}[\rm cm^{3}/s]$ for a DM mass of 10 GeV self-annihilating into the $b\bar{b}$ channel. This constraint is $\sim 10^{3}$ times weaker than that obtained from Milky Way dSphs at the same DM mass (which is currently the most stringent). We showed that a limiting factor in our methods corresponds to the low number of LSBGs considered. In the footprint of our S16A HSC data set, almost 800 LSBGs are identified. However, their distances are well measured for only the 8 objects used in this paper. Given the broad distance distribution of the entire sample, we presented a detailed systematic uncertainties analysis where we studied the impact of uncertainties in the distance measurement in our dark matter limits.

In the high energy limit where the number of detected photon counts is very low, the constraints on the annihilation cross-section scales with the number of objects $N$, instead of $\sqrt{N}$. Indeed, above $\sim 20$ GeV, Fermi-LAT has rarely observed $\gamma$-ray photons in our ROIs. It follows that the number of objects is essential in obtaining constraints on the DM annihilation cross section. Although the current LSBG catalog used in this work contains $\sim 800$ objects in the HSC $200~{}{\rm deg^{2}}$ footprint, this is going to increase to $\sim 6000$ in the final data release of HSC when the HSC survey finish to cover $1400~{}{\rm deg^{2}}$. Moreover, in the near future, LSST will cover over $20000~{}{\rm deg^{2}}$ of sky with the depth of down to $27$ mag, and more than $10^{5}$ LSBGs can be expected to be identified. We found that by means of conservative forecast analysis the expected constraints on the DM cross section using our method will not be stronger than those using dSphs.

However, if the number of spectroscopically followed up LSBGs significantly increases, the constraint using those may reach to ones using dSphs because they are scaled with the number of detected LSBGs in the high DM mass. This is particularly the case in the high energy regime where the effective background is small compared to the expected signal counts.

Acknowledgments

We are grateful to Jenny E. Greene, Johnny P. Greco, Michael A. Strauss and Yoshiyuki Inoue for their careful read of the manuscript and useful comments. This work is in part supported by MEXT Grant-in-Aid for Scientific Research on Innovative Areas (No. 15H05887, 15H05890, 15H05892, 15H05893, JP18H04340, 18H04359, 18J00277). 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.

Appendix

In many cases, DM annihilation signals from objects of interest are expected to be small compared to the astrophysical background/foreground. Therefore, to discuss how the constraints on the annihilation cross-section are expected to increase with the number of objects included in the joint likelihood, we consider the useful definition of $b_{\rm eff}$ introduced in Eq. 4.2. In practice, the quantities $\lambda_{\rm sig}$ and $\lambda_{\rm bkg}$ (necessary to compute $b_{\rm eff}$) can be computed with the use of the gtmodel tool within FermiTools. For $\lambda_{\rm sig}$, we computed infinite-statistics model maps, where the model only contained the DM target in question. As for $\lambda_{\rm bkg}$, we ran gtmodel assuming our best-fit background/foreground model obtained with our maximum-likelihood fitting procedure. Simple arithmetic operations follow to obtain $b_{\rm eff}$.

Reference [71] defines three different regimes in which the behaviour of the DM limits improve differently with additional data taking. First, in the systematics-limited regime ($b_{\rm eff}\gtrsim 2000$), the DM limits are expected to augment incrementally as our knowledge of the astrophysical background improves. DM limits obtained from targets in this regime are not expected to improve dramatically with additional data taking. Second, for objects in the background limited regime ($50\lesssim b_{\rm eff}\lesssim 2000$), the constraints can improve as $\propto\sqrt{N_{\rm obj}}$ (where $N_{\rm obj}$ is the number of objects). Third, for objects in the signal-limited regime ($b_{\rm eff}\lesssim 50$) the DM limits improve with $\propto N_{\rm obj}$. In Figure 6, we show the effective background as a function of DM mass for one of the LSBGs in our study (LSBG-285). As can be seen, $b_{\rm eff}$ is in the signal-limited regime for $M_{\rm DM}\gtrsim 100$ GeV and is out of the systematic-limited regime for all DM masses. Note that for lower masses ($M_{\rm DM}\lesssim 20$ GeV), the DM constraint is expected to scale only with $\sqrt{N_{\rm obj}}$. However, as seen in Figure 5, the constraint becomes more stringent. This seems to oppose the statement above but this is simply due to the lack of the statistics and the level of background is highly dependent on the sky position. The constraint will scale with $\sqrt{N_{\rm obj}}$ once we obtain a good amount of samples. On the contrary, in the signal-limited regime, as the background is subdominant, the fluctuation of the background on the sky position does not affect the constraints significantly. This is why we obtain a clean scaling relation of the DM constraint at high mass scale where the thermal relic scenario is still not excluded.