# diffSph: a Python tool to compute diffuse signals from dwarf
spheroidal galaxies

Martin Vollmann, Finn Welzmüller, Lovorka Gajović  0000-0003-2193-8922
0009-0009-2456-8884 0000-0001-9512-8565

###### Abstract

So far no diffuse emissions in dwarf spheroidal satellites of the Milky Way
have ever been observed. Given that dwarf galaxies are predominantly composed
of Dark Matter, the discovery of these signals could offer valuable insights
into understanding the nature of Dark Matter. We present “diffSph”, a Python
tool which in its present version provides fast predictions of such diffuse
signals in radio frequencies. It also features a very comprehensive module for
the computation of “J” and “D” factors that are relevant for indirect Dark
Matter detection using gamma rays. Routines are coupled to parton-shower
algorithms and Dark Matter halo mass functions from state-of-the-art kinematic
fits. This code is also useful for testing generic hypotheses (not necessarily
associated with any Dark Matter candidate) about the cosmic-ray
electron/positron sources in the dwarf galaxies. The diffSph tool has already
been employed in searches for diffuse signals from dwarf spheroidal galaxies
using the LOw Frequency ARray (LOFAR).

##  1 Introduction

Due to their small star formation rates, Milky Way satellite dwarf spheroidal
galaxies (dSph) are weak emitters, and to this date, no diffuse emission from
these objects has been observed at any frequency. Most of their scarce
luminosity is produced by their member stars instead. Dwarf galaxies are,
however, Dark-Matter (DM) rich environments and in several theories of
particle DM such diffuse emissions are predicted via annihilation or decay
[1]. While the signal strengths are typically weak, they are typically within
the detectable range of current telescopes. Several collaborations have
already conducted dedicated searches for these signals using e. g. gamma-ray
telescopes [2] and found no compelling evidence of any diffuse signal in the
observed targets. Instead, their non-detection results have been used in order
to set exclusion limits on the parameter space of particle DM models.

Gamma rays are the “golden” messenger for indirect DM detection better-
controlled uncertainties, theoretically motivated models (e. g. WIMP),
phenomenologically promising signatures and ever increasing improvements in
sensitivity. Nevertheless, complementary indirect searches using radio
interferometers such as the LOw Frequency ARray (LOFAR) can yield to even more
stringent exclusion bounds on particle DM models [3, 4]. The latter limits
are, however, subject to much larger uncertainties than those encountered in
the gamma-ray counterpart.

Independent of these intrinsic limitations, if diffuse synchrotron-radiation
emission in dSphs is discovered soon, solid theoretical interpretations are
necessary. Therefore, we introduce diffSph, a Python tool that produces signal
templates for diffuse synchrotron emissions within the context of diffusive
cosmic-ray propagation powered by magnetic-field turbulence and for _generic_
hypotheses on the origin of these signals. Special focus is put on the
hypothesis that the cosmic-ray electrons are produced via DM annihilation or
decay but the modularity of diffSph allows for several other hypotheses (not
associated to DM) to be tested.

More specifically, the diffSph package includes several pre-computed objects
that speed up the computations that are relevant for testing several
interesting/popular particle DM scenarios. For instance, diffSph can
efficiently generate synchrotron spectra from decaying or WIMP DM with
particle masses in the 1 GeV - 100 TeV range by evaluating pre-computed kernel
functions. This renders its performance quite efficient without losing much
accuracy.

The diffSph code has also a module (diffsph.profiles) that computes the line-
of-sight and angular integrals that are relevant for indirect DM detection
phenomenology in general. This is linked to several mass models of dwarf
galaxies found in the literature. In particular, this module uses these mass
models as input in its routines for the computation of the so-called
astrophysical $J$\- and $D$-factors, which are essential pieces for indirect
DM detection with gamma rays. The code can be easily accessed and downloaded
via GitHub, e. g.

git clone https://github.com/mertio1/diffsph.git

If the user wants to install diffSph globally, they can do it by going into
diffSph’s main folder and on a terminal typing (the installation can take a
few minutes)

python setup.py bdist_wheel

pip install .

An extensive online documentation is available at
https://mertio1.github.io/diffsph/ and, as a PDF file, at this URL. The main
novelty of diffSph with respect to already existing tools (e. g. [5, 6]), is
that instead of using a method-of-images expansion for the cosmic-ray electron
trasport equation, a Fourier decomposition is employed which turns out to be
more suitable for the computation at hand.

The article is organized as follows: in section 2 we provide a review about
the physical modelling of the diffuse emission in dwarf galaxies. In the next
part (section 3) we describe how diffSph is structured and explain how the
different hypotheses that are incorporated in the code can be tested using it.
In section 4 we show some concrete examples and we then conclude. In the
appendices we include a brief discussion on some asymptotic approximations and
analytical results that are adopted by diffSph.

##  2 Theoretical background

The description of the transport of cosmic-ray electrons (CREs) in dSphs is a
puzzling problem. This is because _no_ magnetic fields in these galaxies have
been observed and the resolution of existing cosmological magneto-hydrodynamic
simulations is not sufficient in order to infer their distribution and
strength in these galaxies. In this work we adopt the arguably simplest
transport model for the CRE, which has become the standard lore in the
literature, see Ref. [7] and references therein. In this framework, the CRE
number density $n_{e}(r,E)$ as a function of the galacto-centric radius $r$
and the CREs’ energy $E$, is obtained by solving a stationary and spherically-
symmetric diffusion-loss equation

$D(E)\nabla^{2}n_{e}+\frac{\partial}{\partial E}[b(E)n_{e}]+q_{e}=0\ ,$ (1)

and a boundary condition $\left.n_{e}(r,E)\right|_{r=r_{h}}=0$, where $r_{h}$
is the _visible halo radius_ parameter. The model is an adaptation of the
Milky Way’s leaky-box CR transport model to the geometry of the dwarf
spheroidal galaxies. In particular, the diffusion and energy-loss coefficients
are respectively given by $D(E)=D_{0}(E/E_{0})^{\delta}$ and
$b(E)=b_{0}\,[1+(B/B_{c})^{2}]E^{2}$, where $\delta$, $D_{0}/E_{0}^{\delta}$,
$b_{0}$, $B_{c}$ are constants111We use $e=4.80326\times 10^{-10}$ Fr,
$m_{e}=9.1093837\times 10^{-28}$ g and $T_{\textrm{CMB}}=2.725$ K.  and $B$ is
an _effective parameter_ that captures the effect of the (randomly oriented)
magnetic field in the galaxy. It is assumed that $B$ is constant within the
$r<r_{h}$ volume.

