Foreground effect on the $J$-factor estimation of ultra-faint dwarf spheroidal galaxies

TL;DR

This paper evaluates how foreground star contamination affects the estimation of the J-factor in ultra-faint dwarf spheroidal galaxies, proposing a data-driven method to improve the robustness of dark matter searches.

Contribution

It introduces a likelihood analysis method that accounts for foreground contamination and demonstrates how to estimate this effect using mock data and future spectroscopic surveys.

Findings

Foreground contamination can be effectively estimated using control regions.

The proposed method improves the stability of J-factor estimation.

Future spectrographs can facilitate better foreground star discrimination.

Abstract

Dwarf spheroidal galaxies (dSphs) are promising targets for the gamma-ray dark matter (DM) search. In particular, DM annihilation signal is expected to be strong in some of the recently discovered nearby ultra-faint dSphs, which potentially give stringent constraints on the $\mathcal{O}(1)$ TeV WIMP DM. However, various non-negligible systematic uncertainties complicate the estimation of the astrophysical factors relevant for the DM search in these objects. Among them, the effects of foreground stars particularly attract attention because the contamination is unavoidable even for the future kinematical survey. In this article, we assess the effects of the foreground contamination on the astrophysical $J$-factor estimation by generating mock samples of stars in the four ultra-faint dSphs and using a model of future spectrographs. We investigate various data cuts to optimize the quality of the data and apply a likelihood analysis which takes member and foreground stellar distributions into account. We show that the foreground star contaminations in the signal region (the region of interest) and their statistical uncertainty can be estimated by interpolating the foreground star distribution in the control region where the foreground stars dominate the member stars. Such regions can be secured at future spectroscopic observations utilizing a multiple object spectrograph with a large field of view; e.g. the Prime Focus Spectrograph mounted on Subaru Telescope. The above estimation has several advantages: The data-driven estimation of the contamination makes the analysis of the astrophysical factor stable against the complicated foreground distribution. Besides, foreground contamination effect is considered in the likelihood analysis.

Tables

Table 1. Table 1: The input parameters of each dSph. The distances from the earth and projected half-light radii are shown by d 𝑑 d and r e subscript 𝑟 𝑒 r_{e}   (McConnachie, 2012 ) . The DM halo and kinematical parameters ρ s subscript 𝜌 𝑠 \rho_{s} , r s subscript 𝑟 𝑠 r_{s} , α 𝛼 \alpha , β 𝛽 \beta , γ 𝛾 \gamma , and β ani subscript 𝛽 ani \beta_{\text{ani}} are determined by fitting the stellar data provided by 200 for Ursa Major II, Coma Berenices, Ursa Major I and 201 ( a ) for Segue 1 under the same procedure as Geringer-Sameth et al. ( 2015 ) . Quoted values are those obtained from the χ 2 superscript 𝜒 2 \chi^{2} minimization. #3 #3 #3 The slopes of the dSph DM profile are known to be poorly constrained, so that the corresponding values in the table should be regarded as one of possible choices to generate merely mock data of UFDs. The J Input subscript 𝐽 Input J_{\text{Input}} shows the J 𝐽 J -factors calculated within an angular radius of 0.5 ​ degree 0.5 degree 0.5\,\text{degree} under the input DM halo parameters and distance.

Model dSph	$d$ [kpc]	$r_e$ [pc]	${\log }_{10}\left(\frac{{\rho }_s}{\left[M_{⨀}/{\text{pc}}^3\right]}\right)$	${\log }_{10}\left(\frac{r_s}{\left[\text{pc}\right]}\right)$	$\alpha$	$\beta$	$\gamma$	$-{\log }_{10}\left(1-{\beta }_{\text{ani}}\right)$	${\log }_{10}\left(\frac{J_{\text{Input}}}{\left[{\text{GeV}}^2/{\text{cm}}^5\right]}\right)$
Ursa Major II	32	149	-0.370	2.62	2.36	3.28	0.0328	-0.975	19.70
Coma Berenices	44	77	-0.283	2.27	2.87	6.79	0.178	-0.894	18.74
Segue 1	23	29	0.306	1.93	0.973	3.94	1.15	-0.00155	20.05
Ursa Major I	97	319	0.587	1.97	2.89	8.04	0.302	-0.625	18.92

Table 2. Table 2: The bulk velocity, the truncation radius and cut conditions for each dSph. The bulk velocity of Ursa Major II, Coma Berenices, and Ursa Major I is from Simon & Geha ( 2007 ) and Segue 1 from Simon et al. ( 2011 ) . The truncation radii are from Geringer-Sameth et al. ( 2015 ) . See the text for more details.

Model dSph	$v_{\text{dSph}}$ [km/s]	$r_{\text{max}}$ [pc]	${\log }_{10}{(g/[\text{cm}/{\text{s}}^2])}_{\text{lower}}$	${\log }_{10}{(g/[\text{cm}/{\text{s}}^2])}_{\text{upper}}$	${\left[\text{Fe/H}\right]}_{\text{lower}}$	${\left[\text{Fe/H}\right]}_{\text{upper}}$
Ursa Major II	-116.5	294	$0.2$	$4.9$	$-4.5$	$-1.5$
Coma Berenices	98.1	238	$0.1$	$4.7$	$-4.3$	$-1.5$
Segue 1	208.5	139	$0.9$	$5.1$	$-6.1$	$-1.2$
Ursa Major I	-55.3	732	$0.0$	$3.7$	$-4.9$	$-0.9$

Table 3. Table 3: The averaged numbers of the member (foreground) stars are given by N Mem subscript 𝑁 Mem N_{\text{Mem}} ( N FG subscript 𝑁 FG N_{\text{FG}} ). The Raw column shows the numbers of the stars after the colour-magnitude cut and the cut of the region of interest. The details of the naive cuts and membership selection are given in the text.

	Condition		Raw		Naive cut		Membership selection
Model dSph	${\theta }_{\text{ROI}}$ [degree]	$i_{\text{max}}$ [mag]	$N_{\text{Mem}}$	$N_{\text{FG}}$	$N_{\text{Mem}}$	$N_{\text{FG}}$	$N_{\text{Mem}}$	$N_{\text{FG}}$
Ursa Major II	0.65	21	80	829	76	75	54	5
		21.5	150	988	141	103	89	4
		22	233	1149	214	132	131	4
Coma Berenices	0.65	21	35	579	34	58	29	2
		21.5	58	743	55	85	44	2
		22	92	898	85	110	66	1
