# The self-confinement of electrons and positrons from dark matter

Marco Regis<EMAIL_ADDRESS>Dipartimento di Fisica, Università di
Torino, Via P. Giuria 1, 10125 Torino, Italy  Istituto Nazionale di Fisica
Nucleare, Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy  Michael
Korsmeier<EMAIL_ADDRESS>The Oskar Klein Centre, Department of
Physics, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden  Gianni
Bernardi  INAF - Istituto di Radioastronomia, via Gobetti 101, 40129 Bologna,
Italy  Department of Physics and Electronics, Rhodes University, PO Box 94,
Grahamstown, 6140, South Africa  South African Radio Astronomy Observatory,
Black River Park, 2 Fir Street, Observatory, Cape Town, 7925, South Africa
Giada Pignataro  Dipartimento di Fisica e Astronomia, Università degli Studi
di Bologna, via P. Gobetti 93/2, 40129 Bologna, Italy  INAF - Istituto di
Radioastronomia, via Gobetti 101, 40129 Bologna, Italy  Javier Reynoso-Cordova
Dipartimento di Fisica, Università di Torino, Via P. Giuria 1, 10125 Torino,
Italy  Istituto Nazionale di Fisica Nucleare, Sezione di Torino, Via P. Giuria
1, 10125 Torino, Italy  Piero Ullio  Scuola Internazionale Superiore di Studi
Avanzati (SISSA), via Bonomea 265, 34136 Trieste, Italy  INFN, Sezione di
Trieste, via Valerio 2, 34127 Trieste, Italy  Institute for Fundamental
Physics of the Universe (IFPU), via Beirut 2, 34151 Trieste, Italy

###### Abstract

Radiative emissions from electrons and positrons generated by dark matter (DM)
annihilation or decay are one of the most investigated signals in indirect
searches of WIMPs. Ideal targets must have large ratio of DM to baryonic
matter. However, such “dark” systems have a poorly known level of magnetic
turbulence, which determines the residence time of the electrons and positrons
and therefore also the strength of the expected signal. This typically leads
to significant uncertainties in the derived DM bounds. In a novel approach, we
compute the self-confinement of the DM-induced electrons and positrons.
Indeed, they themselves generate irregularities in the magnetic field, thus
setting a lower limit on the presence of the magnetic turbulence. We
specifically apply this approach to dwarf spheroidal galaxies. Finally, by
comparing the expected synchrotron emission with radio data from the direction
of the Draco galaxy collected at the Giant Metre Radio Telescope, we show that
the proposed approach can be used to set robust and competitive bounds on WIMP
DM.

##  I Introduction

Dark matter (DM) is a fundamental ingredient in the formation of structures in
our Universe. Yet, its fundamental nature remains elusive. Weakly interacting
massive particles (WIMPs) are one of the most investigated classes of DM
candidates in the literature Silk  _et al._ (2010). The weak interaction
implies that WIMPs in DM halos annihilate in pairs or decay into detectable
species. Thus, information on the WIMP physical properties can be obtained
through indirect searches, namely by studying the astrophysical signals
associated to the annihilation/decay. Among the various channels, a sizable
production of electrons and positrons is a general feature of WIMP models
Regis and Ullio (2008). They in turn generate synchrotron emission in the
magnetized atmosphere of astrophysical structures, which could be observed as
a diffuse radio emission centered on the DM halo distribution Colafrancesco
_et al._ (2006).

The search for faint, diffuse radio emission in targeted galaxies is expected
to undergo a golden era in the forthcoming years. New radio facilities, like
LOFAR, the GMRT, ASKAP and MeerKAT, have started operation. In the next
decade, the SKAO will be the radio telescope providing unmatched sensitivity,
angular resolution and frequency coverage, to search in a deeper and deeper
way for diffuse radio emissions Heald  _et al._ (2020).

It looks appealing to search for WIMP signatures with radio telescopes Regis
_et al._ (2017, 2021). An important obstacle immediately appears, and it is
the presence of radio diffuse emission from “ordinary” astrophysical
processes. How to disentangle it from the possible DM contribution? The
conceptually simplest way to address the issue is to focus on dark, quiescent
systems, where the star formation rate is suppressed, and thus cosmic-ray (CR)
accelerators that would lead to ultra-relativistic electrons and to
synchrotron radiation are absent.

According to numerical simulations, the halos of galaxies are populated by
smaller systems, with many of such subhalos being inefficient in forming
stars. A population of satellites that is observationally detected and that
lacks recent star formation is given by dwarf spheroidal (dSph) galaxies
McConnachie (2012). A key question about a possible synchrotron emission from
dark systems is about the presence of a significant magnetic field. While it
is plausible to assume a coherent component at the level of $\mu$G Regis  _et
al._ (2015), either generated during the history of the object or coming from
the large scale magnetic lines of the host galaxy, the presence of a turbulent
component, in absence of star formation, is more uncertain. The turbulence is
fundamental in setting the confinement scale of the electrons and positrons.
This is a crucial point for the WIMP signal since the objects we just
mentioned are typically quite small, and electrons might travel well outside
the object region before emitting a significant synchrotron radiation.

Such uncertainty on the magnetic properties has led to an uncertainty of
several orders of magnitude in the predicted radio signal in dSphs Regis  _et
al._ (2014), and thus in the bounds on the WIMP annihilation rate Natarajan
_et al._ (2013); Spekkens  _et al._ (2013); Natarajan  _et al._ (2015); Regis
_et al._ (2014, 2017); Kar  _et al._ (2019); Cook  _et al._ (2020); Vollmann
_et al._ (2020); Basu  _et al._ (2021); Gajović  _et al._ (2023). Its
assessment is thus mandatory to determine whether or not dSphs are promising
targets for WIMP radio searches.

In this work, we propose a novel approach, which considers the fact that the
electrons and positrons injected by WIMPs generate themselves magnetic
irregularities and so they induce a certain level of turbulence. We
quantitatively compute such effect and derive under which conditions this
mechanism is relevant to effectively confine $e^{+}-e^{-}$ inside “dark”
systems.

The paper is organized as follows. Section II describes the diffusion and
turbulence equations and their solution (with numerical details in Appendix A
and an approximate analytical solution in Appendix D). The formalism is then
applied to the case of dSph galaxies in Section III, where we also describe
the impact of the magnetic field modeling on the estimate of WIMP-induced
synchrotron radiation. In Section IV, we focus on a specific case, the Draco
dSph, and we derive bounds on WIMP DM by comparing the expected emission with
observations from the GMRT telescope. Details on the DM profile in Draco and
on the radio observations are reported in Appendix C and B. We summarize our
findings in Section V.

##  II Evolution of the electron/positron density and magnetic turbulence

The propagation of CRs in turbulent magnetic fields is an inherently difficult
problem Berezinsky  _et al._ (1990). The charged CRs are deflected in the
magnetic fields. In a constant magnetic field, CRs spiral around the magnetic
field lines. However, as can be shown in linear pertubation theory, the
presence of a (small) turbulent magnetic field on top of the regular field,
randomly changes the angle between the motion of the CR and the regular
magnetic field, effectively leading to the diffusion of the CR along the
magnetic field lines. At the same time, the motion of the charged CRs induces
magnetic fields and stimulates turbulence Cesarsky (1980). The equations
describing the coupling of magnetic fields to the CR plasma are studied in
magnetohydrodynamics. While there are numerical tools (e.g., particle-in-cell
simulations) to solve these equations, the currently available computing power
is not sufficient for the scale of (dwarf) galaxies, given the several orders
of magnitude of separation between the galaxy scale and the CR gyro-radius.
Therefore, we rely on on a more phenomenological-driven method, where the
transport of electrons and positrons is described by diffusion, advection, and
energy losses. Then we assume spherical symmetry. This is a simplifying but
fair approximation, since representative of a model in which a strong gradient
in the CR intensity may manifest in connection to a gradient in the source
function and one in the turbulence distribution. On the other hand, the two
gradients may not be aligned (e.g., a spherical source versus turbulence
flowing along lines of a regular magnetic field being large scale dipole or
toroidal).

With above assumptions, the transport equation is given by Blandford and
Eichler (1987):

$\frac{\partial n_{e}}{\partial t}=\frac{1}{r^{2}}\frac{\partial}{\partial
r}\bigg{[}r^{2}D\frac{\partial n_{e}}{\partial
r}-r^{2}v_{A}n_{e}\bigg{]}+\frac{2v_{A}}{r}\frac{\partial}{\partial
E}\bigg{[}\frac{p}{3}\beta\,c\,n_{e}\bigg{]}-\frac{\partial}{\partial
E}\bigg{[}\dot{E}n_{e}\bigg{]}+q_{\rm CR}\,.$ (1)

Here, $n_{e}$ is the electron (or positron) density in the energy interval
$[E,E+dE]$, $D$ is the spatial diffusion coefficient which may depend on the
radial distance $r$ and energy $E$, $v_{A}$ is the Alfvén velocity that
transports the electrons radially away from the center of the system, $c$ is
the speed of light, $q_{\rm CR}$ is the source term, and $\dot{E}$ describes
the energy losses due to synchrotron emission and inverse Compton
scattering.111We note that there are two different ways to write the advection
term in literature. The divergence term can be partly reshuffled between
electron current (first term) and the adiabatic energy-losses (second term).

In this work we consider electrons and positrons to be injected by the pair
annihilation of DM particles. We assume that there is no other source coming
from astrophysical processes. Thus, the source term is given by:

