{'fraction_non_alphanumeric': 0.030718129618662614, 'fraction_numerical': 0.01704560084455545, 'mean_word_length': 4.740171445462607, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Black hole merger simulations in wave dark matter environments

The interaction of binary black hole mergers with their environments can be studied using numerical relativity simulations. These start only a short finite time before merger, at which point appropriate initial conditions must be imposed. A key task is therefore to identify the configuration that is appropriate for the binary and its environment at this stage of the evolution. In this work we study the behaviour of wave dark matter around equal mass black hole binaries, finding that there is a preferred, quasi-stationary profile that persists and grows over multiple orbits, in contrast to heavier mass dark matter where any overdensity tends to be dispersed by the binary motion. Whilst different initial configurations converge to the preferred quasi-stationary one after several orbits, unwanted transient oscillations are generated in the process, which may impact on the signal in short simulation runs. We also point out that naively superimposing the matter onto a circular binary results in artificially eccentric orbits due to the matter backreaction, which is an effect of the initial conditions and not a signature of dark matter. We discuss the further work required so that comparison of waveforms obtained with environments to vacuum cases can be done in a meaningful way.

In practise, the energy densities required to give significant effects during the inspiral and ringdown parts of the signal are in most cases high relative to the expected astrophysical values [19]. In the case of dark matter (DM), the average galactic densities as measured from observations of galactic rotation curves are at best of the order of M /pc 3 or GeV/cm 3 , with the local density in the Solar neighbourhood of the order ∼ 0.01M /pc 3 or ∼ 0.1GeV/cm 3 [35][36][37][38][39]. In the units of numerical relativity (NR) simulations, which are used to model the merger signals of compact objects, the effect of the density is measured relative to the Schwarzschild radius R s and scales with the compact object's mass. In such units average galactic DM densities are of the order Clearly some enhancement in the density around a black hole (BH) or other compact object relative to this value is required in order for the effect to be above numerical error in a simulation (which is a minimum requirement for it to be modelled and measurable in observational data). Several mechanisms that create DM overdensities around isolated BHs do exist, with one well motivated one being the formation of dark matter spikes [40][41][42][43][44][45]. These arise from the accretion and adiabatic redistribution of particle dark matter in the potential well around black holes as originally suggested by Gondolo and Silk [40]. For wave dark matter , where light bosonic particles form a condensate with astrophysical scale de Broglie wavelengths, similar accretion effects occur for which the resulting profiles have been studied in [67][68][69][70][71][72][73][74][75][76]. The density of the cloud grown via such accretion depends strongly on the asymptotic dark matter environment, but for higher mass candidates gives a power-law enhancement close to the black hole that can be significant. At the other end of the scale, where the wavelength is significantly larger than the black hole, the gradient pressure of the field (sometimes also called the quantum pressure) resists any kind of overdensity and tends to smooth out the profile, suppressing any overdensity.
Another possible enhancement mechanism is the superradiant instability, in which a bosonic field can extract energy and angular momentum from a highly spinning black hole via repeated scattering in the ergoregion (see [77] for a review). Simulations with light massive vector fields suggest these superradiant clouds can grow to be up to ∼ 10% of the mass of the black hole [78], which takes the coefficient in equation (1) to ∼ 10 −5 in the best case, and a combination of both superradiance and accretion may lead to even higher densities [62]. The potential for such bound states to form around BH binaries has been studied in [79], as well as around neutron stars [80][81][82][83].
A key question is whether overdensities that may form around isolated objects persist during a binary merger, which is what will be observed in GW data. For example, particle DM spikes have been shown with N-body simulations to disperse for equal mass mergers, meaning that objects close to merger or with a violent merger history are likely to have lost arXiv:2210.09254v1 [gr-qc] 17 Oct 2022 their DM environment [13,84,85]. For this reason, the key targets for detecting environments are extreme mass ratio inspirals (EMRIs) in LISA [86], where the cloud may still be maintained during the inspiral. Another advantage of EM-RIs is that there is a higher relative impact on the curvature for a given DM density for larger black hole masses, and the dephasing in the frequency of the signal during the inspiral accumulates over many orbits that can be potentially be observed in band or across ground based and space based detectors [87][88][89][90][91][92][93][94]. Studies of the impact of dark matter spikes [25,[27][28][29][95][96][97] and superradiant clouds [30,93,[98][99][100][101][102][103] on EMRIs show that they are potentially detectable with LISA observations.
Beyond the inspiral regime of EMRIs, it is also interesting to consider whether the strong gravitational non-linearities present during a roughly equal mass merger may give rise to distinctive features in the gravitational wave signal or other electromagnetic emissions 1 . These signatures provide information in a different regime to the inspiral and thus their combination could confirm a detection or provide evidence for a particular candidate. Studies of similar mass binaries necessitate the use of NR simulations, in which the Einstein equations are solved numerically for the evolution of the binary and its environment, from some initial state prior to merger until after coalescence and ringdown. Due to the computational expense, such simulations can at best cover the last few (order 10) orbits before the merger. Ultimately the goal is to generate waveform templates for binary mergers that include environmental effects, but a key question that ought to first be answered is whether, at such a late stage in the merger, such an environment will still be present at all, and if so, what spatial configuration it will have. In other words, what is the correct initial data for the matter environment?
In this paper we study this question for the case of wave dark matter accreting onto an equal mass, non spinning BH binary, focussing on the regime where the wavelength of the scalar is of the same order as the Schwarzschild radii of the individual black holes (which is also similar to their separation at the start of an NR simulation). We begin in Sec. III by using toy simulations of fixed BH orbits to show that in this case the accretion of dark matter onto the binary is not disrupted by the orbital motion as it is for higher mass particle candidates. There is instead a quasi-stationary profile that builds up over time, providing a well-motivated initial configuration for such matter in NR simulations.
Having identified a well-motivated profile, we then study the impact of using different initial profiles on the DM evolution in Sec. IV. We study the profile of the matter on the background of a binary merger simulated in full general relativity, initially neglecting the backreaction of the matter onto the binary motion. We compare cases that start with our 1 For example, in cases where high densities are generated during the merger, certain dark matter models may have electromagnetic counterparts arising from self-annihilations, which would provide an alternative way of identifying particular candidates, such as the Peccei-Quinn QCD axion [104][105][106][107].
quasi-stationary profile to more arbitrary configurations such as gaussians, and study the effect of cutting off the accretion at some finite radius. Finally, in Sec. V we use recently developed techniques [108] to construct constraint satisfying initial data, and turn on backreaction to study the effect of superimposing the different matter profiles on a circular vacuum inspiral. Some background information on the set up is given in Sec. II, and throughout this work we use geometric units where G = c = 1.
Our work builds on a number of earlier related investigations into DM environments of compact object binaries. The interaction of a black hole binary and a scalar field environment in the early-inspiral regime, where the separation is large and almost constant with time, has been explored via effective field theory [109,110], weak field approximations [111][112][113][114][115][116][117][118][119]; perturbative schemes [99,101,[120][121][122][123][124][125], and N-body and mesh numerical simulations [29,79,[126][127][128]. The effect of a non-vacuum environment on the post-merger "ringdown" regime, particularly the effect on the frequencies of the characteristic quasi-normal modes has also been explored by a number of authors [129][130][131][132][133]. To model the highly relativistic and dynamical merger one must use full numerical relativity. An axion-like scalar field environment was considered by Yang et al. [134], simulating the effect of a thin shell of matter on the binary merger. Choudhary et al. [135] go further, starting from the last orbit before merger and using a Gaussian as an initial profile for the scalar field. They find a change in the post-merger ringdown signal caused by the increase in effective mass of the final black hole. Ikeda et. al. [79] studied the bound states that may form around binaries via superradiance, and find they can be well described with a perturbative "gravitational molecule" description. Most recently Zhang et al. [136] simulated binary mergers with a spherical scalar field shell, and examined the effect on the gravitational recoil of the binary and the scalar and gravitational radiation. They found that the scalar cloud accelerated the merger, and increased the recoil kick. Related work has also been done in modified gravity in the context of scalar-tensor theories [137][138][139][140][141][142]. We note that in works where the scalar field grows due to superradiance or a non-minimal coupling to gravity the state of an isolated BH often depends only on the properties of the BH, and not the surrounding DM environment. However, we emphasise that similar issues to those identified in this work regarding the ambiguity of the initial state could still arise if the individual scalar clouds interact and form a common cloud prior to the point at which the simulation is started -in most (if not all) cases, no analytic form for a common binary cloud is known.

II. MODELLING WAVE DARK MATTER AROUND BLACK HOLES
For sub eV dark matter the occupation number of the particles in each state is high, with the de Broglie wavelength much larger than the particle separation [48]. We can then treat it as a classical field, in particular, a scalar field for spin-0 bosonic dark matter, which results in wave-like behaviour on astrophysical scales [56] (see [46,48,51] for reviews). Specifically, the system we consider is Einstein gravity with a minimally coupled massive complex scalar field ϕ, described by the action with a simple quadratic potential The dynamics of the scalar field is thus governed by the Klein-Gordon equation on a curved background where µ is a parameter related to the scalar field mass 2 and the field has an associated Compton wavelength λ c = 2π/µ. In this work we consider a regime where the scalar field wavelength is comparable in size to the black hole radius, corresponding to masses of ∼ 10 −9 − 10 −17 eV for BH masses between 10 − 10 9 M . The interaction of such massive scalar fields with isolated black holes has been extensively studied, showing that long-lived scalar clouds can grow around black holes either from simple gravitational accretion from the environment [47,67,68,76], or via the mechanism of superradiance for spinning black holes [77]. In the former case, the solution in the asymptotically flat region far from the black holes is a spatially homogeneous oscillatory solution of the form ϕ = ϕ 0 e −iµt , which describes a fluid of roughly constant density and zero pressure on average -i.e. dark matter. The solution closer to the black hole is described by the Heun functions [69][70][71]143], with characteristic oscillations in the spatial profile on length scales set by the scalar wavelength.
In this paper we go beyond the single black hole spacetime and consider the simplest possible black hole binary with total ADM mass M: two equal mass (M BH ≈ 0.5M) non-spinning black holes with an initial separation of d ≈ 12M on roughly circular orbits. The exact parameters are given in Table I, which result in an inspiral of about 10 orbits before merger with an initial orbital period T ∼ 270M.
We study the interaction between the scalar field and the binary black hole solving the Klein-Gordon equation with different levels of approximation in the background metric: 1.-Fixed orbit simulations -the metric background is the superposition of two isotropic BH solutions, moving on circular orbits. We evolve the scalar field on this background to test whether a dark matter profile accumulates or disperses over time. See Sec. III.
2.-Fully general relativistic evolution with G = 0 -the background is now evolved in full general relativity, but we neglect the backreaction of the matter onto the metric to focus on the impact of the binary on the matter evolution. See Sec. IV.
3.-Fully general relativistic evolution with G = 1 -finally we turn on backreaction to study the impact of the DM environment on the binary motion. See Sec. V.
In each case we take the scalar mass µ = 0.34M −1 , corresponding to a scalar wavelength of around λ c ∼ 18M, slightly larger than the black hole separation but much smaller than the binary period. For black holes in the LVK frequency band, this corresponds to a mass of ∼ 10 −9 eV, while for supermassive binaries detectable with LISA it can go down to ∼ 10 −17 eV for a binary with ADM mass 10 9 M . This choice was motivated by a brief study of the angular momentum flux as described in [144,145], where it was identified as the value that gave the largest exchange of angular momentum with the binary 3 . Our simulations therefore represent a "best case" for the impact of the dark matter on the binary.
Our results are presented in the following sections, with details of the code set up and validation contained in the Appendix.

III. FIXED ORBIT SIMULATIONS: FINDING THE QUASI-STATIONARY PROFILE
During the early stages of a binary merger, the black holes are widely separated and follow approximately Keplerian orbits, with emission of gravitational waves tending to circularise the orbits over time [146][147][148][149]. During this phase, the dark matter profiles of the two BHs will evolve largely independently, generating nontrivial density profiles via accretion. A common DM cloud will also tend to circularise the orbits via the effect of dynamical friction [150]. As the binary separation decreases, at some point the DM clouds will merge and interact 4 , eventually forming a common cloud which is "stirred up" by the binary. It will not, therefore, have a smooth gaussian profile and will carry some angular momentum due to its interaction with the spacetime curvature. To investigate the resulting configuration, we construct an approximate toy model for the late inspiral where we model the gravitational field as a superposition of two isotropic black hole metrics where dΩ 2 = dθ 2 + sin 2 θ dφ 2 and FIG. 1. Evolution of the real part of the scalar field (top) and energy density (bottom) on the fixed orbit binary background after 0, 2, 4 and 6 orbits respectively. The binary generates a scalar cloud that quickly settles into a quasi-stationary spiralling profile that grows in amplitude over time. In the centre around each black hole spikes form in the scalar field amplitude, resulting in an enhanced energy density, with an additional accumulation of matter in the potential well in the middle of the binary. Further out we see regions of both higher and lower density forming, with the patterns on a length scale related to the binary separation and scalar wavelength. Movie can be found in https://youtu.be/XevfJKLO9ec.
is an effective gravitational potential. We impose that the black holes (m 1 = m 2 = M BH ) follow circular Keplerian orbits of radius d/2 and frequency with their centres located at Close to each of the black holes this metric tends to a Schwarzschild metric in isotropic coordinates, while far from the black holes it tends towards a weak field limit. This metric is not a solution of the Einstein equations, but is merely designed to study the way in which a common quasi-stationary profile can form in a period where the orbits are not yet rapidly decaying as at merger. We start the simulations with the binary immersed in a homogeneous scalar field, choosing ϕ(t = 0) = ϕ 0 and Π(t = 0) = −iµϕ 0 , and study the evolution of the scalar field over several orbits. The (real) parameter ϕ 0 controls the asymptotic density of the infinite reservoir, but since we neglect backreaction for this simulation, we can rescale the results to any physical asymptotic density we choose.
We find that the scalar field rapidly accretes from its asymptotic value into a cloud around the two black holes and forms a persistent spiral profile within a few orbits that co-rotates with the binary, see Fig. 1. The scalar field (top panel) is pushed towards large amplitudes, such that the energy density (bottom panel) around and between the black holes increases by several orders of magnitude. Once the profile forms it grows over time homogeneously, fed by the asymptotic reservoir of dark matter imposed at the boundaries.
This quasi-stationary profile can be studied in more detail in Fig. 2, where we plot the evolution of the density and scalar field profiles along the axis of the binary, normalised relative to their central values. Even within the first two orbits, both quantities have already settled into the persistent profile that grows steadily in amplitude over time. We observe density spikes around each of the black holes, an accumulation of scalar matter in the potential well between the black holes, and smaller amplitude density peaks further away on lengthscales that depend on the binary separation and scalar wavelength. The scalar cloud is significantly enhanced compared to a superposition of two of the scalar field profiles found around isolated black holes (like those studied in [64,[67][68][69][70][71]), due to the non-linear effect of the combined gravitational potential.
We note that this persistent profile contrasts with higher mass particle-like dark matter, where dark matter density spikes have been shown to disperse under the influence of the binary motion [13,84,85]. The enhanced density is important both for its effect on the gravitational wave signal and also for FIG. 2. Normalised density and scalar field profiles along half of the fixed orbit binary axis. The black dot illustrates the location of one of the black holes. We see that after a relatively short amount of time (approximately two orbits) the field has gone from a homogeneous profile to a quasi stationary profile peaked around and between the BHs (the zero values within the horizon are due to the choice of gauge and excision conditions for the evolution). This results in a region of enhanced density in the potential well between the binary BHs. potential direct detection for models with standard model couplings [105].
Even with this simple fixed orbit model, the finite size of the numerical domain prevents us from evolving the system for many orbits. However, the persistent and consistently growing scalar field profile shown here allows us to predict the qualitative behaviour of the cloud at late times -it seems that a common spiral shaped cloud would continue to grow until the dark matter reservoir is eventually exhausted, at which point the binary would be left with an isolated spiral shaped cloud that would gradually decay away. Assuming that the reservoir is not exhausted before merger, the cloud should continue to grow and not disperse. 5

IV. GR EVOLUTION WITH G = 0: THE IMPACT OF THE BINARY ON THE MATTER EVOLUTION
We now study the evolution of the wave dark matter on a fully general relativistic binary merger, with the Klein-Gordon equation evolved on a dynamical spacetime described by the Einstein's field equations where R µν is the Ricci tensor, and T µν is the energy momentum tensor of the scalar field. To first isolate the impact of the binary on the matter evolution, we turn off the backreaction by setting Newton's constant G = 0 in Eqn. (10). This ensures that different cloud configurations evolve in the same BBH background, described by the black holes' vacuum trajectories.
The BH initial conditions are Bowen-York data [152][153][154][155][156][157][158][159][160] for the parameters in table I and we solve the Hamiltonian constraint for the correction to the conformal factor arising from the non zero boosts.
The main goal is to compare the evolution of several initial scalar field configurations. First, we observe that the profile found in our toy model with fixed orbits is close to the stationary profile for the fully GR binary solution, as we see very little transient evolution when we impose it on the binary and begin the evolution. This can be contrasted with other possible choices for the initial scalar field profile, such as a Gaussian, where we find significant initial transients before the preferred profile is reached.
The scalar profiles we examine can be classified into two types: extended and isolated clouds. In the former class, the energy density of the cloud reaches the boundary of our simulated domain, allowing for continued accretion from spatial infinity. (This is imposed using extrapolating boundary conditions as described in [161].) Within this class we study: • Fixed Orbit (FO): We take as initial data the resulting quasi-stationary scalar field profile after evolving the binary in the homogeneous dark matter halo for 6 orbits ϕ FO , see Fig. 1. This can be considered the "most correct" initial condition for a scalar field that is still accreting up to the merger.

• Homogeneous (H):
We take the same homogeneous initial profile described in the previous section by setting ϕ(t = 0) = ϕ 0 and Π(t = 0) = −iµϕ 0 . In this profile we still need to grow the quasi-stationary profile, so expect some transient evolution.
In the second class of initial conditions, the scalar cloud has a sufficiently large radius to cover the BBH, but the energy density goes to zero at the boundaries, reproducing an isolated cloud which has exhausted its dark matter reservoir. (For these cases we use Sommerfeld radiative boundary conditions as described in [161].) Here we study two cases: • Fixed Orbit Isolated (FOI): We apply a f (r) ∼ tanh[r − 60M] envelope to the FO scalar field profile, so that ϕ FOI (t = 0) = f (r)ϕ FO (t = 0) and the configuration reproduces an isolated cloud which has exhausted its dark matter reservoir. Again this is a "correct" profile around the black holes, but the cut off we introduce is rather arbitrary, and therefore introduces some transient evolution in the overall shape of the cloud as it settles into a quasi-stationary, isolated profile. • Gaussian (G): We choose a Gaussian profile for ϕ(t = 0) = ϕ 0 exp −r 2 /σ 2 and Π(t = 0) = −iµϕ(t = 0), where r is the distance to the centre of the binary. We use σ = 1/0.03M and choose ϕ 0 so that the total mass of the initial scalar cloud is the same as the isolated fixed orbit cloud described above. We expect transient evolution in this case before we arrive at the quasi-stationary profile.
We plot 2D slices of the initial density configurations perpendicular to the binary orbital axis in Fig. 3.
To quantify the differences we extract the value of the scalar field profile along the line joining the BHs, see the top panel of Fig. 4. In the bottom panel of Fig. 4 we track the value of the scalar field at the centre of the binary, which provides a reasonable indication of the amount of transient evolution. We see as expected that the scalar field needs some time to settle down to the quasi-stationary configuration, with the most transient evolution in the case G and the least in FO. 6 However, 6 Note that the scalar field in the FO and FOI initial data goes to zero inside the horizon. However when we evolve it in the full GR evolution in the moving puncture gauge [162,163] it quickly relaxes to give a continuous non-zero density spike. This difference is mainly due to the different choice of lapse between the fixed orbit metric and the moving punctures gauge used for the full relativistic evolution. One should thus consider the change in the sub-horizon initial scalar field profile to be an artefact of the change in gauge choice. even when starting from very different profiles the scalar field evolves into the same kind of spiral shaped, quasi-stationary, configuration as described in section III within a few orbits. The final state post merger is a single density spike with power law tails, of the kind studied in [64,[67][68][69][70][71]76], with solutions that can be approximated by confluent Heun functions [143].
The main difference between the extended (FO, H) and isolated (FOI, G) clouds is that the former can continue to accrete from infinity, so we see the central amplitude continuing to grow throughout the merger, increasing the density and size of the final scalar cloud. The homogeneous data accretes rapidly at the centre, and in only a few orbits reaches the fixed orbit quasi-stationary configuration. The difference in amplitude between the final cloud from the homogeneous data and FO data is expected given that the fixed orbit data has effectively been accreting for an extra six orbital periods during our fixed-orbit simulation.
For the isolated configurations (FOI and G), the cloud does not continue to grow as there is no reservoir of DM at the boundary. In the case of FOI, because there is no asymptotic density to confine it, the cloud settles into a more distributed (but qualitatively similar) shape -this is a consequence of imposing an arbitrary cut off in the density at finite radius. The Gaussian case, on the other hand, shows strong transient behaviour with several large oscillations around the quasistationary profile, on a timescale roughly corresponding to the period of the binary. After ∼ 2-3 orbits the same configuration is reached -the fact that the final shapes and amplitudes match is a result of starting with the same cloud mass, and implies that the two cases have radiated equal amounts, and been FIG. 4. We plot the |ϕ| profile along the axis of the binary for different initial data profiles. The amplitudes of the fixed orbit isolated (FOI) and Gaussian profiles are normalised so that the integrated mass of the whole cloud is the same. The amplitude of the homogeneous profile is chosen to match the initial asymptotic value of the fixed orbit profile and is set to unity in these plots. The black dot illustrates the location of one of the black holes. The bottom panel plots the evolution of the scalar field amplitude at the centre of the binary with time. Markers illustrate the time of the upper panels. A movie showing the evolution of these profiles can be found in https://youtu.be/wGUUUutuiyU. accreted equally by the BHs, during their transient evolution.
These results demonstrate that within only a few orbits the transients in the cloud will die away and a quasi-stationary configuration will be reached. However, the initial profile can make a significant difference to the transient evolution of the scalar cloud during this time, and therefore may have an impact for shorter simulations. We have not yet considered the backreaction of matter onto the metric, but where this is included it is possible that these transient effects could impact on the black hole trajectories, as well as creating additional radiation of gravitational waves. We will discuss this further in the following section.

V. GR EVOLUTION WITH G = 1: THE IMPACT OF THE MATTER ON THE BINARY EVOLUTION
The ultimate goal of NR simulations with environments is to quantify the impact of the matter on the binary evolution, and resulting gravitational wave emission, so as to learn about the properties of the DM. To do so we must include the backreaction of the evolving matter on the metric background. We therefore restore G = 1 and solve the Hamiltonian and momentum constraints using the novel CTTK method [108]. In particular, we use the hybrid CTTK approach, where we choose a spatially varying mean curvature that depends on the energy density distribution K 2 = 24πGρ, where ρ is the ADM energy density, and solve the constraints obtaining corrections for both the conformal factor and traceless conformal extrinsic curvature, on top of an initial Bowen-York solution for a boosted black hole binary in vacuum. This means the matter environment introduces corrections to the effective initial momenta and bare masses of the black holes. Now that the amplitude of our field has a physical meaning (it is related to the matter density roughly as ρ ∼ ϕ 2 0 ), we quantify this by specifying the relative cloud mass to (vacuum ADM) binary mass ratio.
We superpose our profiles onto the vacuum BH parameters for the masses and momenta, and evolve the metric plus matter to merger, which we find happens faster in the presence of the DM clouds. Some dephasing of the signal compared to the vacuum evolution is expected due to the effects of dynamical friction, radiation of the scalar cloud, and backreaction onto the metric, but the impact we see is unexpectedly large. Examining the black hole trajectories in Fig. 5 reveals that, even for light clouds where M cloud /M ≈ 0.01%, the black holes are pushed into eccentric orbits, giving rise to large deviations from the circular trajectories that we would expect at this late stage of the inspiral and a prompt merger. We find that this effect is proportional to the energy density near the black holes, being the largest for both the extended and isolated fixed orbit profiles (FO and FOI) due to the presence of density spikes near the horizons. For the more artificial initial conditions H and G where the initial energy density near the black holes is smaller, the effect is smaller but still non-negligible. The use of the CTTK method means that in general we obtain larger corrections to the effective initial black hole momenta, and smaller corrections to the effective bare masses, compared to the CTT method. However the unwanted eccentricity is an unavoidable consequence of the non-trivial matter environment, and is particularly significant where there are high energy and momentum flux densities near the black hole horizons, as in our preferred quasi-stationary solution.
Because the deviations are significant, we cannot simply extract waveforms and compare them to the vacuum case to determine how the presence of a scalar cloud impacts the gravitational wave signal, as we would not be able to conclude that such signatures are caused by the presence of the scalar field rather than by the induced (physical but unwanted) eccentricity. Even if one tunes the initial momenta and masses of the black hole punctures so that their initial coordinate velocities and accelerations match the vacuum case, this is not sufficient to remove the effect. Further work is required to remove this effect. One should follow a similar method to those used to obtain low-eccentricity initial data for binary black holes in vacuum spacetimes [164][165][166][167]. In these methods, one measures the trajectory and eccentricity over 2-3 full orbits for an initial choice of input parameters, then uses a high order post-Newtonian [168] approximation or gradient descent to estimate the correction to these input parameters (typically the initial momenta and masses of the black hole punctures) needed to reduce the eccentricity. One then iterates this procedure until the desired eccentricity is achieved. Our case is more complicated because the black holes are perturbed by the backreaction from the matter distribution, which we solve for numerically. However, one could in principle obtain a similar scheme by incorporating the effect of a matter cloud into a new post-Newtonian binary model, something which is beyond the scope of this paper 7 .

VI. DISCUSSION
One of the key challenges in constructing reliable gravitational wave templates for BH binary mergers with environments is establishing the correct initial data for numerical relativity simulations, which begin only a relatively short time before the merger.
Here we conducted numerical simulations of the accretion of wave dark matter around binary black holes in fixed orbits, choosing a scalar mass that gives a large interaction with the BH binary (i.e., the regime where the wavelength of the field is similar to the Schwarzschild radii of the BHs). We found that the scalar field quickly converges to a persistent non-trivial profile, peaked around each black hole, which grows in amplitude as the accretion continues.
We then explored how different choices of initial data affected the evolution of the DM cloud during the binary merger. Our results suggest that the profile found in the fixed orbit simulations is an attractor solution, as other choices of initial data converged to the same distribution over the course of several orbits. This means that for long numerical simulations the precise choice of initial scalar profile may be unimportant, as the field will quickly converge to the quasi-stationary distribution. However for numerical simulations which only consist of a small number of orbits, the transients and the resulting loss of control over the initial data from using a non stationary profile may affect the final results. We also saw that allowing continued accretion from a scalar dark matter reservoir over the course of the merger results in a larger scalar cloud around the final black hole, confirming the importance of the local DM environment for the system's evolution.
Lastly, we obtained constraint satisfying initial data for the scalar profiles obtained from our fixed orbit simulations. To the best of our knowledge, this is the first time that environments with non zero angular momentum have been studied in the initial conditions for DM environments around binaries. We found that naively using the parameters for a quasicircular binary merger in vacuum will not usually produce low eccentricity orbits once matter is added, and because of this one cannot simply compare the resulting waveforms to establish the signatures of dark matter environments. For example, one may find that the merger happens sooner merely because of the increased eccentricity, and not due to effects like dynamical friction and gravitational radiation.
We conclude that a more comprehensive approach is needed to find initial data that both satisfies the Einstein constraints and produces realistic low-eccentricity inspirals, similar to that developed for vacuum binaries [156,157,[164][165][166]. This should ideally incorporate the quasi-stationary profile as part of the process of solving for the initial data, as is done in neutron star simulations [175][176][177], or match to a post-Newtonian model [168][169][170][171][172][173] for the matter and black holes in the late inspiral regime.
Only once the problem of initial conditions is under control can the degeneracies between the effects of a non-trivial dark matter environment and changes in other parameters of the binary be quantified, so that we can examine how to disentangle the two. We may then also extend the parameter space to include unequal mass and spinning black hole binaries, where kicks may occur [26], as well as exploring a wider range of scalar field masses, and the impact of self-interactions.

ACKNOWLEDGEMENTS
We would like to thank Lam Hui and Macarena Lagos for their input and discussions in the early stages of this project. We also thank Giuseppe Ficarra We use the 3 + 1 formalism [152,160,[179][180][181] to foliate the four dimensional spacetime metric into three-dimensional spacelike hypersurfaces indexed by a time coordinate t with line element ds 2 = −α 2 dt 2 + γ i j (dx i + β i dt)(dx j + β j dt), where α is called the lapse, β i the shift and γ i j the spacial metric. The normal to the hypersurface is given by n µ = (−1/α, β i /α). In addition, we decompose the complex scalar field into two real scalar fields ϕ a = (ϕ 1 , ϕ 2 ) as ϕ = ϕ 1 + iϕ 2 , and reduce the second order Klein-Gordon equation to two first order differential equations where Π a is the conjugate momentum of ϕ a , K is the trace of the extrinsic curvature K i j = 1 2α (−∂ t γ i j + D i β j + D j β i ), and D i the covariant derivative associated to the spatial metric γ i j .
For the simple fixed orbit model described in section III we only solve equations (12) and (13) in such a background. We excise (set the evolution variables to zero) within a small region around each black hole centre, inside the horizon, to avoid numerical errors, and evaluate the values of the metric components and their derivatives analytically at each point on the grid. For the relativistic merger, in addition to solving Eqns. (12) and (13), we solve the full Einstein equations numerically. For this we use the CCZ4 formalism [182] with the moving puncture gauge [163,[183][184][185][186]. In both cases we use the open-source numerical relativity code GRCHOMBO [187,188] with adaptive mesh refinement [161].
We use a simulation box length L = 512M and 8 levels of mesh refinement (See Figs. 6 and 7 for convergence tests). Taking advantage of the symmetry in the xy plane we impose reflecting boundary conditions at z = 0, while for the other boundaries we impose either first order extrapolating boundary conditions (matching the first derivative on the exterior ghost cells to that inside the simulation grid) or Sommerfeld boundary conditions. As a diagnostic quantity, we define the effective energy density ρ E measured by timelike observers as with ADM quantities ρ = −n µ n µ T µν and S i = −n a T ia , such that the total energy in a 3D volume is E = ρ E dV .