Segue 1	0.65	21	24	620	22	60	19	1
		21.5	43	748	39	84	34	0
		22	61	953	56	123	49	0
Ursa Major I	0.65	21	42	680	37	32	26	1
		21.5	55	831	48	39	34	1
		22	63	953	56	44	38	1

Table 4. Table 4: The numerical value of the J 𝐽 J -factor obtained by our fit, illustrated in Fig.  3 . We show the averaged median values of log 10 ⁡ ( J / [ GeV 2 / cm 5 ] ) subscript 10 𝐽 delimited-[] superscript GeV 2 superscript cm 5 \log_{10}(J/[\text{GeV}^{2}/\text{cm}^{5}]) (the first values), the averages of the 68% quantiles (the first uncertainties) and the standard deviations of the median values (the second uncertainties).

		Mock Ursa Major II	Mock Coma Berenices	Mock Segue 1	Mock Ursa Major I
$i=21$	KI17	${19.49}_{-0.27}^{+0.24}\pm 0.15$	${18.61}_{-0.66}^{+0.60}\pm 0.44$	${20.07}_{-0.65}^{+0.75}\pm 0.38$	${18.38}_{-0.75}^{+1.51}\pm 0.58$
	Conventional	${19.20}_{-0.24}^{+0.30}\pm 0.20$	${18.65}_{-0.43}^{+0.44}\pm 0.39$	${20.23}_{-1.00}^{+0.95}\pm 0.92$	${18.91}_{-0.83}^{+1.50}\pm 0.68$
	Contaminated	${20.86}_{-0.12}^{+0.12}\pm 0.10$	${20.59}_{-0.15}^{+0.17}\pm 0.18$	${21.59}_{-0.16}^{+0.18}\pm 0.11$	${19.24}_{-0.19}^{+0.18}\pm 0.16$
$i=21.5$	KI17	${19.65}_{-0.19}^{+0.19}\pm 0.12$	${18.72}_{-0.38}^{+0.43}\pm 0.33$	${20.06}_{-0.51}^{+0.64}\pm 0.37$	${18.37}_{-0.74}^{+1.58}\pm 0.58$
	Conventional	${19.08}_{-0.21}^{+0.22}\pm 0.16$	${18.60}_{-0.38}^{+0.41}\pm 0.38$	${20.33}_{-0.69}^{+0.93}\pm 0.59$	${18.95}_{-0.91}^{+1.44}\pm 0.70$
	Contaminated	${20.83}_{-0.11}^{+0.10}\pm 0.08$	${20.69}_{-0.13}^{+0.14}\pm 0.16$	${21.57}_{-0.14}^{+0.14}\pm 0.08$	${19.32}_{-0.17}^{+0.16}\pm 0.14$
$i=22$	KI17	${19.71}_{-0.16}^{+0.16}\pm 0.11$	${18.75}_{-0.28}^{+0.30}\pm 0.27$	${20.06}_{-0.37}^{+0.64}\pm 0.24$	${18.42}_{-0.79}^{+1.65}\pm 0.59$
	Conventional	${19.02}_{-0.19}^{+0.20}\pm 0.15$	${18.50}_{-0.33}^{+0.38}\pm 0.31$	${20.45}_{-0.63}^{+0.86}\pm 0.52$	${18.81}_{-0.88}^{+1.51}\pm 0.76$
	Contaminated	${20.82}_{-0.09}^{+0.08}\pm 0.06$	${20.62}_{-0.11}^{+0.12}\pm 0.13$	${21.68}_{-0.12}^{+0.13}\pm 0.07$	${19.30}_{-0.16}^{+0.15}\pm 0.13$

Equations

$$\Phi (E,\Delta \Omega )=\left[\frac{C⟨\sigma v⟩}{4\pi m_{\text{DM}}^2}\sum _f{b}_f{\left(\frac{d{N}_{\gamma }}{dE}\right)}_f\right]\times J(\Delta \Omega ).$$

$$J(\Delta \Omega )=\left[{\int }_{\Delta \Omega }d\Omega {\int }_{l.o.s.}dl{\rho }^2(l,\Omega )\right].$$

$${\sigma }_{l.o.s}^2(R)=\frac{2}{{\Sigma }_{\ast }(R)}{\int }_R^{\infty }dr\left(1-{\beta }_{\text{ani}}(r)\frac{R^2}{r^2}\right)\frac{{\nu }_{\ast }(r){\sigma }_r^2(r)}{\sqrt{1-R^2/r^2}},$$

$${\sigma }_r^2(r)=\frac{1}{{\nu }_{\ast }(r)}{\int }_r^{\infty }{\nu }_{\ast }(r^{\prime }){\left(\frac{r^{\prime }}{r}\right)}^{2{\beta }_{\text{ani}}}\frac{GM(r^{\prime })}{r^{\prime 2}}d{r}^{\prime }.$$

$${\rho }_{\text{DM}}(r)={\rho }_s(r/r_s{)}^{-\gamma }(1+(r/r_s{)}^{\alpha }{)}^{-(\beta -\gamma )/\alpha },$$

$${\nu }_{\ast }(r)=(3/4\pi r_e^3)(1+(r/r_e{)}^2{)}^{-5/2}.$$

$$-2\ln {\mathcal{L}}_{\text{s}}=-2\sum _i\ln (f_{\text{Mem}}(v_i,R_i)),$$

$$f_{\text{Mem}}(v,R)=2\pi R{\Sigma }_{\ast }(R)C_{\text{Mem}}\mathcal{G}[v;v_{\text{Mem}},{\sigma }_{l.o.s}(R)].$$

$$-2\ln {\mathcal{L}}_{\text{m}}=-2\sum _i\ln (s{f}_{\text{Mem}}(v_i,R_i)+(1-s)f_{\text{FG}}(v_i,R_i)),$$

$$f_{\text{FG}}(v,R)=2\pi R{C}_{\text{FG}}\prod _{j=1}^3\mathcal{G}[v;v_{\text{FG}j},{\sigma }_{\text{FG}j}],$$

$$N_{\text{sg}}^i=T_{\text{obs}}\times {\int }_{\Delta E_i}d{E}^{\prime }\int d{E}_{\gamma }{\mathcal{A}}_{\text{eff}}(E_{\gamma })\Phi (E_{\gamma },\Delta \Omega )\mathcal{R}(E_{\gamma },E^{\prime }).$$

$${\mathcal{L}}_j(⟨\sigma v⟩)$$

