Evolution of diffuse scalar clouds around binary black holes

The use of modern effective field theory techniques has sparked significant developments in many areas of physics, including the study of gravity. Case in point, such techniques have recently been used to show that binary black holes can amplify incident, low-frequency radiation due to an interplay between absorption at the horizons and momentum transfer in the bulk spacetime. In this paper, we further examine the consequences of this superradiant-like mechanism on the dynamics of an ambient scalar field by taking the binary's long-range gravitational potential into account at the nonperturbative level. Doing so allows us to capture the formation of scalar clouds that are gravitationally bound to the binary. If the scalar is sufficiently light, the cloud can be diffuse enough (i.e., dilute while having considerable spatial extent) that it engulfs the binary as a whole. Its subsequent evolution exhibits an immensely rich phenomenology, which includes exponential growth, beating patterns, and the upscattering of bound states into scalar waves. While we find that these effects have negligible influence on the binary's inspiral in the regime wherein our approximations are valid, they offer new, analytic insight into how binary black holes interact with external perturbations. They may also provide useful, qualitative intuition for interpreting the results from future numerical simulations of these complex systems.


I. INTRODUCTION
The details of how binary black holes evolve and coalesce are today well established [1][2][3][4][5][6][7][8]. Indeed, our ability to extract gravitational-wave signals from data gathered by the LIGO Scientific and Virgo Collaborations [9][10][11] is a testament to how accurately we understand these systems according to general relativity. That being said, while we now have a strong grasp of how binary black holes behave when isolated in empty space, their dynamical response to external perturbations is less well understood.
The presence of some matter distribution on top of this spacetime generically introduces multiple new scales into the problem, which can lead to a myriad of new effects. As a first step, recent work in addressing this question has focused on the ramifications of perturbing a binary black hole with a Klein-Gordon field. 1 While a scalar field is a natural starting point due to its technical simplicity, it is also of particular interest on phenomenological grounds, as the existence of new ultralight fundamental fields is a generic prediction of string theory [23][24][25][26][27]. Because these hypothetical particles may couple only very weakly to the Standard Model, prospects for their detection rest on finding novel and innovative probes. In this regard, the advent of gravitational-wave astronomy is particularly opportune [28][29][30][31].
On one end of the spectrum, it is now known that ultralight bosons can rapidly form "clouds" around a rotating black hole due to a superradiant instability, which is most pronounced when the field's Compton wavelength is comparable to the gravitational radius of the hole . If such a cloud forms around a black hole belonging to a binary system, Baumann et al. [54][55][56][57] have shown that this scalar (or vector) cloud can lead to a large dephasing of the gravitational-wave signal when it undergoes a resonant transition. These resonances are excited during the early inspiral phase, when the cloud's characteristic size λ is much smaller than the binary's orbital separation a. For field configurations with even larger spatial extent, Bernard et al. [58] recently reported the existence of global quasinormal modes, which were observed in numerical simulations of a scalar Gaussian pulse of size λ ∼ a scattering off a binary. Lastly, working analytically in the long-wavelength limit (λ a), I recently demonstrated that binary black holes can amplify incident, low-frequency radiation under certain conditions [59].
These advances notwithstanding, there is still much to learn about how binary black holes interact with external fields and whether these interactions can be used to search for new physics. The goal of this present paper is to further our understanding of this problem in the long-wavelength limit by incorporating the effects of the binary's long-range gravitational potential at the nonperturbative level. In doing so, we will be able to study not just the scattering of radiation but also the formation and evolution of scalar clouds diffuse enough to engulf the binary as a whole. As we will show, the same mechanism that led to the amplification of radiation discovered in Ref. [59] triggers, among other things, the exponential growth or decay of these clouds.
At the heart of this phenomenon is the absorptive nature of the binary's constituents: Because an ambient field must always satisfy purely ingoing boundary conditions at the horizons, it is "dragged" alongside the black holes as they orbit one another. Consequently, while a part of this scalar field irreversibly crosses the horizons and deposits its energy into the individual black holes, what remains in the bulk spacetime is agitated by the binary's motion and can either gain or lose momentum as a result. Together, these two effects facilitate a direct exchange of energy and momentum between the binary and the scalar field, which can proceed in either direction depending on their relative phase velocities. Given this description, the close parallel between this process and the reflection of electromagnetic waves off moving conductive boundaries [60][61][62] is unsurprising. The underlying mathematics also shares many features in common with superradiant phenomena in arXiv:2004.03570v1 [hep-th] 7 Apr 2020 general [32,33,63], making it natural to regard this mechanism as a novel variant of superradiance fueled by the binary's orbital motion. Accordingly, we will refer to this mechanism as orbital superradiance.
The remainder of this paper proceeds as follows: In Sec. II, we begin by reviewing the low-energy effective field theory (EFT) constructed in Ref. [59] to describe the propagation of a long-wavelength scalar field on a binary black hole spacetime. Owing to the inherent separation of scales, this problem is analytically tractable and a perturbative solution is obtained in Sec. III. The result is then discussed in two stages: In Sec. IV, we track the evolution of a scalar cloud and calculate the rate at which different bound states grow or decay as a result of orbital superradiance. Then, in Sec. V, we show that the periodic forcing exerted by the binary inevitably converts a fraction of these bound states into outgoing scalar radiation. Also in this section, we reanalyze the effect of mode mixing on the amplification of scalar waves. Finally, a brief summary of our key results is presented in Sec. VI. Supplementing this main line of discussion are Appendices A-D, which contain additional technical details and some of the lengthier derivations.