$q_{\rm CR}(r,E)=\langle\sigma v\rangle_{f}\frac{\rho_{\rm
DM}^{2}(r)}{2\,m_{\rm DM}^{2}}\frac{dN_{e}^{f}}{dE}\;,$ (2)

where $\rho_{\rm DM}$ is the DM energy density, $m_{\rm DM}$ is the DM mass,
$\langle\sigma v\rangle_{f}$ is the velocity-averaged annihilation cross
section into the final state $f$, and ${dN_{e}^{f}}/{dE}$ is the energy
spectrum of the electrons and positrons produced from the annihilation, which
depends on the annihilation channel $f$. We model $\rho_{\rm DM}$ with a
Navarro-Frenk-White (NFW) DM profile Navarro  _et al._ (1996): $\rho_{\rm
NFW}(r)=\rho_{s}\,r_{s}/r/(1+r/r_{s})^{2}$, with two free parameters,
$\rho_{s}$ and $r_{s}$. In Appendix C we discuss also a different choice,
considering a cored profile.

The advection term depends on the Alfvén speed that, for a given magnetic
field, can be estimated from the plasma density $n_{p}$ through $v_{A}\simeq
63\,{\rm km/s}\,\sqrt{10^{-3}{\rm cm^{-3}}/n_{p}}\,B_{0}/\mu{\rm G}$.

Figure 1: Left: Power spectrum of the magnetic turbulence as a function of the
wave number, for a reference DM scenario (detailed in the text) and at three
different distances from the center of the system, 0.2 kpc (dashed), 0.5 kpc
(solid), and 2.0 kpc (dotted). We show the effect of varying the strength of
the regular magnetic field, considering $B_{0}=0.1\,\mu$G (red),
$B_{0}=1.0\,\mu$G (black), and $B_{0}=10\,\mu$G (blue). We remind that
$E\simeq(10^{6}{\rm pc^{-1}}/k)\,(B_{0}/\mu{\rm G})$ GeV. Right: Spatial
diffusion coefficient as a function of energy for the same particle DM
scenario and radial distances of the left panel. We investigate the dependence
on the size and normalization of the DM profile by considering
$\rho_{s}=4\times 10^{7}\,M_{\odot}/{\rm kpc}^{3}$ with $r_{s}=0.1,\,1,\,10$
kpc (red, black, blue) and $\rho_{s}=4\times 10^{8}\,M_{\odot}/{\rm kpc}^{3}$
with $r_{s}=0.1$ kpc (green).

The diffusion of CRs is mediated by scattering processes on the turbulent
magnetic field and we denote the power spectrum of the magnetic turbulence
with $W=W(k,r,t)$. From linear perturbation theory, the relation between the
diffusion coefficient and the turbulence spectrum is given by Berezinsky  _et
al._ (1990):

$D(r,p,t)=\frac{D_{\rm B}(p)4/\pi}{kW(r,k,t)}\;,$ (3)

where $D_{\rm B}(p)=r_{\rm L}(p)c\beta/3$ is the Bohm diffusion coefficient.
For resonant interaction between Alfvén waves and CR scattering, the wave
number $k$ and the momentum $p$ of the CR are related by via the Larmor
radius, $r_{\rm L}(p_{\rm res})=1/k_{\rm res}\simeq 3.3\times 10^{7}\,(E/{\rm
GeV)}\,(\mu{\rm G}/B_{0})$ km.

Turbulence in a galaxy with sizeable star formation rate, such as the Milky
Way, is usually assumed to be injected in the system as an astrophysical
feedback, such as from supernova explosions Blasi (2013); on the other hand,
recent reanalyses - mostly in connection to spectral features in local CR
fluxes - have shown that streaming instabilities of CR themselves may be an
additional and relevant source of turbulence in the Galaxy Blasi  _et al._
(2012); Evoli  _et al._ (2018); Jacobs  _et al._ (2022). Turning to a “dark”
system, like dwarf spheroidal galaxies or isolated DM clumps, in which the
star formation rate is highly suppressed, we discuss here the possibility that
this second effect alone can become sufficiently efficient to give rise to a
diffusive halo. In this scenarios, we assume that turbulence is self-generated
only by the electrons and positrons from DM annihilation. This is clearly a
lower bound for the turbulence (since some rare astrophysical events might
still contribute) and thus a conservative estimate. In particular, if giant
halos around galaxies exist, as possibly suggested by $\gamma$-ray data Gabici
_et al._ (2021), some non-negligible CR and plasma densities might be present
at the dSph location, leading to non-negligible generation of turbulence and
damping. Also, depending on the annihilation final states, there might also be
a significant protons/antiprotons production from DM, which would enhance the
turbulence. Including these additional contributions is beyond the scope of
this work, which focuses on determining an inescapable lower limit. Let us
also stress again that we are assuming spherical symmetry and that the
diffusion is occurring along magnetic field lines that are radially directed.
If the magnetic field configuration of a real target including all these
contributions has a different symmetry, our description has to be accordingly
modified.

We model the cascading of turbulence and the advection with Alfven waves, and
the turbulence power spectrum evolves according to Jones (1993); Miller and
Roberts (1995)

$\frac{\partial W}{\partial t}=\frac{\partial}{\partial
k}\left[D_{kk}(W)\frac{\partial W}{\partial
k}\right]-\frac{1}{r^{2}}\frac{\partial}{\partial r}(r^{2}v_{A}W)+\Gamma_{\rm
CR}(n_{e},k)W\,.$ (4)

The first term on the RHS describes the turbulent cascade by a diffusion in
$k$, the second term accounts for advection with Alfvén velocity $v_{A}$, and
the third term includes the source of turbulence given by resonant streaming
instability. We assume Kolmogorov turbulence, setting the wave number
diffusion coefficient to:

${D_{kk}(W)=c_{k}v_{\rm A}k^{7/2}\sqrt{W}},$

with the numerical factor $c_{k}=0.052$.

The growth rate of the turbulence through streaming instability is given by
Evoli  _et al._ (2018); Jacobs  _et al._ (2022):

$\Gamma_{\rm
CR}=\frac{4\pi\,c\,v_{A}}{3\,k\,W(k)\,B_{0}^{2}/(8\pi)}\left[\beta(p)\,p^{4}\left|\frac{\partial
f}{\partial r}\right|\right]_{p=p_{res}}\,,$ (5)

where $B_{0}$ is the coherent magnetic field of the system,
$4\pi\,p^{2}\,f\,dp=n_{e}\,dE$ and we consider $E\simeq p$ since the electrons
are relativistic. In Eq. (4) we neglected damping. There are different damping
terms discussed in the literature. A nice overview is given for example in
Ref. Recchia  _et al._ (2022). For the systems under consideration, i.e., with
low stellar turbulence, the most important process would be the ion-neutral
damping, due to momentum transfer or charge exchange of ions and neutrons. We
model it adding a term $-\Gamma_{\rm IND}W$ to the RHS of Eq. (4) with
$\Gamma_{\rm IND}$ as described in Ref. Jacobs  _et al._ (2022). By
considering reasonable plasma densities for dwarf galaxies, we found that
damping is not relevant in these systems. Thus, for simplicity, we neglect the
damping term in the rest of the paper.

We solved the two coupled Eqs. (1) and (4) numerically, with a Crank-Nicolson
and explicit scheme, respectively. We refer the reader to Appendix A for
details on the numerical implementation.

Figure 2: Power spectrum of the magnetic turbulence as a function of the wave
number, for the same reference DM scenario of Fig. 1, but with the
annihilation cross section increased by one order of magnitude, $\langle\sigma
v\rangle=10^{-23}\,\rm{cm^{3}/s}$. We show the effect of varying the Alfvén
speed, considering $v_{A}=6.3$ km/s (red), $v_{A}=63$ km/s (black),
$v_{A}=630$ km/s (blue), and $v_{A}=6300$ km/s (green) and report the results
at three different distances from the center of the system, 0.2 kpc (dashed),
0.5 kpc (solid), and 2.0 kpc (dotted).

To understand the level of turbulence that can be provided by DM and the
associated diffusion coefficient, we show some illustrative solutions in Figs.
1 and 2. For these plots, we define a reference model with the following
ingredients: $m_{\rm DM}=100$ GeV, $\langle\sigma
v\rangle=10^{-24}\,\rm{cm^{3}/s}$, annihilation into $b\bar{b}$,
$\rho_{s}=4\times 10^{7}\,M_{\odot}/{\rm kpc}^{3}$ and $r_{s}=1$ kpc (the
latter two numbers are chosen to have a typical DM density of dSph galaxies,
see later the case of Draco). We also set $B_{0}=1.0\,\mu$G and $v_{A}=63$
km/s, unless otherwise specified.

From both figures, one can notice that the level of turbulence generated by DM
increases towards the center. The main reason behind this behaviour is simply
the higher DM density (and thus source term) at smaller radii. In the left
panel of Fig. 1, we see that as $B_{0}$ increases, $W$ decreases. This can be
understood from Eq. (5) with $B_{0}$ being at the denominator. In other words,
the level of turbulence is set by the ratio between the turbulent magnetic
field induced by DM and the coherent magnetic field, and, for a higher value
of the latter, one gets the same ratio (i.e., the same turbulence) only for a
higher DM source term.

