Relativistic drag forces on black holes from scalar dark matter clouds of all sizes

We use numerical simulations of scalar field dark matter evolving on a moving black hole background to confirm the regime of validity of (semi-)analytic expressions derived from first principles for both dynamical friction and momentum accretion in the relativistic regime. We cover both small and large clouds (relative to the de Broglie wavelength of the scalars), and light and heavy particle masses (relative to the BH size). In the case of a small dark matter cloud, the effect of accretion is a non-negligible contribution to the total force on the black hole, even for small scalar masses. We confirm that this momentum accretion transitions between two regimes (wave- and particle-like) and we identify the mass of the scalar at which the transition between regimes occurs.

Gravitational interactions with compact objects are one of the most promising tools for investigating DM properties, since they do not rely on any additional interactions with the Standard Model.Extreme mass-ratio inspirals (EMRIs), in particular, provide an optimal system for studying environmental effects on the inspiral gravitational waveform, since they may complete about 10 4 − 10 5 orbits before merger, meaning that small dephasing effects are integrated over long timescales.Since these kinds of systems typically reside in the center of a galaxy, the smaller object is expected to pass through the DM core, where densities are highest [38][39][40].However, even with next-generation detectors, prospects of observing a signal often rely on enhancements in the density above those in the core, e.g., due to superradiance, accretion of DM spikes, or self-interactions in the DM (see, e.g., [41][42][43][44][45][46][47][48]).
Several effects are expected to give rise to dephasing in a binary's gravitational-wave signal when additional matter is present.A key one is dynamical friction (DF): a gravitational drag force due to an overdensity (a "gravitational wake") that develops behind a massive object as it moves through a medium.First described by Chandrasekhar [49] for a nonrelativistic Newtonian perturber moving through a cloud of noninteracting particles, it was then extended to different media (such as fluids [50][51][52]), different geometries (e.g., spherical or slab-like [53][54][55]), and including relativistic corrections [56][57][58][59].In the context of scalar field DM, the DF was first computed for a nonrelativistic Newtonian perturber [9], and then extended to include velocity dispersion [60,61], self-gravity [38,39] (see also [62,63]) or self-interactions [41,64], for binary systems [65,66], and relativistic perturbers [67,68].It has recently been shown that relativistic effects can play an important role during the final few orbits of an EMRI, producing a detectable effect on the evolution of the binary [69].
Another effect responsible for dephasing in the inspiral of BH binaries is the accretion of matter (and its momentum) onto the BHs as they move relatively to the medium, as originally studied by Bondi [70,71].In most cases it results in a smaller effect than DF (see, e.g., [72]), but it is nevertheless important to characterize it, especially in the case of small environments: while the strength of DF increases with the medium (or, more precisely, the wake) size, the effect of momentum accretion is roughly independent of it, and so becomes more important for smaller environments (wakes).In the context of light scalars this accretion was studied by Unruh [73] (see also [41,42,68,74,75]).
In a previous work [67] some of us characterized the DF effects for large scalar DM clouds, by evolving numerically the scalar field around a BH in uniform linear motion, and (based on fluid-media results [57,58]) suggested a phenomenological model for the relativistic corrections which included a pressure-like and Bondi accretion terms.Soon after, analytic expressions for the drag force on the BH were obtained for a similar setup [68], suggesting that our phenomenological terms could be removed for very light scalars (i.e., γM mc/ℏ ≪ 1, where γ := 1/ √ 1 − v 2 is the Lorentz factor of a BH of mass M and velocity v) by using the analytic expression for Unruh accretion [73], without the need of any free parameter.
In this work, we evolve numerically the scalar field on a moving BH spacetime to confirm that, indeed, the phenomenological relativistic corrections introduced in our previous work can be accounted for by a better modelling of the accretion process.In particular, we find that Unruh accretion captures very well our numerical results for sufficiently light scalars; this is particularly evident for small clouds, which we first simulate in this work, where the accretion of momentum gives a comparable (or dominant) effect to DF.We also derive the drag force on the BH in the geometrical optics limit (for a cloud of particles following timelike geodesics), which is independent of the particles' spin and, in particular, applies to CDM.The (semi-)analytic expressions that we provide are shown to describe well our numerical results in the different regimes.
Hereafter, we adopt the mostly positive metric signature and use geometrized units in which G = c = 1.The scalar field mass will be parametrized by the inverse (Compton) length-scale µ := m/ℏ, and we will often use the dimensionless quantity α s = M µ to present our results.For reference: a value of α s ∼ 10 −2 corresponds to a scalar with m ∼ 10 −12 eV for a solar-mass BH (M ∼ M ⊙ ), or to m ∼ 10 −22 eV for a supermassive BH (M ∼ 10 10 M ⊙ ).

