Constraints on dark matter annihilation from the Event Horizon Telescope Observations of M87$^\star$

The fast developments of radio astronomy open a new window to explore the properties of Dark Matter (DM). The recent direct imaging of the supermassive black hole (SMBH) at the center of M87 radio galaxy by the Event Horizon Telescope (EHT) collaboration is expected to be very useful to search for possible new physics. In this work, we illustrate that such results can be used to detect the possible synchrotron radiation signature produced by DM annihilation from the innermost region of the SMBH. Assuming the existence of a spike DM density profile, we obtain the flux density due to DM annihilation induced electrons and positrons, and derive new limits on the DM annihilation cross section via the comparison with the EHT integral flux density at 230 GHz. Our results show that the parameter space can be probed by the EHT observations is largely complementary to other experiments. For DM with typical mass regions of being weakly interacting massive particles, the annihilation cross section several orders of magnitude below the thermal production level can be excluded by the EHT observations under the density spike assumption. Future EHT observations may further improve the sensitivity on the DM searches, and may also provide a unique opportunity to test the interplay between DM and the SMBH.


Introduction
The nature of Dark Matter (DM) is one of the biggest problems in modern physics and cosmology [1][2][3]. Quite a lot of particle models beyond the standard model were proposed in literature [4][5][6][7]. Some classical candidates are well-motivated, such as weakly interacting massive particles (WIMPs) [7,8], axions or axion-like particles [9][10][11][12], and dark photons [13,14]. Meanwhile, many detecting methods and technologies have also been developed, such as detecting the feasible energy missing in particle collisions [15], observing the scattering signals between DM and detector particles [8,16], and measuring the DM induced products in cosmic rays [17][18][19][20]. Among these candidates, WIMPs with typical masses from sub-GeV to multi-TeV are most widely studied, due partly to the so-called "WIMP miracle" the corect relic density of DM can be naturally obtained in this model. WIMPs can undergo annihilation or decay, yielding standard model elementary particles by multi-channels, such as electrons and positrons. The electrons and positrons can imprint in the electromagnetic spectrum due to radiative processes such as synchrotron, inverse Compton scattering, and/or bremsstrahlung radiation [21][22][23]. There have been considerable efforts to study multi-wavelength electromagnetic emission from DM annihilation in a variety of galaxy clusters and dwarf galaxies, e.g., Refs. [24][25][26][27][28][29][30][31][32][33]. A sizable parameter space of WIMPs was excluded by these observations. The center of a galaxy is expected to be a site with potentially maximum DM concentration, although the detailed DM density profile remains uncertain in the galaxy center due to the effect from the baryon interaction and/or the accretion of the SMBH [34][35][36][37][38]. The previous observations with poor spatial resolutions are difficult to probe the environment parameters close to the black hole, and hence increase the uncertainties of the signal prediction. With the development of technologies of the Very Long Baseline Interferometry (VLBI) [39], the Event Horizon Telescope (EHT) program, a global network of millimeter and sub-millimeter observational facilities, has been established successfully and has imaged the SMBH M87 to the event horizon scale with extremely high angular resolution [40][41][42][43]. Besides crucial tests of the physical laws in the strong gravity field, the EHT observations open a new window to study broad types of new physics models including the particle nature of DM [11,44,45]. The EHT observations enable one to scrutinize the physical and astrophysical processes just surrounding the SMBH, and can thus provide a unique probe of the DM interplay with the black hole.
In this work, we search for radio emission signature induced by the DM annihilation at subparsec scales using the EHT M87 results [40,41]. For the DM density distribution in the center of the galaxy, a spike profile is expected due to the adiabatic accretion of the SMBH [36,37,46], and the number of electrons and positrons annihilated by DM would be enhanced several orders. We will discuss the consequence of the possible radio emission via the synchrotron radiation by high-energy electrons and positrons in the magnetized accretion disk of M87 . Particularly, we find new properties of the relation between the WIMP annihilation cross section and the synchrotron fluxes, due to the annihilation saturation of the density near the horizon. Consequently a finite range on the cross section and mass plane is able to be probed by the data.
The rest of this paper is organized as follows. In Sec. 2, we present the spatial distribution of DM around the SMBH, especially the DM spike profile. In Sec. 3, we review the propagation of electrons and positrons with the synchrotron cooling and advection effect, and then derive the limits on DM parameters using the EHT data. We conclude our work with some discussion in Sec. 4.

