Gamma rays from dark matter spikes in EAGLE simulations

Intermediate Mass Black Holes (IMBHs) with a mass range between $100 \, \text{M}_\odot$ and $10^6 \, \text{M}_\odot$ are expected to be surrounded by high dark matter densities, so-called dark matter spikes. The high density of self-annihilating Weakly Interacting Massive Particles (WIMPs) in these spikes leads to copious gamma-ray production. Sufficiently nearby IMBHs could therefore appear as unidentified gamma-ray sources. However, the number of IMBHs and their distribution within our own Milky Way is currently unknown. In this work, we provide a mock catalogue of IMBHs and their dark matter spikes obtained from the EAGLE simulations, in which black holes with a mass of $10^5 \, \text{M}_\odot/h$ are seeded into the centre of halos greater than $10^{10} \, \text{M}_\odot/h$ to model black hole feedback influencing the formation of galaxies. The catalogue contains the coordinates and dark matter spike parameters for over 2500 IMBHs present in about 150 Milky Way-like galaxies. We expect about $15^{+9}_{-6}$ IMBHs within our own galaxy, mainly distributed in the Galactic Centre and the Galactic Plane. In the most optimistic scenario, we find that current and future gamma-ray observatories, such as Fermi-LAT, H.E.S.S. and CTAO, would be sensitive enough to probe the cross section of dark matter self-annihilation around IMBHs down to many orders of magnitude below the thermal relic cross section for dark matter particles with masses from GeV to TeV. We have made the IMBH mock catalogue and the source code for our analysis publicly available, providing the resources to study dark matter self-annihilation around IMBHs with current and upcoming gamma-ray observatories.


Introduction
The nature of dark matter is one of the most important questions in fundamental physics [1][2][3][4].One of the most popular candidates for dark matter are Weakly Interacting Massive Particles (WIMPs), which naturally arise in well-motivated extensions of the Standard Model of particle physics [5][6][7][8].The thermal production of WIMPs at the measured relic density of Ω DM = 0.26 implies the production of Standard Model particles, including gamma rays, neutrinos, and anti-particles through the self-annihilation cross section Ω DM h 2 ≈ 3 × 10 −27 cm 3 s −1 /⟨σv⟩ [9,10].This canonical annihilation cross section, approximated at the time of chemical decoupling, is applicable primarily for indirect detection scenarios in case of velocity-independent (s-wave) processes, and may not necessarily be probed by direct dark matter searches.The indirect detection of self-annihilating dark matter includes searches for gamma-ray emission from regions with large WIMP number densities n WIMP , as the self-annihilation rate scales with n 2 WIMP .A lot of effort has been put into searching for a dark matter signal in regions with large dark matter number densities, including the Galactic Centre [11][12][13][14][15][16], dwarf galaxies [17][18][19][20] and galaxy clusters [21][22][23].Another very promising target for indirect dark matter searches are the environments of black holes dominated by their gravitational potential.Gondolo and Silk (1999) [24] applied the formalism developed for the adiabatic growth of power-law stellar cusps by Quinlan, Hernquist and Sigurdsson (1995) [25], to the dark matter density profile around the supermassive black hole (SMBH) at the Galactic Centre, Sgr A * .Since the dark matter annihilation rate is proportional to the squared dark matter density ρ 2 DM , these adiabatically grown profiles, dubbed spikes, would lead to a significant enhancement of the gamma-ray signal.Subsequent studies explored the implications of the existence of spikes around SMBHs [26,27], deriving constraints on the dark matter cross section based on the spike around Sgr A * [26,[28][29][30][31][32][33].Zhao and Silk (2005) [34] and Bertone, Zentner, and Silk (2005) [35], proposed spikes around Intermediate Mass Black Holes (IMBHs) [36] as promising targets for indirect dark matter searches.IMBHs cover a mass range between ∼ 10 2 − 10 6 M ⊙ and are expected to form via gravitational runaway [37], the direct collapse of primordial gas in early forming halos [38], as remnants of Population III stars [39,40] or from primordial black holes (PBHs) [41].Potential IMBH candidates have been found in ultraluminous X-ray sources (ULXs) [42].These sources have luminosities above the Eddington limit for compact objects with M ≲ 20 M ⊙ and can therefore not be explained by neutron stars and stellar mass black holes [36].The most well-studied and recognized systems are ω Centauri and 47 Tucanae, for which black hole masses in the order of 10 4 M ⊙ have been found [43,44].However, additional observational data are needed to confirm these measurements [36].The first conclusive evidence for an IMBH is the detection of the gravitational waves from the binary black hole merger event GW190521 with an inferred mass of the remnant black hole of 142 +28 −16 M ⊙ [45].However, the observable black hole mass range of these detectors is currently limited, spanning from a few solar masses up to a couple of hundred solar masses [46].Therefore, current detectors do not allow to probe IMBHs in the M ≳ 10 3 M ⊙ mass range.The upcoming LISA experiment [47] will cover the 10 −4 − 10 −2 Hz frequency range of gravitational waves enabling the detection of IMBHs above the M ≳ 10 3 M ⊙ mass regime [48].The distributions and luminosity of IMBHs and their dark matter spikes have been estimated by Bertone, Zentner and Silk (2005) [35], showing that the spikes may be detected as unidentified point-like gamma-ray sources.Dark matter spikes around IMBHs in the Milky Way would result in tens or more point-like sources with identical energy spectra, which would make them a smoking gun signature for dark matter annihilation [35,49].Early gamma-ray searches from dark matter annihilation around IMBHs have been reported in Aharonian et al. (2008) [50] and Bringmann, Lavalle and Salati (2009) [51].Nowadays, the advancements in cosmological simulations, such as the EAGLE [52] and IllustrisTNG [53] simulations, provide a more comprehensive and refined understanding of the impact of massive black holes and their associated phenomena within their host galaxies.However, it is important to acknowledge the substantial theoretical uncertainties associated with the formation and evolution of IMBHs inherent in these simulations.Recent studies, such as Fujii et al. (2024) [54], have highlighted the complexities of IMBH formation, particularly the intricate baryonic processes that occur on small scales [55][56][57].In their work, Fujii et al. conducted high-resolution hydrodynamic simulations of giant molecular clouds, demonstrating how IMBHs can form from the mergers of very massive stars in dense star clusters, where both the total mass and the metallicity of the cloud play crucial roles.These studies suggest that current cosmological simulations do not fully capture these detailed processes due to their limited spatial and temporal resolutions.Nevertheless, integrating insights from cosmological simulations allows us to build a more comprehensive picture of the distribution of IMBHs and their corresponding dark matter spikes within the Milky Way.In this work, we provide a mock catalogue of IMBHs within Milky Way-like galaxies using the EAGLE simulations [52].The catalogue provides information about the expected number of IMBHs, their spatial distribution and their dark matter spike parameters, including the expected gamma-ray flux from dark matter self-annihilation.The IMBH catalogue, the catalogue of our selection of Milky Way-like galaxies within the EAGLE simulations and the source code are publicly available at [58] and [59].We perform our analysis within the framework of a flat Λ Cold Dark Matter (ΛCDM) cosmology, following the parameters from the Planck mission [9], i.e.Ω m = 0.307, Ω Λ = 0.693, Ω b = 0.04825, h = 0.6777, σ 8 = 0.8288 and n s = 0.9611.The structure of this article is as follows: In Section 2, we introduce the theoretical background of the dark matter spikes and the gamma-ray flux expected from dark matter annihilation.We describe the EAGLE simulations and our selection criteria for Milky Way-like galaxies in Section 3. The properties of selected EAGLE galaxies, the analysis steps to extract the IMBH coordinates and the dark matter spike parameters for each Milky Way-like galaxy are described in Section 4. We present our results and discussion for the detectability of a dark matter annihilation signal around IMBHs in Section 5 & Section 6.Finally, we conclude our findings in Section 7.

Dark matter spikes
We start with describing the theoretical framework of dark matter spikes around intermediate mass black holes.This framework is used to calculate the expected gamma-ray flux from dark matter self-annihilation around IMBHs in Section 4 & 5.

Dark matter density surrounding IMBHs
We follow the theoretical framework of Bertone, Zentner and Silk (2005) [35] to calculate the dark matter density profile surrounding IMBHs.We parametrise the dark matter density profile as follows with the Schwarzschild radius r schw = 2Gm BH /c 2 , the mass of the black hole m BH , the spike radius r sp , the cutoff radius r cut , the dark matter weak cusp density profile ρ wcusp (r), the dark matter spike density profile ρ sp (r) and the dark matter density profile of the host halo ρ halo (r), in which the black hole formed.Therefore, the dark matter density profile around an IMBH consists of four regions: 1.) the region inside 2r schw , in which we assume the dark matter density to be zero, 2.) a weak cusp between 2r schw and the cutoff radius, characterised by r −0.5 due to dark matter s-wave annihilation [60,61], 3.) the region between the cutoff radius and the spike radius, which corresponds to the dark matter spike and 4.) the region outside the spike radius, which follows the dark matter density profile of the host halo.We discuss these individual components and how to compute the dark matter spike parameters in the following.
In N-body cosmological simulations, the Navarro-Frenk-White (NFW) profile has been shown to describe the dark matter density profile of galaxies very well [62].Therefore, we assume that the dark matter density profile of the host halo ρ halo (r) follows a NFW profile [62] given by with the normalisation constant ρ 0 and the scale radius r s .The dark matter density profile ρ NFW (r) is used to calculate the radius of the sphere of gravitational influence r h of the black hole which is given by [63] M (< r h ) = 4π We follow the standard assumption of r sp = 0.2r h to determine the spike radius [32,63,64].
The dark matter spike density ρ sp (r) follows a power law with spike index γ sp and is given by The spike index γ sp is related to the initial power-law index γ of the dark matter density profile of the host halo by [24,25] For the NFW profile with γ = 1 the spike index reduces to γ sp = 7/3.The spike density ρ sp (r) diverges at small radii, however, the self-annihilation of dark matter particles results in a weak cusp for the dark matter density at the saturation radius r sat which is given by with the dark matter mass m χ , the dark matter annihilation cross section times the relative velocity ⟨σv⟩, the age of the universe t 0 and the formation time of the black hole t f [24].Therefore, the saturation radius r sat depends on the assumed dark matter particle.For r ≤ r sat and s-wave annihilation, the dark matter distribution follows [60,61] (2.7) Finally, we define the cutoff radius r cut as For dark matter masses in the GeV to TeV scale and typical dark matter cross sections of ⟨σv⟩ ∼ 10 −26 cm 3 s −1 , the saturation radius r sat is typically in the order of 10 −3 pc, so r cut = r sat .A typical dark matter density profile around an IMBH is shown in Figure 1.In this particular example, we assume a black hole mass of 1.6 × 10 5 M ⊙ , the NFW profile as dark matter halo Figure 1: Dark matter density profile around an IMBH assuming a black hole mass of 1.6 × 10 5 M ⊙ , the NFW profile as dark matter halo profile with ρ 0 = 1.9 GeV cm −3 and r s = 2.3 kpc, a spike radius of r sp = 4.2 pc and cutoff radius of r cut = 2.3 × 10 −3 pc.These values correspond to the parameters from a typical IMBH in our analysis.
profile with ρ 0 = 1.9 GeV cm −3 and r s = 2.3 kpc, a spike radius of r sp = 4.2 pc and cutoff radius of r cut = 2.3 × 10 −3 pc.These values correspond to the parameters from a specific and typical IMBH in our analysis, as we show later in Section 4.