$$\pi (J)$$

Full text

††pubyear: 2017††pagerange: Foreground effect on the $J$-factor estimation of ultra-faint dwarf spheroidal galaxies–C

1 Introduction

Various astrophysical observations such as the dynamics of galaxy clusters (Zwicky, 1933), rotation curves of spiral galaxies (Rubin, Thonnard & Ford, 1978, 1980), and gravitational lensing (McLaughlin, 1999; Łokas & Mamon, 2003; Clowe et al., 2006; Bradac et al., 2006), strongly indicate the existence of dark matter (DM) in the astronomical objects. A recent global fit of the Cosmic Microwave Background (CMB), Large Scale Structure (LSS), and Supernovae (SNe) observations (Ade et al., 2016) reveal that quarter of the total energy of the universe consists of DM. One of the most attractive candidates of DM is weakly interacting massive particle (WIMP), which naturally explains the observed dark matter density with its annihilation channels into lighter standard model particles. Particularly, the WIMP dark matter with $\lesssim{\cal O}(1)$ TeV has drawn attention in the context of the physics beyond the standard model such as supersymmetry (see e.g. Jungman, Kamionkowski & Griest, 1996 also Murayama, 2007; Feng, 2010).

Gamma-ray indirect detection experiment, which aims to observe gamma-rays induced by the DM annihilation, has a strong sensitivity to this ${\cal O}(1)\,\text{TeV}$ WIMP. Among various astronomical objects, dwarf spheroidal satellite galaxies (dSphs) associated with the Milky Way are the ideal targets due to its small distance ($\sim 10\--{\text{a few hundred}}\,\text{kpc}$ from the solar system) and dense DM environment with low astrophysical background (Cholis & Salucci, 2012; Lefranc et al., 2016). However, recent studies show that expected signal flux coming from the dSphs is significantly affected by various uncertainties such as the statistical procedure (Martinez et al., 2009), DM distribution  (Geringer-Sameth et al., 2015; Bonnivard et al., 2015; Hayashi et al., 2016), stellar distribution (Ullio & Valli, 2016), unresolved binary stars ( Mateo, Olszewski & Walker, 2008; Koch et al., 2007; Minor, 2013; Simon & Geha, 2007; Simon et al., 2011; McConnachie & Cote, 2010; Koposov et al., 2011; Kirby et al., 2013; Simon et al., 2015) and foreground contamination (Bonnivard et al., 2016; Ichikawa et al., 2017).

Although future deep spectroscopic surveys would mitigate most of these systematic uncertainties, the foreground contamination remains problematic because the fraction of the foreground would not be suppressed or even become worse in the future observation.

In Ichikawa et al. (2017) (hereafter KI17), we have investigated the effects of the foreground contamination for classical dSphs. We found that, when simple color-magnitude, velocity, and surface-gravity cuts are employed to select the member stars, in which case the foreground fraction is less than 5%, the contamination can lead to an overestimation of the signal flux by a factor of $\sim 3$.

In KI17, we resolve this foreground effect by introducing a multi-component fit in which the distributions of the member stars and foregrounds are mixed.

The situation is more problematic for ultra-faint dSphs (UFDs). The UFDs were discovered after SDSS II and contains smaller number of the stars inside the system. Although recent kinematical analyses (Bonnivard et al., 2015; Ackermann et al., 2015; Geringer-Sameth et al., 2015) indicate that the signal flux coming from the UFDs can be much stronger than that from the classical dSphs, the uncertainties of these signal fluxes are much larger due to the lack of the knowledge of the kinematics inside the system. In particular, a recent study (Bonnivard et al., 2016) reveals that the foreground contamination can significantly affect the estimation of the signal flux by two orders of magnitude at most. Therefore, precise analysis of the foreground effect for the UFDs is required and will play an essential role in the future deeper spectroscopic surveys.

In this paper, we test the foreground effect for the UFDs by generating realistic stellar mock data and applying the likelihood analysis developed in KI17. We also compare the results with those obtained by the other conventional analyses. The organization of this paper is as follows. In Sec. 2, we review the formula of the gamma-ray signal flux and defines the so-called $J$-factor. In Sec. 3, we provide the procedure of our analysis. The results of the fits are given in Sec. 4. Finally, we summarize our discussion in Sec 5.

2 Signal Flux and $J$-factor

The gamma-ray signal flux of DM annihilation stemmed from the dSphs can be expressed by the following formula:

[TABLE]

The coefficient $C$ is 1/2 for Majorana and 1/4 for Dirac dark matter. Dark matter mass is defined by $m_{\text{DM}}$. The product of the total annihilation cross section $\sigma$ and the relative velocity $v$ is averaged with the velocity distribution function (represented by $\langle\dots\rangle$). The branching fraction of the annihilation channel $f$ is denoted by $b_{f}$, while the differential number density of photons from a given final-state $f$ is given by $(dN_{\gamma}/dE)_{f}$.

The factor $J$ after the parenthesis in the right-hand side (so-called $J$-factor) reflects the amount of the squared DM density inside the cone with a solid angle $\Delta\Omega$:

[TABLE]

Here we define the dark matter profile at a distance $l$ and angle $\Omega$ by $\rho(l,\Omega)$. The integration of $l$ is performed along the line-of-sight.

Currently, the most common way of the DM profile estimation is to apply dynamical mass models based on Jeans equations to the line-of-sight velocity of stars in the dSph. Under the assumption of the spherical symmetry for luminous and dark components of the dSph, the projected velocity dispersion at a projected radius $R$ can be written by

[TABLE]

where $r$ denotes the un-projected distance from the centre of the dSph, and $\Sigma_{*}(R)$ is the projected spatial stellar distribution obtained by integrating the stellar distribution $\nu_{*}(r)$ along the projected direction. The anisotropy parameter $\beta_{\text{ani}}$ is defined by $\beta_{\text{ani}}=1-(\sigma^{2}_{\theta}/\sigma^{2}_{r})$ where we define the radial, azimuthal, and polar components of the 3-dimensional dispersion curve as $\sigma_{r}$, $\sigma_{\theta}$, and $\sigma_{\phi}$, respectively, in a spherical coordinate and take $\sigma_{\theta}=\sigma_{\phi}$ for the spherical symmetry.

The radial dispersion curve, $\sigma^{2}_{r}$, is related to the gravitational potential (i.e. dark matter profile) through the spherical Jeans equation (Binney & Tremaine, 2008).