II. EFFECTIVE FIELD THEORY
Our goal is to study the evolution of a real Klein-Gordon field φ(x) around a binary black hole. In realistic scenarios, the energy density in this field is always expected to be sufficiently dilute that its backreaction may be neglected as a first approximation. Even in this test-field limit, however, obtaining a solution is prohibitively difficult for current analytic methods due to the complexity of the spacetime; hence, further assumptions are necessary to render the problem tractable.
In this paper, we restrict our attention to scalar-field configurations whose characteristic length scale λ is much greater than the binary's orbital separation a. Because the individual black holes cannot be resolved by such a longwavelength field, the binary as a whole behaves like an effective point particle that couples to the scalar via a set of multipole moments. This coarse-grained description of the system is essentially an extension of the multipolar post-Minkowskian formalism [64] with the addition of a scalar field.
The interactions between this effective point particle and the fields living in the bulk are encoded in the action [59] where M = M 1 + M 2 is the total mass of the binary. The scalar field φ(τ) ≡ φ( ( (z(τ)) ) ) is to be evaluated at the position of the binary's barycenter, which travels along the worldline z µ (τ) with proper time given by dτ = −g µν dx µ dx ν . In the second term, 2 the composite operators O L (τ) are symmetric and tracefree (STF) tensors localized on the worldline that capture how this point particle's internal degrees of freedom interact with the long-wavelength scalar [65][66][67]. The ellipsis in (2.1) alludes to the presence of analogous composite operators that couple to the gravitational field [68], although these will play no role in our discussion.
Induced multipole moments.-Generically, each of the composite operators in (2.1) can be decomposed into two pieces [69]: , the first term is a permanent multipole that is present if the internal degrees of freedom can directly source the scalar field. Such terms are present, for instance, were the action in (2.1) to describe a binary neutron star system in a scalar-tensor theory of gravity [70][71][72][73][74][75][76]. A binary black hole in general relativity, however, is not generic because stringent no-hair theorems stipulate that Kerr black holes cannot source real scalar fields [77][78][79][80][81][82][83][84][85][86]; , which represents the binary's multipole moments that are induced in response to external perturbations. From now on, we drop the subscript (R) to declutter our notation.
As φ(x) is assumed to behave like a test field, the solution for these induced multipoles is given by linear response theory. With retarded boundary conditions imposed [65,66,69], (2.2) Written in this quantum mechanical language, the expectation value · · · above requires the notion of some Hilbert space for the internal degrees of freedom. When specified correctly, this formalism can be used to systematically incorporate quantum effects like Hawking radiation [87]. That being said, in this work we will be interested only in purely classical observables. It is then unnecessary to specify the density matrix with which this expectation value is taken, as the commutator is simply a c-number when Planck-suppressed terms are neglected [87].
By matching this post-Minkowskian formulation of the binary to a post-Newtonian description valid in its near zone, Ref. [59] showed that the classical solution to (2.2) is in the nonrelativistic, low-frequency limit. In writing this solution, we have chosen coordinates such that the binary's barycenter is at rest at the origin. The motion of the Nth black hole, whose horizon has area A N , is then given by the vector z N (t). Expansion parameters.-At this stage, it is worth enumerating the four different expansion parameters that appear in this EFT. It will be convenient to assume that the two black holes have comparable masses for the purposes of power counting, although this formalism remains valid for arbitrary mass ratios as long as the binary's constituents are widely separated.
The first two expansion parameters come from the post-Newtonian description of the binary in its near zone: The solution in (2.3) is only the leading term in a series organized as an expansion in the binary's orbital velocity v ∼ GM/a and the ratio of timescales GMω. To be specific, the latter is the ratio of the black holes' light-crossing times to the characteristic time scale ω −1 of the scalar. The fact that (2.3) depends only on the areas of the black holes but not on their spins is a consequence of working to leading order in these parameters [59]. The induced multipoles are generated as a result of absorption of the scalar across the horizons [67], and a black hole's absorption cross section in the low-frequency limit is s-wave dominated and equal to its area [88].
The approximation of the entire binary as an effective point particle introduces an additional two expansion parameters: the post-Minkowskian parameter GM/λ and the ratio of length scales a/λ. The first of these characterizes the nonlinearity of a given term in the solution due to self-interactions of the gravitational field. In this work, we work to first order in GM/λ and will treat it nonperturbatively in order to capture the bound states of φ(x), but otherwise we neglect all higher-order corrections in GM/λ.
As for the fourth expansion parameter, we will-rather unusually from the point of view of an EFT-treat terms with different powers of a/λ on equal footing, despite higherorder terms being parametrically suppressed. Doing so will allow us to keep track of the mixing between different angularmomentum modes, which leads to interesting consequences.
These four parameters control different aspects of the perturbative expansion, but enforcing the two conditions v 1 and a/λ 1 is often sufficient to ensure that we are inside the EFT's regime of validity. For a scalar cloud that is gravitationally bound to a binary, its characteristic frequency is set by the scalar field's mass, ω µ, up to some nonrelativistic binding energy ∼ GM µ/λ dB . In the denominator is the de Broglie wavelength λ dB ∼ (GM µ 2 ) −1 , which determines the characteristic length scale of the cloud. This second expression may be used to recast the condition a/λ 1 into an upper bound for the scalar's mass; namely, which follows after using v 2 ∼ GM/a and v 3 ∼ GMΩ. As a rough guide, (2.4) says that a scalar field should have a mass µ 10 −11 (v/0.1)(M /M) eV if it is to engulf a binary of total mass M with orbital velocity v.
An additional upper bound must be established for freelypropagating scalar waves that impinge on the binary with frequency ω and momentum k ∼ 1/λ. For high-momentum modes, we choose to replace the necessary condition a/λ 1 with the sufficient condition aω 1 for simplicity. The latter equivalently reads For low-momentum modes, the quantity a/λ can be arbitrarily small, so now the ultraviolet (UV) cutoff for this EFT is set by another expansion parameter; namely, GMω. Since ω ∼ µ in this limit, the condition GMω 1 is equivalent to µ/Ω v −3 . This upper bound is weaker than that of (2.4), thus we are guaranteed to remain in the EFT's regime of validity when the UV cutoffs in (2.4) and (2.5) are both respected.
Equation of motion.-Extremizing the total effective action with respect to φ(x), we obtain denotes the wave operator on flat space and we have also included the leading contribution from the binary's gravitational potential on the lhs. On the rhs, the induced multipoles O L (t) are given by (2.3); hence, this is a linear, homogeneous differential equation for φ(x).
It is worth remarking that the delta function on the rhs of (2.6) inevitably leads to singularities. Likewise, singularities also arise from the operators O L (t), which are functions of the scalar field and its derivatives evaluated at the origin. These UV divergences, which in this case originate from the point-particle approximation of the binary, are commonplace in EFTs and can be dealt with in the usual way by utilizing a convenient regulator and a renormalisation scheme [89,90]. This procedure turns out to be unnecessary, however, up to first order in perturbation theory.

III. PERTURBATIVE SOLUTION VIA GREEN'S FUNCTION
We obtain an approximate solution to (2.6) by treating the interaction terms on the rhs as small perturbations. Denoting the differential operator on the lhs by D x , this entails looking for a solution of the form φ = φ (0) + φ (1) + · · · , where the zeroth-order piece is an exact solution to the noninteracting theory, D x φ (0) (x) = 0, while the first-order correction is (3.1) The first term is the particular integral sourced 3 by the induced multipoles O where the superscript (0) indicates that the expression in (2.3) is to be evaluated using the zeroth-order solution φ (0) . Meanwhile, the second term in (3.1) is the complementary function, D x φ (1) cf (x) = 0, whose inclusion may be necessary to ensure that the overall solution satisfies our choice of boundary conditions.
In this section, we discuss the three ingredients that make up the particular integral. We begin by writing down the general solution φ (0) to the noninteracting theory, which is then fed into (2.3) to obtain explicit expressions for O L (t). Finally, taking their convolution with the retarded Green's function G(x, x ) produces the end result. 3 More accurately, the terms on the rhs of (2.6) should be viewed as sink terms, since the induced multipoles arise from absorption.

A. The noninteracting theory
To establish some nomenclature and introduce the basis for our perturbative approach, we begin in this subsection with a brief review of the Coulomb functions [91,92].
As the noninteracting theory is time-translation invariant and spherically symmetric, we may look for solutions of the form φ(x) ∝ R(r)Y m (x)e −iωt , where Y m (x) are the usual spherical harmonics. The radial part R(r) must then satisfy the differential equation The resulting set of solutions can be divided into three categories depending on the value of this quantity.
In what follows, we define k ≡ k(ω) as the appropriate root of k 2 , namely Radiation modes.-Let us begin with the case k 2 ≥ 0. Defining ζ −GM µ 2 /k, two linearly-independent solutions to (3.2) are where H ± are Coulomb functions. (We follow the conventions in Ref. [93].) From their asymptotic forms at large r, given in (A1), we can deduce that these solutions correspond to ingoing (−) and outgoing (+) spherical waves. For later purposes, it will also be useful to have defined a particular linear combination of these radiation modes. Let where F is another Coulomb function. This solution describes a superposition of ingoing and outgoing waves in equal measure and is regular at the origin as a result. Additionally, let us define the mode functions which describe the three-dimensional spatial profile of these scalar waves. Yukawa modes.-Solutions to (3.2) for the case k 2 < 0 may now be obtained by analytic continuation. Represented in terms of Whittaker functions [93], they read where the Coulomb phase shift σ (ζ) and the Gamow factor C (ζ) are given by [91,93] For imaginary k defined according to (3.3), the R − (k, r) solution is seen to grow exponentially with r and is therefore unphysical. In contrast, R + (k, r) describes a nonpropagating field profile with characteristic size λ ∼ 1/|k |. Depending on the value of k ∈ iR >0 , this solution can either be singular or regular at r = 0.
The set of singular solutions includes the ω = = 0 mode, which has the asymptotic form at large r up to some constant prefactor. One easily recognizes this as the Yukawa potential sourced by a point charge at the origin, albeit with corrections coming from the gravitational potential of the point mass M. In general, we will refer to this set of singular solutions as the Yukawa modes. Because these modes correspond to having pure imaginary k in the continuous domain iR >0 modulo a discrete set of points to be discussed below, the term continuum states will be used to refer to the combined set of radiation modes and Yukawa modes.
Bound states.-The solution R + (k, r) is regular at the origin when k takes special values such that −iζ = n is an integer and n ≥ + 1. These regular solutions have corresponding frequencies ω = ±E n given by which is reminiscent of the hydrogen bound-state spectrum. Accordingly, this set of regular solutions will be called the bound states. To make the connection to the hydrogen atom even more explicit, let us define R n (r) to be a rescaled version of the solution R + (k, r) when evaluated at ζ = in. 4 In terms of the Whittaker functions, it reads The second line follows from Eq. (13.14.32) of Ref.
[93], which shows that R + and R become proportional to one another when ζ → in. Note also that −ik = GM µ 2 /n in this limit, thus R n (r) is a real-valued function. With this definition, the mode functions are exactly the (orthonormal) hydrogen wavefunctions, albeit with GM µ in place of the fine-structure constant.
General solution.-Because we are treating absorption of the scalar field by the black holes perturbatively via interaction terms, the origin is devoid of sinks or sources in the noninteracting theory. Consequently, the zeroth-order solution φ (0) must be regular at r = 0. This boundary condition precludes the existence of Yukawa modes at this order and, moreover, a net flux of radiation into or out of the origin is also prohibited. The general solution is thus given by the linear combination The first term is the sum over a superposition of ingoing and outgoing waves [the factor of 2 follows from (3.4b) and (3.5)], where a given mode (ω, , m) has an ingoing amplitude specified by the coefficient I > ω m . The ">" symbol is used to emphasize that this function can be chosen without loss of generality to have support only in the domain ω ≥ 0, corresponding to positive-frequency modes. The negative-frequency modes are then automatically taken into account by the complex conjugate (c.c.) terms. Further note that I > ω m must vanish for ω ∈ (−µ, µ) as per the boundary conditions described above. Meanwhile, the second term in (3.12) is the sum over bound states, with a conventional prefactor of 1/ √ 2µ included to render the coefficients c