Gamma-ray flux from dark matter annihilation
The dark matter density profile around IMBHs summarised above is crucial for predicting the expected gamma-ray flux from dark matter annihilation.As presented by Bertone, Zentner and Silk (2005) [35], the gamma-ray flux Φ from dark matter annihilation can be expressed as with the dark matter annihilation spectrum dN/dE and the distance D from the observer to the IMBH.We assumed r cut ≫ r schw and r sp ≫ r cut to simplify the equation.Unlike the work by Bertone, Zentner and Silk (2005) [35], we do not neglect the integral for r < r cut and take the weak cusp into account.Therefore, the gamma-ray flux from dark matter annihilation can be calculated for a specific dark matter particle, the distance to the IMBH and its corresponding dark matter spike parameters.Assuming γ sp = 7/3, we note here that the flux is effectively proportional to ⟨σv⟩ 2/7 m −9/7 χ since the cutoff radius r cut itself is proportional to the dark matter cross section and mass via r cut ∝ ⟨σv⟩ 3/7 m −9/7 χ (see Equations 2.4, 2.6, and 2.8).

EAGLE simulations
The Evolution and Assembly of GaLaxies and their Environments (EAGLE) project [52,65] is a cosmological simulation following the evolution of galaxies in a ΛCDM universe adopting the cosmological parameters advocated by the Planck Collaboration [9].The simulations are performed using a modified version of the N-Body Tree-PM Smoothed Particle Hydrodynamics (SPH) (GADGET-3) code [66].Gravitationally bound structures are identified using the Subfind algorithm [67,68].The EAGLE simulations are calibrated to reproduce the stellar mass function, galaxy sizes, and the galaxy mass-black hole mass relation at z ∼ 0 and include a variety of physical processes, such as star formation and feedback, stellar mass loss, black hole accretion and active galactic nucleus (AGN) feedback.The formation scenarios of black holes ending up in the centre of galaxies, i.e. the remnants of Population III stars, the collapse of cold gas in early-forming and massive halos, or the runaway collisions of stars and stellar mass black holes [69], cannot be resolved by the EAGLE simulations.Therefore, seed black holes of mass 10 5 M ⊙ /h are placed into the centre of halos greater than 10 10 M ⊙ /h if they do not already contain a black hole, following the approach in Springel et al. (2005) [70].The black holes can grow by accreting gas from their surroundings and by merging with other black holes.Their accretion rate is calculated using the Bondi-Hoyle-Lyttleton accretion rate [71] and modified as described by Schaye et al. (2015) [52].At each time step of the simulation, the black holes are manually repositioned by moving them to the location of the neighbouring particle with the lowest gravitational potential, which has a relative velocity less than 25 % of the sound velocity and has a distance to the black hole smaller than three gravitational softening lengths.These conditions ensure that black holes in gas-poor halos do not jump to nearby satellites.The EAGLE simulation used in this work is the reference dataset Ref-L0100N1504 performed in a periodic box with a comoving side length of L = 100 cMpc (comoving Mpc), total number of particles of N = 2 × 1504 3 , initial baryonic particle mass of m g = 1.81 × 10 6 M ⊙ and total dark matter particle mass of m χ = 9.7 × 10 6 M ⊙ .The comoving Plummer-equivalent gravitational softening length is 2.66 ckpc and the maximum physical softening length is 0.7 ckpc.The database is split into the EAGLE galaxy database containing information about halos, galaxies and their merger trees, and the EAGLE particle data containing information about each individual gas, dark matter, star and black hole particle within the simulation * .This work makes use of both the EAGLE galaxy database and the EAGLE particle data.The galaxy database is used to select Milky Way-like galaxies and to determine the formation halo of the black holes, and the particle data is used to extract information about the black holes themselves, such as their mass and their coordinates.

