Contribution of Dark Matter Annihilation to Gamma-Ray Burst Afterglows near Massive Galaxy Centers

The existence of dark matter (DM) is strongly supported by convincing evidence in cosmology and astrophysics (e.g., Bertone et al. 2005). Numerous astronomical processes and evolutions are partly or predominantly affected by DM. However, the nature of DM remains a mystery, and its direct detection has yet to be achieved. The mainstream DM model invokes DM particles, the most promising candidate being the weakly interacting massive particles (WIMP; see, e.g., Jungman et al. 1996; Bertone et al. 2005; Bergström 2012). The WIMP model can be tested by detecting potential signals from WIMP annihilations, which can ultimately produce Standard Model particles, such as neutrinos, photons, electrons, and positrons. Searching for gamma-ray emissions produced by DM annihilation is one method to indirectly detect DM; for this task, nearby galaxy centers and dwarf spheroidal galaxies are appealing targets (e.g., Ackermann et al. 2015; Abdallah et al. 2018; Johnson et al. 2019). Moreover, the radio detection of synchrotron radiation induced by electrons and positrons produced by DM annihilation in galaxies is considered a promising method to constrain the particle nature of DM (e.g., Egorov & Pierpaoli 2013; Storm et al. 2013). However, such emissions are generally weak, and possible detections are limited to galaxies in the local universe.

If significant signals can be powered by electrons produced by DM annihilation, an accelerator should be established in a high-density DM halo. In addition, the number density of system-provided electrons should be comparable to that of DM electrons (DMEs). The jets generated from gamma-ray bursts (GRBs) near galaxy centers are ideal accelerators.

GRBs are the most luminous explosions in the universe. The ultrarelativistic jets in the line of sight launched by newly born magnetars or hyperaccreting black holes can trigger observable GRBs (see the reviews by Mészáros 2006; Zhang 2007, 2018; Liu et al. 2017). The prompt gamma-ray emissions and multiband afterglows of GRBs generally originate from the internal and external shock phases, respectively (e.g., Mészáros & Rees 1993; Rees & Mészáros 1994; Mészáros & Rees 1997; Sari et al. 1998; Piran 2004). Most GRB afterglows can be explained by synchrotron and synchrotron self-Compton (SSC) processes under the conditions of external shocks sweeping the interstellar medium (ISM, or interstellar wind) and accelerating electrons. This explanation has been thoroughly verified by multiwavelength observations of GRB 190114C (MAGIC Collaboration et al. 2019; Fraija et al. 2019; Wang et al. 2019). In the prompt emission phase, the number of DMEs is much less than that of electrons emitted from the central engine, so the electromagnetic radiation of DMEs might be observable only in the afterglow phase.

In this paper, the lightest supersymmetric particle in the WIMP model, namely, the neutralino, is chosen; this particle has four different annihilation channels, i.e., w⁺ w⁻, $b\bar{b}$, μ⁺ μ⁻, and τ⁺ τ⁻. We then study the acceleration and radiation mechanisms of DMEs in external shocks and subsequently predict and analyze the DM contributions to the light curves and spectra of GRB afterglows. In Section 2, we provide a detailed description of the method employed to calculate the light curves and spectra of GRB afterglows, including the effects of DM. In Section 3, the results using different DM parameters are shown. In Section 4, we analyze the locations of GRBs in massive galaxies associated with the effects of DM. The conclusions and a discussion are presented in Section 5.

2.1. Dynamical Evolution of the External Forward Shock

As the GRB ejecta launched from the central engine interact with the circumburst medium, a relativistic shock is generated that propagates through the medium. We adopt an approximate dynamical evolution model to discuss the evolution of the external forward shock (Huang et al. 1999, 2000)

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Gamma }}}{{dm}}=-\displaystyle \frac{{{\rm{\Gamma }}}^{2}-1}{{M}_{\mathrm{ej}}+\epsilon m+2(1-\epsilon ){\rm{\Gamma }}m},\end{eqnarray} \tag{ 1 }$$

where Γ, M_ej, m, and are the Lorentz factor, the ejecta mass, the swept mass from the external medium, and the radiation efficiency of the external shock, respectively. The problem can be solved by introducing two more differential equations

$$\begin{eqnarray}&&\displaystyle \frac{{dm}}{{dR}}=2\pi {R}^{2}(1-\cos {\theta }_{{\rm{j}}}){{nm}}_{{\rm{p}}}\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{dR}}{{dt}}=\beta c{\rm{\Gamma }}({\rm{\Gamma }}+\sqrt{{{\rm{\Gamma }}}^{2}-1}),\end{eqnarray} \tag{ 3 }$$