Under the assumption of constant $\beta_{\text{ani}}$, this dispersion curve can be expressed as (van der Marel, 1994; Mamon & Łokas, 2005)

[TABLE]

Here $G$ is the gravitational constant, and $M(r)$ is the enclosed mass of the spherical dark matter halo: $M(r)\equiv\int^{r}_{0}4\pi r^{\prime 2}{\rho_{\text{DM}}}(r^{\prime})dr^{\prime}$. From Eq. (3) and Eq. (4), we can estimate the DM profile $\rho_{\text{DM}}$ by constructing the dispersion curve $\sigma_{l.o.s.}(R)$ from the observational data.

As we have discussed in KI17, the dominant uncertainty of the signal flux comes from the $J$-factor. This is because while the parenthesis in Eq. (1) is well controlled by the calculation of particle physics, the estimation of the $J$-factor is limited by the number of the kinematical stellar data of the dSphs. Although the uncertainty of the $J$-factor is still under discussion,#1#1#1 This is due to the various biases in the estimation: the statistical procedure, DM halo model, stellar distribution, unresolved binaries and foreground contamination, as reviewed in KI17.

the error bar can be a few orders of magnitude larger for the UFDs, which we focus on in this work. To suppress both the statistical and systematic uncertainties, future deep spectroscopic observation is mandatory.

3 Analysis

In this section, we introduce the mock-based analysis developed in KI17. In our analysis, we first generate realistic mock dSph stellar data including foreground stars. We sample this stellar data by accounting for a spectroscopic capability, which provides realistic mock samples of a future observation. We next attempt to decrease the foreground fraction by imposing a selection rule. In this paper, we consider two approaches: naive cuts and selection by using the membership probability. Finally, we perform the halo profile estimation by using two types of the likelihood functions, which have single and mixed component(s) in their distribution function respectively. In Sec. 4, we will provide the results of the analyses by three combinations of the selections and fits: naive cut + mixed component fit, membership selection + single component fit, and naive cut + single component fit. They correspond to the KI17, conventional and the most naive approaches, respectively.

3.1 Mock dSphs

As models of the mock dSphs, we consider the four UFDs (Ursa Major II, Coma Berenices, Segue 1, and Ursa Major I), in which the observation suggests abundant DM (Hayashi et al., 2016; Bonnivard et al., 2015; Ackermann et al., 2015; Geringer-Sameth et al., 2015). We use the same DM halo profile, the stellar distribution and the domain of parameters scanned as those of Geringer-Sameth et al. (2015) based on the data provided by the kinematical observations (200, ; 201, a)#4#4#4The data of 200 was kindly provided by Josh Simon (private communication). and use the obtained DM profiles for the inputs of the dSph mocks.

In our analysis, the generalized dark matter halo density profile (Hernquist, 1990; Dehnen, 1993; Zhao, 1996) is adopted as the input dark matter profile for the mock data and fit of the likelihood analysis:

[TABLE]

where $r$ denotes the (un-projected) distance from the centre of the dSphs, and parameters $\rho_{s},\,r_{s}$ represent the typical density and scale of the halo respectively, while parameters $\alpha,\,\beta,\,\gamma$ determine the shape of the halo density profile. We also assume Plummer profile (Plummer, 1911) for the member stellar distribution:

[TABLE]

Here $r_{e}$ denotes the projected half-light radius of the dSph and we normalize the stellar distribution $\nu_{*}(r)$ to satisfy $\int 4\pi r^{2}\nu_{*}(r)dr=1$. The input parameters are shown in Table #3.#5#5#5 In order to construct the kinematical data by using the method of Cuddeford (1991) consistently, we set the range of the anisotropy $\beta_{\text{ani}}$ to be $\beta_{\text{ani}}<0$. See, e.g., Ciotti & Morganti (2010) for the limitation of the halo parameters in the analytical solution of the Jeans equation.

The mock stellar data of each dSph is constructed by assigning the colour, chemical abundance, and kinematical information. Synthetic colour-magnitude diagrams are generated by utilizing the PARSEC stellar isochrones (Bressan et al., 2012) to represent observed properties of each dSph. In detail, we first randomly draw a stellar initial mass from the Salpeter initial mass function (IMF). For that mock star, the age is drawn from an uniform distribution in the range 1010.10-1010.12 years, motivated by the fact that the UFDs analysed in this work have been reported to be dominated by an old stellar population (de Jong et al., 2008). Similarly, the value of metallicity ([Fe/H]) is drawn from a Gaussian distribution with (mean, dispersion)$=$($-2.5$, $0.3$), ($-2.5$, $0.3$), ($-2.7$, $0.7$), and ($-2.2$, $0.6$) for Ursa Major II, Coma Berenices, Segue 1, and Ursa Major I, respectively, which are approximately consistent with those estimated by Kirby et al. (2011) and Norris et al. (2010). Based on a theoretical isochrone for the given age and [Fe/H] values obtained above, the absolute magnitude, colour and surface gravity corresponding to the stellar initial mass are assigned. The apparent magnitude and observed colour are then calculated by adopting the distance modulus from 201 (b) and adding typical photometric errors increasing toward fainter magnitudes (0.012 at $i=20.0$ and 0.024 at $i=22.0$) as well as the Galactic extinction from Schlafly & Finkbeiner (2011). At this point, the star is discarded if it is fainter than the i-band limiting magnitude of 22.5. The mock stars are repeatedly generated until the number of member stars brighter than the limiting magnitude estimated by Martin et al. (2008) is reached.

The stellar IMF can deviate from the Salpeter IMF below $M\sim 0.5M_{\odot}$ (see, e.g. the Kroupa (2001) IMF). Magnitude limits ($i=21.0,\,21.5$ and $22.0$) we adopted, however, correspond to stellar masses well above 0.65 $M_{\odot}$, and thus the mock photometric data are not significantly affected by the choice of IMF.

An example of the resulting CMD is shown in Fig 1. To build 50 mock data for each dSph, the whole process is repeated 50 times by adding a Gaussian noise consistent with the uncertainty in the number of member stars estimated by Martin et al. (2008). The position and velocity of each star are assigned consistently with the input dark matter potential using the method of Cuddeford (1991) with the assumptions of the constant velocity anisotropy and spherical distribution.

The non-member stars belonging to the Milky Way galaxy are also included, which are generated by the Besançon model (Robin et al., 2003) as shown by blue points in Fig 1.#6#6#6 We have not included photometric errors for foreground stars, for the number of the foreground stars above the limiting magnitudes is less sensitive to the inclusion of the photometric errors than the case of the dSph member stars.