Milky Way-like galaxy selection
The gamma-ray flux from dark matter self-annihilation is antiproportional to the squared distance from the observer to the IMBH, as shown in Equation 2.9.Therefore, we are particularly interested in IMBHs within our own Milky Way, which thus leads to higher gamma-ray fluxes.In the following, we describe how we select galaxies with properties similar to the Milky Way within the EAGLE simulations.
We derive our selection criteria of Milky Way-like galaxies using the previous works of Callingham et al. (2019) [72], Ortega-Martinez et al. (2022) [73] and Bignone, Helmi and Tissera (2019) [74] as a starting point.Furthermore, based on Wang et al. (2020) [75], we apply an additional requirement on the halo mass M 200 defined as the mass enclosed within a sphere of radius R 200 , which is the radius at which the mean density of the halo is 200 times the critical density of the universe.Wang et al. (2020) found that the mass M 200 of the Milky Way is likely to be in the range 0.5 − 2.0 × 10 12 M ⊙ and we therefore select Milky-Way-like galaxies within EAGLE with this particular mass range.In addition, we require that galaxies have a stellar mass range M * (r < 30 kpc) of 10 10.4 -10 11.2 M ⊙ based on Ortega-Martinez et al. ( 2022) [73].The selection criteria regarding the halo mass and stellar mass are close to those used in Sanderson et al. (2020) [76] based on the FIRE-2 simulation to generate synthetic surveys resembling Gaia DR2 in data structure, magnitude limits, and observational error.Furthermore, we apply selection cuts on the current star formation rate (SFR) of the galaxy and the stellar disk-to-total mass ratio f disk .We use the stellar disk-to-total mass ratio f disk for massive EAGLE galaxies at redshift z = 0 from Proctor et al. ( 2024) [77] who applied Gaussian mixture models to the kinematics of stellar particles and identified the disk, bulge, and intra-halo light (IHL) components of EAGLE galaxies.We follow the selection cuts from Bignone, Helmi and Tissera (2019) regarding the SFR of the galaxy and the stellar disk-to-total mass ratio f disk , i.e. we require the SFR to be in the range 0.1 − 3 M ⊙ yr −1 and f disk > 0.4.Since the Milky Way has not undergone any major mergers within the past couple of billion years [78], we also require that the host halos are relaxed systems, i.e. the distance between the centre of mass and the centre of potential of the galaxy is less than 0.07R 200 and that the halo mass in the substructures is less than 10 % of the halo mass of the galaxy [72].Furthermore, we require that the satellite galaxies are located at a distance of 40 kpc < r ′ < 300 kpc, where r ′ = r(10 12 M/M 200 ) 1/3 is the distance in units of the virial radius [72].This results in satellite distribution similar to the Milky Way.The selection criteria are summarised below.We do not consider halos that have been labelled as spurious within the EAGLE database since these entities should not be considered as genuine galaxies [79].The resulting dataset consists of about 150 Milky Way-like galaxies and about 6300 associated satellites.In the following, we often refer to the central galaxy in these systems as main galaxy and their corresponding satellites as satellite galaxies.In the EAGLE simulations, the main galaxies are classified by SubGroupNumber = 0 and the satellite galaxies by SubGroupNumber > 0.

Properties of selected EAGLE galaxies
We first have a detailed look at our selection of Milky Way-like galaxies within the EAGLE simulations to ensure that the selected main galaxies meet the properties of the Milky Way.
Figure 2 shows the distribution of the main galaxy mass M 200 (left), the SFR (middle) and the galactic disk, bulge and IHL mass fractions f (right).By construction, the halo mass M 200 of the our selection of main galaxies ranges between 0.5 − 2.0 × 10 12 M ⊙ with a median halo mass M200 of 1.25 +0.31 −0.32 × 10 12 M ⊙ .Here and in the following sections, the errors on the median are calculated by determining the 16 th and 84 th percentile of the distribution in order to cover 68 % of the data around the median.The median halo mass M200 agrees very well with the halo mass from a variety of observations using Gaia DR2 data [80,81] in Wang et al. (2020) [75] which results in ∼ 1.2 × 10 12 M ⊙ .The mass distribution in Figure 2 (left) shows a distinct peak at M200 with significantly less galaxies with masses at the lower and upper end of our mass range.This highlights that the majority of our selected EAGLE galaxies are very close to the current estimate of the halo mass of the Milky Way.The SFR distribution in Figure 2 (middle) shows a moderate peak around its median star formation rate SFR of 1.53 +0. 69 −0.78 M ⊙ yr −1 .Chomiuk and Povich (2011) [82] considered different determinations of the SFR in the Milky Way and showed that the SFR converges to 1.9 ± 0.4 M ⊙ yr −1 , which lies within the error range of SFR.Lastly, the galactic disk, bulge and IHL mass fraction distributions in Figure 2 (right) show distinct peaks for each mass fraction.The median values fdisk = 0.61 +0.08 −0.10 , fbulge = 0.30 +0.07 −0.06 and fIHL = 0.08 +0.08 −0.03 indicate that our selection of main galaxies are composed of a prominent disk component with a smaller bulge and an even smaller IHL component.We note that the stellar disk-to-total mass ratio f disk of the Milky Way is measured to be around 0.86 [83] and is therefore significantly higher than f disk of our selection of EAGLE galaxies.The maximum stellar disk-to-total mass ratio f disk of our selection of EAGLE galaxies is about 0.82, meaning that none of the EAGLE galaxies reaches a disk component that is as dominant as the Milky Way disk component.Although this could potentially be a limitation of our analysis, we will show in Section 4.4 that the stellar disk-to-total mass ratio does not seem to have a significant impact on our estimate of the number of IMBHs within our Milky Way.Overall, we conclude that the key properties of our selection of Milky Way-like galaxies within the EAGLE simulations align well with those measured for the Milky Way itself.