where ${\theta }_{{\rm{j}}}$ is the half opening angle of the ejecta, n is the number density of the circumburst medium, t is the time measured in the observer frame, m_p is the proton mass, and $\beta \,=\sqrt{1-1/{{\rm{\Gamma }}}^{2}}$ is the ejecta velocity. Generally, the circumburst medium can be classified into two cases: an ISM and a stellar wind. The mass density of the circumburst medium is a constant in the former case and declines with ∼R⁻² in the latter case. In this work, we consider only the ISM case. We also neglect the evolution of θ_j since numerical simulations showed that sideways expansion is not important (Zhang & MacFadyen 2009).

2.2. Electron Distribution

With the contribution of DMEs considered, the evolution of the shock-accelerated electrons can be expressed as a function of the radius R, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}R}\left(\displaystyle \frac{{{dN}}_{{\rm{e}}}^{{\prime} }}{d{\gamma }_{{\rm{e}}}^{{\prime} }}\right)+\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}{\gamma }_{{\rm{e}}}^{{\prime} }}\left({\dot{\gamma }}_{{\rm{e}}}^{{\prime} }\displaystyle \frac{{{dt}}^{{\prime} }}{{dR}}\displaystyle \frac{{{dN}}_{{\rm{e}}}^{{\prime} }}{d{\gamma }_{{\rm{e}}}^{{\prime} }}\right)={\hat{Q}}_{{\rm{ISM}}}^{{\prime} }+{\hat{Q}}_{{\rm{DM}}}^{{\prime} },\end{eqnarray} \tag{ 4 }$$

where ${{dN}}_{{\rm{e}}}^{{\prime} }/{d\gamma }_{{\rm{e}}}^{{\prime} }$ is the instantaneous electron energy spectrum, ${\gamma }_{{\rm{e}}}^{{\prime} }$ is the Lorentz factor of the shock-accelerated electrons in the comoving frame of the shock, ${dt}^{\prime} /{dR}$ = 1/Γc (with t' being the time in the shock comoving frame and Γ being the Lorentz factor of the external forward shock), ${\dot{\gamma }}_{{\rm{e}}}^{{\prime} }$ is the cooling rate of electrons with the Lorentz factor ${\gamma }_{{\rm{e}}}^{{\prime} }$, and ${\hat{Q}}_{{\rm{ISM}}}^{{\prime} }{dR}$ and ${\hat{Q}}_{{\rm{DM}}}^{{\prime} }{dR}$ represent the injection of electrons from the ISM and DM into the shock during its propagation from R to R + dR.

Here, ${\hat{Q}}_{{\rm{ISM}}}^{{\prime} }=\bar{K}{\gamma }_{{\rm{e}}}^{{\prime} -p}$ with $\bar{K}\approx 4\pi (p-1){R}^{2}n{\gamma }_{{\rm{e}},{\rm{\min }}}^{{\prime} p-1}$ is adopted to describe the injection behavior of newly shocked circumburst medium electrons (CMEs) in the ISM, where p (>2) is the power-law index and ${\gamma }_{{\rm{e}},\min }^{{\prime} }\leqslant {\gamma }_{{\rm{e}}}^{{\prime} }\leqslant {\gamma }_{{\rm{e}},\max }^{{\prime} }$ is adopted for ${\gamma }_{{\rm{e}}}^{{\prime} }$.

The shocked electrons and the magnetic fields share the fractions _e and _B of the thermal energy density in the forward shock downstream. Since DMEs are involved in our work, the minimum Lorentz factor of the shock-accelerated electrons can be expressed as ${\gamma }_{{\rm{e}},{\rm{\min }}}^{{\prime} }={\epsilon }_{{\rm{e}}}\eta ({\rm{\Gamma }}-1)(p-2){m}_{{\rm{p}}}/(p-1){m}_{{\rm{e}}}$, where η (≤1) is the ratio of the number density of CMEs to the number density of all electrons, including both CMEs and DMEs, and m_e denotes the electron masses. The maximum Lorentz factor of electrons is ${\gamma }_{{\rm{e}},\max }^{{\prime} }=\sqrt{9{m}_{{\rm{e}}}^{2}{c}^{4}/8B^{\prime} {e}^{3}}$, with ${B}^{{\prime} }=\sqrt{32\pi {\rm{\Gamma }}({\rm{\Gamma }}-1){{nm}}_{{\rm{p}}}{\epsilon }_{B}{c}^{2}}$ being the magnetic field behind the shock (e.g., Kumar et al. 2012).

The evolution of the DME fluid is described by the following diffusion-loss equation (neglecting reacceleration and advection effects, see, e.g., Colafrancesco et al. 2006; Borriello et al. 2009), i.e.,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{{{dn}}_{{\rm{e}},\chi }}{d{\gamma }_{{\rm{e}},\chi }}={\boldsymbol{\nabla }}\cdot \left[K({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}}){\boldsymbol{\nabla }}\displaystyle \frac{{{dn}}_{{\rm{e}},\chi }}{d{\gamma }_{{\rm{e}},\chi }}\right]\\ \,+\,\displaystyle \frac{\partial }{\partial {\gamma }_{{\rm{e}},\chi }}\left[b({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}})\displaystyle \frac{{{dn}}_{{\rm{e}},\chi }}{d{\gamma }_{{\rm{e}},\chi }}\right]+Q({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}}),\end{array}\end{eqnarray} \tag{ 5 }$$

