Gravitational wave probes on self-interacting dark matter surrounding an intermediate mass black hole

The presence of dark matter overdensities surrounding a black hole can influence the evolution of a binary system. The gravitational wave signals emitted by a black hole binary offer a promising means to probe the dark matter environments near a black hole. The dense region of dark matter can lead to the dephasing of gravitational waveforms, which can be detected by upcoming experiments such as the Laser Interferometer Space Antenna (LISA). The dark matter density profile around the black hole can vary for different dark matter models. Our study specifically investigates the impact of the ultralight self-interacting scalar dark matter (SIDM) on the gravitational wave signals emitted by black hole binaries. A distinctive characteristic of SIDM surrounding a black hole, as opposed to collisionless dark matter, is the formation of a soliton core. We perform a Fisher matrix analysis to estimate the size of the soliton and the corresponding SIDM parameter space that future LISA-like gravitational wave experiments can explore.


I. INTRODUCTION
A wealth of compelling experimental data, such as the galaxy rotation curves and gravitational lensing observations, firmly supports the existence of dark matter (DM) [1][2][3][4][5].In pursuit of unraveling the nature of dark matter, extensive attempts have been undertaken including direct and indirect searches [6][7][8] as well as collider studies [9].However, despite these significant efforts, precise properties of DM such as its mass and interactions remain largely unknown.This enduring challenge of understanding DM stands as one of the longstanding conundrums in modern physics [10].Even though the cold dark matter (CDM), which serves as a fundamental component in the highly successful ΛCDM model, has played a remarkable role in elucidating a large-scale structure of the Universe [11], it encounters challenges such as the too-big-to-fail problem and the core-cusp problem when scrutinizing smaller scales ≲ 1 kpc [12][13][14][15][16][17][18][19][20].The null results from DM search experiments, combined with the presence of those unsettled small-scale structure issues, have spurred investigations into alternative avenues beyond the conventional CDM paradigm.
We study one such possible DM candidate, the ultralight self-interacting scalar dark matter, which is dubbed as SIDM for brevity .A prominent characteristic of such SIDM is a capability to form a stable configuration consisting of condensed bosonic fields, and the pres-ence of such soliton cores can alleviate the tensions on small-scale structures.By analyzing observations which can be affected by the existence of such a soliton, one can expect to extract valuable information about its fundamental properties.As a concrete example for such astrophysical observations, we explore a black hole binary system in the presence of SIDM and corresponding gravitational wave (GW) signals.Our goal is to investigate how future GW satellite experiments can be utilized to probe the properties of SIDM such as its mass and coupling.
For concreteness, our discussions focus on an intermediate mass black hole (IMBH) (typical mass range of order ∼ 10 2 M ⊙ − 10 5 M ⊙ which is the mass situated between stellar-mass black holes and supermassive black holes).The DM concentration in a vicinity of the IMBH is often referred to as "minispike" in contrast to the DM "spike" profile around a supermassive black hole (SMBH) at the center of a galaxy.IMBHs could form, for instance, in the centers of dwarf galaxies and globular clusters, and many IMBHs can well exist inside a single galaxy [45,46].Compared to SMBHs at the centers of galaxies, IMBHs are less likely to experience major mergers because of their association to smaller halos and galaxy masses, and yet surrounding DM profiles can be sufficiently dense to withstand tidal disruptions.The minispikes around IMBHs may hence well persist until the present epoch, which is a topic of active research.As an observational probe of such DM environments, we estimate GW signals originating from a black hole binary involving an IMBH.The GW frequency from the binary system depends on a mass ratio, and a stellar mass black hole orbiting the IMBH is of particular interest because the GW signals fall within the frequency band sensitive to planned space-based interferometers such as LISA (∼ 0.1mHz − 1Hz) [47].
The waveform of GWs is intricately linked to the properties of surrounding DM environments.In our investigation, we specifically focus on studying the effects of SIDM surrounding an IMBH.Notably, the presence of a solitonic core within a minispike marks this scenario as distinct from the vacuum case as well as from the conventional collisionless DM scenarios.The dynamics of a stellar black hole spiraling through such a DM environment are influenced by a dynamical friction and an accretion.These interactions result in accumulated phase shifts in gravitational waveforms, deviating from the predictions of the vacuum case.Such effects on the gravitational waveforms have been investigated actively for the collisionless DM spike profile surrounding a central black hole which an accompanying black hole moves through [48][49][50][51][52][53][54][55][56][57][58][59].More recently, the effect of a SIDM cloud has also been studied when a black hole binary system transverses through it, performing the waveform analysis in the time domain [60].Deviations in the speed of GWs that pass through the SIDM halo have also been studied to probe the mass and self-interaction [61].Further exploration and analysis of various DM models and configurations would provide an additional motivation for future GW experiments, enabling us to gain a deeper understanding of the nature of dark matter.
The aim of our study is to demonstrate that future GW experiments, such as LISA, have the potential to identify SIDM through its distinctive stable self-bounded structure in the central region of a black hole binary system.This study would contribute to shed light on the nature of dark matter and explore the unique signatures it may imprint on GW observations.Our paper is structured as follows.Section II overviews a soliton core formation around a black hole, which sets up the model to be investigated for the GW signal estimation.Section III then describes the evolution of a binary system in the presence of the SIDM soliton core.We discuss how a stellar black hole spiraling around a central IMBH can be affected by a dynamical friction and an accretion.Section IV presents our main results of the gravitational signal analysis.We clarify the SIDM parameter space which the planned satellite experiments such as LISA can probe, along with the accuracy of the DM parameter estimation from a Fisher matrix analysis.Section V is devoted to a discussion/conclusion.