In the right panel, we show the diffusion coefficient obtained from Eq. (3),
for the same reference scenario, and investigating the impact of the size and
normalization of the DM mass density. First we vary $r_{s}$, keeping fixed
$\rho_{s}$ (red, black and blue curves) and see that the confinement becomes
more and more effective (i.e., the diffusion coefficient decreases) as the
source size increases. Note that we do not consider diffusion coefficient
larger than $10^{31}{\rm cm^{2}/s}$. The reason is related to the fact that
Eq. (1) is an effective equation, in particular it is not covariant. This
means that for large $D$ (and fixed size of the diffusion region), it can
describe propagation of particles faster than light. To prevent this issue, we
impose $D\leq 10^{31}{\rm cm^{2}/s}$. For such maximal value of $D$, the
picture is essentially equivalent to a free escape Regis  _et al._ (2015), so
in other words we treated all the cases providing $D>10^{31}{\rm cm^{2}/s}$ as
free escape.

In the right panel of Fig. 1 we also show a case (green line) with a “small”
size of the source, $r_{s}=0.1$ kpc (compared to $\sim$ 1 kpc of a dwarf
galaxy), but with $\rho_{s}$ increased by one order of magnitude with respect
to the reference scenario, $\rho_{s}=4\times 10^{8}\,M_{\odot}/{\rm kpc}^{3}$.
This shows that the effect we are describing in this work, with a focus on
dwarf galaxies, can be relevant in general for dark clumps, that, on the other
hand, must have large DM over-densities.

Fig. 2 shows that the effect of advection is marginal, for the systems under
consideration. It can have an impact only for extremely large values of the
Alfvén speed, $v_{A}\gtrsim 10^{3}$ km/s. This can be understood by estimating
the time-scales associated to the three physical processes included in the
equations:

$\tau_{\rm diff}\simeq 10^{15}\,{\rm s}\left(\frac{L}{{\rm
kpc}}\right)^{2}\frac{10^{28}{\rm cm^{2}/s}}{D}\;,\;\tau_{\rm adv}\simeq
3\times 10^{15}\,{\rm s}\,\frac{L}{{\rm kpc}}\frac{10\,{\rm
km/s}}{v_{A}}\;,\;\tau_{\rm loss}\simeq 4\times 10^{16}\,{\rm s}\,\frac{{\rm
GeV}}{E}\left[1+0.1\left(\frac{B}{\mu{\rm G}}\right)^{2}\right]\,.$ (6)

These estimates tell us that for $D\gtrsim 10^{28}{\rm cm^{2}/s}$, advection
and energy losses are typically subdominant, and that for $D\simeq 10^{29}{\rm
cm^{2}/s}$, which is the case for different benchmark cases shown in our
plots, one needs $v_{A}\sim 10^{3}$ km/s to make advection the dominant
transport mechanism. In Fig. 2, we increased the annihilation cross section to
$\langle\sigma v\rangle=10^{-23}\,\rm{cm^{3}/s}$ in order to have a picture
with lower $D$, i.e., to have scenarios with more relevant impact from
advection.

##  III The case of dwarf spheroidal galaxies and synchrotron radiation

In this Section, we set the DM profile to the case of a prototypical dSph
galaxy. In particular, we consider an NFW profile with $\rho_{s}=4\times
10^{7}\,M_{\odot}/{\rm kpc}^{3}$ and $r_{s}=1$ kpc, that are the best-fit
values found in the case of Draco, see Appendix C. Here Draco is taken as a
benchmark dSph, and considering other classical dSph galaxies would produce
very similar plots. In Figs. 3 and 4 we show again solutions of Eqs. (1) and
(4) (the latter reported in terms of $D$ instead of $W$ using Eq. (3)).

First, we explore the dependence on the strength of the coherent magnetic
field. We assume this magnetic component to be present independently from DM.
Plausible estimates point towards a strength at the level of $\mu$G or
fraction of $\mu$G, see, e.g., discussion in (Regis  _et al._ , 2015). In Fig.
3, we take $\langle\sigma v\rangle=10^{-24}\,\rm{cm^{3}/s}$ (and all the other
parameters as in the reference model detailed above), and see that for
$B_{0}<1\,\mu$G, the self-confinement is important. Indeed for $B_{0}\sim$
fraction of $\mu$G, we get low values of $D$, and electron densities that are
orders of magnitude larger than in the free escape case (as it can be seen by
comparing red and blue curves in the right panel).

In Fig. 4, we set $B_{0}=1\,\mu$G and explore the dependence on the
annihilation cross section. It is clear that for larger $\langle\sigma
v\rangle$ the source term is larger, and so the generation of turbulence is
more effective, leading to a lower diffusion coefficient and larger electron
density. It is important to note that the effect deviates from a simple linear
scaling. There are two saturation effects: for too low level of injection from
DM, one always gets a “free-escape” scenario, whilst on the contrary for very
large injection one ends up in a scenario where the residence time is mostly
determined by advection (see discussion of Eq. (6)), which is independent from
$\langle\sigma v\rangle$.

Figure 3: Spatial diffusion coefficient (left) and equilibrium number density
of electrons (right) as a function of energy for a reference scenario
(detailed in the text) and investigating the dependence on the regular
magnetic field. We show three different radial distances, 0.2 kpc (dashed),
0.5 kpc (solid), and 2.0 kpc (dotted), and three different magnetic field
strengths $B_{0}=0.1\,\mu$G (red), $1\,\mu$G (black), and $10\,\mu$G (blue,
where Eq. (3) would give $D>10^{31}{\rm cm^{2}/s}$).

From the determination of the $e^{+}-e^{-}$ number density $n_{e}$ that we
just discussed, we can now compute the synchrotron emissivity at a given
frequency $\nu$ by folding $n_{e}$ with the total radiative emission power
$P_{syn}$ (Rybicki and Lightman, 1979):

$j_{syn}(\nu,r)=\int dE\,P_{syn}(r,E,\nu)\,n_{e}(r,E)\;.$ (7)

The unpolarized synchrotron power depends on the total magnetic field $B$,
i.e., on the sum of the coherent (that we call $B_{0}$) plus random (that we
call $\delta B$) magnetic fields. We compute the latter from $\delta
B=B_{0}\,\sqrt{\int dk\,W(k)}$. For models in the ball-park of the bounds we
are going to derive in the next Section, $\delta B\ll B_{0}$, and thus
$B\simeq B_{0}$, which is reassuring, since we are working under the quasi-
linear approximation.

It is well known that the synchrotron power roughly scales with $B^{2}$ with a
cutoff for $B\lesssim\nu/\rm{GHz}\,(15\,\rm{GeV}/E)^{2}$. Since $B$ enters in
the determination of $n_{e}$, the dependence of the emissivity from the
magnetic field is not trivial (and depends on $m_{\chi}$, $\langle\sigma
v\rangle$ and on the annihilation channel).

We show some illustrative cases in Fig. 5, taking again a typical dSph, Draco,
and a radio frequency of 650 MHz. We vary the magnetic field strength between
0.1 and 1 $\mu$G, and consider two DM masses, 10 and 100 GeV, and two
annihilation channels, $b\bar{b}$ and $\tau^{+}\tau^{-}$. We set the cross
section so to have non-negligible confinement (more precisely we choose the
value corresponding to the limit, computed in the next Section, for
$B_{0}=1\,\mu$G).

We see that for $m_{\chi}=100$ GeV, i.e., for scenarios efficiently injecting
electrons at energies far from the cutoff of the synchrotron power, an
increase of $B_{0}$ from 0.1 to 1 $\mu$G does not lead to significantly larger
fluxes, or even to lower emissivity, see the $\tau^{+}\tau^{-}$ case. Namely,
the reduction in $n_{e}$ for larger $B_{0}$, due to a smaller confinement
time, is compensating the increase in the synchrotron power, and the
dependence of the signal on the $B_{0}$ strength is mild. On the other hand,
for lower DM masses, the scaling with $B_{0}$ can be dramatic since the peak
of the distribution of the injected $e^{+}-e^{-}$ is now at an energy lower
than the one corresponding to the cutoff in the synchrotron power.

Figure 4: Spatial diffusion coefficient (left) and equilibrium number density
of electrons (right) as a function of energy for a reference scenario
(detailed in the text) and investigating the dependence on the annihilation
cross section. We show three different radial distances, 0.2 kpc (dashed), 0.5
kpc (solid), and 2.0 kpc (dotted), and three different annihilation cross
sections $\langle\sigma v\rangle=10^{-24}\,\rm{cm^{3}/s}$ (black),
$\langle\sigma v\rangle=10^{-23}\,\rm{cm^{3}/s}$ (red), and $\langle\sigma
v\rangle=10^{-22}\,\rm{cm^{3}/s}$ (blue).

##  IV Bounds from radio observations of the Draco dSph galaxy

Our final goal is to derive bounds on WIMP DM by comparing the radio emission
induced by DM in dSph galaxies with observations.

Figure 5: Emissivity at 650 MHz as a function of the radial distance from the
dSph center. In the left panel we consider the $b\bar{b}$ annihilation
channel, and set the annihilation cross section to $\langle\sigma
v\rangle=3.3\times 10^{-25}\,\rm{cm^{3}/s}$ ($1.5\times
10^{-24}\,\rm{cm^{3}/s}$) for $m_{\chi}=10$ GeV (100 GeV). On the right panel
we show the case of annihilation into $\tau^{+}\tau^{-}$ with $\langle\sigma
v\rangle=1.4\times 10^{-25}\,\rm{cm^{3}/s}$ ($2.9\times
10^{-24}\,\rm{cm^{3}/s}$) and $m_{\chi}=10$ GeV (100 GeV).