The energy-loss coefficient receives contributions from two processes:
synchrotron and inverse Compton scattering (ICS) radiation. The former is
proportional to the squared magnetic field $B$ and for the latter we neglect
all radiation fields but the cosmic microwave background (CMB):

$\displaystyle b_{0}$ $\displaystyle=$
$\displaystyle\frac{32\pi^{3}e^{4}k_{B}^{4}T_{\textrm{CMB}}^{4}}{135\hbar^{3}m_{e}^{4}c^{10}}\simeq
2.652\times 10^{-17}\,\textrm{GeV}^{-1}\textrm{s}^{-1}\
,\quad(\texttt{utils.consts.b0})\ ,$ (2) $\displaystyle B_{c}$
$\displaystyle=$
$\displaystyle\sqrt{\frac{8\pi^{3}k_{B}^{4}T_{\textrm{CMB}}^{4}}{15\hbar^{3}c^{3}}}\simeq
3.238\,\mu\textrm{G}\ ,\quad(\texttt{utils.consts.Bc})\ ,$ (3)

the bracketed variable names in the equations above indicate their
corresponding nomenclature in diffSph.

The source term $q_{e}(r,E)$ in Eq. (1) can be factored as

$q_{e}(r,E)=\mathcal{Q}\,\mathcal{R}(r)\,\mathcal{S}(E)\ ,$ (4)

where $\mathcal{R}(r)$, $\mathcal{S}(E)$ are hypothesis-dependent functions;
and $\mathcal{Q}$ is a normalization constant. We consider generic hypotheses,
e. g. WIMP, decaying DM, point source, etc. (see Table 1). In the WIMP
hypothesis case, for example, the normalization constant is given by222If the
variable self_conjugate is set True $\mathcal{Q}=\langle\sigma
v\rangle/(2m_{\rm DM}^{2})$, where $\langle\sigma v\rangle$ is the total
velocity-averaged annihilation cross section times the relative WIMPs’
velocity and $m_{\rm DM}$ is their mass; $\mathcal{R}(r)=\rho^{2}(r)$, where
$\rho(r)$ is the DM halo density and $\mathcal{S}(E)\equiv(1/\langle\sigma
v\rangle)\times\operatorname{d\\!}\,\langle\sigma
v\rangle/\operatorname{d\\!}E$ is the WIMP-model specific differential
electron/positron “yield” per annihilation.

Hypothesis | $\mathcal{Q}$ | $\mathcal{R}(r)$ | $\mathcal{S}(E)$
---|---|---|---
’wimp’ | $\langle\sigma v\rangle/(2m_{\rm DM}^{2})$ | $\rho^{2}(r)$ | $(1/\langle\sigma v\rangle)\times\operatorname{d\\!}\,\langle\sigma v\rangle/\operatorname{d\\!}E$
’decay’ | $\Gamma_{\chi}/m_{\rm DM}$ | $\rho(r)$ | $(1/\Gamma_{\chi})\times\operatorname{d\\!}\,\Gamma_{\chi}/\operatorname{d\\!}E$
’generic’ | $\mathcal{N}$ | $\mathcal{R}_{\star}(r)$ | $E_{0}^{\Gamma-1}/E^{\Gamma}$
Table 1:  Possible hypotheses of diffSph and their associated source-term
functions and parameters.

###  2.1 Fourier-expanded solution

Independent of the hypothesis, Eq. (1) can be solved semi-analytically as a
linear combination of (Fourier) basis functions so that e. g. the brightness
$I_{\nu}(\theta)$ as a function of the observed angular radius $\theta$ can be
expressed as

$I_{\nu}(\theta)=\frac{\mathcal{Q}}{4\pi}\sum_{n=1}^{\infty}h_{n}\times
X_{n}(\nu)\times f_{n}(\theta/\theta_{h})$ (5)

where $R$ is the dSph’s distance and the variable
$\theta_{h}=\arcsin(r_{h}/R)\approx r_{h}/R$ is the angular radius of the
diffuse signal. The basis functions $f_{n}(x)$, which are defined as

$f_{n}(x)=2\,\int_{x}^{1}\frac{\sin(n\pi y)\,{\rm d}y}{\sqrt{y^{2}-x^{2}}}\ ,$
(6)

do not depend on either the CRE sources nor their propagation: they are
universal [3].

The remaining $\mathcal{Q}\times h_{n}\times X_{n}(\nu)$ terms in the formula
carry all the dependence on both the transport and CRE source models.
Concretely, the $n$-th order _halo/bulge factors_ $h_{n}$ are given by

$h_{n}=\frac{2}{r_{h}}\int_{0}^{r_{h}}{\rm d}rr\sin\frac{n\pi
r}{r_{h}}\mathcal{R}(r)\ .$ (7)

and the _spectral coefficients_ $X_{n}(\nu)$ can be written in terms of the
spectral function $X(\nu)$ as $X_{n}(\nu)\equiv X(\nu;\tau_{0}/n^{2})$. The
latter is given in terms of the diffusion time scale parameter
$\tau_{0}=r_{h}^{2}/D_{0}$ by

$X(\nu;\tau_{0})=\frac{2\sqrt{3}e^{3}B}{m_{e}c^{2}}\int_{0}^{\infty}\frac{{\rm
d}z}{z}\mathcal{F}(z)\frac{E_{c}(\nu/z)}{b(E_{c}(\nu/z))}{\rm
e}^{-\eta(E_{c}(\nu/z);\tau_{0})}\int_{E_{c}(\nu/z)}^{\infty}{\rm
d}E\,\mathcal{S}(E){\rm e}^{+\eta(E;\tau_{0})}\ ,$ (8)

where $E_{c}(\nu)=\sqrt{2\pi
m_{e}^{3}c^{5}\nu/(3eB)}=\epsilon_{0}\sqrt{\nu/B}$ ($\epsilon_{0}=5.576$ GeV
$\mu$G1/2 GHz-1/2, utils.consts.eps0) can be understood as the typical energy
of a CRE emitting synchrotron radiation at frequency $\nu$ for a magnetic
field $B$.

The diffSph package includes several built-in templates for the radial
functions $\rho(r)$ and $\mathcal{R}_{\star}(r)$ and for the electron/positron
yields $\mathcal{S}(E)$ (see Tab. 1). While the radial templates
$\mathcal{R}(r)$ are given analytically in Tables 2-3, the $\mathcal{S}(E)$
functions have to be computed numerically in the case of DM annihilation or
decay. In the standard approach, $\mathcal{S}(E)$ is given by