where ${{dn}}_{{\rm{e}},\chi }/d{\gamma }_{{\rm{e}},\chi }$ is the DME equilibrium spectrum, $K({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}})$ is the diffusion coefficient, and $b({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}})$ stands for the energy loss rate. Since GRBs are stellar-scale events, we reasonably assume that the values of the DM densities in these regions are constants. Then the gradients of the DME densities should be considered as 0, so the first term on the right-hand side of the above equation can be neglected. Thus, the steady-state energy spectrum for the DMEs in the ISM can be expressed as (e.g., Colafrancesco et al. 2006; Borriello et al. 2009)

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{{\rm{e}},\chi }({\gamma }_{{\rm{e}},\chi })}{d{\gamma }_{{\rm{e}},\chi }}=\displaystyle \frac{1}{b({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}})}{\int }_{{\gamma }_{{\rm{e}},\chi }}^{{\gamma }_{{\rm{e}},\chi ,\max }}Q(\zeta ,{\boldsymbol{r}})d\zeta ,\end{eqnarray} \tag{ 6 }$$

where ${\gamma }_{{\rm{e}},\chi ,\max }={M}_{\chi }/{m}_{{\rm{e}}}$ is the maximum Lorentz factor of the DMEs, ${M}_{\chi }$ is the DM mass, and $Q({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}})$ is the source item, which can be expressed as

$$\begin{eqnarray}&&Q({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}})=\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{\rho }_{\chi }({\boldsymbol{r}})}{{M}_{\chi }}\right)}^{2} \langle {\sigma }_{A}\upsilon \rangle \displaystyle \frac{{{dN}}_{{\rm{inj}}}}{d{\gamma }_{{\rm{e}},\chi }},\end{eqnarray} \tag{ 7 }$$

where ρ_χ(r) is the DM density profile, $\langle {\sigma }_{A}\upsilon \rangle$ is the annihilation cross section that has a typical value ($3.0\times {10}^{-26}\,{\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}$ is used in this work), and ${{dN}}_{\mathrm{inj}}/d{\gamma }_{{\rm{e}},\chi }$ is the DME injection spectrum, which can be obtained by the Dark SUSY package (e.g., Gondolo et al. 2004; Bringmann et al. 2018). In addition, the energy loss term, $b({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}})$, involving the energy loss of synchrotron radiation and the inverse Compton scattering of the cosmic microwave background (CMB) and starlight photons, which are the faster energy loss processes for driving the equilibrium of DMEs, can be expressed by

$$\begin{eqnarray}\begin{array}{rcl}b({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}}) & = & {b}_{{\rm{IC}}}({\gamma }_{{\rm{e}},\chi })+{b}_{{\rm{Synch}}.}({\gamma }_{{\rm{e}},\chi },{\boldsymbol{r}})\\ & = & {b}_{{\rm{IC}}}^{0}{u}_{{\rm{CMB}}}{E}_{{\rm{e}},\chi }^{2}+{b}_{{\rm{IC}}}^{0}{u}_{{\rm{SL}}}{E}_{{\rm{e}},\chi }^{2}\\ & & +{b}_{{\rm{Synch}}.}^{0}{B}_{{\rm{ISM}}}^{2}({\boldsymbol{r}}){E}_{{\rm{e}},\chi }^{2},\end{array}\end{eqnarray} \tag{ 8 }$$

where the coefficients have the values b_Syn.⁰ = 0.0254 and b_IC⁰ = 0.76 in units of 10⁻¹⁶ GeV⁻¹ (e.g., Colafrancesco et al. 2006; McDaniel et al. 2018), the photon energy densities are ${u}_{\mathrm{SL}}=8\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$ for starlight and ${u}_{{\rm{CMB}}}=0.25\,{\rm{eV}}\,{{\rm{cm}}}^{-3}$ for CMB photons (e.g., Porter et al. 2008; Profumo & Ullio 2010), ${E}_{{\rm{e}},\chi }={\gamma }_{{\rm{e}},\chi }{m}_{{\rm{e}}}{c}^{2}$ is the energy of the DMEs in units of GeV, and B_ISM(r) is the strength of the interstellar magnetic field, which is assumed to be a constant, $\sim 10\,\mu {\rm{G}}$.