II. THEORY
We consider a complex scalar φ minimally coupled to gravity described by the action which results in the Klein-Gordon equation We assume a scalar field dilute enough that its backreaction on the spacetime geometry is negligible at leading order; so, effectively we consider that the scalar field evolves on a (vacuum) Schwarzschild geometry.The scalar field energy-momentum tensor is where

III. COORDINATE SYSTEMS
We start with the Schwarzschild metric in isotropic coordinates ( t, x, ȳ, z), corresponding to the "BH frame", which we then boost by a factor γ in the ∂/∂ x direction; the resulting coordinates (t ′ , x ′ , y ′ , z ′ ) correspond to the "scalar field frame".By adding a spatially-constant shift in the x ′ -coordinate (i.e., x := x ′ − vt ′ ), we obtain a time-invariant metric in which the BH remains at a fixed coordinate position -we call this coordinate system, in which we perform the numerical evolution, the "simulation coordinates" (t, x, y, z).
The 3 + 1 ADM decomposition of the Schwarzschild metric in simulation coordinates is where the lapse, shift and nonzero components of the spatial metric are, respectively, where with r2 := γ 2 x 2 + y 2 + z 2 .While we perform the numerical computations in simulation coordinates, we will present the final gravitational drag forces in the BH frame, where the analytic expressions are more naturally derived.The rate of change of the ADM momentum in the BH frame can be obtained from the one in simulation coordinates using In Ref. [67] we assumed dP ADM t ≈ 0 for small dP ADM x , but did not justify this further.Whilst, as suggested in Ref. [68], this assumption is not necessarily valid, the treatment in Ref. [67] was self-consistent.We provide a clearer and more thorough justification of the assumptions relating to this accretion term in Appendix A.

IV. NUMERICAL FRAMEWORK
Our numerical setup is substantially the same as in Ref. [67].For completeness, we briefly recap the main elements of the methods in this section, and refer the reader to [67] for more information.The technical details on the numerical grid setup are given in Appendix D.
We evolve the scalar field by solving the system on a fixed Schwarzschild geometry, where Π is the conjugate momentum, as defined by Eq. ( 8), and K is the trace of the extrinsic curvature of the background K ij , which in simulation coordinates is simply given by since the metric is time-invariant.We set homogeneous initial conditions for the scalar field across the grid, with Re Π(t = 0) = 0, Re φ(t = 0) = φ 0 , Im Π(t = 0) = µφ 0 , and Im φ(t = 0) = 0. We use an initial amplitude φ 0 = 0.1, but this is an arbitrary choice, since we neglect the backreaction of the field onto the metric and the system is linear, which implies that the final result can be rescaled to different physical densities (assuming that its backreaction remains negligible).

V. GRAVITATIONAL DRAG
To compute the relativistic drag force acting on the moving BH we use the framework developed in Ref. [76], which allows us to find the leading order term in the scalar rest mass density using the test field approximation.The drag force includes both the effects of DF and momentum accretion, and is defined as where t ν µ is the Einstein's pseudotensor of the total spacetime metric g µν , which includes the backreaction from the scalar field (see, e.g., Ref. [76]).The "curvature momentum" P g i (and the force F i ) depend on the slicing of the spacetime Σ o (t) -they are well defined once the observers are specified.
For an asymptotically flat spacetime, the ADM momentum can be decomposed into a curvature part and a scalar field part, where, at leading order in ϵ (with |φ| ∼ ϵ), Here dS j := d 2 x √ σN j , where σ is the determinant of the induced metric on ∂Σ o and N j its outward unit-normal.Thus, at leading order in ϵ, the drag force can be obtained directly from an evolution of the scalar field on a fixed background -sidestepping the actual computation of the backreaction on the metric.That is, Finally, as shown in Appendix A, the last expression implies that the steady-state drag force on the BH frame (i.e., in isotropic coordinates) can be computed from where the right-hand side is to be evaluated numerically in simulation coordinates, with Σ i ⊂ Σ o a 3-dimensional surface outside the horizon, which contains the curvature singularity.

