# New 511 keV line data provides strongest sub-GeV dark matter constraints

Pedro De la Torre Luque<EMAIL_ADDRESS>The Oskar Klein
Centre, Department of Physics, Stockholm University, Stockholm 106 91, Sweden
Polytechnic School, Universidad San Pablo-CEU, CEU Universities. 28668
Boadilla del Monte, Spain  Shyam Balaji<EMAIL_ADDRESS>Physics
Department, King’s College London, Strand, London, WC2R 2LS, United Kingdom
Joseph Silk<EMAIL_ADDRESS>Institut d’Astrophysique de Paris, UMR 7095 CNRS &
Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France  Department
of Physics and Astronomy, The Johns Hopkins University, 3400 N. Charles
Street, Baltimore, MD 21218, U.S.A.  Beecroft Institute for Particle
Astrophysics and Cosmology, University of Oxford, Keble Road, Oxford OX1 3RH,
U.K.

###### Abstract

We study sub-GeV dark matter (DM) particles that may annihilate or decay into
Standard Model (SM) particles producing an electron-positron cascade that
results in positronium bound state formation after energy losses. This
comprises an exotic injection component in the Milky Way that leaves an
imprint in the 511 keV photon line due to the decay of positronium into two
photons. In this work, we use $\sim 16$ yr of SPI spectrometer data from the
INTEGRAL satellite to constrain DM properties. We include three major
novelties in our study: i) we account for positron diffusion and propagation,
as well as positron losses due to annihilation in flight and other energy
losses, ii) we include the free electron density suppression away from the
Galactic plane and iii) we derive limits for decaying DM for the first time
with SPI data. We show that the predicted longitude and latitude profiles
change significantly for different DM masses, contrary to what has previously
been assumed. In addition, we find that the limits derived from this new set
of SPI data are the strongest on sub-GeV DM to date across almost the entire
DM mass range considered (from MeV to a few GeV), excluding cross-sections
down to $10^{-32}$ cm3 s${}^{-1}\,\,(\text{for}\,\,m_{\chi}\sim
1\,\text{MeV})\lesssim\langle\sigma v\rangle\lesssim 10^{-26}$ cm3
s${}^{-1}\,\,(m_{\chi}\sim 5\,\text{GeV})$ and lifetimes up to
$10^{29}\,\textrm{s}\,(m_{\chi}\sim 1\,\text{MeV})\lesssim\tau\lesssim
10^{27}\,\textrm{s}$ ($m_{\chi}\sim 5$ GeV), whilst considering best-fit
cosmic ray (CR) propagation and diffusion parameters. These limits surpass
even the most stringent complementary cosmological and astrophysical limits
over most of the mass range considered.

##  I Introduction

Our Galaxy emits a distinctive $\gamma$-ray emission line at 511 keV that
indicates the presence of positrons and electrons annihilating into two
photons in the interstellar medium (ISM). This line has been measured by
various instruments, such as the SPI spectrometer on the INTEGRAL satellite
Siegert  _et al._ (2016) and the COSI balloon telescope Kierans  _et al._
(2020), and has been reviewed in several works Prantzos  _et al._ (2011);
Kierans  _et al._ (2019); Siegert (2023). The $\gamma$-ray line has two
components, from the galactic disk and bulge respectively. The bulge component
is very narrow, bright and focused around the galactic center (GC), with a
photon flux of about $10^{-3}$ cm-2 sec-1 Siegert  _et al._ (2016). This has
previously been used to set constraints on possible positron sources, since
the line is attributed to the annihilation of (thermal) $e^{+}e^{-}$ into
$\gamma\gamma$ via (para-)positronium formation.

Figure 1: Diffuse $511$ keV line emission longitude profile in the latitude
bin $-10.75^{\circ}<b<10.25^{\circ}$ from DM direct annihilation into
electron-positron pairs for DM masses $m_{\chi}$ ranging from $1$ (cyan), $10$
(green), $100$ (blue) and $1000$ MeV (black) respectively. To facilitate the
comparison with SPI data (shown in red), the emission for each mass is scaled
by the factor indicated in the legend

.

This requires sources injecting positrons in the regions where the emission is
observed, with injection energy that cannot be larger than a few tens of MeV
Beacom and Yuksel (2006); Boehm (2009); Sizun  _et al._ (2007). The
$\gamma$-ray line from the Galactic disk seems compatible with positron
production from the decay of isotopes that are made by stars (see e.g. Refs.
Prantzos  _et al._ (2011); Bartels  _et al._ (2018)), but the source of the
positrons in the bulge is still a mystery (and considered “one of the most
intriguing problems in high energy astrophysics” Prantzos  _et al._ (2011)).
Some recent solutions proposed for the source of bulge positrons include low-
mass $X$-ray binaries Bartels  _et al._ (2018) and neutron star mergers Fuller
_et al._ (2019).

The source of the $511$ keV $\gamma$-ray line from the galactic bulge is still
unclear and could be connected to dark matter (DM), which is expected to be
abundant around the Galactic Center (GC). This possibility has been explored,
at least since Ref. Boehm  _et al._ (2004a). However, not all DM models are
appropriate. For example, DM decay is not compatible with the shape of the
signal, as seen by Refs. Vincent  _et al._ (2012); Ascasibar  _et al._ (2006),
though these works did not use this incompatibility to place strong limits on
DM. Moreover, from studies of the related in-flight positron annihilation and
final state radiation (FSR) (also called inner bremsstrahlung) the positron
injection energy must be around the MeV scale Beacom  _et al._ (2005); Boehm
and Uwer (2006); Beacom and Yuksel (2006); Boehm (2009), which means that if
DM is to account for the emission, the particle mass must be below some tens
of MeV, unless there are some other intermediate steps in the annihilation
process. Such a low DM mass was thought to be inconsistent with cosmological
data, but recent studies have shown that it can be reconciled with
observations if there is also a small amount of neutrino injection from DM
annihilation in the early universe Escudero (2019); Sabti  _et al._ (2020). We
will use these recent findings to motivate study of light DM scenarios that
may produce the 511 keV line.

A seminal work by Ref. Vincent  _et al._ (2012) showed that the amount of DM
annihilation needed to explain the observed 511 keV line emission in the
galactic bulge (using about $8$ years of data collected by INTEGRAL/SPI in an
energy bin of 5 keV width), was parameterized by the approximate relation

$\displaystyle\langle\sigma v\rangle_{511}\simeq 5\times
10^{-31}\Big{(}\frac{m_{\chi}}{3\leavevmode\nobreak\
\text{MeV}}\Big{)}^{2}\,\frac{\text{cm}^{3}}{\text{sec}},$ (1)

for a Navarro-Frenk-White (NFW) profile DM density distribution for the best
fit parameters considered in Ref. Vincent  _et al._ (2012).

However, the previous studies have assumed that positrons follow the parent DM
distribution and have neglected positron propagation in the Galaxy that
approximately follows an energy ($E$)-dependent diffusive motion (with
$D(E)\propto E^{\delta}$, where $D$ is the diffusion coefficient
characterizing their motion and the index $\delta>0$ implies that the higher
the energy of the positrons injected, the faster they move from their
injection point). This assumption was adopted since energy losses at MeV
energies become the dominant process during positron propagation Strong and
Moskalenko (1998); however, diffusion is still very important, and leads to a
$511$ keV profile that changes as a function of the DM mass, as we demonstrate
in Fig. 1.

Remarkably, given the latest data from the AMS-02 experiment Aguilar  _et al._
(2013) on secondary cosmic rays (CRs), the diffusion of CR particles below GeV
scales seems to be faster than expected, such that we obtain better fits to
the data including a change of regime in the diffusion coefficient at low
energies Weinrich, N.  _et al._ (2020); Luque  _et al._ (2021); Luque (2021);
de la Torre Luque  _et al._ (2022). Moreover, the diffusion process is caused
by the interaction of charged particles with turbulence in the interstellar
medium, leading to energy exchange and what is usually called diffusive
reacceleration Heinbach and Simon (1995); Strong and Moskalenko (1998);
Cardillo (2019) (a kind of second order Fermi acceleration), promoting low-
energy particles to higher energies, thus making them diffuse faster. It is
important to realize that some level of reacceleration, besides being favored
by CR analyses Strong and Moskalenko (1998); Luque  _et al._ (2021), must be
present if particles are diffusing off of ISM plasma turbulence. Therefore,
instead of assuming that the 511 keV line follows the distribution of
positrons injected by DM (i.e. the DM distribution), we evaluate the 511 keV
profile from the distribution of propagated (diffuse) positrons injected by
DM. In fact, a few recent works De la Torre Luque  _et al._ (2023a, b, c),
have demonstrated the importance of adding propagation effects in determining
constraints for sub-GeV DM particles.

Here, we perform a detailed evaluation of the $511$ keV line emission from DM
and utilize the non-compatibility of these signals with the observed $511$ keV
line morphology, to constrain the DM annihilation/decay cross-section and
lifetime for the kinematically accessible channels. In doing so, we highlight
three major novelties. i) Firstly, we include the transport of positrons when
calculating the $511$ keV diffuse line emission in the Galaxy. We will in fact
demonstrate that the diffusion of positrons considerably changes their spatial
morphology (since particles of any particular energy propagate at a distinct
rate and also have energy losses that differ depending on the energy of the
particles injected by DM). ii) Secondly, we include in our evaluations the
fact that the free electron density decays very quickly as we move away from
the Galactic Plane. This results in an important suppression of the latitude
profile flux of the $511$ keV line emission. Accounting for this effect is
necessary for any search for positron sources that could produce $511$ keV
radiation and is not limited to sub-GeV DM. iii) Thirdly, we find the limits
obtained from the full SPI dataset are most constraining below DM masses of
$\sim 1$ GeV, a mass range of particular interest given low sensitivity in DM
direct detection experiments to light DM Cirelli  _et al._ (2023); De la Torre
Luque  _et al._ (2023c). Furthermore, we note that this is the first time (as
we have identified), that strong limits on the lifetime of decaying DM have
been derived with this dataset.

##  II Evaluation of the 511 keV signal

It has been widely assumed in the literature that the distribution of diffuse
positrons generated by the decay/annihilation of DM inherits the parent DM
Galactic distribution. However, this does not account for the fact that
positrons can propagate and thereby diffuse large distances from their source,
as pointed out by, e.g. Ref. Jean, P.  _et al._ (2009), thereby changing their
distribution and, eventually, the resulting profile of the observed diffuse
$511$ keV line in the Galaxy. Therefore, in this work, we evaluate the DM-
induced $511$ keV line flux considering, for the first time, the full
propagation of positrons, including all sources of energy losses and other
effects, such as reacceleration Seo and Ptuskin (1994); Osborne and Ptuskin
(1987), triplet pair production Gaggero  _et al._ (2013) and in-flight
annihilation Heitler (1936). In this way, we calculate the steady-state
solution for the distribution of positrons injected by DM (i.e. the diffuse
positrons) that will then thermalize in the ISM and, eventually, form
positronium bound states when finding free electrons in the interstellar gas
(a similar approach to that followed by Ref. Calore  _et al._ (2022)). To do
this whilst simultaneously taking into account that the spatial distribution
of diffuse positrons changes with energy, we integrate the distribution of
diffuse positrons over all possible energies.

Then, we compare our results with more than $16$ years of data collected by
INTEGRAL/SPI, which recently provided the latitudinal and longitudinal
profiles of the diffuse 511 keV line emission with unprecedented quality
Siegert  _et al._ (2019).

Positron propagation: Following a similar procedure as in Ref. De la Torre
Luque  _et al._ (2023c), we first compute the distribution and energy spectra
of the positrons produced by DM annihilation/decay in the Galaxy with a
customized version De La Torre Luque (2023) of the DRAGON2 code Evoli  _et
al._ (2017, 2018), a dedicated CR propagation code designed to self-
consistently numerically solve the full diffusion-advection-loss equation for
the propagation of charged particles in the Galactic environment Ginzburg and
Syrovatskii (1969). We use the same expression for the diffusion coefficient
and propagation parameters (the best-fit values obtained from combined fits to
AMS-02 CR data) as in Ref. De la Torre Luque  _et al._ (2023c), where we refer
the reader for more details. It is important to remark that these simulations
include all the sources of energy losses, inelastic interactions, triplet pair
production Gaggero  _et al._ (2013) and in-flight annihilation of positrons.