When the external forward shock encounters the circumburst medium, the DMEs are accelerated. Here, a simple scenario of Fermi-type shock acceleration is applied in which the energy of the electrons can increase by a factor on the order of Γ during the first shock crossing and increase by a factor of $\bar{g}\sim 2.0$ in subsequent shock-crossing cycles in the frame of the shock downstream (e.g., Gallant & Achterberg 1999; Achterberg et al. 2001). It should be noted that DMEs obtain energy of $\delta ={\epsilon }_{{\rm{e}}}\eta {m}_{{\rm{p}}}(p-2)({\rm{\Gamma }}-1)/{m}_{{\rm{e}}}(p-1){\rm{\Gamma }}$ in the precursor of the shock before acceleration (e.g., Lemoine & Pelletier 2017) due to the interactions of the DMEs with the electromagnetic field (e.g., Kumar et al. 2015). After every cycle, only a fraction of the DMEs can return to the upstream region, while the residual DMEs escape far downstream. According to the results of Monte Carlo simulations of test particles, we assume that the DMEs downstream have the same probability ${\bar{P}}_{\mathrm{ret}}(\sim 0.4)$ to return to the upstream region and that the probability of returning upstream is unity (e.g., Lemoine & Pelletier 2003). Hence, the Lorentz factors and number density spectra of the escaping DMEs after the nth cycle are $\gamma {{\prime} }_{{\rm{e}},\chi }^{N}={\bar{g}}^{N}{\rm{\Gamma }}({\gamma }_{{\rm{e}},\chi }+\delta )$ and ${f}_{\mathrm{esc}}^{N}(\gamma {{\prime} }_{{\rm{e}},\chi }^{N})\,={\bar{P}}_{\mathrm{ret}}^{N}{\bar{g}}^{-N}(1-{\bar{P}}_{\mathrm{ret}}){{dn}}_{{\rm{e}},\chi }({\gamma }_{{\rm{e}},\chi })/k{\rm{\Gamma }}d{\gamma }_{{\rm{e}},\chi }$, respectively. Here, $k\,=({\gamma }_{{\rm{e}},\chi }+\delta )/{\gamma }_{{\rm{e}},\chi }$. By accumulating the number density spectra of the accelerated DMEs, ${{dn}}_{{\rm{e}},\chi }^{{\prime} }({\gamma }_{{\rm{e}},\chi }^{{\prime} })/d{\gamma }_{{\rm{e}},\chi }^{{\prime} }\,={\sum }_{N\,=\,0}^{N=\infty }{f}_{{\rm{esc}}}^{N}({\gamma }_{{\rm{e}},\chi }^{{\prime} })$ is obtained. The fact that this method can obtain the initial spectrum of shock-accelerated CMEs that is generally assumed to have a power-law form illustrates that simplifying the acceleration of the DMEs in this way is feasible.

It should be noted that when the number density of DMEs is much larger than that of CMEs, the value of η is much less than unity, resulting in ${\gamma }_{{\rm{e}},\min }^{{\prime} }\ll {\rm{\Gamma }}$, which conflicts with the increasing energy of DMEs during the first shock crossing. Thus, the shock-induced acceleration of DMEs in this case is inefficient. In other words, not all DMEs are efficiently accelerated. Considering the contributions of nonaccelerated DMEs simultaneously, ${\hat{Q}}_{{\rm{DM}}}^{{\prime} }({\gamma }_{{\rm{e}},\chi }^{{\prime} },R)$ can be derived as

$$\begin{eqnarray}\begin{array}{rcl}{\hat{Q}}_{{\rm{DM}}}^{{\prime} }({\gamma }_{{\rm{e}},\chi }^{{\prime} },R) & = & 4\pi {R}^{2}\left[A\displaystyle \frac{{{dn}}_{{\rm{e}},\chi }^{{\prime} }({\gamma }_{{\rm{e}},\chi }^{{\prime} })}{d{\gamma }_{{\rm{e}},\chi }^{{\prime} }}\right.\\ & & \left.+(1-A)\displaystyle \frac{{{dn}}_{{\rm{e}},\chi }({\gamma }_{{\rm{e}},\chi })}{{\rm{\Gamma }}d{\gamma }_{{\rm{e}},\chi }}\right],\end{array}\end{eqnarray} \tag{ 9 }$$

where A is set as 1 when $\delta \gtrsim 1;$ otherwise, the value of A equals the ratio of the number density of shock-accelerated DMEs to the number density of all DMEs. Here, $\delta \simeq 1$ is treated as the boundary in the case of efficient acceleration between all and some DMEs.

2.3. Afterglow Radiation