B. Induced multipoles for circular orbits
The zeroth-order solution in (3.12) may now be fed into (2.3) to obtain the binary's induced multipoles. To that end, we begin by introducing some compact notation: Let u ≡ (n, , m) collectively refer to the three integers that specify a bound state and, likewise, let w ≡ (ω, , m) refer to the parameters for a given continuum state. We then write u ≡ n, ,m and w ≡ ,m dω 2π (3.13) to denote summing over the bound and continuum states, respectively. In this new shorthand, the zeroth-order solution reads (3.14) Rather than substitute this directly into (2.3) to produce a set of tensorial objects, it is more convenient to work with the components of O L (t) obtained via projection onto a basis of STF tensors. For each , let Y m L denote the basis vectors that generate the spherical harmonics [i.e., Y m (x) = Y m Lx L ] and satisfy the orthogonality relation [94] The 2 + 1 independent degrees of freedom of O L (t) can then be obtained via the projection where the prefactor of −4πi is included purely for convenience.
To reconstruct O L (t), one simply inverts this relation to find At the moment, the formula in (2.3) for these induced multipoles does not make any assumptions about the black holes' trajectories, apart from requiring that |z N (t)| λ. For simplicity, in this paper we will restrict our attention to circular orbits with frequency Ω oriented such that its angular momentum points along the positive z axis. For this configuration, Ref. [59] showed that the components O m (t) are given by where d is the unit vector parallel to z 1 (0) and characterizes the interaction strength between the black holes and the scalar, with r 1 = aM 2 /M and r 2 = −aM 1 /M denoting the displacements of the black holes from their barycenter. The induced multipoles O L (t) are necessarily real by construction; hence, the definition in (3.17) can be used to show that its components must satisfy the constraint which follows from the identity for the complex conjugate of a spherical harmonic; cf. (A12). This motivates writing It is then easy to show that for a real scalar-field solution of the form φ(x) = φ > (x) + c.c., one obtains O > m (t) by evaluating the rhs of (2.3) using only φ > (x) rather than φ(x). In other words, we find after substituting only the positive-frequency part of (3.14) into (3.18). In obtaining this result, we have made use of the identities in (A10) and have defined explicit expressions for which are provided in (A11). The result in (3.22) can now be used to determine the first-order correction φ (1) via the method of Green's functions.

C. Integration contours
The Green's function in (3.1) is defined by the equation and may be written as the inverse Fourier transform (3.25) where the radial part is given by [95,96] G ω (r, r ) = (−2ik) with r > max(r, r ) and r < min(r, r ).
To evaluate the particular integral in (3.1), we first note that the delta function imposes the restrictions r < = r and r > = r, which make performing the integral over x relatively straightforward. Integrating by parts to move the derivatives Now using (3.15), (3.17), and (A9), this can be shown to simplify to after also using the identity n! ≡ n!!(n − 1)!! for the double factorial. Given that (2 )!! = 2 !, an equivalent expression for the particular integral is Both expressions will turn out to be useful in later sections. It remains to perform the integrals over t and ω. Care must be exercised with the latter because the Green's function contains poles at ω = ±E n and branch points at ω = ±µ and at infinity [95]. To proceed, we split the integral over ω into two parts: Its principal value along the real line gives rise to the continuum states, while the bound states are obtained by integrating over closed contours encircling each of the poles; see also Fig. 1.
The complete first-order solution is thus a sum of three parts: The first term contains the bound states, which we study in Sec. IV, while the second contains the continuum states; discussed in Sec. V. Finally, recall from our earlier discussion that the third term is a complementary function that may be required to satisfy boundary conditions.

IV. BOUND STATES
The bound states surrounding a binary evolve in an intricate manner as a result of their interaction with the black holes. In perturbation theory, this evolution can be regarded as being sourced by the induced multipoles O L (t). The first-order result φ (1) b is obtained by performing the ω integral in (3.28) over closed contours that each encircle one of the poles ω = ±E n of the Green's function.
When t > t , these contours should be closed in the lower half of the complex plane; see Fig. 1. In the limit δ → 0 and L → ∞, the integrals along the vertical paths cancel each other, while the integral over the bottom horizontal path vanishes. The residue theorem can then be applied to show that where the poles have been shifted by an amount −i to enforce retarded boundary conditions and the sum over s = ±1 is used to account for both the positive and negative-frequency solutions.
Note also the presence of the step function θ(t − t ), which follows from the fact that the contour should be closed in the upper half of the complex plane when t < t .
To determine the residue of the gamma function at the pole ω = sE n , we use the standard Laurent expansion Γ(− j + z) = (−1) j /( j!z) + O(z 0 ) valid for any nonnegative integer j to obtain Substituting this back into (4.1) and using (3.10) to rewrite the Whittaker function in terms of R n (r), one finds which can be further simplified to The latter expression follows after making two observations: First, the s = −1 terms are exactly the complex conjugates of the s = +1 terms, which one can show by using the identities in (3.20) and (A13) together with the freedom to relabel m → −m as it is being summed over. Second, the lower bound of the t integral at −∞ yields no contribution because of the −i term in the exponent, thus we need only keep track of the result from the upper bound. Discarding this lower bound constitutes no loss in generality, as the freedom to specify initial conditions for the amplitudes of the bound states at some initial time, say t = 0, is provided by our freedom to choose the complementary function φ (1) cf . The result in (4.4) is not yet in a useful form because O (0) m (t ) contains both the positive and negative-frequency parts of the zeroth-order solution, whose separate contributions we would like to make manifest. To do this, we use the decomposition in (3.21) to write The freedom to swap the terms involving O It is now apparent that the full solution φ = φ (0) + φ (1) + · · · for the bound states has the form cf,u + · · · . The time evolution of the amplitude is given to first order by while c (1) cf,u is a constant term coming from φ (1) cf that we tune in order to choose initial conditions.