We simulate the kinematically accessible channels leading to the production of
$e^{\pm}$ final states through the processes:
$\chi\overline{\chi}\rightarrow\mu^{+}\mu^{-},\,\pi^{+}\pi^{-},\,e^{+}e^{-}$
for DM annihilation and
$\chi\rightarrow\mu^{+}\mu^{-},\,\pi^{+}\pi^{-},\,e^{+}e^{-}$ for DM decay.
The injection spectra $dN_{e}/dE_{e}$ for each channel are calculated
following the prescriptions given in Ref. Cirelli  _et al._ (2021). In these
simulations, we compute the signals produced by light DM from $m_{\chi}=m_{i}$
for annihilation and $m_{\chi}=2m_{i}$ for decays, up to $m_{\chi}=5$ GeV
(where $i=e,\mu,\pi$ depending on the final state channel of interest), taking
$3$ masses per decade and a NFW DM distribution Navarro  _et al._ (1996), as a
baseline.

Simulations are performed with a spatial resolution of $\sim 150$ pc, which is
sufficiently fine for comparisons with SPI data (the binning of SPI data
covers around $3\deg$ Siegert  _et al._ (2019), which implies a spatial zone
which is of several hundred parsecs) and for the smooth NFW profile that we
implement. We simulate electron-positron signals from DM annihilation and
decay in the range of kinetic energies from $100$ eV to $10$ GeV, with an
energy resolution of $5\%$. The reason it is necessary to consider such low
energies will become clear below. We have checked that lowering the minimum
energy or the spatial resolution does not impact our results significantly.
However, changing the DM profile to a more cuspy shape does.

511 keV emission:

Given the propagated and integrated distribution of diffuse positrons in the
Galaxy as a function of 3D position $(x,y,z)$ in galactic coordinates,
$\phi_{\rm e}(x,y,z)$, we calculate the emission of $511$ keV $X$-rays from a
given direction by integrating over the line of sight, as

$\frac{d\phi_{\gamma}^{511}}{d\Omega}=2k_{ps}\int
ds\,s^{2}\frac{\phi_{e}(x_{s,b,l},y_{s,b,l},z_{s,b,l})}{4\pi s^{2}}\,,$ (2)

where $k_{ps}=1/4$ is the fraction of positronium decays (corresponding to
(singlet) parapositronium states) contributing to the $511$ keV line signal,
$\phi_{\rm e}=\int\frac{d\phi_{\rm e}}{dE}dE$ is the energy-integrated flux of
positrons, $d\Omega=dldb\cos b$ is the solid angle element being $l$, $b$ and
$s$ the Galactic longitude, latitude and distance $s$ from the Earth
contributing to the signal.

As a second point of originality in our procedure, we consider that the
probability of producing positronium away from the Galactic Plane (where the
probability of encountering free electrons quickly decreases) should depend
not only on the positron density injected from DM, but also on the ambient
free electron density. Hence, we apply a scaling relation to the $511$ keV
profiles, following the vertical distribution of free electron density in the
Galaxy. This is also motivated by the fact that SPI observations of the
latitude profile seem to show a similar variation of the $511$ keV emission as
that of the free electron density (see section VI of the Appendix VI for a
more detailed discussion). We remark that this scaling leads, in fact, to a
more conservative but realistic evaluation of the predicted $511$ keV
emission. For the distribution of free electrons, we make use of the NE2001
model Cordes and Lazio (2003a, b), which contains the modelling of several
Galactic sub-structures and is derived from a combined fit of pulsar
dispersion measurements, temporal and angular broadening of radio pulses, and
intensity measurements. In this work we consider the vertical dependence of
the electron density, since the free electron density is always much higher
than the diffuse positron density in the Galactic disk (i.e. there will be
always free electrons available to form positronium states in the disk of the
Galaxy and the probability of producing these states will mainly depend on the
density of positrons). This has usually been assumed true in all the Galaxy
(see, e.g. Refs. Keith and Hooper (2021); Calore  _et al._ (2022)), however a
significantly different rate of positronium production is expected out of the
Galactic plane, because of the fast reduction of free electron density when
moving away from the Galactic plane. As we discuss below, while this scaling
has an appreciable effect in the predicted latitude profiles, it leads to
essentially no effect on the predicted longitude profiles of the $511$ keV
line emission.

In particular, we adopt the scaling predicted by the NE2001 model from the
most relevant structures on a Galactic scale, and consider a correction
proposed in Ref. Gaensler  _et al._ (2008), where the authors found that the
height of the thick disk roughly doubles to $\simeq 2$ kpc: i) a thick disk
with large Galactocentric scale height (whose scaling follows a
$\exp[-|z|/H_{1}]$ relation, where $z$ is the Galactocentric height and the
scale $H_{1}$ was found to be $1$ kpc), ii) a thin, annular disk in the inner
Galaxy (whose scaling decays as $\exp[-|z|/H_{2}]$, with $H_{2}\simeq 140$
pc), and iii) a GC component. Therefore, we take this scaling to be maximal at
the center of the Galaxy and exponentially decays like
$\exp[-|z|/H_{1}]+\exp[-|z|/H_{2}]$ at a height $z$, as shown in Fig. 3 of the
AppendixVI.

##  III Sub-GeV dark matter constraints from SPI data

Figure 2: Comparison of the $95\%$ confidence bounds on annihilating (left
panel) or decaying (right panel) DM derived in this work (thick dashed for SPI
Lat and solid lines for SPI Long, respectively) with other existing
constraints. We show the CMB bounds from Slayter Slatyer (2016) and Lopez-
Honorez et al. Lopez-Honorez  _et al._ (2013) (dot-dashed line), the bounds
from XMM Newton De la Torre Luque  _et al._ (2023c) (dotted line), the
previous 511 keV line bounds from Vincent et. al Vincent  _et al._ (2012) (red
dashed line) and from Voyager 1 De la Torre Luque  _et al._ (2023c) (dot-
dashed line). In both panels, we show the $\pi^{+}\pi^{-}$ (blue),
$\mu^{+}\mu^{-}$ (red) and $e^{+}e^{-}$ (green) channels, respectively.

Contrary to what has been assumed in the past, we show in Fig. 1 how the
predicted longitude profiles for the diffuse $511$ keV line change for various
DM masses. For concreteness, we show the longitude profile predicted from
direct annihilation to $e^{+}e^{-}$ for DM masses of $1$, $10$, $100$ and
$1000$ MeV, compared to SPI data. As we see, not only does the normalization
of the signal grow as the DM mass decreases, but also the generated profile is
peakier around the GC. This is due to the fact that the lower the energy of
the injected positrons (i.e. the lower the mass of the DM particle injecting
positrons), the smaller the distance they have to travel before thermalizing
or getting annihilated in flight (an effect that also depends on the gas
density, which is larger towards the GC). As indicated by the
$\chi^{2}$-profile ($\chi^{2}$ vs DM mass) shown in Fig. 4 of Appendix VII,
the lower the mass of the particle, the closer to the shape of the profile
observed by SPI (i.e. lower $\chi^{2}$), however, we can see that a NFW
profile is not able to generate a profile as peaky as the SPI measurements
indicate for DM masses above the electron mass. We have checked that other
typical DM density distributions do not lead to a satisfactory fit to the data
either, unless we invoke a very peaked distribution, such as the Moore profile
Moore  _et al._ (1999) (with slope $\gamma=1.5$), as illustrated in Fig. 5
(App. VIII). This fact also implies that the most constraining data points to
derive DM bounds are those away from the central longitudes. We show similar
plots for the annihilation in all the channels explored here
($\mu^{+}\mu^{-}$, $\pi^{+}\pi^{-}$ and direct $e^{+}e^{-}$) in Fig. 6 of the
Appendix IX, where we also show the predicted latitude profiles for different
DM masses. In the case of the predicted latitude profiles, we remark that: i)
Considering the scaling of the free electron density leads to a predicted
latitude profile very similar to the one measured by SPI, even for high
($\sim$ GeV) DM masses (where one would have a positron distribution much more
spread out around the Galaxy without the scaling applied), ii) the profiles
are relatively similar for different masses (although their normalization is
higher for lower DM masses) due to this scaling, and iii) It is very important
to consider the annihilation in flight of positrons, which does affect the
positron distribution above or below the Galactic plane, but it can reduce
significantly the emission in the plane (the lower the DM mass the more
important this effect becomes). In addition, as confirmed already by other
studies Boehm and Fayet (2004); Vincent  _et al._ (2012), we find that
decaying DM produces a very flat profile, compared to the observations, for
all the masses (see Fig. 7 of App. IX), what can also be seen from the high
values of $\chi^{2}$ obtained in their fit to the SPI longitude and latitude
profiles (see Fig. 4, of App. VII).

Then, we show the bounds derived from SPI data and compare them to other
existing independent DM constraints in the literature in Fig. 2, for DM
annihilation (left panel) and decay (right panel). To derive these limits, we
use a simple $\chi^{2}$ fit to the SPI data, using the curve-fit Python
package, and take the $2\sigma$ error obtained. In both cases, annihilation
and decay, we show the constraints derived in this work from the SPI latitude
profile as a dashed curve labelled “SPI Lat” and the SPI longitude profile as
a solid curve labelled “SPI Long”. We disseminate the channel of annihilation
or decay by colour based on the final state $e^{+}e^{-}$ (green),
$\mu^{+}\mu^{-}$ (red) and $\pi^{+}\pi^{-}$ (blue). In the left panel, we also
compare with the limits obtained by calculating the effect of the $e^{\pm}$
injection from DM annihilation on the CMB anisotropies, from Slatyer Slatyer
(2016) and Lopez-Honorez et al. Lopez-Honorez  _et al._ (2013), which are
shown as dot-dashed lines in Fig. 2 and are labelled “CMB.” These limits are
based on the fact that the injection of charged particles by DM during the
cosmic Dark Ages must increase the residual ionization fraction, thus altering
the anisotropies of the CMB. The CMB limits are not very strong for decaying
light DM, unlike the case of annihilating DM. This is because DM clusters at
redshifts $z\lesssim 100$, which greatly increases the DM annihilation rate
and therefore the $e^{\pm}$ injection, while the decay rate stays the same.
Since the CMB bounds are weak for decaying DM, we show the next strongest
bound which comes from Ref. De la Torre Luque  _et al._ (2023c) labelled
“Voyager-1” and shown as a dot-dashed curve in the right panel of Fig. 2. This
limit comes from the fact that decaying DM results in an $e^{\pm}$ cascade
that can be directly observed in space-based probes such as Voyager 1, without
being suppressed by the solar modulation effect of the heliosphere, due to the
satellite’s distance from the solar system.

Then, we show the strong bounds from Refs. De la Torre Luque  _et al._
(2023c), labelled as “XMM-Newton” in both panels, shown as dotted curves.
These limits are derived from the fact that $e^{\pm}$ products resulting from
an exotic DM injection (through annihilation or decay) may produce
bremsstrahlung radiation and upscatter the low-energy Galactic photon fields
via the inverse Compton process generating a broad emission from $X$-ray to
$\gamma$-ray energies observable in experiments such as Xmm-Newton which have
measured $X$-rays in galactocentric rings around the GC. These limits include
best-fit CR propagation and diffusion parameters and represented (till this
analysis) the strongest astrophysical constraints for this mass range of DM of
$1$ MeV to a few GeV and even surpass cosmological bounds across a wide range
of masses as well.

In addition, we show the bounds derived by Ref. Vincent  _et al._ (2012), also
using $511$ keV data from SPI after a shorter period ($\sim 8$ yr) of
operation. Our limits are stronger because the previous constraint used the
spectrum of the emission observed around the GC at 511 keV, instead of the
full profile. This, added to the fact that we find that the the most
constraining observations are those from high longitudes, explains the
strength of our constraints in comparison to the previous limit.

In fact, we see that the SPI longitude profile provides the strongest limit on
annihilating DM, excluding cross-sections of $\langle\sigma v\rangle\simeq
10^{-32}\,\textrm{cm}^{3}\textrm{s}^{-1}$ at masses of around an MeV,
increasing to around the thermal relic values of $\langle\sigma v\rangle\simeq
10^{-26}\,\textrm{cm}^{3}\textrm{s}^{-1}$ at masses over several GeV. This SPI
longitude profile only gets surpassed by Xmm-Newton $X$-ray bounds at masses
$\gtrsim 700$ MeV. The CMB bounds only become more stringent than the SPI
latitude profile limits at masses in excess of a few tens of GeV. Strikingly,
limits derived from the SPI latitude profile are only marginally weaker (less
than an order of magnitude across the entire DM mass range) than the SPI
longitude profile limits. This was also found in Ref. Calore  _et al._ (2022).
Furthermore, this bound also surpasses the other DM limits across most of the
mass range considered, except XMM-Newton at masses over a few hundreds of MeV
and CMB bounds at slightly higher mass. Finally, we also compare our limits
with those derived from Ref. Boehm  _et al._ (2004b) from FSR emission.