The flux density measured by a telescope can be estimated as

$S_{th}(\nu,\theta_{0})=\int
d\phi\,d\theta\,\sin\,\theta\,\mathcal{G}(\theta,\phi,\theta_{0})\int
ds\,\frac{j_{syn}(\nu,r(s,\theta,\phi))}{4\pi}\;,$ (8)

where $s$ labels the coordinate along the line of sight, $\phi$ and $\theta$
are the angular coordinates, $\theta_{0}$ is the direction of observation, and
$\mathcal{G}$ gives synthesized beam of the telescope.

In the evaluation of Eq. (8), we need to specify five ingredients, two of them
related to the WIMP microscopic properties, i.e., the particle mass $m_{\chi}$
and the annihilation rate $\langle\sigma v\rangle$, and three related to the
dSph, namely, the DM spatial distribution, the regular magnetic field, and the
plasma density (the latter to provide the Alfven velocity, see above). We
derive bounds in the plane $\langle\sigma v\rangle$ versus $m_{\chi}$.

We constrain the DM profile $\rho$ via Jeans analysis and using the dispersion
velocities of the dSph stellar component. In particular, we consider the
parametric form of the NFW profile Navarro  _et al._ (1996), and derive the
likelihood $\mathcal{L}_{Jeans}$ associated to the parameters
$(\rho_{s},r_{s})$ describing the distribution. For more details, see Appendix
C, where we also consider the case of a cored profile described by the Burkert
parametrization.

For what concerns the regular magnetic field and plasma density, we cannot
proceed with a data-driven determination since only upper limits can be found
in the literature. Thus we bracket the uncertainty by considering two
different scenarios for each quantity, one where we assume values smaller than
in the host galaxy (i.e., our Milky Way) but not too far-away, taking
$B_{0}=1\,\mu$G and $n_{g}=10^{-3}\,{\rm cm}^{-3}$, and one where instead we
take a significantly smaller estimate but still above the cosmological values,
i.e., $B_{0}=0.1\,\mu$G and $n_{g}=10^{-5}\,{\rm cm}^{-3}$.

With this modeling at hand we can now evaluate $S_{th}$ for a given dSph. In
particular, we compute it for the Draco dSph and compare to recent data
collected with the upgraded Giant Metrewave Radio Telescope (uGMRT) at 550-750
MHz. For more details on the observations, see Appendix B.

Assuming a Gaussian likelihood, we can compare expected diffuse signals and
observed data through Regis  _et al._ (2015):

$\mathcal{\tilde{L}}_{diff}=e^{-\chi^{2}/2}\;\;\;{\rm
with}\;\;\;\chi^{2}=\frac{1}{N_{pix}^{FWHM}}\sum_{i=1}^{N_{pix}}\left(\frac{S_{th}^{i}-S_{obs}^{i}}{\sigma_{rms}^{i}}\right)^{2}\;,$
(9)

where $i$ denotes the pixel in the Draco radio image, $S_{obs}^{i}$ is the
observed flux density, $\sigma_{rms}^{i}$ is the r.m.s. error, $N_{pix}$ is
the total number of pixels in the area under investigations, that is chosen to
be a circle of $30^{\prime}$ of radius, and $N_{pix}^{FWHM}$ is the number of
pixels within the uGMRT synthesized beam.

We add to the theoretical term $S_{th}$ a spatially flat term $S_{flat}$ that
is included in the fit to account for a possible offset in the zero-level
calibration of the map, and then we define a likelihood $\mathcal{L}_{diff}$
which depends only on the DM parameters by profiling out $S_{flat}$ from the
likelihood in Eq. (9).

Then, combining the Jeans analysis and uGMRT data, we can define, at any given
mass $m_{\chi}$, a global likelihood:

$\mathcal{L}(\langle\sigma
v\rangle,\rho_{s},r_{s})=\mathcal{L}_{diff}(\langle\sigma
v\rangle,\rho_{s},r_{s})\times\mathcal{L}_{Jeans}(\rho_{s},r_{s})\;.$ (10)

To derive the bounds, we assume that $\lambda_{c}\langle\sigma
v\rangle=-2\ln[\mathcal{L}(\langle\sigma
v\rangle,\rho_{s}^{\,lbf},r_{s}^{\,lbf})/\mathcal{L}(\langle\sigma
v\rangle^{b.f.},\rho_{s}^{\,gbf},r_{s}^{\,gbf})]$ follows a
$\chi^{2}$-distribution with one d.o.f. and with one-sided probability given
by $P=\int^{\infty}_{\sqrt{\lambda_{c}}}d\chi\,e^{-\chi^{2}/2}/\sqrt{2\,\pi}$,
where $\langle\sigma v\rangle^{b.f.}$ denotes the best-fit value for the
annihilation rate at that specific WIMP mass. The superscript $gbf$ indicates
the global best-fit (i.e., taking $\langle\sigma v\rangle=\langle\sigma
v\rangle^{b.f.})$), whilst $lbf$ denotes the best-fit of $\rho_{s}$ and
$r_{s}$ for that given $\langle\sigma v\rangle$. The 95% C.L. upper limit on
$\langle\sigma v\rangle$ at mass $m_{\chi}$ is obtained from
$\lambda_{c}=2.71$.

Results are shown in Fig. 6 with black lines for three different annihilation
channels, $b\bar{b}$ (left), $\tau^{+}\tau^{-}$ (right, dashed) and
$\mu^{+}\mu^{-}$ (right, dotted). The reference scenario assumes
$B_{0}=1\,\mu$G and $n_{p}=10^{-3}\,{\rm cm}^{-3}$ (i.e., $v_{A}=63$ km/s). In
the left panel, we show the effect of varying the magnetic field and advection
velocity, by considering $B_{0}=0.1\,\mu$G and $n_{p}=10^{-5}\,{\rm cm}^{-3}$
(we focus on the $b\bar{b}$ channel for the sake of brevity). To understand
the values of $v_{A}$, we remind that $v_{A}\simeq 63\,{\rm
km/s}\,\sqrt{10^{-3}{\rm cm^{-3}}/n_{p}}\,B_{0}/\mu{\rm G}$.

Figure 6: Bounds at 95% C.L. on the annihilation rate as a function of the DM
mass. Our reference scenario has $B_{0}=1\,\mu$G and $n_{p}=10^{-3}\,{\rm
cm^{-3}}$ (which implies $v_{A}\simeq 63$ km/s). The colored bands show the
impact of the uncertainty associated to the determination of the DM density,
assumed to follow an NFW profile. Blue is for $b\bar{b}$, red for
$\tau^{+}\tau^{-}$ and green for $\mu^{+}\mu^{-}$. In the left panel, we show
also alternative scenarios, where $B_{0}=0.1\,\mu$G and/or
$n_{p}=10^{-5}\,{\rm cm^{-3}}$.

The left panel shows the main achievement of this work. We devised a self-
consistent method to derive bounds on WIMP DM from radio observations of dSph.
While in the current state-of-the-art one deals with several orders of
magnitude of uncertainty for the derived bounds (see, e.g., Fig. 5 in (Regis
_et al._ , 2014)), which are due to the unknown description of the dSph
interstellar medium, in this work we show there is an “irreducible” bound,
that is only a factor $\mathcal{O}(1)$ uncertain (for masses above 20-30 GeV).

The two ingredients we need to include in the model (and that are not data-
driven, like on the contrary the parameters of the DM profile) are the
advection and the regular magnetic field. The impact of advection is limited.
The value of $B_{0}$, instead, significantly affects both the synchrotron
power and the confinement time, but the effects on the signal go in the
opposite directions, partially canceling each other. Indeed, the larger is
$B_{0}$ the larger is $P_{syn}$ (thus enhancing the signal), but also the
larger is $D$ (which depletes the signal), as we already discussed in the
previous Section. This implies that the different scenarios differ only by a
factor $\mathcal{O}(1)$ in the expected signal and so in the bound.

The only exception to this argument is for scenarios with low $B_{0}$, low
$m_{\chi}$ and soft channels of annihilation. Indeed, in this case the peak of
the synchrotron power corresponds to an energy well above the peak of the
$e^{+}-e^{-}$ spectrum induced by WIMP annihilations. We show it in the case
of $b\bar{b}$, where this effect is maximal, being the spectrum of emission
quite soft, whilst for leptonic channels the spectrum is harder and the effect
is more limited, also at low masses.

For a fair comparison of Fig. 6 with bounds in other analyses of Draco, let us
mention that the determination of the DM profile described in Appendix C leads
to a J-factor at $0.5^{\circ}$ which is in agreement with recent estimates
Petac  _et al._ (2018); Alvarez  _et al._ (2020), but a factor of $\sim 3$
lower than in the computation of (Bonnivard  _et al._ , 2015), which has been
widely used in the literature. To understand the impact of the uncertainty in
the determination of the DM profile, we include the colored bands in Fig. 6.
They show bounds derived with a procedure that is slightly different from the
one outlined above, namely, the band reports the bounds obtained by taking all
the $\rho_{s}$ and $r_{s}$ value within their 95% C.L. region (derived from
the Jeans analysis).

##  V Conclusions