VI. ANALYTIC EXPRESSIONS
The DF acting on a point-like Newtonian perturber moving at nonrelativistic velocities through a scalar field cloud was first derived in Ref. [9].These expressions were extended to the case of a BH moving at relativistic speeds in Ref. [68], including also the drag force from accretion of momentum.We summarize here the key results, which we will validate against our simulations.
There are two important dimensionless parameters in this problem: the ratio of the BH size to the reduced (relativistic) Compton wavelength of the scalars and the ratio of the characteristic scattering radius to the reduced de Broglie wavelength These parameters control the wave effects, respectively, in the accretion and scattering processes [68], with the field behaving as particles in the semi-classical limits ᾱs ≫ 1 or β ≫ 1, for large azimuthal numbers ℓ ≫ 1 [77]. 1nother relevant dimensionless parameter is the ratio which characterizes the radius r of the cloud (or, more precisely, of the wake) in units of the de Broglie wavelength λ dB .The drag force due to accretion of momentum is independent of r and becomes increasingly important (as compared to DF) for Λ ≲ 1.Although we can find (semi-)analytic expressions for the steady-state drag force in all regimes, they only take a simple closed form in particular limiting cases; all the expressions in this section are given in the BH frame (isotropic coordinates) and have the form (20) with ρ the asymptotic rest-mass density of the medium, and the coefficients D and A characterizing the contribution from DF and accretion, respectively.Let us start with very light scalars (ᾱ s ≪ 1).
where Ψ is the digamma function.For smaller scalar field clouds, but still much larger than the BH size (r/M ≫ 1), we can also find a simple closed-form expression in the wave limit β ≪ 1: e πβ πβ sinh(πβ) A Unruh , (22) where Cin(z) := z 0 (1 − cos t)dt/t is the cosine integral.In both expressions the terms in the first line are an extension to relativistic velocities of the DF expression derived in Ref. [9], while the term in the second line is due to (Unruh) accretion of momentum. 2.
For heavier scalars ᾱs ≫ 1, the field behaves as collisionless particles and one recovers geodesic results (cf.Appendix B); using the WKB approximation one can show that, for large clouds Λ ≫ β, the drag force is 2 The momentum accretion was actually derived in Ref. [68], but it can be obtained directly from the expression for the accretion rate derived by Unruh in [73] 0.012 0.018 0.024 21) D(ᾱ s 1, Λ β) + A particle , Eq. ( 23) overestimates them for the other cases where αs ≲ 1 (middle and bottom panels).Interestingly, the numerical results for these masses seem to be well described by the combination D( ᾱs ≪ 1, Λ ≫ 1) + A( ᾱs ≫ 1) (shown as black dashed curves), that is, the accretion is in the particle limit, whilst the dynamical friction is in the wave limit.
where the critical impact parameter, show the analytic expression for the total force in Eq. ( 22).We see that the wave-like expressions for the accretion and dynamical friction are an excellent fit.
and χ(v) is a general-relativistic correction to DF given in Appendix B. As expected, the contribution of DF to the total drag force is the same as in Eq. ( 21), in the semi-classical limit β ≫ 1, modulo general-relativistic corrections contained in χ(v), which originate from particles with impact parameters b ≳ b cr .For small clouds we cannot find an analytic expression for the drag force, since the integral (B2) for the DF contribution must be evaluated numerically.