The Besançon model generates a set of stars without positional information of individual stars for a given direction. In this paper, we have used the same set of the foreground stars (generated by the Besançon model) for all the 50 mock data, while generated the random distributions of the stars on the sky plane in individual mocks.

We note that an appearent sharp cut-off at $(g-i)<0.3$ seen only for the foreground stars reflects different assumptions about the stellar evolution for horizontal branch stars, which is highly model dependent. Fraction of such blue stars in both the mock member and foreground stars are minor compared to red-giant stars or the main-sequence stars in the mock data.

3.2 Spectrograph

In our analysis, we adopt the same detector capability in KI17 (see Table 3 in KI17). The observing parameters are based on the capability of the Prime Focus Spectrograph (PFS) attached to 8.2 m Subaru telescope. PFS is the next generation spectrograph of the SuMiRe project (Takada et al., 2014; Sugai et al., 2015; Tamura et al., 2016) and the science operation is planned to start around $2019\--2020$. The key advantages of PFS are its large field-of-view ($\sim 1.38^{\circ}$ diameter), 2394 fibers, and the wide wavelength coverage ($380\--1260$ nm) mounted on the large aperture telescope.

One of the main targets is the classical dSphs (Fornax, Sculptor, Draco, Ursa Minor, and Sextans), for which line-of-sight velocities of stars are measured with a precision dv of $\sim 3\,\text{km/s}$ down to magnitudes deeper than $i\sim 21$ covering a wide area well beyond their tidal radii.

The unique capability of PFS has also an advantage in observing ultra-faint dwarf galaxies, increasing the sample size by a factor of 2 or more and simultaneously covering the target galaxy and the foreground/background Milky Way stars. The latter aspect is crucial in efficiently taking the effect of contaminating stars into account as in the analysis presented later in this paper.

To take the spectroscopic capabilities into account, we smear the

mock velocity, surface gravity ($\log g$) and metallicity ([Fe/H]) data with widths corresponding to the expected measurement errors, 3 km/s, 0.5 dex and 0.5 dex, respectively

and select the stars which locate at $r<d\sin\theta_{\text{ROI}}$, reflecting the limitation of the region of interest. Here $d$ denotes the distance of each dSph and $\theta_{\text{ROI}}$ is the angular radius of the region of interest.

The depth of the survey depends on the exposure time. We adopt three cases of the upper bound of the magnitude ($i_{\text{max}}=21,\,21.5$ and $22$). In the first case, we demonstrate the current sensitivity reach.#7#7#7 Since the size of the UFDs is smaller than that of the classical dSphs, the kinematical data provided by the current observations is deeper than the classical dSphs ($i\sim 19.5$).

The second case ($i_{\text{max}}=21.5$) is for a deeper survey with an integration time of several nights. The third case is for an ultimate reach.

3.3 Data selection

Before the likelihood analysis, the foreground contamination in the mock data can be largely reduced by using the information of its position, velocity, surface gravity, metallicity, and colour-magnitude. We here adopt two approaches to the data reduction: naive cut approach and more sophisticated membership selection.

3.3.1 Naive cut

In this approach, we impose the cuts of the velocity, surface gravity, metallicity, and colour-magnitude on the dataset and optimize them by (roughly) tuning the boundaries of the cuts by eye. The velocity cut is a $\pm 60$ [km/s] range from each bulk velocity $v_{\text{dSph}}$. The lower and upper bounds of the surface gravity $g$ and metallicity [Fe/H] are given in Table 2 for each dSph, while the region of the colour-magnitude diagram is shown in Fig. 1. Note that we choose these boundaries to include most of the stars in clumps. Although harder cuts can be imposed to reduce the fraction of the foreground stars, the cut eliminates scattered member stars and the reconstructed velocity distribution can be distorted and derive a bias of the halo estimation. We provide the numbers of the member and foreground stars after the cuts in the ‘Naive cut’ column in Table 3.

3.3.2 Membership selection

The latter strategy utilizes the membership probability of each star. The membership probability is defined by the probability to find a member star at a given position, velocity, surface gravity, and, metallicity. We calculate this membership probability by a conventional approach given by Walker et al. (2009). In the calculation, the distributions of the foreground stars are also taken into account. The distributions of the velocity, surface gravity, and metallicity except for the foreground velocity distribution are modeled by $R$-independent single Gaussians, while the foreground velocity distribution is fixed without free parameters and the spatial distributions are more generally parametrized. The detail of this process is given in Walker et al. (2009) and the appendix of KI17. We select stars within $95$% confidence level of the membership probability. We provide the numbers of the member and foreground stars after this selection in the ‘Membership selection’ column in Table 3. Compared with the case of the ‘Naive cut’, much higher purities of the data is obtained by this procedure, while some fraction of the member star is eliminated. We here stress that this approach assumes the constant velocity dispersion in the membership probability assignment. The member stars eliminated in the selection are mostly due to the constant velocity-dispersion bias and therefore can affect the estimation of the $J$-factors.

3.4 Kinematical fit

In this section, we provide two types of the analysis for the kinematical fit. The first one is the single component fit, in which all the data is regarded as member star. The second one is the mixed component approach developed in KI17, in which the member and foreground distribution are simultaneously fitted. For both the fits, we apply the unbinned likelihood analysis to the halo estimation.

3.4.1 Single Component fit

The single component fit is performed by assuming that the data used for the fit contains only member stars, which implies that the likelihood function is given by

[TABLE]

where $f_{\text{Mem}}(v,R)$ is the distribution function of the member stars. The index $i$ runs all the stars in the mock data set. We assume that the velocity distributions of the member stars can be approximated by a single Gaussian and hence the distribution functions can be expressed as

[TABLE]

Here ${\cal G}[x;\,\mu,\sigma]$ denotes the Gaussian distribution of a variable $x$ with a mean value $\mu$ and a standard deviation $\sigma$. We note that the parameter $v_{\text{Mem}}$ represents the bulk velocity of the dSph and mostly converges to the input bulk velocity $v_{\text{dSph}}$. The distribution functions are normalized by $C_{\text{Mem}}$ to satisfy $\int^{r_{\text{ROI}}}_{0}dR\int^{v_{\text{upper}}}_{v_{\text{lower}}}dvf_{\text{Mem}}(v,R)=1$ where $r_{\text{ROI}}\equiv d\,\sin\theta_{\text{ROI}}$.