A. Mode mixing
We are now in a position to discuss the physical implications of this result. Substituting the expression for O which is written in terms of the energy differences and the matrix elements Note the extra normalization factor of 1/(2µ) in the second line of (4.9) has been included because the sum over continuum states w is dimensionful; cf. (3.13). Written in this way, V uu and V uw both have dimensions of energy.
We have chosen to call the objects in (4.11) matrix elements because-in analogy with quantum mechanics-they characterize the overlap between different modes as a result of their interaction with the binary. From the general form in (4.9), it is clear that even if a given mode u has zero amplitude c u (0) = 0 at an initial time t = 0, the presence of another bound state u u will seed the growth of u as long as V uu 0. Likewise, energy from ingoing radiation can also be captured and converted into bound states, as the second line in (4.9) demonstrates. This appearance of mode mixing is unsurprising in the context of a binary, as the underlying spacetime is not time-translation invariant nor is it rotationally symmetric. 5 Accordingly, the modes u ≡ (n, , m) and w ≡ (ω, , m) do not remain eigenstates of the interacting theory.
Actually, this last statement needs refining because the matrix elements vanish when at least one of + m or + m are odd. This property can be traced back to the spherical harmonics in (4.11), which are evaluated using the unit vector d that is confined to be in the z = 0 plane. As was already pointed out in Ref. [59], the vanishing of these matrix elements has a simple physical interpretation: Modes with + m 2Z correspond to field profiles that are concentrated away from the z = 0 plane and are therefore unappreciable in the neighborhood of the binary. As a result, these modes are effectively blind to the presence of the black holes and are conserved (at this order in perturbation theory), while modes with + m ∈ 2Z interact with the binary and get mixed into one another.
Let us highlight another consequence of mode mixing: As the solution in (4.9) shows, a given mode u will oscillate not just at its natural frequency E n but also at the secondary frequencies |E n − ∆ uu |, which results in a beating pattern when viewed in the time domain. Crucially, note that some these secondary frequencies are much greater than the scalar field's mass µ. If a particle were to have an energy given by one of these high frequencies, we would expect it to escape the gravitational potential of the binary and travel off to infinity. Indeed, the same thing happens in the case of a long-wavelength scalar field, as we show later in Sec. V.

B. Growth rates
For now, we turn our attention to yet another prediction of (4.9). In addition to mode mixing, this solution also reveals that bound states can either grow or decay. To see this, let us suppose that c Because the constant term may be removed by an appropriate choice of c (1) cf,u , the physical effect of the diagonal element V uu is to cause the bound state to grow at the rate It is worth remarking that the linear growth in (4.12) is an approximation that is valid only at early times t 1/Γ u . Once Γ u t becomes of order unity, terms of the form ∼ (Γ u t) p that appear at higher orders in perturbation theory all become relevant. One might naturally expect that resumming these polynomials to all orders will lead to an exponentially-growing solution c u (t) ∝ exp(Γ u t) and, indeed, this turns out to be the case. The details of this resummation procedure, while interesting on theoretical grounds, have been relegated to Appendix B as they will not be relevant to this paper's main line of discussion. In what follows, it will suffice to work with the linear approximation in (4.12).
Written out explicitly, the growth rate for the u ≡ (n, , m) mode reads (4.14) For reasons that have already been discussed, this growth rate vanishes when + m is odd. When nonvanishing, a given mode will grow if 0 < E n < mΩ and will decay otherwise. Accordingly, the bound states of the noninteracting theory turn into quasibound states once their interaction with the binary is taken into account.

C. Backreaction and energy extraction
In many well-motivated scenarios [24][25][26][27], an ultralight scalar couples only very weakly to the Standard Model. As such, if a binary black hole is enveloped by a cloud of quasibound states, we should not expect to observe the evolution of this scalar cloud directly and must infer its presence through more indirect means. One possibility is by examining the way it affects the binary's orbital evolution. Because any momentum transferred from the black holes to the scalar field must be accompanied by an appropriate backreaction of the scalar onto the black holes, 6 any secular increase in the energy E b of the bound states must have been extracted from the energy stored in the binary.
To gain a sense for how much energy is extracted via this process, let us consider a simplified scenario in which only a single modeû ≡ (n,ˆ ,m) is populated initially. Time-averaged over a period much longer than the other timescales in the problem, the rate at which energy is extracted from the binary into this bound state is given by [42][43][44]  u | 2 is the total (initial) energy in the scalar cloud. Note that the timeaveraging procedure eliminates any contribution from mode mixing at this order.
If Γû > 0, the growth of this mode extracts energy from the binary and causes it to inspiral more rapidly than it would in pure vacuum. On the other hand, if Γû < 0, this decaying mode injects energy into the orbit while it is being absorbed and will decelerate the inspiral as a result. Whether either of these effects leave an observable imprint depends on the magnitude of (4.15). A useful measure is to compare it to the energy flux of gravitational waves emitted by the binary, given to leading order by [98] F = 32 5 for the case of circular orbits. Taking the ratio of (4.15) to (4.16), we find where E orb = GM 1 M 2 /2a is the magnitude of the binary's orbital energy and 6 This backreaction can be understood in terms of a force that the scalar exerts onto the black holes [67]. 7 At leading order, a number of results from the literature on the superradiant instability of a single Kerr black hole can be adapted to the study of orbital superradiance because the mathematics describing quasibound states in the long-wavelength limit is identical.
is the rate at which the orbit shrinks due to gravitational radiation.
A cursory glance at (3.19) and (4.14) will confirm that, for fixed component masses M 1 and M 2 , the scalar cloud's growth rate is largest when both black holes are spherical. In this limit, the ratio Γ u /Γ GW can be expressed as a function of three dimensionless quantities: the binary's orbital velocity v, the symmetric mass ratio ν M 1 M 2 /M 2 , and the ratio µ/Ω that relates the scalar's mass to the binary's orbital frequency. An explicit formula is provided in Appendix D, where we also argue that the precise value of ν has little effect on our conclusions. For this reason, we consider only equal-mass binaries in what follows.
The value of Γ u /Γ GW as a function of µ/Ω is shown in Fig. 2 for an equal-mass binary composed of spherical black holes. A value of v = 0.1 has been chosen for the orbital velocity, which is large enough that it is at the limit of validity of the post-Newtonian expansion. As the curves in Fig. 2 would all move downwards for smaller values of v, they represent the largest-possible rates that we can reliably calculate using this EFT.
For small values of µ, the u = (2, 1, 1) mode has the largest growth rate, which reaches a maximum of Γ 211 /Γ GW ∼ 2 × 10 −22 before it crosses the threshold at µ Ω (note E n µ) and turns into a decaying mode. Above this threshold, the (3, 2, 2) mode takes over as the fastest-growing mode until it too becomes a decaying mode at the next threshold µ 2Ω. This pattern continues for increasing values of µ, with ( + 1, , ) being the fastest-growing mode when ( − 1) µ/Ω . The overall trend in Fig. 2 clearly shows that the maximum value Γ u can attain decreases rapidly as µ increases. Although the EFT breaks down as we approach µ/Ω ∼ v −2 , this trend strongly suggests that orbital superradiance is always grossly inefficient. Consequently, the exponential growth of a long-wavelength scalar cloud is unlikely to leave any measurable impact on the evolution of a binary black hole.
In contrast, the decay rates can become much larger as µ increases (see the right panel of Fig. 2), but their observational viability rests on Mû being comparable to E orb . The fact that the growth rates are so small implies that clouds with such high densities are unlikely to have formed dynamically around binary black holes that start off in pure vacuum.
Having said that, there may be a possibility that other processes could generate these clouds, particularly in alternative theories of gravity wherein the scalar field is nonminimally coupled to matter. For instance, does core collapse of a massive star into a black hole remnant leave behind an appreciable scalar cloud? If so, could successive supernova events in a stellar binary lead to a pair of black holes enveloped by a common scalar cloud (assuming an optimal value for µ)? It has been shown that a large amount of scalar radiation can be produced during core collapse in a certain class of scalar-tensor theories [99][100][101][102], although current numerical methods are unable to determine if a scalar cloud can develop around a black hole remnant [103]. Exploring these open questions presents an exciting opportunity for future work.