VII. NUMERICAL RESULTS
We focus our simulations on the boundaries between different regimes of cloud sizes Λ and scalar masses α s , validating the analytic expressions derived in Ref. [68] for light scalars (ᾱ s ≪ 1), and in Appendix B for heavy scalars (ᾱ s ≫ 1).These expressions contain no free parameters.Our numerical results for the drag force are obtained extracting all necessary quantities from our simulations to evaluate Eq. ( 16).The masses we consider are in the range α s ∈ [0.025, 1] and, so, our simulations could not really probe ᾱs ≫ 1.The numerical limitation to probe this limit has to do with the increasingly large frequency of the scalar field oscillations, which becomes increasingly harder to resolve. Figure 1 shows our numerical results along with the plots of the analytic expressions discussed in the previous section for large clouds with different scalar masses (µr, α s ) ∈ {(45, 0.05), (300, 0.5), (300, 1)}.We confirm that for the lightest scalar considered, α s = 0.05 (top panel), the total drag force on the BH is correctly accounted for by the expression in Eq. (21).We can see from the middle and bottom panels that, at larger masses, the (wave-like) Unruh accretion significantly overestimates the total force.In these cases, a better approximation for the drag force from accretion is given by A particle , the second term in Eq. ( 23), although formally it was derived only in the ᾱs ≫ 1 limit.Adding the DF from particle DM (general-relativistic) deflection slightly overestimates the force (as can be seen from the red dotted lines in the plots) for the heavier scalars considered α s = 0.5 , 1, but we expect it to become more accurate for heavier scalars (α s ≫ 1).For intermediate scalar masses (α s ∼ 1), the best description of the drag force on the BH seems to be provided by a combination of the DF expression derived for ᾱs ≪ 1 with the (collisionless particle) accretion expression derived for ᾱs ≫ 1 (black dashed lines).Note that the Unruh accretion is fundamentally wavelike and it only includes the ℓ = 0 mode, which explains why it becomes an increasingly worse description of the numerical results as we approach ᾱs ∼ 1; this is mainly due to the no inclusion of higher ℓ modes, which become effective precisely at ᾱs ∼ 1.On the other hand, the DF expression derived for ᾱs ≪ 1 includes higher ℓ modes and, so, also captures the particle regime.In fact, the only difference between D(ᾱ s ≪ 1, β ≫ 1) and D(ᾱ s ≫ 1) are general-relativistic corrections to DM particle deflection, which we speculate to become effective only at ᾱs ≫ 1.Our results show that that the transition from wave to particle-like regime occurs for scalar masses α s ∼ O(0.1).However, for the DF this transition appears to be much more gradual than for accretion, and for α s ∼ 1 some (general-relativistic) particle effects are still not active.
Figure 2 shows the comparison of our numerical results with the analytic expression (22) for small clouds Λ ≲ 1 in the wave regime β ≪ 1.We confirm that this analytic expression provides an excellent description of our simulations for small clouds (µr = 2.5) with scalar masses α s = 0.05 and 0.025.For these small clouds the contribution from accretion to the drag force is comparable or more important than DF, and it is then crucial to model it well.While our analytic expressions (based on Unruh accretion) fit well the numerical results for Λ ∼ 1, the phenomenological ones used in Ref. [67] (based on Bondi accretion) are not good enough (see Appendix C).