II. SELF-INTERACTING DARK MATTER AROUND A BLACK HOLE
In this work, we consider a real scalar field ϕ SIDM model based on the following action: where m is a mass of the scalar field, and λ denotes a coupling constant responsible for a quartic self-interaction.
We assume that λ is positive so that the interaction is repulsive which can counterbalance the gravity to form a soliton inside a halo.The volume factor d 4 x √ −g is invariant under a general coordinate transformation where g = det(g µν ) denotes a determinant of the Friedmann-Lemaître-Robertson-Walker (FLRW) metric in Newtonian gauge Here, a(t) is a cosmological scale factor, and Φ denotes a scalar perturbation neglecting the contribution of an anisotropic stress tensor.We set c = ℏ = 1 in this section for the notational brevity.
In the nonrelativistic limit where the momentum of particles is of negligible magnitude, |⃗ p| ≪ m, we can rewrite the real scalar field ϕ in terms of a slowly varying complex field ψ and a fast-varying phase e ±imt : The time and momentum variances of ψ are much smaller than m, i.e. ψ/ψ ≪ m and ∇ψ/ψ ≪ m.Substituting Eqs. ( 2) and (3) into Eq.( 1), and neglecting fast oscillatory terms e ±2imt that average out to zero, we have where we have kept the terms up to the first order, and the superscript "•" denotes a differentiation with respect to time t.The equation of motion for ψ gives the nonlinear Schrödinger equation in the FLRW spacetime, where H = ȧ/a denotes a Hubble parameter.The nonlinear term arises from the self-interaction potential Φ self defined by On galactic scales, we assume that the system of scalar fields has been completely decoupled from the background evolution of the Universe, in which case we can neglect the scale factor and the Hubble parameter to rewrite Eq. ( 5) as When a high density medium of scalar fields condensates into the lowest momentum state, it behaves like a single macroscopic fluid exhibiting a superfluidity.The manyparticle system of Eq. ( 7) that involves a multitude of scalar fields is generically difficult to solve.Therefore, we pursue the macroscopic analysis of the superfluid [62][63][64][65][66], by using a mean-field approximation where ψ(⃗ r, t) denotes a wave function of the condensate, and δ ψ(⃗ r, t) is a small perturbation of the system.Then the density of the condensate is given by n(⃗ r, t) = | ψ(⃗ r, t)| 2 .This approach leads to the reduction of the many-body problem to one single field ψ(⃗ r, t) averaging out the effects of all other particles.In this description, we can decompose the condensate wave function using the Madelung transformation [67] where ρ(⃗ r, t) = n(⃗ r, t)m and s denote an amplitude and a phase respectively.Plugging Eqs. ( 8) and ( 9) into Eq.( 7) and taking imaginary and real parts result in hydrodynamical equations where we have defined represents a quantum pressure arising from the fact that a system cannot be of infinitesimal size due to the Heisenberg uncertainty principle.We will work in the regime where the quantum pressure is negligibly smaller than the self-interaction potential We consider a halo system comprised of a black hole (BH) of mass M BH at the center surrounded by the condensate of scalar fields (i.e.dressed BH).In the Newtonian limit, the metric perturbation can be written as where Φ BH = − rS 2r denotes a BH potential with Schwarzschild radius r S = 2GM BH , and Φ dress represents a potential generated by a gravitational interaction of the scalar fields satisfying a Poisson equation When a repulsive pressure due to the self-interaction balances out a gravity, the system can form a stable configuration with vanishing velocity ⃗ v = 0, in which case Eq. ( 11) reduces to Taking a divergence on both sides and assuming a spherically symmetric solution for ρ, we find where r c denotes a characteristic radius of a soliton defined by Changing a variable x = r/r c gives a spherical Bessel equation of order ℓ = 0, where there are two independent solutions, j ℓ=0 = sin x/x and n ℓ=0 = cos x/x.The final density profile of a soliton is given by a linear combination of these two solutions [25,27,41,68] where ρ sin and ρ cos are constants to be determined.The second term dominates for a small r where the gravitational pull from the black hole is significant.Outside of a soliton (r > r c ), on the other hand, the scalar fields start to behave like the collisionless CDM.To facilitate a smooth transition of a density profile from the outer to the inner region, we will use a continuity condition and a mass conservation to match profiles in adjoint sectors.
First, in the outermost region of halo, we consider the Navarro-Frenk-White (NFW) density profile [69]: where the concentration parameter c ≡ r vir /r sc is a ratio between a virial radius and a scale radius, and the halo density is ρ crit ∆ at the virial radius.We take ∆ = 200 and adopt the scaling relation between the concentration and the halo mass (c − M relation) given by Ref. [70] for concreteness.The critical density of the Universe today ρ crit = 2.775h 2 × 10 11 M ⊙ and the reduced Hubble constant h = H 0 /(100 km/s/Mpc) with H 0 = 67.4km/s/Mpc inferred from the CMB [11].Second, as they approach closer to the center, due to the influence of the central BH, the scalar fields are redistributed to form a spike profile [71,72] where ρ sp denotes a density at a reference radius r sp , and γ sp is a slope of the spike.Assuming an adiabatic growth of the spike in the presence of the central BH with an initial slope of γ i , it acquires a final slope of γ sp = (9−2γ i )/(4−γ i ).For instance, γ i = 1 for the NFW and then γ sp = 7/3.The two undetermined parameters ρ sp and r sp are obtained by a continuity condition and a mass conservation where r min = 3r S denotes an innermost stable circular orbit (ISCO), and the upper limit of the integral 5r sp is empirically obtained such that the resulting integral gives twice the mass of a BH [49]. 1 The integrand ρ DM is defined by Third, since the innermost region of halo is described by the soliton configuration in Eq. ( 19), we again use the continuity condition and the mass conservation to fix the parameters ρ sin and ρ cos , To sum up, the final density profile of a halo is composed of three layers Figure 1 shows the density profile of a halo around a 10 4 M ⊙ BH in the CDM model and SIDM model with different characteristic soliton radii for illustration.One can find that densities of soliton cores in the SIDM model are smaller than the density in the CDM model.This is because a greater repulsive self-interaction leads to a larger soliton core r c which reduces densities of solitons.In Sec.III, we will utilize this profile to compute a dynamical friction and an accretion rate for a stellarmass BH moving through an IMBH minispike.
In addition to the solitonic profile, another characteristic feature of SIDM to be compared with the collisionless DM is the effective sound speed due to the pressure caused by the self-interaction.Our analysis employs an effective model where the SIDM system reacts to a moving object by generating phononlike sound waves as follows.
Perturbations of the condensate are manifested by the excitation in the SIDM superfluid.They stand for sound waves that travel through the condensate.To obtain a sound speed at which it travels [65,66], we substitute Eq. ( 8) into Eq.( 7) with neglecting the contributions of the gravity and the perturbation where g = 3λ 4m 2 .The stationary solution which represents a ground state of the condensate takes the following form: where µ is a coefficient to be determined.Plugging this ansatz into Eq.( 29) and neglecting the kinetic term give The excited states can be found from seeking a solution of the form where Plugging this ansatz into Eq.( 29), and taking the Fourier transformation give Combining these equations gives the dispersion relation in k-space, where the sound speed is defined by The sound speed originates from the self-interaction, and it changes the energy spectrum of exited states.
In the small k limit (long-range), the first term in Eq. ( 36) dominates This implies that the repulsive interaction (λ > 0) results in an oscillating solution.However, the attractive interaction (λ < 0) gives either a growing or a decaying mode, indicating that it cannot form a stable condensate 2 .We also utilize this derived sound speed in Sec.III in studying the evolution of stellar-mass BHs traversing through a soliton core.