The synchrotron radiation power at the frequency $\nu ^{\prime}$ can be obtained by (Rybicki & Lightman 1979)

$$\begin{eqnarray}&&{P}_{{\rm{syn}}}^{{\prime} }({\nu }^{{\prime} })=\displaystyle \frac{\sqrt{3}{e}^{3}{B}^{{\prime} }}{{m}_{{\rm{e}}}{c}^{2}}{\int }_{{\gamma }_{{\rm{e}},{\min }}^{{\prime} }}^{{\gamma }_{{\rm{e}},{\max }}^{{\prime} }}\left(\displaystyle \frac{{{dN}}_{{\rm{e}}}^{{\prime} }}{d{\gamma }_{{\rm{e}}}^{{\prime} }}\right)F\left(\displaystyle \frac{{\nu }^{{\prime} }}{{\nu }_{{\rm{c}}}^{{\prime} }}\right)d{\gamma }_{{\rm{e}}}^{{\prime} },\end{eqnarray} \tag{ 10 }$$

where ${\nu }_{{\rm{c}}}^{{\prime} }=3{{eB}}^{{\prime} }{\gamma }_{{\rm{e}}}^{{\prime} 2}/4\pi {m}_{{\rm{e}}}c$ and $F({\nu }^{{\prime} }/{\nu }_{{\rm{c}}}^{{\prime} })=({\nu }^{{\prime} }/{\nu }_{{\rm{c}}}^{{\prime} }){\int }_{{\nu }^{{\prime} }/{\nu }_{{\rm{c}}}^{{\prime} }}^{+\infty }{K}_{5/3}(x){dx}$, where ${K}_{5/3}(x)$ is the modified Bessel function of order 5/3. The number density of SSC seed photons can be written as (e.g., Fan et al. 2008)

$$\begin{eqnarray}&&{n}_{{\rm{ph}}}^{{\prime} }({\nu }^{{\prime} })\approx \displaystyle \frac{1}{4\pi {R}^{2}}\displaystyle \frac{1}{{ch}{\nu }^{{\prime} }}{P}_{{\rm{syn}}}^{{\prime} }({\nu }^{{\prime} }),\end{eqnarray} \tag{ 11 }$$

and thus, the SSC power at the frequency ${\nu }_{\mathrm{ic}}^{{\prime} }$ is (e.g., Blumenthal & Gould 1970)

$$\begin{eqnarray}\begin{array}{rcl}{P}_{{\rm{SSC}}}^{{\prime} }({\nu }_{{\rm{ic}}}^{{\prime} }) & = & \displaystyle \frac{3{\sigma }_{{\rm{T}}}{ch}{\nu }_{{\rm{ic}}}^{{\prime} }}{4}{\int }_{{\nu }_{min}^{{\prime} }}^{{\nu }_{max}^{{\prime} }}\displaystyle \frac{{n}_{{\rm{ph}}}^{{\prime} }({\nu }^{{\prime} })d{\nu }^{{\prime} }}{{\nu }^{{\prime} }}\\ & & \times {\int }_{{\gamma }_{{\rm{e}},min}^{{\prime} }}^{{\gamma }_{{\rm{e}},max}^{{\prime} }}\displaystyle \frac{y(q,g)}{{\gamma }_{{\rm{e}}}^{{\prime} 2}}\displaystyle \frac{{{dN}}_{{\rm{e}}}^{{\prime} }}{d{\gamma }_{{\rm{e}}}^{{\prime} }}d{\gamma }_{{\rm{e}}}^{{\prime} },\end{array}\end{eqnarray} \tag{ 12 }$$

where $y(q,g)=2q\,\mathrm{ln}q+(1+2q)(1-q)+8{q}^{2}{g}^{2}(1-q)/(1+4{qg})$, q = w/4g(1−w), $g={\gamma }_{{\rm{e}}}^{{\prime} }h{\nu }^{{\prime} }/{m}_{{\rm{e}}}{c}^{2}$, and $w=h{\nu }_{\mathrm{ic}}^{{\prime} }/{\gamma }_{{\rm{e}}}^{{\prime} }{m}_{{\rm{e}}}{c}^{2}$. The observed spectral flux is (e.g., Granot et al. 1999)

$$\begin{eqnarray}&&{F}_{{\nu }_{\mathrm{obs}}}=\displaystyle \frac{1+z}{4\pi {D}_{{\rm{L}}}^{2}}\int {\int }_{(\mathrm{EATS})}P^{\prime} ({\nu }^{{\prime} }){D}^{3}d{\rm{\Omega }},\end{eqnarray} \tag{ 13 }$$

