\section{Introduction}

The nature of Dark Matter (DM) is one of the biggest problems in modern physics and cosmology \cite{Rubin:1980zd, Clowe:2006eq, Larson:2010gs}. Quite a lot of particle models beyond the standard model were proposed in literature \cite{Jungman:1995df, Bertone:2004pz, Feng:2010gw, Roszkowski:2017nbc}. Some classical candidates are well-motivated, such as weakly interacting massive particles (WIMPs) \cite{Roszkowski:2017nbc, Liu:2017drf}, axions or axion-like particles \cite{Peccei:1977hh, Weinberg:1977ma, Davoudiasl:2019nlo, Yuan:2020xui}, and dark photons \cite{An:2020jmf, Caputo:2021eaa}. Meanwhile, many detecting methods and technologies have also been developed, such as detecting the feasible energy missing in particle collisions \cite{Slatyer:2017sev}, observing the scattering signals between DM and detector particles \cite{Liu:2017drf,Aprile:2020tmw}, and measuring the DM induced products in cosmic rays \cite{Porter:2011nv, Ge:2020yuf, Cao:2020bwd, Chen:2020gcl}. 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~\cite{Profumo:2010ya, McDaniel:2017ppt, McDaniel:2018vam}.
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.~\cite{Yuan:2009yy, Huang:2011xr, Bringmann:2012vr, Bergstrom:2013jra, Ackermann:2015zua, Slatyer:2015jla, Li:2015kag, Abdallah:2016ygi, Li:2018kgy,Cang:2020exa}. 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 \cite{Ferrarese:2004qr,Vasiliev:2008uz,Gondolo:1999ef, Gnedin:2003rj, Fields:2014pia}.
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) \cite{Middelberg:2008qc}, 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$^\star$ to the event horizon scale with extremely high angular resolution \cite{Akiyama:2019cqa, Akiyama:2019bqs, Akiyama:2019fyp, Akiyama:2019eap}. 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 \cite{Davoudiasl:2019nlo, Chen:2019fsq, Chen:2021lvo}. 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$^\star$ results~\cite{Akiyama:2019cqa, Akiyama:2019bqs}. For the DM density distribution in the center of the galaxy, a spike profile is expected due to the adiabatic accretion of the SMBH \cite{Gondolo:1999ef, Ullio:2001fb, Gnedin:2003rj}, 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$^\star$. 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.~\ref{spike}, we present the spatial distribution of DM around the SMBH, especially the DM spike profile. In Sec.~\ref{electron propagation around SMBH}, 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.~\ref{summary}.

\section{Density profile of DM halo with SMBH}\label{spike}

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~\cite{Navarro:1995iw, Navarro:1996gj}. It is parameterized as
\begin{align}\label{NFW_profile}
\rho_{\mathrm{NFW}}(r)=\frac{\rho_0}{\left(r / r_{0}\right)\left(1+r / r_{0}\right)^{2}},
\end{align}
where $r_0$ is the scale radius and $\rho_0$ is the normalization constant. Eq.~(\ref{NFW_profile}) shows a $r^{-1}$ density cusp at the center.