From the right panel of Fig. 2, we see once again that the SPI longitude
profile provides the strongest limit on decaying DM, excluding lifetimes of
$\tau\simeq 10^{29}\,\textrm{s}$ at masses of around an MeV decreasing to
around values of $10^{27}\,\textrm{s}$ at masses of over a GeV. This SPI
longitude profile gets surpassed by Xmm-Newton $X$-ray bounds at masses
$\gtrsim 1$ GeV. Limits derived from the SPI latitude profile are about an
order of magnitude weaker than the SPI longitude profile limits across the
entire mass range. The next strongest limit, coming from Voyager 1 is
significantly weaker (about three orders weaker than the SPI longitude profile
limits). There are also cosmological bounds for decaying DM, such as in Refs.
Liu  _et al._ (2023a, b), but we omit these since they are not competitive.
This illustrates the power of the $511$ keV line in constraining sub-GeV DM.
We remark that including the free electron density scaling changes our
latitude bounds by a factor of $2$-$3$, leading to a more conservative
estimate.

##  IV Discussion and conclusion

In this work, we have explored the possibility of sub-GeV dark matter (DM)
particles that annihilate or decay into standard model (SM) particles,
resulting in an electron-positron cascade that eventually forms positronium
bound states after thermalization. Approximately one-fourth of these bound
states decay into photons with an energy of 511 keV, which can be detected by
the SPI spectrometer on board the INTEGRAL satellite. We have used the SPI
longitude and latitude profiles to derive new limits on the DM properties,
such as the mass, the annihilation cross-section, and the decay lifetime. We
have improved previous analyses by taking into account two important effects
that have not hitherto been considered: the diffusion and propagation of
positrons in the interstellar medium, and the reduction of the free electron
density away from the Galactic plane. These effects have a significant impact
on the shape and intensity of the 511 keV photon line, and therefore on the
resulting DM constraints. Additionally, we have constrained decaying DM with
this SPI dataset for the first time. We have found that our limits are the
strongest on sub-GeV DM across almost the entire DM mass range considered
(from MeV to a few GeV), ruling out thermally-averaged cross-sections of
around $10^{-32}$ cm3s${}^{-1}\lesssim\langle\sigma v\rangle\lesssim 10^{-26}$
cm3s-1 and decay lifetimes of $10^{29}\,\textrm{s}\lesssim\tau\lesssim
10^{27}\,\textrm{s}$ respectively) by including best-fit cosmic ray (CR)
propagation and diffusion parameters. Our limits are among the strongest and
most robust DM constraints to date, surpassing the most stringent
complementary cosmological and astrophysical limits over majority of the mass
range considered.

However, it is important to comment on the uncertainties affecting our
results. First, the total systematic uncertainties of the SPI data points are
not totally under control, since they depend on some template models for the
shape of the signal. This could result in extra systematic uncertainties
($\mathcal{O}(10\%)$) not included in the public data. We would expect
increasing the systematic uncertainty in our fits to weaken our limits at most
by a factor of two, leaving our main conclusions unchanged. Currently, the
details of the propagation of CRs beyond the our local environment are not
known. There are recent indications that diffusion is not uniform throughout
the Galaxy Amenomori  _et al._ (2021); Cao  _et al._ (2023); Collaboration
(2023); De La Torre Luque  _et al._ (2023), and theoretically there is nearly
no doubt that one should expect inhomogeneous (and probably anisotropic)
propagation of charged particles in the Galaxy Schroer  _et al._ (2022);
Reichherzer  _et al._ (2022); Cerri  _et al._ (2017). However, we only have
access to CR data in the Solar System, which is limited to the GeV to TeV
range Balaji (2023). This also implies that our diffusion models are still
quite uncertain at sub-GeV energies, which forces us to use extrapolations of
our best-fit (spatially uniform) models above the GeV in this study. On top of
this, strong convection or winds near the GC could also significantly change
the spectrum and profiles of positrons and electrons. However, a correct
estimation of such effects is far from being possible with the current data. A
study of the impact of the uncertainties associated with the propagation
parameters is beyond the scope of this paper, but will be addressed in a
future work. An evaluation of the diffuse $511$ keV line emission from
Galactic positron sources will help strengthen the robustness of these limits
and deepen our understanding of this mysterious emission. But what this work
clearly demonstrates is that diffusion should be taken into account in our
future studies, since it plays a crucial role in shaping the spectrum and
profile of this line. Data from the future COSI satellite Siegert  _et al._
(2020), with better spatial resolution, would also provide a better insight of
the processes and structures producing the $511$ keV line in the Milky Way.

In addition, possible tests of this class of DM scenarios, by characterizing
the diffuse MeV-scale $\gamma$-ray emission from the halo of the Milky Way,
with proposed $\gamma$-ray telescopes such as AMEGO McEnery  _et al._ (2019)
or e-ASTROGAM De Angelis  _et al._ (2018) would provide useful background
models. These could be used to further strengthen bounds on sub-GeV DM.

##  V Acknowledgements

We are grateful to Pierluca Carenza, Marco Cirelli, Daniele Gaggero, Dan
Hooper, Jordan Koechler, Tim Linden, Pierrick Martin and Thomas Siegert for
useful discussions. SB is supported by the STFC under grant ST/X000753/1. PDL
is supported by the European Research Council under grant 742104 and the
Swedish National Space Agency under contract 117/19. This project used
computing resources from the Swedish National Infrastructure for Computing
(SNIC) under project Nos. 2021/3-42, 2021/6-326, 2021-1-24 and 2022/3-27
partially funded by the Swedish Research Council through grant no. 2018-05973.

## References

  * Siegert  _et al._ (2016) T. Siegert _et al._ , Astron. Astrophys. 586, A84 (2016), arXiv:1512.00325 [astro-ph.HE] .
  * Kierans  _et al._ (2020) C. A. Kierans _et al._ , Astrophys. J. 895, 44 (2020), arXiv:1912.00110 [astro-ph.HE] .
  * Prantzos  _et al._ (2011) N. Prantzos _et al._ , Rev. Mod. Phys. 83, 1001 (2011), arXiv:1009.4620 [astro-ph.HE] .
  * Kierans  _et al._ (2019) C. A. Kierans _et al._ , “Positron annihilation in the galaxy,” (2019), arXiv:1903.05569 [astro-ph.HE] .
  * Siegert (2023) T. Siegert, Astrophysics and Space Science 368 (2023), 10.1007/s10509-023-04184-4.
  * Beacom and Yuksel (2006) J. F. Beacom and H. Yuksel, Phys. Rev. Lett. 97, 071102 (2006), arXiv:astro-ph/0512411 .
  * Boehm (2009) C. Boehm, New Journal of Physics 11, 105009 (2009).
  * Sizun  _et al._ (2007) P. Sizun, M. Cassé, S. Schanne, and B. Cordier, “Constraints on the injection energy of positrons in the galactic centre region,” (2007), arXiv:astro-ph/0702061 [astro-ph] .
  * Bartels  _et al._ (2018) R. Bartels _et al._ , Mon. Not. Roy. Astron. Soc. 480, 3826 (2018), arXiv:1803.04370 [astro-ph.HE] .
  * Fuller  _et al._ (2019) G. M. Fuller _et al._ , Phys. Rev. Lett. 122, 121101 (2019), arXiv:1811.00133 [astro-ph.HE] .
  * Boehm  _et al._ (2004a) C. Boehm _et al._ , Phys. Rev. Lett. 92, 101301 (2004a), arXiv:astro-ph/0309686 .
  * Vincent  _et al._ (2012) A. C. Vincent, P. Martin, and J. M. Cline, JCAP 04, 022 (2012), arXiv:1201.0997 [hep-ph] .
  * Ascasibar  _et al._ (2006) Y. Ascasibar _et al._ , Mon. Not. Roy. Astron. Soc. 368, 1695 (2006), arXiv:astro-ph/0507142 .
  * Beacom  _et al._ (2005) J. F. Beacom, N. F. Bell, and G. Bertone, Phys. Rev. Lett. 94, 171301 (2005), arXiv:astro-ph/0409403 .
  * Boehm and Uwer (2006) C. Boehm and P. Uwer, “Revisiting bremsstrahlung emission associated with light dark matter annihilations,” (2006), arXiv:hep-ph/0606058 [hep-ph] .
  * Escudero (2019) M. Escudero, JCAP 02, 007 (2019), arXiv:1812.05605 .
  * Sabti  _et al._ (2020) N. Sabti _et al._ , JCAP 01, 004 (2020), arXiv:1910.01649 .
  * Strong and Moskalenko (1998) A. W. Strong and I. V. Moskalenko, The Astrophysical Journal 509, 212 (1998).
  * Aguilar  _et al._ (2013) M. Aguilar _et al._ (AMS Collaboration), Phys. Rev. Lett. 110, 141102 (2013).
  * Weinrich, N.  _et al._ (2020) Weinrich, N. _et al._ , A&A 639, A131 (2020).
  * Luque  _et al._ (2021) P. D. L. T. Luque _et al._ , JCAP 07, 010 (2021), arXiv:2102.13238 [astro-ph.HE] .
  * Luque (2021) P. D. L. T. Luque, JCAP 11, 018 (2021), arXiv:2107.06863 .
  * de la Torre Luque  _et al._ (2022) P. de la Torre Luque _et al._ , JCAP 07, 008 (2022), arXiv:2202.03559 [astro-ph.HE] .
  * Heinbach and Simon (1995) U. Heinbach and M. Simon, Astrophys. J. 441, 209 (1995).
  * Cardillo (2019) M. Cardillo, Galaxies 7, 49 (2019), arXiv:1904.12476 .
  * De la Torre Luque  _et al._ (2023a) P. De la Torre Luque, S. Balaji, and P. Carenza, (2023a), arXiv:2307.13731 [hep-ph] .
  * De la Torre Luque  _et al._ (2023b) P. De la Torre Luque, S. Balaji, and P. Carenza, (2023b), arXiv:2307.13728 [hep-ph] .
  * De la Torre Luque  _et al._ (2023c) P. De la Torre Luque, S. Balaji, and J. Koechler, (2023c), arXiv:2311.04979 [hep-ph] .
  * Cirelli  _et al._ (2023) M. Cirelli _et al._ , JCAP 07, 026 (2023), arXiv:2303.08854 .
  * Jean, P.  _et al._ (2009) Jean, P. _et al._ , A&A 508, 1099 (2009).
  * Seo and Ptuskin (1994) E.-S. Seo and V. S. Ptuskin, ApJ 431, 705 (1994).
  * Osborne and Ptuskin (1987) J. Osborne and V. Ptuskin, in  _International Cosmic Ray Conference_ , Vol. 2 (1987) p. 218.
  * Gaggero  _et al._ (2013) D. Gaggero, P. D. Serpico, and V. Bonnivard, in  _33rd International Cosmic Ray Conference_ (2013) p. 0361.
  * Heitler (1936) W. Heitler,  _The quantum theory of radiation_ , International Series of Monographs on Physics, Vol. 5 (Oxford University Press, Oxford, 1936).
  * Calore  _et al._ (2022) F. Calore _et al._ , Phys. Rev. D 105, 063026 (2022), arXiv:2112.08382 [hep-ph] .
  * Siegert  _et al._ (2019) T. Siegert _et al._ , A& A 627, A126 (2019).
  * De La Torre Luque (2023) P. De La Torre Luque, “Dragon2_optimized_dm&antinuclei,”  (2023).
  * Evoli  _et al._ (2017) C. Evoli _et al._ , JCAP 02, 015 (2017), arXiv:1607.07886 [astro-ph.HE] .
  * Evoli  _et al._ (2018) C. Evoli _et al._ , JCAP 07, 006 (2018), arXiv:1711.09616 [astro-ph.HE] .
  * Ginzburg and Syrovatskii (1969) V. L. Ginzburg and S. I. Syrovatskii,  _Topics in Astrophysics and Space Physics, New York: Gordon and Breach, 1969_ (Gordon & Breach Publishing Group, 1969).
  * Cirelli  _et al._ (2021) M. Cirelli _et al._ , Phys. Rev. D 103, 063022 (2021), arXiv:2007.11493 [hep-ph] .
  * Navarro  _et al._ (1996) J. F. Navarro, C. S. Frenk, and S. D. White, Astrophys. J. 462, 563 (1996), arXiv:astro-ph/9508025 .
  * Cordes and Lazio (2003a) J. M. Cordes and T. J. W. Lazio, “Ne2001. ii. using radio propagation data to construct a model for the galactic distribution of free electrons,” (2003a), arXiv:astro-ph/0301598 [astro-ph] .
  * Cordes and Lazio (2003b) J. M. Cordes and T. J. W. Lazio, “Ne2001.i. a new model for the galactic distribution of free electrons and its fluctuations,” (2003b), arXiv:astro-ph/0207156 [astro-ph] .
  * Keith and Hooper (2021) C. Keith and D. Hooper, Phys. Rev. D 104, 063033 (2021).
  * Gaensler  _et al._ (2008) B. M. Gaensler _et al._ , Publications of the Astronomical Society of Australia 25, 184–200 (2008).
  * Slatyer (2016) T. R. Slatyer, Phys. Rev. D 93, 023527 (2016), arXiv:1506.03811 [hep-ph] .
  * Lopez-Honorez  _et al._ (2013) L. Lopez-Honorez _et al._ , JCAP 07, 046 (2013), arXiv:1303.5094 [astro-ph.CO] .
  * Moore  _et al._ (1999) B. Moore _et al._ , The Astrophysical Journal 524, L19–L22 (1999).
  * Boehm and Fayet (2004) C. Boehm and P. Fayet, Nucl. Phys. B 683, 219 (2004), arXiv:hep-ph/0305261 .
  * Boehm  _et al._ (2004b) C. Boehm, T. A. Enßlin, and J. Silk, Journal of Physics G: Nuclear and Particle Physics 30, 279–285 (2004b).
  * Liu  _et al._ (2023a) H. Liu _et al._ , Phys. Rev. D 108, 043530 (2023a), arXiv:2303.07366 [astro-ph.CO] .
  * Liu  _et al._ (2023b) H. Liu _et al._ , Phys. Rev. D 108, 043531 (2023b), arXiv:2303.07370 [astro-ph.CO] .
  * Amenomori  _et al._ (2021) M. Amenomori _et al._ (Tibet ${\mathrm{AS}}_{\gamma}$ Collaboration), Phys. Rev. Lett. 126, 141101 (2021).
  * Cao  _et al._ (2023) Z. Cao _et al._ , Physical Review Letters 131 (2023), 10.1103/physrevlett.131.151001.
  * Collaboration (2023) I. Collaboration, Science 380, 1338 (2023), https://www.science.org/doi/pdf/10.1126/science.adc9818 .
  * De La Torre Luque  _et al._ (2023) P. De La Torre Luque _et al._ , A&A 672, A58 (2023).
  * Schroer  _et al._ (2022) B. Schroer _et al._ , MNRAS 512, 233 (2022), https://academic.oup.com/mnras/article-pdf/512/1/233/42894576/stac466.pdf .
  * Reichherzer  _et al._ (2022) P. Reichherzer _et al._ , MNRAS 514, 2658–2666 (2022).
  * Cerri  _et al._ (2017) S. S. Cerri _et al._ , Journal of Cosmology and Astroparticle Physics 2017, 019–019 (2017).
  * Balaji (2023) S. Balaji, Phys. Lett. B 845, 138157 (2023), arXiv:2211.03807 [astro-ph.HE] .
  * Siegert  _et al._ (2020) T. Siegert _et al._ , The Astrophysical Journal 897, 45 (2020).
  * McEnery  _et al._ (2019) J. McEnery _et al._ , “All-sky medium energy gamma-ray observatory: Exploring the extreme multimessenger universe,” (2019), arXiv:1907.07558 [astro-ph.IM] .
  * De Angelis  _et al._ (2018) A. De Angelis _et al._ , Journal of High Energy Astrophysics 19, 1–106 (2018).
  * Ackermann  _et al._ (2017) M. Ackermann _et al._ , The Astrophysical Journal 840, 43 (2017).
  * Gondolo and Silk (1999) P. Gondolo and J. Silk, Physical Review Letters 83, 1719–1722 (1999).
  * Lacroix (2018) T. Lacroix, Astron. Astrophys. 619, A46 (2018), arXiv:1801.01308 [astro-ph.GA] .
  * Balaji  _et al._ (2023) S. Balaji _et al._ , JCAP 08, 063 (2023), arXiv:2303.12107 [hep-ph] .
  * Spekkens and Giovanelli (2005) K. Spekkens and R. Giovanelli, Astron. J. 129, 2119 (2005), arXiv:astro-ph/0502166 .
  * Walker  _et al._ (2012) S. A. Walker _et al._ , Mon. Not. Roy. Astron. Soc. 424, 1826 (2012), arXiv:1205.2276 [astro-ph.CO] .
  * Foëx  _et al._ (2012) G. Foëx _et al._ , Astronomy & Astrophysics 546, A106 (2012).
  * Vegetti and Vogelsberger (2014) S. Vegetti and M. Vogelsberger, Mon. Not. Roy. Astron. Soc. 442, 3598 (2014), arXiv:1406.1170 [astro-ph.CO] .
  * Navarro  _et al._ (2004) J. F. Navarro _et al._ , Mon. Not. Roy. Astron. Soc. 349, 1039 (2004), arXiv:astro-ph/0311231 .