In this work we investigated the generation of magnetic turbulence due to the
injection of electrons and positrons from DM annihilation. In objects where
the star formation rate is significant, this mechanism is typically
subdominant with respect to the magnetic irregularities given by, e.g.,
supernova explosions. On the other hand, we found that in “dark” and DM dense
systems, it can be the driving mechanism, and it can lead to a non-negligible
self-confinement of the electrons and positrons.

We quantitatively assessed the effect by solving numerically a system of two
coupled differential equations that describe the evolution of the electron
number density and of the turbulence.

As a concrete application, we derived solutions in the case of dSph galaxies,
which are one of the prime targets in DM indirect searches. After determining
the equilibrium density of electrons and positrons, we computed the associated
synchrotron radiation. This emission peaks in the radio regime and the search
for faint, diffuse radio sources is expected to significantly progress in the
forthcoming years thanks to SKAO and its precursors.

Our findings suggest that the self-generated turbulence implies a non-
negligible lower limit to the residence time of the electrons in dSphs. The
dependence of this effect on the strength of the coherent magnetic field
$B_{0}$ partially cancels the dependence of the synchrotron power on $B_{0}$,
implying low uncertainty in the predicted DM-induced synchrotron emission.
This allows to derive robust bounds, something that represents a significant
improvement with respect to the current state-of-the-art where bounds are
instead heavily dependent on the magnetic assumptions (see, e.g., discussion
in (Regis  _et al._ , 2014)).

Finally, we used data of the Draco dSph collected at the uGMRT at 550-750 MHz,
and compared them to the expected DM synchrotron emission computed within the
self-confinement framework. We found that the proposed approach can already
allow to set robust and competitive bounds on WIMP DM, summarized in Fig. 6.
Near future radio observations with the SKAO and its precursors have the
capability to significantly further push down the derived constraints.

## Acknowledgements

We would like to thank C. Evoli and S. Recchia for useful discussions, and the
Institute for Fundamental Physics of the Universe (IFPU), Trieste, where the
early phase of this work was carried out, for the hospitality and support.

MR and JR acknowledge support from the project “Theoretical Astroparticle
Physics (TAsP)” funded by the INFN, from ‘Departments of Excellence 2018-2022’
grant awarded by the Italian Ministry of Education, University and Research
(miur) L. 232/2016 and from the research grant ‘From Darklight to Dark Matter:
understanding the galaxy/matter connection to measure the Universe’ No.
20179P3PKJ funded by miur. MK is supported by the Swedish Research Council
under contracts 2019-05135 and 2022-04283 and the European Research Council
under grant 742104. PU is supported by the research grant “The Dark Universe:
A Synergic Multi-messenger Approach” number 2017X7X85K under the program PRIN
2017 funded by the The Italian Ministry of Education, University and Research
(MIUR), and by the European Union’s Horizon 2020 research and innovation
program under the Marie Sklodowska-Curie grant agreement No 860881-HIDDeN.

## References

  * Silk  _et al._ (2010) J. Silk _et al._ ,  _Particle Dark Matter: Observations, Models and Searches_, edited by G. Bertone (Cambridge Univ. Press, Cambridge, 2010).
  * Regis and Ullio (2008) M. Regis and P. Ullio, Phys. Rev. D 78, 043505 (2008), arXiv:0802.0234 [hep-ph] .
  * Colafrancesco  _et al._ (2006) S. Colafrancesco, S. Profumo, and P. Ullio, Astron. Astrophys. 455, 21 (2006), arXiv:astro-ph/0507575 .
  * Heald  _et al._ (2020) G. Heald _et al._ (SKA Magnetism Science Working Group), Galaxies 8, 53 (2020), arXiv:2006.03172 [astro-ph.GA] .
  * Regis  _et al._ (2017) M. Regis, L. Richter, and S. Colafrancesco, JCAP 07, 025 (2017), arXiv:1703.09921 [astro-ph.HE] .
  * Regis  _et al._ (2021) M. Regis _et al._ , JCAP 11, 046 (2021), arXiv:2106.08025 [astro-ph.HE] .
  * McConnachie (2012) A. W. McConnachie, Astron. J. 144, 4 (2012), arXiv:1204.1562 [astro-ph.CO] .
  * Regis  _et al._ (2015) M. Regis, L. Richter, S. Colafrancesco, S. Profumo, W. J. G. de Blok, and M. Massardi, Mon. Not. Roy. Astron. Soc. 448, 3747 (2015), arXiv:1407.5482 [astro-ph.GA] .
  * Regis  _et al._ (2014) M. Regis, S. Colafrancesco, S. Profumo, W. J. G. de Blok, M. Massardi, and L. Richter, JCAP 10, 016 (2014), arXiv:1407.4948 [astro-ph.CO] .
  * Natarajan  _et al._ (2013) A. Natarajan, J. B. Peterson, T. C. Voytek, K. Spekkens, B. Mason, J. Aguirre, and B. Willman, Phys. Rev. D 88, 083535 (2013), arXiv:1308.4979 [astro-ph.CO] .
  * Spekkens  _et al._ (2013) K. Spekkens, B. S. Mason, J. E. Aguirre, and B. Nhan, ApJ 773, 61 (2013), arXiv:1301.5306 [astro-ph.CO] .
  * Natarajan  _et al._ (2015) A. Natarajan, J. E. Aguirre, K. Spekkens, and B. S. Mason, arXiv e-prints , arXiv:1507.03589 (2015), arXiv:1507.03589 [astro-ph.CO] .
  * Kar  _et al._ (2019) A. Kar, S. Mitra, B. Mukhopadhyaya, T. R. Choudhury, and S. Tingay, Phys. Rev. D 100, 043002 (2019), arXiv:1907.00979 [astro-ph.HE] .
  * Cook  _et al._ (2020) R. H. W. Cook, N. Seymour, K. Spekkens, N. Hurley-Walker, P. J. Hancock, M. E. Bell, J. R. Callingham, B. Q. For, T. M. O. Franzen, B. M. Gaensler, L. Hindson, M. Johnston-Hollitt, A. D. Kapińska, J. Morgan, A. R. Offringa, P. Procopio, L. Staveley-Smith, R. B. Wayth, C. Wu, and Q. Zheng, ”Mon. Not. Roy. Astron. Soc.” 494, 135 (2020), arXiv:2003.06104 [astro-ph.GA] .
  * Vollmann  _et al._ (2020) M. Vollmann, V. Heesen, T. W. Shimwell, M. J. Hardcastle, M. Brüggen, G. Sigl, and H. J. A. Röttgering, ”Mon. Not. Roy. Astron. Soc.” 496, 2663 (2020), arXiv:1909.12355 [astro-ph.HE] .
  * Basu  _et al._ (2021) A. Basu, N. Roy, S. Choudhuri, K. K. Datta, and D. Sarkar, ”Mon. Not. Roy. Astron. Soc.” 502, 1605 (2021), arXiv:2101.04925 [astro-ph.HE] .
  * Gajović  _et al._ (2023) L. Gajović _et al._ , (2023), arXiv:2303.12155 [astro-ph.CO] .
  * Berezinsky  _et al._ (1990) V. S. Berezinsky, S. V. Bulanov, V. A. Dogiel, and V. S. Ptuskin,  _Astrophysics of cosmic rays_ , edited by V. L. Ginzburg (1990).
  * Cesarsky (1980) C. J. Cesarsky, Ann. Rev. Astron. Astrophys. 18, 289 (1980).
  * Blandford and Eichler (1987) R. Blandford and D. Eichler, Phys. Rept. 154, 1 (1987).
  * Navarro  _et al._ (1996) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996), arXiv:astro-ph/9508025 [astro-ph] .
  * Blasi (2013) P. Blasi, Astronomy and Astrophysics Reviews 21, 70 (2013), arXiv:1311.7346 [astro-ph.HE] .
  * Blasi  _et al._ (2012) P. Blasi, E. Amato, and P. D. Serpico, Phys. Rev. Lett. 109, 061101 (2012), arXiv:1207.3706 [astro-ph.HE] .
  * Evoli  _et al._ (2018) C. Evoli, P. Blasi, G. Morlino, and R. Aloisio, Phys. Rev. Lett. 121, 021102 (2018), arXiv:1806.04153 [astro-ph.HE] .
  * Jacobs  _et al._ (2022) H. Jacobs, P. Mertsch, and V. H. M. Phan, JCAP 05, 024 (2022), arXiv:2112.09708 [astro-ph.HE] .
  * Gabici  _et al._ (2021) S. Gabici, S. Recchia, F. Aharonian, and V. Niro, Astrophys. J. 914, 135 (2021), arXiv:2101.05016 [astro-ph.HE] .
  * Jones (1993) T. W. Jones, ApJ 413, 619 (1993).
  * Miller and Roberts (1995) J. A. Miller and D. A. Roberts, ApJ 452, 912 (1995).
  * Recchia  _et al._ (2022) S. Recchia, D. Galli, L. Nava, M. Padovani, S. Gabici, A. Marcowith, V. Ptuskin, and G. Morlino, Astron. Astrophys. 660, A57 (2022), arXiv:2106.04948 [astro-ph.HE] .
  * Rybicki and Lightman (1979) G. B. Rybicki and A. P. Lightman,  _Radiative processes in astrophysics_ (1979).
  * Petac  _et al._ (2018) M. Petac, P. Ullio, and M. Valli, JCAP 12, 039 (2018), arXiv:1804.05052 [astro-ph.GA] .
  * Alvarez  _et al._ (2020) A. Alvarez, F. Calore, A. Genina, J. Read, P. D. Serpico, and B. Zaldivar, JCAP 09, 004 (2020), arXiv:2002.01229 [astro-ph.HE] .
  * Bonnivard  _et al._ (2015) V. Bonnivard _et al._ , Mon. Not. Roy. Astron. Soc. 453, 849 (2015), arXiv:1504.02048 [astro-ph.HE] .
  * Vollmann (2020) M. Vollmann, arXiv e-prints , arXiv:2011.11947 (2020), arXiv:2011.11947 [astro-ph.HE] .
  * Perley and Butler (2017) R. A. Perley and B. J. Butler, ApJS 230, 7 (2017), arXiv:1609.05940 [astro-ph.IM] .
  * Intema  _et al._ (2009) H. T. Intema, S. van der Tol, W. D. Cotton, A. S. Cohen, I. M. van Bemmel, and H. J. A. Röttgering, A&A 501, 1185 (2009), arXiv:0904.3975 [astro-ph.IM] .
  * Offringa  _et al._ (2014) A. R. Offringa, B. McKinley, N. Hurley-Walker, F. H. Briggs, R. B. Wayth, D. L. Kaplan, M. E. Bell, L. Feng, A. R. Neben, J. D. Hughes, J. Rhee, T. Murphy, N. D. R. Bhat, G. Bernardi, J. D. Bowman, R. J. Cappallo, B. E. Corey, A. A. Deshpande, D. Emrich, A. Ewall-Wice, B. M. Gaensler, R. Goeke, L. J. Greenhill, B. J. Hazelton, L. Hindson, M. Johnston-Hollitt, D. C. Jacobs, J. C. Kasper, E. Kratzenberg, E. Lenc, C. J. Lonsdale, M. J. Lynch, S. R. McWhirter, D. A. Mitchell, M. F. Morales, E. Morgan, N. Kudryavtseva, D. Oberoi, S. M. Ord, B. Pindor, P. Procopio, T. Prabu, J. Riding, D. A. Roshi, N. U. Shankar, K. S. Srivani, R. Subrahmanyan, S. J. Tingay, M. Waterson, R. L. Webster, A. R. Whitney, A. Williams, and C. L. Williams, MNRAS 444, 606 (2014), arXiv:1407.1943 [astro-ph.IM] .
  * Walker  _et al._ (2015) M. G. Walker, E. W. Olszewski, and M. Mateo, MNRAS 448, 2717 (2015), arXiv:1503.02589 [astro-ph.GA] .
  * Read  _et al._ (2018) J. I. Read, M. G. Walker, and P. Steger, Mon. Not. Roy. Astron. Soc. 481, 860 (2018), arXiv:1805.06934 [astro-ph.GA] .
  * Binney (1980) J. Binney, Monthly Notices of the Royal Astronomical Society 190, 873 (1980), https://academic.oup.com/mnras/article-pdf/190/4/873/18194217/mnras190-0873.pdf .
  * Binney and Mamon (1982) J. Binney and G. A. Mamon, Monthly Notices of the Royal Astronomical Society 200, 361 (1982), https://academic.oup.com/mnras/article-pdf/200/2/361/9403099/mnras200-0361.pdf .
  * Flewelling  _et al._ (2020) H. A. Flewelling, E. A. Magnier, K. C. Chambers, J. N. Heasley, C. Holmberg, M. E. Huber, W. Sweeney, C. Z. Waters, A. Calamida, S. Casertano, X. Chen, D. Farrow, G. Hasinger, R. Henderson, K. S. Long, N. Metcalfe, G. Narayan, M. A. Nieto-Santisteban, P. Norberg, A. Rest, R. P. Saglia, A. Szalay, A. R. Thakar, J. L. Tonry, J. Valenti, S. Werner, R. White, L. Denneau, P. W. Draper, K. W. Hodapp, R. Jedicke, N. Kaiser, R. P. Kudritzki, P. A. Price, R. J. Wainscoat, S. Chastel, B. McLean, M. Postman, and B. Shiao, ApJS 251, 7 (2020), arXiv:1612.05243 [astro-ph.IM] .
  * Navarro (1996) J. F. Navarro, Symposium - International Astronomical Union 171, 255–258 (1996).
  * Burkert (1995) A. Burkert, The Astrophysical Journal 447, L25 (1995).
  * Read and Steger (2017) J. I. Read and P. Steger, MNRAS 471, 4541 (2017), arXiv:1701.04833 [astro-ph.GA] .
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, PASP 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Genina  _et al._ (2020) A. Genina, J. I. Read, C. S. Frenk, S. Cole, A. Benítez-Llambay, A. D. Ludlow, J. F. Navarro, K. A. Oman, and A. Robertson, MNRAS 498, 144 (2020), arXiv:1911.09124 [astro-ph.GA] .
  * Reynoso-Cordova  _et al._ (2022) J. Reynoso-Cordova, M. Regis, and M. Taoso, JCAP 10, 038 (2022), arXiv:2203.13735 [astro-ph.GA] .