VIII. DISCUSSION
Future detections of EMRIs by space-based detectors like LISA, TianQin and Taiji [78][79][80][81] will open up new windows on BH environments.It has been suggested that superradiant clouds may be detectable by LISA [82-  20)-( 23).In the case of a scalar field wavelength comparable to the size of the BH and smaller and very small scalar field clouds (shaded region on the bottom left of this graph), there is no valid analytic expression for the DF force.Here we can only calculate the drag force on the BH numerically, but as one can see from our results, the analytic expressions give a good order of magnitude estimate even in this regime.
In this work we have calculated numerically the force resulting from DF and accretion on a BH boosted through a uniformly dense scalar field medium, and validated the analytic expressions derived from first principles in their different regimes of validity.In particular, we have shown that the total relativistic drag force has the form given in Eq. (20), where the coefficients D and A contain, respectively, the contributions from DF and momentum accretion.The different regimes for which analytic expression are known are summarized in Fig. 3.As discussed in Ref. [68] and confirmed here, the additional "pressure" correction that we had considered in our previous work [67] is not necessary when the accretion is correctly accounted for (by using Unruh's expression instead of Bondi's), which becomes particularly evident for the smaller clouds (Λ ≲ 1) we have considered here.Interestingly, we found that at intermediate masses ᾱs ≲ 1 the accretion process transitions between a wave description to a particle one (though the latter expression is formally derived in the ᾱs ≫ 1 limit).The general-relativistic effects on particle DF are expected to become effective only at ᾱs ≫ 1, larger than the masses considered in our simulations.
Having established the methods for extracting and quantifying the relativistic drag forces on BHs moving though scalar field DM clouds, our simulations can now be extended to more complex cases, such as those including DM self-interactions or BH spin, and other fundamental fields such as massive bosons of spin 1 and 2. Another interesting extension of our work is to go beyond the test field approximation.All our results were obtained under the assumption that the diluted cloud has a mass much smaller than one of the BH.If that assumption holds, the subleading corrections can be captured by a similar approach to the one used here (but more complicated).Systems for which that assumption does not hold need a different (semi)analytic approach (e.g., [39,98]).momentum is negligible compared to the one in the spatial part, and so we ignored the contribution of dP ADM t /dt in Eq. (7).While this assumption is not formally valid, we provide here a more careful treatment which confirms that the results presented in Ref. [67] are self-consistent, and thus correct.
Consider three concentric 2-dimensional spheroidal surfaces around the BH: the outer surface ∂Σ o , where the ADM mass and overall change in momentum is defined, the inner surface ∂Σ i , which is the last well-resolved surface in our simulations before the BH event horizon, and the BH horizon itself ∂Σ BH .The rate of change in the scalar field momentum in the region between ∂Σ BH and ∂Σ o is where dV := d 3 x √ −g/α.The flux across the BH horizon in our coordinates is always zero, so we can discard the second term on the right-hand side.In the steady state the cloud profile is approximately stationary, so the left-hand term should be small (but not negligible). 3 In the simulation coordinates, there is always a (small) steady state growth in the momentum of the cloud of Ṁ v, where Ṁ is the steady-state increase in the BH mass over time measured in its rest frame with respect to proper time in that frame.Whilst in that frame this causes an increase in the BH mass and momentum over time, in our coordinates it manifests as the field collecting outside the horizon.Therefore we have that We also need to note that the force acting on the BH is not simply dP ADM i /dt ≈ − ∂Σo dS j αT j i , unless dP φ i /dt = 0. 4 We need to isolate dP g i /dt (see Eq. ( 11)), using the 3 In practice it is large initially, but we want to exclude these changes, which are related to the transient growth in Σo−Σ BH dV αT 0 i that occurs in our simulations while dynamically evolving towards the steady state.The cloud profile rapidly settles into the steady-state one at smaller r, and more gradually spreads outwards.Whilst there should also be a change in the source term once the outer profile settles into the steady state, we have seen in previous work that this volume integral is dominated by the part at smaller r, and so in practice it also settles quickly to its steady state value [76]. 4The rate of change of the ADM momentum gives the total force acting on the system (BH + scalar) inside the volume of interest Σo.For a cloud with compact support (and using a sufficiently large ∂Σo), we would have dP ADM i /dt = 0, meaning that the total force acting on an isolated system vanishes, but the force acting on the BH may still be nonzero if it exchanges momentum with the scalar.
fact that in the steady regime in the simulation coordinates.Thus, in these coordinates the force on the BH is For the energy fluxes, the analysis is the same.However, since the background spacetime has a time-like Killing vector, the (volume) source term is exactly zero.Thus, in simulation coordinates, in the steady state we have that where the Ṁ /γ 2 follows from the fact that P g t = M/γ in the simulation coordinates, plus the transformation to the BH frame proper time that appears in Ṁ .We can see that, using Eq. ( 7) to transform the ADM momentum to the BH rest frame, the terms in Ṁ cancel exactly such that In the BH frame the steady state is such that dP φ x /d t = 0, so In practice in our simulations we cannot resolve the field close to the horizon, but only down to ∂Σ i .However, it is easy to check that for a ∂Σ i close enough to the BH event horizon (e.g., Fig. 5 of Ref. [67]), we have ∂Σi That is, the left-hand side is approximately independent of the choice of Σ i provided it is close to the horizon.This is the quantity that we extract from our simulations and plot as the total BH drag force, in both the current and previous work [67].
It is straightforward to find the critical impact parameter b cr (v), which separates between the particles that are absorbed and the ones that are only deflected by the BH; this is given in Eq. (24).The drag force due to the accretion of particles is simply given by The contribution to the DF from the particles that are deflected is found from introducing an infrared cutoff b max in the impact parameter of the beam, and where π − 2ϕ ∞ is the deflection angle, which is given by [100] (B3) with K and F , respectively, the complete and the incomplete elliptic integrals of the first kind, with elliptic modulus The parameter μ := M/ l (with l the latus rectum) and the eccentricity ẽ satisfy (B5) The general-relativistic effects on the deflection are nonnegligible for particles with impact parameter b ≫ b cr /v.Thus, for large clouds b max /λ dB ≫ β (≫ 1), one can find a semi-analytic expression with the (numerically-evaluated) factor and is accurate to < 1% for all v.For a cutoff in the radius r instead (see, e.g., Ref. [9]) the total gravitational drag force on BHs moving through collisionless particle-like media (with no velocity dispersion) is given by Eq. ( 23).