3.4.2 Mixed Component fit

In the mixed component fit, the stellar distribution is considered to be the sum of the foreground and member star distribution. The likelihood function ${\cal L}_{\text{m}}$ is defined by introducing the membership fraction parameter $s$ as follows

[TABLE]

where $f_{\text{FG}}(v,R)$ is the distribution function of the foreground stars. We model the foreground distribution function by the production of the three Gaussians, corresponding to the foreground thin disc, thick disc, and halo components:

[TABLE]

with $v_{\text{FG}j}$, $\sigma_{\text{FG}j}$ ($j=1,\,2,\,3$) being parameters of the distribution. Here we assume that the parameters $\sigma_{\text{FG}j}$ are independent of $R$ in contrast to the dispersion of the member star. The constant $C_{\text{FG}}$ denotes the normalization factor to satisfy $\int^{r_{\text{ROI}}}_{0}dR\int^{v_{\text{upper}}}_{v_{\text{lower}}}dvf_{\text{FG}}(v,R)=1$.

For the sake of the convergence of the mixed component fit, we constrain the parameters $v_{\text{FG}j}$, and $\sigma_{\text{FG}j}$ by using the data in the control region (i.e., the region in which the number of the member stars is negligible). In KI17, we deduce the foreground velocity distribution by using the data out of the region of the velocity cut and interpolate it to the signal region. For the UFD case, on the other hand, since the bulk velocities of these dSphs are not as large as that of classical dSphs, the foreground estimation by using the control region in the velocity distribution does not efficiently work. In the case of the classical dSphs considered in KI17, the bulk velocity of the member stars is largely different from the Milky Way stellar halo component, which makes the velocity cut work most efficiently. As Fig. 2 shows, on the contrary, the distribution of the foreground stars overlaps with the distribution peak of the UFD, hence the significant amount of the foreground stars contaminate the signal region, which prevent using the same method in KI17.

Instead, we define the control regions in the distribution of the spatial position by setting an annulus centred at each dSph galaxy from the radius of the signal region to the PFS threshold, $\theta=0.65^{\circ}$.

For precise determination of the DM profile and increasing the value of the $J$-factor, we need to take a large signal region as possible, though the signal region must be included in PFS threshold and the range for the control region must be reserved.

Here, the radii of the signal regions are chosen to be $2r_{e},\,4r_{e},\,4r_{e}$, and $r_{e}$ for Ursa Major II, Coma Berenices, Segue 1, and Ursa Major I respectively, based on their half-light radii $r_{e}$. We should bear in mind that the radius of the signal region is not optimized. Instead, we tried to take the radius so that its angular diameter is about a half of FoV.

When we perform a fit to the control region, we take into account the effect of the thin and thick disc components of the foreground stars in addition to the halo component, because the disc components remain after surface gravity and metallicity cuts in case of UFDs. This contrasts to the case of classical dSphs, where the foreground stars mainly belong to the halo component after the naive cut. In order to represent the three foreground components we assume the foreground distribution can be expressed by a sum of three Gaussian functions. We first perform fits by the three Gaussian model for control region data on which colour-magnitude, $i_{\text{max}}$ and ROI cuts are imposed, and obtain the best-fitting values and standard deviations of each Gaussian. Then we perform secondary fits for the control region data with all naive cuts (colour-magnitude, $i_{\text{max}}$, ROI, surface gravity and [Fe/H]) imposed on, using the best-fitting Gaussians achieved in the first fit as the priors. Here we obtain the best-fitting values and standard deviations of $v_{\text{FG}j}$, $\sigma_{\text{FG}j}$, which are defined as $v_{\text{FG0}j}$, $\sigma_{\text{FG0}j}$, $dv_{\text{FG}j}$, and $d\sigma_{\text{FG}j}$ respectively. Finally we use this information as a prior for $v_{\text{FG}j}$, $\sigma_{\text{FG}j}$ by multiplying $\prod_{j=1,2,3}{\cal G}[v_{\text{FG}j};\,v_{\text{FG0}j},\,dv_{\text{FG}j}]\,{\cal G}[\sigma_{\text{FG}j};\,\sigma_{\text{FG0}j},\,d\sigma_{\text{FG}j}]$ to the likelihood function ${\cal L}$ in Eq.(9).

3.4.3 Fit algorithm

The likelihood function (multiplied by the foreground priors for the mixed component fit) are searched by performing the Metropolis-Hastings algorithm (Metropolis et al., 1953; Hastings, 1970) of the Markov Chain Monte Carlo (MCMC) method. The parameter set of the single component fit consists of the five free parameters of the dark matter halo ($\rho_{s},\,r_{s},\,\alpha,\,\beta,\,\gamma$), one velocity anisotropy parameter $\beta_{\text{ani}}$ and one nuisance parameters ($v_{\text{Mem}}$), while the mixed component fit also has the other seven nuisance parameters ($s$, $\,v_{\text{FG}j}$, $\sigma_{\text{FG}j}$). In the MCMC method, the halo parameters are searched under the flat/log-flat priors within the range of $-4<\log_{10}(\rho_{s}/[\text{M}_{\bigodot}/\text{pc}^{3}])<4$, $-2<\log_{10}(r_{s}/[\text{kpc}])<5$, $0.5<\alpha<3$, $3<\beta<10$, $0<\gamma<1.2$ and $-1<\log_{10}(1-\beta_{\text{ani}})<1$.

3.5 Strategy

Using 50 mocks for each case ($i_{\text{max}}=21,\,21.5$, and 22), we test three types of the $J$-factor estimation: the method of KI17 (naive cut + mixed component fit), Conventional analysis (membership selection + single component fit), and Contaminated fit (naive cut + single component fit). We here stress that in the KI17 approach, the velocity distribution of the foreground is parametrized by the fit and therefore the error bar of the $J$-factor involves the uncertainty of the foreground distribution, while we fix the spatial stellar distributions of member and foreground stars in the likelihood. This contrasts with the Conventional approach in which a fixed model of the foreground velocity distribution and parametrized spatial distributions are used in the selection.

4 Results