$\mathcal{S}(E)=\sum_{I}\textrm{BR}_{I}\frac{\operatorname{d\\!}N_{I}}{\operatorname{d\\!}E}\
,$ (9)

where BR stands for _branching ratio_ , i. e. BR${}_{I}=\langle\sigma
v\rangle_{I}/\langle\sigma v\rangle$, ($\sum_{I}$BR${}_{I}=1$), and the index
$I$ runs over primary final-state particle/antiparticle pairs (referred to as
_channels_ in diffSph). These depend on the mass of the DM particle and are
readily available in specialized software tools like DarkSUSY [8], MicrOMEGAS
[9] and PPPC [10]. In this work, we use the (electroweak corrected) results by
PPPC for the following possible primary final states:

$\displaystyle\texttt{channel}=\\{e^{+}e^{-}\ (\texttt{`ee'}\textrm{ or
}\texttt{`ee\\_MC'}),\mu^{+}\mu^{-}\ (\texttt{`mumu'}),{}$
$\displaystyle{}\tau^{+}\tau^{-}\ (\texttt{`tautau'}),\nu\bar{\nu}\
(\texttt{`nunu'}),$ $\displaystyle q\bar{q}\ (\texttt{`qq'}),c\bar{c}\
(\texttt{`cc'}),b\bar{b}\ (\texttt{`bb'}),t\bar{t}\ (\texttt{`tt'}),{}$
$\displaystyle{}W^{+}W^{-}\ (\texttt{`WW'}),ZZ\ (\texttt{`ZZ'}),hh\
(\texttt{`hh'})\\}\ ,$ (10)

where $q$ refers to any neutral combination of the light quarks $u$, $d$, $s$
and their antiparticles. In the $e^{+}e^{-}$ case, there are two
possibilities. The channel = ‘ee’ scenario corresponds to the
$\mathcal{S}(E)=\delta(E-m_{\rm DM})$ hypothesis, whereas channel = ‘ee_MC’
incorporates soft/collinear radiation effects which are obtained numerically
using the aforementioned software packages.

##  3 Architecture of diffSph

Fig. 1 shows what the main capabilities of diffSph 1.0. are and how these are
interconnected. The chief purpose of the code is obtaining fast predictions
for synchrotron fluxes in dwarf galaxies once the hypothesis is given. This
can be done using the pyflux module. Users can also exclude hypotheses given
some experimental data using the limits module.

[width=.8]figs/mindmap_diffsph.pdf

Figure 1:  Code architecture of diffSph. Main modules are represented as white
blobs with blue borders and their extensions (red blobs) show the associated
top-level functions. Sub-modules are shown in blue while the input parameters
in gray. Rectangular nodes correspond to outsourced inputs and their use is
optional.

###  3.1 Predicted fluxes with pyflux

The pyflux module calls python functions that are contained in the profiles
and spectra submodules in order to obtain diffuse emission fluxes for a given
object and a set of parameters. For example, if the user wants to obtain the
brightness $I_{\nu}(\theta)$ of a given galaxy galaxy, at the frequency nu (in
GHz) and as a function of angle theta (in arcmin) in the WIMP hypothesis, the
pyflux module will compute and add all necessary series factors in Eq. (5) for
that particular galaxy. The user has to specify which radial template rad_temp
(Einasto, NFW, etc.) and the reference from which the relevant parameters
($\rho_{s}$, $r_{s}$, etc.) should be extracted. The current version includes
mass models from Refs. [11, 12, 13, 14, 15] (see Table 4). Half-light radii
and distances are obtained from [16, 17, 18].

The default transport-model parameter set of diffSph is inspired by CR
propagation fits from the Milky Way such as [19]: $B=2\,\mu$G,
$D_{0}=2.77\times 10^{28}$ cm2/s and 333However, note that other authors [20,
21, 22, 23] obtain $\delta\simeq 0.5$ consistent with Kraichnan diffusion.
$\delta=0.35$. Note that these parameters are chosen for the sake of
concreteness as no measurement of any of these exist and that the predictions
that are obtained by using them are fairly optimistic. As a concrete example,
the diffSph prediction for the brightness (in Jy/sr) of a 50 GeV WIMP
annihilating into $\mu^{+}\mu^{-}$ (i. e. best-fit model for the Galactic-
Center-Excess proposed in Ref. [24]) at 150 MHz that is emitted from Segue 1
at $\theta=0$ and adopting a Herquist/Diemand/Zhao (HDZ) mass model with
parameters fitted by Ref. [25], is obtained by typing the following command:

In[ ]: |  from diffsph import pyflux as pf
---|---
| pf.synch_brightness(
|  theta = 0., nu = .15, galaxy = ’Segue I’, rad_temp = ’HDZ’, ref =
’1408.0002’,
|  hyp = ’wimp’, sv = 3.e-26, mchi = 50., channel = ’mumu’)
Out[ ]: | 1.1055442452868214

####  3.1.1 Determination of $N$

The solution given in Eq. (5) is an infinite series on Fourier-like modes
$n=0,1,2,\ldots$ However, diffSph computes this up to a maximum value $n=N$.
In order to find a prescription for the determination of $N$ that is
reasonably large, there are two aspects that need to be considered.

First, as usual in Fourier any decomposition, the larger the mode index $n$
the better the spatial resolution of the solution. As a rule of thumb, $N\sim
1/\delta_{\textrm{res}}$ where $\delta_{\textrm{res}}\sim\Delta
r/r\simeq\Delta\theta/\theta$ is the targeted angular (radial) resolution of
the theoretical prediction (5). On the hand, regardless of the $r$\- or
$\theta$-resolution arguments, a much stronger condition that the solutions
have to fulfill is that these should convergence on a global level. More
specifically, $N$ has to be large enough that the total flux density as
integrated over the complete galaxy leads to a convergent result. This is
given by

$S_{\nu}=\frac{1}{R^{2}}\frac{r_{h}^{2}\,\mathcal{Q}}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}h_{n}\times
X_{n}(\nu)+\mathcal{O}\left(r_{h}^{4}/R^{4}\right)\ .$ (11)

We thus introduce two variables, a Boolean one high_res, which is set False by
default and a floating one accuracy (in %). Given these two, diffSph
determines $N$ by evaluating the ratio $|s_{m+1}/S_{m}|$ for ever increasing
values of $m$, where $s_{m}$ is defined as

$s_{m}=h_{m}X_{m}\ ,$ (12)

and, depending on the high_res variable, $S_{m}$ is given by