##  VI Height scale of the free electron density

As explained above, in order to consider that different regions of the Galaxy
will have a different rate (probability) of positronium formation, that also
depends on the distribution of the ambient (free) electrons, we consider a
scaling of the $511$ keV line emission following the NE2001 model Cordes and
Lazio (2003a, b) for the free electron density distribution. In particular, we
consider that the disk of the Galaxy will always contain a much larger density
of free electrons in the ISM gas than diffuse positrons (i.e. there will be
always more electrons available than positrons to form positronium states,
thus being “saturated”, and the probability of positronium formation will
mainly depend on the diffuse positron distribution). However, this assumption
must be broken well above and below the disk, since the density of free
electrons falls very quickly as the galactic height, $z$, is increased, as
shown by Fig. 3. In particular, in this figure we depict the electron density
distribution, at the center of the Galaxy, as a function of the Galactic
height, normalized to the electron density at the GC position. This is the
scale factor that we apply to our evaluations of the $511$ keV line. This
scaling is particularly important for very low DM masses, when the positrons
lose energy very fast in the disk (but not outside the disk, since the density
of gas producing these losses is much lower) and the fraction of positrons
annihilated in flight (directly proportional to the gas density where the
positrons propagate) becomes relevant.

This is also motivated by the fact that, as we observe from the SPI
measurements of the latitude profile, the intensity of the $511$ keV line
decays by around one order of magnitude at $10^{\circ}$ above or below the
disk. Assuming the the height of the halo is around $4$-$10$ kpc, this means
that with every $\sim 0.36$-$0.9$ kpc, the intensity of the line decays by one
order of magnitude, similar to the factor $\gtrsim 6$ that the free electron
density decays every $\sim$kpc above or below the disk. While the predicted
latitude profiles from DM appear to be quite flat (compared to the SPI data)
for all the masses explored here, using this scaling we find that the
predicted latitude profiles become much more similar to the shape of the SPI
measurements. However, we remark that the longitude profiles predicted are
roughly unaffected by this scaling.

Figure 3: Vertical profile of the free electron density normalized at the
Galactic Center as a function of height above the Galactic Plane $z$. We take
this distribution from the NE2001 model Cordes and Lazio (2003a, b), as
indicated in the text.

##  VII $\chi^{2}$ vs dark matter mass from the fits of SPI data

In this appendix, we show the $\chi^{2}$ values obtained from the fits of the
predicted $511$ keV latitudinal and longitudinal profiles to the SPI
measurements ($\sim 16$ yr), as a function of the DM mass. The $\chi^{2}$
value denotes how similar the predicted profile (for the best-fit annihilation
rate) is to the one observed by SPI. Fig. 4 shows the $\chi^{2}$ dependence
with the DM mass for decay and annihilation, for the longitude (left panel)
and latitude (right panel) profiles. First, we remark that DM decay does not
lead to any good fit of the SPI measurements at all (except at high masses for
the latitude profile). In turn, DM annihilation leads to profiles of the $511$
keV emission that more closely resemble the data. However, while the longitude
profiles favors low DM masses, the latitude profiles start to worsen below
$\sim 10$ MeV. This misalignment in preferred mass between the latitude and
longitude profile may indicate that DM is not a favoured explanation for the
observed 511 keV line emission. We remind the reader that these results
strongly depend on the DM density distribution, which as we mentioned earlier,
we set to follow a NFW distribution, with density at Solar system position of
$0.4$ GeV/cm3.

Figure 4: Reduced $\chi^{2}$ ($\chi^{2}$ over number of SPI data points) as
function of DM mass $m_{\chi}$ for annihilation and decay. The left panel
shows the values obtained for the analysis of the SPI longitude profile and
the right panel the results for the SPI latitude profile. The channels are
shown in the same colour as Fig. 2.

##  VIII 511 keV emission profiles from common Galactic dark matter density
distributions

In this appendix, we report a comparison of the longitude (left panels of Fig.
5) and latitude (right panels of Fig. 5) profiles predicted for DM
annihilation, assuming different standard Galactic DM distributions, for DM
masses of $10$ and $100$ MeV. In particular, we show the predictions obtained
for the NFW profile (cyan lines), a contracted NFW profile (c-NFW), like the
one that provides the best-fit to the Galactic Center Excess (GCE) Ackermann
_et al._ (2017), which has a slope of $\gamma=1.25$ (green lines), a Moore
profile Moore  _et al._ (1999) , which has a slope of $\gamma=1.5$ (blue
lines), and a cuspy profile, as an example of a distribution that accounts for
spikes formed around the central black hole Sgr A*, taken from Refs. Gondolo
and Silk (1999); Lacroix (2018); Balaji  _et al._ (2023) (black lines).

We defer the full study of the compatibility of the different DM distributions
with the SPI observations to a future work. However, this comparison allows us
to see how the profiles are sensitive to the different choice of the Galactic
DM distribution. In particular, we see that profiles that are not very steep
(with slope $\gamma\lesssim 1.25$) are not able to produce a good simultaneous
fit to the longitudinal and latitudinal SPI profiles for DM mass exceeding an
MeV. However, a steeper profile would be able to provide relatively
satisfactory fits for MeV masses. However, we remark that a high slope, like
that adopted in the Moore profile, is already in tension with some
astrophysical observations. This is because it can predict too much DM in the
inner regions of galaxies and galaxy clusters to comply with rotation curve
observations Spekkens and Giovanelli (2005), $X$-ray observations of galactic
clusters Walker  _et al._ (2012), gravitational lensing observations Foëx  _et
al._ (2012); Vegetti and Vogelsberger (2014) and N-body simulations of
$\Lambda$CDM halos Navarro  _et al._ (2004). In general, this leads to the
conclusion that DM cannot be the dominant source of positrons in the Galaxy,
but, instead, it is expected that a combination of different sources is
producing the peaked $511$ keV emission profiles that SPI observes.

Figure 5: Comparison of the predicted diffuse $511$ keV line emission from a
low-mass DM particle of mass $m_{\chi}=10$ MeV (top panels) and $m_{\chi}=100$
MeV (bottom panels) annihilating into electron-positron pairs for different DM
density profiles.The left panels are the predicted longitude profiles, for
latitude bin $-10.75^{\circ}<b<10.25^{\circ}$ while the right panels are the
latitude profile, for a longitude bin of $-5.25^{\circ}<l<3.75^{\circ}$. We
fix the annihilation rate to the values shown in the plots and scale the
profiles for different masses to facilitate the comparison with the SPI data.
We show the NFW profile (cyan), contracted NFW (green), Moore (blue) and cuspy
DM profiles.

##  IX Annihilation and decay profiles for different channels