IMBH number distribution and coordinates
Given that major mergers of black holes, i.e. mergers of black holes with similar mass, can lead to the disruption of the dark matter spike [35], we only consider unmerged IMBHs in this analysis.Furthermore, we exclude black holes with a mass m BH > 10 6 M ⊙ to stay within the mass range of IMBHs.For each Milky-Way like galaxy in the EAGLE simulations, we determine its GroupNumber and assign IMBHs with the same GroupNumber to the galaxy.In total, we find about 2500 IMBHs of which about 2000 IMBHs are associated with the main galaxies and the remaining ∼ 500 IMBHs are associated with satellite galaxies.The number distribution of IMBHs, i.e. the number of IMBHs within a galaxy N BH versus the number of galaxies N g with N BH , is shown in Figure 3 (left).We distinguish between the number distribution of all IMBHs, i.e.IMBHs associated with main or satellite galaxies (labelled as 'M+S' in the following), and the number distribution of IMBHs associated with the main galaxies only (labelled as 'M' in the following), thus excluding IMBHs associated with satellite galaxies.The median number of IMBHs ÑBH,M+S within our selection of galaxies is 15 +9 −6 , i.e. we expect about 15 IMBHs within a Milky Way-like galaxy and its corresponding satellite galaxies.Considering only IMBHs associated with the main galaxy, we find a median number of IMBHs of ÑBH,M = 12 +8 −6 .We further note that ∼ 20 % of satellite galaxies that contain at least one star particle, i.e. excluding dark matter halos, are hosting at least one IMBH.These facts make not only the Milky Way itself an interesting target for IMBHs searches but also underscores the importance of its satellite galaxies.One should notice that the nonobservation of IMBHs with masses larger than ∼ 10 5 M ⊙ in Milky Way satellite galaxies is not in tension with the above-mentioned prediction given the present associated uncertainties.However, the estimated number of IMBHs within the Milky Way differs significantly from the 101 ± 22 IMBHs found by Bertone, Zentner and Silk (2005) [35].We discuss potential causes for these differences in Section 6. Next, we identify the spatial distribution of IMBHs within their main galaxies, focusing on their radial distribution and their density distribution in galactic coordinates.For each main galaxy, we determine the coordinates of the IMBHs for a coordinate system with its origin at the galactic centre.We define the galactic centre of each galaxy as the centre of potential, which is given in the EAGLE galaxy database by the CentreOfPotential x, CentreOfPotential y and CentreOfPotential z in the SubHalo table.The coordinates of the IMBHs are given in the EAGLE particle database by the Coordinates parameter.As we are interested in the IMBH coordinates in the galactic frame, we rotate the coordinate systems so that the galaxy angular momentum vector (or spin vector), given by Spin x, Spin y and Spin z, is aligned with the z-axis of our coordinate system.This step ensures The distributions for all IMBHs, i.e.IMBHs associated with main or satellite galaxies (M+S), and for the IMBHs associated with the main galaxies only (M) is shown.
that the disk of the galaxy is located at a galactic latitude of 0 • .Afterwards, we calculate the distance of the IMBH to the galactic centre d GC and to the Sun d Sun , and the corresponding galactic coordinates, i.e. galactic latitude b and longitude l.We rescale the distance between the Galactic Centre and the Sun d ′ GS based on the halo mass M 200 of the main galaxy using d ′ GS = d GS (M 200 /10 12 M ⊙ ) 1/3 with d GS = 8.33 kpc [84].Figure 3 shows the cumulative radial distribution for all IMBHs (M+S) and for IMBHs associated with the main galaxies only (M).In both distributions, the IMBHs are concentrated towards the centre of their main galaxy.As expected, comparing the two cumulative radial distributions indicates that the IMBHs in the satellite galaxies are located at larger distances to the galactic centre.The median distance of all IMBHs to the galactic centre is dGC,M+S = 94 +124 −71 kpc and the median distance of the IMBHs associated with the main galaxy is dGC,M = 69 +122 −50 kpc.About 80% of all IMBHs and about 86% of the IMBHs associated with the main galaxies are within a distance of ∼ 200 kpc to the galactic centre.A 2D map of the positions of IMBHs (M+S) in galactic coordinates is shown in Figure 4.For the 2D map, the catalogue was upsampled by randomly generating an angle ϕ r between 0 and 2π, and adding ϕ r to the azimuthal angle of the IMBHs in the reference frame with its origin at the galactic centre.The resulting coordinates are added to the IMBH catalogue and the process is repeated 100 times.This way, we achieve an upsampling of the IMBH coordinates by a factor 100 under the assumption that the distribution of IMBHs is independent of the azimuthal angle.The azumithal angle distribution of the original IMBH catalogue was found to be uniformly distributed, which justifies our upsampling method.The probability density function (PDF) is calculated using a Gaussian kernel density estimation (KDE) with Scott's rule as bandwidth selection method [85] and the Haversine metric for distance calculation [86].The larger the PDF value in a given region, the higher the IMBH number density in that region.The contours correspond to the integral of the PDF for a given sky area such that they contain 10%, 20%, 30%, 40% and 50% of the total number of IMBHs.The PDF shows that the IMBHs are not uniformly distributed in the galaxy.Instead, they are concentrated towards the centre of the galaxy and along the galactic plane.We further discuss the consequences of this distribution for the detectability of a dark matter annihilation signal in Section 6.
If not explicitly stated otherwise, we make use of all IMBHs (M+S) for our calculations in the following Sections.Figure 4: 2D map of IMBHs associated with main or satellite galaxies (M+S) in galactic coordinates.The PDF is calculated using a gaussian kernel density estimation and the contours correspond to the integral of the PDF for a given sky area such that they contain 10%, 20%, 30%, 40% and 50% of the total number of IMBHs.Only 1 % of the upsampled IMBH coordinates are depicted here (see text for details).