$\displaystyle
S_{m}=\left\\{\begin{array}[]{c}\displaystyle\sum_{n=1}^{m}\frac{(-1)^{n-1}}{n}s_{n}\quad,\quad\texttt{high\\_res
= True}\\\
\displaystyle\sum_{n=1}^{m}(-1)^{n-1}s_{n}\quad,\quad\texttt{high\\_res =
False}\end{array}\right.\ .$ (15)

The algorithm (pyflux.which_N) then stops when the ratio is smaller than the
chosen value (default value is 1 %) of the accuracy variable. The value of $m$
at which the algorithm stops will be the chosen $N$. A key feature of diffSph
is that while this algorithm is run, the code stores the values of all series
coefficients $\\{s_{m}\\}_{m=1}^{N}$ in the cache folder, so that the
computations are not carried out all over again when the same set of model
parameters are inputted.

###  3.2 The limits module

Besides generating flux predictions via the pyflux module, diffSph offers the
possibility of producing 2$\sigma$ exclusion limits (or limit forecasts) on
the relevant signal-strength parameters of any hypothesis (e. g.
$\langle\sigma v\rangle$ in the case of the WIMPs) given a non-detection. This
is done by the limits module. There are two possibilities which we illustrate
(in the context of the WIMP hypothesis) below.

####  3.2.1 sigmav_gausslim

This is the strategy that was employed in Refs. [4, 26]. In this approach, we
use the fact that the leading-mode shape function $f_{n=1}(\theta/\theta_{h})$
defined in (6), resembles a Gaussian function

$I_{\nu}=a_{\nu}\operatorname{e\\!}^{-\frac{\theta^{2}}{2\sigma_{h}^{2}}}\ ,$
(16)

where444This expression is obtained by minimizing
$\mathcal{F}(\sigma_{h})=\int_{0}^{\theta_{h}}\operatorname{d\\!}\theta\theta[\frac{2\theta_{h}^{2}}{\pi\sigma_{h}^{2}}\operatorname{e\\!}^{-\theta^{2}/(2\sigma_{h}^{2})}-f_{0}\left(\frac{\theta}{\theta_{h}}\right)]^{2}$
$\sigma_{h}\equiv 2\,\theta_{h}/5$. For a given $\sigma_{h}$, this Gaussian
template can be used as a fake source or, in the context of a likelihood-ratio
test, as the tested model. The non-observation of diffuse signals translates
into a limit on the amplitude $a_{\nu}$ (in Jy/beam) which diffSph can
reinterpret as exclusion limits on the annihilation cross section of WIMPs.

####  3.2.2 Forecasts with sigmav_limest

Another unique capability of diffSph is that of producing estimated upper
limits on the relevant signal-strength parameters ($\langle\sigma v\rangle$,
$\Gamma_{\chi}$, etc.) of the diffuse emission given the root-mean-square
noise of the observation. This is computed using a method devised in Ref.
[26], where it is shown that

$\mathcal{Q}_{\textrm{lim}}\approx
20.6\frac{I_{\nu}^{\textrm{rms}}\Omega_{\textrm{beam}}}{\sum_{n=1}^{N}(-1)^{n-1}h_{n}X_{n}(\nu)}\frac{2\,r_{h}}{\textrm{HFD}}\sqrt{\frac{\Omega(\theta_{h})}{\Omega_{\textrm{beam}}}}\
,$ (17)

where $\mathcal{Q}_{\text{lim}}$ is the estimated limit for
$\mathcal{Q}=\langle\sigma v\rangle/(2m_{\rm DM})$, $\Gamma_{\chi}/m_{\rm
DM}$, etc; HFD stands for _Half-Flux Diameter_ and
$\Omega(\theta_{h})=\pi\theta_{h}^{2}$.

###  3.3 Submodules

As depicted in Fig. 1, the top-level functions that are defined in the main
pyflux and limits modules are constructed using the profiles, spectra and
utils sub-modules. We briefly describe these below

####  3.3.1 Halo/bulge profiles module

The profiles module includes halo/bulge factor ($h_{n}$) calculators (see Eq.
(7)) for generic $\mathcal{R}(r)$ functions. We provide several benchmark
templates (see Tables 2-3) for these functions as well. These are
characterized by a set of parameters, e. g. $\rho_{s}$, $r_{s}$, etc., that
can be provided by the user. For example, the following command outputs the
$h_{n=3}$ factor in GeV2/cm5 for an Einasto profile with the given parameters
and assuming DM annihilation.

| Navarro-Frenk-White | Isothermal | Hernquist-Diemand-Zhao
---|---|---|---
rad_temp | NFW | SIS | HDZ
$\rho(r)$ | $\displaystyle\frac{\rho_{s}r_{s}}{r(1+r/r_{s})^{2}}$ | $\displaystyle\frac{\sigma_{v}^{2}}{2\pi Gr^{2}}$ | $\displaystyle\frac{\rho_{s}r_{s}^{\gamma}}{r^{\gamma}(1+r^{\alpha}/r_{s}^{\alpha})^{(\beta-\gamma)/\alpha}}$

| Einasto | Pseudo Isothermal | Burkert | cored NFW
---|---|---|---|---
rad_temp | Enst | ps_ISO | Bkrt | cNFW
$\rho(r)$ | $\displaystyle\rho_{s}\operatorname{e\\!}{}^{-\frac{2}{\alpha}\left(\frac{r^{\alpha}}{r_{s}^{\alpha}}-1\right)}$ | $\displaystyle\frac{\rho_{s}}{1+r^{2}/r_{s}^{2}}$ | $\displaystyle\frac{\rho_{s}}{(1+r/r_{s})(1+r^{2}/r_{s}^{2})}$ | $\displaystyle\frac{\rho_{s}}{r_{c}/r_{s}+r/r_{s}(1+r/r_{s})^{2}}$

Table 2:  Dark Matter profile templates available in diffSph. | Plummer | Point Source | Constant
---|---|---|---
rad_temp | Plmm | PS | Const
$\mathcal{R}_{\star}(r)$ | $\displaystyle\frac{3}{4\pi r_{\star}^{3}}\frac{1}{(1+r^{2}/r_{\star}^{2})^{5/2}}$ |  $\displaystyle\frac{1}{4\pi r^{2}}\delta(r)$ or $\displaystyle\frac{1}{(2\pi r_{\star}^{2})^{3/2}}\operatorname{e\\!}^{-\frac{r^{2}}{2r_{\star}^{2}}}$ | $\frac{3}{4\pi r_{\star}^{3}}$
Table 3:  Radial function templates for generic hypotheses in diffSph.

In[ ]: |  from diffsph.profiles import hfactors
---|---
| hfactors.halo_factor(
|  n = 3, rh = 0.8, rad_temp = ’enst’, rs = 0.5, rhos = 2.4, hyp = ’wimp’)
Out[ ]: | 3.263539652523281e+23

Yet another unique feature of diffSph is the outsourcing of the halo
parameters, which are obtained from state-of-the-art catalogues from stellar
kinematics (see Table 4). These can be accessed (default option) if the
Boolean variable manual is set True. For instance,

In[ ]: |  from diffsph.profiles import massmodels
---|---
| massmodels.h(
|  n = 3, galaxy = ’Ursa Minor’, ratio = 0.8, rad_temp = ’NFW’,
|  hyp = ’wimp’, manual = False, ref = ’Martinez’)
Out[ ]: | 4.289075479711235e+23

In the previous example we introduced ratio variable, which is defined as the
ratio $r_{h}/r_{\star}$, where $r_{\star}$ is the half-light radius of the
galaxy. The profiles module also includes brightness, emissivity and flux-
density halo/bulge factor calculators in the Regime “A”, “B” and “C”
approximations introduced in Ref. [26] (see e. g. Eqs. (26), (27)). In the
same way as before, users can decide whether to use mass models shown in Tab.
4 or to input the parameters manually. Particularly interesting for prompt-
emission signals (e. g. gamma rays) are the (differential and solid-angle
integrated) $J$\- and $D$-factor calculators (which are equivalent to the
regime “A” halo/bulge brightness and flux-density factors respectively).

Name | Nr. of modelled galaxies | Models | Reference
---|---|---|---
’1309.2641’ | 19 | (c)NFW, Burkert, Einasto | [11]
’1408.0002’ | 20 | HDZ | [12]
’1608.01749’ | 1 (Draco) | HDZ | [13]
’1706.05481’ | 4 | HDZ | [14]
Table 4:  dSph mass models linked to the diffSph profiles module.

####  3.3.2 The spectra module

This module allows user to compute the spectral function $X(\nu)$
(spectra.synchrotron.X_gen) for a generic (user-given) source function
$\mathcal{S}(E)$. This is done by performing the following integral
numerically

$X(\nu;\tau_{0})=\int_{0}^{\infty}\operatorname{d\\!}E\,\mathcal{S}(E)\hat{X}(E;\tau_{0})\
,$ (18)

where $\hat{X}(E;\tau_{0})$ is the universal spectral “kernel” function (see
appendices). If the user wishes to consider a power law spectrum such as the
’generic’ entry in Table 1, they should rather use pre-computed spectral
function spectra.synchrotron.X.

Pre-computed spectral functions for DM searches are also available. These are
obtained from Monte Carlo parton shower simulations (in the case of DM
annihilation or decay) from PPPC [10].

$X(\nu;\tau_{0})=\frac{4\sqrt{3}e^{3}BE_{\nu}}{m_{e}c^{2}b(E_{\nu})}\mathcal{M}_{K,k\,m_{\rm
DM}/2}\left(\frac{2E_{\nu}}{k\,m_{\rm
DM}c^{2}},\eta(E_{\nu});\delta\right)\quad,\quad k=1,2\ ,$ (19)

where $k=1$ for decaying DM and $k=2$ if WIMPs are considered. The function
$\mathcal{M}_{K,M}(\xi,\eta,\delta)$ is given in the appendices.

####  3.3.3 Utilities

The utils modules contains a comprehensive library of parameters that are
relevant for solving the diffusion-loss equation. In particular all unit
conversions factors and fundamental constants that are used by all modules and
submodules of diffSph can also be found there (see utils.consts).

Further useful tools such as calculators of the signal’s Half Flux Diameter
(HFD) and Full Width at Half Maximum (FWHM) are also implemented. These are
particularly interesting for the estimation of the limits that is mentioned in
3.2.2 and discussed in Ref. [26].

##  4 Examples

The diffSph tool can create signal templates for prompt (gamma rays) and
synchrotron (radio) emission from DM annihilation as shown in Fig. 2.

[width=.85]figs/promptcolor.png

[width=.85]figs/synchcolor.png

Figure 2: Diffuse emission profiles for Draco in the WIMP hypothesis for gamma
rays (upper plot) and radio waves (lower plot).

More specifically, by specifying a particular halo model for a given dwarf
galaxy (e. g. Draco as in the figure), users can easily generate differential
$J$ and $D$ factors as functions of the galaxy’s galactocentric angular
distance. In the upper figure such $J$-factor profiles are shown. Two
parametrizations for the DM halo are considered. On the left hand side the
generalized NFW with the best-fit parameters from Ref. [25] while on the right
hand side, the same template is plotted but using an NFW profile with model
parameters extracted from Ref. [11].

In the lower panels of Fig 2 we show the synchrotron emission profile that
corresponds to the aforementioned DM halo models of the Draco galaxy. While
the $J$ factors only depend on the DM halo mass function, the synchrotron maps
shown in the figure as specific to the frequency of 150 MHz and a 10 GeV WIMP
annihilating primarily into a $\mu^{+}\mu^{-}$ state with the default
propagation parameters.

Another use case of diffSph is the estimation of upper limits on the WIMP
annihilation cross-section or the decay rate of decaying DM candidates given a
noisy radio signal.

[width=.9]figs/sigmav_cvni.pdf

[width=.9]figs/decay_cvni.pdf

Figure 3: CVnI limit estimates on the annihilation cross-section (left panel)
and decay rate (right panel) of DM for several primary final states from non-
observation of synchtrotron radiation.

For that purpose, users can use the sigmav_limest() and decay_rate_limest()
functions contained in the limits module. Fig. 3 shows for several primary
final states the estimated upper limits that can be put on these quantities
for a noisy signal of 100 $\mu$Jy/beam with a beam size of 20” at 150 MHz.

In particular, the following example code can be used to reproduce the
estimated annihilation cross-section and decay rates of DM particles as a
function of their mass for $\mu^{+}\mu^{-}$ final states that shown in red in
Fig. 3:

from diffsph import limits as lim
import numpy as np

mass_grid = np.logspace(0, 4) # WIMP masses

sigmav = [lim.sigmav_limest( nu = .150, rms_noise = 100, beam_size = 20,
galaxy = ’Canes Venatici I’, rad_temp = ’HDZ’, ratio = 1, D0 = 1e27, delta =
’kol’, B = 1, mchi = m, channel = ’mumu’, self_conjugate = True, manual =
False, high_res = False, accuracy = 1, ref = ’1408.0002’ ) for m in mass_grid]
# calculating cross-sections

Gamma = [lim.decay_rate_limest( nu = .150, rms_noise = 100, beam_size = 20,
galaxy = ’Canes Venatici I’, rad_temp = ’HDZ’, ratio = 1, D0 = 1e27, delta =
’kol’, B = 1, mchi = m, channel = ’mumu’, self_conjugate = True, manual =
False, high_res = False, accuracy = 1 ref = ’1408.0002’ ) for m in mass_grid]
# calculating decay rates

###  4.1 Generic power-law spectrum

Another gap that is filled by diffSph is considering generic (non-DM) CRE
sources in dSphs. These have been largely neglected in the literature. In Fig.
4 we show $X(\nu)$ in the case that such power-law CRE injection hypotheses
are assumed.

As expected, the resulting synchrotron spectra also enjoy power-law
dependences $X(\nu)\propto\nu^{-s}$. It is also well known that the spectral
indices $s$ are related to the exponent $\Gamma$ in the CRE source function.
The way $s$ is related to $\Gamma$ is regime dependent, though. In the fast
diffusion situation (regime “C”), the CRE distribution $n_{e}$ is $\propto
E^{-(\Gamma+\delta)}$, which means that $s=(\Gamma+\delta-1)/2$ [27]. In the
opposite limit, where energy losses dominate, $n_{e}\propto E^{-(1+\Gamma)}$
and $s=\Gamma/2$. This is precisely what is observed in all panels of the
figure.

[width = .9]figs/powspecs.pdf

Figure 4:  Spectral function $X(\nu)$ for generic power-law dependent CRE
source functions for several combinations of the $\Gamma$ and $\delta$
parameters.

##  5 Conclusions

We introduced diffSph, a computational assistant (in Python) for the search of
diffuse emissions in dwarf spheroidal galaxies. Besides the most basic
features such as computing spectra and signal profiles given some hypothesis
about the CRE sources, this tool includes several additional use cases. In
particular, the code can be used to obtain constraints on annihilating WIMPs
or decaying DM models. Non-DM scenarios, e. g. point source with power-law CRE
spectra, are also included in diffSph. The halo library is also quite
comprehensive and can be used in order to compute model-specific $J$\- and $D$
factors for the models contained in the dSph catalogue.

The code’s excellent performance is due to a novel analytical approach that
was introduced in Ref. [3] and the incorporation and automatization of a
number of built-in benchmark electron spectra (from e. g. parton shower
algorithms) and spatial distribution functions. In future versions, modules
for ICS and prompt emissions for the benchmark hypotheses presented here will
also be included. We also aim at also enlarging the number of built-in models
and improve further the code’s efficiency and user experience.

This is an open source project licensed under an MIT License.

## Acknowledgements

LG acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum
Universe” – 390833306. MV would like to thanks Kai Otto for his very useful
programming hacks.

##  Appendix A Numerical prescriptions

In contrast to the computation of Eq. (7), which is rather straightforward
once $\mathcal{R}(r)$ is given; the spectral function (8) is much more
involved and requires a dedicated treatment. As a first step, we introduce the
hypothesis independent function $\hat{X}(E)$, which defined in such a way that
$X(\nu)$ can be expressed as
$X(\nu)=\int_{0}^{\infty}\operatorname{d\\!}E\,\mathcal{S}(E)\hat{X}(E)$.
Namely,

$\hat{X}(E)=\frac{4\sqrt{3}e^{3}BE_{\nu}}{m_{e}c^{2}b(E_{\nu})}\mathcal{M}\left(\frac{E_{\nu}}{E},\eta(E_{\nu}),\delta\right)=X_{0}\frac{B}{1+(B/B_{c})^{2}}\frac{1}{E_{\nu}}\mathcal{M}\left(\frac{E_{\nu}}{E},\eta(E_{\nu}),\delta\right)\
,$ (20)

where the constant $K_{0}$ is given by

$X_{0}=\frac{4\sqrt{3}e^{3}}{m_{e}c^{2}b_{0}}\simeq
0.0354\,\textrm{erg}\,\textrm{GHz}^{-1}\,\mu\textrm{G}^{-1}\,\textrm{GeV}\
,\quad(\texttt{utils.consts.X0})\ .$ (21)

The “master” function $\mathcal{M}$ is defined in terms of the pitch-angle
averaged synchrotron-radiation function $\mathcal{F}(y)$ as follows

$\mathcal{M}(\xi,\eta;\delta)=\int_{\xi}^{\infty}{\rm
d}z\,\mathcal{F}(z^{2})\exp\left[-\eta\left(z^{1-\delta}-\xi^{1-\delta}\right)\right]\
,$ (22)

where the $\mathcal{F}$ function reads

$\mathcal{F}(y)=y^{2}\left(K_{4/3}(y)K_{1/3}(y)-\frac{3}{5}\left[K_{4/3}^{2}(y)-K_{1/3}^{2}(y)\right]\right)\
,$ (23)

and $K_{1/3}(y)$, $K_{4/3}(y)$ are the modified Bessel functions with indices
1/3 and 4/3 respectively. As evident by Eq. (20), once the master function is
known for all (positive) values of $\xi$ and $\eta$ (and, of course, $\delta$)
we can determine $\hat{X}$ for generic parameter sets. Nevertheless, the
kernel $\mathcal{M}$ is still a very complicated function. Its numerical
computation can take very long and for some choices of the parameters,
standard integration algorithms fail to converge. In order to overcome this,
the following approximations are very useful and have been implemented in
diffSph.

$\displaystyle M(\xi,\eta,\delta)$ $\displaystyle\to$
$\displaystyle\int_{\xi}^{\infty}{\rm d}x\,\mathcal{F}(x^{2})\quad,\ \eta\to
0\quad(\eta\ll\xi^{-(1-\delta)})$ $\displaystyle M(\xi,\eta,\delta)$
$\displaystyle\to$ $\displaystyle\int_{\xi}^{\infty}{\rm d}x\,f_{0}\times
x^{2/3}\exp[-\eta(x^{1-\delta}-\xi^{1-\delta})]=f_{0}\Gamma\left(\frac{5}{3(1-\delta)}\right)\frac{\eta^{-\frac{5}{3(1-\delta)}}{\rm
e}^{\eta\xi^{1-\delta}}}{1-\delta}\quad,$ $\displaystyle{}\
\eta\to\infty\quad(\eta\gg(1/3)^{-(1-\delta)}\quad,\quad\eta^{1-\delta}\ll(1/3)^{1-\delta})$
$\displaystyle M(\xi,\eta,\delta)$ $\displaystyle\to$
$\displaystyle\int_{u}^{\infty}{\rm
d}x\,\mathcal{F}(x^{2})\delta(\eta(x^{1-\delta}-\xi^{1-\delta}))=\frac{\xi^{\delta}}{(1-\delta)\eta}\mathcal{F}(\xi^{2})\quad,\
\eta\to\infty\quad(\eta\gg\xi^{-(1-\delta)})$

In practice, however, the way diffSph evaluates $\mathcal{M}$ is by
interpolating pre-computed values in the logarithmic spaced grid
$\log\xi\in[-15,+1.2]$, $\log\eta\in[-3,+13]$ and for selected values of the
$\delta$ parameter: $\delta=1/6$, $2/6$ (’kol’), $3/6$ (’kra’) and $4/6$,
where ’kol’ and ’kra’ stand for the standard Kolmogorov [28] and Kraichnan
[29] CR propagation scenarions respectively. The user can certainly turn this
option off by setting the Boolean ’fast_comp’ variable in
(spectra.synchrotron.htX()) to 0. The ranges were chosen in a way that the
extrapolations behave as expected from known asymptotic expressions.

###  A.1 Spectra from parton showers

The functions $\mathcal{M}$ and $\hat{X}$ are by construction hypothesis
independent. Therefore, diffSph users can compute synchrotron spectra for
basically any ansatz for $\mathcal{S}(E)$ by using Eq. (18) (see
spectra.synchrotron.X_gen()). However, for those users interested in testing
specific particle DM scenarios in which the DM annihilates or decays primarily
into two-body standard-model final states ($e^{+}e^{-},b\bar{b},\ldots$),
diffSph contains a number of pre-computed spectral functions $X_{k;K,m_{\rm
DM}}(\nu)$, where $K=e^{+}e^{-},b\bar{b},\ldots$ and the index $k$ specifies
the process: $k=1$ for decay and $k=2$ for annihilation. These are obtained
using the analogue of Eq. (20), where instead of the master function
$\mathcal{M}$ we use

$\mathcal{M}_{K,k\,m_{\rm DM}/2}(\xi,\eta;\delta)=\int_{0}^{1}{\rm
d}x\,\mathcal{M}\left(\frac{\xi}{x},\eta,\delta\right)\,s_{K,k\,m_{\rm
DM}/2}(x)\quad,\quad k=1,2\ ,$ (24)

where $x=E/(k\,m_{\rm DM}/2)$, $s_{K,k\,m_{\rm DM}/2}(x)=(k\,m_{\rm
DM}/2)\times\mathcal{S}_{K,k\,m_{\rm DM}/2}(E)$ and $\mathcal{S}_{K,k\,m_{\rm
DM}/2}(E)$ is the injection function which is obtained from parton shower
algorithm assuming a two-body final state $K$.

###  A.2 Spectra from power-law CRE distribution

For a generic power-law CRE distribution
$\mathcal{S}^{0}(E)=E_{0}^{\Gamma-1}/E^{\Gamma}$, we define

$\displaystyle\mathcal{M}_{\Gamma}(\xi,\eta,\delta)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{\operatorname{d\\!}x}{x^{\Gamma}}\,\mathcal{M}\left(\frac{\xi}{x},\eta,\delta\right)=\xi^{-(\Gamma-1)}\int_{0}^{\infty}\operatorname{d\\!}\xi^{\prime}\,\xi^{\prime\Gamma-2}\mathcal{M}\left(\xi^{\prime},\eta,\delta\right)\
,$ (25) $\displaystyle=$
$\displaystyle\frac{1}{\xi^{\Gamma-1}}\int_{0}^{\infty}\operatorname{d\\!}z\,\mathcal{F}(z^{2})\operatorname{e\\!}{}^{-\eta\,z^{1-\delta}}f_{0}(z;\eta,\delta)\
,$

where
$f_{0}(z;\eta,\delta)=\int_{0}^{z}\operatorname{d\\!}\xi\,\xi^{\Gamma-2}\operatorname{e\\!}{}^{\eta\,\xi^{1-\delta}}$
can be obtained analytically.

##  Appendix B Analytic limits

Depending on how large or small the propagation parameter $\tau_{0}$ is with
respect to the typical energy-loss time scale $t_{\rm loss}\simeq E_{\rm
max}/b(E_{\rm max})$, one can neglect either the diffusion-coefficient term
(regime “A”) or all energy losses (regime “C”) in Eq. (1) [3]. In the
intermediate case, where none of the terms in the diffusion-loss equation can
be neglected, one can make the approximation that all terms with $n>1$ in Eq.
(5) are negligible. Regardless of the regime, Eq. (5) can be further
simplified to

$I_{\nu}^{a}\approx\frac{\mathcal{Q}^{a}}{4\pi}\,\hat{H}_{\rm
reg.}^{a}(\theta)\times X^{a}_{n=1}(\nu)\ ,$ (26)

where

$\hat{H}_{\rm
RA}^{a}(\theta)=\sum_{n=1}^{\infty}\,h^{a}_{n}f_{n}\left(\frac{\theta}{\theta_{h}}\right)\
,\ \hat{H}_{\rm
RB}^{a}(\theta)=h^{a}_{1}f_{1}\left(\frac{\theta}{\theta_{h}}\right)\ {\rm
and}\ \hat{H}_{\rm
RC}^{a}(\theta)=\sum_{n=1}^{\infty}\frac{1}{n^{2}}\,h^{a}_{n}f_{n}\left(\frac{\theta}{\theta_{h}}\right)\
.$ (27)

Alternative (but equivalent) expressions can be found in Ref. [3]. These
expressions are implemented in the profiles module.

## References

  * [1] Ya. B. Zeldovich, A. A. Klypin, M. Yu. Khlopov, and V. M. Chechetkin.  Astrophysical constraints on the mass of heavy stable neutral leptons.  Sov. J. Nucl. Phys., 31:664–669, 1980.
  * [2] M. L. Ahnen et al.  Limits to dark matter annihilation cross-section from a combined analysis of MAGIC and Fermi-LAT observations of dwarf satellite galaxies.  JCAP, 1602(02):039, 2016, 1601.06590.
  * [3] Martin Vollmann.  Universal profiles for radio searches of dark matter in dwarf galaxies.  arXiv e-prints, page arXiv:2011.11947, November 2020, 2011.11947.
  * [4] L. Gajović et al.  Weakly interacting massive particle cross section limits from LOFAR observations of dwarf spheroidal galaxies.  Astron. Astrophys., 673:A108, 2023, 2303.12155.
  * [5] Alex McDaniel, Tesla Jeltema, Stefano Profumo, and Emma Storm.  Multiwavelength Analysis of Dark Matter Annihilation and RX-DMFIT.  JCAP, 09:027, 2017, 1705.09384.
  * [6] Moritz Hütten, Céline Combet, and David Maurin.  CLUMPY v3: $\gamma$-ray and $\nu$ signals from dark matter at all scales.  Comput. Phys. Commun., 235:336–345, 2019, 1806.08639.
  * [7] Sergio Colafrancesco, S. Profumo, and P. Ullio.  Detecting dark matter WIMPs in the Draco dwarf: A multi-wavelength perspective.  Phys. Rev. D, 75:023513, 2007, astro-ph/0607073.
  * [8] Torsten Bringmann, Joakim Edsjö, Paolo Gondolo, Piero Ullio, and Lars Bergström.  DarkSUSY 6 : An Advanced Tool to Compute Dark Matter Properties Numerically.  JCAP, 07:033, 2018, 1802.03399.
  * [9] G. Belanger, F. Boudjema, P. Brun, A. Pukhov, S. Rosier-Lees, P. Salati, and A. Semenov.  Indirect search for dark matter with micrOMEGAs2.4.  Comput. Phys. Commun., 182:842–856, 2011, 1004.1092.
  * [10] Marco Cirelli, Gennaro Corcella, Andi Hektor, Gert Hutsi, Mario Kadastik, Paolo Panci, Martti Raidal, Filippo Sala, and Alessandro Strumia.  PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection.  JCAP, 03:051, 2011, 1012.4515.  [Erratum: JCAP 10, E01 (2012)].
  * [11] Gregory D. Martinez.  A robust determination of Milky Way satellite properties using hierarchical mass modelling.  MNRAS, 451(3):2524–2535, August 2015, 1309.2641.
  * [12] Alex Geringer-Sameth, Savvas M. Koushiappas, and Matthew Walker.  Dwarf galaxy annihilation and decay emission profiles for dark matter experiments.  Astrophys. J., 801(2):74, 2015, 1408.0002.
  * [13] Koji Ichikawa, Miho N. Ishigaki, Shigeki Matsumoto, Masahiro Ibe, Hajime Sugai, Kohei Hayashi, and Shun-ichi Horigome.  Foreground effect on the $J$-factor estimation of classical dwarf spheroidal galaxies.  Mon. Not. Roy. Astron. Soc., 468(3):2884–2896, 2017, 1608.01749.
  * [14] Koji Ichikawa, Shun-ichi Horigome, Miho N. Ishigaki, Shigeki Matsumoto, Masahiro Ibe, Hajime Sugai, and Kohei Hayashi.  Foreground effect on the $J$-factor estimation of ultrafaint dwarf spheroidal galaxies.  Mon. Not. Roy. Astron. Soc., 479(1):64–74, 2018, 1706.05481.
  * [15] Miguel A. Sanchez-Conde, Mirco Cannoni, Fabio Zandanel, Mario E. Gomez, and Francisco Prada.  Dark matter searches with Cherenkov telescopes: nearby dwarf galaxies or local galaxy clusters?  JCAP, 12:011, 2011, 1104.3530.
  * [16] Alan W. McConnachie.  The observed properties of dwarf galaxies in and around the Local Group.  Astron. J., 144:4, 2012, 1204.1562.
  * [17] Nicolas F. Martin, Jelte T. A. de Jong, and Hans-Walter Rix.  A comprehensive Maximum Likelihood analysis of the structural properties of faint Milky Way satellites.  Astrophys. J., 684:1075–1092, 2008, 0805.2945.
  * [18] Beth Willman, Marla Geha, Jay Strader, Louis E. Strigari, Joshua D. Simon, Evan Kirby, Nhung Ho, and Alex Warres.  Willman 1—A Probable Dwarf Galaxy with an Irregular Kinematic Distribution.  AJ, 142(4):128, October 2011, 1007.3499.
  * [19] Troy A. Porter, Gudlaugur Johannesson, and Igor V. Moskalenko.  The GALPROP Cosmic-ray Propagation and Nonthermal Emissions Framework: Release v57.  Astrophys. J. Supp., 262(1):30, 2022, 2112.12745.
  * [20] Michael Korsmeier and Alessandro Cuoco.  Testing the universality of cosmic-ray nuclei from protons to oxygen with AMS-02.  Phys. Rev. D, 105(10):103033, 2022, 2112.08381.
  * [21] Felix Kahlhoefer, Michael Korsmeier, Michael Krämer, Silvia Manconi, and Kathrin Nippel.  Constraining dark matter annihilation with cosmic ray antiprotons using neural networks.  JCAP, 12(12):037, 2021, 2107.12395.
  * [22] Yoann Génolini, Mathieu Boudaud, Marco Cirelli, Laurent Derome, Julien Lavalle, David Maurin, Pierre Salati, and Nathanael Weinrich.  New minimal, median, and maximal propagation models for dark matter searches with Galactic cosmic rays.  Phys. Rev. D, 104(8):083005, 2021, 2103.04108.
  * [23] Carmelo Evoli, Giovanni Morlino, Pasquale Blasi, and Roberto Aloisio.  AMS-02 beryllium data and its implication for cosmic ray transport.  Phys. Rev. D, 101(2):023013, 2020, 1910.04113.
  * [24] Mattia Di Mauro and Martin Wolfgang Winkler.  Multimessenger constraints on the dark matter interpretation of the Fermi-LAT Galactic center excess.  Phys. Rev. D, 103(12):123005, 2021, 2101.11027.
  * [25] Alex Geringer-Sameth, Savvas M. Koushiappas, and Matthew Walker.  DWARF GALAXY ANNIHILATION AND DECAY EMISSION PROFILES FOR DARK MATTER EXPERIMENTS.  The Astrophysical Journal, 801(2):74, mar 2015.
  * [26] M. Vollmann, V. Heesen, T. Shimwell, M. J. Hardcastle, M. Brüggen, G. Sigl and H. Röttgering,  Radio constraints on dark matter annihilation in Canes Venatici I with LOFAR.  Mon. Not. Roy. Astron. Soc. 496 (2020) no.3, 2663-2672 doi:10.1093/mnras/staa1657.
  * [27] George B. Rybicki and Alan P. Lightman.  Radiative processes in astrophysics.  A Wiley Interscience Publication, 1979.
  * [28] A. Kolmogorov.  The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds’ Numbers.  Akademiia Nauk SSSR Doklady, 30:301–305, January 1941.
  * [29] Robert H. Kraichnan.  Inertial-Range Spectrum of Hydromagnetic Turbulence.  Phys. Fluids, 8:1385–1387, 1965.