Density profile of DM halo with SMBH
We assume that the global DM halo density profile of M87 galaxy is described by the Navarro-Frenk-White (NFW) profile, which is expected to be a universal density profile for cold DM halo in hierarchically clustering Universes [47,48]. It is parameterized as where r 0 is the scale radius and ρ 0 is the normalization constant. Eq. (2.1) shows a r −1 density cusp at the center. However, as proposed by Gondolo and Silk [36], the adiabatic accretion of DM onto the SMBH in the galactic center would significantly enhance the DM density around it and form a spiky structure in the DM density distribution. Further studies showed that some dynamical effects, such as the mergers of black holes and the interactions with the baryonic matter, may shallow or even disrupt the spike [37,49,50]. Nevertheless, we assume that the DM spike had formed around M87 and was survived from such dynamical effects. In this scenario, the spike density profile reads [49]: where δ is the slope of the original density profile, and δ sp = (9 − 2δ)/(4 − δ). We take δ = 1 to match the NFW profile. R sp is the spike radius, which is given as (2.3) In the above equation, M BH is the mass of black hole, α δ is a normalization coefficient which characterize the mass ratio of DM spike and black hole. Noticed that the depletion of DM density in the inner region of spike due to DM annihilations set a saturate DM density ρ sat = m χ / σv t BH . where t BH is the age of the black hole, and m χ and σv are repectively the mass and thermally averaged annihilation cross section of DM particle.
Using ∂ρ χ /∂t = − σv ρ 2 χ /m χ , one obtains spike density within the region R Sch ≤ r < R sp [36]: where R Sch is the Schwarzschild radius. Based on above analysis, when taking into account the effects of accretion of SMBH and DM annihilations, one should modify the NFW profile to the following piecewise distribution [51,52]: (2.5) According to the EHT results, the mass of M87 is 6.4×10 9 M , and the corresponding Schwarzschild radius is R Sch = 5.9 × 10 −4 pc. Taking α δ = 0.1, we fix r 0 = 20 kpc for the halo (as for the Milky Way), and normalize ρ 0 ≈ 2.5 GeV cm −3 as Ref. [52] done. Then, we could calculate R sp 220 pc throught Eq.(2.3). Assuming the age of M87 is t BH = 10 9 yr and using Eq. (2.5), we displayed the DM density profile of M87 in Fig. 1. In this plot, the corresponding Schwarzschild radius of M87 and the areas covered by EHT observations are also shown.

Propagation of the Electrons and Positrons
The distribution of electrons and positrons from the DM annihilation can been obtained via solving the propagation equation in the presence of synchrotron radiation and advection, which reads [51-53] where f e (r, p) is the equilibrium distribution function of electrons and positrons in momentum space at radius r and momentum p, assuming a steady-state. The first term on the left hand side describes the spatial diffusion, with D(r, p) being the diffusion coefficient. The second and third terms are the advection current of accretion flow and the energy gain of electrons and positrons caused by adiabatic compression, with v(r) = −c R Sch /r is the radial infall velocity of electrons and positrons in the accretion flow. Finally, the last term describe the energy loss due to radiative processes. We here consider the stationary state solution of the transport equation (3.1) in the spherical symmetry geometry. The source function q(r, p) is the DM annihilation rate in momentum space, which is related to the annihilation rate Q(r, E) in energy space through following equation: where σv is the thermally averaged annihilation cross-section, and dN inj e ± ,i /dE is the e ± injection spectrum through annihilation channels i with branching ratios BR i . In this paper, we respectively consider annihilation channels i = e + e − , µ + µ − , τ + τ − , bb with 100% branching ratios, and extract corresponding e ± injection spectrum from the DarkSUSY or PPPC packages [54]. The electron and positron energy spectrum is then written as In order to appropriately describe the propagation of electrons/positrons produced by DM annihilations with taking into account the influence of SMBH, one need to consider two propagation zones separately. Firstly, for r outside the accretion radius r acc , the only energy loss mechanism for electron/positrons are radiative processes discussed above. While for r ≤ r acc , one has more complicate picture due to the existence of accretion flow towards central SMBH. As a consequence, at a given injection radius R inj , corresponding energy transfer of electron/positrons is governed by the competition between two physical processes, i.e., energy loss due to conventional radiative processes and the energy gain due to adiabatic compression along the plasma accretion flow. For the ultra-relativistic electrons/positrons with their radiative energy loss is dominated by synchrotron emission, transport equation (3.1) yields an integral analytic solution [51]: where we choose accretion radius r acc = 2GM/v 2 f ∼ 0.04 pc for the typical galactic wind velocity v f 500 − 700 km s −1 . p inj ≡ p inj (R inj ; r, p) is the injection momentum of an electron/positron is injected at R inj (≥ r) and arriving at r with momentum p. The adiabatic compression in the advection process leads to velocity field of the accretion flow, with the momentum gain rateṗ ad . The associated characteristic curves of p inj (R inj ; r, p) are obtained by solving following differential equation [51,52]: dp dr =ṗ ad (r, p) +ṗ(r, p) v(r) ṗ ad (r, p) +ṗ syn (r, p) v(r) . (3.5) In above equation, approximation is due to the fact that the total radiative momentum lossṗ(r, p), is dominated by synchrotron radiation, anḋ Ep m 2 e c 3 . (3.6) Following method introduced in Refs. [51,52], for a homogeneous magnetic field in the galactic halo, the analytical solution of p inj (R inj ; r, p) for Eq. (3.5) is given by It should be emphasize that the denominator in Eq. (3.7) may vanish and even become negative, thus leads to nonphysical values for injection momentum. To be specific, this depends on the efficiency of the accretion flow and characterizes the region of the injection parameters (R inj , p inj ) for a given arrival point (r, p). In practice, p inj holds positive for We then use the value of min[r acc , R 0 inj ] as an effective upper bound for the integration in Eq. (3.4). While for lower bound for the integration, we take r = r ISCO with r ISCO being the radius of the Innermost Stable Circular Orbit (ISCO), which indicates the innermost radius to maintain the stable orbit in the equatorial plane for a time-like test particle under small perturbations (see Appendix for details).