Figure 6: Comparison of the predicted longitude (left panels) and latitude
(right panels) of the $511$ keV emission from DM annihilation for different DM
masses, through the $\mu^{+}\mu^{-}$ (top panels), $\pi^{+}\pi^{-}$ (middle
panels) and direct $e^{+}e^{-}$ (bottom panels) channels. The lowest mass for
the case of the $\pi^{+}\pi^{-}$ annihilation channel is actually
$m_{\chi}=2m_{\pi}\simeq 139$ MeV and the minimum energy shown in the
$\mu^{+}\mu^{-}$ is $m_{\chi}=2m_{\mu}\simeq 105$ MeV. We fix the annihilation
rate to the values shown in the plots and scale the profiles for different
masses to facilitate the comparison with the SPI data. The DM profiles are
shown in the same colours and the left and right panels are the same latitude
bins as Fig. 5.

In this appendix, we report the predicted latitude and longitude profiles of
the $511$ keV emission from DM decay and annihilation. In particular, we show
in Fig. 6 a comparison of profiles obtained assuming DM annihilation through
the $\mu^{+}\mu^{-}$ (top panels), $\pi^{+}\pi^{-}$ (middle panelbos) and
direct $e^{+}e^{-}$ (bottom panels) channels for different DM masses (from
twice the mass of the final state, around $100$ MeV, to $5$ GeV, in the case
of the $\pi^{+}\pi^{-}$ and $\mu^{+}\mu^{-}$ channels, and from $1$ MeV to $1$
GeV for the direct $e^{+}e^{-}$ channel). The left panels illustrate the
predicted longitude profiles of the $511$ keV emission and the right panels
the predicted latitude profiles. These comparisons allow us to see that the
data points corresponding to the higher longitude and latitude angles (in
absolute value) are the most constraining ones, which explains why the bounds
derived from this dataset are much stronger than those derived from Ref.
Vincent  _et al._ (2012) (besides the fact that these measurements are
reported after double the duration of event collection). As we see, the
predicted longitude profiles show a clear trend: the lower the mass of the DM
particle, the more concentrated the $511$ keV emission around the GC. This is
expected, since a lower DM mass means a lower injection energy for the
positrons, and the lower the energy of the positrons, the slower they diffuse.
Hence lower energy positrons have a distribution more similar to the DM
distribution (NFW in these calculations) and lose their energy closer to their
injection point. As we pointed out in the previous section, the predictions
from a NFW DM distribution do not lead to a $511$ keV emission longitude
profile as peaky as that observed by the SPI data even for the low DM masses
that are accessible in the direct $e^{+}e^{-}$ channel. In turn, as we see
from the predicted latitude profiles, masses of around $100$ MeV reproduce the
observed SPI latitude profile quite well. We remark that this is due to the
scaling that accounts for the decay of the density of free electrons away from
the GC, applied for first time in this work. Without accounting for this
effect, we observe that the latitude profiles become much more diffuse (i.e.
flatter profiles) than what is observed. In addition, for the latitude
profiles, we see that the same trend wherein lower DM mass leads to a profile
more concentrated around the Galactic Plane is observed above some tens of
MeV. However, once the energy loss (mainly Coulomb and ionization energy
losses) becomes the dominant process during the positron propagation, the
$511$ keV emission at the Galactic plane starts to decay quickly, given that
these particles suffer from severe energy losses and annihilation in flight
losses (since this effect is proportional to the gas density, it does not
affect the positrons above and below the Galatic plane). We have investigated
that the minimum energy employed in these simulations ($100$ eV) is low enough
that our results do not change appreciably for a minimum energy of $1$ eV.

Finally, we show in Fig. 7 the predicted profiles from decay of DM for masses
of $1$ MeV to $1$ GeV, in the direct $e^{+}e^{-}$ channel. As already observed
in the past, and discussed above, we observed that both the latitude and
longitude profiles of the predicted $511$ keV emission are much flatter (for
every mass) than the ones measured by SPI, disfavouring DM decay as the main
producer of the emission around the GC. This is simply due to the fact that
the total positron emission from decay is proportional to the DM density,
while in the case of annihilation the emission is $\propto\rho_{DM}^{2}$.

Figure 7: Similar to Fig. 6 but for the case of DM decay into the direct
$e^{+}e^{-}$ channel.


# Proximity induced superconductivity in a topological insulator

Philipp Rüßmann<EMAIL_ADDRESS>Institute of Theoretical Physics
and Astrophysics, University of D-97074 Würzburg, Germany  Peter Grünberg
Institut and Institute for Advanced Simulation, Forschungszentrum Jülich and
JARA, D-52425 Jülich, Germany  Stefan Blügel  Peter Grünberg Institut and
Institute for Advanced Simulation, Forschungszentrum Jülich and JARA, D-52425
Jülich, Germany

###### Abstract

Interfacing a topological insulator (TI) with an $s$-wave superconductor (SC)
is a promising material platform that offers the possibility to realize a
topological superconductor through which Majorana-based topologically
protected qubits can be engineered. In our computational study of the
prototypical SC/TI interface between Nb and Bi2Te3, we identify the benefits
and possible bottlenecks of this potential Majorana material platform.
Bringing Nb in contact with the TI film induces charge doping from the SC to
the TI, which shifts the Fermi level into the TI conduction band. For thick TI
films, this results in band bending leading to the population of trivial TI
quantum-well states at the interface. In the superconducting state, we uncover
that the topological surface state experiences a sizable superconducting gap-
opening at the SC/TI interface, which is furthermore robust against
fluctuations of the Fermi energy. We also show that the trivial interface
state is only marginally proximitized, potentially obstructing the realization
of Majorana-based qubits in this material platform.

##  I Introduction

The promise to realize topologically protected qubits based on the non-Abelian
exchange statistics of Majorana zero modes that can be nucleated in
topological superconductors sparked a lot of research interest in the field in
the past years Fu and Kane (2008); Cook and Franz (2011); Alicea  _et al._
(2011); Sato and Ando (2017); Flensberg  _et al._ (2021). Despite the ongoing
effort in increasing, for example, the decoherence times of qubits, the entire
field of quantum computing still faces material optimization challenges of the
different pursued qubit architectures and Majorana-based platforms are no
exception de Leon  _et al._ (2021).

TIs are a very versatile material and recent advances in fabrication
techniques open up new ways to use their unique properties for novel
functionalities in devices Hasan and Kane (2010); Breunig and Ando (2022). For
example, recently a TI-based Josephson junction was fabricated that
demonstrated the usability of TIs in the field of superconducting transmon
qubits Schmitt  _et al._ (2022). The fundamental property of the
superconducting proximity effect in SC/TI heterostructures has been a major
research focus in the past years, which was pursued both experimentally Xu
_et al._ (2015, 2014a); Stolyarov  _et al._ (2021); Flötotto  _et al._ (2018);
Xu  _et al._ (2014b); Bai  _et al._ (2020); Schmitt  _et al._ (2022) and
theoretically Fu and Kane (2008); Chiu  _et al._ (2016); Park  _et al._
(2020); Legg  _et al._ (2022). For instance, using scanning tunneling
microscopy Xu  _et al._ (2015, 2014a); Stolyarov  _et al._ (2021), angle-
resolved photoemission Flötotto  _et al._ (2018); Xu  _et al._ (2014b), or
transport experiments Bai  _et al._ (2020), it was demonstrated that TIs of
the (Bi,Sb)2(Te,Se)3 family can indeed be proximitized by an interface to a
superconductor. Elemental $s$-wave superconductors like Nb Schmitt  _et al._
(2022); Flötotto  _et al._ (2018) or Pb Stolyarov  _et al._ (2021) but also
two-dimensional van der Waals superconductors like NbSe2 Xu  _et al._ (2015,
2014a, 2014b) or PdTe2 Park  _et al._ (2020); Bai  _et al._ (2020) have been
used successfully to open a superconducting gap in the electronic structure of
a TI. It was demonstrated that the proximity effect decays exponentially in
the topological surface state (TSS) with its distance to the SC/TI interface
(i.e. with the TI film thickness) Flötotto  _et al._ (2018); Xu  _et al._
(2014a); Park  _et al._ (2020). Further analysis even revealed possible
signatures of Majorana zero modes in form of zero-bias voltage peaks inside
supercurrent vortices of the SC/TI heterostructures Xu  _et al._ (2015,
2014a).

Comparing different SC/TI interfaces, it is clear that the details of the
interface chemistry strongly influence the proximity effect. This is reflected
(i) in the different extent to which the quantum-well states and the TSS of
thin TI films are gapped out upon contact with a SC and (ii) in how quickly
the proximity-induced gap in the different states decay from the SC/TI
interface Xu  _et al._ (2014b); Flötotto  _et al._ (2018); Park  _et al._
(2020). These differences in the proximity effect were attributed to the
degree of hybridization of the respective TI wave functions with the
superconductor and the resulting possibility for Cooper pair tunneling into
the TI Xu  _et al._ (2014b); Park  _et al._ (2020). Only in the recent past,
it was realized that at the SC/TI contact charge transfer and band bending
effects play an important role for the proximity effect and the opening of a
topological gap Stolyarov  _et al._ (2021); Legg  _et al._ (2022). It was
understood that the effect of charge transfer can even be beneficial in
proximitized TI nanowires due to the induced inversion symmetry breaking, as
it can enhance the sub-band splitting of TI states required for the
realization of topological superconductivity in this geometry Legg  _et al._
(2022).

Despite the numerous research on the proximity effect in SC/TI
heterostructures, the interplay of the microscopic details of the interface,
the local electronic structure, and the factors controlling the proximity
effect in this material platform remains elusive, and to date no topological
qubit exists. Thus, to fill this gap, a deeper understanding of the physics of
the SC/TI interface is obligatory, which we aim at in this study.

The focus of this computational study is the induced proximity effect in the
electronic structure of the TI Bi2Te3 in contact to the $s$-wave
superconductor niobium (Nb). For the TI we consider thin films of Bi2Te3 of
different thicknesses between two and ten quintuple layers (QLs), which
corresponds to thicknesses of $\approx 2-10$ nm. The SC contact is modelled
using six layers of Nb(111) in order to build a minimal structural model for
the SC/TI interface. Our calculations are based on density functional theory
(DFT) and the Bogoliubov-de Gennes (BdG) method to include the effect of
superconductivity together with a realistic description of the electronic
structure Rüßmann and Blügel (2022a). Further calculational details can be
found in the “Methods” section.

The calculational setup is shown in Fig. 1(a) where the “free” (i.e. the QL
adjacent to the vacuum region) and the “contact” QL (i.e. the QL in contact to
the superconductor) are highlighted. Throughout this work we discuss the
electronic structure on both free and contact sides of the TI film. In
general, in the TI/SC heterostructure a TSS will appear on the free side at
the interface between TI and vacuum, and also as a contact TSS at the
interface between the TI and the SC which can hybridize with the states from
the SC. Chemically, the SC/TI interface is a metal-semiconductor contact where
charge transfer and band bending effects can lead to the occupation of trivial
interface states in the accumulation region Bahramy  _et al._ (2012);
Michiardi  _et al._ (2022). Our work discusses these different states that
appear at the SC/TI contact and uncovers the interplay of wave function
localization, superconducting proximity effect which may allow to further
optimize the SC/TI interface using elemental SCs in the future.

##  II Results

###  II.1 Normal state properties of the Bi2Te3/Nb interface

Figure 1(b) shows an exemplary normal state density of states (DOS) without
superconductivity for a three QL thick Bi2Te3 film in contact to Nb. The DOS
is integrated over the regions of the free QL, the contact QL and the Nb
region. The DOS in the free QL hardly changes compare to the DOS in
freestanding Bi2Te3 (not shown). In contrast, the QL in contact with Nb
clearly hybridizes with the states originating from the Nb film, introducing
additional states that are particularly visible in the region of the TI bulk
band-gap around $E-E_{\mathrm{F}}=-0.3\mathrm{\,eV}$. Contacting Nb with
Bi2Te3 leads to a charge transfer from Nb to the TI of around $0.15$ electrons
($\approx 0.87$ electrons per nm2), resulting in an upward shift of the Fermi
level. This Fermi-level shift is obvious from the band structure shown in Fig.
1(c) where the color code indicates the localization of the states. The Fermi
level resides inside the conduction band of the TI as shown schematically in
Fig. 1(d). Furthermore, the hybridization of the TSS in the contact QL with
states from the Nb film strongly affects the dispersion of the surface state,
which is visible in an upward shift and an opening of the Dirac cone formed by
the TSS.