where "EATS" is the equal-arrival time surface corresponding to the same observer time, ${\nu }^{{\prime} }=(1+z){\nu }_{\mathrm{obs}}/D$ (with D being the Doppler factor of the emitter), ${D}_{{\rm{L}}}$ is the luminosity distance in the standard ΛCDM cosmology model (${{\rm{\Omega }}}_{M}=0.27$, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.73$, and ${H}_{0}=71\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$), and z is the redshift of the burst.

Following the above method, we can calculate the light curves and spectra of GRB afterglows, including the effects of DM. The universal parameters of the external shock model are set as the isotropic kinetic energy ${E}_{{\rm{k}},{\rm{i}}{\rm{s}}{\rm{o}}}={10}^{53}\,{\rm{e}}{\rm{r}}{\rm{g}}$, ${\epsilon }_{{\rm{e}}}={10}^{-1}$, ${\epsilon }_{B}={10}^{-3}$, ${\rm{\Gamma }}=200$, p = 2.3, ${\theta }_{{\rm{j}}}=0.1$, and z = 1.

Figure 1(a) shows the light curves of GRB afterglows including and excluding the effects of DM with $n=0.1\,{\mathrm{cm}}^{-3}$. DM annihilation occurs in the $b\bar{b}$ channel, and the DM particle mass is ${M}_{\chi }=10\,\mathrm{GeV}$. The solid and dashed lines represent the light curves and spectra with ${\rho }_{\chi }=5.0\times {10}^{4}\,\mathrm{GeV}\,{\mathrm{cm}}^{-3}$ and 0, respectively. The magenta and blue lines correspond to the synchrotron emission in the optical R band and $0.3\mbox{--}10\,\mathrm{KeV}$, respectively. The red lines correspond to the synchrotron + SSC emissions in the range of 0.1–1 GeV, and the green lines correspond to the SSC emissions in the range of 0.3–1 TeV. The fluxes in the ranges of 0.3–10 keV and 0.3–1 TeV are multiplied by 10 and 100, respectively. In this case, the annihilation DM effects are adequately displayed in the multiband light curves of GRB afterglows. Figure 1(b) shows the synchrotron + SSC spectra in the same situation as Figure 1(a). The green, red, and blue lines represent the spectra at ${10}^{3}\,{\rm{s}}$, ${10}^{4}\,{\rm{s}}$, and ${10}^{5}\,{\rm{s}}$, respectively. Once the effects of DM are considered, in addition to fluence increasing, the shapes of the spectra also change.

Zoom In Zoom Out Reset image size

Figure 2 shows the synchrotron + SSC spectra of GRB afterglows with different parameters, i.e., DM density, particle DM mass, annihilation channel, and CME density. Figure 2(a) displays the synchrotron + SSC spectra of GRB afterglows with different n at $t={10}^{4}\,{\rm{s}}$ after the GRB is triggered. The DM density is set to $5\times {10}^{4}\,\mathrm{GeV}\,{\mathrm{cm}}^{-3}$ or $0\,\mathrm{GeV}\,{\mathrm{cm}}^{-3}$, and the DM particle mass is taken as 10 GeV. The red, blue, green, and magenta lines represent $n=1.0\,{\mathrm{cm}}^{-3}$, ${10}^{-1}\,{\mathrm{cm}}^{-3}$, ${10}^{-2}\,{\mathrm{cm}}^{-3}$, and ${10}^{-3}\,{\mathrm{cm}}^{-3}$, respectively. The solid and dashed lines represent the cases including and excluding the contribution of the DMEs, respectively. Different CME densities result in different spectral shapes, and the inclusion of DMEs further makes the variety of the spectra. For a lower CME density, the DME effect is more significant.

Zoom In Zoom Out Reset image size

The spectra in four different annihilation channels are presented in Figure 2(b). It is obvious that the fluxes in the ${w}^{+}{w}^{-}$ and $b\bar{b}$ annihilation channels are higher than those in the ${\tau }^{+}{\tau }^{-}$ and ${\mu }^{+}{\mu }^{-}$ annihilation channels.

Figure 2(c) shows the synchrotron + SSC spectra of GRB afterglows for different DM densities and DM masses at $t={10}^{4}\,{\rm{s}}$. The values of the DM density are set to $5\times {10}^{4}\,\mathrm{GeV}\,{\mathrm{cm}}^{-3}$ and $5\times {10}^{5}\,\mathrm{GeV}\,{\mathrm{cm}}^{-3}$, and the values of DM particle mass are taken as 10 GeV and 100 GeV. The DM annihilation channel is assumed to be the $b\bar{b}$ channel, and the CME number density is $n=1\,{\mathrm{cm}}^{-3}$. Figure 2(d) presents the afterglow spectra for the same parameters of Figure 2(c), except the parameters $n={10}^{-3}\,{\mathrm{cm}}^{-3}$, and ${\rho }_{\chi }$ = $5\,\times {10}^{3}\,\mathrm{GeV}\,{\mathrm{cm}}^{-3}$ and $5\times {10}^{4}\,\mathrm{GeV}\,{\mathrm{cm}}^{-3}$. The black solid lines in Figures 2(c) and (d) denote the spectra excluding the contributions of the DMEs. The flux values of the black lines are slightly higher than those of the red solid lines (corresponding to lower DM densities and high DM particle masses) because ${\gamma }_{{\rm{e}},\min }$ decreases slightly when the effects of DM are considered. According to all the lines, a higher DM density and lower DM particle mass can dramatically increase the fluxes. As shown by the dashed blue line in Figure 2(d), the flux exhibits prominent improvements in the low- and high-energy parts due to the contributions of inefficiently accelerated DMEs and the increased number of seed photons, respectively.