Fig. 3 shows the results of these three approaches, namely, the method of KI17, the Conventional analysis, and the Contaminated fit by blue, orange, and green bars, respectively. Here we give the averaged median values of $\log_{10}(J/[\text{GeV}^{2}/\text{cm}^{5}])$ for each fit by the dots. The lighter error bars show the averages of the widths of the $68\,\%$ quantiles, while the darker ones show the square roots of the $68\,\%$ quantiles and the standard deviations of the median values, written in an additional way to the lighter ones. The grey dashed lines show the input values. For each dSph, three bars with the same colours correspond to the case of $i_{\text{max}}=21,\,21.5$, and $22$ with $\theta_{\text{ROI}}=0.65$ respectively, from the left. All $J$-factors are calculated within an angular radius of $0.5\,\text{degree}$ (i.e., $\Delta\Omega=2.4\times 10^{-4}\,\text{sr}$), which is the standard size for the $J$-factor calculation. We here choose the distance from the centre of the dSph to the outermost observed member star $r_{\text{max}}$ as the most conservative radius given by  Geringer-Sameth et al. (2015).

In the Contaminated analysis (green bars in Fig 3), the overestimation of the $J$-factor becomes more than an order of the magnitude. This is because the dispersion curve inflates due to the foreground contamination which mainly locate at the outer region with a large velocity dispersion ($\sim 30\--40$ km/s). Since the fraction of the foreground contaminating stars is more than 50 %, the overestimation is much larger than that of the classical case.

On the other hand, the $J$-factors seem to be successfully reproduced by the Conventional approach (orange bars in Fig 3). However, since this approach assumes a constant velocity dispersion in its membership selection, the dispersion curve after the selection tends to be more or less constant as a function of the radius, which lead to a small bias to the $J$-factor estimation. We provide the typical uncertainty in the dispersion curve and its sample-to-sample scatter for the 50 mocks in the left column of Fig. 4. The red lines show the median value of the dispersion curves obtained by the fit of the Conventional approach (averaged by the 50 mocks), while the green band shows the (averaged) 68 % quantile. The median values of the dispersion curves also fluctuate sample by sample, reflecting the quality of the sample. We show this fluctuation by the orange shaded regions which are obtained by the square root sum of the standard deviation of the median values of the 50 mocks and the 68 % quantiles. The input dispersion curves are also shown by the grey dashed lines.

For the Ursa Major II, the figure shows that the dispersion curve is flatter than that of the input, which makes the $J$-factor underestimated. This fact seems to be caused by the constant velocity dispersion bias. This effect becomes more significant for a larger size of stellar data, as can be seen in the three orange bars. Meanwhile, since the changes of the dispersions curves of the other dSphs are not as large as the Ursa Major II case, the effect of the bias is not seen.

We also note the results for the Ursa Major I case. Although the number of the stars in Ursa Major I does not significantly differ from the other dSphs (see Table 3), both the Conventional and KI17 approaches cannot determine the $J$-factor as precisely as those of the other dSphs. It seems to originate in the fact that the inner part of the dark matter profile is not well determined. We left this analysis to future work.

KI17 approach (blue bars in Fig. 3) also provides successful $J$-factor estimations. For the Ursa Major II case, the $J$-factor is getting converged to the input value for deeper observations, while it cannot be seen in the conventional approach. For other dSphs, the result of our analysis is compatible with the conventional one, which reflects the fact that the velocity dispersion at the inner part is more or less constant.

We also give the distribution of the dispersion curve obtained by the fit in the right column of Fig. 4. Compared with the distribution of the Conventional approach, the width of the 68 % quantile is larger in the outer region, while it becomes smaller in the inner region (except for the Ursa Major I case). The results of the $J$-factor estimation implies that the width in the inner region preferentially affects the uncertainties of the $J$-factors. We also note that the median dispersion curve of the Ursa Major II case successfully follows the input curve at $R\sim 250$ pc, in contrast with that of the Conventional approach. The sensitivity of the indirect dark matter detection observing gamma-rays from the dSphs depends directly on the median values and uncertainties of $J$-factors. See Appendix B for those who are interested in this dependence in a concrete example.

5 Summary

In this paper, we have investigated the effect of the foreground contamination on the estimation of astrophysical factor, using the mock kinematical data of the four representative ultra-faint dwarf spheroidal galaxies. This is because we cannot completely distinguish the foreground stars from the dSph’s member stars even if imposing several data cuts. We have adopted our developed fitting analysis, KI17, utilizing the future spectroscopic survey, PFS. Such a multi-object spectrograph with large field of view enables us to observe numerous number of stellar spectra required for the KI17 analysis.

For comparison, we have performed three types of the $J$-factor estimation: the KI17 methods, the Conventional analysis and the Contaminated fit. As the result of the analysis, the $J$-factor value estimated by the Contaminated analysis is up to a few hundred times larger than the input value and its confidence interval is significantly small, because all stellar data after naive cut are regarded as member star even including the foreground contamination. On the other hand, the KI17 and Conventional analysis can reproduce the input $J$-factor value within $1\sigma$ confidence levels except for Ursa Major II.

For the case of Ursa Major II, which has the non-flat velocity dispersion curve, the Conventional approach underestimates the $J$-factor value with respect to the input value. This seems to originate from the assumption of the constant velocity dispersion at the membership selection. Accordingly, the member stars in the outer region tend to be rejected.

The likelihood function of the KI17 method includes the information of both the foreground stars and the member ones together with the parameters describing their distribution functions, and the properties of the foreground stars are roughly determined by the photometric and spectroscopic observations of the stars in the control region. This method allows us to treat correctly and statistically the effect of the foreground contamination for the observational data. Moreover, the method can provide the validation of the assumption concerning the velocity dispersion curve in the Conventional analysis. Therefore, our statistical method should become powerful tool for the $J$-factor estimate of the MW dSphs in the PFS-era.

Acknowledgments

We are grateful to the referee for her/his careful reading of our paper and thoughtful comments. We would like to give special thanks to Josh Simon and Marla Geha for giving us the kinematic data of UFD galaxies. This work is supported in part by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT), Japan, No. 16H01090 (for K. H.), No. 15H05889, No. 25105011 (for M. I.), No. 16H02176, No. 26104009 (for S. M.), No. 17H02878 (for S. M. and M. I.), No. 17K14249 (for M. N. I.), and No. 18J21186 (for S. H.). The work of K.I. is also supported by the JSPS Research Fellowships for Young Scientists. Finally, Kavli IPMU is supported by World Premier International Research Center Initiative (WPI), MEXT, Japan.

Appendix A Table of obtained $J$-factors

We show the numerical values of the $J$-factors and its errors obtained by our fits in Table 4.

Appendix B Implication to gamma-ray detections