Figure 1: Calculational setup and normal-state electronic structure of a
topological insulator – superconductor (SC/TI) heterostructure. (a)
Calculational setup of $n$ quintuple layers (QLs) Bi2Te3 ($n=2-10$) in contact
with 6 layers Nb(111). Indicated are the “contact” and the “free” QLs of the
Bi2Te3 film. (b) Density of states integrated in the regions of the free QL,
the contact QL and the Nb region for the 3 QL SC/TI heterostructure shown in
(a). (c) Electronic band structure of the SC/TI heterostructure where the
localization of the wave functions is given as red (localized in Bi2Te3) or
blue (localized in Nb) colors. The Dirac point of the topological surface
state (TSS) on the free Bi2Te3 side is seen around
$E\approx-0.45\mathrm{\,meV}$ at the $\Gamma$-point. (d) Schematic
representation of the changes to the TSS upon bringing the TI in contact to
Nb.

###  II.2 Superconducting proximity effect in ultra-thin SC/TI
heterostructures

We start the investigation of the superconducting proximity effect in the
topological insulator with the thinnest TI film we consider here. For two QLs
Bi2Te3 the free and the contact side are not completely decoupled as seen, for
example, in the sizable hybridization-induced gap at the Dirac points in the
TSS of freestanding TIs in the ultra-thin limit Zhang  _et al._ (2010). Figure
2 shows the density of states in the superconducting state of the two QLs of
Bi2Te3 and the Nb regions in a narrow energy window around the Fermi level.
The superconducting $s$-wave pairing in Nb leads to the well-known opening of
a gap in the density of states of the Nb layers. The principal coherence peak
of Nb is used to identify the size of the superconducting gap in the
superconductor region (see Fig. 2c) which is called $\Delta_{0}$. The rich
structure in the electronic bands at the complex SC/TI interface, together
with the hybridization of the states and their partial localization in Nb or
the TI region, weaken the superconductivity at the interface and lead to a
soft superconducting gap where the DOS does not abruptly drop to zero for
energies $|E-E_{\mathrm{F}}|<\Delta_{0}$. A second gap at roughly half the
size of intrinsic Nb gap (highlighted with green dashed lines in Fig. 2b)
arises from the contact QL of the Bi2Te3 film. This proximity induced gap is
sizable due to the strong overlap of the wave functions in the interface
regions. Finally, in the DOS of the free QL a third gap is evident which
reaches only a fifth of the value of the Nb gap (orange dashed lines). We
attribute this to the surface states arising on the free side that only has an
exponential tail of its wave functions extend into the superconductor region.
Therefore, the proximity induced gap in the free side is significantly reduced
compared to the contact side.

Further insides into the proximity-induced gap on both sides of the TI film
can be gained from an analysis of the superconducting band structure. Figure 3
shows the band structure of 2 QLs Bi2Te3 on Nb around the Fermi energy. The
color code indicates the localization of the wave functions and panels (b-j)
show magnified views into selected intersections of bands with the Fermi
energy where the superconducting energy gap opens up. We recover the three
identified gap sizes from the DOS plots in Fig. 2 in states that are,
respectively, mainly localized in Nb (blue bands in Fig. 3b,c), states that
are localized in the contact QL (panel e) or states mainly localized in the
free side of the Bi2Te3 thin film (panels g,h). We point out that the states
at (d) and (e,f) arise from the TSS whereas the states at (h,i,j) are bulk-
like quantum-well states of the thin TI films. In agreement with the DOS
calculations discussed before, the induced gap in the TSS states in particular
in $\overline{\Gamma K}$ is significantly smaller than the gap in the quantum-
well states that couple much more strongly to the Nb states (panels i and j).
Furthermore, it is noteworthy that the induced gap in the TSS in the
orthogonal directions $\overline{\Gamma K}$ and $\overline{\Gamma M}$ differ
considerably (cf. Fig. 3f,g). The two states have a different Fermi velocity
which is inversely proportional to the density of states at the Fermi energy
$N(E_{\mathrm{F}})=\int_{\bm{k}\in\mathrm{FS}}(\hbar|\bm{v}_{\mathrm{F}}(\bm{k})|)^{-1}\,\mathrm{d}^{2}k$.
Consequently the superconducting gap size
$\Delta\sim\exp\bigl{(}-\frac{1}{\lambda N(E_{\mathrm{F}})}\bigr{)}$, that
follows from the BCS theory Bardeen  _et al._ (1957) with the electron-phonon
coupling constant $\lambda$, is expected to be larger for the crossing of the
TSS with $E_{\mathrm{F}}$ in $\overline{\Gamma M}$ compared to the
$\overline{\Gamma K}$ direction. Nevertheless, the localization of the wave
function (the state in $\overline{\Gamma M}$ has a longer wave function tail
into the Nb region) will also play a significant role in the formation of the
differently sized induced gap in the TSS and cannot be neglected.

Figure 2: Superconducting density of states for 2 QLs Bi2Te3 on Nb(111). (a-c)
integrated DOS in the free QL, the contact QL and the Nb region, respectively.
Three distinct gap sizes seen in the DOS are indicated with dashed vertical
lines which highlight the superconducting gap of Nb (red, $E=\pm\Delta_{0}$)
the principal induced gap in the contact quintuple layer (green, $E=\pm
0.53\Delta_{0}$) and the induced gap in the free QL (orange, $E=\pm
0.21\Delta_{0}$). Shown is only the particle component of the particle-hole
space used in the BdG formalism. Figure 3: Superconducting gap in different
bands of a 2QL Bi2Te3 / Nb(111) heterostructure. (a) Normal state band
structure around the $\Gamma$ point. (b-j) Gap in the superconducting state of
the bands around the intersections with the Fermi level marked in (a). The
three dashed horizontal lines indicate the features in the density of states
highlighted in Fig. 2. Note that the superconducting proximity gap in the
topological surface state is shown in the red bands for $\overline{\Gamma K}$
on the free side in (d) and for $\overline{\Gamma M}$ on the free and the
contact sides in (e) and (f), respectively.

###  II.3 Band bending and superconducting proximity effect in thick TI / SC
hetrostructures