IMBH mass and formation redshift
In order to calculate the dark matter spike parameters for each individual IMBH, we need to know the mass m BH and the formation redshift z f of the IMBH.We extract the mass and formation time of the IMBHs at redshift z = 0 in the EAGLE particle database from the BH Mass and BH FormationTime parameters.Figure 5 shows the mass and formation redshift distribution of IMBHs within our selection of Milky Way-like galaxies.They are determined by extracting the distribution for each individual galaxy and then calculating the mean of the distributions per bin.Therefore, these distributions represent the average distributions one would expect for a Milky Way-like galaxy.The large majority of IMBHs have a mass close to the initial seeding mass of 10 5 M ⊙ /h = 1.48 × 10 5 M ⊙ .The median mass of IMBHs is 1.49 +0.14 −0.02 • 10 5 M ⊙ and the distribution rapidly decreases for increasing black hole mass.Since we excluded black holes that encountered any merger during their lifetime, the mass accretion of IMBHs is purely caused by the accretion of gas from their surroundings.The formation redshift distribution of IMBHs is shown in Figure 5 (right).The median formation redshift of IMBHs is 2.78 +2.06 −2.01 with the largest number of IMBHs having formed at a low redshift with z ≲ 1.This meets our expectation since the seed black holes are placed into the centre of halos with a total mass larger than 10 10 M ⊙ /h.The number of halos with a mass larger than 10 10 M ⊙ /h increases with increasing evolution time of the universe, resulting in a higher number of IMBHs forming at a lower redshift.Note that the low formation redshift and the small growth of the black holes is a consequence of the limited resolution of the EAGLE simulation.As briefly discussed in Section 1, IMBHs are expected to form via processes that take place at much higher redshifts, i.e z ≳ 10 [37][38][39]41].However, due to the black hole seeding mechanism applied in EAGLE, most of the black holes in the simulation appear at z ≲ 5. We therefore assume here that the black holes actually formed at higher redshifts and then grew until the seeding mass of 10 5 M ⊙ /h was reached.Thus, the majority of the growth of the black holes and consequently the formation of the dark matter spikes takes place before the black holes have reached the EAGLE seeding mass.

Correlation between N BH and galaxy properties
We calculate the correlation between the properties of the main galaxies and the number of IMBHs N BH in each galaxy in order to investigate potential indicators for the presence of IMBHs.We calculate the correlation coefficient c i between the number of IMBHs N BH and the galaxy properties as follows: whereas cov(X, Y ) is the covariance of X and Y , σ X is the standard deviation of X and i = {M 200 , SFR, f disk }. Figure 6 shows the number of all IMBH N BH in each galaxy versus the mass of the galaxy M 200 (left), the SFR (middle) and the stellar disk-to-total mass ratio f disk (right).The strongest correlation is observed between the number of IMBHs N BH and the galaxy mass M 200 with a correlation coefficient of c M 200 = 0.51.Whereas galaxies with a mass M 200 ≲ 8 × 10 11 M ⊙ contain no more than ∼ 20 IMBHs, galaxies with M 200 ≳ 1.2 × 10 12 M ⊙ can contain up to ∼ 40 IMBHs.Both the star formation rate and the stellar disk-to-total mass ratio f disk show only a very mild correlation with the number of IMBHs N BH with c SFR = 0.11 and c f disk = −0.07,respectively.We discuss the consequences of these results in more detail in Section 6.

Dark matter spike parameters
For each IMBH, we calculate the dark matter spike parameters r sp , ρ(r sp ) and r cut .Therefore, it is necessary to determine the host halo in which the IMBH formed.We extract the formation redshift z f of each IMBH at z = 0 as described in the previous section and determine the closest redshift z c with z c ≤ z f for which a snapshot in the EAGLE simulations is available.In this snapshot, we identify the IMBH based on its ParticleIDs and identify its host galaxy based on its GroupNumber and SubGroupNumber.We extract the dark matter density profile ρ host (r) of the host galaxy using the mass profile M (< r) information, given by the ApertureSize and Mass DM parameters from the Aperture table of the EAGLE galaxy database.We fit the dark matter density profile to the NFW profile as described in Equation 2.2 and obtain the normalisation constant ρ 0 and the scale radius r s using the least squares method [87].We use ρ 0 and r s to calculate the radius of gravitational influence r h using Equation 2.3.We apply the IMBH mass at z = 0 for Equation 2.3, assuming an adiabatic growth of the IMBH since its formation time.Finally, the spike radius r sp is calculated by r sp = 0.2r h [63] and the spike density ρ(r sp ) by evaluating the NFW profile at r sp .In some cases, the IMBH is not assigned to any halo at z c .In this case, we use the dark matter density profile of the closest halo at z c with a mass larger than the required mass to form a IMBH, i.e.M > 10 10 M ⊙ /h, to calculate the spike parameters.Therefore, we assume that the IMBH has formed in the next closest halo with sufficient mass.Furthermore, IMBHs assigned to a spurious halo at z c or containing zero Mass DM values from the Aperture table are not considered.The spike radius and spike density distribution are shown in Figure 7.The median spike radius is 3.95 +1.81 −1.27 pc and the median spike density is 1.19 +2.57−0.80 • 10 3 GeV cm −3 .

Results
In this section, we present our results for the detectability of a gamma-ray signal from dark matter annihilation around IMBHs in the Milky Way.Assuming a dark matter mass m χ and cross section ⟨σv⟩, we calculate the cutoff radius r cut using Equation 2.8. Figure 8 shows the distribution of the cutoff radius for m χ = 500 GeV and ⟨σv⟩ = 3 × 10 −26 cm 3 s −1 , assuming the bb−annihilation channel.The median cutoff radius is 2.13 +0.35 −0.43 • 10 −3 pc.The larger the assumed dark matter cross section ⟨σv⟩, the larger the cutoff radius as more self-annihilation events occur and the saturation of the dark matter spike is reached at larger radii.The expected gamma-ray flux from dark matter annihilation is calculated by implementing the distance D to the IMBH, the dark matter mass m χ , cross section ⟨σv⟩ and spike parameters r sp , ρ(r sp ) and r cut into Equation 2.9.We compute the integrated luminosity of IMBHs by calculating the number of IMBHs N BH (Φ) that surpass a certain flux threshold Φ(E > E th ), assuming a typical energy threshold of (a) E th = 100 GeV for Imaging Atmospheric Cherenkov Telescopes (IACTs) and (b) E th = 100 MeV for space-based gamma-ray observatories.Figure 9 (left) shows the average integrated luminosity over all Milky Waylike galaxies for dark matter masses between 0.5 TeV and 1.5 TeV, fixed annihilation cross section of ⟨σv⟩ = 3 × 10 −26 cm 3 s −1 , the bb-channel and E th = 100 GeV.For the range of dark matter masses and energy threshold chosen here, the integrated luminosity increases for increasing dark matter mass.This is due to the dark matter annihilation spectrum being integrated from E th to m χ , resulting in an increasing number of photons for higher dark matter mass.The average H.E.S.S. flux sensitivity for the H.E.S.S. galactic plane survey is ∼ 5 × 10 −13 cm −2 s −1 [88] and depicted by the grey vertical line in Figure 9. Independent of the dark matter masses considered here, all IMBH fluxes are expected to surpass the H.E.S.S. sensitivity, making IMBHs promising targets for ground-based gamma-ray observatories.In order to test the limits of our analysis, we lowered the dark matter cross section until only ∼ 2.3 IMBHs would exceed the H.E.S.S. sensitivity.This number is motivated by the 90 % confidence level of a non-detection from Poisson statistics.We find that ∼ 2.3 IMBHs exceed the H.E.S.S. sensitivity for a dark matter mass of m χ = 500 GeV at a cross section of ∼ ⟨σv⟩ = 7 × 10 −37 cm 3 s −1 , see Figure 9 (right).However, this cross section cannot be directly translated to an upper cross section limit that H.E.S.S. would be able to probe because (a) H.E.S.S. does not have a full sky coverage and (b) H.E.S.S. does not reach the flux sensitivity of the galactic plane survey in all of its observations.We discuss these limitations in more detail in the next section. .The Fermi-LAT flux sensitivity at the Galactic Centre is lower than at higher galactic latitudes and longitudes due to the higher gamma-ray background from the diffuse galactic emission.
All IMBH fluxes are expected to surpass the Fermi-LAT sensitivity at large galactic longitudes and latitudes independent of the dark masses chosen here.Additionally, the majority of IMBH fluxes are expected to surpass the Fermi-LAT Galactic Centre sensitivity.The fact that Fermi-LAT's flux sensitivity is high enough for both small and large galactic coordinates to detect the majority of IMBHs indicates that the instrument has the potential in detecting a gamma-ray signal from dark matter annihilation around IMBHs independent of the IMBH sky position.Similarly to Figure 9 (right), we lower the dark matter cross section until ∼ 2.3 IMBHs would exceed the Fermi-LAT flux sensitivity at the Galactic Centre.We find this to be the case at ⟨σv⟩ ∼ 10 −32 cm 3 s −1 .The corresponding integrated luminosity is shown in Figure 10 (right).Unlike ground-based gamma-ray observatories, Fermi-LAT provides data for the full gamma-ray sky and is therefore not limited to a specific sky region, such as the Galactic Plane.We discuss the consequences of this fact in more detail in the next section.Furthermore, we investigate the impact of different dark matter density profiles on our results.In addition to the NFW profile, we assume that the dark matter density of the IMBH formation halos follows a cored density profile.Since we find that the EAGLE data does not provide a sufficient resolution in order to properly test the cored profile, we refer the interested reader to the results in Appendix A. A comparison of different dark matter annihilation channels and their impact on the expected gammm-ray flux around IMBHs is discussed in Appendix B.  6 Discussion