Appendix C: Comparison to previous results
Here we discuss the differences between the analytic expressions used in this work and the phenomenological ones (based on fluid-like results) used to fit the numerical results in our previous work [67].There we fit our results using an expression inspired by Bondi accretion [70,71], which describes the accretion of a collisional fluid onto a BH as it moves through the medium and has the form with λ an O(1) constant factor and c s the asymptotic sound speed of the medium.The error bars in our numerical results, due to the long-lived oscillations in the scalar field, were not small enough to reach any definite conclusions on the appropriateness of the description in terms of Bondi accretion.In particular, considering v ≫ c s , we found that any value of λ ∈ [0, 1] was a good fit to our results.The increase in the overall (DF + accretion) drag force that we observed compared to the analytic expressions was then attributed to a pressure-like correction to the DF expression D(ᾱ s ≪ 1, Λ ≫ 1) in Eq. ( 21) (note that we also did not include the Lorentz factor γ in Λ and β).Soon after, analytic expressions for the relativistic drag force on BHs moving through light scalars were derived in Ref. [68], which found results consistent with ours for the two lightest cases (α s = 0.05 and 0.2) we had considered.However, in that reference it was argued that, while Bondi accretion is a good approximation for collisional fluids, The red dotted curve represents the expression in Eq. ( 44) of Ref. [67], including the "pressure correction" and using Bondi accretion, and the black dashed curves show the expression for the drag force using Unruh accretion (rather than Bondi) without any phenomenological correction.We see that whilst the previously suggested pressure correction fits the simulation data well at large r, it does not provide a good fit for smaller clouds, whereas the expression for Unruh accretion is a good fit in both regimes (for small αs).
the momentum accretion of a light (wave-like) scalar field is better described by an expression derived by Unruh in Ref. [73], which reads This adds a non-negligible contribution to the force at relativistic velocities, which removes the need for the (ad-hoc) pressure correction in the DF.It has also the advantage of being fully derived from first principles, containing no free parameters.As shown in Fig. 4, the force on the BH from Unruh accretion (red dotted line), unlike the one from Bondi's (green dashed line), is negligible compared to DF at low velocities, but has a larger effect at high velocities.The Unruh accretion term is particularly relevant for small clouds (Λ ≲ 1).In that regime it is also not equivalent to the pressure term that was suggested in Ref. [67].This is illustrated in Fig. 5, where the red dotted curves show the expression in Eq. ( 44) of Ref. [67], which includes the "pressure correction" factor P corr = 1 + κp/ρ (with p a scalar field effective pressure and κ another fudge factor) and the Bondi-type accretion.Whilst it provides a good fit for a large cloud Λ ≫ 1 (top panel), when Λ ≲ 1 (bottom panel) the expressions in this paper, which account for the accretion correctly, clearly give a better fit.