We now turn our attention to the limit of thick TI films in contact with the
$s$-wave superconductor. Figure 4 shows the band structure of the 10 QL Bi2Te3
/ Nb(111) heterostructure where the contribution of different QLs from the
free side (QL#1) to the contact side with Nb (QL#10) are shown in (a). On the
free side the Fermi level resides within the bulk band-gap below the onset of
the conduction band marked with a purple arrow and only the TSS bridges the
gap to the valence band. On the contact side on the other hand the charge
transfer from the Nb contact shifts the Fermi level shifts upwards, which
fills a state at the bottom of the conduction band (highlighted with the
purple arrow). This effect is particularly strong in the topological insulator
because of the low density of states in the bulk band-gap that can screen the
excess charge. As shown in the calculated shift of the highest lying core
levels of the Bi and Te atom in Fig. 4(b), the screening of the charge doping
by the Nb contact decays over $\sim 4$ QLs after which the Fermi level also
approaches the bottom of the conduction band (as seen in the middle panel of
(a)).

Next, we discuss the influence of a variation of the Fermi level on the band
bending and the resulting superconducting gap in the electronic structure
shown schematically in Fig. 5. This models the effect of intrinsic charge
doping and allows to investigate the effect of charge puddles on the
electronic structure and the proximity effect in SC/TI heterostructures.
Figure 6 shows the electronic structure in the normal state (a-c) compared to
the superconducting state (d-f) for three locations of the Fermi level: (i)
$E_{\mathrm{F}}$ inside the TI valence band, (ii) $E_{\mathrm{F}}$ inside the
TI bulk band-gap as in Fig. 4, and (iii) $E_{\mathrm{F}}$ inside the TI
conduction band. Changing the location of the Fermi level in the TI region
shifts the bands of the TI relative to the Nb bands which influences the
hybridization of the TI states with the Nb states. Additionally, having the TI
Fermi level located inside the valence or conduction band leads to higher
density of states in the TI region when it comes in contact to the Nb so that
the charge doping of from the Nb contact is screened efficiently and less band
bending occurs.

The superconducting band structures are displayed in Fig. 6(d-f) for the three
Fermi levels in the TI are shown in terms of the integrated spectral density
(shown at $k_{x}<0$) that shows the particle-like components of the
Bogoliubov-quasiparticles in $\overline{\Gamma K}$-direction. We further
visualize the contributions to the anomalous density along the same path in
the Brillouin zone ($k_{x}>0$) that indicates the amount to which different
bands are proximitized. Here different bands pick up superconductivity due to
the varying degree of tunneling of Cooper pairs from the superconductor into
the respective TI bands because of the localization of the wavefunctions and
their orbital character. Comparing the calculations for different Fermi
levels, we find only minor changes in the superconducting gap of the bands
that derive from the Nb layers ($|k_{x}|>0.15\,\mathrm{\AA}$) and in the TSS
at the contact side (green ellipses). We conclude that the contact TSS
hybridizes strongly with the superconductor and thus is proximitized
efficiently with an induced gap size of roughly half of the gap of the
superconductor. On the other hand the TSS on the free side is not proximitized
as it is decoupled from the Nb contact through the thick TI film due to its
exponential decay from the free surface (see also Fig. A.1 in the supplemental
material).

Interestingly, for the Fermi level inside the TI bulk band-gap shown in Fig.
6(e), the filled state at the bottom of the conduction band is localized in
the interface region with a maximum in the second QL from the Nb contact with
only a small fraction of the wave function extending into the Nb region. This
effectively decouples the interface state (highlighted with the dashed
ellipses) from the superconductor so that the resulting proximity effect is
hardly visible. Finally, the conduction band states of the TI, that cross the
Fermi level in Fig. 6(f) in the region highlighted by the black ellipse, are
located within the bulk of the TI film and thus are not proximitized at all.

Figure 4: Band bending in $\sim 10\,\mathrm{nm}$ thick Bi2Te3 on Nb. (a) Band
structure contributions of first, fifth and tenth (i.e. contact) quintuple
layer (QL) which shows the band bending as a result of screening of the excess
charge transferred from the Nb contact to the topological insulator
(highlighted with purple arrows indicating the bottom of the conduction band
around $\overline{\Gamma}$). (b) Core level shifts of the highest lying core
state in the TI region (Te $4d$ and Bi $5d$). Shown are the contributions of
the 3 distinct Te layer, the Bi layers as well as the average value over the
QLs. Figure 5: Schematic illustration of the charge-transfer-induced band
bending in the TI region at different positions of the Fermi level located in
the bulk band gap, the valence band (VB) or the conduction band (CB). Note
that the electron transfer from Nb to the TI always leads to an upward shift
of the Fermi level at the SC/TI contact which is screened more efficiently
when the DOS in the TI region is higher (i.e. when $E_{\mathrm{F}}$ lies in VB
or CB). Figure 6: Superconducting band structure of 10 QL Bi2Te3 / Nb(111) for
varying Fermi level in the bulk region of the TI changing from (a) Fermi level
inside the bulk TI valence band, (b) inside the bulk band-gap and (c) inside
the conduction band. (d-f) zoomed views of (a-c) along $\overline{\Gamma K}$
comparing the spectral density ($k_{x}<0$) and the anomalous density
($k_{x}>0$). The purple arrows (green ellipses) in (d-f) highlight the
topological surface state on the free (contact) side of the SC/TI
heterostructure.

###  II.4 Decay of the superconducting proximity effect in the TSS

Finally, we discuss the decay of the superconducting proximity effect in the
TI film. We focus on the results from the 10 QL thick TI film as the thickest
film we considered in this work. Figure 7 shows the decay of the order
parameter of superconductivity (i.e. the anomalous density $\chi$) from the Nb
interface into the TI film. Shown are the result for the three positions of
the Fermi level in the bulk TI region discussed in the previous section.
Starting from the Nb layers on the left-hand side the order parameter quickly
drops and characteristic oscillations in the atoms within the quintuple layers
are seen. In the vacuum region the order parameter then drops exponentially as
only the exponential tail of the TI wave functions penetrate into the vacuum.

In order to quantify the decay of the order parameter we fit a polynomial
decay $f(z)=\alpha z^{-\beta}$ which we find to match our data much better
than an exponential decay. In the fit we exclude the first two QLs as well as
the vacuum region in order to analyze the long-distance limit without
exaggerating effects from the region where band bending plays a major role.
When the bulk TI Fermi level resides in the bulk band-gap or at the top of the
valence band we find a fast decay ($\sim z^{-\beta}$) of the proximity effect
in the superconducting order parameter with exponents $\beta=2.68$ and
$\beta=2.76$, respectively. On the other hand the Fermi level lies within the
conduction band of the TI, the decay of the order parameter approaches the
slow decay of $\sim 1/z$ known from metal/SC interfaces Deutscher and de
Gennes (1969); Falk (1963). We attribute the faster decay of the order
parameter for the Fermi level residing in the bulk band-gap or at the edge of
the valence band with the significantly reduced number density of states in
these circumstances which make the SC/TI interface behave rather like a
semiconductor/SC interface than a metal/SC interface.

Figure 7: Decay of the superconducting order parameter (i.e. the average
anomalous density) for 10 QL thick Bi2Te3 in contact to Nb. The shaded regions
indicate the Nb and the vacuum regions and the vertical dotted lines
highlights the beginning of the quintuple layers. The dashed lines show fits
with functions $f(z)=\alpha z^{-\beta}$ where the fit exponent $\beta$ is
given in the figure legend.

##  III Discussion

Comparing our results to existing data in the literature, we find that the
slow polynomial decay of the proximity effect in the order parameter is well
in line with previous model-based studies of the proximity effect in SC/TI
heterostructures Chiu  _et al._ (2016). Our results show that the localization
of the wave functions in the SC/TI heterostructure determines the possibility
for Cooper pair tunneling into the respective states. The TSS on the free side
away from the SC contact is localized at the interface to the vacuum and only
has an exponentially decaying tail of its wave function into the TI.
Consequently the proximity gap in the free TSS shows an exponential decay with
the TI thickness which is in line with earlier observations of the decaying
proximity gap as seen for example in experiments Xu  _et al._ (2015, 2014a);
Stolyarov  _et al._ (2021); Flötotto  _et al._ (2018); Xu  _et al._ (2014b).

On the other hand, the absence of the proximity effect in some quantum-well
states, which we similarly attribute to a negligible hybridization with Nb, is
not always found experimentally Flötotto  _et al._ (2018). However, this
apparent discrepancy could be an artifact of the structural model we chose in
our computational approach. In general, the SC/TI interface between TIs and
elemental superconductors is amorphous and poly-crystalline SC phases would be
a more realistic description of the interface which is however computationally
not feasible. A proper treatment would thus encompass extremely large unit
cells with thick amorphous superconducting regions that show some intermixing
at the SC/TI interface. The defect-induced broadening of the states and the
disordered nature of the interface can result in a quasi-continuum of states
of SC states throughout the Brillouin zone with a slower decay into the
surface TI region. Such a situation could therefore offer more possibilities
for hybridization with the TI electronic structure which would result in more
efficient Cooper pair tunneling into the TI and thus a larger induced
proximity gap in the TI bands compared to our present results.

In contrast to amorphous interfaces with elemental SCs, the interfaces between
vdW materials can in principle be much cleaner Novoselov  _et al._ (2016),
which will have consequences for the charge doping at the SC/TI interface. One
can expect that the presence of the vdW gap between the different 2D materials
(TI and SC) will lessen the charge doping which was also not reported in
earlier studies using PdTe2 and NbSe2 to proximitize a TI Xu  _et al._ (2014a,
b, 2015); Park  _et al._ (2020). Depending on the fabrication process a higher
interface quality interface is however not a given as seen in the self-formed
PdTe2/(Bi,Sb)2Te3 heterostructure where considerable intermixing is present at
the interface Bai  _et al._ (2020).

In future studies, we plan to investigate the effect of intermixing and
diffusion of superconductor atoms into the surface QLs of the TI and
scrutinize its effect on the proximity effect. Such a dusting with a moderate
concentration of superconducting atoms is believed to stabilize the proximity
effect as shown in the very long observation of the superconducting gap
amorphous atomic Pb layer on Bi2Te3 around Pb islands Stolyarov  _et al._
(2021). We furthermore plan to probe the influence of buffer layers that can
modify (e.g. partially screen) the charge doping from SC contact to the TI.

In conclusion, our results contain mixed messages for the ability to use
elemental SC/TI interfaces as a platform for Majorana-based topological
qubits. Encouraging is the sizable SC gap-opening we find in the contact TSS
(in the order of half of the gap size of Nb), which is even robust against
variations of the Fermi level in the TI region. Therefore, the proximity
effect in the TSS in contact with the SC is also expected to be robust to
fluctuations in the Fermi level, as they occur in the presence of charge
puddles, which are a common problem in TI materials Beidenkopf  _et al._
(2011); Skinner  _et al._ (2012); Bömerich  _et al._ (2017). Moreover, the
charge transfer and the resulting band bending in thick TI films contacted
with Nb are comparable with recent model calculations that showed the
beneficial effect of metallization for the technologically important geometry
of a proximitized TI nanowire Legg  _et al._ (2022). Detrimental for Nb/Bi2Te3
as a platform for topological superconductivity could be the trivial quantum-
well states in the TI, which can become occupied due to band bending or
intrinsic doping. Depending on their localization at the SC/TI interface, the
Cooper pair tunneling into these states is strongly suppressed, which leads to
very small SC gap-openings in these trivial states. Our findings may help in
understanding and optimizing the interface between TIs and superconductors as
a suitable material platform in the field of topologically protected quantum
computing.

##  IV Methods

The density functional theory (DFT) results of this work are produced using
the full-potential relativistic Korringa-Kohn-Rostoker Green’s function method
(KKR) Ebert  _et al._ (2011) as implemented in the JuKKR code package The
JuKKR developers (2021). The local density approximation (LDA) is used to
parameterize the exchange correlation functional Vosko  _et al._ (1980). We
use $80\times 80$ $k$-points in the Brillouin zone integration during self-
consistency, spherical screening clusters of $1$ nm radius around each atom
(containing $120-180$ sites), an $\ell_{max}=2$ cutoff in the angular momentum
expansion of the Green’s function with an exact partitioning of the unit cell
into atomic cells Stefanou  _et al._ (1990); Stefanou and Zeller (1991), and
correct the truncation error arising from the finite $\ell_{max}$ cutoff using
Lloyd’s formula Zeller (2004). A Lloyd’s formula-like renormalization of the
energy integration weights is used to shift the Fermi level of the topological
insulator with respect to the Nb superconductor when the position of the Fermi
level is changed in the calculation. The TI/Nb(111) system in this study is
calculated in a symmetric TI(111)/Nb(111)/TI(111) geometry in order to
counteract spurious shifts in the electronic structure arising from charge-
transfer-induced dipole moments. In order to facilitate the calculations we
construct the TI/Nb interface starting from the experimental crystal structure
of Bi2Te3 (111) of 4.4 Å Atuchin  _et al._ (2012). Then 6 layers of Nb(111)
scaled ($-6\%$ compared to Nb equilibrium lattice constant) to the same in-
plane lattice constant as the TI film are placed in a distance of 2.6 Å (i.e.
the TI inter-layer distance). Because the focus of our investigation is the
induced superconductivity in the electronic structure of the TI film we
neglect the effect of lattice relaxations at the interface.

The superconducting properties are calculated based on the Bogoliubov-de
Gennes (BdG) extension to the KKR method Csire  _et al._ (2015); Rüßmann and
Blügel (2022a) which is abbreviated as “KKR-BdG” in this work. In the KKR-BdG
method, the Bogoliubov-de Gennes Hamiltonian (given here in atomic units)

$H_{\mathrm{BdG}}(\bm{x})=\left(\begin{array}[]{cc}-\nabla^{2}-E_{\mathrm{F}}+V(\bm{x})&\mathcal{D}(\bm{x})\\\
\mathcal{D}^{*}(\bm{x})&\nabla^{2}+E_{\mathrm{F}}-V^{*}(\bm{x})\end{array}\right)$
(1)

is solved self-consistently within the framework of density functional theory.
Here, $\bm{x}$ denotes the position in space, $V(\bm{x})$ is the normal state
electronic potential and $\mathcal{D}(\bm{x})$ is the superconducting pairing
potential responsible for the formation of Cooper pairs. Following the
parametrization of the exchange-correlation functional introduced by Suvasini
et al. Suvasini  _et al._ (1993), the superconducting pairing potential in
atom $i$ can be written as
$\mathcal{D}_{i}(\bm{x})=\lambda_{i}\chi_{i}(\bm{x})$ with a set of atom-
dependent semi-phenomenological coupling constants $\lambda_{i}$ and the
anomalous density $\chi_{i}(\bm{x})$. In analogy to the BCS theory of
superconductivity, the coupling constants $\lambda_{i}$ can be interpreted as
the electron-phonon coupling, which results in the attractive interaction that
form the Cooper pairs in conventional superconductors. We neglect the
intrinsic electron-phonon coupling in the TI region and model the
superconducting state with coupling constants

$\lambda_{i}=\left\\{\begin{array}[]{rl}\lambda_{0}&\mathrm{if}\
i\in\mathrm{sites\ of\ Nb\ atoms}\\\ 0&\mathrm{else}\end{array}\right.$ (2)

where $\lambda_{0}$ is a positive real-valued constant. This is motivated by
the smallness of the electron-phonon coupling in 3D TIs of the Bi2Te3-family
that is found to be in the order of 0.1 Heid  _et al._ (2017) whereas in Nb it
is an order of magnitude larger Bauer  _et al._ (1998). The pairing potential
in atom $i$ is calculated self-consistently using

$\mathcal{D}_{i}(\bm{x})=\lambda_{i}\chi_{i}(\bm{x})$ (3)

with the anomalous density $\chi$ that is the order parameter of the
superconductivity. This results in a finite intrinsic superconducting pairing
potential in the Nb layers and vanishing intrinsic pairing potential in the TI
layers. More details of the KKR-BdG method can be found in Ref. Rüßmann and
Blügel, 2022a. In this work, we scale the coupling constant such that the
superconducting gap in Nb is of the order of 1 mRy to facilitate the
resolution of the proximity-induced superconducting gaps in Bi2Te3. The
induced gaps can then be scaled with respect to the experimental Nb gap size
of 0.12 mRy (1.6 meV) Odobesko  _et al._ (2019), where the scaling constant
for the smaller gaps Bi2Te3 should remain the same. We also varied the model
for the coupling constant to stay constant throughout the heterostructure with
the same value in the TI region and the SC region. This however does not
affect the outcome of the results (i.e. the band structures and the decay of
the anomalous density stay unaffected), which highlights the resilience of our
conclusions with respect to this approximation.

The series of DFT calculations in this study are facilitated through the
AiiDA-KKR plugin Rüßmann  _et al._ (2021, 2021) to the AiiDA infrastructure
Huber  _et al._ (2020). This has the advantage that the data provenance is
automatically stored in compliance to the FAIR principles of open data
Wilkinson  _et al._ (2016). The complete data set that includes the full
provenance of the calculations is made publicly available in the materials
cloud repository Talirz  _et al._ (2020); Rüßmann and Blügel (2022b).

## Data Availability

The data was generated using the JuKKR code (https://jukkr.fz-juelich.de) via
the AiiDA-KKR plugin Rüßmann  _et al._ (2021, 2021) to the AiiDA
infrastructure Huber  _et al._ (2020) and are publicly available on the
materialscloud archive Talirz  _et al._ (2020); Rüßmann and Blügel (2022b).

## Code Availability

The source codes of the AiiDA-KKR plugin Rüßmann  _et al._ (2021, 2021) and
the JuKKR code The JuKKR developers (2021) are published as open source
software under the MIT license at https://github.com/JuDFTteam/aiida-kkr and
https://iffgit.fz-juelich.de/kkr/jukkr, respectively.

## Acknowledgements

We thank Achim Rosch and Henry Legg for fruitful discussions about the effect
of band bending in the SC/TI heterostructure. We acknowledge support by the
Joint Lab Virtual Materials Design (JL-VMD) and we are grateful for computing
time granted by the JARA Vergabegremium and provided on the JARA Partition
part of the supercomputer CLAIX at RWTH Aachen University (project number
jara0191). This work was furthermore funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – Cluster of Excellence Matter and Light for Quantum
Computing (ML4Q) EXC 2004/1 – 390534769, and we thank the Bavarian Ministry of
Economic Affairs, Regional Development and Energy for financial support within
High-Tech Agenda Project “Bausteine für das Quantencomputing auf Basis
topologischer Materialien mit experimentellen und theoretischen Ansätzen“.

## Competing Interests

The authors declare no competing financial or non-financial interests.

##  Appendix A Supplement

###  A.1 Wave function localization

Figure A.1 shows the superconducting band structure of 10QL Bi2Te3
proximitized by an interface to Nb(111) which is compared to the localization
of the wave functions throughout the heterostructure. We find that the states
that originate from the Nb superconductor and penetrate only little into the
TI region ($|k_{x}|>0.1\,\mathrm{\AA}$) have a large superconducting gap. The
topological surface that is located on the free side (highlighted by purple
arrows) on the other hand has no wave function overlap with the Nb region and
therefore is not proximitized. The interface states that form due to charge-
transfer induced band bending at the TI/Nb interface also only has a marginal
extend of the wave function with the Nb region and therefore only picks up a
tiny superconducting gap as discussed in the main text.

Figure A.1: Superconducting band structure (a) and wave function localization
(b) of 10QL Bi2Te3 proximitized by an interface with Nb(111). The band
structure in (a) is integrated over the full heterostructure and the wave
function localization in (b) is shown along the blue dash-dotted line in (a).
The red and grey dashed lines in (b) separate the regions of the
superconductor ($z<3\,\mathrm{\AA}$), the TI film ($3\,\mathrm{\AA}\leq z\geq
105\,\mathrm{\AA}$) and the vacuum ($z>105\,\mathrm{\AA}$) in the
heterostructure. The purple arrows at $k_{x}\approx 0.08\,\mathrm{\AA}$
highlight the topological surface state on the free side. The green arrows at
$k_{x}\approx 0.04\,\mathrm{\AA}$ indicate the interface states.

## References

  * Fu and Kane (2008) L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * Cook and Franz (2011) A. Cook and M. Franz, Phys. Rev. B 84, 201105(R) (2011).
  * Alicea  _et al._ (2011) J. Alicea, Y. Oreg, G. Refael, F. von Oppen, and M. P. A. Fisher, Nature Phys. 7, 412 (2011).
  * Sato and Ando (2017) M. Sato and Y. Ando, Rep. Prog. Phys. 80, 076501 (2017).
  * Flensberg  _et al._ (2021) K. Flensberg, F. von Oppen, and A. Stern, Nat. Rev. Mater. 6, 944–958 (2021).
  * de Leon  _et al._ (2021) N. P. de Leon, K. M. Itoh, D. Kim, K. K. Metha, T. E. Northup, H. Paik, B. S. Palmer, N. Samarth, S. Sangtawesin, and D. W. Steuerman, Science 372, eabb2823 (2021).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Breunig and Ando (2022) O. Breunig and Y. Ando, Nat. Rev. Phys. 4, 184 (2022).
  * Schmitt  _et al._ (2022) T. W. Schmitt, M. R. Connolly, M. Schleenvoigt, C. Liu, O. Kennedy, J. M. Chávez-Garcia, A. R. Jalil, B. Bennemann, S. Trellenkamp, F. Lentz, E. Neumann, T. Lindström, S. E. de Graaf, E. Berenschot, N. Tas, G. Mussler, K. D. Petersson, D. Grützmacher, and P. Schüffelgen, Nano Lett. 22, 2595 (2022).
  * Xu  _et al._ (2015) J. Xu, M. Wang, Z. L. Liu, J. Ge, X. Yang, C. Liu, Z. A. Xu, D. Guan, C. L. Gao, D. Qian, Y. Liu, and Q. Wang, Phys. Rev. Lett. 114, 017001 (2015).
  * Xu  _et al._ (2014a) J.-p. Xu, C. Liu, M.-x. Wang, J. Ge, Z.-l. Liu, X. Yang, and Y. Chen, Phys. Rev. Lett. 112, 217001 (2014a).
  * Stolyarov  _et al._ (2021) V. S. Stolyarov, S. Pons, S. Vlaic, S. V. Remizov, D. S. Shapiro, C. Brun, S. I. Bozhko, T. Cren, T. V. Menshchikova, E. V. Chulkov, W. V. Pogosov, Y. E. Lozovik, and D. Roditchev, J. Phys. Chem. Lett. (2021), 10.1021/acs.jpclett.1c02257.
  * Flötotto  _et al._ (2018) D. Flötotto, Y. Ota, Y. Bai, C. Zhang, K. Okazaki, A. Tsuzuki, T. Hashimoto, J. N. Eckstein, S. Shin, and T.-c. Chiang, Sci. Adv. 4, eaar7214 (2018).
  * Xu  _et al._ (2014b) S.-Y. Xu, N. Alidoust, I. Belopolski, A. Richardella, C. Liu, M. Neupane, G. Bian, S.-h. Huang, R. Sankar, C. Fang, B. Dellabetta, W. Dai, Q. Li, M. J. Gilbert, F. Chou, N. Samarth, and M. Z. Hasan, Nature Phys. 10, 943 (2014b).
  * Bai  _et al._ (2020) M. Bai, F. Yang, M. Luysberg, J. Feng, A. Bliesener, G. Lippertz, A. A. Taskin, J. Mayer, and Y. Ando, Phys. Rev. Mat. 4, 094801 (2020).
  * Chiu  _et al._ (2016) C.-K. Chiu, W. S. Cole, and S. Das Sarma, Phys. Rev. B 94, 1 (2016).
  * Park  _et al._ (2020) K. Park, G. Csire, and B. Újfalussy, Phys. Rev. B 102, 134504 (2020).
  * Legg  _et al._ (2022) H. F. Legg, D. Loss, and J. Klinovaja, Phys. Rev. B 105, 155413 (2022).
  * Rüßmann and Blügel (2022a) P. Rüßmann and S. Blügel, Phys. Rev. B 105, 125143 (2022a).
  * Bahramy  _et al._ (2012) M. Bahramy, P. King, A. de la Torre, J. Chang, M. Shi, L. Patthey, G. Balakrishnan, P. Hofmann, R. Arita, N. Nagaosa, and F. Baumberger, Nat. Commun. 3, 1159 (2012).
  * Michiardi  _et al._ (2022) M. Michiardi, F. Boschini, H.-H. Kung, M. X. Na, S. K. Y. Dufresne, A. Currie, G. Levy, S. Zhdanovich, A. K. Mills, D. J. Jones, J. L. Mi, B. B. Iversen, P. Hofmann, and A. Damascelli, Nat. Commun. 13, 3096 (2022).
  * Zhang  _et al._ (2010) Y. Zhang, K. He, C.-Z. Chang, C.-L. Song, L.-L. Wang, X. Chen, J.-F. Jia, Z. Fang, X. Dai, W.-Y. Shan, S.-Q. Shen, Q. Niu, X.-L. Qi, S.-C. Zhang, X.-C. Ma, and Q.-K. Xue, Nature Phys. 6, 584–588 (2010).
  * Bardeen  _et al._ (1957) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957).
  * Deutscher and de Gennes (1969) G. Deutscher and P. de Gennes,  _Superconductivity: Proximity Effects_ , edited by R. D. Parks, Vol. 2 (Marcel Dekker, New York, 1969).
  * Falk (1963) D. S. Falk, Phys. Rev. 132, 1576 (1963).
  * Novoselov  _et al._ (2016) K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. C. Neto, Science 353, 6298 (2016).
  * Beidenkopf  _et al._ (2011) H. Beidenkopf, P. Roushan, J. Seo, L. Gorman, I. Drozdov, Y. S. Hor, R. J. Cava, and A. Yazdani, Nature Phys. 7, 939–943 (2011).
  * Skinner  _et al._ (2012) B. Skinner, T. Chen, and B. I. Shklovskii, Phys. Rev. Lett. 109, 176801 (2012).
  * Bömerich  _et al._ (2017) T. Bömerich, J. Lux, Q. T. Feng, and A. Rosch, Phys. Rev. B 96, 075204 (2017).
  * Ebert  _et al._ (2011) H. Ebert, D. Ködderitzsch, and J. Minár, Rep. Prog. Phys. 74, 096501 (2011).
  * The JuKKR developers (2021) The JuKKR developers, “The Jülich KKR Codes,” (2021), https://jukkr.fz-juelich.de.
  * Vosko  _et al._ (1980) S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).
  * Stefanou  _et al._ (1990) N. Stefanou, H. Akai, and R. Zeller, Comput. Phys. Commun. 60, 231 (1990).
  * Stefanou and Zeller (1991) N. Stefanou and R. Zeller, J. Phys.: Cond. Matter 3, 7599 (1991).
  * Zeller (2004) R. Zeller, J. Phys.: Condens. Matter 16, 6453 (2004).
  * Atuchin  _et al._ (2012) V. Atuchin, T. Gavrilova, K. Kokh, N. Kuratieva, N. Pervukhina, and N. Surovtsev, Solid State Comm. 152, 1119 (2012).
  * Csire  _et al._ (2015) G. Csire, B. Újfalussy, J. Cserti, and B. Győrffy, Phys. Rev. B 91, 165142 (2015).
  * Suvasini  _et al._ (1993) M. B. Suvasini, W. M. Temmerman, and B. L. Györffy, Phys. Rev. B 48, 1202 (1993).
  * Heid  _et al._ (2017) R. Heid, I. Y. Sklyadneva, and E. V. Chulkov, Sci. Rep. 7, 1095 (2017).
  * Bauer  _et al._ (1998) R. Bauer, A. Schmid, P. Pavone, and D. Strauch, Phys. Rev. B 57, 11276 (1998).
  * Odobesko  _et al._ (2019) A. B. Odobesko, S. Haldar, S. Wilfert, J. Hagen, J. Jung, N. Schmidt, P. Sessi, M. Vogt, S. Heinze, and M. Bode, Physical Review B 99, 115437 (2019).
  * Rüßmann  _et al._ (2021) P. Rüßmann, F. Bertoldo, J. Bröder, J. Wasmer, R. Mozumder, J. Chico, and S. Blügel, Zenodo (2021), https://doi.org/10.5281/zenodo.3628251.
  * Rüßmann  _et al._ (2021) P. Rüßmann, F. Bertoldo, and S. Blügel, npj Comput. Mater. 7, 13 (2021).
  * Huber  _et al._ (2020) S. P. Huber, S. Zoupanos, M. Uhrin, L. Talirz, L. Kahle, R. Häuselmann, D. Gresch, T. Müller, A. V. Yakutovich, C. W. Andersen, F. F. Ramirez, C. S. Adorf, F. Gargiulo, S. Kumbhar, E. Passaro, C. Johnston, A. Merkys, A. Cepellotti, N. Mounet, N. Marzari, B. Kozinsky, and G. Pizzi, Sci. Data 7, 300 (2020).
  * Wilkinson  _et al._ (2016) M. D. Wilkinson, M. Dumontier, I. J. Aalbersberg, G. Appleton, M. Axton, A. Baak, N. Blomberg, J.-W. Boiten, L. B. da Silva Santos, P. E. Bourne, J. Bouwman, A. J. Brookes, T. Clark, M. Crosas, I. Dillo, O. Dumon, S. Edmunds, C. T. Evelo, R. Finkers, A. Gonzalez-Beltran, A. J. G. Gray, P. Groth, C. Goble, J. S. Grethe, J. Heringa, P. A. C. ’t Hoen, R. Hooft, T. Kuhn, R. Kok, J. Kok, S. J. Lusher, M. E. Martone, A. Mons, A. L. Packer, B. Persson, P. Rocca-Serra, M. Roos, R. van Schaik, S.-A. Sansone, E. Schultes, T. Sengstag, T. Slater, G. Strawn, M. A. Swertz, M. Thompson, J. van der Lei, E. van Mulligen, J. Velterop, A. Waagmeester, P. Wittenburg, K. Wolstencroft, J. Zhao, and B. Mons, Sci. Data 3, 160018 (2016).
  * Talirz  _et al._ (2020) L. Talirz, S. Kumbhar, E. Passaro, A. V. Yakutovich, V. Granata, F. Gargiulo, M. Borelli, M. Uhrin, S. P. Huber, S. Zoupanos, C. S. Adorf, C. W. Andersen, O. Schütt, C. A. Pignedoli, D. Passerone, J. VandeVondele, G. P. Thomas C. Schulthess, Berend Smit, and N. Marzari, Sci. Data 7, 299 (2020).
  * Rüßmann and Blügel (2022b) P. Rüßmann and S. Blügel, “Materials Cloud Archive 2022.XXX,”  (2022b), to be added during review.