Since many parameters related to jets are uncertain and the ISM and DM can prominently influence the light curves and spectra of GRB afterglows, it is difficult to distinguish the distributions of DM unless the properties of the circumburst medium can be constrained by observations. For a full description of the medium properties of GRB afterglows, the location of GRBs in the galaxy should be considered, as discussed below.

To simplify, we assume a DM halo with the virial mass ${M}_{\mathrm{vir}}$ in the virial radius ${R}_{\mathrm{vir}}$. The mean density is equal to the virial overdensity ${{\rm{\Delta }}}_{\mathrm{vir}}$ times the mean background density ${\rho }_{{\rm{u}}}$, given by (e.g., Bullock et al. 2001)

$$\begin{eqnarray}&&{M}_{\mathrm{vir}}\equiv \displaystyle \frac{4\pi }{3}{{\rm{\Delta }}}_{\mathrm{vir}}{\rho }_{{\rm{u}}}{R}_{\mathrm{vir}}^{3}.\end{eqnarray} \tag{ 14 }$$

${{\rm{\Delta }}}_{\mathrm{vir}}=200$ is generally preferred, which is independent of cosmology. Thus, the corresponding virial radius and virial mass are ${R}_{200}$ and ${M}_{200}$, respectively. Furthermore, the concentration parameter (e.g., Bullock et al. 2001) can be defined as

$$\begin{eqnarray}&&{c}_{200}\equiv \displaystyle \frac{{R}_{200}}{{r}_{{\rm{c}}}},\end{eqnarray} \tag{ 15 }$$

which is also related to the virial mass (e.g., Buote et al. 2007), i.e.,

$$\begin{eqnarray}&&{c}_{200}=\displaystyle \frac{6.9}{(1+z)}{\left(\displaystyle \frac{{M}_{200}}{{10}^{14}{M}_{\odot }}\right)}^{-0.178}.\end{eqnarray} \tag{ 16 }$$

Moreover, the virial mass can be derived by integrating different DM density profiles in the viral radius R₂₀₀,

$$\begin{eqnarray}&&{M}_{200}=4\pi {\int }_{0}^{{R}_{200}}\rho (r){r}^{2}\,{dr}.\end{eqnarray} \tag{ 17 }$$

The DM density of the galaxies is a function of radius. The general Navarro–Frenk–White (NFW) density profile reads (Navarro et al. 1997)

$$\begin{eqnarray}&&{\rho }_{\mathrm{NFW}}(r)=\displaystyle \frac{{\rho }_{{\rm{c}}}}{{\left(r/{r}_{{\rm{c}}}\right)}^{{\gamma }_{{\rm{c}}}}{\left[1+(r/{r}_{{\rm{c}}})\right]}^{3-{\gamma }_{{\rm{c}}}}},\end{eqnarray} \tag{ 18 }$$

where ${\rho }_{{\rm{c}}}$ and ${r}_{{\rm{c}}}$ are the characteristic density and radius, respectively. Here, we consider four values of the slope: ${\gamma }_{{\rm{c}}}=1.0$, 1.1, 1.2, and 1.3, respectively. The profile proposed by Moore et al. (1999) corresponding to the slope ${\gamma }_{{\rm{c}}}=1.5$ is derived according to N-body simulations, which suggests a steeply rising DM density toward the galaxy center. In addition, the Einasto law is written as (Merritt et al. 2006)

$$\begin{eqnarray}&&{\rho }_{\mathrm{Einasto}}(r)={\rho }_{{\rm{c}}}\exp [-{d}_{N}({\left(r/{r}_{{\rm{c}}}\right)}^{1/N}-1)],\end{eqnarray} \tag{ 19 }$$

where ${d}_{N}\approx 3N-1/3+0.0079/N$ for $N\gtrsim 0.5$.

