The First Robust Evidence Showing a Dark Matter Density Spike Around the Supermassive Black Hole in OJ 287

Black hole dynamics suggests that dark matter would redistribute near a supermassive black hole (SMBH) to form a density spike. However, no direct evidence of a dark matter density spike around an SMBH has been identified. In this Letter, we present the first robust evidence showing a dark matter density spike around an SMBH. We revisit the data of the well-known SMBH binary OJ 287 and show that the inclusion of the dynamical friction due to a dark matter density spike around the SMBH can satisfactorily account for the observed orbital decay rate. The derived spike index γsp=2.351−0.045+0.032 gives an excellent agreement with the value γ sp = 2.333 predicted by the benchmark model assuming an adiabatically growing SMBH. This provides a strong verification of the canonical theory suggested two decades ago modeling the gravitational interaction between collisionless dark matter and SMBHs.


Introduction
In the past few decades, various studies showed that a supermassive black hole (SMBH) can alter the nearby dark matter density distribution to form a density spike (Young 1980;Gondolo & Silk 1999;Merritt 2004a;Gnedin & Primack 2004;Merritt 2004b;Sadeghian et al. 2013;Nampalliwar et al. 2021).The dark matter density would be steepened within the sphere of influence of a SMBH due to the conservation of angular momentum and radial action.Consequently, we expect that such a high dark matter density near a SMBH would significantly trigger the rate of dark matter annihilation to give a strong emission of high-energy gamma-rays (Bertone et al. 2002;Gnedin & Primack 2004;Fields et al. 2014;Shapiro & Shelton 2016).However, no strong gamma-ray emission has been detected so far near any SMBH, including the SMBH in our Galaxy (Sgr A*) (Fields et al. 2014).
Recently, Chan & Lee (2023) claim that the data of two nearby black hole low-mass X-ray binaries, A0620-00 and XTE J1118+480, might reveal the existence of dark matter density spikes around their black holes.The dynamical friction exerted by the dark matter density spikes can satisfactorily explain the abnormal large orbital decay rates of the companion stars in these two binaries.However, the evidence of any dark matter density spike around a SMBH is definitely lacking, even though many recent studies are still modelling the black hole binary inspirals with the existence of dark matter density spikes (Yue & Cao 2019;Tang et al. 2021;Dai et al. 2022;Becker et al. 2022;Li et al. 2022;Qunbar & Stone 2023).In particular, we expect that the gravitational waves (GWs) emitted by supermassive black hole binaries (SMBHBs) can help reveal the properties of SMBHs and verify our theoretical understanding about these exotic systems, although the emission of these low-frequency GWs can only be detected by future space GW interferometers.Moreover, confirming the existence of a dark matter density spike can also help constrain the parameters of dark matter annihilation or the decay rate (Gondolo & Silk 1999;Kar et al. 2023).
In this letter, we present the first evidence of the existence of a dark matter density spike around a SMBH.We revisit the data of the well-known SMBHB OJ 287, which consists of a secondary SMBH with mass m BH ≈ 1.5 × 10 8 M ⊙ orbiting a primary SMBH with mass M BH ≈ 1.8 × 10 10 M ⊙ .We show that the energy loss rate due to GW emission is significantly smaller than the observed energy loss rate, with a discrepancy of more than 4.3σ.With the assumption of a dark matter density spike around the primary SMBH, the dynamical friction can provide the extra energy loss rate which can satisfactorily account for the observed energy loss rate.The derived spike index gives an excellent agreement with the predicted value based on the adiabatically growing SMBH model (Gondolo & Silk 1999;Fields et al. 2014).Therefore, this is the strongest evidence so far to reveal the existence of dark matter density spike around a SMBH.