##  Appendix A Numerical implementation

In this Appendix, we report details on the numerical solution of Eqs. (1) and
(4).

###  A.1 Electron/positron density

Equation (1) has been finite-differenced by means of the Crank-Nicolson
scheme. It has built on the case without advection, described in the Appendix
of (Regis  _et al._ , 2015) (see also Vollmann (2020) for a semi-analytical
solution). The full discretized equation is given by:

$\displaystyle\frac{n^{t+1}_{i,j}-n^{t}_{i,j}}{\Delta t}$
$\displaystyle=\frac{1}{2}\bigg{[}L^{t}_{i,j}+L^{t+1}_{i,j}\bigg{]}$ (11)

where we defined $n(t_{t},r_{i},\tilde{E}_{j})=n^{t}_{i,j}$ and consider a
linear grid in space $r(r_{i})=r_{i}$ and a logarithmic grid in energy
$\tilde{E}_{j}=\log\Big{(}E_{j}\Big{)}$. The operators in Eq. (11) are given
by (where each line corresponds to each of the five different terms on the RHS
of Eq. (1)):

$\displaystyle L^{t}_{i,j}=$ $\displaystyle\frac{1}{r_{i}^{2}(\Delta
r)^{2}}\left[\Big{(}r_{i+\frac{1}{2}}\Big{)}^{2}D_{i+\frac{1}{2},j}\Big{(}n^{t}_{i+1,j}-n^{t}_{i,j}\Big{)}\right.$
(12)
$\displaystyle\qquad\qquad\qquad\qquad-\left.\Big{(}r_{i-\frac{1}{2}}\Big{)}^{2}D_{i-\frac{1}{2},j}\Big{(}n^{t}_{i,j}-n^{t}_{i-1,j}\Big{)}\right]$
(13) $\displaystyle-v_{A}\bigg{[}\frac{2}{r_{i}}n^{t}_{i,j}+\frac{1}{\Delta
r}\Big{(}n^{t}_{i,j}-n^{t}_{i-1,j}\Big{)}\bigg{]}$ (14)
$\displaystyle-\frac{1}{E_{j}\Delta\tilde{E}}\bigg{[}\left(\frac{dE}{dt}(E_{j+1})\right)n^{t}_{i,j+1}-\left(\frac{dE}{dt}(E_{j})\right)n^{t}_{i,j}\bigg{]}$
(15) $\displaystyle+\Big{(}q^{(\rm CR)}\Big{)}^{t}_{i,j}.$ (16)

having defined $r_{i\pm\frac{1}{2}}=(r_{i}+r_{i\pm 1})/2$.

As mentioned, we choose a linear grid in $r$. This is because it gives a
faster convergence than a log grid. We checked the results are however
independent from this choice. We set $L=5\,r_{s}$ (where $r_{s}$ is the scale
radius of the DM profile) and define

$\displaystyle\Delta_{r}$ $\displaystyle=\frac{L}{n_{r}-0.5}\qquad,\qquad
r_{i}=\Big{(}i+\frac{1}{2}\Big{)}\Delta_{r}\quad\bigg{|}i\in\\{0,1,...,n_{r}-1\\}$

We impose the following boundary conditions:

$\displaystyle\partial_{r}n\Big{(}r=0,E,t\Big{)}$
$\displaystyle=0\qquad\qquad\qquad\longrightarrow$ $\displaystyle
n^{t}_{-1,j}$ $\displaystyle=n^{t}_{0,j}$ (17) $\displaystyle
n\Big{(}r=L,E,t\Big{)}$ $\displaystyle=0\qquad\qquad\qquad\longrightarrow$
$\displaystyle n^{t}_{n_{r}-1,j}$ $\displaystyle=0$ (18) $\displaystyle
n\Big{(}r,E=E_{\rm min},t\Big{)}$
$\displaystyle=0\qquad\qquad\qquad\longrightarrow$ $\displaystyle n^{t}_{i,0}$
$\displaystyle=0$ (19) $\displaystyle n\Big{(}r,E=E_{\rm max},t\Big{)}$
$\displaystyle=\text{continuous}\qquad\longrightarrow$
$\displaystyle\mathrm{d}n^{t}_{i,n_{E}-1}$
$\displaystyle=\mathrm{d}n^{t}_{i,n_{E}-2}$ (20)

Note that we include a fictitious value in $r$ below grid,
$r_{-1}=-\Delta_{r}/2$. In this way the boundary condition is exactly at
$r=0$. This is one of the motivations for choosing a linear grid since with a
log grid one can only impose the boundary condition at a small $r\neq 0$.