However, as proposed by Gondolo and Silk~\cite{Gondolo:1999ef}, 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 \cite{Merritt:2002vj,Gnedin:2003rj,Bertone:2005hw}. Nevertheless, we assume that the DM spike had formed around M87$^\star$ and was survived from such dynamical effects. In this scenario, the spike density profile reads~\cite{Merritt:2002vj}:
\begin{equation}\label{spike_profile}
\rho_{\rm sp}(r)=\rho_0\left(\frac{R_{\rm sp}}{r_0}\right)^{-\delta}\left(\frac{R_{\rm sp}}{r}\right)^{\delta_{\rm sp}},
\end{equation}
where $\delta$ is the slope of the original density profile, and $\delta_{\rm sp}=(9-2\delta)/(4-\delta)$. We take $\delta = 1$ to match the NFW profile. $R_{\rm sp}$ is the spike radius, which is given as
\begin{equation}\label{spike_radius}
R_{\rm sp}=\alpha_\delta r_0 \left(\frac{M_{\rm BH}}{\rho_0 r_0^3}\right)^{\frac{1}{3-\delta}}.
\end{equation}
In the above equation, $M_{\rm BH}$ is the mass of black hole, $\alpha_{\delta}$ 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 $\rho_{\rm sat}= m_\chi/\langle\sigma v\rangle t_{\rm BH}$. where $t_{\rm BH}$ is the age of the black hole, and $m_{\chi}$ and $\langle\sigma v\rangle$ are repectively the mass and thermally averaged annihilation cross section of DM particle. Using $\partial \rho_{\chi}/\partial t= - \langle\sigma v\rangle \rho_{\chi}^2/m_{\chi}$, one obtains spike density within the
region $R_{\rm Sch} \leq r<R_{\rm sp}$~\cite{Gondolo:1999ef}:
\begin{align}\label{spike_density}
\rho_\chi(r)=\frac{\rho_{\rm s p}(r) \rho_{\rm s a t}}{\rho_{\rm s p}(r)+\rho_{\rm sat}},
~\end{align}
where $R_{\rm 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 \cite{Aloisio:2004hy, Lacroix:2016qpq}:
\begin{align}
\label{profile}
\rho_\chi(r)=\left\{\begin{array}{ll}
0 & \qquad r<R_{\rm Sch},\\
\frac{\rho_{\rm s p}(r) \rho_{\rm s a t}}{\rho_{\rm s p}(r)+\rho_{\rm sat}} & \qquad R_{\rm Sch} \leq r<R_{\rm sp}, \\
\rho_{\mathrm{NFW}}(r)
& \qquad r \geq R_{\rm sp}.
\end{array}\right.
\end{align}

\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.8\linewidth]{M87_spike.pdf}
\end{center}
\caption{Dark Matter Density Profile of M87$^\star$. Density distribution for an NFW halo profile (black dashed) with $\delta=1$ and the same halo with a spike in the central region for a non-annihilating DM (red dashed). Also shown are DM with $m_\chi=100$ GeV and thermally average annihilation cross section $\langle\sigma v\rangle= 10^{-31}, 10^{-29}, 10^{-27}$, and $10^{-25}~\mathrm{cm}^{3}~{\rm s}^{-1}$ (blue, green, orange and purple solid line), respectively. The dot-dashed black vertical line shows corresponding Schwarzschild radius of M87$^\star$, and the gray shadowed region corresponding to the areas covered by EHT observation. In this plot, the lifetime of SMBH is chosen with $10^{9}$ years and $R_{\rm sp}\simeq 220$ pc.}
\label{Fig.spike}
\end{figure}

According to the EHT results, the mass of M87$^\star$ is $6.4 \times 10^{9} ~M_{\odot}$, and the corresponding Schwarzschild radius is $R_{\rm Sch}=5.9\times 10^{-4}$~pc. Taking $\alpha_{\delta}=0.1$, we fix $r_0=20$~kpc for the halo (as for the Milky Way), and normalize $\rho_{0} \approx 2.5~\mathrm{GeV}~\mathrm{cm}^{-3}$ as Ref.\cite{Lacroix:2016qpq} done. Then, we could calculate $R_{\rm sp}\simeq 220$~pc throught Eq.(\ref{spike_radius}). Assuming the age of M87$^\star$ is $t_{\rm BH}=10^9$ yr and using Eq.~(\ref{profile}), we displayed the DM density profile of M87$^\star$ in Fig.~\ref{Fig.spike}. In this plot, the corresponding Schwarzschild radius of M87$^\star$ and the areas covered by EHT observations are also shown.

\section{Synchrotron emission induced by DM annihilations at the center of M87}\label{electron propagation around SMBH}

\subsection{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  \cite{Aloisio:2004hy, Regis:2008ij, Lacroix:2016qpq}
\begin{equation}
\label{propagation}
    -\frac{1}{r^{2}} \frac{\partial}{\partial r}\left[r^{2} D \frac{\partial f_e}{\partial r}\right]+v \frac{\partial f_e}{\partial r}-\frac{1}{3 r^{2}} \frac{\partial}{\partial r}\left(r^{2} v\right) p \frac{\partial f_e}{\partial p}+\frac{1}{p^{2}} \frac{\partial}{\partial p}\left(\dot{p} p^{2} f_e\right)=q(r, p)
\end{equation}
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\sqrt{R_{\rm 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~(\ref{propagation}) 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:
\begin{equation}
q(r,p) = \frac{c}{4\pi p^2}Q(r,E)=\frac{c}{4\pi p^2} \frac {\langle\sigma v\rangle\rho_\chi^2(r)}{2 m^2_{\chi}} \sum_i {\rm BR}_i\frac{d N^{\rm inj}_{e^\pm, i}}{d E}(E),
\end{equation}
where $\langle\sigma v\rangle$ is the thermally averaged annihilation cross-section, and $d N^{\rm inj}_{e^\pm,i}/dE$ is the $e^\pm$ injection spectrum through annihilation channels $i$ with branching ratios ${\rm BR}_i$. In this paper, we respectively consider annihilation channels $i=e^{+}e^{-},~\mu^{+}\mu^{-}, ~\tau^{+}\tau^{-},~b\bar{b}$ with $100\%$ branching ratios, and extract corresponding $e^\pm$ injection spectrum from the {\tt DarkSUSY} or {\tt PPPC} packages~\cite{Cirelli:2010xx}. The electron and positron energy spectrum is then written as
\begin{equation}\label{number density}
n_{e}(r, E)=\frac{4 \pi p^{2}}{c} f_{e}(r, p).
\end{equation}

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_{\rm acc}$, the only energy loss mechanism for electron/positrons are radiative processes discussed above. While for $r\leq r_{\rm 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_{\rm 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~(\ref{propagation}) yields an integral analytic solution~\cite{Aloisio:2004hy}:
\begin{equation}\label{propagation_solution}
f_{e}(r, p)= \int_{r}^{r_{\mathrm{acc}}}\frac {Q_{i}\left(R_{\mathrm{inj}},p_{\mathrm{inj}}\right)}{v(R_{\rm inj})}\left(\frac{R_{\mathrm{inj}}}{R_{\mathrm{Sch}}}\right)^{\frac{5}{2}}\left(\frac{p_{\mathrm{inj}}}{p}\right)^{4} \mathrm{~d} R_{\mathrm{inj}},
\end{equation}
where we choose accretion radius $r_{\rm acc}=2GM/v^2_f\sim 0.04$ pc for the typical galactic wind velocity $v_f \simeq 500-700~{\rm km}~{\rm s}^{-1}$. $p_{\rm inj}\equiv p_{\rm inj}\left(R_{\rm inj}; r, p\right)$ is the injection momentum of an electron/positron is injected at $R_{\rm inj}~(\geq 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 $\dot{p}_{\rm{ad}}$. The associated characteristic curves of $p_{\rm inj}\left(R_{\rm inj}; r, p\right)$ are obtained by solving following differential equation~\cite{Aloisio:2004hy,Lacroix:2016qpq}:
\begin{equation}\label{momentum gain}
\frac{d p}{d r}=\frac{\dot{p}_{\rm ad} (r, p)+\dot{p} (r, p)}{v(r)}\simeq \frac{\dot{p}_{\rm ad} (r, p)+\dot{p}_{\rm syn} (r, p)}{v(r)}.
\end{equation}
In above equation, approximation is due to the fact that the total radiative momentum loss $\dot{p} (r, p)$, is dominated by synchrotron radiation, and
\begin{align}
\dot{p}_{\rm ad} (r, p) &=-\frac{1}{3r^2}p\frac{\partial}{\partial r}\left[r^2 v(r)\right], \nonumber\\
\dot{p}_{\rm s y n} (r, p)&=\frac{4}{3} \sigma_{T} \frac{B^2(r)}{8 \pi} \frac{E p}{m_e^2 c^3}.
\end{align}
Following method introduced in Refs.~\cite{Aloisio:2004hy,Lacroix:2016qpq}, for a homogeneous magnetic field in the galactic halo, the analytical solution of $p_{\mathrm{inj}}\left(R_{\mathrm{inj}} ; r, p\right)$ for Eq.~(\ref{momentum gain}) is given by
\begin{equation}\label{injection_momentum}
    p_{\mathrm{inj}}\left(R_{\mathrm{inj}} ; r, p\right)=p\left[\frac{k_{0} R_{\mathrm{Sch}}^{-\frac{1}{2}}}{c} R_{\mathrm{inj}}^{\frac{3}{2}} p\left(\frac{r}{R_{\mathrm{inj}}}-1\right)+\left(\frac{R_{\mathrm{inj}}}{r}\right)^{\frac{1}{2}}\right]^{-1}.
\end{equation}
It should be emphasize that the denominator in Eq.~(\ref{injection_momentum}) 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_{\rm inj}, p_{\rm inj})$ for a given arrival point $(r, p)$. In practice, $p_{\rm inj}$ holds positive for $R_{\rm inj} < R^0_{\rm inj}$ with
\begin{equation}\label{Rinj0}
R_{\mathrm{inj}}^{0}=r+\frac{c}{k_{0} p}\left(\frac{r}{R_{\mathrm{Sch}}}\right)^{-\frac{1}{2}}.
\end{equation}
We then use the value of ${\rm min}[r_{\rm acc}, R^0_{\rm inj}]$ as an effective upper bound for the integration in Eq.(\ref{propagation_solution}). While for lower bound for the integration, we take $r=r_{\rm ISCO}$ with $r_{\rm 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).

\subsection{Synchrotron flux density due to DM annihilations}

The power of synchrotron emission per unit frequency $\nu$ emitted by an electron of given energy $E$ and pitch angle $\theta$, in a magnetic field $B$, is written as~\cite{Ghisellini_2013}
\begin{align}
P_{\mathrm{syn}}(\nu, E_e, B, \theta) &=\frac{\sqrt{3} e^{3} B \sin \theta}{m_{\mathrm{e}} c^{2}} F\left(\nu/\nu_{c}\right), \nonumber\\
F(\nu/\nu_c) &\equiv \frac{\nu}{\nu_c}\int_{\nu/\nu_c}^{\infty} d y K_{5 / 3}(y).
\end{align}
In above equation, $F(\nu/\nu_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 $\nu_c=3e B E^2 \sin \theta/(4\pi m^3_e c^5)$. Taking average over pitch angles, one arrived
\begin{align}
\langle P_{\rm syn}\rangle(\nu, E, B) =\frac{1}{2}\int_{0}^{\pi} d \theta \sin \theta P_{\rm syn}(\nu, E_e, B, \theta),
\end{align}

Equipped with the $e^\pm$ energy spectrum $n_e(r,E)$ in Eq.~(\ref{number density}), the emissivity of synchrotron radiation is given by convoluting $n_e$ with the averaged power $\langle P_{\rm syn}\rangle$~\cite{Ghisellini_2013}:
\begin{align}
j_{\rm s y n}(\nu, r)=2 \int_{m_{e}}^{M_{\chi}} dE \langle {P}_{\rm s y n} \rangle(\nu, E, B) n_e(r,E),
\end{align}
where the factor 2 takes into account electrons and positrons. The integrated flux density measured by a detector can be estimated as
\begin{equation}\label{flux density}
S_{\rm syn}(\nu)= \int d\Omega_{\rm obs}\int_{l.o.s}dI_{\rm syn}(\nu),
\end{equation}
where $d\Omega_{\rm obs}$ is the element of the solid angle subtended by the image plane, $dI_{\rm syn}(\nu)$ is the differential specific intensity of synchrotron emission, and the integration is performed along the line of sight. The specific intensity $I_{\rm syn}$ and emissivity $j_{\rm s y n}$ are related by standard radiative transfer equation~\cite{Rybicki:2004hfl}:
\begin{align}
\frac{d I_{\rm syn}(\nu, s)}{d s}=-\alpha(\nu, s) I_{\rm syn}(\nu, s)+\frac{j_{\rm syn}(\nu, s)}{4 \pi},
\end{align}
where the increment $ds$ along a line of sight, and $\alpha(\nu, 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.~\cite{Aloisio:2004hy}, the synchrotron self-absorption is negligible if only the synchrotron losses are considered which is the same for our case. Then Eq.~(\ref{flux density}) reduces to $I_{\rm syn}(\nu)=\int ds j_{\rm syn}(\nu, s)/4 \pi$.

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_{\rm syn}(\nu)/\nu^3$ is frame-invariance due to the Liouville theorem, which builds connection between the observed and the emitted specific intensity as
\begin{equation}\label{redsfhit}
I_{\rm obs}\left(\nu_{\rm obs}\right) = \left(\frac{\nu_{\rm obs}}{\nu_{\rm em}}\right)^3I_{\rm em}\left(\nu_{\rm em}\right)= g^3 I_{\rm em}\left(\nu_{\rm em}\right),
\end{equation}
where $g=\nu_{\rm obs}/\nu_{\rm 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~\cite{Akiyama:2019cqa, Akiyama:2019fyp}, 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=\sqrt{1-R_{\rm Sch}/r}$. With this expression, substituting Eq.~(\ref{redsfhit}) into Eq.~(\ref{flux density}), one obtains observed total flux density. Here the integration of $d\Omega_{\rm obs}$ is performed over the EHT image field.

\subsection{Constraints on DM parameters from the EHT observations}\label{results and discussion}

In this work we use the observational data given in Refs.~\cite{Akiyama:2021qum, Akiyama:2021tfw}, which adopted array of the EHT to map out M87$^\star$ at a resolutions of $\sim 120~\mu$as, corresponding to a spatial resolution of $\sim 8 R_{\rm Sch}$. And the Schwarzschild radius of SMBH is estimated as $R_{\rm{Sch}}=5.9\times 10^{-4}~{\rm pc}$, corresponding to angular size $7.3~\mu {\rm as}$. And the total flux density in the image field of 120 $\mu {\rm 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$^\star$ to be $B \sim 1-30~\mathrm{G}$~\cite{Akiyama:2021tfw}, 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_{\rm{Sch}}$, which fully covers the integrated range in our calculation.

\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{sed.pdf}
\caption{Spectra of synchrotron emission from DM annihilation, for $b\bar{b}$ and $\mu^{+}\mu^{-}$ channels and magnetic field of 1~G and 30~G (solid lines). The DM mass and cross section are chosen as $m_{\chi}=100$~GeV and $\langle\sigma v\rangle = 10^{-30}~ \mathrm{cm}^{3}~\mathrm{~s}^{-1}$. For comparison, the contribution of regular NFW profile for $\mu^{+}\mu^{-}$ channel are also shown (dot-dashed lines).
The vertical dot-dashed  line corresponds to the observational frequency, i.e., 230 GHz.}
\label{sed}
\end{figure}

The spectral energy distributions (SEDs) of the synchrotron emission from the DM annihilation around the SMBH can be derived, as shown in Fig.~\ref{sed} 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.~\ref{Fig.spike}. 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 $b\bar{b}$ annihilation channel than those in the leptonic channels such as $\mu^+\mu^-$, 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.~\ref{sed}. 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.~(\ref{Rinj0}) 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.~\ref{sed}.

\begin{figure*}[htbp]
\centering
\includegraphics[width=0.49\textwidth]{flux1G.pdf}
\includegraphics[width=0.49\textwidth]{flux30G.pdf}
\caption{The predicted flux of different cross section for four channels. The left/right panels are fluxes predicted by magnetic field at 1~G/30~G, respectively. And the cross-section used in left panel and right panel are $\langle\sigma v\rangle = 10^{-32}~ \mathrm{cm}^{3}~\mathrm{~s}^{-1}$ and $\langle\sigma v\rangle = 10^{-30}~\mathrm{cm}^{3} \mathrm{~s}^{-1}$. The red line is $b\bar{b}$, red line is $e^{+}e^{-}$, green line is $\mu^{+}\mu^{-}$, gray line is $\tau^{+}\tau^{-}$, and black dot-dashed line is the observed flux.}
\label{flux}
\end{figure*}

We can then set excluded regions on the $\langle\sigma v\rangle-m_\chi$ plane by comparing the expected flux with the observed flux around M87$^\star$, which by presented as the black dot-dashed line in Fig.~\ref{flux}. 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 $\langle\sigma v\rangle-m_\chi$ plane and calculate the exclusion regions for DM masses from sub-GeV to TeV. The results are presented in Fig.~\ref{bounds}, 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; \cite{Bergstrom:2013jra}), the Fermi-LAT (green line; \cite{Ackermann:2015zua}), H.E.S.S. (purple line; \cite{Abdallah:2016ygi}), and the energy injection into CMB plasma during the dark ages (red line; \cite{Slatyer:2015jla}).

\begin{figure*}[htbp]
\centering
\includegraphics[width=0.49\textwidth]{mumu.pdf}
\includegraphics[width=0.49\textwidth]{tautau.pdf}
\includegraphics[width=0.49\textwidth]{bb.pdf}
\includegraphics[width=0.49\textwidth]{ee.pdf}
\caption{Limits on Dark Matter Annihilation. Constraints on the WIMP pair annihilation cross section as a function of the WIMP mass at $\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$, $b\bar{b}$, $e^{+}e^{-}$ channels by EHT M87$^\star$ 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; \cite{Bergstrom:2013jra}), the Fermi-LAT (green line; \cite{Ackermann:2015zua}), H.E.S.S. (purple line; \cite{Abdallah:2016ygi}), CMB (red line; \cite{Slatyer:2015jla}) and thermal relic abundance (black dot-dashed).}
\label{bounds}
\end{figure*}

\section{Discussion}\label{summary}

In this paper, we investigated the DM-induced synchrotron emission in M87$^\star$. 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 $\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$, $b\bar{b}$ 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 $\langle\sigma v\rangle \sim 10^{-34}-10^{-27}{\rm cm}^3~{\rm s}^{-1}$ is more severely constrained. On the other hand, the spike density would be depleted for sufficiently large $\langle\sigma v\rangle$ and reaches saturation, resulting in no constraints when $\langle\sigma v\rangle$ is larger than $\sim 10^{-27}{\rm cm}^3~{\rm 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 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.

\section*{Acknowledgements}

We thank Yifan Chen, Yi-Zhong Fan, Lei Feng, Lingyao Kong, Xu Pan, Chi Tian and Yue-Lin Sming Tsai for helpful discussion and valuable comments in various aspects. R.D. is supported in part by the National Key R\&D Program of China (2021YFC2203100). Z.Q.S is supported by the National Natural Science Foundation of China (NSFC) under Grants No. U1738210 and No. 12003074; X.H. and Q.Y. are supported by Chinese Academy of Sciences and the Program for Innovative Talents and Entrepreneur in Jiangsu.


\section{Introduction}
    Since the 1960s, superconducting quantum interference devices (SQUIDs) have been the benchmark for sensitive magnetic field measurement at low frequencies~\cite{Ryhanen1989-ci,Pizzella2001-pl}.
    However, developments in the spin-exchange relaxation-free (SERF) regime over the past two decades have allowed optically-pumped magnetometers (OPMs) to achieve a similar, if not better sensitivity compared to SQUIDs \citep{kominis_subfemtotesla_2003}, with a number of distinct advantages.
    SQUIDs require cryogenic cooling which increases costs, limits the device portability, and increases the separation between the device and object of interest.
    OPMs offer reduced separation distances over SQUIDs, which is of particular interest for magnetoencephalography (MEG)~\cite{Sander2012-kv,Colombo2016-hg,Sheng2017-sp,Limes2020-eg}, where OPM-MEG is able to offer an improved spatial~\cite{Tan1990-yx,Boto2016-ya}, and  spatio-temporal~\cite{gialopsou_improved_2021} resolution.

    These advantages can be similarly applied to the diagnostic measurement of batteries through non-invasive current mapping~\cite{Hu2020-kn,bason_non-invasive_2021}, and for fundamental physics research~\cite{Pustelny2013-dq,Abel2019-tm}.
    While OPMs offer many advantages over SQUIDS, a number of challenges remain, especially around the construction of dense arrays.
    Commercially available OPMs currently focus on single-sensor designs, and when multiple units are operated in close proximity issues arise from radio frequency (RF) cross-talk, as well as high sensor surface temperatures.
    To achieve dense OPM array systems, these problems must be addressed.

    Here we present a magnetometer system designed with an emphasis on scalability and modularity, both for easy adaptability of sensor locations and sensing regimes.
    The magnetometer is constructed from individual modules, used to configure the layout of the OPM system.

    Additionally, optical elements inside modules can be individually replaced or reconfigured to allow implementation of different sensing schemes.

    To illustrate the adaptability of the magnetometer system, two configurations of the instrument operating in the ultra-sensitive SERF regime are presented. A single sensor configuration was used to capture alpha-rhythms from the brain of a human participant, demonstrating suitability for high sensitivity measurements. Additionally, a two-sensor gradiometer with superior sensitivity performance through common noise rejection was characterised.

    The scalability of the modular OPM system, along with the flexibility to incorporate a variety of magnetometry schemes, make this system an ideal development platform for the realisation of dense sensor arrays.

\section{Method} \label{sec.Method}
    \begin{figure}
        \centering
        \includegraphics[width=\columnwidth]{images/Modular_OPM_2x2.png}
        \caption{An example of multi-axis magnetic gradiometer connecting 4 sensor modules in $2\times 2$ array configuration with light source and beam distribution modules. The laser beam is conditioned with lens, linear polariser (LP) and quarter-waveplate ($\lambda/4$). After passing through vapour cells, the beam power is monitored with photodiodes (PD). A non-polarising $50:50$ beam splitter ($50:50$ BS) and a mirror ensure that the probe light from a single beam is delivered to all sensors.}
        \label{fig:2x2}
   \end{figure}

    Our system utilises \SI[product-units = repeat]{45 x 45 x 40}{\milli\metre} cubic modules which inter-connect to produce a desired OPM system layout. An example of $2\times 2$ sensor array is shown in Fig.~\ref{fig:2x2}.
    The enclosure of each module is 3D-printed from polycarbonate, allowing for design flexibility while also providing excellent mechanical rigidity with a \SI{139}{\degree C} glass transition temperature~\cite{Shamim2014-zu}.
    The system is assembled from three types of modules: light source, sensor, and beam distribution modules. Light source modules introduce and condition laser light for use in the magnetometer system. Light can be internally generated by a vertical cavity surface-emitting laser (VCSEL), or brought from external sources via optical fibres.

            \begin{figure}
            \centering
            \includegraphics[width=1\columnwidth]{images/Modular_OPM_grad.png}
            \caption{

            A schematic (top) and a photograph (bottom) of the single-axis gradiometer modular OPM consisting of a light source module and two sensor modules. The dimensions of each module are $W\times D\times H=\SI[product-units = repeat]{45 x 45 x 40}{\milli\metre}$. Modules can be snapped together to produce different array configurations.

            }
            \label{fig:optics}
        \end{figure}

    Sensor modules contain optical elements, vapour cells, magnetic field coils and photodiodes required for magnetic field measurements.
    Finally, beam distribution modules are used for light delivery, allowing for greater flexibility in positioning of the sensors.
    Due to the adaptability of its optics packages, measurement schemes and array configurations, the magnetometer can be tailored to a large variety of applications.

    The sensor module design presented in this paper operates in the high-sensitivity SERF regime, requiring large atomic densities in a low magnetic field environment~\citep{allred_high-sensitivity_2002}.

    A gradiometer (Fig.~\ref{fig:optics}) is constructed using one light source module and two sensor modules in a similar configuration to~\cite{nardelli_conformal_2020}.
    Our modular system can further be extended to a multi-axis gradiometer instrument as shown in Fig.~\ref{fig:2x2}.

    The SERF sensor modules utilise the Hanle resonance (zero-field resonance) as a measurement scheme~\cite{Dupont-Roc1969-yj}, relying on the monitoring of the transmitted power of a resonant circularly polarised laser beam probing a vapour cell.

    When a magnetic field component that is perpendicular to the laser beam direction, $B_{x}$ (respectively $B_{z}$), is scanned through zero field, atoms are pumped into a dark state and the magnitude of the transmitted light follows a Lorentzian distribution with peak transmission at $B_{x}=0$ (respectively $B_{z}=0$). By applying a modulating field $B_{x,RF}\cos{(\omega t)}$ across the zero-field resonance peak, and passing the photodiode response through a lock-in amplifier (LIA), a dispersion signal is produced of the form~\cite{cohen-tannoudji_diverses_1970,tierney_optically_2019}
    \begin{equation}
        V(B_x)= A\frac{\gamma B_x \tau}{1 +(\gamma B_x \tau)^2},
    \end{equation}
    where $\gamma$ is the gyromagnetic ratio for the active atomic species,
    $\tau$ is the relaxation time of the spin-polarised atomic population, and $A$ is a proportionality constant (with units \SI{}{\volt/\pico\tesla}) determined by calibrating the sensor output.
    The LIA also reduces 1/f noise, improving the signal-to-noise ratio of the resonance signal~\cite{kastler_hanle_1973}. A similar argument can be made for the other orthogonal component $B_z$.

    A sensor module contains three orthogonal coil pairs positioned around a vapour cell, used to apply offset and modulation fields. The field $B_{x}$ is modulated at frequency $\omega=\SI{926}{\hertz}$ with a peak-to-peak of ~ \SI{60}{\nano\tesla} at the centre of the vapour cell.
    To achieve sufficient atomic density for the SERF regime, the vapour cells are heated using alternating currents at \SI{81}{\kilo\hertz} and \SI{200}{\kilo\hertz}, driven through twisted manganin heating wires wrapped around the vapour cells. These AC frequencies are chosen to be orders of magnitude higher than the rubidium atomic response bandwidth to avoid broadening of the measured resonance. The vapour cell housing is thermally insulated to minimise heating power and surface temperature of sensor module. In doing so, a stand-off distance of \SI{8.5}{\milli\metre} from the end of the sensor module housing to the cell centre was achieved.

    Light from an external cavity diode laser, tuned to the \ensuremath{{}^{87}\mathrm{Rb}}{} D1 transition (\SI{795}{\nano\metre}), is coupled into a polarisation maintaining (PM) optical fibre and input to the light-source module of the sensor (Fig.~\ref{fig:optics}).
    The linearly polarised beam is then collimated to a $1/e^2$ diameter of \SI{3}{\milli\metre} before being circularly polarised and split with a 50:50 non-polarising beamsplitter cube.
    The reflected beam is directed through a \SI[product-units = repeat]{5 x 5 x 5}{\milli\metre} blown-glass \ensuremath{{}^{87}\mathrm{Rb}}{} cell in sensor module 1, containing \SI{50}{\torr} of N$_{2}$ and \SI{650}{\torr} of Ne as buffer gases \footnote{Blown-glass \ensuremath{{}^{87}\mathrm{Rb}}{} cells were procured from Twinleaf LLC.}. In sensor module 2, an identical vapour cell is probed with the beam coming from the transmitted arm of the beamsplitter. Photodiodes are placed behind the vapour cells to monitor the probing beam power. Their signals are passed through external trans-impedance amplifiers and fed to a pair of synchronised LIAs before being recorded with a digital acquisition card (DAQ) \footnote{LabJack T7 pro, LabJack Corporation, Co, USA}.

    Sensitivity measurements were taken inside a four-layer mumetal cylindrical shield \footnote{MS-2 mu-metal shield, Twinleaf LLC} with an internal magnetic field magnitude below \SI{10}{\nano\tesla} and magnetic noise of \SI{10}{\femto\tesla/\sqrt{\hertz}}. The offset fields were applied using the sensor module's internal coils, and modulation fields were generated using field coils built into the shield. The LIA output voltages were calibrated by applying well-defined $B_x$ fields.

    To demonstrate suitability of the modular OPM system for MEG applications, a single-sensor magnetometer configuration was operated to capture human brain alpha-rhythms~\cite{Cohen1968-nf}. The setup is similar to the three-module gradiometer configuration shown in Fig. \ref{fig:optics}, with the difference that only one sensor module is present. The MEG measurements were taken inside a three-layer mumetal cylinder, with an internal diameter of \SI{50}{\centi\metre} and axial length of \SI{100}{\centi\metre}. The single-sensor module OPM was placed at the centre of the cylinder, supported by a wooden arm which the participant was able to rest their head against. A prerecorded voice prompted the participant to open and close their eyes every \SI{12}{\second} for \SI{150}{\second} to stimulate an alpha-band response when their eyes were shut (photic blocking)~\cite{Cohen1972-jw}. The optimum placement of the sensor over the primary visual cortex (Oz position) was measured according to the standard $10-10$ system~\cite{Koessler2009-pc}. The analogue signal from the OPM-MEG system was sampled at \SI{200}{\hertz} using the same DAQ setup as previously described.

\section{Results and discussion}
    \begin{figure}
            \centering
            \includegraphics[width=\columnwidth]{images/Noise_Floor_Figure.png}
            \caption{The \SIrange{1}{200}{\hertz} frequency compensated linear spectral density of the demodulated signals from the sensor module 1 (blue), the sensor module 2 (red) and from the balanced photodiodes gradiometer configuration (yellow). The gradiometer configuration removes common-mode noise between the two sensors, which can be seen by the noise reduction at \SI{50}{\hertz}, and the broader peak in the \SIrange{55}{75}{\hertz} band. Sensor module 1 has a noise-floor of \SI{65}{\femto\tesla/\sqrt{\hertz}}, sensor module 2, \SI{83}{\femto\tesla/\sqrt{\hertz}}, and the gradiometer configuration \SI{47}{\femto\tesla/\sqrt{\hertz}} in the \SIrange{5}{45}{\hertz} band.
            }
            \label{fig:ffts}
        \end{figure}

  The sensitivity measurements were taken over a period of \SI{60}{\second} during which the magnetometer signals from sensor modules 1 and 2 were recorded.
  Data were analysed for each sensor separately, and as a software gradiometer (sensor module 1 -- sensor module 2).
   The frequency spectra in Fig.\,\ref{fig:ffts} show noise floors of sensor modules 1 and 2 to be \SI{65}{\femto\tesla/\sqrt{\hertz}} and \SI{83}{\femto\tesla/\sqrt{\hertz}}, respectively, in the \SIrange{5}{45}{\hertz} band.
   As both the \ensuremath{{}^{87}\mathrm{Rb}}{} atoms and the LIA act as a low-pass filter, the linear density spectrum has been divided by the transfer function of each sensor module~\cite{colombo_four-channel_2016}. The transfer function follows a first-order low-pass filter from \SI{1}{\hertz} to \SI{250}{\hertz}, and the bandwidths of the sensor modules were found to be \SI{213}{\hertz} for module 1, and \SI{219}{\hertz} for module 2.
   Using the gradiometer configuration the sensitivity improves to \SI{47}{\femto\tesla/\sqrt{\hertz}}, exploiting the suppression of common light intensity noise. In addition, the gradiometer configuration offers a common-mode noise suppression ratio of 6.1 on the \SI{50}{\hertz} powerline peak compared to module 2. When considering the inter-cell distances, the gradiometer measurement translates to a magnetic field gradient sensitivity of \SI{10.4}{\femto\tesla/\centi\metre/\sqrt{\hertz}} in the frequency range between \SI{5}{\hertz} and \SI{45}{\hertz}.
   For this measurement the \ensuremath{{}^{87}\mathrm{Rb}}{} vapour cells were operated at temperatures of \SI{129}{\degree C} and \SI{124}{\degree C} in sensor modules 1 and 2, respectively.
   We note that the optimal sensitivity of \ensuremath{{}^{87}\mathrm{Rb}}{} vapour in the SERF regime occurs at \SI{150}{\degree C}~\cite{Shah2007-sa}, leading to a factor 4 increase in vapour density~\cite{Seltzer2008-kh}, which would require an alternative to the polycarbonate-printed material.
   In addition to the potential sensitivity increase from the \ensuremath{{}^{87}\mathrm{Rb}}{} density, further gains can be made in the low-frequency domain by implementing active feedback on vapour cell temperature and laser power.

   These adaptions, currently in development, coupled with the sensor design presented here, will result in a fully modular OPM system, re-configurable by the user.

    \begin{figure}
        \centering
        \includegraphics[width=\columnwidth]{images/Window_Hann_summer_reduced_range.png}
        \caption{Spectrogram of the \SIrange{8}{12}{\hertz} alpha-band response to a participant opening/closing their eyes. The overlaid white trace denotes the audio cue trigger signal from the stimulus PC. High values of the trigger signal indicate the time when the participant is instructed to close their eyes, the low value is when the participant has their eyes open. The spectrogram colour scheme details the high periods of activity in yellow, and the low periods in green.}
        \label{fig:spectrogram}
    \end{figure}

   Data from the alpha-rhythm measurement with a single sensor module were analysed in the frequency domain, with the spectrogram shown in Fig.~\ref{fig:spectrogram}. We see an increase in the activity in the \SIrange{8}{12}{\hertz} alpha-band when the participant's eyes are closed, with a strong correlation with the TTL cue. Although we already show that the current modular sensor is suitable for MEG, future improvements of the signal-to-noise ratio can be made by suppressing common-mode noise using a gradiometer.

  We have designed a modular OPM system and have demonstrated its operation in two configurations, one of which was used to record MEG signals from a human participant.
  The ability of this highly-customisable OPM system to implement new measurement schemes, array configurations, dual-scheme OPM arrays, or multi-modal quantum sensor systems, offers a unique combination of adaptability and the capability for in-field use.
  This modular system is an important development tool for the advancement of dense OPM arrays in medical and industrial fields, providing a rapid testing platform to benchmark OPM system architectures.

\section*{Additional material}
    The participant gave written informed consent to take part, after explanation of the procedure and purpose of the experiment.
    The study has been approved by the Brighton and Sussex Medical School Research Governance and Ethics Committee (ER/BSMS3100/1).

    The data and code that support the findings of this study are available from the authors upon reasonable request.

\bibliographystyle{ieeetr}