We show that the existence of clouds of ultralight particles surrounding black holes during their cosmological history as members of a binary system can leave a measurable imprint on the distribution of masses and orbital eccentricities observable with future gravitationalwave detectors. Notably we find that for nonprecessing binaries with chirp masses inlineformulammlmath displayinlinemmlmrowmmlmi mathvariantscriptMmmlmimmlmommlmommlmn10mmlmnmmlmsubmmlmrowmmlmiMmmlmimmlmrowmmlmrowmmlmo stretchyfalsemmlmommlmrowmmlmsubmmlmrowmmlmathinlineformula formed exclusively in isolation largerthanexpected values of the eccentricity i.e. inlineformulammlmath displayinlinemmlmrowmmlmiemmlmimmlmommlmommlmsupmmlmrowmmlmn10mmlmnmmlmrowmmlmrowmmlmommlmommlmn2mmlmnmmlmrowmmlmsupmmlmrowmmlmathinlineformula at gravitationalwave frequencies inlineformulammlmath displayinlinemmlmrowmmlmsubmmlmrowmmlmifmmlmimmlmrowmmlmrowmmlmiGWmmlmimmlmrowmmlmsubmmlmommlmommlmsupmmlmrowmmlmn10mmlmnmmlmrowmmlmrowmmlmommlmommlmn2mmlmnmmlmrowmmlmsupmmlmtext mmlmtextmmlmtext mmlmtextmmlmiHzmmlmimmlmrowmmlmathinlineformula would provide tantalizing evidence for a new particle of mass between inlineformulammlmath displayinlinemmlmrowmmlmo stretchyfalsemmlmommlmn0.5mmlmnmmlmommlmommlmn2.5mmlmnmmlmo stretchyfalsemmlmommlmommlmommlmsupmmlmrowmmlmn10mmlmnmmlmrowmmlmrowmmlmommlmommlmn12mmlmnmmlmrowmmlmsupmmlmtext mmlmtextmmlmtext mmlmtextmmlmieVmmlmimmlmrowmmlmathinlineformula in nature. The predicted evolution of the eccentricity can also drastically affect the inband phase evolution and peak frequency. These results constitute unique signatures of boson clouds of ultralight particles in the dynamics of binary black holes which will be readily accessible with the Laser Interferometer Space Antenna as well as future midband and decihertz detectors.
['Cosmology', 'Astrophysics', 'and Gravitation', 'General Relativity and Quantum Cosmology', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena', 'High Energy Physics - Phenomenology', 'High Energy Physics - Theory']
Signatures of Ultralight Bosons in the Orbital Eccentricity of Binary Black Holes
{'Signatures of ultralight bosons in the orbital eccentricity of binary black holes': "Mateja Boˇskovi'c, 1 Matthias Koschnitzke, 1, 2 and Rafael A. Porto 1  \n1  \nDeutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany 2 II. Institut fur Theoretische Physik, Universitat Hamburg,  \nLuruper Chaussee 149, 22761 Hamburg, Germany  \nWe show that the existence of clouds of ultralight particles surrounding black holes during their cosmological history as members of a binary system can leave a measurable imprint on the distribution of masses and orbital eccentricities observable with future gravitational-wave detectors. Notably, we find that for nonprecessing binaries with chirp masses M ≲ 10 M ⊙ , formed exclusively in isolation, larger-than-expected values of the eccentricity, i.e. e ≳ 10 -2 at gravitational-wave frequencies f GW ≃ 10 -2 Hz, would provide tantalizing evidence for a new particle of mass between [0 . 5 , 2 . 5] × 10 -12 eV in nature. The predicted evolution of the eccentricity can also drastically affect the in-band phase evolution and peak frequency. These results constitute unique signatures of boson clouds of ultralight particles in the dynamics of binary black holes, which will be readily accessible with the Laser Interferometer Space Antenna, as well as future mid-band and Deci-hertz detectors.  \nIntroduction. The birth of gravitational-wave (GW) science [1] heralds a new era of discoveries in astrophysics, cosmology, and particle physics [2]. Measuring the properties of GW signals with current and future observatories, such as the Laser Interferometer Space Antenna (LISA) [3], the Einstein Telescope (ET) [4] and Cosmic Explorer (CE) [5], as well as other Mid-band [6] and Decihertz detectors [7], not only will unravel the origins of binary black hole (BBH) mergers, it also opens the possibility to discover (very-weakly-coupled) ultralight particles that are ubiquitous in theories of the early universe [812]. Notably, the mass, spin alignment, and eccentricity are expected to be correlated with formation channels, isolated or dynamical , e.g. [13-31]; whereas boson clouds (or 'gravitational atoms' [8, 9]), formed around black holes via superradiance instabilities [32-36], can produce a large backreaction on the orbital evolution. Following analogies with atomic physics [37], the cloud may encounter Landau-Zener (LZ) resonances [38], or ionization effects [39-41]. The presence of a cloud then leads to large finite-size effects [37, 42], floating/sinking orbits [38], as well as other sharp features [40], that become unique signatures of ultralight particles in the BBH dynamics.  \nFor the most part, up until now backreaction effects have been studied under the simplified assumption of planar, quasi-circular orbits. The reason is twofold [37]. First, several formation scenarios lead to spins that are parallel to the orbital angular momentum [18]. Second, the decay of eccentricity through GW emission in vacuum [43, 44] is expected to have circularized the orbit in the late stages of the BBH dynamics. We retain here the former but relax the latter assumption. As we shall see, adding eccentricity not only introduces a series of overtones [41, 45, 46], it can also have a dramatic influence in the orbital dynamics as the cloud transits a LZ-type transition. Although the strength of the new resonances is proportional to the eccentricity, depending on their position and nature (floating or sinking), a small departure from circularity can lead to transitions that not only would deplete the cloud, but also  \ninduce a rapid growth of eccentricity toward a large critical (fixed-point) value: e cr ∈ [0 . 3 , 0 . 6]. As measurements of the eccentricity are correlated with formation channels, the predicted increase can impact the inferred binary's origins. Measurements of larger-than-expected eccentricities would then provide strong evidence for the existence of a new ultralight particle in nature. In particular, because of the critical fixed point, a fraction of the BBHs undergo a rapid growth of orbital eccentricity to a common value. As a result, the distribution of masses and eccentricities may feature a skewed correlation by the time they reach the detector's band. Furthermore, for chirp masses M < 10 M ⊙ and spin(s) aligned with the orbital angular momentum-expected to exclusively form in the field-the presence of a boson cloud at earlier times can shift a fraction of the population toward values of e ≳ 10 -2 at 10 -2 Hz, readily accessible to LISA [3]. Furthermore, the GW-evolved eccentricity may also be within reach of the planned mid-band [6] or Deci-hertz [7] observatories. For all such events, a new ultralight boson of mass [0 . 5 , 2 . 5] × 10 -12 eV forming a cloud and decaying through a LZ-type transition prior to detection, may be the ultimate culprit.  \nThe more drastic evidence is given when the resonant transition occurs in band with measurable frequency evolution. A plethora of phenomena are discussed in [37, 38] for circular orbits. In addition to overtones, the increase in eccentricity would imply that higher harmonics become more relevant, which in turn affects the peak frequency of the GWs, even for floating orbits. We point out here some salient features and elaborate further on the details elsewhere [46].  \nThe gravitational atom. Ultralight particles of mass µ can form a cloud around a rotating black hole of mass M , via superradiant instabilities [8, 9]. The typical mass of the (initial) cloud scales as M c , 0 /M ≃ α , whereas its typical size is r c ≃ r g α 2 , with r g ≡ GM c 2 , and  \nα = GMµ ℏ c ≃ 0 . 1 ( M 15 M ⊙ ) ( µ 10 -12 eV ) . (1)  \nThe (scalar) cloud evolves according to a Schrodingerlike equation [47, 48], with eigenstates | a ⟩ ≡ | n a l a m a ⟩ , and ( n, l, m ) the principal, orbital and azimuthal angular momentum, 'quantum numbers'. (For vector clouds [38, 48, 49], we must include the total angular momentum.)  \nThe energy eigenvalues of the cloud scale as ϵ nlm = µ ( 1 -α 2 2 n 2 + f nl α 4 + h nl ˜ amα 5 ) , with ˜ a the dimensionless spin of the black hole, see [48]. At saturation, we have ˜ a ≃ α , whereas the combined system black hole plus cloud may still be rapidly rotating. One of the main difference w.r.t. ordinary atoms, however, is the presence of a decay/growing time, Γ -1 nlm ∝ µα 4 l +5 , for a given eigenstate [9, 47, 48, 50]. The (scalar) cloud may be populated by the dominant growing mode, | 211 ⟩ , or an excited state, | 322 ⟩ . Depending on α , they may be robust to GW emission (from the cloud itself) on astrophysical scales [9, 51-55]. They can also deplete through resonant transitions in binaries [37, 38], as we discuss here. In what follows we work with G = ℏ = c = 1 units.  \nGravitational collider goes eccentric. Following [37, 38] we consider a boson cloud around a black hole of mass M in a bound orbit with a companion object of mass M ⋆ , with q ≡ M ⋆ /M the mass ratio. The coordinates are centered at the black hole plus cloud system, with R ⋆ the radial distance to the perturber, and φ ⋆ the azimuthal angle. We consider planar motion with the spin parallel to the orbital angular momentum, with the orbit described by the semi-major axis a and the eccentricity e , while φ ⋆ corresponds to the true anomaly. We take the orbital frequency to be positive such that the two, co-rotating and the counter-rotating, orientations are identified by ˙ φ ⋆ = s | ˙ φ ⋆ | , with s = ± 1.  \nThe gravitational perturbations of the companion induce mixing of the atomic levels. For a perturber outside of the cloud R ⋆ ≫ r c the off-diagonal matrix elements of the Hamiltonian, ⟨ a | V ⋆ | b ⟩ , are given by a multipole expansion that can be written as an harmonic series [37, 38]  \n⟨ a | V ⋆ | b ⟩ l ⋆ = ∑ | m ⋆ |≤ l ⋆ η ( m ⋆ ) ab e -im ⋆ φ ⋆ , (2)  \n̸  \n̸  \nwith η ( m ⋆ ) ab ∝ R -( l ⋆ +1) ⋆ . The matrix elements obey selection rules which determine possible transitions, which we refer as hyperfine (only ∆ m = 0), fine (∆ ℓ = 0 , ∆ n = 0), and Bohr (∆ n = 0), respectively [37, 38].  \n̸  \nFor illustrative purposes, we consider a two-level model. The Hamiltonian equation is given by  \ni ( ˙ c a ˙ c b ) =    -∆ ϵ 2 η 0 ( R ⋆ R 0 ) -( l ⋆ +1) e i ∆ mφ ⋆ c . c . ∆ ϵ 2 -i Γ b    ( c a c b ) , (3)  \nwith ∆ m ≡ m b -m a , ∆ ϵ ≡ ϵ b -ϵ a the energy split, Γ b the width of the decaying mode, and η 0 the value of the perturbation at a reference point R 0 . Furthermore, since  \n(vacuum) GW emission is expected to reduce the initial eccentricity prior to encountering the resonant transition, and for the purpose of analytical understanding, in what follows we describe the orbital evolution in the Hamiltonian, H , of (3) using a small-eccentricity approximation, 1  \nφ ⋆ ≃ ϑ +2 e sin ϑ, R ⋆ ≃ a (1 -e cos ϑ ) , (4)  \n˙ ϑ ≡ s Ω , Ω = √ M (1 + q ) /a 3 , (5)  \nin terms of ϑ , the mean anomaly [56], and apply the Jacobi-Anger expansion into Bessel functions. Hence,  \nH = D + ∞ ∑ k = -∞   η k e i ( k +∆ m ) ϑ η k e -i ( k +∆ m ) ϑ   , (6)  \nD =   -∆ ϵ 2 ∆ ϵ 2 -i Γ b   , η k ∼ η 0 f 2 3 ( l ⋆ +1) e | k | | k | ! , f ≡ Ω Ω 0 ,  \nwhere we traded distance for orbital frequency. The case ( e, Γ b ) = 0 was studied in [38]. The slow GW-induced evolution of the orbital frequency, Ω( t ) ≃ Ω 0 + γ 0 t with γ 0 = 96 5 qM 5 / 3 Ω 11 / 3 0 (1+ q ) 1 / 3 , leads to a LZ transition [57, 58] between the energy levels. The transition is triggered for  \nΩ 0 = s ∆ ϵ ∆ m > 0 . (7)  \nThis value, dictated by the spectrum of the cloud, will serve as our reference point in the evolution of the binary.  \nIgnoring backreaction effects (see below), the LZ solution is controlled by the parameter z 0 ≡ η 2 0 / ( γ 0 | ∆ m | ), which determines the adiabaticity of the transition. As famously demonstrated in [57, 58], starting in the far past from the | a ⟩ state, in the limit 2 πz 0 ≫ 1, the eigenstate of the (time-dependent) Hamiltonian yields a complete population transfer into the (decaying) | b ⟩ mode in the far future. As it turns out, although the solution changes at finite time, controlled by the parameter v 0 ≡ Γ b / √ γ 0 | ∆ m | , the asymptotic properties of the system (ignoring backreaction) are remarkably robust against the value of the decaying width for the | b ⟩ state [59].  \nFor eccentric orbits, the evolution in (6) also features a transition at Ω 0 (for k = 0). However, it introduces a series of overtones  \nΩ k = f k Ω 0 , f k = ∆ m ∆ m + k , k ∈ Z . (8)  \nProvided each k -resonance is sufficiently narrow, we can ignore the other ( k ' = k ) terms in (6). As in [38], we can  \n̸  \nlinearize the orbital evolution near the transition, Ω( t ) = Ω k + f ( e ) γ k t , where f ( e ) = 1+ 73 e 2 24 + 37 e 4 96 (1 -e 2 ) 7 / 2 and γ k ≡ γ 0 f 11 / 3 k , such that the LZ solution now depends on the modified z k ≡ η 2 k f ( e ) γ k | ∆ m + k | and v k ≡ Γ b √ f ( e ) γ k | ∆ m + k | , respectively.  \nOrbital backreaction. Dissipative effects, such as GW emission, from the binary [38, 60] or the cloud itself [51-55], ionization [39-41, 60], and decay widths [59, 61, 62], strongly influence the LZ phenomenology, and vice versa. We focus here on the prevailing case of two-body GW emission, with the companion outside of the cloud, thus focusing on (hyper)fine resonances, combined with a two-level LZ transition into a decaying mode.  \nThe orbital dynamics is governed by flux-balance equations at infinity [38, 41, 44, 46, 62, 63], and at the black hole's horizon [8, 64, 65]:  \n˙ E o + ˙ E c + ˙ M = F GW ≡ -32 f ( e ) 5 M 5 q 2 ( q +1) a 5 , (9)  \n˙ L o + s ( ˙ L c + ˙ S ) = T GW ≡ F GW Ω g ( e ) f ( e ) , (10)  \n˙ M = 2Γ b E c( b ) , ˙ S = 2Γ b L c( b ) , (11)  \nwith g ( e ) = 1+ 7 e 2 8 (1 -e 2 ) 2 , and ˙ M , ˙ S the change of mass and spin due to the decay of the | b ⟩ state onto the black hole. 2 The orbital energy and angular momentum are given by E o = -M 2 q 2 a and L 2 o = ( M 5 q 3 )(1 -e 2 ) 2( q +1) | E o | , while for the cloud is a sum over the populated states, E c( i ) ≡ ( M c , 0 /µ ) ϵ i | c i | 2 , and similarly for L c( i ) with ϵ i → m i .  \nThe above equations can then be rewritten as  \nd Ω dt = rγ 0 f 11 / 3 f ( e ) , (12)  \nr ≡ ˙ E o F GW = 1 -b sgn( s ∆ m ) f -11 / 6 √ f ( e ) γ 0 | ∆ m + k | d | c a | 2 dt , (13)  \nde 2 dt = 2 3 f 8 / 3 γ 0 Ω 0 f ( e ) √ 1 -e 2 × (14) [ r ( f -√ 1 -e 2 ) -f + g ( e ) f ( e ) ] ,  \nin terms of the orbital parameters, where  \nb ≡ 3 M c , 0 M | ∆ m | f -3 / 2 | ∆ m + k | -1 / 2 (1 + q ) 1 / 3 αq ( M Ω 0 ) 1 / 3 Ω 0 √ γ 0 , (15)  \nparameterises the backreaction effects on the orbit due to the cloud. It is worth emphasising that the above equations apply to generic (planar) motion, regardless  \nof the value for the eccentricity. As we shall see, even for small initial conditions, the orbit is affected by large backreaction effects due to the presence of the cloud.  \nAs anticipated by the analysis in [38] for the case of circular orbits (which we encourage the reader to consult for further details), 'effective' LZ parameters emerge: ζ k ( t ) ≡ z k /r ( t ) and w k ( t ) ≡ v k / √ r ( t ), making it a fully nonlinear system. We can nonetheless estimate the value of the energy-momentum transfer near the resonance by self-consistently solving the condition ζ k = z k /r k ( ζ k ). For moderate-to-large population transfer ( ζ k ≳ 1), we find the limiting results:  \n(16)  \nr k ≃ 2 ( 1 + √ 1 -sgn( s ∆ m ) b k z k v k ) , ( w k ≫ ζ k )  \nr k ≃ ( 1 -sgn( s ∆ m ) b k 4 √ z k ) -1 , ( w k ≪ ζ k ) -1 .  \nAs discussed in [38], the orbital evolution branches into either floating ( r ≃ 0), for s ∆ m < 0, or sinking orbits ( r ≳ 1), for s ∆ m > 0. However, except for the trivial case when ζ k ≪ 1, due to the nonlinear nature of the problem the transfer of energy and angular momentum from the cloud to the orbit does not simply reduce to the quest for adiabaticity of the LZ transition, not even for w k ≪ ζ k . For instance, for extreme cases, with z k ≫ 1, the (unperturbed) transition spreads over long time scales, ∆ t LZ ≃ 4 √ z k /γ k [66], which in turn reduces the orbital impact, as we see in (16). As it turns out, in the large backreaction scenario, the sweet spot for floating orbits occurs when b k ≫ √ z k . Even though, due to the properties of the LZ solution, a strong decay width ( w k ≫ ζ k ) does not alter this picture, the impact on the orbit evolution as well as the population transfer becomes suppressed by 1 /v k , as shown in (16). On the other hand, for the sinking case, the largest values of r k are obtained for nonadiabatic transitions.  \nEccentric fixed point. For the GW-dominated epochs, with r ≃ 1, the leading order term in (14) vanishes, and the first contribution is at O ( e 2 ). Likewise for the k = 0 (main) resonance, for which the first term is ∝ ( r 2 -11 3 ) e 2 . As a result, the eccentricity is damped unless the orbit gets a large kick ( r ≳ 7 . 3). As the influence of the cloud increases, the RHS of (14) asymptotes (modulo a positive prefactor) to ( f k -1)( r -1), in which case it enters at leading order. Moreover, the differences in the GWfluxes in (9) and (10) generate a distinction between the early and late resonances. In the floating case, with r ≃ 0, the eccentricity grows for the early resonances ( f k < 1) and decays for the late ones ( f k > 1). This can be understood by noticing that, when ˙ E o ≃ 0, we have ˙ L o ∝ ( Ω -Ω 0 ΩΩ 0 + O ( e 2 ) ) , and using d ( L 2 o ) ∝ -d ( e 2 ) the eccentricity grows for Ω k < Ω 0 and decays whenever Ω k ≥ Ω 0 . This trend is reversed in the sinking case.  \nBecause of the changes in the evolution of the eccen-  \n<!-- image -->  \nFIG. 1. BBH eccentricities at f GW = 10 -2 Hz ( left ), evolved with a uniformly distributed q ∈ [0 . 1 , 1] and a boson cloud on the heavier black hole. The pale blue dots account for the values without a cloud [15]. (BBHs with e ≲ 10 -6 are not shown.) Cumulative effect, i.e. the ratio of binaries with eccentricities above a given value e 0 , ( right ), with (solid) and without (dotted) a cloud, both at 10 -2 Hz (blue) and 1 Hz (pink), respectively.  \n<!-- image -->  \ntricity across different resonances, it is instructive to look at the opposite limit e → 1. In that case, the RHS of (14) becomes ∝ r -1 (1 -e ) 3 . Let us consider the case of a floating orbit. Since the sign of de dt is positive for Ω k < Ω 0 , but turns negative when the eccentricity approaches e ≃ 1, this implies the existence of a critical 'attractor' fixed point, e cr , given by the condition g ( e cr ) /f ( e cr ) = f k [cf. (14)]. For instance,  \ne cr = { 0 . 46 , 0 . 35 , 0 . 29 } , for | ∆ m | = { 1 , 2 , 3 } , (17)  \nwith k = -1. Similarly, an unstable fixed point develops for the earlier and main sinking resonances.  \nFor the case of floating orbits (with s ∆ m < 0), if the backreaction is sufficiently effective to enforce r k ≃ 0 while the eccentricity approaches the critical point, one can then estimate the floating time ∆ t FL ≃ b k / √ γ k , leftover population | c a ( ∞ ) | 2 ≲ r k , and notably the growth of the eccentricity upon exiting the resonant transition,  \ne fin ≃ e cr √ 1 -e -C k , with C k ∼ √ γ k Ω k b k . (18)  \nAlthough we have used a small-eccentricity approximation to describe the initial stages of the cloud's evolution in (6), we have demonstrated through numerical studies that the behavior described above remains valid for generic (planar) orbits. See [46] and the appendices for details.  \nThe cloud's eccentric fossil. As it was argued in the literature [15-17], the distribution of masses and eccentricities observed with LISA can in principle distinguish between formation channels. However, the contrast between vacuum evolution and the large eccentricities produced by the cloud's resonant transition can lead to dramatic changes in the expected evolution of the system. As a proof of concept, we take the stellar-mass BBH population studied in [15], with chirp masses M ≲ 10 M ⊙ , expected to form exclusively in isolation, and with the spins aligned with the orbital angular momentum. As a  \nconsequence, the assumption of equatorial (uninclined) motion may be implemented without loss of generality (as done in [15]).  \nWe consider clouds of ultralight bosons of mass between 10 -13 and 10 -11 eV, surrounding black holes in co-rotating orbits. Superradiance may then excite the | 322 ⟩ state which, depending on the parent black hole's mass and birth orbital frequency, will experience a series of (hyper)fine transitions. 3 To illustrate the distinct physical effects, and following [14, 15], we consider a birth orbital frequency (for the cloud+BBH system) at Ω ini /π ≃ 10 -4 Hz, 4 and evolve, using the peak GW frequency [68]  \nf GW ≃ Ω π (1 + e ) 1 . 1954 (1 -e 2 ) 3 / 2 , (19)  \nuntil f GW = 10 -2 and 1 Hz. The final distribution is shown on the left panel of Fig. 1. While some of the BBHs experience an early overtone of the hyperfine transition, the majority are affected by the fine overtones instead. The BBHs then float over a period of time while increasing the orbital eccentricity. Moreover, the cloud typically either terminates there or decays later at the k = 0 resonance. Depending on the parameters, the ultimate decay may decrease the eccentricity or have a small impact on the orbit. As a result, a wedge-type distribution emerges, with the heavier black holes (within each wedge) subject to the largest increase in eccentricities.  \nThe cumulative effect is shown on the right panel of  \nFIG. 3. Evolution of the peak frequency through a resonant transition (at t -t res = 0) in the LISA band of a GW170809like event, compared to the evolution without the cloud.  \n<!-- image -->  \nFIG. 2. Percentage of binaries with eccentricities above 0.01 at f GW = 10 -2 Hz for different values of µ .  \n<!-- image -->  \nFig. 1, where a significant fraction of the population (with different parent masses) is affected by the resonances, yielding values of the eccentricities at 1Hz that may be within reach of mid-band and Decihertz detectors. As the value of µ increases (decreases) the location of the wedge in the distribution moves toward lower (higher) masses. The dependence on the value of µ for this population of BBH (with M ≲ 10 M ⊙ ), reaching e ≳ 10 -2 at f GW = 10 -2 Hz, is shown in Fig. 2.  \nEccentric in band. Because of the connection to formation channels, we discussed a sub-population of BBHs. However, similar conclusions apply to black holes with higher masses [46]. For instance, for a GW170809type event [69] ( M≃ 24 M ⊙ ), with a parent black hole M ≃ 20 M ⊙ (carrying the cloud) and a (heavier) companion M ⋆ ≃ 40 M ⊙ , we find f GW ≳ 10 -2 Hz at the k = -1 fine transition for µ ≃ 1 . 5 × 10 -12 eV. The BBH reaches the resonance and floats, with approximately constant orbital frequency for about six years, while the eccentricity increases from e ≲ 10 -2 to e ≃ 0 . 1, and likewise the peak GW frequency grows, while the cloud depletes. The resulting frequency evolution till merger, which is distinct from the growth of the eccentricity that may occur due to other astrophysical mechanisms [70, 71], 5 is displayed in Fig. 3. In addition to the notable features, higher harmonics would also become more relevant as the eccentricity increases [72]. As for the case of large tidal Love numbers [37, 42], 6 new dedicated templates will be needed to search for these phenomena in the GW data.  \nConclusions. We have shown that the presence of a boson cloud surrounding a black hole in a binary system can impact the distribution of masses and eccentricities observable with GW detectors. We have also  \nfound that a greater-than-expected value of the eccentricity, e ≳ 10 -2 at GW frequencies f GW ≃ 10 -2 Hz, develops for a (sub-)population of isolated stellar-mass BBHs (with M ≲ 10 M ⊙ ), right at the heart of the LISA band. Likewise, these BBHs will decay through GW emission to values of the eccentricities, i.e. e ≃ 10 -4 -10 -3 , within experimental reach of mid-band [6] and Decihertz [7] detectors. The observation of such GW signals would then provide tantalizing evidence for the existence of an ultralight particle of mass between 5 × 10 -13 and 2 . 5 × 10 -12 eV in nature. Furthermore, we have also shown that inband resonance transitions are possible, yielding dramatic changes in the GW frequency evolution, constituting yet another smoking-gun signature of the imprint of a boson cloud in the BBH dynamics.  \nThere are several venues for further exploration. First, unlike Bohr-type resonances, we have concentrated here on (co-rotating) hyperfine and fine transitions which occur outside of the cloud (for the range of parameters we considered), and therefore are not subject to ionization/dynamical friction [39, 41]. Preliminary studies suggest that a similar increase of eccentricity occurs for certain type of Bohr transitions at higher frequencies, which would put them within reach of the ET and CE detectors [46], but a more in-depth study is needed to take all relevant effects into account. Second, although generic, we have considered the case of uninclined orbits. This is justified for the populations of BBHs we considered here (formed in isolation with spins parallel with the angular momentum). However, to encompass also dynamically-formed systems, we must add inclination and new (off-plane) transitions [41]. While our results remain unchanged for quasi-planar motion, we also expect similar conclusions to apply for inclined orbits. (In fact, as shown in [74, 75], resonant transitions tend to equatorialize the orbit.) Finally, identical results can be drawn also for neutron star/black hole binaries. For instance, those formed in isolation have a parent black hole with mass near M ≃ 7 M ⊙ , and likewise in binaries with negligible eccentricity at f GW ≃ 1Hz [13, 25, 31].  \nThe presence of a boson cloud would then also lead to larger-than-expected eccentricities, providing additional circumstantial evidence for a new ultralight particle in nature.  \nAcknowledgements. We thank G.M. Tomaselli for informative discussions and the authors of [74] for sharing a draft of their related work. The work of MB and MK is supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC 2121 'Quantum Universe' - 390833306. MB and RAP are supported in part by the ERC Consolidator Grant 'Precision Gravity: From the LHC to LISA' provided by the European Research Council (ERC) under the European Union's H2020 research and innovation programme (grant No. 817791).", 'Appendix A: Gravitational atom': "The range of ultralight masses that we probe yield a (very) high 'axion decay constant' ( f a ), generically suppressing self-interactions [76-79] and coupling to other species [80-84]. In addition, for the black hole masses we studied, with 5 M ⊙ ≲ M ≲ 50 M ⊙ , we arrive at small-tomoderate values of α ≲ 0 . 25, keeping relativistic corrections to the hydrogenic states small [55, 60, 85, 86]. In this regime, the eigenvalues are given by [37, 48]  \nϵ nlm = µ [ 1 -α 2 2 n 2 -( 1 8 n + 6 2 l +1 -2 n ) α 4 n 3 + 16 2 l (2 l +1)(2 l +2) ˜ amα 5 n 3 ] . (A1)  \nFurthermore, from Detweiler's approximation [47, 48] (with our sign convention)  \n-Γ nlm ≃ 2˜ r + C nl g lm α 4 l +5 ( m Ω H -ω nlm ) , (A2)  \nwith g lm ≡ ∏ l k =1 [ k 2 (1 -˜ a 2 ) + (˜ am -2 r + ω ) 2 ] , C nl ≡ 2 4 l +1 ( n + l )! n 2 l +4 ( n -l -1)! ( l ! (2 l )!(2 l +1)! ) 2 , ˜ r + = 1 + √ 1 -˜ a 2 , M Ω H = ˜ a/ (2˜ r + ); whereas for the cloud itself decaying into GWs, we use the (non-relativistic) approximation in [51].", 'Appendix B: Tidal interactions': "For equatorial orbits and resonances triggered away from the cloud, the tidal interaction is given by [37, 38]  \n⟨ a | V ⋆ | b ⟩ ≡ ∞ ∑ l ⋆ =2 ∑ | m ⋆ |≤ l ⋆ η ( ⋆ ) ab e -im ⋆ φ ⋆ , (B1)  \nη ( ⋆ ) ab = -qα r R -( l ⋆ +1) ⋆ 4 π 2 l +1 ∣ ∣ ∣ Y ∗ ( ⋆ ) ( π 2 , φ ⋆ )∣ ∣ ∣ I r I Ω  \nI r ≈ ∫ d rr 2 ˆ R b ˆ R a r l ⋆ (B2)  \nc ⋆ , ∞  \nI Ω ≡ ∫ d Ω Y ∗ a ( θ, ϕ ) Y ( ⋆ ) ( θ, ϕ ) Y b ( θ, ϕ ) , (B3)  \nwhere ( ⋆ ) ≡ ( l ⋆ , m ⋆ ), r ≡ r/r c , R ⋆ ≡ R ⋆ /r c , ˆ R c = r 3 / 2 c R c is the (dimensionless) hydrogenic radial wavefunction, Y lm is the spherical harmonic. We leave to [46] the discussion on resonances 'inside the cloud', including dipolemediated transitions [60, 86, 87].  \nThe Jacobi-Anger expansion applied to (B1) [using (4)],  \ne ± i ∆ m ( ϑ +2 e sin( ϑ )) = ∞ ∑ k = -∞ ( ± 1) k J k (2 e ∆ m )e i ( k ± ∆ m ) ϑ , (B4)  \ncan be applied to the off-diagonal terms of the Hamiltonian [cf. (3) in the main text]. Using the properties of the Bessel function 7 , the tidal perturbation [cf. (6) of the main text] becomes  \nη ( ⋆ ) ab,k = η ( ⋆ ) ab, 0 f 2 3 ( l ⋆ +1) (∆ me ) | k | | k ! | ( 1 + ( l ⋆ +1) k 2∆ m ) + O ( e | k | +1 )  \n( ⋆ ) ab, 0 = -qα r c ( r c a 0 ) l ⋆ +1 4 π 2 l ⋆ +1 ∣ ∣ ∣ Y ∗ ( ⋆ ) ( π 2 , ϑ ) ∣ ∣ ∣ I r I Ω ,  \n, η (B5)  \n̸  \nwhere a 0 = [ M (1+ q ) / Ω 2 0 ] 1 / 3 . This interaction is nonzero only if the selection rules are satisfied [37, 38]: -m a + m ⋆ + m b = 0, l b + l ⋆ + l a = 2 p , | l a -l b | ≤ l ⋆ ≤ l a + l b . Furthermore, for equatorial orbits, for even (odd) l ⋆ only the spherical harmonics even (odd) in m ⋆ = 0 are nonzero.", 'Appendix C: Atomic resonances': 'As the populated state has a maximal azimuthal number m max = n -1, at the (hyper)fine resonances, it can only transition into states with lower m . Such transitions are only possible on co-rotating orbits, where they obey  \nFIG. 4. Orbital resonance frequencies of possible (fine) transitions for the | 322 ⟩ -component of the cloud. The blue and red arrows represent vacuum and floating BBH evolution, respectively; whereas the grey arrows point to various (rapidly decaying) modes that ultimately deplete the cloud. The thinner the line the smaller the cloud gets.  \ns ∆ m < 0, yielding floating-type motion 8 . From the selection rules for the tidal interactions, the | 211 ⟩ state has only one hyperfine transition to | 21 -1 ⟩ , while the only possible fine transition, to | 200 ⟩ , can only occur inside the cloud. Furthermore, for small values of α , we find that the floating time of the hyperfine transition would take longer than a Hubble time, preventing them to reach the LISA band. At the same time, (barring a precise fine tuning of the birth frequency of the BBH+cloud) for large values of α we expect the | 211 ⟩ component of the cloud to decay through its own GW emission before reaching the resonant transition (see also [62, 88]). On the other hand, the (longer-lived) | 322 ⟩ state may experience various types of resonances. In contrast to early and late resonances with k = ± 1, all of the k = 0 hyperfine transitions to | 32 m ⟩ happen at the same frequency (∆ m drops out of the ratio). The dominant main (hyperfine) transition is the one to | 320 ⟩ , with l ⋆ = 2 [as the | 322 ⟩ → | 32 -2 ⟩ resonance can only be mediated by the hexadecapole ( l ⋆ = 4), making it extremely weak and nonadiabatic]. Transitions to the | 32 ± 1 ⟩ states are not possible for equatorial orbits. The fine resonances from the excited state are (octopolar) to the | 31 -1 ⟩ , | 311 ⟩ and (quadrupolar) | 300 ⟩ states, in that order in the frequency domain. We show a succession of transitions in Fig. 4. We ignored the (sinking) highl ⋆ Bohr resonances that can overlap with the range of frequencies that we consider here, since they do not significantly affect the dynamical evolution. We postpone the general analysis of Bohr transitions to [46].', 'Linear solution': 'The solution of the (linear) LZ transition with ˙ Ω k ≃ f ( e ) γ k ≃ const, including a decaying width, is given  \nby [59, 61, 62] (see also [38, 58])  \n| c a | 2 = exp ( -v k τ -π 2 z k )∣ ∣ ∣ D iz k ( e i 3 π 4 ( τ -iv k ) )∣ ∣ ∣ 2 (D1) | c b | 2 = exp ( -v k τ -π 2 z k ) z k ∣ ∣ ∣ D iz k -1 ( e i 3 π 4 ( τ -iv k ) )∣ ∣ ∣ 2 ,  \nwhere D are parabolic cylinder functions [58], and we introduced a dimensionless time  \nτ ≡ t √ | ∆ m + k | f ( e ) γ k . (D2)  \nNotice that the physical impact of γ k (associated with the timescale of GW radiation) on the LZ dynamics is reflected via the ratio with the other two relevant parameters, for instance, the ratio with the (time) scale associated with the perturbation, ( η/ √ γ k ) 2 , measures the adiabaticity of the transition, while Γ b / √ γ k describes the impact of of the decaying width relative to the dynamical time. Remarkably, ignoring backreaction effects, even though the solution changes with respect to the nondecaying case, the transition probability at infinity, given by | c a ( ∞ ) | 2 = e -2 πz k , turns out to be independent of v k . 9  \nThe presence of a decay width, however, tends to smooth the LZ transition, by transferring a fraction of the initial state earlier than the case with v k = 0, and also by damping late-time oscillations. Consequently, for large widths, d dτ | c a | 2 peaks earlier than the v k = 0 case. In the limit v k ≫ z k , the LZ solution acquires a simple form [61, 62]  \n| c a | 2 = exp [ -2 z k ( arctan ( τ v k ) + π 2 )] , (D3)  \n| c b | 2 = | c a | 2 z k τ 2 + v 2 k . (D4)  \nIn this regime, d dτ | c a | 2 peaks at τ max ≃ -z k v k , and the width of the transition roughly scales as  \n∆Ω √ | ∆ m + k | f ( e ) γ ≃ 2 v k √ 1 + 2 z 2 k . (D5)  \nIn the two-level system with a decaying mode, we can relate the total energy and the angular momentum of the | a, b, BH ⟩ state as [cf. (9)-(11) in the main text]  \n( ˙ L c + ˙ S ) = ( ˙ E c + ˙ M ) ∆ m ∆ ϵ × ρ ϵ ( t ) ρ m ( t ) (D6)  \nρ x ≡ x a ∆ x d dt | c a | 2 -Γ a | c a | 2 -∑ ˆ Γ a | c a | 2 +( a → b ) , (D7)  \nwhere x ≡ { ϵ, m } , and the ˆ Γ represents other sources of dissipation, such as GW emission from the cloud [9, 5155] or ionization [39-41]. From the Schrodinger equation [cf. (3)], we find  \nρ x = d dt | c a | 2 +2 | c a | 2 Γ a +2 x b ∆ x ∑ ( ˆ Γ a | c a | 2 + ˆ Γ b | c b | 2 ) ,  \nand a similar equation applies with a ↔ b . After superradiance saturates the growth of the cloud in the | a ⟩ state, we have Γ a ≃ 0. Moreover, for the floating time of the resonances we consider here, we can ignore other sources of dissipation during the LZ transition (setting ˆ Γ ≃ 0) [46]. Hence, ρ x = ( d/dt ) | c a | 2 , which implies ρ ϵ /ρ m = 1, yielding the expression in (12)-(14) of the main text.', 'Nonlinear backreaction effects': "Because of the backreaction on the orbital frequency, which in turns controls the LZ transition, the problem becomes nonlinear, and it depends on the level-occupancy through the derivative of the parent state occupancy [cf. (13) in the main text]. We plot the this derivative for the linear LZ problem in Fig. 5. In the parts of the parameter space where it has compact support, the backreaction simply renormalizes the parameters ( z k , v k ) → ( ζ k , w k ) [cf. (14) of the main text] near the maximum value, as  \nζ k = η 2 k | ∆ m + k | ˙ Ω k = η 2 k | ∆ m + k | γ k f ( e ) r = z k r k ( ζ k ) , (D8)  \nand with the energy transfer itself depending on the leveloccupancy of the cloud. The solution can, nonetheless, be found self-consistently in terms of the relevant parameters.  \nThis mapping provides a useful indicator, even in the regime of strong backreaction, to evaluate whether, e.g., floating can occur, by indicating the breakdown of the linearization of the full problem. From the self-consistency condition (D8) the 'true' time scale for the nonlinear problem follows T = √ r k τ . The LZ function derivatives  \nFIG. 5. (Negative) derivative of the parent state occupancy, evaluated for a linear LZ transition. Brackets in the legend correspond to ( z k , v k ). Note the decrease, and shift to the left, of the amplitude, as well as the widening of the function and damping of the late-time oscillations, as the ratio v k /z k increases.  \n<!-- image -->  \ncan now be expressed in terms these quantities,  \nd | c i | 2 dτ = √ z k ζ k d | c i | 2 d T , (D9)  \nleading to  \n-d | c a | 2 dτ ∣ ∣ ∣ τ = -ζ k w k = √ z k ζ k ψ k ( ζ k , w k ) , (D10) ψ k ≡ 1 2 e -1 2 ( πζ k ) [ 2 w k ∣ ∣ ∣ ∣ C H iζ ( C w k ∆ w k 2 )∣ ∣ ∣ ∣ 2 -( C H 1 -iζ ( C ∗ ∆ w ∗ k 2 ) H iζ ( C ∆ w k 2 ) +c . c . ) ] ,  \nwhere C = 1 + i , ∆ w k = w k (1 -iζ k ) and H is the Hermite function. In what follows, for the sake of notation brevity, we take the small-eccentricity approximation, with f ( e ) ≃ 1. This is also justified by the fact that the eccentricity before the LZ transition is typically small across the parameter space, which we use to evaluate the type of transition the cloud will experience.", 'Negligible-decay regime': 'The ψ k function simplifies in the asymptotic regime  \nlim ζ → 0 ψ k ∼ √ πζ , & lim ζ →∞ ψ k ∼ 1 4 √ ζ ( w k ≪ ζ k ) ,  \nin which  \nζ k ∼ z k ( 1 -sgn( s ∆ m ) b k 4 √ z k ) , (D11)  \nand the value of r k given in (16) of the main text. As advertised, the energy transfer for the adiabatic floating  \n( s ∆ m< 0), is extremised by large values of b k with moderate z k parameters. In contrast, sinking ( s ∆ m > 0) is consistent with large adiabaticity only for moderate backreaction b k < 4 √ z k . Furthermore, the large backreaction limit ( b k ≫ 1) of the weakly-adiabatic regime  \nζ k ∼ sgn( s ∆ m ) ( z k b 2 k π ) 1 / 3 , r k ∼ sgn( s ∆ m )( πb 2 k z 2 k ) 1 / 3 ,  \nis only possible for sinking orbits. In such scenarios, the strong backreaction then mostly leads to nonadiabatic sinking transitions, with a potentially significant impact on the orbit', 'Strong-decay regime': '̸  \nIn general, the impact of the decay width depends on the dynamical timescale of the LZ transition and should be compared with the strength of the coupling. Moreover, for eccentric orbits, even if z 0 > v 0 at f 0 = 1, this hierarchy may be reversed for f k = 1, due to the e 2 | k | suppression in z k . Similarly to the case of negligible decay, we can also obtain approximate relations in the limit described by (D3), yielding  \nψ k ∼ 2 ζ k e -ζ k [ π -2 arctan( ζ k )] w k ( ζ 2 k +1) , (D12)  \nr k ∼ 1 ± 2 b k z k e -ζ k [ π -2 arctan( ζ k )] v k ( ζ 2 k +1) . (D13)  \nIn general, this regime shares various qualitative behavior as in the w k ≪ ζ k case, but with the adiabaticity gain/loss and orbital impact suppressed by the ratio b k /v k . In particular, large adiabaticity is consistent only for floating orbits, where we have  \nζ k ≃ √ 2 e √ b k z k v k , r k = z k ζ k → 0 . ( b k ≫ 1)  \nIn contrast, for sinking orbits, the strong backreaction requires a small population transfer, and we find  \nζ k ≃ v k 2 b k , r k ≃ 2 b k z k v k ≫ 1 . ( b k ≫ 1)  \nNotice that, somewhat counter-intuitively, a strongdecay width not only does not necessarily imply the total depletion of the cloud, instead it suppresses r k , hence the ability of the system to float, resulting in a lesser amount of the cloud being transferred to the decaying mode (see (E9) below).', 'Appendix E: Floating': "For values of the energy transfer r k ≲ 0 . 2, the growth of the orbital frequency is sufficiently suppressed to allow  \nfor the possibility of floating. In that case, the growth of (initially small) eccentricity is given, in units of the dynamical time introduced in (D5), by  \n̸  \ne ( t ) ≃ √ e 2 in + I e τ (1 -f k ) , ( f k = 1)  \n(E1) e ( t ) ≃ e in exp { -11 6 I e τ } , ( f k = 1) (E2) I e ≡ 2 3 √ γ 0 Ω 0 f 5 / 6 k | ∆ m + k | 1 / 2 .  \nThe critical points of the evolution of the eccentricity depend both on | ∆ m | and | k | . For the first few values of | k | they are described by the polynomial fit e cr = g 0 + g 1 k + g 2 k 2 , where g 0 = { 0 . 3 , 0 . 2 , 0 . 16 } , g 1 = { 0 . 18 , 0 . 16 , 0 . 14 } , g 2 = { 2 , 1 . 6 , 1 . 3 } × 10 -2 , all for | ∆ m | = { 1 , 2 , 3 } in respective order.  \nTo calculate the evolution of the eccentricity towards the fixed point, we change the time variable in the evolution equations [cf. (14)], from t to e 2 ( t ), yielding  \n∫ e 2 fin e 2 in d ( e 2 ) √ 1 -e 2 [ g ( e ) f ( e ) -f k ] = I e b k ( 1 -| c a ( ∞ ) | 2 ) , (E3)  \nThe result can then be approximated by (18) in the main text, where  \nC k = c ∆ m,k [ I e b k ( 1 -| c a ( ∞ ) | 2 ) + e 2 in 1 -f k ] , (E4)  \nand c ∆ m, -1 = { 2 . 37 , 2 , 1 . 37 } for ∆ m = {-1 , -2 , -3 } . From the value of e 2 fin in (E3), we can also estimate the duration of the floating period  \n∆ τ Fl = 1 I e ∫ e 2 fin e 2 in d ( e 2 ) 1 √ f ( e ) √ 1 -e 2 [ g ( e ) f ( e ) -f k ] , ≃ C k ˜ c ∆ m,k I e , c ∆ m,k ˜ c ∆ m,k = 1 √ f ( e cr ) . (E5)  \nAn exemplary parameter space of final eccentricity and floating time is shown in Fig. 6.  \nStrong floating provides a distinct phase of the nonlinear LZ transition, during which typically most of the population transfer occurs. From (13) of the main text, we have  \nd | c a | 2 dτ = r ( τ ) -1 b k , (E6)  \nyielding a linear-in-time decay of the population during floating, and a transfer of population given by  \n| c a ( ∞ ) | 2 ≃ | c a ( τ < Fl) | 2 -(1 -r k ) ∆ τ Fl b k . (E7)  \nIn the strong decay regime, the condition r ( t ) → 0 at  \n<!-- image -->  \nFIG. 6. Final eccentricity ( left ) and floating time ( right ), for the | 322 ⟩ → | 31 -1 ⟩ transition with k = -1, assuming initial conditions that lead to a (long-lasting) floating orbit ( r -1 ≃ 0) and e in ≪ 1.  \n<!-- image -->  \nthe resonance is necessary at each point in time, but not sufficient to guarantee a steady floating-type period. In addition, there must be enough of the cloud left to sustain a small r ( t ). Following [62], one can estimate the sufficient condition by considering the minimum amount of cloud needed to 'startjump' a floating period. Applying the linear LZ solution (D3) in (13) of the main text, we find  \n| c a | 2 min ≃ v k √ f ( e in ) 2 b k z k ( e in ) (1 -r k ) . (E8)  \nThe left hand side can be interpreted as the minimal amount of cloud needed to start floating at a particular resonance. In turn, if the right hand side is larger than one, floating cannot start. The same condition can also be used to estimate the amount of cloud left when floating stops, by matching into the linear LZ solution backwards , from the end of the floating time,  \n| c a ( ∞ ) | 2 ≃ v k √ f ( e fin ) 2 b k z k ( e fin ) (1 -r k ) , (E9)  \nwhich then becomes the portion of the cloud surviving after floating stops. Strong decay and small z k could in principle interrupt the floating period and leave a moderate amount of the cloud intact. 10 However, at the overtones we have z k ∼ e | 2 k | , which is increasing during floating. Hence, as the eccentricity approaches the critical point, for instance at the k = -1 overtone, the value of z k increases by a factor ( e cr /e in ) 2 ≃ 10 2 -10 3 , thus significantly extending the floating period, and reducing the amount of cloud left after the transition, in comparison with the na¨ıve estimate in (E8).  \nIn general, the equations in (E3), (E5) and (E9) must be solved self-consistently in order to determine the end state of floating. As an estimate, we may apply (E1) to (E9), and assuming e in ≪ 1, we have  \n∆ τ Fl b k / √ f ( e cr ) ≃ ( x -1 2 x + √ ( x +1) 2 4 x 2 -λ x ) , (E10) λ = v k 2 b k z k ( e in ) , x = I e b k (1 -f k ) e 2 in .  \nNotice that for the k = -1 overtone, the dependence on e in drops out from the ratio λ/x . In this case λ ≪ x , and we find the longest periods of floating, ∆ τ Fl ≃ b k / √ f ( e cr ), and largest depletion of the | c a ( ∞ ) | 2 ≲ r k . Depending on the resonances and the parameter space, such hierarchy may be also valid for higher overtones.", 'Appendix F: Numerical validation': "Our semi-analytic approach, described in Apps. D and E, rely in part on the small initial eccentricity of the BBHs prior to entering the transition region. We have validated this by numerically solving the Schrodinger equation for arbitrary (planar) orbits, coupled with the energy-momentum balance equations [cf. (3),(9)-(11) of the main text], for a number of representative examples, including a broad inital eccentricity range 10 -2 ≲ e in ≲ 0 . 5. For this purpose, we used the generic description of the orbital evolution, R ⋆ = a (1 -e cos u ), as a function of the eccentric anomaly u , which is related to the mean anomaly via Kepler's equation u -e sin( u ) = ϑ (e.g. [56]). We used the NDSolve routine in Mathematica [89], monitoring the violation of unitarity in the Γ b = 0 regime, residual for Γ b ≥ 0, as well as the state occupation numbers (D1). We have also checked the consistency of our results by comparing with the values in [48] for ( e, Γ b ) = 0.  \nIn general, the numerical results are broadly consis-  \n<!-- image -->  \nFIG. 7. Numerical evolution of the examples described in App. F. The figures correspond to the same z k but different values for b k , with the left plot having a stronger backreaction parameter (a few times larger) than on the right. The curves show the eccentricity evolution for vanishing (pink) and large decay v k ≫ z k (blue, dot-dashed), both yielding floating orbits ( r k ≃ 0), except for the case of a large decay with small backreaction (blue curve on the right panel has r k ≃ 0 . 8). The semi-analytic solution (black, dashed) [cf. (14) in the main text], is in remarkable agreement within its regime of validity ( r k ≃ 0).  \n<!-- image -->  \ntent with the analytical arguments. As expected, the largest discrepancy occurs for the estimates of r ( t ) and | c a ( ∞ ) | 2 . For the floating regime, the analytic results tend to overestimate their respective values by a factor of few to an order of magnitude. Hence, our conclusions based on the analytic approximations err on the side of caution . In Fig. 7 we present (left) the numerical eccentricity evolution of a two-level system with ∆ m = 2 for values { z 0 = 4 , v 0 = 2 , b 0 = 200 , e 0 = 10 -2 } , at the main resonance. For these parameters, we have { z -1 = 7 · 10 -3 , v -1 = 3 . 4 , b -1 = 450 , I e = 2 · 10 -4 } at the k = -1 overtone, for which one finds a strongdecay regime: v -1 /z -1 = 512. We start the numerical evaluation before the first overtone, following also the case v -1 = 0. Away from the resonance, the orbital evolution closely follows standard GW evolution in vacuum. The approximations from Sec. D then correctly predict the strong floating that occurs (in both setups) at the frequency given by (8). Notice that, modulo small deviations, both cases follow the r k ≃ 0 prediction for the eccentricity growth in (E3). Furthermore, broadly consistent with our estimates, the left-over occupancy of the inital state is given by | c a ( ∞ ) | 2 ≃ 10 -4 and | c a ( ∞ ) | 2 ≃ 10 -3 , with and without the decay, respectively. Finally, reducing the backreaction to b 0 = 25 (Fig. 7, right), the decaying case does not develop floating, as r k ≃ 0 . 8, although as expected the eccentricity will still grow slightly and the transfer of population is somewhat increased compared to the linear LZ transition. In contrast, the non-decaying case still exhibits a strong floating, follows the predicted growth of the eccentricity, and transfers the parent state up to | c a ( ∞ ) | 2 ≃ 5 × 10 -3 .", 'Appendix G: Eccentricity distribution': "To evaluate the distribution of eccentricities in the population from [15], we order (in frequency) all possible hyperfine and fine resonances with the correspond-  \ning overtones, up to | k | = 5. (A subset is displayed in Fig. 4.) We then calculate the corresponding parameters ( z k , v k , b k ) for every resonance the cloud may encounter. The condition that the initial state had enough time to reach the resonance is imposed, constrained by the lifetime imposed by its own GW emission. We evolve in vacuum via [43, 44] between transitions. At every resonance, we estimate r k as explained in App. D, and from there the floating time, eccentricity growth, and left-over cloud, as explained in App. E. We use a cutoff at 10 -3 of the initial M c, 0 /M , to estimate when the cloud's influence on the orbital dynamics becomes negligible.  \nFor the BBHs shown in Fig. 1 of the main letter, only a few percent experience a floating transition | 322 ⟩ → | 31 -1 ⟩ with k < -1, while 25 % float at the k = -1 resonance, making it the dominant one. Approximately 15 % of the population float at the k = 0 transition. Only a few experience resonances to | 311 ⟩ or | 300 ⟩ states. The black holes with the largest masses, hence largest values of α for fixed µ , can float at the | 320 ⟩ hyperfine transition, where of the order of 5 % (each), see an earlier resonance and increase eccentricity there; or see the main or k = 1 resonance and decrease eccentricity instead. We do not find a significant impact for overtones with | k | > 4. We show in the left panel of Fig. 8 the distribution of different k 's experienced by the BBH population we studied in the main text. We show in the right panel the distribution of initial eccentricities before floating starts. As a consistency check, notice that for most of the resonances we have e in ≲ 0 . 1, which further support the validity of the small-eccentricity approximation. In Fig. 9 we plot the eccentricity evolution for three representative cases drawn from the population in Fig. 1 of the main text.  \n<!-- image -->  \nFIG. 8. Number of floating resonances for specific k ( left ) and binned by initial eccentricity ( right ) for the floating orbits the population shown in Fig. 1 experiences.  \n<!-- image -->  \nFIG. 9. Eccentricity evolution (in terms of the GW frequency) for representative examples from Fig. 1 in the main text. Solid lines show the evolution with the cloud, while dashed lines represents the standard vacuum evolution.  \n<!-- image -->", 'Birth frequency': "Similar to Fig. 1 in the main letter, we show in Fig. 10 the eccentricity distribution at f GW = 10 -2 Hz, but with a lower birth frequency for the BBH+cloud system, at Ω ini /π = 5 × 10 -5 Hz. Although the number of larger-than-expected eccentricity points is somewhat lower, the plot illustrates the robustness of our predictions against changes in the initial conditions. In general, a large portion of BBHs that merge within a Hubble time will be born with orbital frequencies between [0 . 5 , 1] × 10 -4 Hz [14]. Moreover, BBHs that are formed through a common-envelope mechanism that can bring the binary closer, e.g. [67], at the same time may produce a (younger) secondary black hole that can ultimately carry the boson cloud, pushing the cloud's birth frequency towards higher values, thus avoiding altogether ( k = 0) hyperfine transitions that can decrease the eccentricity. Overall, our findings demonstrate that a skewedtype distribution and larger-than-expected eccentricities are robust predictions of the presence of boson clouds in the BBH dynamics.", 'Appendix H: Background (black hole) parameters': "Finally, let us comment on the backreaction on the black hole mass and spin. We will implement this in detail in a forthcoming paper. However, it is worth emphasising already that this effect does not change the main message of the letter: fixed points in the eccentricity evolution and the imprint in its final distribution at observable GW frequencies. We provide here a brief argument to support this claim.  \nFirst of all, the changes in α (through the growth of the black hole mass) move the position of the resonances towards higher frequencies, which may help in sustaining the floating condition. Regarding the changes in spin, while they depend on the sign of m b , most of transitions tend to decrease the value of the parent black hole's spin, which in turn affects the superradiant condition. As a result, the | 322 ⟩ state may develop a decay width. However, even in the most extreme scenario, and assuming the entire cloud is instantly reabsorbed, we would still have Γ | 322 ⟩ ≃ α 13 (with ℓ = 2 for the excited state). As a consequence, as we show in Fig. 11, the decay time turns out to be longer than the floating time for all of the fine transitions. This effectively allows us to ignore the resulting width of the | 322 ⟩ state in the determination of the growth of eccentricity.  \nFor the hyperfine transitions, on the other hand, the situation is more subtle, since the floating times are in principle longer. In that case, one has to take into account that the decrease of the spin does not occur instantaneously, but on the time-scales of the LZ transition governing the population of the | b ⟩ state. Preliminary studies indicate that the qualitative results for the hyperfine resonance are not significantly modified. In particular, since hyperfine transitions fall typically in the weak-decay regime, the floating condition is expected to be generally robust. This expectation is also in agreement with the results of [62], where the hyperfine transition | 211 ⟩ → | 21 -1 ⟩ was considered including the evolution of the background parameters. The authors did not find a significant disruption of the floating mechanism  \n<!-- image -->  \nFIG. 10. Same as Fig. 1 in the main letter, but with Ω ini /π = 5 × 10 -5 Hz as the birth orbital frequency for the BBH+cloud system.  \n<!-- image -->  \nFIG. 11. Ratio of floating and (excited-state) decay time for the k = -1 overtone of the fine transition, assuming the parent black hole's spin is reduced to ˜ a = 0 . 9˜ a s , with ˜ a s = 2 α the value at saturation.  \n<!-- image -->  \n- [1] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), Open Data from the Third Observing Run of LIGO, Virgo, KAGRA, and GEO, Astrophys. J. Suppl. 267 , 29 (2023), arXiv:2302.03676 [gr-qc].\n- [2] R. Alves Batista et al. , EuCAPT White Paper: Opportunities and Challenges for Theoretical Astroparticle Physics in the Next Decade, (2021), arXiv:2110.10074 [astro-ph.HE].\n- [3] K. G. Arun et al. (LISA), New horizons for fundamental physics with LISA, Living Rev. Rel. 25 , 4 (2022), arXiv:2205.01597 [gr-qc].\n- [4] M. Maggiore et al. , Science Case for the Einstein Telescope, JCAP 03 , 050, arXiv:1912.02622 [astro-ph.CO].\n- [5] D. Reitze et al. , Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO, Bull. Am. Astron. Soc. 51 , 035 (2019), arXiv:1907.04833 [astro-ph.IM].\n- [6] S. Baum, Z. Bogorad, and P. W. Graham, Gravitational Wave Science in the Mid-Band with Atom Interferometers, (2023), arXiv:2309.07952 [gr-qc].\n- [7] S. Kawamura et al. , Current status of space gravitational wave antenna DECIGO and B-DECIGO, PTEP 2021 , 05A105 (2021), arXiv:2006.13545 [gr-qc].\n- [8] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper,  \ndue to the absorption of mass and angular momentum. Notice, in addition, that for the excited state we have Γ | 322 ⟩ ≪ Γ | 211 ⟩ , which further supports this conclusion. We postpone a full analysis to [46].  \n- and J. March-Russell, String Axiverse, Phys. Rev. D 81 , 123530 (2010), arXiv:0905.4720 [hep-th].\n- [9] A. Arvanitaki and S. Dubovsky, Exploring the String Axiverse with Precision Black Hole Physics, Phys. Rev. D 83 , 044026 (2011), arXiv:1004.3558 [hep-th].\n- [10] D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643 , 1 (2016), arXiv:1510.07633 [astro-ph.CO].\n- [11] M. Demirtas, C. Long, L. McAllister, and M. Stillman, The Kreuzer-Skarke Axiverse, JHEP 04 , 138, arXiv:1808.01282 [hep-th].\n- [12] V. M. Mehta, M. Demirtas, C. Long, D. J. E. Marsh, L. McAllister, and M. J. Stott, Superradiance in string theory, JCAP 07 , 033, arXiv:2103.06812 [hep-th].\n- [13] K. Belczynski, V. Kalogera, and T. Bulik, A Comprehensive study of binary compact objects as gravitational wave sources: Evolutionary channels, rates, and physical properties, Astrophys. J. 572 , 407 (2001), arXiv:astroph/0111452.\n- [14] I. Kowalska, T. Bulik, K. Belczynski, M. Dominik, and D. Gondek-Rosinska, The eccentricity distribution of compact binaries, Astronomy &amp; Astrophysics 527 , A70 (2011).\n- [15] K. Breivik, C. L. Rodriguez, S. L. Larson, V. Kalogera, and F. A. Rasio, Distinguishing Between Formation  \n- Channels for Binary Black Holes with LISA, Astrophys. J. Lett. 830 , L18 (2016), arXiv:1606.09558 [astroph.GA].\n- [16] A. Nishizawa, E. Berti, A. Klein, and A. Sesana, eLISA eccentricity measurements as tracers of binary black hole formation, Phys. Rev. D 94 , 064020 (2016), arXiv:1605.01341 [gr-qc].\n- [17] A. Nishizawa, A. Sesana, E. Berti, and A. Klein, Constraining stellar binary black hole formation scenarios with eLISA eccentricity measurements, Mon. Not. Roy. Astron. Soc. 465 , 4375 (2017), arXiv:1606.09295 [astroph.HE].\n- [18] C. L. Rodriguez, M. Zevin, C. Pankow, V. Kalogera, and F. A. Rasio, Illuminating Black Hole Binary Formation Channels with Spins in Advanced LIGO, Astrophys. J. Lett. 832 , L2 (2016), arXiv:1609.05916 [astro-ph.HE].\n- [19] C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee, and F. A. Rasio, Post-Newtonian Dynamics in Dense Star Clusters: Highly-Eccentric, Highly-Spinning, and Repeated Binary Black Hole Mergers, Phys. Rev. Lett. 120 , 151101 (2018), arXiv:1712.04937 [astro-ph.HE].\n- [20] C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee, K. Kremer, F. A. Rasio, J. Samsing, C. S. Ye, and M. Zevin, Post-Newtonian Dynamics in Dense Star Clusters: Formation, Masses, and Merger Rates of HighlyEccentric Black Hole Binaries, Phys. Rev. D 98 , 123005 (2018), arXiv:1811.04926 [astro-ph.HE].\n- [21] M. E. Lower, E. Thrane, P. D. Lasky, and R. Smith, Measuring eccentricity in binary black hole inspirals with gravitational waves, Phys. Rev. D 98 , 083028 (2018), arXiv:1806.05350 [astro-ph.HE].\n- [22] L. Randall and Z.-Z. Xianyu, Eccentricity without Measuring Eccentricity: Discriminating among Stellar Mass Black Hole Binary Formation Channels, Astrophys. J. 914 , 75 (2021), arXiv:1907.02283 [astro-ph.HE].\n- [23] X. Fang, T. A. Thompson, and C. M. Hirata, The Population of Eccentric Binary Black Holes: Implications for mHz Gravitational Wave Experiments, Astrophys. J. 875 , 75 (2019), arXiv:1901.05092 [astro-ph.HE].\n- [24] I. M. Romero-Shaw, P. D. Lasky, E. Thrane, and J. C. Bustillo, GW190521: orbital eccentricity and signatures of dynamical formation in a binary black hole merger signal, Astrophys. J. Lett. 903 , L5 (2020), arXiv:2009.04771 [astro-ph.HE].\n- [25] M. A. Sedda, Dissecting the properties of neutron star - black hole mergers originating in dense star clusters, Commun. Phys. 3 , 43 (2020), arXiv:2003.02279 [astroph.GA].\n- [26] H. Glanz and H. B. Perets, Common envelope evolution of eccentric binaries, Monthly Notices of the Royal Astronomical Society 507 , 2659 (2021), arXiv:2105.02227 [astro-ph.SR].\n- [27] M. Zevin, I. M. Romero-Shaw, K. Kremer, E. Thrane, and P. D. Lasky, Implications of Eccentric Observations on Binary Black Hole Formation Channels, Astrophys. J. Lett. 921 , L43 (2021), arXiv:2106.09042 [astro-ph.HE].\n- [28] A. Gualandris, F. M. Khan, E. Bortolas, M. Bonetti, A. Sesana, P. Berczik, and K. Holley-Bockelmann, Eccentricity evolution of massive black hole binaries from formation to coalescence, Monthly Notices of the Royal Astronomical Society 511 , 4753-4765 (2022).\n- [29] M. Garg, S. Tiwari, A. Derdzinski, J. G. Baker, S. Marsat, and L. Mayer, The minimum measurable eccentricity from gravitational waves of LISA massive black  \n- hole binaries, (2023), arXiv:2307.13367 [astro-ph.GA].\n- [30] P. Saini, Resolving the eccentricity of stellar mass binary black holes with next generation ground-based gravitational wave detectors 10.1093/mnras/stae037 (2023), arXiv:2308.07565 [astro-ph.HE].\n- [31] R. Dhurkunde and A. H. Nitz, Search for eccentric NSBH and BNS mergers in the third observing run of Advanced LIGO and Virgo, (2023), arXiv:2311.00242 [astro-ph.HE].\n- [32] Y. B. Zel'Dovich, Generation of Waves by a Rotating Body, Soviet Journal of Experimental and Theoretical Physics Letters 14 , 180 (1971).\n- [33] Y. B. Zel'Dovich, Amplification of Cylindrical Electromagnetic Waves Reflected from a Rotating Body, Soviet Journal of Experimental and Theoretical Physics 35 , 1085 (1972).\n- [34] W. H. Press and S. A. Teukolsky, Floating Orbits, Superradiant Scattering and the Black-hole Bomb, Nature 238 , 211 (1972).\n- [35] W. E. East, Massive Boson Superradiant Instability of Black Holes: Nonlinear Growth, Saturation, and Gravitational Radiation, Phys. Rev. Lett. 121 , 131104 (2018), arXiv:1807.00043 [gr-qc].\n- [36] R. Brito, V. Cardoso, and P. Pani, Superradiance: New Frontiers in Black Hole Physics, Lect. Notes Phys. 906 , pp.1 (2015), arXiv:1501.06570 [gr-qc].\n- [37] D. Baumann, H. S. Chia, and R. A. Porto, Probing Ultralight Bosons with Binary Black Holes, Phys. Rev. D 99 , 044001 (2019), arXiv:1804.03208 [gr-qc].\n- [38] D. Baumann, H. S. Chia, R. A. Porto, and J. Stout, Gravitational Collider Physics, Phys. Rev. D 101 , 083019 (2020), arXiv:1912.04932 [gr-qc].\n- [39] D. Baumann, G. Bertone, J. Stout, and G. M. Tomaselli, Ionization of gravitational atoms, Phys. Rev. D 105 , 115036 (2022), arXiv:2112.14777 [gr-qc].\n- [40] D. Baumann, G. Bertone, J. Stout, and G. M. Tomaselli, Sharp Signals of Boson Clouds in Black Hole Binary Inspirals, Phys. Rev. Lett. 128 , 221102 (2022), arXiv:2206.01212 [gr-qc].\n- [41] G. M. Tomaselli, T. F. M. Spieksma, and G. Bertone, Dynamical friction in gravitational atoms, JCAP 07 , 070, arXiv:2305.15460 [gr-qc].\n- [42] H. S. Chia, T. D. P. Edwards, D. Wadekar, A. Zimmerman, S. Olsen, J. Roulet, T. Venumadhav, B. Zackay, and M. Zaldarriaga, In Pursuit of Love: First Templated Search for Compact Objects with Large Tidal Deformabilities in the LIGO-Virgo Data, (2023), arXiv:2306.00050 [gr-qc].\n- [43] P. C. Peters and J. Mathews, Gravitational radiation from point masses in a Keplerian orbit, Phys. Rev. 131 , 435 (1963).\n- [44] P. C. Peters, Gravitational Radiation and the Motion of Two Point Masses, Phys. Rev. 136 , B1224 (1964).\n- [45] E. Berti, R. Brito, C. F. B. Macedo, G. Raposo, and J. L. Rosa, Ultralight boson cloud depletion in binary systems, Phys. Rev. D 99 , 104039 (2019), arXiv:1904.03131 [grqc].\n- [46] M. Boskovic, M. Koschnitzke, and R. A. Porto, In preparation, .\n- [47] S. L. Detweiler, KLEIN-GORDON EQUATION AND ROTATING BLACK HOLES, Phys. Rev. D 22 , 2323 (1980).\n- [48] D. Baumann, H. S. Chia, J. Stout, and L. ter Haar, The Spectra of Gravitational Atoms, JCAP 12 , 006,\n- arXiv:1908.10370 [gr-qc].\n- [49] W. E. East, Superradiant instability of massive vector fields around spinning black holes in the relativistic regime, Phys. Rev. D 96 , 024004 (2017), arXiv:1705.01544 [gr-qc].\n- [50] S. R. Dolan, Instability of the massive Klein-Gordon field on the Kerr spacetime, Phys. Rev. D 76 , 084001 (2007), arXiv:0705.2880 [gr-qc].\n- [51] H. Yoshino and H. Kodama, Gravitational radiation from an axion cloud around a black hole: Superradiant phase, PTEP 2014 , 043E02 (2014), arXiv:1312.2326 [gr-qc].\n- [52] R. Brito, V. Cardoso, and P. Pani, Black holes as particle detectors: evolution of superradiant instabilities, Class. Quant. Grav. 32 , 134001 (2015), arXiv:1411.0686 [gr-qc].\n- [53] A. Arvanitaki, M. Baryakhtar, and X. Huang, Discovering the QCD Axion with Black Holes and Gravitational Waves, Phys. Rev. D 91 , 084011 (2015), arXiv:1411.2263 [hep-ph].\n- [54] R. Brito, S. Ghosh, E. Barausse, E. Berti, V. Cardoso, I. Dvorkin, A. Klein, and P. Pani, Gravitational wave searches for ultralight bosons with LIGO and LISA, Phys. Rev. D 96 , 064050 (2017), arXiv:1706.06311 [grqc].\n- [55] N. Siemonsen, T. May, and W. E. East, Modeling the black hole superradiance gravitational waveform, Phys. Rev. D 107 , 104003 (2023), arXiv:2211.03845 [gr-qc].\n- [56] S. Tremaine, Dynamics of Planetary Systems (2023).\n- [57] L. Landau, Zur theorie der energieubertragung. ii (1932).\n- [58] C. Zener, Nonadiabatic crossing of energy levels, Proc. Roy. Soc. Lond. A 137 , 696 (1932).\n- [59] V. M. Akulin and W. P. Schleich, Landau-zener transition to a decaying level, Phys. Rev. A 46 , 4110 (1992).\n- [60] R. Brito and S. Shah, Extreme mass-ratio inspirals into black holes surrounded by scalar clouds, Phys. Rev. D 108 , 084019 (2023), arXiv:2307.16093 [gr-qc].\n- [61] N. V. Vitanov and S. Stenholm, Pulsed excitation of a transition to a decaying level, Phys. Rev. A 55 , 2982 (1997).\n- [62] T. Takahashi, H. Omiya, and T. Tanaka, Evolution of binary systems accompanying axion clouds in extreme mass ratio inspirals, Phys. Rev. D 107 , 103020 (2023), arXiv:2301.13213 [gr-qc].\n- [63] T. Takahashi, H. Omiya, and T. Tanaka, Axion cloud evaporation during inspiral of black hole binaries: The effects of backreaction and radiation, PTEP 2022 , 043E01 (2022), arXiv:2112.05774 [gr-qc].\n- [64] G. Ficarra, P. Pani, and H. Witek, Impact of multiple modes on the black-hole superradiant instability, Phys. Rev. D 99 , 104019 (2019), arXiv:1812.02758 [gr-qc].\n- [65] L. Hui, Y. T. A. Law, L. Santoni, G. Sun, G. M. Tomaselli, and E. Trincherini, Black hole superradiance with dark matter accretion, Phys. Rev. D 107 , 104018 (2023), arXiv:2208.06408 [gr-qc].\n- [66] N. V. Vitanov, Transition times in the Landau-Zener model, Phys. Rev. A 59 , 988 (1999), arXiv:quantph/9811066.\n- [67] K. Belczynski et al. , Evolutionary roads leading to low effective spins, high black hole masses, and O1/O2 rates for LIGO/Virgo binary black holes, Astron. Astrophys. 636 , A104 (2020), arXiv:1706.07053 [astro-ph.HE].\n- [68] L. Wen, On the eccentricity distribution of coalescing black hole binaries driven by the Kozai mechanism in globular clusters, Astrophys. J. 598 , 419 (2003), arXiv:astro-ph/0211492.  \n- [69] B. P. Abbott et al. (LIGO Scientific, Virgo), GWTC1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9 , 031040 (2019), arXiv:1811.12907 [astro-ph.HE].\n- [70] L. Randall and Z.-Z. Xianyu, An Analytical Portrait of Binary Mergers in Hierarchical Triple Systems, Astrophys. J. 864 , 134 (2018), arXiv:1802.05718 [gr-qc].\n- [71] L. Randall and Z.-Z. Xianyu, Observing Eccentricity Oscillations of Binary Black Holes in LISA, (2019), arXiv:1902.08604 [astro-ph.HE].\n- [72] D. Wadekar, J. Roulet, T. Venumadhav, A. K. Mehta, B. Zackay, J. Mushkin, S. Olsen, and M. Zaldarriaga, New black hole mergers in the LIGO-Virgo O3 data from a gravitational wave search including higher-order harmonics, (2023), arXiv:2312.06631 [gr-qc].\n- [73] B. Su, Z.-Z. Xianyu, and X. Zhang, Probing Ultralight Bosons with Compact Eccentric Binaries, Astrophys. J. 923 , 114 (2021), arXiv:2107.13527 [gr-qc].\n- [74] G. M. Tomaselli, T. F. M. Spieksma, and G. Bertone, The resonant history of gravitational atoms in black hole binaries, (2024), arXiv:2403.03147 [gr-qc].\n- [75] G. M. Tomaselli, T. F. M. Spieksma, and G. Bertone, The legacy of boson clouds on black hole binaries, (2024), arXiv:2407.12908 [gr-qc].\n- [76] H. Yoshino and H. Kodama, The bosenova and axiverse, Class. Quant. Grav. 32 , 214001 (2015), arXiv:1505.00714 [gr-qc].\n- [77] A. Gruzinov, Black Hole Spindown by Light Bosons, (2016), arXiv:1604.06422 [astro-ph.HE].\n- [78] M. Baryakhtar, M. Galanis, R. Lasenby, and O. Simon, Black hole superradiance of self-interacting scalar fields, Phys. Rev. D 103 , 095019 (2021), arXiv:2011.11646 [hepph].\n- [79] H. S. Chia, C. Doorman, A. Wernersson, T. Hinderer, and S. Nissanke, Self-interacting gravitational atoms in the strong-gravity regime, JCAP 04 , 018, arXiv:2212.11948 [gr-qc].\n- [80] J. a. G. Rosa and T. W. Kephart, Stimulated Axion Decay in Superradiant Clouds around Primordial Black Holes, Phys. Rev. Lett. 120 , 231102 (2018), arXiv:1709.06581 [gr-qc].\n- [81] M. Boskovic, R. Brito, V. Cardoso, T. Ikeda, and H. Witek, Axionic instabilities and new black hole solutions, Phys. Rev. D 99 , 035006 (2019), arXiv:1811.04945 [gr-qc].\n- [82] H. Fukuda and K. Nakayama, Aspects of Nonlinear Effect on Black Hole Superradiance, JHEP 01 , 128, arXiv:1910.06308 [hep-ph].\n- [83] T. F. M. Spieksma, E. Cannizzaro, T. Ikeda, V. Cardoso, and Y. Chen, Superradiance: Axionic couplings and plasma effects, Phys. Rev. D 108 , 063013 (2023), arXiv:2306.16447 [gr-qc].\n- [84] Y. Chen, X. Xue, and V. Cardoso, Black Holes as Neutrino Factories, (2023), arXiv:2308.00741 [hep-ph].\n- [85] E. Cannizzaro, L. Sberna, S. R. Green, and S. Hollands, Relativistic Perturbation Theory for Black-Hole Boson Clouds, Phys. Rev. Lett. 132 , 051401 (2024), arXiv:2309.10021 [gr-qc].\n- [86] F. Duque, C. F. B. Macedo, R. Vicente, and V. Cardoso, Axion Weak Leaks: extreme mass-ratio inspirals in ultralight dark matter, (2023), arXiv:2312.06767 [gr-qc].\n- [87] S. L. Detweiler and E. Poisson, Low multipole contributions to the gravitational selfforce, Phys. Rev. D 69 ,  \n084019 (2004), arXiv:gr-qc/0312010.  \n- [88] Y. Cao and Y. Tang, Signatures of ultralight bosons in compact binary inspiral and outspiral, Phys. Rev. D 108 ,  \n- 123017 (2023), arXiv:2307.05181 [gr-qc].\n- [89] W. R. Inc., Mathematica, Version 13.2, champaign, IL, 2022."}