Appendix D: Numerical set-up
The analytic expression we wish to verify is only valid for small β, so we take α s = 0.05 and 0.025, although we also compare to previous results for α s = 0.5.We use the fixed background feature of the open source numerical relativity code GRChombo [101,102] to solve Eqs. ( 8) and ( 9) on a fixed metric background in the boosted isotropic coordinates described above.The scalar field is evolved using the method of lines with its spatial derivatives evaluated using fourth-order finite difference stencils, a fourth-order Runge-Kutta time integration, and a hierarchy of grids with 2:1 resolution.The value of the metric and its derivatives are calculated locally from the analytic expressions at each point.
Due to the causal structure of the BH in isotropic Schwarzschild coordinates, the metric naturally imposes ingoing boundary conditions at the horizon.At the outer boundary, we implement nonzero, time oscillating boundary conditions for the scalar field by setting the field to be spatially constant in the radial direction, extrapolating from the values within the numerical domain.This simulates the effects of a roughly constant energy density, but can introduce unphysical effects in very long simulations.These effects can be easily identified by varying the domain size.We find that these start to appear on a timescale of roughly twice the size of the numerical domain (in code units), which ultimately limits the time for which the growth of the cloud can be studied.
We choose the size of the simulation domain L such that the approximate boundary conditions that we use do not spoil the solution near the BH in the time necessary to reach a steady state.Therefore, we take L = 4096M , with 9 (2 : 1) refinement levels and the coarsest level having 128 3 grid points.We take advantage of the symmetries of the problem in Cartesian coordinates to reduce the domain to 64 2 ×128 points.This ensures that we maintain a spatial resolution of dx = 0.0625M around the horizon of the BH at r ∼ M/2 and that the BH horizon lies well within the area covered by the finest level.In order to resolve the temporal oscillations of the field sufficiently on the coarsest level, we use a time step dt coarse such that we have at least 32 time steps per period of oscillation, i.e.T = 2π/µ > 32 dt coarse .
The simulations used to extract the numerical results in this work are for the most part (with the addition of α s = 0.025) the same as the ones in Ref. [67].Therefore we do not repeat these here, but instead refer to Appendix A.2 in that paper, where we showed convergence tests in

FIG. 1 .
FIG.1.Our numerical results for large clouds Λ ≫ 1 (shown by the vertical error bars) considering three scalar field masses αs = 0.05, 0.5 and 1, from top to bottom.The dot-dashed and dotted curves represent, respectively, the analytic expressions for the total force in the light(21) and heavy (23) scalar limits.Unruh accretion reproduces the numerical results in the case of αs ≪ 1 (top panel), but overestimates them for the other cases where αs ≲ 1 (middle and bottom panels).Interestingly, the numerical results for these masses seem to be well described by the combination D( ᾱs ≪ 1, Λ ≫ 1) + A( ᾱs ≫ 1) (shown as black dashed curves), that is, the accretion is in the particle limit, whilst the dynamical friction is in the wave limit.

FFIG. 2 .
FIG.2.Our numerical results for small clouds Λ ≲ 1 (indicated by the vertical error bars) for the scalar masses αs = 0.05 (top panel) and αs = 0.025 (bottom panel).The dashed curves show the analytic expression for the total force in Eq. (22).We see that the wave-like expressions for the accretion and dynamical friction are an excellent fit.

FIG. 3 .
FIG. 3. Summary of the different DF and accretion regimes, showing the regions of validity of the expressions in Eqs.(20)-(23).In the case of a scalar field wavelength comparable to the size of the BH and smaller and very small scalar field clouds (shaded region on the bottom left of this graph), there is no valid analytic expression for the DF force.Here we can only calculate the drag force on the BH numerically, but as one can see from our results, the analytic expressions give a good order of magnitude estimate even in this regime.

FIG. 4 .
FIG.4.The different contributions to the total drag force on the BH, comparing the addition of Bondi accretion (as used in our previous work) and Unruh accretion (as in[68]) given by the green dashed and red dotted lines, respectively.Comparing with the DF force alone (grey solid line) we can see that Bondi accretion has a non-negligible contribution to the force at low velocities, while Unruh accretion only becomes important at very relativistic speeds.All three cases are computed with the same parameters αs = 0.05 and µr = 45.We set λ = 0.5 for the O(1) factor in the expression for Bondi accretion in Eq. (C1).

45 FFFIG. 5 .
FIG. 5. Our numerical results (shown as vertical error bars)for αs = 0.05 in the case of a large cloud (µr = 45) on the top panel, and a small cloud (µr = 2.5) on the bottom panel.The red dotted curve represents the expression in Eq. (44) of Ref.[67], including the "pressure correction" and using Bondi accretion, and the black dashed curves show the expression for the drag force using Unruh accretion (rather than Bondi) without any phenomenological correction.We see that whilst the previously suggested pressure correction fits the simulation data well at large r, it does not provide a good fit for smaller clouds, whereas the expression for Unruh accretion is a good fit in both regimes (for small αs).