Number of IMBHs
Our analysis provides the number and distribution of IMBHs in a Milky Way-like galaxy under the assumption that a 10 10 M ⊙ /h halo is populated by one IMBH with 10 5 M ⊙ /h that subsequently grows through accretion.This black hole formation scenario in the EAGLE simulations is motivated by two key factors: firstly, the resolution of the EAGLE simulations is insufficient to accurately represent the actual formation processes of black holes.Secondly, the feedback resulting from the growth of these black holes plays a pivotal role in the formation of galaxies, influencing star formation in massive galaxies and altering the gas profiles of their host halos [52,65].In this minimal scenario for IMBH formation, the average number of IMBHs within a Milky Way-like galaxy and its corresponding satellite galaxies is 15 +9 −6 , which is significantly lower than the 101 ± 22 IMBHs found by Bertone, Zentner and Silk (2005) [35].This difference is likely caused by the different black hole seeding mechanisms considered in our analysis.The EAGLE simulations seed a black hole in halos with masses greater than 10 10 M ⊙ /h.On the other hand, the work of Bertone, Zentner and Silk (2005) follows the procedure of Koushiappas, Bullock, and Dekel (2004) [90], in which the seeding of black holes is based on a critical halo mass as a function of, among others, the redshift and gas temperature.This leads to host halo masses down to ∼ 10 7 M ⊙ , which allows black holes to form in halos up to three orders of magnitudes smaller than in our analysis.The effect is indirectly illustrated in the redshift distribution shown in Figure 5 (right).Whereas the majority of the black holes in our analysis are seeded at low redshifts with z ≲ 5, the black holes in Bertone, Zentner and Silk (2005) are seeded at redshift z ≳ 10.This is due the fact that the halos acquire more mass over time and reach the required (lower) seeding mass sooner in the work of Bertone, Zentner and Silk (2005).Despite these technical arguments, the approach we employ with the EAGLE simulations offers a more timely understanding of the formation and evolution of galaxies and is validated against the latest observations.While the Bertone, Zentner and Silk (2005) includes smaller progenitor halos in their semi-analytical models, these models lacked the dynamic environmental effects and feedback mechanisms that are now known to significantly impact galaxy formation and black hole growth.EAGLE, on the other hand, represent state-of-the-art cosmological simulations that have been calibrated against a range of observables in our Universe.These simulations include, among others, detailed modelling of star formation, stellar evolution, metal enrichment, supernova feedback, and AGN feedback.These processes are crucial as they influence the thermodynamic properties of the interstellar and intergalactic medium, and hence the formation and evolution of galaxies and black holes within these environments.Finally, while we are confident that our results represent a more robust estimate based on our current understanding of galaxy evolution, we acknowledge that they are not the only possible predictions.The seeding mechanism used in EAGLE is supported by a wide range of observations [65] but does not exclude other plausible scenarios or models.Furthermore, it is important to note that the number of IMBHs N BH within our selection of Milky Way-like galaxies scatters significantly between the individual galaxies, varying from a minimum of 1 IMBH per galaxy up to a maximum of 43 IMBHs per galaxy.This number strongly correlates with the galaxy mass M 200 as can be seen in Figure 6 (left).Other massrelated parameters, such as the total or star mass of the galaxy show a similar correlation with the number of IMBHs.A precise measurement of M 200 of the Milky Way is rather challenging and different approaches can lead to significantly different results, see Wang et al. (2020) [75] for a detailed comparison.Our choice of the M 200 -range to select Milky Way-like galaxies within EAGLE is well motivated and aligns with the range of the actual M 200 measurements of the Milky Way but future, more precise constrains of the Milky Way mass will improve our predictions of the number of IMBHs and make them more robust.The actual stellar disk-to-total mass ratio f disk of the Milky Way is measured to be around 0.86, which is significantly higher than f disk of most of our selected EAGLE galaxies.However, as shown in Figure 6 (right), the correlation between f disk and the number of IMBHs is very mild with a correlation coefficient of c f disk = −0.07.The number of IMBHs within the Milky Way predicted by our analysis is therefore expected to be only mildly affected by the difference between f disk of our selected EAGLE galaxies and the true f disk value of the Milky Way.Lastly, we do not find a strong correlation between N BH and the SFR with a correlation coefficient of c SFR = 0.11 as shown in Figure 6 (middle).This indicates that the presence of IMBHs barely promotes the total SFR of our selection of EAGLE galaxies.

