The effect of wave dark matter on equal mass black hole mergers

For dark matter to be detectable with gravitational waves from binary black holes, it must reach higher than average densities in their vicinity. In the case of light (wave-like) dark matter, the density of dark matter between the binary can be significantly enhanced by accretion from the surrounding environment. Here we show that the resulting dephasing effect on the last ten orbits of an equal mass binary is maximized when the Compton wavelength of the scalar particle is comparable to the orbital separation, $2\pi/\mu\sim d$. The phenomenology of the effect is different from the channels that are usually discussed, where dynamical friction (along the orbital path) and radiation of energy and angular momentum drive the dephasing, and is rather dominated by the radial force (the spacetime curvature in the radial direction) towards the overdensity between the black holes. Whilst our numerical studies limit us to scales of the same order, this effect may persist at larger separations and/or particle masses, playing a significant role in the merger history of binaries.

Introduction.-Gravitational-waveobservations provide a unique window that can be used not only to infer the astrophysical properties of black holes (BHs), but also to gather information about the environments they live in.The presence of matter around BHs during a binary merger event results in modifications to the trajectories, which in turn changes the gravitational-wave signal in a characteristic way [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19].Environments may arise from standard matter, such as accretion discs, or from dark matter (DM).In this work we focus on the latter case.
The DM energy densities required to give significant effects on the signal are high relative to the expected average galactic values, with the latter determined by large scale observations [20][21][22][23][24]. Therefore the impact of such effects may be expected to be small [6].However, average galactic densities describe DM on large scales only, and its distribution on small scales (in particular the parsec and below scales relevant for astrophysical BHs) is not well constrained [25].There exist several mechanisms that could create DM overdensities around isolated BH.One well known possibility is the superradiant instability, in which a bosonic field extracts energy and angular momentum from a highly spinning black hole via repeated scattering in the ergoregion [26][27][28][29][30][31] (see [32] for a review).Another more prosaic effect is simply the accretion of dark matter in the potential well around BHs, which results in the formation of "dark matter spikes" [33] (a combination of both superradiance and accretion may lead to even higher densities [34]).Such spikes were originally proposed in the context of WIMP-like dark matter, but in general their profile is a power law with an exponent that depends on the effective equation of state of the dark matter [35][36][37][38][39][40][41].However, they also occur for low mass, wave-like DM candidates, with a form that is dependent on the relative Compton wavelength of the DM particle and the black hole horizon [42][43][44][45][46][47].In both cases, the DM density near the BHs depends on the asymptotic dark matter environment and on the particle properties.* josu.aurrekoetxea@physics.ox.ac.ukFIG. 1. Dephasing in the coalescence time ∆t c for a 10-orbit binary for different scalar masses µ and initial densities ρ 0 .Here tc is the merger (coalescence) time in the absence of a dark matter cloud.The effect is maximized for µ ≈ 0.45M −1 , corresponding to a Compton wavelength of the dark matter particle that is comparable to the initial separation of the orbit However, a key question is whether these overdensities around isolated objects persist during a binary merger event.
In this Letter, we build on our study to better understand how generic such an effect is, and to properly quantify the impact that the DM has on the binary close to merger.We focus on the effect of wave DM on equal mass BH mergers, and in particular its dependence on the mass of the scalar particle.We simulate a 10-orbit binary black hole in an initially homogeneous dark matter environment starting from initial conditions satisfying the Hamiltonian and momentum constraints.We identify the decay of the orbit (and, as a consequence, dephasing of the gravitational wave signal) as being a direct result of the scalar cloud.Our key results are illustrated in Fig. 1, where we show the dephasing is maximized when the mass of the scalar particle is such that its Compton wavelength is comparable to the initial separation of the orbit λ c = 2π/µ ∼ d 0 .In addition, we are able to quantify the different channels that contribute to the dephasing in our scenario, finding the dominant effect to be driven not by radiation or dynamical friction drag forces, as are often discussed, but rather the attraction of the binary to the central overdensity.
Key background and physical setup.-Weconsider a minimally coupled massive complex scalar field Φ described by the action with a quadratic potential V (Φ) = µ 2 Φ * Φ/2, where µ is the scalar field mass.The dynamics of the scalar field is given by the Klein-Gordon equation coupled to gravity In the case of a single BH immersed in a reservoir of such scalar DM, the stationary solution near the black hole is described by the Heun functions [43,[110][111][112], with a power law envelope and characteristic oscillations in the spatial profile on length scales set by the scalar wavelength.In the case of a binary no analytic form for a stationary state is known, but simulations [104] using the numerical codes GRCHOMBO [113] and GRDZHADZHA [114] have indicated that for a range of initial configurations and within a few orbits, the scalar matter evolves into a persistent quasi-stationary profile with density spikes near the black holes and an overdensity between them.
Ideally we would set this "natural" quasi-stationary DM configuration as an initial condition, and study the impact the cloud has on the binary merger using general relativity.However, even if an analytic form was known, a consistent solution of the GR constraints would lead to changes to the initial effective masses and momenta of the black holes for different densities and profiles, making comparisons of the subsequent evolutions difficult to interpret.In particular, it is hard to know if the additional dephasing is arising due to matter effects or due to the increased initial eccentricity of the orbits.One can mitigate this by applying eccentricity reduction schemes to the initial data, but the fact that the clouds can be very dense We refer to these as the "cloud" and "no cloud" cases, respectively.Simulation movie [109].
near the horizon makes this challenging as the eccentricity is extremely sensitive to small changes.
In this Letter we take a simpler approach.Given the short relaxation timescale of the cloud (∼ 2 orbits), compared to the timescale of the merger we are simulating (∼ 10 orbits), we start all simulations from a homogeneous configuration with fixed initial density ρ 0 = µ 2 φ 2 0 and allow the cloud to build up dynamically during the simulation.To do so, we choose homogeneous initial conditions for the real and imaginary components of the scalar field, Φ = (φ 0 , 0) and ∂ t Φ = (0, µφ 0 ).As we vary µ (which gives us different cloud configurations, Fig. 2), we adjust φ 0 so as to we keep the initial density ρ 0 unchanged, thus the initial trajectories of the binary as a result of solving the initial data will be the same in all cases, which allows comparison between different masses (and different scalar cloud profiles) in a more controlled way.See the Supplemental Material for more details about the numerical implementation of the evolution equations and initial data.
Dephasing of the binary.-Westudy the tensor gravitational-wave modes emitted by the binary black hole extracting the Newman-Penrose scalar Ψ 4 with tetrads proposed by [115], projected into spin-weight −2 spherical harmonics, Ȳ lm dΩ, where dΩ = sin θ dθ dϕ is the area element on the S 2 unit sphere.The merger or coalescence time for our 10-orbit binary in vacuum is tc ≈ 2000 M, defined as when |ψ 22 | peaks, which is the dominant mode.This is also the case for small initial densities ρ 0 since there is less backreaction of the matter on the binary metric.For a given density, a smaller effect is also seen for masses µ ≪ M −1 , as the FIG. 3. Dephasing of the gravitational wave signal due to the accretion, dynamical friction and emission of wave dark matter around a binary black hole merger.Top panel is the real part of the ψ 22 mode, whilst mid and bottom panels are its modulus and phase, respectively.The black solid line corresponds to (µ, φ 0 ) = (68, 50) × 10 −4 , which we refer to as "no cloud", as the Compton wavelength of the scalar field is much larger and we do not efficiently excite a DM cloud.The blue solid line corresponds to (µ, φ 0 ) = (4300, 0.79) × 10 −4 , and causes a ∆t c /t c ≈ 10% dephasing of the merger time.The initial densities for both these cases is ρ 0 ≈ 10 −9 M −2 .Movie in [109].
Compton wavelength λ c ≫ d and the cloud is not efficiently excited, see Fig. 2. We refer to the case of small µ as the "no cloud" configuration (it has the same initial, non zero density, but no structure forms around the binary), and the higher µ case as "with cloud".A typical result can be seen in Fig. 3.
As expected, the presence of wave dark matter around the binary results in a dephasing of the gravitational-wave signal ∆t c ≡ t c − tc , which is caused by effects like accretion and dynamical friction from the cloud.We will give explanations and order of magnitude estimates for the various effects in the following section.
In Fig. 1 we compare the dephasing for different DM masses µ ∈ {0.0068, 0.86}M −1 , corresponding to wavelengths λ c ∈ {924, 7}M that span a range above and below the initial binary separation d 0 ≈ 12 M.We find that the dephasing is maximized for µ ≈ 0.45M −1 , corresponding to If the mass is smaller µ < 0.45M −1 , the cloud quickly becomes suppressed and the dephasing becomes negligible, but we note that larger separations earlier in the lifetime of the binary may support clouds at smaller masses.If the mass is larger µ > 0.45M −1 , the dephasing is smaller but remains significant, and we still find an efficient excitation of the cloud.On the one hand, this is not so surprising -even at our highest mass, we are still in a regime where µ ≈ M −1 , and so as the merger radiates gravitational waves and inspirals in, it eventually approaches an orbital separation comparable to λ c .On the other hand, in other studies (see e.g.[84]) one of-ten finds that the behaviour at this limit is already reasonably well described by the particle limit, and so we might have expected to see a greater dissipation of the cloud and suppression of the effect.The fact that this is not the case implies that the mass does not need to be very finely tuned for the effects to be significant, and motivates a more detailed study to find the boundary between the wave and particle regimes.
We also vary the asymptotic energy density to find the value at which the dephasing is detectable in our simulations during the last 10 orbits, which gives an indication of the value required for effects to be significant at merger.See the conclusion section for these values in physical units.
Quantification of the causes of the dephasing.-Toquantify the origins of the dephasing we identify the changes in the energy, angular and radial momentum of the binary that relate to the presence of the cloud of matter.We follow the approach of [84,[116][117][118], and define a current J µ = ξ ν T µ ν in the direction ξ ν and associated charge and flux where N i is the outward normal direction to the surface that bounds the volume.If ξ ν is a Killing vector then ∇ µ J µ = 0 and the change in charge is balanced by a flux through a surface.When this is not the case the conservation laws require an additional "source" term describing the exchange of the charge between matter and curvature.It is this quantity that corresponds to gravitational forces in the Newtonian limit, and that quantifies the way in which momentum is extracted from the binary by the matter 1 .
The explicit expressions for the energy, angular and radial momentum charges, sources and fluxes in terms of the ADM variables are given in the Supplemental Material.In Fig. 4 we plot the time integration of each of these quantities and verify their agreement.In each case the black line should equal the sum of the blue and red lines and provides a check on the error.It is the red line that quantifies the exchange of the relevant charge from the binary to the matter, and therefore that drives the dephasing.
In the top panel of Fig 4, we see that the energy of the matter within the volume increases, due to the flux of matter energy across the outer surface -this is simply reflecting the fact that the central cloud density grows over time due to accretion from the environment.The increase is partially offset by a negative source term, which is mainly driven by the accretion of energy into the BHs, increasing their masses by approximately 1% over the course of the merger (we can check that this agrees to the change in their measured masses from the apparent horizon finder).After the merger, the energy in the −1

DM cloud
Source from/to curvature Flux from/to infinity −3 Conservation law for the energy, angular (corotating) and radial (inward) momentum of the matter in a sphere of radius 40M around the binary that contains the main cloud overdensity (Fig. 2).The black line shows the change in each respective charge for the cloud over time and is the sum of the other two lines.The red line describes the time integration of the exchange of the charge from the curvature to the matter.The blue line is (minus) the total flux of the charge into the outer bounding surface.
cloud around the remnant decreases slightly as some is accreted onto the remnant, but it does not completely dissipate.
In the second panel, the angular momentum held in the matter cloud initially increases as the curvature of the binary "stirs up" the cloud during the transient phase, then reaches a reasonably steady state during which the rate of extraction of angular momentum from the spacetime curvature (the stirring) is balanced by its flux out of the outer surface.The result is that the angular momentum of the binary is decreased, and carried away by scalar waves.We can view the source from/to curvature as a dynamical friction effect -the extraction of the angular momentum of the spacetime of the binary by the matter.After merger there is an increased flux of radiation from the outer surface of the volume which carries away all the angular momentum built up in the cloud during the inspiral.
In the third panel, the radial momentum held in the matter cloud is tracked.Here, we see that the matter cloud overall gains some inward radial momentum.However, the accretion of inward radial momentum from the outer surface is roughly balanced by its loss into curvature -partly as a result of the binary accreting the ingoing momentum, and partly as a result of it being attracted to the central overdensity.This gives the binary an inward pull that accelerates during the final plunge.
These measurements allow us to quantify which of the three contributions accounts for most of the dephasing.For circular orbits the angular momentum is given by J = Mr 2 (2π/T ) and very roughly we can say that the total dephasing ∆t c /t c ∼ n∆T 0 /(nT 0 ), where n is the number of orbits and T 0 is the pe-riod of the first orbit.The changes in T can be estimated in terms of the total changes of mass, angular momentum and radius (expressed in terms of the final amplitudes of the red lines in Fig. 4) as where we have used ∆r/r 0 ≈ ∆P tot r t2 c /(r 0 M 0 tc ) in the last term, obtained from integrating the inward radial force on each BH over the inspiral time tc .Using that r 0 = d 0 /2, n ≈ 10, tc ≈ 2000M and the values observed in Fig. 4, we estimate O(0.1)%, O(0.01)% and O(10)% for each term, respectively.We conclude that the radial force towards the overdensity is the dominant cause of the dephasing.
Conclusion.-Using general relativistic simulations of a binary accreting dark matter, we have shown that the dephasing in the gravitational-wave signal of an equal mass black hole merger is maximized when the Compton wavelength of the dark matter particle is comparable to the orbital distance of the binary, 2π/µ ∼ d.We need the mass of the scalar to be sufficiently large for a central overdensity to build up -low mass scalars suppress structure on smaller scales than their Compton wavelength.Converting into physical units, the optimal scalar mass to induce dephasing in the last 10 orbits of an equal mass binary with total mass M is then µ ≈ 5 × 10 −17 M/10 6 M ⊙ −1 eV, which can result in a 10% dephasing during the last 10 orbits of the binary (taking the blue line in Fig. 1) for asymptotic densities around the BH of This is high relative to the average DM density, but our measured dephasing is only over a short period of the binary's lifetime (∼ 10 orbits), and has a cumulative effect.Therefore smaller densities could give sufficient dephasing to be detectable assuming the effect is triggered at larger separations (which would also allow lower mass candidates to contribute to the effect), or if observations can happen over a longer time frame (e.g. by combined LISA/LVK observations).The simulations in this work demonstrate that accumulation of wave-like dark matter between the binary could have a significant effect on the merger history of binaries, unlike in particle cases where the dark matter tends to disperse.As recently suggested in [74], it could even go some way to explaining the final parsec problem.In particular, we highlight the importance of considering the radial force arising from any central overdensity that forms, in addition to the radiation of waves carrying angular momentum and energy away from the binary.As noted above, the effects remain significant even at the higher end of the masses that we can probe in our simulations, at which µM ≈ 1.Further investigations should be made to determine the point at which particle-like behaviour takes effect, and to study the importance of the relativistic features in our simulations such as the presence of black hole horizons.

NUMERICAL IMPLEMENTATION, DIAGNOSTICS AND CONVERGENCE TESTS
We evolve the gravity sector solving the Einstein field equations for a line element that we decompose in the usual ADM form [119] ds 2 = −α 2 dt 2 + γ i j (dx i + β i dt)(dx j + β j dt), (7) where γ i j is the three-dimensional spatial metric that we decompose into a conformal factor and a conformally related metric γ i j = γi j /χ.The lapse and shift gauge functions α and β i determine the choice of spatial hyperslicings and their coordinates, which in numerical relativity are dynamically determined.The extrinsic curvature tensor K i j = (2D (i β j) − ∂ t γ i j )/2α is decomposed into a trace K and a traceless part A i j , i.e.K i j = A i j + (1/3)Kγ i j .We evolve use the CCZ4 formulation [120] and the moving puncture gauge [121][122][123][124] with GRCHOMBO [113,125,126].We solve the Hamiltonian constraint using Bowen-York initial data using the CTTK hybrid method [127] (see table I for the binary parameters).In the homogeneous limit, this reduces to choosing the trace of the extrinsic curvature tensor K 2 = 24πGρ and solving for a correction of the conformal factor χ sourced by the traceless components A i j .We choose K < 0 so that the scalar field is initially decaying, which chooses the more conservative impact on the merger.The value is small and the effect of either choice on the overall trends observed is not significant.We use a simulation box length L = 512M and 9 levels of mesh refinement (See Fig. 5 for convergence tests).We impose reflecting boundary conditions at z = 0, while for the other boundaries we impose either zeroth order extrapolating boundary conditions (matching the values on the exterior ghost cells to the value (radially directed) inside the simulation grid) or Sommerfeld boundary conditions.
Following the approach of [116,117], the energy momentum tensor is decomposed into the energy density, momentum density and stress-energy density measured by normal observers as The conservation law is written as where these correspond to the change in charge in the cloud, the exchange between matter from/to curvature, and the flux to/from infinity.The energy and angular momentum currents are related to the time-like, angular and radial directions J F t = N i β i (αρ − β j S j ) + α(β k S i k − αS i (11) F {φ ,r} = −N i β i S j ξ j {φ ,r} + αN i S i j ξ j {φ ,r} , where N i = (x, y, z)/r is the normalised radial unit vector, with s i = (x, y, z)/r and N i = s i / (γ jk s j s k ).The source terms where ∂ t γ i j = −2αK i j + D i β j + D j β i , and both ∂ t α and ∂ t β i are given by our moving puncture gauge conditions.The quantities ∂ µ ξ ν are all zero except ∂ x ξ y = 1 and ∂ y ξ x = −1 We also track the flux through inner surfaces that move together with the black holes, which introduces additional advection terms to the flux where N BH i is defined above.

FIG. 2 .
FIG.2.Cloud density for two values of µ: one large with λ c ≈ d 0 (right) in which case the binary obtains an enhanced density cloud, and one smaller, λ c ≫ d 0 (left) in which case the pressure coming from the long wavelength of the collective excitations of the field prevents a high density cloud from forming.We refer to these as the "cloud" and "no cloud" cases, respectively.Simulation movie[109].
FIG.5.Convergence testing the gravitational-wave phase evolution for the largest dephasing case: (µ, φ 0 ) = (4300, 2.5) × 10 −3 , so that the initial density is ρ 0 ≈ 10 −8 M −2 .We use three different resolutions with N 3 number of grid points on the coarsest level.The decrease in the error is consistent with 1 st order.
0, −y, x, 0) and ξ ν r = −(0, x, y, z)/r.The respective charges Q and fluxes F in terms of ADM quantities are then