Combining the above equations, one can obtain the characteristic radius and density of the DM density profiles with the virial mass ${M}_{200}={10}^{12}\,{M}_{\odot }$. The results are shown in Figure 3. The blue and green solid lines represent the Einasto law and Moore law, respectively. The red lines represent the NFW formula, where the solid, dashed, dotted, and dashed–dotted lines denote the NFW profile with the following four values of the inner slope: 1.0, 1.1, 1.2, and 1.3, respectively. The gray shaded region indicates the DM density range of $\sim {10}^{3}-{10}^{5}\,\mathrm{GeV}\,{\mathrm{cm}}^{-3}$. The vertical black dashed line denotes the location 1 pc from the center of the galaxy. According to the magnetic field model used in the galaxy center (Aloisio et al. 2004), ${B}_{\mathrm{ISM}}$ increases rapidly with decreasing radius within 1 pc, so the value of ${B}_{\mathrm{ISM}}$ assumed above would be inappropriate within 1 pc. Finally, the brown shaded region represents the possible locations of GRBs with significant contributions from DM, i.e., several to tens of parsecs from the center of a galaxy with virial mass $\sim {10}^{12}\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

Moreover, the effects of DM spatial diffusion (e.g., Colafrancesco et al. 2006), reacceleration, and advection (e.g., Strong et al. 2007) might be important at the galaxy scale or beyond, which would decrease the DME densities around GRBs. These effects are neglected for our calculations.

In this paper, we study the contributions of DM annihilation to GRB afterglows. If GRBs occur at distances of several to tens of parsecs from the centers of massive galaxies, the effects of DMEs should be significant. Moreover, the influences of different DM particle masses and annihilation channels on GRB afterglows are adequately reflected in the flux changes.

To date, no GRB has been discovered at a distance of tens of parsecs from a galaxy center, but one may expect that these events could be detected in the future because the density of stars near the galaxy center is generally higher than that in other regions within the host galaxy. Observations show that the majority of stars formed at distances of a few parsecs from the center of the Milky Way are older than 5 Gyr (e.g., Blum et al. 2003), and one pulsar, PSR ${\rm{J}}1745-2900$, was detected ∼0.12 pc from Sagittarius ${{\rm{A}}}^{\star }$ (Eatough et al. 2013). Of course, the dense gas and bright galaxy center will trouble the observations of those distant GRBs.

In addition to DMEs, gamma-ray emission is another final product of DM annihilation due to ${\pi }^{0}$ decay. The gamma-ray flux produced by ${\pi }^{0}$ decay from a galaxy with ${M}_{\mathrm{vir}}={10}^{12}\,{M}_{\odot }$ at z = 1 is far below the observational threshold and hence can be neglected. However, if this galaxy is located near the Milky Way, the gamma-ray flux can be detected. Therefore, we use observations of a neighboring galaxy, M31, which has a virial mass of $\sim {10}^{12}\,{M}_{\odot }$, to test our model under the assumption that the gamma-ray emission near the center of M31 originated from DM annihilation. In the observation, a gamma-ray excess with a flux of $(5.6\pm 0.6)\times {10}^{-12}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ was detected in the energy range from 0.1 to $100\,\mathrm{GeV}$ within the central M31 region of $\sim 5\,\mathrm{kpc}$ (Ackermann et al. 2017). We can calculate the gamma-ray flux by the formula

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d{{\rm{\Phi }}}_{\gamma }({E}_{\gamma })}{{{dE}}_{\gamma }} & = & \displaystyle \frac{1}{{D}_{A}^{2}}\displaystyle \frac{\langle {\sigma }_{A}\upsilon \rangle }{2{m}_{\chi }^{2}}\displaystyle \frac{{{dN}}_{\gamma }({E}_{\gamma })}{{{dE}}_{\gamma }}{E}_{\gamma }^{2}\\ & & \times {\int }_{{r}_{\min }}^{{r}_{\max }}{\rho }_{\chi }^{2}(r){r}^{2}{dr},\end{array}\end{eqnarray} \tag{ 20 }$$

where ${{dN}}_{\gamma }({E}_{\gamma })/{{dE}}_{\gamma }$ is the gamma-ray spectrum corresponding to ${M}_{\chi }=100\,\mathrm{GeV}$ and the $b\bar{b}$ annihilation channel and ${r}_{\max }\,=5\,\mathrm{kpc}$, and D_A is the distance of M31. The NFW profile with ${\gamma }_{{\rm{c}}}=1.3$ is adopted. By integrating the energy range from 0.1 to $100\,\mathrm{GeV}$, the gamma-ray flux is $5.7\,\times {10}^{-13}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$. This flux is less than that observed in M31, indicating that our model is self-consistent.

This work was supported by the National Natural Science Foundation of China under grants 11822304, 11890692, 11773007, and U1531130, and the Guangxi Science Foundation under grant 2018GXNSFFA281010.