###  A.2 Magnetic turbulence

Equation (4) is solved using the explicit scheme.
First we define $Q^{(W)}=\Gamma_{\rm CR}(n)\,W/v_{A}$. Then, since all terms
on the RHS of Eq. (4) are proportional to $v_{A}$, the stationary solution is
independent of $v_{A}$, and by redefining $\tau=v_{A}k\,t$, we can write Eq.
(4) as:

$\frac{\partial W}{\partial\tau}=c_{k}\frac{1}{k}\frac{\partial}{\partial
k}\left[k^{\frac{7}{2}}\sqrt{W}\frac{\partial W}{\partial
k}\right]-\frac{1}{k}\frac{1}{r^{2}}\frac{\partial}{\partial
r}(r^{2}W)+\frac{1}{k}Q^{(W)}.$ (21)

The replacement $\tau=v_{A}k\,t$ and its dependence on $k$ might be unexpected
at first glance. We note that this replacement does not allow to easily
reconstruct the solutions for finite times $t$. But in the steady state
solution ($t\rightarrow\infty$ and $\tau\rightarrow\infty$), where the LHS of
the equation disappears, the solution for $W$ is identical. We verified this
explicitly for a few examples. The replacement allows faster convergence at
small $k$.

The full discretized equation is given by (where each line corresponds to each
term in Eq. (21)) :

$\displaystyle W^{t+1}_{i,j}$ $\displaystyle=W^{t}_{i,j}$ (22)
$\displaystyle+\frac{c_{k}\Delta\tau}{k_{j}^{2}(\Delta
l)^{2}}\left[\Big{(}k_{j+\frac{1}{2}}\Big{)}^{\frac{5}{2}}\sqrt{\left(\frac{W^{t}_{i,j}+W^{t}_{i,j+1}}{2}\right)}\Big{(}W^{t}_{i,j+1}-W^{t}_{i,j}\Big{)}-\Big{(}k_{j-\frac{1}{2}}\Big{)}^{\frac{5}{2}}\sqrt{\left(\frac{W^{t}_{i,j}+W^{t}_{i,j-1}}{2}\right)}\Big{(}W^{t}_{i,j}-W^{t}_{i,j-1}\Big{)}\right]$
(23)
$\displaystyle-\frac{\Delta\tau}{k_{j}}\bigg{[}\frac{2}{r_{i}}W^{t}_{i,j}+\frac{1}{\Delta{r_{i}}}\Big{(}W^{t}_{i,j}-W^{t}_{i-1,j}\Big{)}\bigg{]}$
(24) $\displaystyle+\frac{\Delta\tau}{k_{j}}\Big{(}Q^{(W)}\Big{)}^{t}_{i,j}.$
(25)

where we have $W(\tau_{t},r_{i},l_{j})=W^{t}_{i,j}$, and used a linear grid in
space $r(r_{i})=r_{i}$ (the same as for the density equation (11)) and a log
grid in the wave number $k(l_{j})=k_{j}=\exp(l_{j})$. We note that the wave
number is related to energy since $E/{\rm GeV}=9.2\times
10^{5}\,(B_{0}/\mu{\rm G})({\rm pc}^{-1}/k)$. However, we solve Eq. (4)
extending the gird in energy to avoid impacts from the boundary conditions
(see below) imposed on $W$.

We defined $k_{j\pm\frac{1}{2}}=\exp[(l_{j}+l_{j\pm 1})/2]$.

Finally, we impose the following boundary conditions:

$\displaystyle\partial_{k}W\Big{(}r,k=k_{\rm min},t\Big{)}$ $\displaystyle=0$
$\displaystyle\longrightarrow$ $\displaystyle W^{t}_{i,0}$
$\displaystyle=W^{t}_{i,1}$ (26) $\displaystyle W\Big{(}r=L,k,t\Big{)}$
$\displaystyle=0$ $\displaystyle\longrightarrow$ $\displaystyle
W^{t}_{n_{r}-1,j}$ $\displaystyle=0$ (27) $\displaystyle W\Big{(}r,k=k_{\rm
max},t\Big{)}$ $\displaystyle=0$ $\displaystyle\longrightarrow$ $\displaystyle
W^{t}_{i,n_{k}}$ $\displaystyle=0$ (28)
$\displaystyle\partial_{r}W\Big{(}r=0,k,t\Big{)}$ $\displaystyle=0$
$\displaystyle\longrightarrow$ $\displaystyle W^{t}_{i,-1}$
$\displaystyle=W^{t}_{i,0}$ (29)

We apply the Neumann boundary condition Eq. (26) because the asymptotic
solution for $k_{\rm min}\rightarrow 0$ in the diffusion-only scenario is a
constant value $W$. However, this asymptotic behavior is broken by the
advection term as can be seen for example in Fig. 1. The solution goes quickly
to zero as $k$ becomes smaller than the resonant value corresponding to
$E=m_{\rm DM}$. Since we set $k_{\rm min}$ much smaller than this value the
boundary condition becomes irrelevant. Also for the boundary condition at
$k_{\rm max}$ in Eq. (28) we could have imposed the asymptotic behavior of the
diffusion-only equation. On the other hand, we explicitly checked that this
boundary condition affects the solution only down to $k_{\rm
max}/\mathcal{O}(1)$, and we chose $k_{\rm max}$ corresponding to an energy
much smaller than the ones used to derive the DM bounds. Thus, also the
boundary condition at $k_{\rm max}$ becomes irrelevant and do not impact our
results.

For the results shown in this paper, we set $E_{\rm min}=0.1$ GeV, $E_{\rm
max}=m_{\rm DM}$, $k_{\rm min}=100\,{\rm pc}^{-1}$ and $k_{\rm
max}=10^{8}\,{\rm pc}^{-1}$.

##  Appendix B Radio observations and data reduction

Figure 7: Image of the field centred on the Draco dSph galaxy with the uGMRT
at 650 MHz. The angular resolution is $8^{\prime\prime}\times
3^{\prime\prime}$ and the rms noise is 10 $\mu$Jy beam-1.

Observations of the Draco dSph galaxy were carried out with the uGMRT at band
4. Details are reported in Table 1. Observations were preceded by a 20 minute
observation of the bandpass calibrator 3C 286, assumed to be 21.4 Jy at 550
MHz, with a spectral index $\alpha=-0.85$ (Perley and Butler, 2017). Data
reduction was carried out with the Source Peeling and Atmospheric Modeling
(SPAM, Intema  _et al._ , 2009) pipeline. Given the wide (200 MHz) band,
observations were split in four sub-bands of 50 MHz each and calibrated
separately. The pipeline carries out the excision of radio frequency
interference, followed by initial calibration and imaging of each sub-band
individually. Images were further deconvolved individually in order to obtain
a sky model for selfcalibration. Calibration solutions were specifically
determined for bright, offending sources that were, therefore, accurately
subtracted. The calibrated sub-bands were eventually jointly imaged with
WSClean v3.1 (Offringa  _et al._ , 2014), using a robust = 0 weighting scheme,
yielding to a $8^{\prime\prime}\times 3^{\prime\prime}$ synthesized beam. The
deconvolved image is shown in Fig. 7, with a 10 $\mu$Jy beam-1 rms noise.

The high resolution image was used to model all the compact sources brighter
than 50 $\mu$Jy beam-1 which were subtracted from the visibilities. Residual
visibilities were then imaged at lower resolution $32^{\prime\prime}\times
15^{\prime\prime}$ and few more sources were identified and subtracted to
obtain a final 32 $\mu$Jy beam-1 rms noise for the image that has been used to
derive the bounds shown in the main text.

Table 1: Observation details. Target | RAJ2000 | DECJ2000 | Frequency range | Channel width | Integration time
---|---|---|---|---|---
|  |  | (MHz) | (kHz) | (hours)
Draco | $17^{\rm h}20^{\rm m}12.401$ | $+57^{\circ}20^{\prime}55^{\prime\prime}$ | $550-750$ | 97.656 | 3.1

##  Appendix C Jeans analysis of the Draco dSph galaxy

As mentioned on Section III, in order to constrain the spatial distribution of
the DM density, we performed a Jeans analysis using Draco’s stellar
kinematics. We used the line-of-sight (LOS) velocities sample from Ref. Walker
_et al._ (2015). In addition to LOS velocities, this dataset contains
information on effective temperatures, metallicities and surface gravities of
stars, which has been used to select the final sample of LOS velocities,
following the same procedure as in Ref. Read  _et al._ (2018). The LOS
projected spherical radial velocity is Binney (1980); Binney and Mamon (1982)

$\sigma_{\rm{LOS}}^{2}(R)=\frac{2}{\Sigma(R)}\int_{R}^{\infty}\left(1-\beta\frac{R^{2}}{r^{2}}\right)\nu\sigma_{r}^{2}\frac{rdr}{\sqrt{r^{2}-R^{2}}},$
(30)

where $\Sigma(R)$ is the projected luminous surface density, $\beta$ is the
orbital velocity anisotropy, $\nu(r)$ is the stellar density and
$\sigma_{r}(r)$ is the radial velocity distribution which can be obtained
solving the spherical Jeans equation

$\frac{1}{\nu(r)}\frac{\partial}{\partial
r}(\nu(r)\sigma_{r}^{2})+\frac{2\beta(r)\sigma_{r}^{2}}{r}=-\frac{GM(<r)}{r^{2}},$
(31)