In this appendix, we discuss the sensitivity of the future indirect dark matter detection observing gamma-rays from dSphs, in particular, focusing on how it depends on the median values and uncertainties of J-factors obtained in section 4. In order to make our analysis concrete, we consider the Cherenkov Telescope Array project (Cherenkov Telescope Array Observatory gGmbH, 2017) and the wino DM as a gamma-ray observatory and a DM candidate, respectively.

Wino DM is one of the most attractive WIMP candidates, where it is introduced in the supersymmetric extension of the standard model as the superpartner of the neural weak gauge boson. The wino DM attracts many attentions at present, because it is predicted in the anomaly-mediated supersymmetry breaking scenario explaining the Higgs boson mass of 125 GeV as well as non-observation of new physics signals at collider experiments (Ibe et al., 2012). Its mass is predicted to be O(1) TeV and its annihilation cross section is boosted by the so-called the Sommerfeld effect (Hisano et al., 2004), so that the indirect dark matter detection is expected to work very efficiently to detect the wino DM. We refer Lefranc et al. (2016) and Cirelli et al. (2011) for the branching fraction of each annihilation channel and the corresponding fragmentation function $(dN_{\gamma}/dE)_{f}$, respectively, to calculate the photon flux in Eq. (1).

CTA will provide significantly improved sensitivity to WIMP DM with its high angular resolution and its wide energy coverage. Here we briefly review the CTA analysis in Lefranc et al. (2016). The number of photon count in $i$-th energy bin $N^{i}_{\gamma}$ is given by $N^{i}_{\gamma}=N^{i}_{\text{sg}}+N^{i}_{\text{bg}}$. Here $N^{i}_{\text{sg}}$ denotes the number of the signal photon in the $i$-th energy bin $\Delta E_{i}$, given by

[TABLE]

Here, $T_{\text{obs}}$ is the observation time for each UFD, which is assumed to be $T_{\text{obs}}$ = 50 hours as a benchmark. The differential flux $\Phi$ is the one obtained by Eq. (1), The effective area $\mathcal{A}_{\text{eff}}(E_{\gamma})$ and the energy resolution of $\mathcal{R}(E_{\gamma},E^{\prime})\equiv\mathcal{G}[E^{\prime};E_{\gamma},\sqrt{8\ln(2)}\delta_{\text{res}}(E_{\gamma})]$ are given in Cherenkov Telescope Array Observatory gGmbH (2017). The background rate, $\nu_{i}$, is also given in Cherenkov Telescope Array Observatory gGmbH (2017), from which we obtain the number of the background photon $N_{\text{bg}}^{i}$ by multiplying $T_{\text{obs}}$. In our analysis, the size of each energy bin is the same as that of the background rate in the reference.

The likelihood function of the indirect detection for a specific dSph (indexed with $j$) is given by

[TABLE]

where $N^{i}_{\text{obs}}$ denotes the photon number observed in $i$-th energy bin. To estimate the mean sensitivity of the annihilation cross section for the null observation, we assume $N^{i}_{\text{obs}}={N^{i}_{\text{bg}}}$. The factor $\pi(J)$ represents the uncertainty of the $J$-factor which is discussed in this paper. Total likelihood function is given by the product of each likelihood function $\mathcal{L}(\langle\sigma v\rangle)=\prod_{j}\mathcal{L}_{j}(\langle\sigma v\rangle)$. Then we perform the statistical test with the condition $\chi^{2}\equiv-2\ln[\mathcal{L}(\langle\sigma v\rangle)/\mathcal{L}_{\text{max}}]<2.71$ (95% confidence level). The sensitivity line, i.e. the prospected upper limit on $\langle\sigma v\rangle$, is obtained by $\chi^{2}=2.71$ as the function of $M_{\text{DM}}$.

The upper panel of Fig. 5 shows the sensitivity lines of photon cross section $\langle\sigma v\rangle_{\gamma\gamma}$ for the case of $i_{\text{max}}=21.5$. For the illustration, we convert the sensitivity lines $\langle\sigma v\rangle$, obtained by $\chi^{2}=2.71$, into $\langle\sigma v\rangle_{\gamma\gamma}=b_{\gamma\gamma}\langle\sigma v\rangle$.#8#8#8 For the ”pure” photon model where $b_{\gamma\gamma}^{\text{pure}}=1$, for instance, its sensitivity $\langle\sigma v\rangle_{\gamma\gamma}^{\text{pure}}$ is obtained by $(1+b_{\gamma Z}^{\text{Wino}}/2b_{\gamma\gamma}^{\text{Wino}})\langle\sigma v\rangle_{\gamma\gamma}^{\text{Wino}}$, where $\langle\sigma v\rangle_{\gamma\gamma}^{\text{Wino}}$ is the sensitivity of Wino DM in Fig. 5. This reinterpretation is verified because the continuum spectrum of Wino DM barely affect to its sensitivity lines.

The sensitivity line obtained by the Contaminated approach is about 100 times severer than the other methods. The sensitivity lines simply reflect the estimated value of the $J$-factor, which show the importance of a careful estimation of the $J$-factors, since otherwise the dark matter model will be constrained too aggressively. The Conventional and the KI17 approaches avoid such a problem and, when all of the four dSphs are taken into account in the likelihood, there is no significant difference between them at the level of the present observational depth. We here note that in case of the observation of Segue 1 only, even which is the most promising target due to its large $J$-factor, the difference between the Conventional and the KI17 approaches is more significant. The difference will also appear when the $J$-factors are estimated at deeper observation, as can be expected from the $J$-factors in Fig. 3. In the lower panel of Fig. 5 we show the improvement of the sensitivity with the KI17 method by increasing the observation depth. From the view point of the thermal Wino dark matter search where its mass is predicted to be about 3 TeV  (Hisano et al., 2007), it will be crucial to choose the observational depth at around $i_{\text{max}}=21.5$, as can be seen in the panel.

Appendix C Effect of the truncation radius

The truncation radius dependence of the $J$-factor is shown in Fig. 6 for each $i_{\text{max}}=21,21.5$, or $22$ at the four dSphs. The figure shows that the both Conventional and KI17 analyses reproduce the input $J$-factor well for any $R_{\text{trunc}}$. #9#9#9The only exception is seen for the case of Ursa Major II at the Conventional analysis (top right panel). Moreover, it can also be seen that the value of $J$-factor is not alter significantly when $r_{\text{max}}$ is similar to or larger than $r_{\text{max}}$. We hence set the truncation radius of the dark matter halo to $r_{\text{max}}$ in Table 2 in our analysis of this paper.