Spatial distribution of IMBHs
We find that the IMBHs are not uniformly distributed in the galaxy, but that they are concentrated towards the centre of the galaxy and along the Galactic Plane.We suppose that the IMBHs do not follow a uniform distribution due to the manual repositioning of black holes applied in the EAGLE simulations.Although the positions and trajectories of massive black holes are affected by dynamical friction in reality, EAGLE lacks the resolution to capture this effect.Instead, the effect of dynamical friction is modelled by manually repositioning black holes to the location of the neighbouring particle with the lowest gravitational potential at each time step of the simulation [52], as already discussed in Section 3.1.Without repositioning, black hole growth is negligible and SMBHs do not end up in the centre of massive galaxies [91] which would be in strong contradiction with observations.The manual repositioning of black holes applied in the EAGLE simulations therefore seems to capture the effect of dynamical friction reasonably well and is considered to be the main explanation for the IMBH distribution along the Galactic Plane.For a detailed study of the effect of black hole repositioning in galaxy formation simulations, we refer the reader to Bahé et al. (2022) [91].
Additionally, we suppose that the rather late seeding of the black holes at low redshifts is further contributing to the spatial distribution of the black holes along the Galactic Plane.Due to their injection at low redshifts, the black holes are more prone to interaction with the Galactic Plane which therefore makes them more spatially correlated to the baryonic matter in the Galactic Plane.Furthermore, we observe that the majority of the IMBHs are located within the main galaxy, although a considerable number is also found in the associated satellite galaxies.On average 15 +9 −6 IMBHs are distributed within the main galaxy and its corresponding satellites, and 12 +8 −6 IMBHs are present in the main galaxy only, indicating that 3 +2 −2 of the 15 +9 −6 IMBHs are distributed within satellite galaxies.It is also interesting to note that, within our selection of Milky Way-like galaxies, about 20 % of the satellite galaxies with at least one star particle contain at least one IMBH.For the Milky Way, more than 60 satellite galaxies within ∼ 400 kpc have been observed so far [92].This makes not only the Milky Way itself but also its satellite galaxies an interesting target for IMBH searches.Given the distribution of IMBHs within both the main galaxy and its satellite galaxies, ground-based gamma-ray observatories, which do not have full sky coverage, should therefore focus on observations of the Galactic Centre, the Galactic Plane, and satellite galaxies to maximise the number of IMBHs within their field of view.Fortunately, all current generation ground-based gamma-ray observatories, i.e.H.E.S.S., MAGIC and VERITAS, have a Galactic Centre survey and observations of many satellite galaxies in their program [13,17,88,[93][94][95][96].Due to the high gamma-ray fluxes that are expected from the dark matter annihilation around IMBHs, we do not expect the flux sensitivity to be the limiting factor for the detectibility of a potential gamma-ray signal with ground-based observatories.Instead, the low expected number of IMBHs within the Milky Way is likely to limit the detectibility of a potential gamma-ray signal.Integrating the PDF from Figure 4 over the sky region of the H.E.S.S. galactic plane survey [88], i.e. for −65 • < l < 110 • and |b| < 3 • , and excluding the IMBHs within satellite galaxies, we expect about N BH,HESS = 0.6 +0. 4  −0.3IMBHs within the field of view.The upcoming Cherenkov Telescope Array Observatory (CTAO) [97] is planning to observe the galactic plane covering a larger sky area with |l| < 90 • and |b| < 6 • , resulting in an expected number of IMBHs within the field of view of N BH,CTAO = 1.1 +0.8  −0.6 .The H.E.S.S. extragalactic survey (HEGS) [98] covering a set of extragalactic observations will increase the sky area covered by H.E.S.S. and therefore improve the chances of detecting an IMBH signal although the average flux sensitivity is not expected to be as high as for the galactic plane survey.Other gamma-ray observatories, such as HAWC [99], LHAASO [100] and Fermi-LAT [101], are able to observe (almost) the full sky and are therefore very well suited for the detection of a gamma-ray signal from dark matter annihilation around IMBHs.

Gamma-ray flux detectability of IMBHs
The integrated luminosity functions of IMBHs shown in Figure 9 and Figure 10 indicate that a potential gamma-ray signal from dark matter annihilation around IMBHs is detectable with current and future gamma-ray observatories.These objects are expected to appear as unidentified, point-like sources with identical energy spectra.Depending on the experiment, such analyses should be able to probe a wide range of dark matter masses, ranging from sub-GeV to tens of TeV, and cross sections down to ⟨σv⟩ ∼ 7 × 10 −37 cm 3 s −1 in the most optimistic scenario.To our knowledge, these limits would be the most stringent limits on dark matter annihilation cross sections for dark matter masses in the range of ∼ GeV − TeV.However, we emphasize that the search for dark matter spikes around IMBHs is affected by significant uncertainties.The predicted number and spatial distribution of IMBHs in the Milky Way depends on the assumed formation scenario and seeding mechanism applied within a simulation [49].Furthermore, our limits on the dark matter cross section depend on the expected number of IMBHs within a galaxy, which is strongly correlated to the mass M 200 of the galaxy.This introduces an additional systematic uncertainty that can only be reduced by determining the mass of the Milky Way more precisely in the future.Moreover, the IMBH formation scenarios have an impact on the dark matter distribution around the IMBHs at z = 0.The formation of IMBHs at high redshift is still subject of current research and the details remain to be fully understood [36].The dark matter cross section an analysis is able to probe is directly influenced by those uncertainties.Deriving strict limits associated with high confidence levels is therefore challenging and the limits provided in this article should be treated with caution.However, by integrating state-of-the-art cosmological simulations and implementing recent measurements of the Milky Way, we argue that our analysis offers a novel and greatly improved approach on identifying dark matter spikes around IMBHs in comparison to previous studies.Here, we have covered the case of dark matter annihilation into gamma rays, but the analysis can be easily extended to other indirect detection channels, such as the detection of neutrinos [102] with the IceCube [103] and KM3NeT [104] experiments, using the IMBH catalogue and the source code provided in this work.Furthermore, IMBHs could potentially emit X-ray and radio emissions from the accretion of matter in the accretion disk.Gaggero et al. (2017) [105] and Scarcella et al. (2021) [106] investigated the multi-wavelength detectability of primordial and astrophysical black holes, respectively.They found that black holes concentrated in the denser central regions of the Galaxy are more likely to accrete gas and produce detectable emissions compared to those located in the less dense outer regions.However, the high gas density region along the Galactic Disk spans only a few ∼ 0.1 deg in Galactic latitude.Although we find that the IMBHs in our analysis are mainly distributed towards the Galactic Centre and along the Galactic Disk, the concentration of these objects within |b| ≲ 0.3 deg is very low.It is therefore unlikely to find an IMBH within this region and we do not expect the IMBHs from our analysis to be strong multiwavelength emitters.We postpone a more detailed analysis of the multiwavelength detectability of these objects to future work.Our IMBH catalogue and the catalogue of our selection of Milky Way-like galaxies within EAGLE is publicly available at [59].The galaxy catalogue contains, among others, the mass parameters, the SFR and the stellar disk-to-total mass ratio of each individual galaxy.The IMBH catalogue contains, among others, the coordinates, mass, formation redshift and spike parameters for each individual IMBH.Each column of the catalogues is described in detail in Table 2 & 3. We also provide separate files for which we calculated the gamma-ray fluxes for different dark matter masses m χ and annihilation cross sections ⟨σv⟩.The columns of these files are described in Table 4.The source code used to generate the IMBH catalogue and the gamma-ray fluxes is publicly available at [58].It provides a detailed description of the analysis steps and can be used to generate the IMBH catalogue for different EAGLE simulations.

Conclusions & Future Work
We presented a mock catalogue of IMBHs and their dark matter spikes in Milky Way-like galaxies derived from the EAGLE simulations.The catalogue contains the coordinates, mass, formation redshift and dark matter spike parameters for each individual IMBH.On average, our selection of Milky Way-like galaxies contains 15 +9 −6 IMBHs, primarily distributed towards the Galactic Centre and along the Galactic Plane.We demonstrated that the gamma-ray flux from dark matter annihilation around IMBHs should be detectable by both current and forthcoming gamma-ray observatories, including H.E.S.S, Fermi-LAT and the upcoming Cherenkov Telescope Array Observatory (CTAO).Depending on the experiment, such analyses should be able to probe a wide range of dark matter masses, ranging from sub-GeV to tens of TeV, and cross sections down to ⟨σv⟩ ∼ 7 × 10 −37 cm 3 s −1 in the most optimistic scenario.To the best of our knowledge, these limits would be the most stringent constraints on dark matter annihilation cross sections for dark matter masses in the range of ∼ GeV − TeV.However, it is crucial to consider the limitations of the EAGLE simulations in capturing the complex processes involved in IMBH formation due to their spatial and temporal resolutions, as highlighted by recent studies of IMBH formation [54][55][56][57].These limitations introduce systematic uncertainties to our analysis that may have a substantial impact on our results.The IMBH catalogue and the source code used to generate the catalogue and calculate the gamma-ray fluxes are publicly available at [58] and [59].The source code provides a detailed description of the analysis steps and can be used to generate the IMBH catalogue for different EAGLE simulations.In future work, we aim to apply to the IMBH mock catalogue on data of the H.E.S.S. Galactic Plane survey, the H.E.S.S. extragalactic survey and on observations of satellite galaxies to search for a potential gamma-ray signal from dark matter annihilation around IMBHs by investigating the unidentified point sources within these observations.If none of the sources matches the expected dark matter annihilation spectrum, we will provide upper limits on the dark matter cross section.number of IMBHs that can be detected by ground-based gamma-ray observatories is limited by the number of IMBHs within the field of view of the experiment, as we discuss in more detail in Section 6.

B Dark matter annihilation channels
In the main text we considered dark matter annihilation to occur via the bb-and τ − τ +channel.In the following, we briefly investigate the effect of different annihilation channels on the gamma-ray flux from dark matter self-annihilation around IMBHs. Figure 12 (left) shows the gamma-ray spectra per dark matter annihilation for m χ = 500 GeV and the bb-, τ − τ + -, W − W + -and ZZ-channels.Generally, the bb-, W − W + -and ZZ-channels result in fairly similar gamma-ray spectra with the most significant differences in the high energy regime (≳ 100 GeV).All spectra decrease with increasing energy.The τ − τ + -channel leads to a lower gamma-ray emission compared to the other channels for E ≲ 60 GeV and to a significantly higher gamma-ray emission for E ≳ 60 GeV. Figure 12 (right) shows the integrated luminosity of IMBHs for E th = 100 GeV, ⟨σv⟩ = 3 × 10 −26 cm 3 s −1 and m χ = 500 GeV for the same four annihilation channels as in Figure 12 (left).The τ − τ + -channel results in the highest gamma-ray fluxes, followed by the W − W + -, ZZ-and bb-channel.This meets our expectation since we consider a energy threshold of 100 GeV and the τ − τ + -channel leads to the highest number of gamma-ray photons per annihilation for E ≳ 60 GeV.Again, the bb-, W − W + -and ZZ-channels lead to fairly similar integrated luminosities as expected from the gamma-ray spectra.Considering a fixed flux threshold of Φ(E > 100 GeV) = C Galaxy catalogue of selected Milky Way-like galaxies

10 − 7 10 −5 10 −
i) Halo mass M 200 range: 0.5 − 2.0 × 10 12 M ⊙ ii) Stellar mass M * (r < 30 kpc) range: 10 10.4 -10 11.2 M ⊙ iii) Star formation rate range: 0.1-3 M ⊙ yr −1 iv) Stellar disk-to-total mass ratio f disk larger than 0.4 v) Distance between the centre of mass and the centre of potential of the galaxy is less than 0.07R 200 vi) Total mass in substructures is less than 10 % of the total mass of the galaxy Selection criteria for satellite galaxies: i) Distance from halo centre r in the range: 40 kpc < r ′ < 300 kpc with r ′ = r(10 12 M ⊙ /M 200 ) 1/3