Synchrotron flux density due to DM annihilations
The power of synchrotron emission per unit frequency ν emitted by an electron of given energy E and pitch angle θ, in a magnetic field B, is written as [55] In above equation, F (ν/ν c ) is the standard function depicting the spectral behavior of synchrotron radiation, K 5/3 is the modified Bessel function of order 5/3, and critical synchrotron frequency is defined as ν c = 3eBE 2 sin θ/(4πm 3 e c 5 ). Taking average over pitch angles, one arrived where dΩ obs is the element of the solid angle subtended by the image plane, dI syn (ν) is the differential specific intensity of synchrotron emission, and the integration is performed along the line of sight. The specific intensity I syn and emissivity j syn are related by standard radiative transfer equation [56]: where the increment ds along a line of sight, and α(ν, s) is the absorption coefficient. The absorption process is dominated by the synchrotron self-absorption, which is expected only to be important at low frequencies. As pointed out in Ref. [51], the synchrotron selfabsorption is negligible if only the synchrotron losses are considered which is the same for our case. Then Eq. (3.12) reduces to I syn (ν) = dsj syn (ν, s)/4π. We next briefly discuss the issue of gravitational redshift. The synchrotron emission produced at the innermost spike region will be influenced by the gravitational potential, the modification of photon frequency due to relativistic Doppler effect and gravitational redshift must be taken into account. For this purpose, notice that quantity I syn (ν)/ν 3 is frame-invariance due to the Liouville theorem, which builds connection between the observed and the emitted specific intensity as where g = ν obs /ν em is the redshift factor which encode both relativistic Doppler effects and gravitational redshift (see Appendix for detailed calculation). In principle, the realistic astrophysical black hole should be straightforward to extent type which is parameterized by its mass M and the dimensionaless spin parameter a = Jc 2 /M . Since current EHT resolution still can not fully determine the spin parameter of SMBH [40,42], we here simply consider Schwarzschild case to derive a conservative limit, and it is straightforward to extent our discussion to Kerr one follow from the expression of redshift factor in Appendix. The redshift factor for Schwarzschild black hole takes a simple form as g = 1 − R Sch /r. With this expression, substituting Eq. (3.14) into Eq. (3.12), one obtains observed total flux density. Here the integration of dΩ obs is performed over the EHT image field.