III. BLACK HOLE BINARY EVOLUTION
In the previous section, we have discussed the DM halo which surrounds a BH with a mass m 1 (≡ M BH ).When a stellar-mass BH with a mass m 2 is caught gravitationally by such a BH, the formed binary system can produce the observable GWs.The GWs from such a binary system carry the information of a surrounding DM halo; thus one can probe the property of the DM environment with these GWs [48,49].In this work, we focus on the case with q ≡ m 2 /m 1 ≪ 1 which can emit GWs in the frequency range of LISA.
With the density profile of the DM halo, one can derive dynamical equations of the binary system.The equation of motion in a radial direction is where r 2 is the distance between a stellar-mass BH and the barycenter of the system, r is the relative distance 2 See Refs.[26,38] for the case of the attractive interaction to form a stable soliton.
of two objects, ω s is the angular velocity, and M 1 is the total mass inside the orbit given by which includes the contributions from the central BH and the surrounding DM halo.
Since the mass ratio q is much less than 1, the orbit can be regarded as a quasicircular orbit and the DM halo can be considered as unperturbed. 3In this approximation, both ṙ2 and r2 vanish, and we obtain where M ≡ M 1 + m 2 is a total mass of the system.
In addition to the gravitational pull of DM inside the orbit, the DM halo can also affect the evolution of the binary system by a dynamical friction and an accretion.Since we consider a small stellar-mass BH orbiting a large central BH with its own unperturbed DM halo, only the dynamical friction and the accretion affecting the dynamics of the stellar-mass BH are taken into account.
The combined effects of GW emissions, the dynamical friction, and the accretion decrease the orbital energy and reduce the relative distance of the binary.The evolution of relative distance reads as where µ ≡ M 1 m 2 /M denotes a reduced mass, and the backreaction force due to the emission of GWs is given by [76] where v ≡ rω s is the velocity of a stellar-mass BH.
The dynamical friction of a stellar-mass BH orbiting inside the SIDM medium takes the form [77] where the function I is defined by where Λ ≡ vt/r min with t denoting the time for which the stellar-mass BH has traveled.The Coulomb logarithm ln Λ incorporates the minimum and maximum impact parameters of the stellar mass BH.The Mach number M ≡ v/c s parametrizes how fast the stellar-mass BH moves through the medium with respect to the sound speed c s defined by Eq. ( 37) with substituting the density profile with Eq. ( 28).For the supersonic expression, in the limit M ≫ 1, we recover the conventional steady state collisionless dark matter result I → ln(vt/r min ) with the replacement of vt → r max .
Here we use the Coulomb logarithm ln Λ = ln m 1 /m 2 which is commonly adopted in the literature [55, 78, 79] 4 .
In the subsonic limit M ≪ 1, on the other hand, the dynamical friction exerted by the collisionless medium is larger than the SIDM partly due to the self-interacting repulsive force that hinders the accumulation of the DM in the wake.The F Ac denotes the accretion drag force where μ = ṁ2 (1 + q) −2 , and ṁ2 denotes the accretion rate.For the SIDM medium, we adopt the Bondi-Hoyle-Lyttleton accretion rate [51,[82][83][84][85] while for the collisionless medium we use the form In the transonic regime v ∼ c s ≪ c, for example, the accretion rate in the SIDM medium is much larger than the collisionless case Such an enhancement can arise from "funneling effect", in which the particles tend to funnel into inward-spiraling trajectories [85].The collisions and interactions among 4 Such a choice of ln Λ would suffice for our purpose of demonstrating the potential gravitational signals, partly because of its weak logarithmic dependence and also other relevant uncertainties such as the modeling of GW detector noise and sensitivity.The exact value of the Coulomb term Λ is not well known and is usually obtained numerically in the literature.For instance, r min can be estimated by demanding that the numerically calculated force should match the analytically estimated one [80].
The value of rmax can be, for instance, the radius under the influence of BH gravity (Roche radius or Hill radius) or the soliton radius.Different authors use different values for the Coulomb term, and we refer the readers to the existing literature (e.g.[42,55,60,77,80,81]) for more detailed discussions on the dynamical friction in the presence of the gaseous medium and collisional DM.
the particles redistribute their angular momenta, and, combined with the gravitational influence of the BH, the particles tend to follow more radial trajectories toward the accreting BH.The relative magnitude of the accretion drag force and the dynamical friction for the collisionless medium is and for the SIDM case, in the supersonic limit v ≫ c s , In comparison to the binary in the vacuum, the binary in the DM halo can be also affected by the dynamical friction and accretion.These effects accelerate the loss of orbit energy of the binary; thus the number of orbital cycles before the coalescence N cyc gets smaller.
In Fig. 2, we show the evolution of forces contributed by GWs, dynamical frictions, and accretions respectively where τ ≡ t| r=rISCO − t is the time to ISCO, and benchmark parameters of the binary system are given by the set 1 in Table I.One can find that, in the final stage of inspiral which is relevant for GW observations, the force due to the emission of GWs is dominant for both CDM and SIDM halos; thus the orbital decay of the binary is still mainly determined by the emission of GWs, rather than by the dynamical friction or accretion.Figure 3 demonstrates the corresponding difference of N cyc compared to the binary in the vacuum, ∆N cyc = N cyc (vacuum) − N cyc (with DM halo).One can find that ∆N cyc is primarily contributed by the dynamical friction (represented by dot-dashed lines), while the contribution from the accretion (dashed lines) is subdominant.Even though the force due to the GW emission is significant and hence can affect how quickly the stellar-mass BH approaches toward the central BH, its contribution to the accumulated phased shift in the waveform is negligible.

IV. GRAVITATIONAL WAVE ANALYSIS
The waveform of GWs produced from the inspiral of a binary system is by [76] h where D L is a luminosity distance to a source, M c = µ 3/5 M 2/5 denotes a chirp mass, t ret = t − D L /c is a retarded time, ι is the angle between an orbital angular momentum axis of a binary and a direction to a detector, f is a frequency of GWs which is given by f = ω GW /(2π) with ω GW = 2ω s , and Ψ is a phase of GWs. Figure 4 shows the difference in the phase of GWs, ∆Ψ = Ψ(vacuum) − Ψ(with DM halo) (55) for different DM parameter sets compared with the vacuum case, and the parameters of the binary system are given by the set 1 in Table I.One can find the phase of the GW waveform is affected, the so-called dephasing effect, by the DM environment.Over the inspiral period, this phase shift effect can accumulate and result in observable signals of sufficient magnitude to probe the characteristics of the dark matter environment.Even though the dephasing effect is calculated with a numerical approach in our quantitative discussions, it would be informative to provide an analytical assessment to see how it is affected by our model parameters.The conservation of energy gives where the orbital energy of the binary is and the power of energy loss due to the GW emission, the dynamical friction, and the accretion are respectively given by To make a rough estimation, we consider the supersonic regime v ≫ c s , neglecting the dependence of a sound speed in the F DF and F Ac .We also assume that ρ halo is a constant.With using v = rω s and Eq. ( 41), we can rewrite Eq. ( 56) as with Here the dimensionless quantities y and x are defined as x ≡ τ ω GW,ISCO , with The difference in the phase of GWs can be written as The LISA is in a heliocentric orbit and consists of an equilateral triangle formed by three spacecrafts, each separated by a distance of 2.5 million kilometers from one another.The center of mass for the constellation, known as the guiding center, is in a circular orbit at 1 AU and 20 • behind the Earth.Choosing the polar coordinate system with the Sun at its origin, the strain of GWs in a detector is given by [86] where F + and F × are the detector response functions, ϑ and φ are the latitude and longitude of the binary in the polar coordinate system, and χ is the polarization angle.
Here, ∆t is the delay between the arrival time of GWs at the Sun and the arrival time at the detector, which is given by The detector response functions are given by with Here, α = 2πt/yr + α 0 is the orbital phase of the guiding center, and β = 2πn/3 + β 0 (with n = 0, 1, 2 for three spacecrafts) is the relative phase of the spacecraft within the constellation.The parameters α 0 and β 0 give the initial ecliptic longitude and orientation of the constellation.
In the limit of a large signal-to-noise ratio (SNR), the posterior probability distribution of the source parameters can be approximated by a multivariate Gaussian distribution centered around the true values.The corresponding covariance can be estimated by the inverse of the Fisher information matrix.For a network including N independent detectors, the Fisher matrix can be written as where d is given by S 2 (f ) , . . ., hN (f ) and θ denotes the parameter vector with true value θ.
Here S i (f ) is the noise power spectral density for the ith detector and hi (f ) is the Fourier transformation of the time domain signal.The bracket operator (A, B) for two functions A(t) and B(t) is defined as where A * is its complex conjugate.The total SNR is given by (d, d).
The root-mean-squared (1σ) errors of parameters and the correlation coefficients among parameters can be estimated by the inverse of the Fisher matrix Σ = Γ −1 , which can be written as and For the binary with SIDM halo, the parameter vector is θ = {r c ; m 1 , m 2 , D L , ι, χ, ϑ, φ, ϕ ISCO , t ISCO }, where ϕ ISCO and t ISCO are the source phase and time respectively at ISCO.We consider the binary system with masses given in Table I, {ι, χ, ϑ, φ} are set to be π/4, {ϕ ISCO , t ISCO } are set to be 0, and D L is set to vary the SNR. Figure 5 shows the probability distribution of {r c , m 1 , m 2 } with the fiducial parameter values θ given by set 1 in Table I and r c = 100r S .The tilded parameter is defined as the value with respect to the fiducial value adopted in the Fisher matrix analysis θ = θ/ θ.The dotted vertical lines in the one-dimensional likelihood function marginalized over the other parameters indicate the 1σ interval for the parameters.Two-dimensional ellipses are the 1-, 2-, and 3-σ confidence contour plots.The correlation values defined by Eq. ( 79) are also shown in the figure and r c is highly correlated with m 1 and anticorrelated with m 2 .This correlation arises because the increase in r c results in the decrease in the soliton density amplitude, and such an effect on the GW signals can be compensated by the increase (decrease) of the central (accompanying) BH mass.0 .9 9 9 9 5 1 .00 0 0 0 1 .00 0 0 5 0 .98 0 0 .9 9 0 1 .00 0 1 .0 1 0 1 .0 2 0 0 .9 9 9 9 6 0 .9 9 9 9 8 1 .00 0 0 0 1 .00 0 0 2 1 .00 0 0 4 0 .9 9 9 9 5 1 .00 0 0 0 1 .00 0 0 5 0 .9 9 9 9 6 0 .9 9 9 9 8 1 .00 0 0 0 1 .00 0 0 2 1 .00 0 0 4 Focusing on the DM related parameters, the error of r c (∝ √ λ/m 2 ) is shown in Fig. 6.The existence of the soliton core is the characteristic feature of SIDM, and the SIDM parameter space for which the error on r c is small enough is of our particular interest.For the parameter range of our interest, the stellar BH is at a radius around O(1 ∼ 10)r S a half year before it reaches the ISCO.For instance, for {m 1 /M ⊙ , m 2 /M ⊙ } = {10 4 , 1} illustrated in Fig. 6, a half year observation of LISA corresponds to receiving the signals when the stellar BH moves though the SIDM soliton from r ∼ 18r S to r ∼ 3r S .The stellar BH is hence well inside the soliton core during the observation, and a bigger soliton core density leads to a bigger dynamical friction and accretion.We can hence expect that the signal can be more sensitive to the parameter space which leads to a bigger soliton core density.Our modeling of the DM profile around a central BH in Sec.II has led to a bigger soliton density for a smaller soliton radius and the accuracy is indeed better for a smaller soliton size r c . Figure 6 shows that a half year of LISA observation is expected to measure the characteristic radius of soliton r c within 1% (10%) accuracy when r c ≲ 444r S (r c ≲ 2700r S ).The binary with a smaller total mass M and a larger mass ratio q increases the expected accuracy for a given r c /r S because it results in a bigger dephasing effect.
We can also argue how the black hole masses affect the estimated parameter accuracy illustrated in Fig. 6 by looking at their effects on the dephasing.The analytical expression of Eq. ( 68) also helps in understanding the behaviors of Fig. 6.The difference in dephasing with different values of m 2 is partly due to the difference in the dynamical friction which becomes bigger for a bigger m 2 when the central IMBH mass is fixed.The dynamical friction is also sensitive to the velocity of the stellar mass BH and equivalently to the frequency of the orbit.Such a dependence of the dynamical friction on the frequency is significant as inferred from v −2 dependence in Eq. (44).Consequently, the dephasing for the m 1 = 10 5 M ⊙ example becomes smaller compared with the m 1 = 10 4 M ⊙ examples in Fig. 6 because a larger central mass leads to a larger initial orbital radius and slower orbital motion.
Based on the results in Fig. 6, we can obtain a detectable region (orange hatched) in the parameter space FIG. 7. The detectable region (orange hatched) in the parameter space of {m, λ} from a half year LISA observation of the binary system by requiring σr c /rc < 1.The dashed purple line corresponds to the limit rc ≫ λ dB , and the bullet cluster constraint σ/m ≲ 1cm 2 /g is represented by the dashed blue line.The dashed green line represents the Jeans scale bound rJ ≲ 1 kpc.Further information is provided in the main text.
of {m, λ} by requiring σ rc /r c < 1 in Fig. 7.The dashed purple line corresponds to the limit where the characteristic soliton scale is larger than the de Brogile wavelength r c ≫ λ dB to avoid the effect of quantum pressure [42].The dashed blue line corresponds to the bullet cluster constraint σ/m ≲ 1cm 2 /g.To form enough subgalaxy scale structures, we demand a Jeans length should not exceed 1 kpc scale r J ≲ 1 kpc (dashed green line) [24].The Jeans length of the SIDM model with a quartic self-interaction is determined by a characteristic soliton scale r J = r c (which is obtained by balancing the self-interacting repulsive force and the gravitational attraction force) [24,40].
Therefore, our gravitational probes on the SIDM model will be able to shed light on the uncharted area of the parameter space that has not been surveyed by other measurements.However, the limiting factor of our study is the presence of a degeneracy between two parameters λ and m arising from their underlying connection to r c (∝ √ λ/m 2 ).To break this degeneracy, one would need other probes which possess different dependence on the DM parameters, for instance, the observables which are not directly linked to the height of a soliton profile.Such an inquiry is left for future work.

V. CONCLUSIONS
We have studied how SIDM can be probed by forthcoming GW experiments.It would be ideal to identify a region with a large DM density for effective DM exploration, and the overdense region surrounding a BH pro-vides a promising environment for studying DM properties.Through the Fisher matrix analysis, we have demonstrated that SIDM scenarios can be distinguished from the vacuum scenario, and we clarified the previously unexplored SIDM parameter space which can be probed by the planned GW experiments.The future space-borne GW detectors hence would offer a promising means to study and identify DM, which can complement ongoing terrestrial DM search experiments.While we focused on the LISA specifications in our analysis, it has been proposed that space-based GW detectors, such as LISA and Taiji, can form an observation network [87].With the aid of such a network in the future, the observable parameter space can be further enlarged.The anticipated data from future GW experiments hold great promise in shedding light on the nature of DM and providing insights into the properties of SIDM.
Before concluding our discussions, it is worth mentioning that a more precise estimation of the dynamical friction and accretion processes in the presence of SIDM around a BH deserves further exploration beyond the simplified estimation used in our current modeling.Additionally, investigating the DM profile of SIDM around a BH, which extends beyond the scope of our modeling presented in this paper, represents an intriguing topic for future research.For instance, in our study, we employed the Schrödinger-Poisson equations to describe the DM profile around the BH, treating it as a source of fixed, external gravitational potential.This approach provides a reasonable starting point for exploring GW probes, especially considering that the dephasing of the waveform predominantly originates from the radius farther away from the Schwarzschild radius due to the larger number of orbital cycles experienced by the stellar BH.However, further investigations with more precise treatments, such as utilizing the Einstein-Klein-Gordon equations, are warranted and left for future work.

FIG. 1 .
FIG. 1.The DM density profile around a 10 4 M⊙ black hole in the CDM model and SIDM model with different soliton characteristic radii.

34 FIG. 2 .
FIG. 2. The evolution of the forces contributed by GWs, dynamical frictions, and accretions in the CDM model and the SIDM model with different soliton characteristic radii.τ is the time to reach the ISCO.The force due to the gravitational wave emission is dominant, and three GW curves are indistinguishable among different parameter sets in this figure (represented by the top dotted curve).

9 FIG. 3 .
FIG. 3.The difference in the number of orbital cycles compared with the binary in the vacuum, ∆Ncyc = Ncyc(vacuum)−Ncyc(with DM halo), as a function of the time τ to reach the ISCO.

FIG. 4 .
FIG.4.The dephasing of the gravitational waveform for our example parameter sets with respect to the vacuum case, ∆Ψ = Ψ(vacuum) − Ψ(with DM halo) as a function of the GW frequency.The vertical dashed line denotes the frequency at the ISCO.The corresponding τ (the time to reach the ISCO) at each frequency is also shown.

1 FIG. 6 .
FIG.6.1σ uncertainty of the SIDM soliton radius σr c as a function of rc (x axis is scaled to the Schwarzschild radius and y axis is scaled to a given fiducial value in the Fisher matrix analysis).The SNR is set to 100 and 10 for illustration and a half year of observation by LISA is assumed.To illustrate the soliton size of interest, the gray dotted lines indicate the relative accuracy of 0.1 and 0.01.

TABLE I .
The benchmark parameter sets we use in the analysis.The mass m2 changes with time due to the accretion, and the values here are given when r = rISCO.