Figure 2 :
Figure 2: Number of our selection of Milky Way-like galaxies within the EAGLE simulations N g versus the galaxy mass M 200 (left), the star formation rate (SFR) (middle) and the galactic disk, bulge and intra-halo light (IHL) mass fractions f (right).The red lines correspond to the median values and the red shaded regions to the errors on the median, which are calculated by determining the 16 th and 84 th percentile of the distribution.

Figure 3 :
Figure 3: Number distribution (left) and cumulative radial distribution (right) of IMBHs.The distributions for all IMBHs, i.e.IMBHs associated with main or satellite galaxies (M+S), and for the IMBHs associated with the main galaxies only (M) is shown.

Figure 6 :
Figure 6: Number of IMBHs N BH in each galaxy versus the mass of the galaxy M 200 (left), the SFR (middle) and the stellar disk-to-total mass ratio f disk (right).

8 Φ 3 Figure 11 :
Figure11: Left: Example of a dark matter halo profile obtained from the EAGLE simulations.The red line represents the best fit model for the NFW profile and the yellow line the best fit model for the cored profile in which the core index γ c is a free parameter.The inset plot shows a zoom of the ∼ 1 kpc region.The bottom panel shows the residuals R for both models.Right: Integrated luminosity of IMBHs obtained under the assumption of a NFW profile, a cored profile with γ c as free parameter, a cored profile with fixed γ c = 0.9 and a cored profile with fixed γ c = 0.3.A dark matter particle with ⟨σv⟩ = 3 × 10 −26 cm 3 s −1 and m χ = 500 GeV, and the bb-channel was assumed.The grey dashed line depictes the average H.E.S.S. flux sensitivity for the H.E.S.S. Galactic Plane survey[88].

Table 1 :
[77]mn descriptions for the galaxy catalogue of selected Milky Way-like galaxies Integer identifier of this galaxy within its FoF halo at z = 0.The condition subgroup number = 0 selects central galaxies.Star formation rate of the galaxy fdisk -Stellar disk-to-total mass ratio of the galaxy (values from[77]) [77]ge-Stellar bulge-to-total mass ratio of the galaxy (values from[77]) fihl -Stellar IHL-to-total mass ratio of the galaxy (values from[77])n sat -Number of satellite galaxies associated with the galaxy n sat stars -Number of satellite galaxies with at least one star particle associated with the galaxy