V. OUTGOING RADIATION
The periodic forcing that a binary exerts onto a surrounding cloud inevitably leads to a fraction of the scalar field being upscattered and ejected as outgoing radiation. Additionally, ingoing radiation can scatter off this binary and undergo amplification when given the right initial conditions. In our perturbative approach, both of these phenomena are encoded in the principal value of the ω integral in (3.29). Using (3.6) to replace the Whittaker function with the radial solution R + (k, r), the result is where the first-order correction to the outgoing amplitude for a given mode w ≡ (ω, , m) is It is worth highlighting that the solution in (5.1) is not of the form φ = φ > + c.c.; hence, A (1) w implicitly accounts for both positive and negative-frequency modes. Because the overall solution must be real, these coefficients will have to satisfy a constraint analogous to (3.20). Combining (A18a) with the freedom to relabel ω → −ω and m → −m as they are being integrated and summed over, respectively, we arrive at the constraint E SW at which energy is carried away from the binary in the form of scalar waves, it is useful to first recast the zeroth-order solution (3.14) into a form similar to (5.1). The identity in (A18b) can be used to rewrite the sum over continuum states as where I w ≡ I ω m = I > ω m + e πζ (−1) +m I > * −ω −m implicitly accounts for the ingoing amplitudes of both positive and negative-frequency modes, and can be seen to satisfy the requisite constraint Taken in combination with (5.1), the full solution φ = φ (0) + φ (1) + · · · for the continuum states has the form where R w I w + A w is the total outgoing amplitude, with A w given to first order in the interactions in (5.2a). Now integrating the (t, r) component of the scalar's energymomentum tensor over a spherical shell of radius r and taking the limit r → ∞, the time-averaged power loss is given by the difference between the energy flux flowing into and out of the system, These quantities have simple expressions when integrated over all time:

A. Ejection of bound states
To better understand the physical implications of (5.1), let us start by supposing-as we did in Sec. IV C-that there is no ingoing radiation and only a single bound stateû is populated at zeroth order. In this case, the outgoing radiation we compute represents the portion of the scalar cloud that is being ejected out of the system. The energy flux for this process is given to leading order by where the delta function in the denominator is associated with the integral over all time, dt ≡ 2πδ(0). Because this formula is quadratic in A (1) w , it is formally of second order in the interactions, thus . E SW will generally be smaller than the rate . E b at which energy extracted from the binary fuels the growth of bound states; cf. (4.15). However, this hierarchy inverts when the scalar field is sufficiently light, as we will show.
Since ζ and, consequently, σ (ζ) are real when k 2 > 0 [93], taking the absolute square of (5.2a) yields (5.9) One might recognise S (ζ) ≡ |s (ζ)| 2 as the Sommerfeld enhancement factor [104][105][106], which may be rewritten as after using standard identities for the gamma function [104]. On the mathematical level, this factor arises naturally in our calculations because the interaction terms involve evaluating derivatives of Coulomb functions at the origin. We run into difficulties, however, when attempting to assign to this factor its usual physical interpretation. We elaborate further in later parts of this section. When only a single bound stateû is populated, a close inspection of (3.22) reveals that mode mixing will generate continuum states with frequencies in a discrete set given by ω m = |En + (m −m)Ω|. As we argued at the end of Sec. IV A, these continuum states are generated alongside newly populated quasibound states u û, which oscillate at the same set of frequencies |E n − ∆ uû | ≡ |En + (m −m)Ω|. The subset of these continuum states with ω m > µ are radiation modes that propagate to infinity. Substituting (3.22) into (5.8), the power in these radiation modes is where k m ≡ k(ω m ) and likewise ζ m ≡ ζ(ω m ). (The intermediate steps for this general type of calculation are presented in Appendix C.) As it stands, this result is not particularly illuminating. It is perhaps most instructive to compare (5.11) to (4.15), in which case one finds Although this is a sum over infinitely many modes, it suffices to keep only the lowest few values of to obtain a good approximation because the prefactor |Y m (d)/(2 +1)!!| 2 decays rapidly like ∼ 1/(2 +1)!. To proceed, let us begin by analyzing the limiting behavior of a given term in (5.12) when k 2 m > 0 and ζ m 1. The sum of Sommerfeld factors reads 8 which has the asymptotic form ∼ 2π|ζ | 2 +1 when ζ → ∞. While it may be natural to want to think of this as describing the usual Sommerfeld enhancement for low-momentum modes (recall ζ ∝ 1/k), in the present context there is no analogous process that occurs on flat space, since bound states cannot form in the absence of the binary's gravitational potential. With this in mind, the Sommerfeld factors appearing in (5.12) are perhaps best regarded as simply an inevitable part of the result rather than engendering any kind of enhancement.
Taking the limit ζ m 1, the corresponding term in (5.12) reduces to since ω m µ in this case. Power counting reveals that this term scales with the EFT's expansion parameters as ∼ v 5 (a/λ dB ) 2 +1 (m − µ/Ω). Accordingly, for clouds witĥ m ∼ O(1) and µ/Ω v −2 , cf. (2.4), the rate at which energy is carried away by low-momentum radiation is parametrically suppressed relative to . E b . Indeed, each term in (5.12) is a monotonically-increasing function of k m , so most of the energy is carried away in highmomentum modes (ζ m 1). In this limit, S (ζ) ∼ 1, thus the corresponding term in (5.12) reduces to To assess the typical size of this term, it is instructive to express it in terms of the ratio f µ/Ω. Also using the definition ω m = |En + (m −m)Ω| and approximating En µ, where ∆m = m −m. This term scales as ∼ v 2 +6 when f , m, and ∆m are all of order unity, meaning the rate at which energy is carried away by high-momentum radiation is-in this case-also parametrically suppressed relative to . E b . However, when f 1, the 1/ f prefactor enhances this term such that power loss to radiation can become significant in comparison to . E b . More precisely, a given high-momentum mode will extract energy from the cloud at a rate greater than (5.17) One may conclude from this simple scaling analysis that scalar clouds cannot form dynamically around binary black holes when the scalar field's mass µ is sufficiently light, as the rate at which the cloud is depleted via scalar radiation is greater than the rate at which it grows due to orbital superradiance. As a rough guide, this occurs when µ 10 −19 (v/0.1) 3 (M /M) eV.
Gravitational waves.-It is worth briefly remarking that a scalar cloud will also emit gravitational waves due to the oscillatory nature of its backreaction onto the spacetime. When only the single bound stateû is populated (and further assuminĝ =m for simplicity), the energy flux of gravitational waves emitted by the cloud is [42,44] .
where Cû (< 1) is some dimensionless prefactor whose exact form will not be important to us, but note C n00 = 0. Comparing this with (4.15), we find