The supermassive black hole binary OJ 287
OJ 287 is a well-known SMBHB as it has been studied for more than a century (Sillanpää et al. 1988;Valtonen & Lehto 1997;Valtonen et al. 2008Valtonen et al. , 2010;;Dey et al. 2018Dey et al. , 2019;;Laine et al. 2020;Komossa et al. 2023;Titarchuk, Seifina & Shrader 2023;Valtonen et al. 2023;Zwick & Mayer 2023;Martinez 2023).This SMBHB contains one secondary SMBH orbiting another very massive primary SMBH so that it is one of the most exotic systems observed in our universe.This binary emits a large amount of X-ray and radio radiation and we expect that it could also emit a huge amount of GW energy.Based on the analysis of the accurately extracted starting epochs of ten optical outbursts of OJ 287 between the years of 1912-2016, an accurate orbit of the secondary SMBH in OJ 287 can be determined (Dey et al. 2018).This study has provided the most robust results and obtained accurate orbital parameters with very small uncertainties for OJ 287 (see Table 1 for the essential orbital parameters).
Although there are some recent studies claiming that the mass of the primary SMBH is only M BH ∼ 10 8 M ⊙ (Komossa et al. 2023) and the mass of the secondary SMBH is smaller by 20% (Titarchuk, Seifina & Shrader 2023), these results are respectively based on the analysis on a particular outburst observation in 2022 (Komossa et al. 2023) and the comparative study focusing on X-ray data only (Titarchuk, Seifina & Shrader 2023).Overall speaking, the orbital parameters given in Dey et al. (2018), which have followed the outburst data in the past 104 years, are still the most comprehensive and robust results for OJ 287.In the followings, we will base on these robust orbital parameters shown in Table 1 to perform our analysis.
As orbital precession exists in OJ 287, the orbit of the secondary SMBH can be described by the relative distance r(φ) between the primary and secondary SMBHs: where a is the orbital semi-major axis, e is the eccentricity, and α = ∆Φ/2π is the precession phase angle.
One of the most intriguing properties of OJ 287 is that the orbit of the secondary SMBH is shrinking.The orbital period decay rate is Ṗ = −(0.00099± 0.00006) (Dey et al. 2018), which means that two SMBHs will merge together after about 12000 years.The expected reason for this orbital shrinking is that energy is given out continuously due to GW emission (Valtonen et al. 2008;Dey et al. 2018).Previous studies following pulsar binaries have shown that the orbital shrinking rate of pulsars agrees with the predicted rate based on GW emission (Weisberg & Taylor 2005).Based on General Relativity, the energy loss rate due to GW emission can be analytically given by (Peters & Mathews 1963;Maggiore 2007;Tang et al. 2021;Li et al. 2022 where µ = m BH M BH /(m BH + M BH ) and M = m BH + M BH are the reduced mass and total mass of SMBHB OJ 287 respectively.
On the other hand, one can convert the orbital period decay rate to the total energy loss rate theoretically.Although the secondary SMBH is orbiting with a large precession angle, we can still apply the Keplerian relation P ∝ a 3/2 , though the proportionality constant is larger due to the precession effect.Also, since the General Relativistic terms only contribute about 1% of the total mechanical energy E and the orbital period decay rate is smaller than 0.1%, the total mechanical energy of the orbital motion can be well-approximated by the conventional Newtonian expression E = −GMµ/2a.Therefore, we can write the total energy loss rate in terms of the orbital period decay rate: (3) By using Eqs.( 2) and ( 3) with M BH = (1.8348± 0.0008) × 10 10 M ⊙ , m BH = (1.5013± 0.0025) × 10 8 M ⊙ , P = 12.067 ± 0.007 years, Ṗ = 0.00099 ± 0.00006, and e = 0.657 ± 0.001 obtained in Dey et al. (2018), and taking a = 1.72 × 10 17 cm based on orbital analysis (Valtonen & Lehto 1997;Laine et al. 2020;Martinez 2023), we get Ė = −(3.66± 0.24) × 10 41 W and ĖGW = −(2.62± 0.02) × 10 41 W. This gives a 4.3σ discrepancy between Ė and ĖGW , which shows a large tension between observations and theoretical prediction assuming solely GW emission.Therefore, we expect that there must exist another important energy loss mechanism in OJ 287.

Dynamical friction of the dark matter density spike
Based on the results of numerical simulations, the density profile of a massive halo formed by collisionless dark matter would follow the Navarro-Frenk-White (NFW) density profile (Navarro et al. 1996): where ρ s and r s are scale density and scale radius respectively.This profile is commonly modeled as the galactic dark matter density profile, including our Galaxy and the M31 galaxy (Sofue 2015).Nevertheless, the galactic central SMBH would re-distribute dark matter to form a dark matter density spike around the SMBH due to conservation of angular momentum and radial action (Gondolo & Silk 1999).Outside the spike region r ≥ r sp , the dark matter density would follow back to the global NFW density profile.To summarize, the dark matter density around the primary SMBH can be described by the following spike model (with General Relativistic correction) (Sadeghian et al. 2013;Eda et al. 2015;Tang et al. 2021;Capozziello, Zare & Hassanabadi 2023;John, Leane & Linden 2023): where For the benchmark model suggesting an adiabatic growth of SMBH, one can relate the spike index γ sp with the power-law index of the dark matter density outside the spike γ: γ sp = (9 − 2γ)/(4 − γ) (Gondolo & Silk 1999;Fields et al. 2014;Eda et al. 2015).As the NFW profile suggests γ = 1 for r ≪ r s , the adiabatic growth model predicts γ sp = 7/3 ≈ 2.333, which gives a very high dark matter density near a SMBH.Note that we did not include the effect of dark matter annihilation in our analysis.A large rate of dark matter annihilation would reduce the dark matter density spike to the so-called annihilation plateau density (Fields et al. 2014).As the age of the SMBH, mass of dark matter particles, and the annihilation cross section are unknown, we neglect the annihilation effect in our analysis.
In the followings, we describe a theoretical framework to model the unknown parameters ρ s , r s , ρ sp and r sp by the known parameters M BH and m BH .First of all, by considering 43 galaxy-scale strong gravitational lenses, there is an empirical relation between the mass of galactic SMBHs and the total dynamical mass of galaxies (Bandara, Crampton & Simard 2009): Surprisingly, this empirical relation gives an excellent agreement with the later simulation result M BH ∝ M 1.55±0.05tot (Booth & Schaye 2010) and it is consistent with the results for elliptical galaxies M BH ∝ M 1.6 tot (Bogdán & Goulding 2015).Based on the empirical relation in Eq. ( 6) with the corresponding uncertainties, we get M tot = 2.21 +3.67  −1.06 × 10 14 M ⊙ for the galaxy hosting OJ 287.By taking this total dynamical mass as the virial mass, we can calculate the virial radius by r 3 200 = 3M tot /4πρ 200 .Here, ρ 200 is defined as By adopting the values of the redshift of OJ 287 z = 0.306 (Benitez & Dultzin-Hacyan 1996), the cosmological density parameters Ω m = 0.315±0.007and Ω Λ = 0.685±0.007,and the Hubble constant H 0 = 67.4± 0.5 km/s/Mpc from Planck's observation (Planck Collaboration 2020), we get r 200 = 3.54 +1.36 −0.69 × 10 24 cm for the galaxy hosting OJ 287.
Furthermore, using the mass-concentration relation of cosmological structures, we can get the concentration parameter c 200 from the total dynamical mass M tot .The empirical mass-concentration relation based on the lensing data of galaxies and galaxy clusters can be written as (Xu et al. 2021) where h = 0.674, C 0 = 5.119 +0.183 −0.185 , γ c = 0.205 +0.010 −0.010 and log 10 (M 0 ) = 14.083 +0.130 −0.133 .Again, by including the corresponding uncertainties, we get c 200 = 4.17 +0.99  −0.55 for OJ 287, which gives r s = r 200 /c 200 = 8.49 +1.01 −0.62 × 10 23 cm.Furthermore, since we know M tot and r s , we can get ρ s = 6.84 +4.23  −1.83 × 10 −26 g cm −3 .Based on the standard spike model, the spike radius r sp is empirically defined by r sp = 0.2r in , where r in is the radius of influence (Fields et al. 2014;Eda et al. 2015;Kavanagh et al. 2020).The radius of influence can be determined by (Merritt 2004a,b;Eda et al. 2015;Kavanagh et al. 2020): Therefore, we get the following analytic relation (Eda et al. 2015;Kavanagh et al. 2020;Mukherjee et al. 2023) Also, using Eq. ( 5) and considering at r = r sp , we have The above two relations Eq. ( 10) and Eq. ( 11) can connect ρ sp and r sp individually with the spike index γ sp .
Since the dark matter density is extremely high near the primary SMBH, the effect of dynamical friction would also be very large.In fact, the effect of dynamical friction of a dark matter density spike has been theorized for a long time.Most theoretical studies have included the effect of dynamical friction in modelling black hole mergers (Dai et al. 2022;Becker et al. 2022;Qunbar & Stone 2023;Mukherjee et al. 2023).The energy loss rate due to dynamical friction is given by (Chandrasekhar 1943;Yue & Cao 2019) where ln Λ ≈ ln M BH /m BH is the Coulomb Logarithm (Kavanagh et al. 2020), ξ(σ) ≈ 1 is a numerical factor depending on the dark matter velocity dispersion σ, and v = (GM BH /p) 1/2 (e 2 − 1) + 2[1 + e cos(1 − α)φ] is the orbital velocity, with p = a(1 − e 2 ) (Tang et al. 2021).If dark matter particles follow a Maxwellian distribution, the numerical factor can be described by ξ(σ) = erf(X) − 2Xe −X 2 / √ π, where X = v/ √ 2σ (Merritt 2013).Assuming the dark matter velocity dispersion is close to the velocity dispersion in the galactic bulge, we can get σ ∼ 440 km/s by using the SMBH-velocity dispersion relation M BH = 10 8.32±0.05M ⊙ (σ/200 km/s) 5. 64±0.32 (McConnell & Ma 2013).Since v ∼ (GM BH /p) 1/2 ∼ 50000 km/s, we can get X ∼ 80, which gives ξ(σ) almost equal to 1. Hence, by using the dark matter density spike expression, we can get the average energy loss rate for one period due to dynamical friction (including precession effect): As we mentioned that ρ sp and r sp depend on the spike index, from the above equation, we can see that the average energy loss rate due to dynamical friction depends on the spike index γ sp only.By writing the total energy decay rate Ė = ĖGW + ĖDF , we can constrain the value of the spike index γ sp .Including all of the uncertainties of the parameters and empirical relations, we can obtain the range of γ sp by adopting Ė = −(3.66± 0.24) × 10 41 W. As ĖDF depends sensitively on γ sp , from Fig. 1, we get a very narrow range of γ sp = 2.351 +0.032 −0.045 , which gives an excellent agreement with the canonical model prediction γ sp = 2.333.This means that the inclusion of the dynamical friction due to dark matter density spike can satisfactorily account for the discrepancy in the energy loss rate.Note that stellar heating effect near the primary SMBH might drive the spike index down to a smaller value (Gnedin & Primack 2004;Merritt 2004b).Nevertheless, if the central total stellar mass near the primary SMBH of OJ 287 is ∼ 0.1% of M BH , the heating time required to change the spike index would be longer than 17 Gyr (Merritt 2004b).Therefore, the stellar heating effect might not be significant for OJ 287.

Discussion
In this letter, we have shown that the observed energy loss rate (i.e. the period decay rate) of OJ 287 is much larger than the predicted energy loss rate solely due to GW emission, with a discrepancy of 4.3σ.Nevertheless, by adding the energy loss rate of the dynamical friction due to the dark matter density spike around the primary SMBH, it can satisfactorily account for the discrepancy and reproduce the observed energy loss rate, with a narrow range of the spike index γ sp = 2.351 +0.032 −0.045 .Surprisingly, this spike index gives an excellent agreement with the canonical model prediction γ sp = 2.333 based on the adiabatic SMBH growing model and the standard global NFW dark matter distribution predicted from numerical simulations.These provide a consistent picture to describe the period decay rate of OJ 287 and reveal the first robust evidence of the existence of a dark matter density spike around a SMBH.On the other hand, Alachkar, Ellis & Fairbairn (2023) recently followed the orbital dynamics to constrain the dark matter spike mass of OJ 287.Our result shows that the dark matter spike mass is about 0.1% of the primary SMBH mass, which is consistent with the constraints obtained in Alachkar, Ellis & Fairbairn (2023) (< 3% of the primary SMBH mass).
Note that the spike index relation γ sp = (9−2γ)/(4−γ) is derived based on the adiabatic growth model of a single SMBH (Gondolo & Silk 1999).For SMBHBs like OJ 287, such a relation might not be applicable.Therefore, the agreement between our constrained γ sp and the canonical value 2.333 may be just a coincidence only.According to numerical simulation results, SMBHBs can scatter dark matter particles and decrease the density in the inner regions (Merritt et al. 2002;Merritt 2013).Therefore, the spike index for OJ 287 would be smaller than the expected value.However, there are some mechanisms which can replenish the inner regions with dark matter particles and stars.For example, Zhao, Haehnelt & Rees (2002) show that the efficient randomization of the orbits can provide a replenishment of the dark matter loss near the SMBHs.Also, Beraldo e Silva et al. (2023) show that the orbits of dark matter particles and stars could be destabilized and brought to the inner galactic region with the crossing of the main bar resonances.These mechanisms provide some possibilities for replenishing dark matter particles to the inner region of OJ 287.Therefore, the large γ sp constrained in our study might still be possible.Hence, our result may provide an important clue to understand the complicated interactions and feedbacks between dark matter and SMBHBs.
In fact, Chan & Lee (2023) have shown the possible existence of dark matter density spikes around stellar-mass black holes.Our results here might support the idea of theoretical prediction that dark matter density spikes might exist around most of the black holes, including intermediate-mass black holes (Lacroix & Silk 2018;Chan 2018) and SMBHs (Gondolo & Silk 1999).In particular, the galactic SMBHs could be good targets for us to study the properties of dark matter, such as the annihilation and decay constraints.Besides, taking the SMBH in our galaxy as an example, dynamical studies using the stars orbiting the SMBH (Sgr A*) can be another way to examine the existence of any density spike around the SMBH (Chan, Lee & Yu 2022;John, Leane & Linden 2023;Shen et al. 2024).Future accurate observations of the stars orbiting the SMBH can help verify the dark matter density spike model and constrain the properties of dark matter.Moreover, future low-frequency GW observations in space can further examine OJ 287 and other similar SMBHBs in our universe.The pattern of the GW signals can reveal the structure of SMBHBs and the SMBHB inspiral process (Hannuksela, Ng & Li 2020;Zhao et al. 2023).The low-frequency GW data can provide the final smoking-gun evidence to verify our result and our understanding of the interactions between dark matter and SMBHs.

Fig. 1 .
Fig. 1.-The shaded region bounded by the black lines is the total energy decay rate constrained by observations (the thick black line indicates the average value).The shaded region bounded by the red lines is the predicted energy decay rate combining the effects of GW and dynamical friction (DF) of the dark matter density spike (the thick red line indicates the average value).The green dotted line indicates the energy decay rate due to GW emission only.The blue vertical dashed line indicates the spike index γ sp = 2.333 predicted by the adiabatic SMBH growing model.Dynamical friction in a dark matter density spike with spike index γ sp = 2.351 +0.032 −0.045 can account for the large orbital period decay rate measured by Dey et al. (2018).