Constraints on DM parameters from the EHT observations
In this work we use the observational data given in Refs. [57,58], which adopted array of the EHT to map out M87 at a resolutions of ∼ 120 µas, corresponding to a spatial resolution of ∼ 8R Sch . And the Schwarzschild radius of SMBH is estimated as R Sch = 5.9 × 10 −4 pc, corresponding to angular size 7.3 µas. And the total flux density in the image field of 120 µas is constrained to be 0.6 Jy and distance to us is about 16.7 Mpc. We follow the EHT collaboration to set the magnetic field strength surrounding M87 to be B ∼ 1 − 30 G [58], which is estimated using the reconstruction of Faraday rotation as well as the image's brightness and flux density based on a one-zone emission model. Higher values of the magnetic field predict too small electron temperature to explain the brightness temperature in the EHT image and are thus disfavored by observations. The estimate of the magnetic field strength is expected to hold at least up to 6R Sch , which fully covers the integrated range in our calculation. The spectral energy distributions (SEDs) of the synchrotron emission from the DM annihilation around the SMBH can be derived, as shown in Fig. 2 for some examples. As a comparison, we also plot the SEDs from the NFW halo profile, which are much lower than those from the spike profile. This can also be expected according to the density profile shown in Fig. 1. Two processes shape the number density of electrons: the energy gain through the adiabatic compression and the energy loss through the synchrotron emission. At high energy end, the energy loss plays an important role and therefore the radiation spectrum is contributed dominantly by the electrons from annihilation directly. However, the low-energy electrons can acquire some energy when they fall into the black hole before the gain and loss are balanced, so the SED, in the low energy range, would be harder than that without convection. Since there are more low-energy electrons in the bb annihilation channel than those in the leptonic channels such as µ + µ − , the synchrotron flux at 230 GHz, at which the energy gain through the adiabatic compression would be more important, in the former case would be larger, as illustrated in Fig. 2. And when the magnetic field is weaker, the density can be contributed by the electrons generated in a larger space as indicated by the upper bound in Eq. (3.8) and more low-energy electrons could gain the energy during the infall as well. This larger number density of electrons would compete with the low emissivity per electron in small magnetic field and a larger synchrotron flux would be produced at 230 GHz. Therefore, the predicted flux become stronger for weaker magnetic fields, as shown in Fig. 2. We can then set excluded regions on the σv − m χ plane by comparing the expected flux with the observed flux around M87 , which by presented as the black dot-dashed line in Fig. 3. Since the backgrounds around the SMBH are extremely complicated, we only set conservative limits by requiring the expected fluxes not to exceed the observed one. We performed a grid scan in the σv − m χ plane and calculate the exclusion regions for DM masses from sub-GeV to TeV. The results are presented in Fig. 4, where the light and dark blue regions respectively correspond to limits based on the magnetic fields of 1 G and 30 G. It is interesting to note that, the constraints for 1 G magnetic field are stronger than those for 30 G magnetic field. This is mainly because the DM-induced electron density is effectively higher due to the infall of electrons from a larger space volume for a smaller cooling in a weaker magnetic field. For comparison, we also plot the 95% upper limits from the AMS (blue line; [27]), the Fermi-LAT (green line; [28]), H.E.S.S. (purple line; [31]), and the energy injection into CMB plasma during the dark ages (red line; [29]).

Discussion
In this paper, we investigated the DM-induced synchrotron emission in M87 . Assuming a spike profile of DM in the close vicinity of the SMBH, we derive new limits on the DM annihilation cross section using the EHT observations at 230 GHz, for µ + µ − , τ + τ − , bb and e + e − channels for typical mass regions of WIMP DM. Since the spike density in the inner regions would be further enhanced for relatively small annihilation cross section, the expected flux would also be amplified. As a consequence, unlike other existing limits which are usually more sensitive to the large annihilation cross sections, in our case the small annihilation cross section within σv ∼ 10 −34 −10 −27 cm 3 s −1 is more severely constrained. On the other hand, the spike density would be depleted for sufficiently large σv and reaches saturation, resulting in no constraints when σv is larger than ∼ 10 −27 cm 3 s −1 .
Moreover, due to a non-monotonic dependence of the injection spectrum on the magnetic field, more stringent constraints are given for a lower magnetic field. If the DM spike does exist, our results are more sensitive in the small cross section region, which are complementary to other experiments. The DM distribution without enhancement from the spiky structure predicts much lower synchrotron emission flux, and  Constraints on the WIMP pair annihilation cross section as a function of the WIMP mass at µ + µ − , τ + τ − , bb, e + e − channels by EHT M87 result, the dark and light blue contour limits assume magnetic field are 30 G and 1 G, respectively. Other constraints, as the legend applies to all panels, from the AMS (blue line; [27]), the Fermi-LAT (green line; [28]), H.E.S.S. (purple line; [31]), CMB (red line; [29]) and thermal relic abundance (black dot-dashed).
thus no effective constraint can be obtained. The forthcoming EHT observations of the SMBH in our Milky Way is expected to improve significantly the results of the current work due to its proximity. In addition, observations at multiple frequencies should also be very helpful in better identifying the DM signal or constraining the model parameters.
frequency of a test-particle at the emission radius r e , Ω = u φ /u t , in the free-fall model, can be written as [60], Ω(r, θ) = 2a 2 M 3 r (r 2 + a 2 ) 2 − a 2 (r 2 − 2M r + a 2 ) sin 2 θ . (4. 3) The source emission is assumed to be monochromatic and isotropic, with a Gaussian intensity. Using the normalization condition g µν u µ em u ν em = −1, we have and therefore, g = −g tt − 2g tφ Ω − g φφ Ω 2 1 + λΩ . (4.5) Here, λ = k φ /k t is a constant of the motion along the photon path. For Schwarzschild black hole with a = 0, Eq. (4.5) simplified to g = 1 − R Sch /r. For Kerr black hole, the radius of ISCO has expression [61]: where the ∓ sign is for test particles in prograde and retrograde orbits, respectively, relative to the spin of the BH and we the coefficients C 1,2 are given as C 2 = 3a 2 + C 2 1 . (4.8) Again, taking a = 0, one obtains r ISCO = 3R Sch for Schwarzschild black hole.