B. Superradiant scattering
We now turn our attention to a different scenario in which the zeroth-order solution is given by a steady stream of ingoing radiation peaked at the single modeŵ ≡ (ω,ˆ ,m). This corresponds to making the choice where Φŵ is in general some complex-valued coefficient with dimensions of energy. It has previously been shown that this ingoing wave can extract energy from the binary's orbital motion and undergo amplification under the right conditions [59]. In this subsection, we extend the results of Ref. [59] in several directions. 9 This is a conservative upper bound that does not take into consideration the possibility that B 2 ˆ can vanish. As an example, if an equal-mass binary is surrounded by a cloud comprised of only thê u = (2, 1, 1) mode, one finds B 01 = 0; hence, energy is predominantly radiated away in the = 1 modes. In this case, . E SW becomes larger than . E b only when µ/Ω v 8 .
Amplification factor.-From (3.22), we learn that a single ingoing modeŵ will scatter into multiple outgoing modes with frequencies in a discrete set given by ω m = |ω + (m −m)Ω|. Included in this spectrum is the original mode with frequencŷ ω ≡ ωm, which typically comprises the majority of the outgoing energy flux as a result of it interfering with ingoing radiation. Explicitly, one expands (5.7) to find Now substituting (3.22) into the above formula and making use of the symmetries in (5.3) and (5.5), we obtain (5.21) to first order in the interactions. (The details of this calculation are presented in Appendix C.) As the Gamow factor C (ζ) ≡ s * (ζ)/(2 + 1)!! when ζ ∈ R, the above expression is already real and can be further simplified to read To gain a sense of how much energy is exchanged during this scattering process, it is instructive to compare (5.22) to the total flux of ingoing radiation . E in SW ; an expression for which is also derived in Appendix C, cf. (C14). This dimensionless ratio defines the total amplification factor Z . E SW / . E in SW , which is given by in the case of a single ingoing modeŵ. For high-momentum modes withk GM µ 2 (ζ 1), the Sommerfeld factor Sˆ (ζ) and the exponential e −2πζ both reduce to unity such that we recover the result in Ref. [59]. This limiting behavior signifies that the binary's long-range gravitational potential has negligible influence on the amplification or absorption of high-momentum modes. For more general values ofk, (5.10) may be used to show that which has the asymptotic form ∼ |ζ | 2 +1 e −2π |ζ | when ζ → ∞. This leads to an exponential suppression of the energy carried away by low-momentum modes. Interestingly, one might naively expect that the appearance of S (ζ) in (5.2a) should lead to Sommerfeld enhancement, but this is directly contradicted by the result in (5.24), which we interpret as follows: In the classical analog of Sommerfeld enhancement, we imagine a stream of particles impinging on a star of radius R [106]. In the absence of gravity, the geometric cross section σ 0 = πR 2 of this star provides a measure of the fraction of particles that collide into it and are subsequently absorbed. However, the actual cross section σ for this interaction can be much larger, especially for particles with low momenta, because the star's attractive gravitational potential is able to pull in particles that have impact parameters greater than R .
With this picture in mind, one should now expect no analogous enhancement to occur in the present scenario. In our setup, the ingoing modeŵ is a spherical wave that is already directed straight at the origin; hence, the presence of the binary's gravitational potential does nothing to affect the amount of radiation that reaches it. While this argues for the lack of Sommerfeld enhancement, it remains to explain the suppression of low-momentum modes observed in (5.24). Although the physical origin of this suppression is still not fully understood, the most likely explanation is that it is due to the conversion of radiation modes into bound states [cf. the second line in (4.9)], which is enhanced at low momenta. A full quantitative analysis is needed to validate this interpretation, although such a task is beyond the scope of this present paper.
Putting these conceptual issues aside, let us discuss the likelihood of observing this energy exchange between the binary and the scalar. As we discussed earlier in Sec. IV A, the particular geometry of the binary prevents it from interacting with any long-wavelength mode whose angular momentum is such that + m 2Z. For the remaining modes, amplification occurs if 0 <ω <mΩ, in which case the binary loses energy and inspirals more rapidly as a result. Otherwise, there is a net absorption of the scalar by the binary, which then gains energy and experiences a slowing down of its inspiral. The feasibility of observing either of these effects depends on the magnitude of . E SW when compared to the outgoing flux of gravitational radiation F. As a rough estimate, we should expect to observe the influence of this energy exchange on the orbital motion only if the ratio is not too much smaller than unity. The amplification factor Z is shown as a function of the ingoing frequency for different values of the scalar field's mass in Fig. 3. As we did in Sec. IV C, we assume an equal-mass binary composed of spherical black holes traveling with orbital velocity v = 0.1. The same reasoning as before justifies limiting ourselves to this specific case: First, the curves in Fig. 3 would all move downwards for smaller values of v, so once again they represent the largest possible values that can be reliably calculated using this EFT. Moreover, the precise value of the binary's mass ratio has little effect on our overall conclusions, as Appendix D argues.
In the left panel of Fig. 3, we see that amplification is most pronounced for the = m = 1 mode, which reaches a maximum value Z ∼ 4 × 10 −10 when ω = 3Ω/4. As the higher = m modes are less efficiently amplified, the overall trend suggests that this orbital superradiant mechanism continues to become increasingly insignificant even for large frequencies ω Ω/v beyond the EFT's regime of validity. Given the smallness of Z and the unlikelihood that the ingoing flux . E in SW of scalar waves would match or exceed the outgoing gravitational-wave flux F in realistic astrophysical scenarios, we deduce that the amplification of long-wavelength scalar fields is observationally inaccessible.
In contrast, the right panel of Fig. 3 demonstrates that absorption continues to become more efficient as ω increases; naturally prompting us to ask: Is there any regime (possibly at some frequency ω Ω/v outside the EFT's regime of validity) in which |Z | and . E in SW are both large enough that they can leave a measurable imprint on the evolution of the binary? On a more theoretical level, it is also interesting to ask: What is the maximum amount of radiation that can be absorbed by a binary black hole in a given time? Because the amplification factor is bounded from below (Z ≥ −1), there are two possibilities for what might occur in the high-frequency regime: Either Z gradually tends to a minimum value (meaning absorption would be most pronounced at high frequencies) or it has a turning point (i.e., there is a critical frequency beyond which absorption becomes less efficient again). The fact that moving black holes can amplify high-frequency radiation [107] (via what is essentially the slingshot effect) is a hint that the latter may be more likely. These questions point to potential directions for future work.
Secondary modes.-To complete our discussion on the scattering of scalar waves, we ought to discuss the additional energy that is carried away by the secondary modes w ŵ generated through mode mixing. Although their contribution to the energy flux is typically subleading because they first appear in . E SW at second order in the interactions, the energy they carry can exceed that of the original mode ifω is sufficiently small, as we now show.
Let us denote this O(A 2 ) correction to the energy flux as .

E
(2) SW . The calculation is almost identical to that in Sec. V A and the end result is found to be SW , we simply divide one by the other to find Observe that (5.28) has the same mathematical structure as (5.12); hence, the analysis will proceed in a largely similar fashion. First, recall that while sums of this kind are to be taken over infinitely many modes, a good approximation can be obtained by keeping only the lowest few values of , since the higher multipoles are factorially suppressed. Next, the fact that each term in (5.28) is a monotonically-increasing function of k m signifies that most of the energy carried away will be in the form of high-momentum modes (ζ m 1). Therefore, let us concentrate on a given term in (5.28) and suppose that ζ m 1. In terms of the dimensionless ratiô f ω/Ω, this term reads where ∆m = m −m. Whenf ,m, and ∆m are all of order unity, this term scales as ∼ v 2 +6 and is thus parametrically suppressed. However, if insteadf 1 (and necessarily µ/Ω <f ifŵ is to be a radiation mode), this term can become arbitrarily large. As a result, the energy carried away in a secondary mode of frequency ω m will dominate over the energy carried by the original modeŵ whenω/Ω v 2 +6 ≤ v 6 .
This phenomenon is particularly interesting ifŵ is a counterrotating mode satisfyingmΩ −ω < 0, since in this case we would predict an amplification factor Z < 0 when truncating to O(A). However, if this ingoing mode hasω v 6 Ω, the outgoing energy flux is dominated by the O(A 2 ) term, which is positive definite; cf. (5.26). Thus, we learn that energy is always extracted from the binary during this kind of scattering process if the ingoing frequency is low enough.
To be clear, while SW can be very large relative to SW , its magnitude is still small in comparison to . E in SW and further decreases asω → 0; meaning the actual amount of energy that a low-frequency, counter-rotating mode extracts from a binary is always negligible. Nonetheless, this calculation illustrates the kinds of rich physics that can arise as a consequence of mode mixing.

VI. DISCUSSION
While we now have a comprehensive picture of how binary black holes evolve when in isolation, questions about their dynamical response to external perturbations are stimulating an emerging area of active research. The benefits to be reaped from this enterprise are twofold: First, the study of how ultralight fields influence the orbital evolution of a binary provides us with the prospect of using gravitational-wave detectors as tools to search for new physics. Second, even in the absence of a discovery, this kind of theoretical work offers new insight into general relativity and the properties of gravitational systems.
In this paper, we tracked the evolution of a long-wavelength scalar field living on a fixed binary black hole background. The interplay between absorption at the horizons and momentum transfer in the bulk gives rise to a novel energy-extraction mechanism, which we have herein dubbed "orbital superradiance." It was previously shown that this mechanism can lead to the amplification of incident, low-frequency radiation [59]. The main novelty in this work is a nonperturbative treatment of the binary's long-range gravitational potential, which facilitates an extension of the results in Ref. [59] to include the formation and evolution of bound states.
The key takeaway is as follows: Consider for simplicity an incident spherical wave that is peaked at some frequencyω and has angular momentum specified by the integers (ˆ ,m). Three effects are triggered when this scalar wave scatters off a binary black hole. First, the ingoing mode is reflected back out with an amplitude that is either amplified or reduced. Amplification occurs if it corotates with the binary at an angular phase velocityω/m smaller than the binary's orbital frequency Ω, while a net absorption of the mode occurs otherwise. Second, owing to the inherent lack of symmetries in this system, multiple secondary modes are generated during scattering, which propagate outwards at certain frequencies given by ω m = |ω + (m −m)Ω|. Third, a fraction of this ingoing wave is captured and converted into (quasi)bound states. It is worth noting that the binary's "dumbbell" geometry establishes a selection rule that requires the integerˆ +m to be even if any of these effects are to take place, since field configurations that defy this condition are unappreciable in the neighborhood of the binary. Likewise, only secondary outgoing modes and bound states with angular momenta satisfying + m ∈ 2Z are generated during this scattering process.
The bound states that form around a binary subsequently evolve in an intricate manner due to a combination of three effects: First, orbital superradiance drives each bound state to either grow or decay exponentially, depending on its angular momentum. Second, the underlying geometry of the spacetime allows different modes to mix into one another, causing each bound state to exhibit a beating pattern by virtue of oscillating at multiple frequencies. For a scalar field of mass µ, the bound states with angular momentum ( , m) oscillate at frequencies given approximately by ω m | µ + (m − m )Ω|. Third, a fraction of these bound states are inevitably upscattered and ejected out of the system as scalar waves. (There is also a concomitant emission of gravitational waves from the cloud, although this is typically a subleading effect.) The rate at which outgoing scalar radiation depletes the energy in the scalar cloud can exceed the growth rates of the bound states when µ is considerably smaller than Ω, in which case scalar clouds can no longer form dynamically around the binary.
All of these effects illustrate the rich phenomenology that can arise in systems with horizons (or dissipative channels, more broadly) when time-translation invariance and rotational symmetry are weakly broken. 10 Furthermore, the calculations underpinning these results demonstrate the usefulness of modern EFT methods in understanding the dynamics of complex systems with multiple hierarchies of scales. Unfortunately, they predict that orbital superradiance is grossly inefficient: The energy extracted from a binary black hole to amplify incident scalar radiation or to fuel the growth of bound states is always negligible for systems within the EFT's regime of validity. Moreover, the trends in Figs. 2 and 3 suggest that this conclusion extends also to scalar fields with higher frequencies or larger masses, as long as the binary is in its early inspiral phase. 11 While this is certainly disappointing from an observational standpoint, our results still constitute useful information about which effects play an important role during a binary's lifetime. Besides, orbital superradiance is expected to be just one of many phenomena that arise when a binary is perturbed by an external field. There is still much to do before a comprehensive 10 In the sense that both symmetries are restored in the EFT upon removal of the interaction terms, which we treated perturbatively. 11 Our perturbative approach breaks down when the binary is closer to merger, although the way the growth rates and amplification factors scale with the orbital velocity suggests that they may become appreciable in this regime. While this may be true, the binary does not remain in this stage for long and the usual version of black hole superradiance quickly takes over once the binary coalesces. survey of all the effects that can occur is in hand. Natural next steps include relaxing some of the assumptions made in this paper. In particular, truncating to leading order in perturbation theory led us to neglect any interaction between the scalar field and the spins of the individual black holes. However, it has been shown that ambient matter generically exerts a "gravitational Magnus force" on spinning black holes [108]; hence, the presence of an external scalar field is likely to have an effect on the precession of the binary's orbital plane. Additionally, it is conceivable that generalizing to the case of eccentric orbits or including gravitational radiation from the binary will also teach us something new about how these systems interact with external fields. That being said, it is important to temper our expectations for observing any effect we study in this long-wavelength limit, since the large separations of scales inherent in this regime typically lead to strong powerlaw suppression by the EFT's expansion parameters.
Indeed, all work to date (including the upward trends in the right panels of Figs. 2 and 3) point to the likelihood of scalar fields with higher frequencies or larger masses having a more dramatic impact on the orbital evolution of a binary black hole; especially when resonant excitations can occur [56,58]. Even so, we believe that studies in the long-wavelength limit will continue to be of value moving forward. The fact that we have a strong analytic handle on the problem in this regime can be used to gain physical intuition for better interpreting the results of numerical simulations, which may be the only recourse in certain scenarios. Of particular interest, for example, is the case of a scalar field whose characteristic size is comparable to the binary's orbital separation. The general ideas and techniques found in this paper could also find applications in other branches of physics that involve open systems wherein one or more spacetime symmetries are weakly broken.

ACKNOWLEDGMENTS
It is a pleasure to thank Vitor Cardoso, Anne-Christine Davis, Eugene Lim, and Ulrich Sperhake for stimulating discussions.

APPENDIX A: PROPERTIES OF THE RADIAL SOLUTIONS
This Appendix provides a collection of useful identities for the radial solutions to (3.2).
Limiting forms.-The identities in this first part have all been reproduced or adapted from Ref. [93]. At large distances (r → ∞), the R ± solutions have the asymptotic forms If instead r → 0, the behavior of the radial solutions around the origin may be inferred from the limiting forms of the Whittaker functions. The solutions that are regular at the origin are all proportional to whereas the irregular solutions are proportional to It is also useful to know the limiting behavior of these solutions for small and large values of ζ. For ζ → 0 with k held fixed (corresponding to a removal of the gravitational potential), one has where h ± are the spherical Hankel functions while j is the spherical Bessel function of the first kind. Instead taking the low-momentum limit k → 0 (i.e., ζ → ∞ with M and µ held fixed), one recovers the usual Bessel functions: Note that the prefactor multiplying R + on the lhs of (A7) has exactly the same k dependence as the outgoing amplitude A Since both R (k, r) and R n (r) are proportional to the Whittaker function M, a good starting point is which follows from (A3) and the identity ∂ L x L = !. The definitions in (3.6) and (3.10) can then be used to show that where the coefficients on the rhs are given by Complex conjugates.-Several instances in the main text exploit identities for the complex conjugates of the mode functions to obtain simplified expressions. These identities are derived here. Combining the well-known identity with the fact that R n (r) is a real function tells us that the complex conjugate of a bound-state mode function is As for the continuum states, the identity [109] [ can be used to show that These can be written in a more useful form by utilizing the circuital relations [109,110] in conjunction with the identity k * (ω) ≡ −k(−ω), which is a consequence of the definition in (3.3). Only the results for R + and R are relevant for physical solutions. They are: [R + (k, r)] * ≡ e +iπ( +iζ) R + (k, r), where we writek ≡ k(−ω) andζ ≡ ζ(−ω) as shorthand.
Combined with (A12), the complex conjugates of the continuum-state mode functions are

APPENDIX B: RESUMMATION AND LATE-TIME BEHAVIOR OF BOUND STATES
The growth of Γ u t in (4.12) invalidates our naive perturbative approach once it becomes of order unity, even though the expansion parameters enumerated in Sec. II all remain small. This kind of secular growth turns out to be generic in any system with an interaction Hamiltonian that persists for all times [111]. For the scenario studied in this paper, this late-time breakdown of perturbation theory poses no threat because the binary will typically coalesce well before Γ u t ∼ 1. That being said, on theoretical grounds, it is interesting to explore how we might obtain an approximate solution to (2.6) that remains valid at late times. The general results may find application in studies of other open systems whose lifetimes exceed 1/Γ u .
The key is to carefully resum the dominant polynomial behavior ∝ t p at each order p in perturbation theory while neglecting subleading terms. Included in this set of terms we will neglect are backreaction effects from the outgoing radiation and Yukawa modes (see Sec. V), because they contribute to c u (t) beginning only at second order. Additionally, higherorder corrections to the formula for the induced multipoles in (2.3), which are suppressed by extra powers of v and GMω, can also be neglected. For added simplicity, we will also assume no ingoing radiation in this Appendix.
With these considerations in mind, resummation amounts to looking for a solution to (2.6) of the form φ(x) = ∞ p=0 φ (p) (x), where each term in this series is sourced by the previous term via the iteration This is, of course, simply a generalization of (3.1). Accordingly, the coefficients for the bound states at order p + 1 are given by after suitably generalizing (4.9). Rather than perform this string of integrals, the trick is to now differentiate twice to obtain In the same way that (B1) is an iterative solution to the equation of motion in (2.6), this set of differential equations in (B3) can be viewed as establishing an iterative method (assuming V is suitably small) for solving the master equation 12 This equation is strongly reminiscent of time-dependent perturbation theory in quantum mechanics, albeit with two key differences: First, the terms on the rhs can be regarded as arising from some interaction Hamiltonian for the system, whereby u|H int |u ∝ −iV uu e −i(m−m )Ωt . Given that the diagonal elements V uu are real, the prefactor of −i indicates that H int is not Hermitian-a necessary condition for this system 12 With the benefit of hindsight, this master equation can be seen to follow more easily from substituting the ansatz φ(x) ∝ u [c u (t)ψ u (x)e −iE n t + c.c.] directly into (2.6). However, doing so obscures the fact that the bound states alone are not a complete solution to the problem. As we discussed in Sec. V A, the production of outgoing radiation is inevitable in this system. to exhibit nonunitary evolution. The other key difference is that (B4) is clearly a set of second-order, rather than first-order, differential equations; reflecting the relativistic nature of this system.
To perform the requisite resummation, we now treat the terms involving the diagonal elements V uu on the rhs of (B4) nonperturbatively. Moving them over to the lhs, the master equation may be rewritten as where is independent of c u and can therefore be considered a source term. To solve this equation, we begin by noting that it is of the form where γ ± are the two zeros of the characteristic polynomial γ 2 ± − (2iE n − V uu )γ ± + 2iE n Γ u . As V E n , it suffices to use the approximate solutions where we defineū ≡ (n, , −m) as shorthand. For later purposes, it will also be useful to define γ u (γ + − γ − )/2. Now choosing boundary conditions such that c u (0) = c (B9) This is only a formal solution because J u depends on the other bound states u u, whose solutions are also given by (B9). To obtain an explicit result, we iterate (B9) in powers of V/γ 1. Starting with c u (t) = c (0) u e Γ u t + O(V/γ), after one iteration we find whereẼ u = E n + iΓ u is the complex frequency of the quasibound state and∆ uu is defined in the same way as ∆ uu in (4.10a) except with E n →Ẽ u . Having carefully resummed the leading polynomial growth to all orders, this solution is valid at late times t 1/Γ u while still being organized as a perturbative expansion in the small parameter V/γ. This resummed solution also brings with it a new prediction: Let us denote the fastest-growing mode by u and its growth rate by Γ . The sum over u in (B10), which quantifies the leading effects due to mode mixing, then tells us that all modes satisfying V uu 0 will grow at the same rate Γ at late times, even if they were initially decaying. 13 While this is a nontrivial result for the system of equations under study, it is irrelevant in the case of a scalar cloud around a binary black hole because even the shortest e-folding time 1/Γ is always orders of magnitude greater than the orbital decay timescale 1/Γ GW . It would therefore be interesting to explore if there are other open systems that could survive long enough to exhibit this universal growth rate at late times.

APPENDIX C: SCALAR-WAVE FLUX
Multiple instances in Sec. V call for the evaluation of the power radiated to infinity in scalar waves. For the sake of efficiency, we here derive a general formula for the energy flux when it is sourced by a single bound stateû or a single ingoing radiation modeŵ. In either of these cases, the components of the induced multipoles take on the general form w . Written out explicitly, the total energy lost to scalar radiation is First order.-Let us begin by evaluating the term in (C5) that is linear in A To proceed, we deduce from the definition of s (ζ) in (5.2b) that s (−ζ) = s * (ζ)e πζ when ζ ∈ R (k 2 > 0). Combined with the freedom to relabel ω → −ω and m → −m as they are being integrated and summed over, respectively, we find (C8) Of course, this general result is valid only for induced multipoles sourced by a single bound stateû or a single ingoing radiation modeŵ. If the former, we would have I * w = 0 ∀ w; meaning A (1) w contributes linearly to the energy flux only if it interferes with ingoing radiation. Performing the sum over w in (C8), we pick up a single nonvanishing contribution at the frequency ωm ≡ω, which yields .
Second order.-In the absence of ingoing radiation, the leading contribution to the energy flux is quadratic in A (1) w . Taking the absolute square of (C4), one finds A (1) w 2 = S (ζ)k 2( +1) (2π) 2 δ(ω − ω m )δ(0)|o m | 2 + δ(ω + ω −m )δ(0)|o ,−m | 2 . (C10) To arrive at this result, we use the fact that the cross terms proportional to δ(ω−ω m )δ(ω+ω −m ) may be discarded because they have nonoverlapping support. 14 It then follows that after also using the freedom to relabel m → −m. To simplify this result one step further, we note that the product ωk 2 +1 is invariant under the transformation ω → −ω whereas ζ changes sign; hence, we may easily perform the integral over ω to obtain E in SW is often a useful measure. For a single ingoing modeŵ, the latter is given by substituting (5.19) into (5.7b). After neglecting cross terms that involve products of delta functions with nonoverlapping support, the end result is where the exponential arises from the identity in (5.5).

APPENDIX D: DIFFERENT MASS RATIOS
For a binary composed of spherical (sph) black holes, the growth rate of the u ≡ (n, , m) mode may be written as in terms of the dimensionless parameters v, ν M 1 M 2 /M 2 , and f µ µ/Ω. Similarly, the amplification factor for the w ≡ (ω, , m) mode is In this case, the expression depends on four dimensionless quantities: v, ν, f µ , and f ω ω/Ω; and note ζ ≡ −v 3 ( f 2 ω / f 2 µ − 1) −1/2 . In both formulas, the effect of the symmetric mass ratio ν enters via the same function n (ν) 8ν 2 (1 + which is normalized such that n (1/4) = 1. The additional prefactor of 1/ν in (D1) causes the ratio Γ u /Γ GW to diverge in the limit ν → 0. This singularity is unrelated to Γ u and is due entirely to Γ GW . Physically, it is reflecting the fact that the timescale over which the orbit shrinks becomes infinite in the limit of a test particle around a host black hole. It follows that Γ u and Z are both proportional to n (ν); hence, our discussion will focus purely on this function's properties.
It presents three different classes of behavior depending on the value of . When = 0, the largest value of n 0 (ν) = 2 − 4ν in the domain ν ∈ (0, 1/4] coincides with ν = 0. This behavior has a simple physical interpretation: For a binary with fixed total mass M, a smaller symmetric mass ratio leads to a larger combined area for the black hole's horizons. In other words, that n 0 (ν) is maximized when ν = 0 is simply corroborating the fact that absorption is more efficient when there is a larger horizon area. Note, however, that the value of this function only changes by a factor of 2 in the domain ν ∈ (0, 1/4]. When = 1, one finds n 1 (ν) = 16ν 2 , which is maximized when the binary's components have equal masses (ν = 1/4). For ≥ 2, this function always has a maximum somewhere in the domain ν ∈ (0, 1/4]. Numerically, we find max n (ν) ∼ exp(1.4 − 2.0 log − 1.5) when 1, which can be quite a large number. For instance, max n (ν) ∼ 10 9 when = 20.
What does this mean for the conclusions in the main text? Given that n 1 (ν) is maximized for equal-mass binaries, the growth rates for the = 1 modes (shown in Fig. 2) and the amplification factors for the same modes (shown in Fig. 3) are indeed the largest values possible within the EFT's regime of validity. For larger values of , carefully selecting an optimal value for ν can enhance the growth rates and amplification factors relative to the equal-mass case, but this enhancement grows exponentially with at best, which is still no match for the factorials in the denominators of (D1) and (D2). Consequently, the general trend remains unchanged: The maximum value that Γ u or Z can attain decreases rapidly as we increase µ or ω.