where $M(<r)$ the mass content up to a radius $r$.

Besides the kinematic data of the stars, a photometric sample is needed in
order to estimate the stellar density distribution in Draco. It is taken from
Ref. Flewelling  _et al._ (2020). Then we model the stellar surface brightness
density as a sum of three (projected) Plummer spheres Binney and Mamon (1982)

Figure 8: Posterior distributions of the DM parameters
$\log_{10}(\rho_{s}/(M_{\odot}/\rm{kcp}^{3}))$ and $\log_{10}(r_{s}/\rm{kpc})$
from equation (33). Left panel shows the results of the MCMC for the NFW
profile. Right panel shows the results in the Burkert case.

$\nu(r)=\sum_{j=1}^{3}\frac{3M_{j}}{4\pi
a_{j}^{3}}\left(1+\frac{r^{2}}{a_{j}^{2}}\right)^{-5/2}\;\;\;,\;\;\;\Sigma(R)=\sum_{j=1}^{3}\frac{M_{j}}{\pi
a_{j}^{2}}\left(1+\frac{R^{2}}{a_{j}^{2}}\right)^{-2},$ (32)

with $M_{j}$ and $a_{j}$ being parameters that we include in the likelihood
and estimate from the data. The contribution to the enclosed mass in Eq. (31)
is computed through the integration of the stellar density, whilst the DM part
is modeled either through a NFW Navarro (1996) or a Burkert Burkert (1995)
profile,

$\rho_{\rm{NFW}}(r)=\frac{\rho_{s}}{\left(\frac{r}{r_{s}}\right)\left(1+\frac{r}{r_{s}}\right)^{2}}\;\;\;,\;\;\;\rho_{\rm{Burkert}}(r)=\frac{\rho_{s}}{\left(1+\frac{r}{r_{s}}\right)\left(1+\left(\frac{r}{r_{s}}\right)^{2}\right)}.$
(33)

In addition to the kinematic and photometric data, we also make use of higher
velocity moments called Virial Shape Parameters (VSPs) which in principle can
help alleviate the mass-anisotropy degeneracy Read and Steger (2017), and are
given by:

$v_{s1}=\frac{2}{5}\int_{0}^{\infty}GM\nu(r)\left(5-2\beta(r)\right)\sigma_{r}^{2}rdr=\int_{0}^{\infty}\Sigma(R)\left\langle
v_{\rm{LOS}}^{4}\right\rangle RdR,$ (34a)
$v_{s2}=\frac{4}{35}\int_{0}^{\infty}GM\nu(r)\left(7-6\beta(r)\right)\sigma_{r}^{2}r^{3}dr=\int_{0}^{\infty}\Sigma(r)\left\langle
v_{\rm{LOS}}^{4}\right\rangle R^{3}dR,$ (34b)

being the RHS of the equation obtained trough the kinematic data and the LHS
the prediction of the model. With all these ingredients we perform an MCMC
analysis using the python package emcee Foreman-Mackey  _et al._ (2013)
implemented through the non-parametric code Gravsphere Read and Steger (2017)
and the python wrapper pyGravSphere Genina  _et al._ (2020). For more details
on the method, see Ref. Reynoso-Cordova  _et al._ (2022). Our results for the
posterior distributions of the DM parameters
$\log_{10}(\rho_{s}/[M_{\odot}/\mathrm{kpc}^{3}])$ and
$\log_{10}(r_{s}/[\mathrm{kpc}])$ are shown in Fig. 8.

In Fig. 9, we show a comparison between the bounds derived in Fig. 6, adopting
an NFW profile, with the ones obtained for a Burkert profile. The difference
between the two cases is limited. This can be understood from: a) the J-factor
over $30^{\prime}$ of the two models is very similar; b) the size of the
Burkert core is significantly smaller than the diffusive region (and thus it
has a limited impact in reducing the turbulence generation); c) the spatial
diffusion is partially washing out the central cusp of the NFW profile.

Figure 9: Similar to Fig. 6 but comparing the bounds derived with the NFW and
Burkert profiles.

##  Appendix D Estimate of the diffusion coefficient with analytical
approximations

In this Appendix we aim at deriving a simple (and clearly approximated)
expression for the diffusion coefficient induced by the self-generation
mechanism described in this work.

By neglecting advection and taking the stationary limit of Eq. (4), one
obtains

$\frac{\partial}{\partial k}\left[D_{kk}(W)\frac{\partial W}{\partial
k}\right]=-\Gamma_{\rm CR}(n)W\;.$ (35)

Let us now consider spectral and spatial features of the equilibrium
distribution of the $e^{+}-e^{-}$ to be separable and assume a power-law for
the energy spectrum, namely:

$n_{e}=10^{-12}\,\rm{GeV^{-1}\,cm^{-3}}\,\tilde{N}_{0}\,\left(\frac{E}{\rm
GeV}\right)^{-3+\delta}\,\tilde{f}(r/r_{*})\;,$ (36)

where $\tilde{N}_{0}$ is just a normalization constant, $\delta$ gives the
spectral index of the energy distribution of $n_{e}$ and $\tilde{f}$ describes
its spatial profile, that we will write in terms of a typical scale $r_{*}$
(which could be, e.g., the half-light radius or the NFW scale radius), where
it is normalized $\tilde{f}(r=r_{*})=1$. Putting together Eqs. (35) and (36),
taking the term $p^{4}\,f$ at the resonance, and after some simple numerical
evaluations, one gets

$-\frac{\partial}{\partial k}\left[k^{3}(k\,W)^{\gamma}\frac{\partial
W}{\partial k}\right]=3.3\times
10^{-8}\,\tilde{N}_{0}\,\left[\frac{(10^{-6}\,\rm{pc})^{-1}}{k}\right]^{\delta}\left(\frac{B_{0}}{\mu\rm{G}}\right)^{-3+\delta}\left|\frac{\partial\tilde{f}}{\partial
r/(100\,\rm{pc})}\right|\;.$ (37)

By describing the turbulence with $W=W_{0}\,k^{\beta}$, and plugging it into
Eq. (37), it follows $\beta=-1-\delta/(\gamma+1)$ and the solution reads:

$W=W_{0}\,k^{-1-\delta/(\gamma+1)}\;\;\;\rm{with}\;\;\;W_{0}=\left(\frac{3.3\times
10^{-8}}{(1-\delta)\,[1+\delta/(\gamma+1)]}\,\frac{\tilde{N}_{0}}{(10^{-6}\,\rm{pc})^{\delta}}\left(\frac{B_{0}}{\mu\rm{G}}\right)^{-3+\delta}\left|\frac{\partial\tilde{f}}{\partial
r/(100\,\rm{pc})}\right|\right)^{1/(1+\gamma)}\;.$ (38)

Notice that, as already discussed in the main text, the turbulence scales
inversely with the strength of the coherent magnetic field, and directly with
the spatial derivative of the distribution function. Finally, from Eq. (3),
namely $D\simeq c\,r_{L}/[3\,k_{res}\,W(k_{res})]$, we arrive at the
expression of the diffusion coefficient:

$D=3.3\times 10^{22}\rm{cm^{2}/s}\,\left(\frac{3.3\times
10^{-8}\,\tilde{N}_{0}}{(1-\delta)\,[1+\delta/(\gamma+1)]}\,\left|\frac{\partial\tilde{f}}{\partial
r/(100\,\rm{pc})}\right|\right)^{-\frac{1}{1+\gamma}}\,\left(\frac{E}{\rm
GeV}\right)\,\left(\frac{B_{0}}{\mu\rm{G}}\right)^{\frac{3-\delta}{1+\gamma}-1}\,\left(\frac{k_{res}}{(10^{-6}\,\rm{pc})^{-1}}\right)^{\frac{\delta}{\gamma+1}}\;.$
(39)

Taking now a typical example, i.e., $\delta=0$, and $\gamma=1/2$ (Kolmogorov),
the estimate of the diffusion coefficient becomes

$D=3\times
10^{27}\rm{cm^{2}/s}\,\left(\tilde{N}_{0}\,\left|\frac{\partial\tilde{f}}{\partial
r/(100\,\rm{pc})}\right|\right)^{-2/3}\,\left(\frac{E}{\rm
GeV}\right)\,\left(\frac{B_{0}}{\mu\rm{G}}\right)\;.$ (40)

From the above expression, it looks clear that for magnetic fields at $\mu$G
level and GeV particles, the effect of self-confinement can be relevant if the
$e^{+}-e^{-}$ energy density is $\gtrsim
10^{-12}\,\rm{GeV/cm^{3}}=10^{-3}\,\rm{eV/cm^{3}}$ and with a significant
spatial variation of the distribution.

In the above estimate, we neglected advection. This is justified only for
$D\gg 10^{27}\rm{cm^{2}/s}$ (see main text), otherwise the full equation has
to be solved. In order to have an estimate of $q_{CR}$ needed to have a
significant turbulence production and so of the DM annihilation rate needed,
one can solve Eq. (2) by again taking the stationary limit and neglecting
advection and energy losses (since in the physical cases considered here the
associated time-scales are longer than the diffusion time-scale), i.e.,
solving $D\,\partial n_{e}/\partial r=-\int_{0}^{r}dr^{\prime}\,r^{\prime
2}\,q_{\rm CR}(r^{\prime})/r^{2}$. It can be easily verified that this
approach leads to similar results as the ones presented in the main text with
a more sophisticated numerical treatment.