Previous works have argued that future gravitationalwave detectors will be able to probe the properties of astrophysical environments where binaries coalesce including accretion disks but also dark matter structures. Most analyses have resorted to a Newtonian modeling of the environmental effects which are not suited to study extrememassratio inspirals immersed in structures of ultralight bosons. In this Letter we use relativistic perturbation theory to consistently study these systems in spherical symmetry. We compute the flux of scalar particles and the rate at which orbital energy is dissipated via gravitational radiation and depletion of scalars i.e. dynamical friction. Our results confirm that the Laser Interferometer Space Antenna will be able to probe ultralight dark matter structures by tracking the phase of extrememassratio inspirals.
['Cosmology', 'Astrophysics', 'and Gravitation', 'General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena']
ExtremeMassRatio Inspirals in Ultralight Dark Matter
{'No Header': '3  \n5', 'Extreme-mass-ratio inspirals in ultra-light dark matter': "Francisco Duque, 1, 2 Caio F. B. Macedo, 3 Rodrigo Vicente, 4 and Vitor Cardoso 2, 5, 6  \n1 Max Planck Institute for Gravitational Physics (Albert Einstein Institute) Am Muhlenberg 1, D-14476 Potsdam, Germany 2 CENTRA, Departamento de F'ısica, Instituto Superior T'ecnico - IST,  \nUniversidade de Lisboa - UL, Avenida Rovisco Pais 1, 1049 Lisboa, Portugal  \nFaculdade de F'ısica, Campus Salin'opolis, Universidade Federal do Par'a, 68721-000, Salin'opolis, Par'a, Brazil  \n4 Institut de Fisica d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology,  \nCampus UAB, 08193 Bellaterra (Barcelona), Spain  \nNiels Bohr International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark 6 Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan  \n(Dated: September 23, 2024)  \nPrevious works have argued that future gravitational-wave detectors will be able to probe the properties of astrophysical environments where binaries coalesce, including accretion disks, but also dark matter structures. Most analyses have resorted to a Newtonian modelling of the environmental effects, which are not suited to study extreme-mass-ratio inspirals immersed in structures of ultralight bosons. In this letter, we use relativistic perturbation theory to consistently study these systems in spherical symmetry. We compute the flux of scalar particles and the rate at which orbital energy is dissipated via gravitational radiation and depletion of scalars, i.e. dynamical friction. Our results confirm that the Laser Interferometer Space Antenna will be able to probe ultra-light dark matter structures by tracking the phase of extreme-mass-ratio inspirals.  \nIntroduction. Progress of gravitational-wave (GW) astronomy requires modeling increasingly complex sources. Extreme-mass-ratio inspirals (EMRIs) - to be observed by the space-borne Laser Interferometer Space Antenna (LISA) mission [1] - correspond to binaries where a stellar-mass compact object orbits around a supermassive black hole (SMBH) with /greaterorsimilar 10 5 M /circledot . EMRIs perform ∼ 10 5 orbital cycles while in the LISA frequency band [24]. They hold a vast scientific potential [5-7], but require the construction of extremely accurate waveforms [8, 9].  \nEMRIs may evolve in non-vacuum environments, such as active galactic nuclei disks [10-19], and should be accreting matter throughout the evolution. The motion of the smaller body generates trailing density wakes in the surrounding medium, acting as dynamical friction that accelerates the inspiral [20, 21]. Environmental effects occur in accretion disks [22, 23] or in galactic dark matter (DM) structures [24-30]. Studies on their relevance for GW astronomy often resort to Newtonian approximations [31-37], but EMRIs inherently probe strong-field gravity [38-44]. Recently, a fully-relativistic framework was developed to study GW emission by EMRIs in spherically symmetric, but otherwise generic, spacetimes [4553]. Relevant differences were found in the GW flux with respect to vacuum for EMRIs immersed in a fluid-like DM halo [45, 46].  \nIn this Letter, we extend the analysis to axionic environments of ultra-light (scalar) fields. New degrees of freedom with mass m Φ ∼ 10 -6 -10 -23 eV appear naturally in extensions of the Standard Model [54-56] and have been suggested as a component of DM [57-59]. They can form self-gravitating structures [60-63] (known as boson stars or DM solitons) which describe well the core of DM halos when the bosons are 'fuzzy' (i.e., their mass is m Φ ∼ 10 -22 -10 -19 eV) [64-69]. Here, their de Broglie wavelength is /lessorsimilar kpc, and they behave as waves  \nat these (and smaller) length scales, alleviating tensions with some galactic observations [70-72]. When their Compton wavelength is comparable to the size of a spinning BH, they can efficiently extract rotational energy from it via superradiance [73-77], and condensate into boson clouds [73, 78-82]. Other mechanisms to form boson condensates include accretion [69, 83, 84] or dynamical capture from a DM halo [85].  \nTheory. Consider a minimally coupled complex scalar Φ described by the action (we use the metric signature ( -+ ++) and geometrized units G = c = 1)  \nS = ∫ d 4 x √ -g ( R 16 π -1 2 ∂ µ Φ ∂ µ Φ ∗ -1 2 V (Φ ∗ Φ)+ L m ) , (1)  \nwith R the Ricci scalar, V a self-interaction potential, and L m describes the matter sector, which interacts with the scalar only via the minimal coupling to the metric g µν . This theory has a global U(1) symmetry with associated Noether current and charge, respectively,  \nJ µ Q = 2 g µν Im[Φ ∗ ∂ ν Φ] , Q = ∫ d 3 x √ -gJ t Q . (2)  \nThe action above results in the Einstein-Klein-Gordon (EKG) system  \nG µν = 8 πT µν , ✷ g Φ = ∂V ∂ Φ ∗ , (3)  \nwhere ✷ g ≡ ∂ µ [ √ -g ∂ µ ( · )] / √ -g , plus a set of equations for the matter sector. Here, G µν is the Einstein tensor and T µν is the total energy-momentum tensor  \nT µν ≡ T Φ µν + T m µν , (4)  \nwith T Φ µν = ∇ ( µ Φ ∗ ∇ ν ) Φ -g µν 2 ( ∇ α Φ ∗ ∇ α Φ+ V ) the scalar field contribution. We focus on sufficiently dilute environments such that V ≈ µ 2 | Φ | 2 , with µ = m Φ / /planckover2pi1 the inverse Compton wavelength of the scalar.  \nEMRI system. The matter sector is the smaller object in the EMRI, modeled by a point particle of mass m p with Lagrangian density  \nL m = m p ∫ d τ 1 √ -g δ (4) [ x -x p ( τ )] , (5)  \nwhere x p ( τ ) is the particle worldline parameterized by its proper time τ . The corresponding energy-momentum tensor is then T µν m = m p ∫ d τ 1 √ -g δ (4) [ x -x p ( τ )] ˙ x µ p ˙ x ν p , where ˙ x µ p ≡ d x µ p / d τ . Treating the secondary as a point particle is a good approximation since the scalar field is light enough to be insensitive to the internal structure of these objects. More precisely, we consider scalars with de Broglie wavelength much larger than the size of the secondary object [38-40].  \nBackground solution. The particle acts as a source of perturbations to a given background solution { ̂ g µν , ̂ Φ } of the EKG system (3), with T m µν = 0. The background spacetime is taken to be static and spherically symmetric  \nd ̂ s 2 ≡ ̂ g µν d x µ d x ν ≈ -A ( r )d t 2 + d r 2 B ( r ) + r 2 dΩ 2 , (6)  \nwith dΩ 2 = d θ 2 + sin 2 θ d ϕ 2 the line element on the 2sphere. The associated background field has the form  \n̂ Φ ≈ φ 0 ( r ) e -iωt , ω ∈ R . (7) This background is used to describe a BH of mass M BH dressed with a scalar configuration with mass M Φ ≡ ∫ d 3 x √ -g T Φ tt /A , which is bound either by the BH gravity ( M BH /greatermuch M Φ ) or by its own self-gravity ( M BH /lessmuch M Φ ) [68]. These two limits correspond, respectively, to boson clouds and DM solitons doped with a parasitic BH. From no-hair results [86, 87], static and spherically symmetric solutions cannot exist, and the scalar field must be swallowed by the BH. However, for ultra-light scalars the accretion timescale can be much larger than the EMRI inspiral (or than an Hubble time [69]), and Eqs. (6) and (7) are a good approximation to the background solution.  \nLinear perturbations. The point particle sources small perturbations to the background solution, resulting in the spacetime metric g µν and scalar field Φ:  \ng µν ≈ ̂ g µν + q δg µν , Φ ≈ ̂ Φ+ q δ Φ , (8) accurate up to O ( q 2 ), with q ≡ m p /M BH . Spherical symmetry allows us to decompose any perturbation into irreducible representations of SO(3), grouped into polar and axial sectors, depending on their behaviour under parity transformations [88-93]. We focus on the polar sector, which dominates GW emission and couples scalar to gravitational perturbations. In the Supplemental Material (SM), we describe how to reduce the system to a set of 'wave-like' partial differential equations, which we solve numerically using techniques previously applied in EMRI modelling [46, 94].  \nAt leading order in q , the particle follows geodesics of the background spacetime ̂ g µν . For simplicity, we focus  \non quasi-circular EMRIs, with particle worldline x µ p ( t ) = ( ε p /A p , r p , π/ 2 , Ω p t ) , with A p ≡ A ( r p ), and where ε p ≡ A p / √ A p -r 2 p Ω 2 p is the particle energy per unit mass and Ω p ≡ √ A ' ( r p ) / 2 r p the orbital angular frequency [95].  \nAdiabatic inspiral and fluxes. As the secondary orbits the central BH, it transfers part of its energy to the scalar environment and to GWs. For small q , the energy dissipated over an orbit is much smaller than the orbital energy, and the secondary evolves adiabatically over a succession of geodesics [8, 96]. To follow its evolution, we only need to compute fluxes in the gravitational and scalar sectors using the linearized EKG system.  \nWe compute energy fluxes at the leading order O ( q 2 ) using T Φ µν and Isaacson's effective energy-momentum tensor T g µν [97],  \n˙ E S L = -σ L lim r → r L r 2 √ AB ∫ dΩ T S tr , (9)  \nwith sector S = { Φ , g } , location L = { H, ∞} (at horizon and spatial infinity, respectively), and σ H, ∞ = -, +. Note that the energy of a non-relativistic scalar field environment is M Φ ≈ Qω , and that the dissipation process results in ˙ M Φ ≈ ˙ Qω . Thus, the fluxes of Noether charge  \n˙ Q L = σ L lim r → r L r 2 √ A B ∫ dΩ J r Q , (10)  \nare necessary to compute the rate of change of the particle energy, ˙ ε p ≡ ∑ S,L,/lscript /lscript ˙ ε S p , L , with /lscript the azimuthal number and [98]  \n/lscript ˙ ε Φ p , L ≈ -/lscript ˙ E Φ L -/lscript ˙ Q L ω m p , /lscript ˙ ε g p , L = -/lscript ˙ E g L m p . (11)  \nThe flux expressions in terms of the perturbation asymptotic amplitudes are given in the SM.  \nFuzzy DM soliton. Let us consider first isolated selfgravitating solutions of the EKG system [60, 63, 86]. We focus on Newtonian boson stars (NBSs), where | ω -µ | /lessmuch µ [98-101], which may describe the inner core of DM halos [57, 58]. They follow the relation [98, 99, 102]  \nM Φ ≈ 9 × 10 8 M /circledot ( 0 . 1 kpc R 98 ) m -2 22 , (12)  \nwith m 22 ≡ m Φ / 10 -22 eV and R 98 the radius enclosing 98% of M Φ (reference values for a Milky-Way-like galaxy [103]). Their oscillating frequency is ω/µ ≈ 1 -0 . 163 ( M Φ µ ) 2 [98, 104], where µ ≈ 4 . 79 yr -1 m 22 , and M Φ µ ∼ 10 -3 m 22 ( M Φ / 10 9 M /circledot ). The metric functions are A ≈ B ≈ 1 + 2 U NBS ( r ), where U NBS is the gravitational potential of the NBS, satisfying ∇ 2 U NBS = µ 2 φ 2 0 , NBS . Accurate fits for U NBS and the radial distribution φ 0 , NBS are known [98] (for details see SM).  \nThe dynamics of NBSs in the non-relativistic regime was studied in Ref. [98], including dissipation via scalar field depletion by quasi-circular binaries with Ω p /lessmuch µ .  \nThere, it was shown that the scalar fluxes are well described by analytic expressions [cf. Eqs. (147)-(149), or our SM]. Here, for configurations with M Φ µ /lessorsimilar 0 . 25, our numerical results for the scalar field fluxes - obtained from solving the linearized fully relativistic EKGsystem agree with these expressions, with relative errors /lessorsimilar 10%, within their range of validity. Even for more compact configurations, they capture the correct dependence of the fluxes on the orbital frequency.  \nFor /lscript -mode scalar field perturbations to be depleted the orbital frequency must satisfy /lscript Ω p > µ -ω , otherwise they are confined. If the orbital dynamics is solely determined by the NBS gravity, it can be shown that Ω p /lessorsimilar µ -ω : the orbital frequency is always too small to generate dipolar ( /lscript = 1) radiation. However, dipolar depletion is possible for more compact relativistic configurations, or in a more realistic description of fuzzy DM cores that includes a central SMBH. Any EMRI observable by LISA (and most near-future GW detectors) must have an orbital frequency Ω p /greatermuch µ ∼ 10 -7 m 22 Hz, thus sourcing relativistic perturbations. Such orbital frequencies cannot be supported by the NBS gravity alone, but they are allowed in compact boson stars [94], or in NBSs hosting a parasitic massive BH, to which we now turn.  \nFuzzy DM soliton hosting a parasitic BH. Large galaxies (like the Milky Way) are expected to host a SMBH [105, 106]. To describe the actual galactic center, we consider a SMBH at the center of a NBS. As mentioned, any spherical scalar configuration around a BH is a transient. However, for an ultra-light boson mass m 22 /lessorsimilar 800, a SMBH with M BH /lessorsimilar 10 6 M /circledot at the center of a soliton core of a DM halo with mass M halo /lessorsimilar 10 15 M /circledot takes more than an Hubble time to accrete 10% of the soliton mass [69]. We can then safely neglect accretion during the observation timescale of an EMRI evolving in these environments.  \nAs shown in SM, if M Φ /greatermuch M BH and M Φ µ /lessmuch 1, the background spacetime can be described by  \nA ≈ ( 1 -2 M BH r ) e 2 U NBS ( r ) , B ≈ 1 -2 m ( r ) r , (13)  \nφ 0 ≈ φ 0 , NBS ( r ) e -2 iµM BH log ( 1 -2 M BH r ) , (14)  \nm ≡ M BH +4 πµ 2 ∫ r 2 M d r ' r ' 2 [ φ 0 , NBS ( r ' )] 2 , (15)  \nwith U NBS and φ 0 , NBS the (isolated) NBS profiles. Note that the central BH dominates the dynamics close to the horizon, where m ( r ) ≈ M BH and U NBS ( r ) ≈ -1 2 M Φ /R 98 /lessmuch 1 [69]. There, the background spacetime reduces to Schwarzschild with a redshifted time coordinate controlled by the NBS compactness, A ≈ (1 -2 M BH /r ) (1 -M Φ /R 98 ). This behaviour was noted for other extended matter distributions around BHs [45, 107, 108] and, as we show in SM, it happens generically when the typical length scale of a Newtonian matter distribution is much larger than the BH size. Far from the  \nwith  \nFIG. 1. Energy loss rate in the gravitational ( | /lscript ˙ ε g p , ∞ | ) and scalar sectors ( | /lscript ˙ /epsilon1 Φ p, ∞ | ) for an EMRI in a NBS with M Φ µ = 0 . 23 hosting a parasitic BH of mass M BH = 0 . 02 M Φ ; the soliton size is R 98 ∼ 8 × 10 3 M BH . Scalar modes are well described by low-frequency analytic expressions from [98] in their regime of validity, and the gravitational contributions follow the vacuum analytic expressions for quadrupole and octupole (all shown as black dots).  \n<!-- image -->  \nhorizon, the NBS gravity dominates the dynamics and, from U NBS /lessmuch 1, we recover the linear superposition of the two Newtonian potentials A ≈ 1+2( U NBS -M BH /r ). The innermost stable circular orbit (ISCO) frequency is Ω ISCO /µ ≈ ( M BH µ ) -1 (1 -M Φ / 2 R 98 ) / 6 √ 6 ( /greatermuch 1 for ultra-light fields), implying that EMRIs in this background can generate relativistic perturbations and dipolar depletion while in the LISA frequency band.  \nThe contributions of scalar field depletion and gravitational radiation to the energy loss rate of an EMRI are shown in Fig. 1, as a function of the distance to the central BH. For illustration purposes, we considered a soliton of mass M Φ µ = 0 . 23 hosting a parasitic SMBH with M BH = 0 . 02 M Φ ; this would accrete the soliton on a timescale ∼ 10 5 yr m -1 22 [69], much larger than the inspiralling timescale or the period of observation for an EMRI in LISA and a boson with m 22 /lessorsimilar 10 2 .  \nFor low-frequency EMRIs, the fluxes are well described by the analytical predictions both in the scalar [98] and gravitational sectors, apart from the resonant behavior in the gravitational sector for Ω p /lessorsimilar ( µ -ω ) //lscript , which is a relativistic feature [94]. The energy dissipated through scalar depletion for orbital frequencies Ω p /greatermuch µ ∼ 10 -7 m 22 Hz can be much larger than the low-frequency predictions, since relativistic corrections become active. The contribution from dipolar scalar field depletion which dominates the energy loss in the scalar channel for all orbital radii where it is active - shows a slightly different behaviour than higher/lscript modes, possibly due to its exceptional character (there is no gravitational radiation for /lscript = 1).  \nFluxes at the BH horizon are several orders of magnitude smaller than at infinity. For low orbital frequencies Ω p /lessorsimilar 16 . 3 //lscript ( M Φ µ ) 2 m 22 nHz the normal modes of  \nthe bosonic configuration are resonantly excited by the EMRI [94, 98, 101]. Even though these scalar perturbations are bound to the soliton, the corresponding GWs carry information about the global structure of the background. The stirring of the bosonic profile by the EMRI at such frequencies also induces resonant accretion, but its contribution to the energy loss rate is orders of magnitude smaller than the flux of energy in GWs at infinity (see Fig. 2 of SM).  \nBoson cloud. Finally, we apply our framework to an EMRI evolving in a spherical (ground-state) boson cloud, a solution of the EKG system where the scalar field structure is sustained by the BH gravity alone. For diluted clouds one may neglect its backreaction on the geometry. The background spacetime is then taken to be (approximately) described by the Schwarzschild metric A ≈ B ≈ 1 -2 M BH /r , whereas the scalar field profile and oscillating frequency are [73, 78-80, 82, 109]  \nφ 0 ≈ Ce -M BH µ 2 r e -2 iµM BH log ( 1 -2 M BH r ) , (16)  \nSuperradiance provides an effective mechanism for formation of these structures [73, 81], which can lead to the growth of clouds with M Φ ≈ 0 . 1 M BH on short timescales [76, 77]. This formation mechanism requires the seed BH to be spinning, and superradiant clouds are born with angular momentum. Thus, they cannot possibly be spherically symmetric. Even though BH spin is key for the growth of these structures, it plays a less crucial role in some of its physics. Spherical boson clouds can be grown via other mechanisms, like accretion [69, 83, 84] or dynamical capture from a DM halo [85].  \nand ω/µ ≈ 1 -1 2 ( M BH µ ) 2 , with C ≡ √ M Φ πM BH ( M BH µ ) 2 .  \nFigure 2 shows the dependence of scalar and GW emission on the orbital radius for an EMRI with M BH µ = 0 . 08 and M Φ = 0 . 1 M BH . Our results are in qualitative agreement with those of Ref. [108]. The small discrepancies are due to differences in the scalar field profile, φ 0 , from the Newtonian approximation employed (c.f., (16)), and the relativistic (numerical) profile used in Ref. [108]. The effect from scalar depletion dominates over the one of gravitational radiation-reaction for sufficiently large orbital radii. Our results suggest that dissipation via dipolar scalar depletion is especially sensitive to the scalar field profile of the background (compare Figs. 1 and 2). Observational prospects. The dephasing of an EMRI gravitational waveform can be roughly estimated by δψ/ (2 π ) ∼ Ω p T 2 obs / ( πT Φ ), where T Φ ∼ /epsilon1 p / ˙ /epsilon1 p [31]. In Fig. 3, choosing T obs = 4yr, we show the region in the parameter space where the dephasing of signal due to the fuzzy DM soliton environment is larger than one cycle. This value is typically used as a rule of thumb for measurability of beyond-vacuum GR effects with LISA [31]. We estimate ˙ ε p by fitting the numerical data in the LISA frequency band, in the range M Φ µ = [0 . 1 , 0 . 25] and M BH /M Φ = [0 . 01 , 0 . 04]. For r p < 15 M BH , the fit underestimates the flux, so our results for the de-  \nFIG. 2. Energy loss rate due to the fluxes at infinity in the gravitational ( | /lscript ˙ ε g p , ∞ | ) and scalar sectors ( | /lscript ˙ /epsilon1 Φ p, ∞ | ) for an EMRI evolving in a boson cloud with M BH µ = 0 . 08 and M Φ = 0 . 1 M BH ; the cloud size is R 98 ∼ 6 × 10 2 M BH . For diluted clouds the test field approximation is valid and the energy dissipated via scalar depletion is linear in M Φ .  \n<!-- image -->  \nFIG. 3. Dephasing of GWs from an EMRI in a DM soliton. Each line stands for a one-cycle dephasing. Above the line the dephasing is larger than one cycle.  \n<!-- image -->  \nive. For a fixed mass-ratio q , the dephasing due to the scalar environment surpasses the criteria threshold provided that M Φ µ and M BH µ are large enough. While some cosmological simulations [6567] indicate that M Φ µ ∼ 10 -3 for low-redshift halos of mass M h ∼ 10 12 M /circledot , there is considerable dispersion in the literature (most probably due to tidal stripping events in the halos history [110]) allowing for M Φ µ /lessorsimilar 0 . 02 √ m 22 . This makes the effect potentially observable for sufficiently compact field configurations of scalars with m Φ /greaterorsimilar 10 -20 eV. For boson clouds, we use directly the numerical results and find that for clouds with mass M Φ /M BH /greaterorsimilar 0 . 01, any scalar field mass M Φ µ /greaterorsimilar 0 . 08 leads to a 1 cycle dephasing for a mass-ratio q /greaterorsimilar 5 × 10 -5 .  \nDiscussion. Our results highlight the importance of relativistic modelling of environmental effects. Although we recover several features previously discussed in a New-  \ntonian context [10, 28, 34, 36, 82, 109], others can only be captured within a fully relativistic framework, such as the behavior of the dipolar mode or the resonant accretion by the central BH.  \nComparison of Figs. 1 and 2 suggests future GW observations can probe and discriminate between different ultra-light dark matter distributions, as the energy flux in the scalar sector is sensitive to the profile of the background scalar distribution. Although we focused the discussion on EMRIs for LISA, our results can be relevant for SMBH binaries in the Pulsar Timing Array (PTA) frequency-band ( ∼ 10 -9 -10 -7 Hz), as the dissipation via scalar depletion dominates over GW radiation reaction for large orbital separations [111-118]. Such prospects motivate further explorations, including a systematic analysis of different ultra-light DM configurations. This could assist in finding semi-analytic descriptions for the energy loss via scalar depletion (i.e. dynamical friction) that would facilitate the implementation of these environmental effects in state-of-the-art EMRI waveform models, like the modular FEW package [8, 119, 120]. Another important extension is the study of more generic orbits, as dynamical friction does not necessarily circularize the orbits. Non-vacuum environments are known to lead to eccentricity excitation [36, 121-123], and DM distributions at galactic centres can induce orbital precession that competes with vacuum GR [124, 125]. Finally, our results can impact binary merger rates, similarly to EMRI formation assisted by gas flows in accretion disks, a mechanism which has been suggested to dominate the fraction of detectable EMRIs by LISA [11].  \nAcknowledgements. We are thankful to Richard Brito for very fruitful discussions. F.D. acknowledges financial support provided by FCT/Portugal through grant No. SFRH/BD/143657/2019. C.F.B.M. would like to thank Funda¸c˜ao Amazˆonia de Amparo a Estudos e Pesquisas (FAPESPA), Conselho Nacional de Desenvolvimento Cient'ıfico e Tecnol'ogico (CNPq) and Coordena¸c˜ao de Aperfei¸coamento de Pessoal de N'ıvel Superior (CAPES) - Finance Code 001, from Brazil, for partial financial support. R.V. is supported by grant no. FJC2021-046551I funded by MCIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR. R.V. also acknowledges the support from the Departament de Recerca i Universitats from Generalitat de Catalunya to the Grup de Recerca 'Grup de F'ısica Te'orica UAB/IFAE' (Codi: 2021 SGR 00649). We acknowledge support by VILLUM Foundation (grant no. VIL37766) and the DNRF Chair program (grant no. DNRF162) by the Danish National Research Foundation. V.C. is a Villum Investigator and a DNRF Chair. V.C. acknowledges financial support provided under the European Union's H2020 ERCAdvanced Grant 'Black holes: gravitational engines of discovery' grant agreement no. Gravitas-101052587. Views and opinions expressed are however those of the author only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 101007855 and No 101131233.  \n- [1] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, et al. , arXiv e-prints (2017), arXiv:1702.00786 [astro-ph.IM].\n- [2] T. Hinderer and E. E. Flanagan, Phys. Rev. D 78 , 064028 (2008), arXiv:0805.3337 [gr-qc].\n- [3] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, Phys. Rev. D95 , 103012 (2017), arXiv:1703.09722 [gr-qc].\n- [4] L. Barack et al. , Class. Quant. Grav. 36 , 143001 (2019), arXiv:1806.05195 [gr-qc].\n- [5] P. Amaro-Seoane et al. (LISA), Living Rev. Rel. 26 , 2 (2023), arXiv:2203.06016 [gr-qc].\n- [6] P. Auclair et al. (LISA Cosmology Working Group), Living Rev. Rel. 26 , 5 (2023), arXiv:2204.05434 [astroph.CO].\n- [7] K. G. Arun et al. (LISA), Living Rev. Rel. 25 , 4 (2022), arXiv:2205.01597 [gr-qc].\n- [8] S. A. Hughes, N. Warburton, G. Khanna, A. J. K. Chua, and M. L. Katz, Phys. Rev. D 103 , 104014 (2021), [Erratum: Phys.Rev.D 107, 089901 (2023)], arXiv:2102.02713 [gr-qc].\n- [9] B. Wardell, A. Pound, N. Warburton, J. Miller,\n- L. Durkan, and A. Le Tiec, Phys. Rev. Lett. 130 , 241402 (2023), arXiv:2112.12265 [gr-qc].\n- [10] E. Barausse, V. Cardoso, and P. Pani, Phys. Rev. D 89 , 104059 (2014), arXiv:1404.7149 [gr-qc].\n- [11] Z. Pan and H. Yang, Phys. Rev. D 103 , 103018 (2021), arXiv:2101.09146 [astro-ph.HE].\n- [12] Z. Pan, Z. Lyu, and H. Yang, Phys. Rev. D 104 , 063007 (2021), arXiv:2104.01208 [astro-ph.HE].\n- [13] A. Derdzinski and L. Mayer, (2022), 10.1093/mnras/stad749, arXiv:2205.10382 [astro-ph.GA].\n- [14] R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. D 102 , 043015 (2020), arXiv:2004.08342 [astro-ph.HE].\n- [15] R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 896 , L44 (2020), arXiv:2006.12611 [astro-ph.HE].\n- [16] R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 125 , 101102 (2020), arXiv:2009.01075 [gr-qc].\n- [17] B. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 892 , L3 (2020), arXiv:2001.01761 [astro-ph.HE].\n- [18] H. Tagawa, Z. Haiman, and B. Kocsis, Astrophys. J. 898 , 25 (2020), arXiv:1912.08218 [astro-ph.GA].\n- [19] H. Tagawa, B. Kocsis, Z. Haiman, I. Bartos, K. Omukai, and J. Samsing, Astrophys. J. 908 , 194 (2021), arXiv:2012.00011 [astro-ph.HE].\n- [20] S. Chandrasekhar, Astrophysics Journal 97 , 255 (1943).  \n- [21] E. C. Ostriker, Astrophys. J. 513 , 252 (1999), arXiv:astro-ph/9810324.\n- [22] N. I. Shakura and R. A. Sunyaev, Astron. Astrophys. 24 , 337 (1973).\n- [23] M. A. Abramowicz and P. C. Fragile, Living Rev. Rel. 16 , 1 (2013), arXiv:1104.5499 [astro-ph.HE].\n- [24] V. C. Rubin, J. Ford, W. K., and N. Thonnard, The Astrophysical Journal 238 , 471 (1980).\n- [25] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462 , 563 (1996), arXiv:astro-ph/9508025.\n- [26] L. Hernquist, The Astrophysical Journal 356 , 359 (1990).\n- [27] E. Barausse, V. Cardoso, and P. Pani, Phys. Rev. D 89 , 104059 (2014).\n- [28] V. Cardoso and A. Maselli, Astron. Astrophys. 644 , A147 (2020), arXiv:1909.05870 [astro-ph.HE].\n- [29] L. Zwick, P. R. Capelo, and L. Mayer, Mon. Not. Roy. Astron. Soc. 521 , 4645 (2023), arXiv:2209.04060 [gr-qc].\n- [30] G. Caneva Santoro, S. Roy, R. Vicente, M. Haney, O. J. Piccinni, W. Del Pozzo, and M. Martinez, Phys. Rev. Lett. 132 , 251401 (2024), arXiv:2309.05061 [gr-qc].\n- [31] B. Kocsis, N. Yunes, and A. Loeb, Phys. Rev. D 84 , 024032 (2011), arXiv:1104.2322 [astro-ph.GA].\n- [32] B. J. Kavanagh, D. A. Nichols, G. Bertone, and D. Gaggero, Phys. Rev. D 102 , 083006 (2020), arXiv:2002.12811 [gr-qc].\n- [33] A. Coogan, G. Bertone, D. Gaggero, B. J. Kavanagh, and D. A. Nichols, Phys. Rev. D 105 , 043009 (2022), arXiv:2108.04154 [gr-qc].\n- [34] P. S. Cole, G. Bertone, A. Coogan, D. Gaggero, T. Karydas, B. J. Kavanagh, T. F. M. Spieksma, and G. M. Tomaselli, Nature Astron. 7 , 943 (2023), arXiv:2211.01362 [gr-qc].\n- [35] M. Garg, A. Derdzinski, L. Zwick, P. R. Capelo, and L. Mayer, Mon. Not. Roy. Astron. Soc. (2022), 10.1093/mnras/stac2711, arXiv:2206.05292 [astroph.GA].\n- [36] G. M. Tomaselli, T. F. M. Spieksma, and G. Bertone, (2023), arXiv:2305.15460 [gr-qc].\n- [37] L. Berezhiani, G. Cintia, V. De Luca, and J. Khoury, (2023), arXiv:2311.07672 [astro-ph.CO].\n- [38] R. Vicente and V. Cardoso, Phys. Rev. D 105 , 083008 (2022), arXiv:2201.08854 [gr-qc].\n- [39] D. Traykova, K. Clough, T. Helfer, E. Berti, P. G. Ferreira, and L. Hui, Phys. Rev. D 104 , 103014 (2021), arXiv:2106.08280 [gr-qc].\n- [40] D. Traykova, R. Vicente, K. Clough, T. Helfer, E. Berti, P. G. Ferreira, and L. Hui, Phys. Rev. D 108 , L121502 (2023), arXiv:2305.10492 [gr-qc].\n- [41] N. Speeney, A. Antonelli, V. Baibhav, and E. Berti, Phys. Rev. D 106 , 044027 (2022), arXiv:2204.12508 [grqc].\n- [42] L. Speri, A. Antonelli, L. Sberna, S. Babak, E. Barausse, J. R. Gair, and M. L. Katz, Phys. Rev. X 13 , 021035 (2023), arXiv:2207.10086 [gr-qc].\n- [43] A. Derdzinski, D. D'Orazio, P. Duffell, Z. Haiman, and A. MacFadyen, Mon. Not. Roy. Astron. Soc. 501 , 3540 (2021), arXiv:2005.11333 [astro-ph.HE].\n- [44] H. Tanaka, T. Takeuchi, and W. R. Ward, Astrophys. J. 565 , 1257 (2002).\n- [45] V. Cardoso, K. Destounis, F. Duque, R. P. Macedo, and A. Maselli, Phys. Rev. D 105 , L061501 (2022), arXiv:2109.00005 [gr-qc].\n- [46] V. Cardoso, K. Destounis, F. Duque,  \n- R. Panosso Macedo, and A. Maselli, Phys. Rev. Lett. 129 , 241103 (2022), arXiv:2210.01133 [gr-qc].\n- [47] K. S. Thorne and A. Campolattaro, 'Non-Radial Pulsation of General-Relativistic Stellar Models. I. Analytic Analysis for L ¿= 2,' Astrophysical Journal, vol. 149, p.591 (1967).\n- [48] S. Detweiler and L. Lindblom, Astrophys. J. 292 , 12 (1985).\n- [49] Y. Kojima, Phys. Rev. D 46 , 4289 (1992).\n- [50] G. Allen, N. Andersson, K. D. Kokkotas, and B. F. Schutz, Phys. Rev. D 58 , 124012 (1998), arXiv:grqc/9704023.\n- [51] O. Sarbach and M. Tiglio, Phys. Rev. D 64 , 084016 (2001), arXiv:gr-qc/0104061.\n- [52] K. Martel and E. Poisson, Phys. Rev. D 71 , 104003 (2005).\n- [53] K. Martel, Phys. Rev. D 69 , 044025 (2004).\n- [54] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38 , 1440 (1977).\n- [55] S. Weinberg, Phys. Rev. Lett. 40 , 223 (1978).\n- [56] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D 81 , 123530 (2010), arXiv:0905.4720 [hep-th].\n- [57] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D 95 , 043541 (2017), arXiv:1610.08297 [astro-ph.CO].\n- [58] L. Hui, Ann. Rev. Astron. Astrophys. 59 , 247 (2021), arXiv:2101.11735 [astro-ph.CO].\n- [59] M. I. Khlopov, B. A. Malomed, and I. B. Zeldovich, Monthly Notices of the Royal Astronomical Society 215 , 575 (1985).\n- [60] D. J. Kaup, Phys. Rev. 172 , 1331 (1968).\n- [61] E. Seidel and W.-M. Suen, Phys. Rev. Lett. 72 , 2516 (1994), arXiv:gr-qc/9309015.\n- [62] J.-w. Lee and I.-g. Koh, Phys. Rev. D 53 , 2236 (1996), arXiv:hep-ph/9507385.\n- [63] S. L. Liebling and C. Palenzuela, Living Rev. Rel. 15 , 6 (2012), arXiv:1202.5809 [gr-qc].\n- [64] D. F. Torres, S. Capozziello, and G. Lambiase, Phys. Rev. D 62 , 104012 (2000), arXiv:astro-ph/0004064.\n- [65] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Phys. 10 , 496 (2014), arXiv:1406.6586 [astro-ph.GA].\n- [66] H.-Y. Schive, M.-H. Liao, T.-P. Woo, S.-K. Wong, T. Chiueh, T. Broadhurst, and W. Y. P. Hwang, Phys. Rev. Lett. 113 , 261302 (2014), arXiv:1407.7762 [astroph.GA].\n- [67] J. Veltmaat, J. C. Niemeyer, and B. Schwabe, Phys. Rev. D 98 , 043509 (2018), arXiv:1804.09647 [astroph.CO].\n- [68] J. Barranco, A. Bernal, J. C. Degollado, A. DiezTejedor, M. Megevand, D. Nunez, and O. Sarbach, Phys. Rev. D 96 , 024049 (2017), arXiv:1704.03450 [grqc].\n- [69] V. Cardoso, T. Ikeda, R. Vicente, and M. Zilh˜ao, Phys. Rev. D 106 , L121302 (2022), arXiv:2207.09469 [gr-qc].\n- [70] B. Moore, Nature (London) 370 , 629 (1994).\n- [71] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490 , 493 (1997), arXiv:astro-ph/9611107.\n- [72] D. H. Weinberg, J. S. Bullock, F. Governato, R. Kuzio de Naray, and A. H. G. Peter, Proc. Nat. Acad. Sci. 112 , 12249 (2015), arXiv:1306.0913 [astroph.CO].\n- [73] R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906 , pp.1 (2015), arXiv:1501.06570 [gr-qc].\n- [74] W. H. Press and S. A. Teukolsky, Nature (London) 238 , 211 (1972).\n- [75] S. R. Dolan, Phys. Rev. D 76 , 084001 (2007), arXiv:0705.2880 [gr-qc].\n- [76] W. E. East and F. Pretorius, Phys. Rev. Lett. 119 , 041101 (2017), arXiv:1704.04791 [gr-qc].\n- [77] W. E. East, Phys. Rev. Lett. 121 , 131104 (2018), arXiv:1807.00043 [gr-qc].\n- [78] S. L. Detweiler, Phys. Rev. D 22 , 2323 (1980).\n- [79] V. Cardoso and S. Yoshida, JHEP 07 , 009 (2005), arXiv:hep-th/0502206.\n- [80] A. Arvanitaki and S. Dubovsky, Phys. Rev. D 83 , 044026 (2011), arXiv:1004.3558 [hep-th].\n- [81] R. Brito, V. Cardoso, and P. Pani, Class. Quant. Grav. 32 , 134001 (2015), arXiv:1411.0686 [gr-qc].\n- [82] D. Baumann, H. S. Chia, J. Stout, and L. ter Haar, JCAP 12 , 006 (2019), arXiv:1908.10370 [gr-qc].\n- [83] L. Hui, D. Kabat, X. Li, L. Santoni, and S. S. C. Wong, JCAP 06 , 038 (2019), arXiv:1904.12803 [gr-qc].\n- [84] K. Clough, P. G. Ferreira, and M. Lagos, Phys. Rev. D100 , 063014 (2019), arXiv:1904.12783 [gr-qc].\n- [85] D. Budker, J. Eby, M. Gorghetto, M. Jiang, and G. Perez, (2023), arXiv:2306.12477 [hep-ph].\n- [86] R. Ruffini and J. A. Wheeler, Phys. Today 24 , 30 (1971).\n- [87] J. D. Bekenstein, Phys. Rev. Lett. 28 , 452 (1972).\n- [88] T. Regge and J. A. Wheeler, Phys. Rev. 108 , 1063 (1957).\n- [89] F. J. Zerilli, Phys. Rev. D 2 , 2141 (1970).\n- [90] N. Sago, H. Nakano, and M. Sasaki, Phys. Rev. D 67 , 104017 (2003), arXiv:gr-qc/0208060.\n- [91] S. L. Detweiler and E. Poisson, Phys. Rev. D 69 , 084019 (2004), arXiv:gr-qc/0312010.\n- [92] K. Martel and E. Poisson, Phys. Rev. D 71 , 104003 (2005), arXiv:gr-qc/0502028.\n- [93] DLMF, 'NIST Digital Library of Mathematical Functions,' http://dlmf.nist.gov/, Release 1.0.26 of 2020-0315, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.\n- [94] C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Phys. Rev. D 88 , 064046 (2013), arXiv:1307.4812 [gr-qc].\n- [95] V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Phys. Rev. D 79 , 064016 (2009), arXiv:0812.1806 [hep-th].\n- [96] L. Barack and A. Pound, Rept. Prog. Phys. 82 , 016904 (2019), arXiv:1805.10385 [gr-qc].\n- [97] R. A. Isaacson, Phys. Rev. 166 , 1272 (1968).\n- [98] L. Annulli, V. Cardoso, and R. Vicente, Phys. Rev. D 102 , 063022 (2020), arXiv:2009.00012 [gr-qc].\n- [99] P.-H. Chavanis, Phys. Rev. D 84 , 043531 (2011), arXiv:1103.2050 [astro-ph.CO].\n- [100] F. Kling and A. Rajaraman, Phys. Rev. D 96 , 044039 (2017), arXiv:1706.04272 [hep-th].\n- [101] J. H.-H. Chan, S. Sibiryakov, and W. Xue, JHEP 08 ,  \n- 045 (2023), arXiv:2304.13054 [astro-ph.CO].\n- [102] M. Boˇskovi'c, F. Duque, M. C. Ferreira, F. S. Miguel, and V. Cardoso, Phys. Rev. D 98 , 024037 (2018), arXiv:1806.07331 [gr-qc].\n- [103] I. De Martino, T. Broadhurst, S. H. H. Tye, T. Chiueh, and H.-Y. Schive, Phys. Dark Univ. 28 , 100503 (2020), arXiv:1807.08153 [astro-ph.GA].\n- [104] F. S. Guzman and L. A. Urena-Lopez, Phys. Rev. D 69 , 124033 (2004), arXiv:gr-qc/0404014.\n- [105] J. Kormendy and D. Richstone, Ann. Rev. Astron. Astrophys. 33 , 581 (1995).\n- [106] J. Kormendy and L. C. Ho, Ann. Rev. Astron. Astrophys. 51 , 511 (2013), arXiv:1304.7762 [astro-ph.CO].\n- [107] E. Figueiredo, A. Maselli, and V. Cardoso, Phys. Rev. D 107 , 104033 (2023), arXiv:2303.08183 [gr-qc].\n- [108] R. Brito and S. Shah, Phys. Rev. D 108 , 084019 (2023), arXiv:2307.16093 [gr-qc].\n- [109] D. Baumann, H. S. Chia, R. A. Porto, and J. Stout, Phys. Rev. D 101 , 083019 (2020), arXiv:1912.04932 [grqc].\n- [110] H. Y. J. Chan, E. G. M. Ferreira, S. May, K. Hayashi, and M. Chiba, Mon. Not. Roy. Astron. Soc. 511 , 943 (2022), arXiv:2110.11882 [astro-ph.CO].\n- [111] G. Agazie et al. (NANOGrav), Astrophys. J. Lett. 952 , L37 (2023), arXiv:2306.16220 [astro-ph.HE].\n- [112] A. Afzal et al. (NANOGrav), Astrophys. J. Lett. 951 , L11 (2023), arXiv:2306.16219 [astro-ph.HE].\n- [113] G. Agazie et al. (NANOGrav), Astrophys. J. Lett. 951 , L8 (2023), arXiv:2306.16213 [astro-ph.HE].\n- [114] J. Antoniadis et al. (EPTA, InPTA:), Astron. Astrophys. 678 , A50 (2023), arXiv:2306.16214 [astro-ph.HE].\n- [115] C. Smarra et al. (European Pulsar Timing Array), (2023), arXiv:2306.16228 [astro-ph.HE].\n- [116] J. Antoniadis et al. , (2023), arXiv:2306.16227 [astroph.CO].\n- [117] A. Ghoshal and A. Strumia, (2023), arXiv:2306.17158 [astro-ph.CO].\n- [118] J. Ellis, M. Fairbairn, G. Hutsi, J. Raidal, J. Urrutia, V. Vaskonen, and H. Veermae, (2023), arXiv:2306.17021 [astro-ph.CO].\n- [119] M. L. Katz, A. J. K. Chua, L. Speri, N. Warburton, and S. A. Hughes, Phys. Rev. D 104 , 064047 (2021), arXiv:2104.04582 [gr-qc].\n- [120] L. Speri, M. L. Katz, A. J. K. Chua, S. A. Hughes, N. Warburton, J. E. Thompson, C. E. A. ChapmanBird, and J. R. Gair, (2023), arXiv:2307.12585 [gr-qc].\n- [121] R. P. Nelson, Astron. Astrophys. 443 , 1067 (2005), arXiv:astro-ph/0508486.\n- [122] V. Cardoso, C. F. B. Macedo, and R. Vicente, Phys. Rev. D 103 , 023015 (2021), arXiv:2010.15151 [gr-qc].\n- [123] M. Rahman, S. Kumar, and A. Bhattacharyya, (2023), arXiv:2306.14971 [gr-qc].\n- [124] K. Destounis, A. Kulathingal, K. D. Kokkotas, and G. O. Papadopoulos, Phys. Rev. D 107 , 084027 (2023), arXiv:2210.09357 [gr-qc].\n- [125] K. Destounis, F. Angeloni, M. Vaglio, and P. Pani, (2023), arXiv:2305.05691 [gr-qc].", "Supplemental Material to 'Axion weak leaks: extreme mass-ratio inspirals in ultralight dark matter'": "Francisco Duque, 1, 2 Caio F. B. Macedo, 3 Rodrigo Vicente, 4 and Vitor Cardoso 2, 5  \n1 Max Planck Institute for Gravitational Physics (Albert Einstein Institute) Am Muhlenberg 1, D-14476 Potsdam, Germany 2 CENTRA, Departamento de F'ısica, Instituto Superior T'ecnico - IST,  \nUniversidade de Lisboa - UL, Avenida Rovisco Pais 1, 1049 Lisboa, Portugal  \n3 Faculdade de F'ısica, Campus Salin'opolis, Universidade Federal do Par'a, 68721-000, Salin'opolis, Par'a, Brazil  \n4 Institut de Fisica d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology,  \nCampus UAB, 08193 Bellaterra (Barcelona), Spain  \n5 Niels Bohr International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark  \nIn this Supplemental Material, we give a detailed derivation of the evolution equations for the gravitational and scalar field linear perturbations in sphericallysymmetric, but otherwise generic, backgrounds. Our procedure is based on previous works on relativistic stars [1, 2], and our own studies of extreme mass-ratio inspirals (EMRIs) immersed in dark-matter (DM) halos constituted by an anisotropic fluid [3, 4]. We also explicitly express the different fluxes in terms of asymptotic amplitudes of the field perturbations. Finally, we provide more details on the background (ultralight DM) solutions considered in our work: self-gravitating DM solitons (either in isolation, or hosting a parasitic black hole), and boson clouds sustained by black holes (BHs).", 'I. A CRASH COURSE ON BH PERTURBATION THEORY': 'The natural framework to study EMRIs is BH perturbation theory [5-7]. Consider a secondary compact object of mass m p sourcing perturbations to a particular (non-vacuum) background solution { ( ̂ M , ̂ g µν ) , ̂ Φ } of the Einstein-Klein-Gordon (EKG) system, representing a massive BH of mass M BH with a scalar field environment of mass M Φ ,  \ng µν = ̂ g µν + q δg µν + O ( q 2 ) , (1) Φ = Φ+ q δ Φ+ O ( q 2 ) . (2)  \n̂ Our perturbation scheme is built over the mass-ratio q ≡ m p /M BH .', 'I.1. Background solution': 'Wefocus on spherically-symmetric backgrounds, which in spherical coordinates have line element  \nd ̂ s 2 ≡ ̂ g µν d x µ d x ν ≈ -A ( r )d t 2 + d r 2 B ( r ) + r 2 dΩ 2 , (3)  \nwhere dΩ ≡ d θ 2 + sin 2 θ d ϕ 2 is the line element on the 2-sphere, with the background scalar field  \n̂ Φ ≈ φ 0 ( r ) e -iωt , (4)  \nwhere ω is (in general) a complex-valued parameter found by solving a boundary value problem, imposing regularity at the BH horizon (or at the origin) and at spatial infinity. When taking the backreaction into account, an exact solution of this form with Im ω = 0 cannot possibly exist, since no-hair theorems prevent these spacetimes to be exactly static. Nevertheless, the previous expressions can still describe dynamical solutions for an interval of time ∆ t /lessorsimilar τ backg ≡ ( M BH /M Φ ) | Im ω | -1 (see, e.g., Refs. [8, 9]). Although not needed, for simplicity, in this work we focus on scalar fields that are light enough ( µM BH /lessmuch 1), so that accretion is suppressed [9] and the inspiral happens on a timescale τ insp /lessmuch min ( τ backg , | Im ω | -1 ) . For such light fields, we can effectively use Im ω ≈ 0 to compute the energy and angular momentum loss rate in EMRIs.  \n/negationslash', 'I.2. Multipolar expansion': "The spherical symmetry of the background solution guarantees that any linear perturbation can be decomposed into irreducible representations of SO(3). For the scalar field, this corresponds to the standard expansion in (scalar) spherical harmonics [10]  \nδ Φ = e -iωt r ∑ /lscript,m φ /lscriptm + ( t, r ) Y /lscriptm ( θ, ϕ ) , (5)  \nδ Φ ∗ = e iωt r ∑ /lscript,m φ /lscriptm -( t, r ) Y /lscriptm ( θ, ϕ ) . (6)  \nwhere ∑ /lscript,m ≡ ∑ ∞ /lscript =0 ∑ /lscript m = -/lscript . We treat the δ Φ and δ Φ ∗ perturbations as independent and factor out the dependence on the background frequency ω , which as we will see below will allow us to obtain linearized equations of motion whose coefficients do not depend on time. Although in the main text we consider only scalar field structures diluted enough, so that self-interactions are negligible, here for generality we include a U(1)-invariant self-interaction potential V . So, it is also useful to expand  \nit, and its first derivative, as  \nV ≈ ̂ V + ∑ /lscript,m δV /lscriptm ( t, r ) Y /lscriptm ( θ, ϕ ) , (7)  \nU ≡ ∂V ∂ | ̂ Φ | 2 ≈ ̂ U + ∑ /lscript,m δU /lscriptm ( t, r ) Y /lscriptm ( θ, ϕ ) . (8) For diluted scalar field structures, we have simply  \n̂ V ≈ µ 2 φ 2 0 , ̂ U ≈ µ 2 , (9) δV /lscriptm = φ 0 r ( φ /lscriptm + + φ /lscriptm -) , δU /lscriptm ≈ 0 , (10)  \nδg axial = ∑ /lscript,m √ 2 /lscript ( /lscript +1) r [ i h /lscriptm 1 c /lscriptm -h /lscriptm 0 c 0 /lscriptm + √ ( /lscript +2)( /lscript -1) 2 h /lscriptm 2 d /lscriptm ] , (11)  \nδg polar = ∑ /lscript,m AH /lscriptm 0 a 0 /lscriptm -i √ 2 H /lscriptm 1 a 1 /lscriptm + 1 B H /lscriptm 2 a /lscriptm + √ 2 K /lscriptm g /lscriptm + √ 2 /lscript ( /lscript +1) r ( h ( e ) /lscriptm 1 b 1 /lscriptm -ih ( e ) /lscriptm 0 b 0 /lscriptm ) + (√ ( /lscript +2)( /lscript +1) /lscript ( /lscript -1) 2 f /lscriptm -/lscript ( /lscript +1) √ 2 g /lscriptm ) G /lscriptm . (12)  \nWe omit the dependence on ( t, r, θ, ϕ ) to avoid cluttering, but the mode perturbations h /lscriptm 1 , h /lscriptm 0 , . . . are functions of only ( t, r ), while a 0 /lscriptm , a /lscriptm , . . . are the ten tensor spherical harmonics depending only on ( θ, ϕ ), which are orthonormal on the 2-sphere, as defined in Ref. [6].  \nIn a similar fashion, the energy-momentum tensor can be expanded as  \nT = ∑ /lscript,m A 0 /lscriptm a 0 /lscriptm + A 1 /lscriptm a 1 /lscriptm + A /lscriptm a /lscriptm + B 0 /lscriptm b 0 /lscriptm + B /lscriptm b /lscriptm + Q 0 /lscriptm c 0 /lscriptm + Q /lscriptm c /lscriptm + D /lscriptm d /lscriptm + F /lscriptm f /lscriptm + G /lscriptm g /lscriptm . (13)  \nFor a given source the expansion coefficients can be obtained by projecting the energy-momentum tensor on the respective spherical harmonic, e.g., A 0 /lscriptm = ( a 0 /lscriptm , T ), with ( · , · ) the inner product on the 2-sphere  \n( R /lscript ' m ' , S /lscriptm ) = ∫ S 2 dΩ ( R /lscript ' m ' µν ) ∗ S /lscriptm λρ η µλ η νρ , (14)  \nwhere  \nη µν = diag ( -1 , 1 , r 2 , r 2 sin 2 θ ) . (15)", 'I.3. Gauge invariance': "The theory of general relativity (GR) is invariant under diffeomorphisms, which translates into the gauge invariance of linear perturbations:  \nx µ → x ' µ = x µ + q ξ µ , δg µν → δg ' µν = δg µν -2 ∇ ( µ ξ ν ) , δ Φ → δ Φ ' = δ Φ -ξ µ ∂ µ Φ ,  \n̂ where ξ µ is the vector field generating the (infinitesimal) diffeomorphism; we can use its four components to impose some coefficients in the expansion of δg µν to vanish. Note that ξ µ can also be expanded in the polar and axial components  \nξ µ polar = ∑ /lscript,m ( -1 A ξ /lscriptm t , Bξ /lscriptm r , 0 , 0 ) Y /lscriptm + ξ /lscriptm Ω r 2 sin θ (0 , 0 , sin θ ∂ θ Y /lscriptm , ∂ ϕ Y /lscriptm ) , ξ µ axial = ∑ /lscript,m ξ /lscriptm ax r sin θ (0 , 0 , ∂ ϕ Y /lscriptm , -sin θ ∂ θ Y /lscriptm ) ,  \nwhere the ξ /lscriptm t , ξ /lscriptm r , ξ /lscriptm Ω , ξ /lscriptm ax are only functions of ( t, r ). Then, in terms of tensor spherical harmonics, we find  \nThe gravitational perturbations can be expanded in the basis of the (ten) tensor spherical harmonics, which can be grouped into polar (electric/even) or axial (magnetic/odd) types, depending on their behavior under parity transformations ( θ, ϕ ) → ( π -θ, π + ϕ ) [6],  \n2 ∇ ξ = (2 ∂ t ξ t -A ' Bξ r ) a 0 -i √ 2 ( ∂ r ξ t + ∂ t ξ r -A ' A ξ t ) a 1 + ( 2 ∂ r ξ r + B ' B ξ r ) a -i √ 2 /lscript ( /lscript +1) r ( ξ t + ∂ t ξ Ω ) b 0 + √ 2 /lscript ( /lscript +1) r ( ∂ r ξ Ω + ξ r -2 r ξ Ω ) b + √ 2( /lscript +1) /lscript ∂ t ξ ax c 0 -i √ 2( /lscript +1) /lscript ( ∂ r ξ ax -ξ ax r ) c + i √ 2 ( /lscript +2)( /lscript +1)( /lscript -1) r ξ ax d + √ 2 ( /lscript +2)( /lscript +1) /lscript ( /lscript -1) r 2 ξ Ω f + √ 2 r 2 (2 rBξ r -/lscript ( /lscript +1) ξ Ω ) g , (16)  \nwhere the prime denotes a derivative with respect to r and, henceforth, we omit the ( /lscript, m ) indices to avoid cluttering. If we pick ξ µ judiciously, we can eliminate one component of the metric perturbations in the axial sector and three in the polar sector. In this work, we adopt the Regge-Wheeler gauge, where the gauge is fixed by setting to zero all terms involving angular derivatives of the highest order [5, 11]  \nh 2 = h ( e ) 0 = h ( e ) 1 = G = 0 . (17)  \nHowever, for /lscript ≤ 1 this choice does not completely fix the gauge, because some of the tensor spherical harmonics are identically zero: more precisely, b 0 = b = c 0 = c = 0 for /lscript = 0, and d = f = 0 for /lscript ≤ 1. The /lscript ≤ 1 modes do not contribute to the radiative degrees of freedom of the gravitational field, and in vacuum can be removed by a gauge transformation [5, 12]. However, in the presence of an environment, they carry physical significance through the matter sector and generate dipolar scalar depletion.", 'I.4. Point-particle source': "The secondary body is modelled as a point-particle with energy-momentum tensor  \nT µν p = m p r 2 sin θ √ B A dt p dτ dx µ dt dx ν dt × δ [ r -r p ( t )] δ [ θ -θ p ( t )] δ [ ϕ -ϕ p ( t )] , (18)  \nWe take the body to be on circular geodesics of the background at a constant radius r p ,  \nr p ( t ) = r p , θ p ( t ) = π 2 , ϕ p ( t ) = Ω p t . (19)  \nThe orbital angular frequency is Ω p = √ A ' p / 2 r p , where A p ≡ A ( r p ) [13]. These orbits admit two conserved quantities: the energy and angular momentum per unit rest mass, respectively, given by [13]  \nε p = A p dt p dτ = A p √ A p -r 2 p Ω 2 p , (20)  \nl p = r 2 p dϕ p dτ = r 2 p Ω p √ A p -r 2 p Ω 2 p . (21)  \nFor circular orbits, the tensor harmonics expansion (13) for the energy-momentum tensor of the secondary greatly simplifies. We have A /lscriptm = A (1) /lscriptm = B /lscriptm = Q /lscriptm = 0, while the non-vanishing coefficients can be expressed in terms of the orbital parameters as follows  \nA 0 /lscriptm = m p √ ABε p r 2 Y /star /lscriptm δ r p , B 0 /lscriptm = m p i √ ABl p r 3 √ ( n +1) δ r p ∂ φ Y /star /lscriptm , Q 0 /lscriptm = -m p √ ABl p r 3 √ ( n +1) δ r p ∂ θ Y /star /lscriptm , G /lscriptm = m p l 2 p √ AB r 4 √ 2 ε p δ r p Y /star /lscriptm , D /lscriptm = m p il 2 p √ AB ε p r 4 √ 2 n ( n +1) δ r p ∂ θφ Y /star /lscriptm , F /lscriptm = m p l 2 p √ AB r 4 2 ε p √ 2 n ( n +1) δ r p ( ∂ φφ -∂ θθ ) Y /star /lscriptm , (22)  \nwhere n = /lscript ( /lscript +1) / 2 -1, Y /star /lscriptm = Y /star /lscriptm ( π/ 2 , 0) and δ r p = δ ( r -r p ).", 'II. EVOLUTION EQUATIONS': 'We now have all the necessary ingredients to obtain the equations of motion to evolve the perturbations.', 'II.1. /lscript ≥ 2': "We start with the /lscript ≥ 2 modes. We will only focus on the polar sector, where gravitational perturbations couple to matter ones. First, redefine  \nH 0 = K + r A S , (23)  \nH 1 = r A ˜ H 1 . (24)  \nThis choice is inspired in studies of perturbations around relativistic stars [2] and the factoring of r and A  \ncaptures the asymptotic behavior at, respectively, large distances and near the BH horizon. Then, we define  \nE µν = G (1) µν -8 π ( T Φ(1) µν + T p µν ) , (25)  \nand follow the same steps as in Ref. [4]  \n- · E θθ -E ϕϕ / sin 2 θ gives an algebraic relation for H 2  \nH 2 = K + r A S -8 πr 2 √ n ( n +1) / 2 F /lscriptm . (26)  \n∂ ˜ H 1 ∂t = √ A B ∂S ∂r ∗ + A r S + AA ' r K -8 π A 2 r 2 φ ' 0 ( φ + + φ -) -4 π A √ n ( n +1) / 2 (2 A + rA ' ) F /lscriptm . (27)  \n- · E tt -AB E rr gives the first second-order 'wavelike' equation for K  \n-∂ 2 K ∂ 2 t + ∂ 2 K ∂ 2 r ∗ +2 √ AB r ∂K ∂r ∗ -[ A r 2 ( /lscript ( /lscript +1) + 2 B -2 + 2 r B ' ) -2 A ' B r +8 πω 2 φ 2 0 +8 πABφ ' 2 0 ] K = [ 2 B r +2 B ' -2 A ' B A +8 π r A ω 2 φ 2 0 +8 π B r φ ' 2 0 ] S +16 π AB r 2 φ ' 0 ( φ + + φ -) -8 πAδV -8 π A 0 /lscriptm -8 π AB r √ n ( n +1) / 2 ∂ F /lscriptm ∂r + 4 π √ n ( n +1) / 2 [ 2 rA ' B -A ( /lscript ( /lscript +1) + 2 + 2 B +2 rB ' ) -8 πr 2 ω 2 φ 2 0 -8 πr 2 ABφ ' 2 0 +8 πr 2 A ̂ V ] F /lscriptm (28)  \n- · E θθ + E ϕϕ / sin 2 θ gives another second-order 'wavelike' equation for S  \n-∂ 2 S ∂ 2 t + ∂ 2 S ∂ 2 r ∗ -[ A r 2 ( /lscript ( /lscript +1) + 2 B -2 + r 2 B ' ) + A ' B 2 r -8 πω 2 φ 2 0 +8 πABφ ' 2 0 +8 πA ̂ V ] S = A r [ 2 A ' B r -A ' B ' -2 A '' B ] K + 8 π r 2 [( A 2 ( B ' -2 B r ) +2 AA ' B ) φ ' 0 +2 A 2 Bφ '' 0 ] ( φ + + φ -) +8 π √ 2 A 2 r G /lscriptm -8 πAr √ n ( n +1) / 2 ∂ 2 F /lscriptm ∂t 2 + 4 π AB √ n ( n +1) / 2 (2 A + rA ' ) ∂ F /lscriptm ∂r + 4 π √ n ( n +1) / 2 A r [ 2 rA ' B + r 2 B A ' 2 A -A ( /lscript ( /lscript +1) -4 + 4 B ) + 16 πr 2 ω 2 φ 2 0 -16 πr 2 ABφ ' 2 0 -16 πr 2 A ̂ V ] F /lscriptm (29)  \n- · Finally, the Klein-Gordon perturbed up to first-order gives  \n-∂ 2 φ + ∂ 2 t + ∂ 2 φ + ∂ 2 r ∗ +2 i ω ∂φ + ∂t + [ ω 2 -A r 2 /lscript ( /lscript +1) -AB ' + A ' B 2 r -A ̂ U ] φ + -8 πABφ ' 2 0 ( φ + + φ -) -r A φ 0 δU = -2 irωφ 0 ∂K ∂t -i r 2 A ωφ 0 ∂S ∂t + i r 2 A √ B A ωφ 0 ∂ ˜ H 1 ∂r ∗ -[ r 2 A ω 2 φ 0 + rBφ ' 0 ( r A ' 2 A -r B ' 2 B -2 ) -r 2 Bφ '' 0 ]( S + A r K ) + i r 2 ω 2 A [( 6 B r + B ' -3 A ' B A ) φ 0 +4 Bφ ' 0 ] H 1 + i 4 πr 3 ωφ 0 √ n ( n +1) / 2 ∂ F /lscriptm ∂t -4 πr 3 ABφ ' 0 √ n ( n +1) / 2 ∂ F /lscriptm ∂r -4 π r 2 AB √ n ( n +1) / 2 [( 4 + r B ' B ) φ ' 0 +2 rφ '' 0 ] F /lscriptm , (30)  \n- · E rθ gives a relation for ∂ ˜ H 1 /∂t  \nThe equation for φ -is similar but with the complex conjugate coefficients multiplying the perturbations.", 'II.2. /lscript = 1': "As mentioned above, for /lscript = 1, our gauge choice in Eq. (17) is not completely fixed because some of the spherical harmonics are identically 0. In particular, g = 0 so we can fix one more gravitational perturbation function to 0. We follow the original choice of Zerilli [5] and pick K = 0. Then  \n- · E tr yields an algebraic relation for H 1  \nH 1 = -1 AB + rAB ' +4 πr 2 ( A ̂ V -ω 2 φ 2 0 + ABφ ' 2 0 ) [ rA ∂H 2 ∂t + i 4 πωA ( φ 0 + rφ ' 0 ) ( φ + -φ -) -i 4 πω √ A B rφ 0 ( ∂φ + ∂r ∗ -∂φ -∂r ∗ ) -4 πArφ ' 0 ( ∂φ + ∂t + ∂φ -∂t ) . ] (31)  \n- · E rθ gives a relation between ∂H 0 /∂r ∗ and the other variables which can then be substituted in E rr to obtain an algebraic relation for H 0  \nH 0 = 1 2 A ( B -1) -rA ' B +8 πr 2 ω 2 φ 2 0 [ 2 rB ∂H 1 ∂t + ( 2 A ( B -1) + rA ' B +8 πr 2 ̂ V -8 πr 2 ω 2 φ 0 ) H 0 -8 π [ ( rω 2 φ 0 + ABφ ' 0 ) ( φ + + φ -) + r √ ABφ ' 0 ( ∂φ + ∂r ∗ + ∂φ -∂r ∗ ) + irωφ 0 ( ∂φ + ∂t + ∂φ -∂t )] (32)  \nThe two equations above can then be substituted in E tt to obtain an equation for ∂H 2 /∂r ∗  \n∂H 2 ∂r ∗ = -1 r √ A B [ ( 1 + B + rB ' +4 πr 2 Bφ ' 2 0 ) H 2 + ( B -1 + rB ' +4 πr 2 ( ̂ V + Bφ ' 0 )) H 0 ] + 4 π r √ AB [ ( rω 2 φ 0 -ABφ ' 0 ) ( φ + + φ -) + irωφ 0 ( ∂φ + ∂t + ∂φ -∂t ) + √ ABrφ ' 0 ( ∂φ + ∂r ∗ + ∂φ -∂r ∗ )] + 8 πr √ AB A 0 , (33)  \n- · Finally, the perturbed Klein-Gordon equation in this gauge choice is written as  \n-∂ 2 φ + ∂t 2 + ∂ 2 φ + ∂ 2 r ∗ +2 i ω ∂φ + ∂t + [ ω 2 -A r 2 /lscript ( /lscript +1) -AB ' + A ' B 2 r -A ̂ U ] φ + -rAδUφ + = r 2 √ ABφ ' 0 ∂H 0 ∂r ∗ + irω √ B A φ 0 ∂H 1 ∂r ∗ + r 2 √ ABφ ' 0 ∂H 2 ∂r ∗ -i r 2 ωφ 0 ∂H 0 ∂t -rBφ ' 0 ∂H 0 ∂t -i r 2 ωφ 0 ∂H 2 ∂t -rω 2 φ 0 H 0 + iω [ φ 0 ( 2 B + r A ( AB ' -A ' B ) ) +2 Brφ ' 0 ] H 1 + [( 2 AB + r 2 ( AB ' + A ' + B ) ) φ ' 0 + rABφ '' 0 ] H 2 , (34)  \nand the same equation for φ -but with complex conjugated coefficients.", 'III. BACKGROUND SOLUTIONS': "In this section, we provide more details on the different backgrounds used for the metric and the scalar field  \ndistribution. First, let us rewrite the g rr metric function  \nin terms of the Misner-Sharp mass function  \nB ( r ) = 1 -2 m ( r ) r . (35)  \nm ' = 2 πr 2 ( 1 -2 m r ) ( | φ ' 0 | 2 + µ 2 + ω 2 A 1 -2 m r | φ 0 | 2 ) , (36)  \n( 1 -2 m r ) φ '' 0 + √ 1 -2 m r A 1 r 2 ( r 2 √ A ( 1 -2 m r ) ) ' φ ' 0 = ( µ 2 -ω 2 A ) φ 0 . (38)  \nlog ' ( A ) = 2 m r 2 ( 1 -2 m r ) +2 πr ( | φ ' 0 | 2 -µ 2 -ω 2 A 1 -2 m r | φ 0 | 2 ) , (37)", 'III.1. Newtonian Boson Stars': "Newtonian Boson Stars (NBS) are localized, lowcompactness, self-gravitating solutions of the system above, where the scalar field is necessarily weak, i.e. | φ 0 , NBS ( r ) | /lessmuch 1. When this is verified, the background metric function is approximately [14]  \nA ( r ) ≈ 1 + 2 U NBS ( r ) , (39)  \nwhere | U NBS ( r ) | /lessmuch 1 is the Newtonian gravitational potential of the NBS. The frequency of the fundamental mode is given by  \nω = µ -γ , (40)  \nwhere 0 < γ /lessmuch µ . As we will see below, Eqs. (36)-(38) can be reduced to the Schrodinger-Poisson system for φ 0 , NBS and U NBS , respectively. This forms an eigenvalue problem for γ where one needs to impose regularity of φ 0 , NBS and U NBS at the origin and demand asymptotic flatness at large distances. It is possible to show that for large r , the scalar field falls off according to a Yukawa potential φ 0 , NBS ∼ e -√ 2 µγr /r , while the gravitational potential behaves as | U NBS | ∼ M Φ /r , where M Φ is the total mass of the scalar configuration  \nM Φ = 4 πµ 2 ∫ ∞ 0 r ' 2 | φ 0 ( r ' ) | 2 . (41)  \nRef. [14] obtained semi-analytic descriptions of the field and gravitational potential accurate up to 1 % with respect to numerical solutions. We reproduce it here for completeness  \nU NBS = µ 2 M 2 Φ f ( x ) , (42) f = a 0 +11 a 0 r 1 x + ∑ 9 i =2 a i x i -x 10 ( x + r 1 ) 11 , x = µ 2 M Φ r , r 1 = 1 . 288 , a 0 = -5 . 132 , a 2 = -143 . 279 , a 3 = -645 . 326 , a 4 = 277 . 921 , a 5 = -2024 . 838 , a 6 = 476 . 702 , a 7 = -549 . 051 , a 8 = -90 . 244 , a 9 = -13 . 734 .  \nand  \nφ 0 = µ 2 M 2 Φ g ( x ) , (43) g = e -0 . 570459 x ∑ 8 i =0 b i x i + b f x 9 . 6 ( x + r 2 ) 9 , x = µ 2 M Φ r , r 2 = 1 . 182 , b 0 = 0 . 298 , b 1 = 2 . 368 , b 2 = 10 . 095 , b 3 = 12 . 552 , b 4 = 51 . 469 , b 5 = -8 . 416 , b 6 = 54 . 141 , b 7 = -6 . 167 , b 8 = 8 . 089 , b f = 0 . 310 . (44)  \nand the frequency of the fundamental mode is  \nγ ≈ 0 . 163 M 2 Φ µ 3 . (45)  \nAs mentioned in the main text, NBSs exhibit a scale invariant mass-radius relation  \nM Φ ≈ 9 × 10 8 M /circledot ( 0 . 1 kpc R 98 ) m -2 22 , (46)  \nwhere we recall m 22 ≡ m Φ / 10 -22 eV, R 98 is the radius enclosing 98% of M Φ .  \nAn important property of NBSs is that near their center, in particular for radius r /lessorsimilar 0 . 1 R 98 , both the scalar field profile and the gravitational potential are approximately constant  \nφ 0 , NBS (0) ≈ 0 . 0662 M 2 Φ µ 3 , (47)  \nU NBS (0) ≈ -1 2 M Φ R 98 . (48)  \nAll the features discussed are illustrated in Fig. 1 of Ref. [14].", 'III.2. Parasitic Black Holes': "Now we want to place a BH in the center of a NBS. As discussed in the main text, no-hair theorems [15, 16] prevent the existence of a static, spherically-symmetric solution in this setup. Consequently, the scalar field distribution will be accreted by the BH and Im( ω ) = 0. However,  \n/negationslash  \nNumerical Relativity results supported by analytic estimates have shown that when the mass of the NBS is much larger than that of the BH, i.e., M BH /lessmuch M Φ , and the star is Newtonian, i.e., µM Φ /lessmuch 1, the accretion timescale is larger than the Hubble time [9]. In particular [9]  \nIm( ω ) ∼ 10 -1 µ 3 M 2 Φ , (49)  \nand the time τ 10% necessary for the BH to accrete 10% of the initial mass of the NBS, M Φ , 0 , is [9]  \nτ 10% 10 Gyr /lessorsimilar 6 . 6 ( 10 10 M Φ , 0 ) 5 ( 10 -22 eV µ ) 6 . (50)  \nTherefore, for the inspiral timescale of an EMRI we can take Im( ω ) ≈ 0.  \nThe hierarchy of length scales determined by the event horizon and the NBS's radius allows us to separate the problem into two distinct regions. Close to the event horizon, the gravitational field is dominated by BH, while at far-away distances, gravity is dictated by the self-gravity of the NBS. In between, there is a transition radius where the two gravitational fields are comparable. It is useful to define the influence radius r i as the location where they match  \nM BH r i = | U NBS ( r i ) | ⇔ ⇔ r i = M BH | U NBS ( r i ) | ∼ 2 M BH M Φ R 98 . (51)  \nSince M BH /lessmuch M NBS and the star is Newtonian, so R NBS /greatermuch M Φ , the transition radius satisfies 2 M BH /lessmuch r i /lessmuch R NBS , meaning it is located far-away from the BH horizon but still deep inside the NBS.", 'III.2.1. Inner region ( r /lessmuch r i )': "In this region, we take the test-field approximation and consider that the scalar field does not back react on the spacetime metric. The ordinary differential equation (ODE) for the mass function, Eq. (36), becomes simply  \nm ' ( r ) = 0 ⇔ m ( r ) = const = M BH , (52)  \nwhere we impose that at the BH horizon this constant is the BH mass.  \nHaving this, Eq. (37) reduces to  \nlog ' ( A ) = 2 M BH r 2 ( 1 -2 M BH r ) ⇔ ⇔ log( A ) = log ( 1 -2 M BH r ) +Υ ⇔ ⇔ A ( r ) = ( 1 -2 M BH r ) e Υ , (53)  \nwhere Υ is a constant which will be fixed by demanding continuity of A with the outer region, where the NBS self-gravity dominates over the BH potential.  \nSimilarly, the Klein-Gordon equation is solved by [14]  \nφ 0 ( r ) = C ( 1 -2 M BH r ) e -2 iµM BH (54)  \nwhere C is a constant, which will also be fixed by demanding continuity of φ 0 with the outer region.  \nIII.2.2. Outer region ( r /greatermuch r i )  \nWe now focus on the region at large radius, far-away from the radius of influence of the BH. Here, the scalar field is non-relativistic and the radial derivatives of φ 0 are subleading in Eqs. (36) and (37). Then, let us take the following ansatz for A and φ 0  \nA = ( 1 -2 M BH r ) e 2 U ( r ) , (55)  \nwhere U ( r ) ∼ O ( | φ 0 | 2 ) is a function to be determined, which satisfies | U ( r ) | /greatermuch M BH /r in this region. Plugging these in Eqs. (36)-(38) and keeping only the leading terms in 1 /r and | φ 0 | 2  \nφ 0 = ( 1 -2 M BH r ) -2 iµM BH ˜ φ 0 ( r ) , (56)  \nm ( r ) = M BH +4 πµ 2 ∫ r r i dr ' r ' 2 | φ 0 ( r ' ) | 2 , (57)  \n˜ φ '' 0 + 2 r ˜ φ ' 0 -2 µ ( µU + γ ) ˜ φ 0 = 0 . (59)  \nU ' = 4 πµ 2 r 2 ∫ r r i dr ' r ' 2 | φ 0 ( r ' ) | 2 , (58)  \nHere, we directly integrated the ODE for m ( r ) and imposed continuity with respect to mass function in the inner region obtained above (Eq. (52)).  \nDifferentiating the ODE for U with respect to r  \nU '' + 8 πµ 2 r 3 ∫ r r i dr ' r ' 2 | φ 0 ( r ' ) | 2 = 4 πµ 2 r 2 | φ 0 | 2 ⇔ ⇔ U '' + 2 r U ' = 4 πµ 2 r 2 | φ 0 | 2 , (60)  \nwhich we identify as the Poisson equation for the gravitational potential U being sourced by the scalar field distribution. Therefore, we identify U with the NBS potential U NBS described in the previous section. Since 2 M BH /lessmuch r i /lessmuch R NBS , for r ∼ r i the potential U NBS is already approximately constant. Therefore, demanding continuity of A ( r ) between the inner and outer region, we can match the constant Υ in the inner solution (53) with the central value of the NBS's gravitational potential, arriving at Υ = 2 U NBS (0).  \nFinally, we identify the Eq. (59) as the (timeindependent) Schrodinger equation for ˜ φ 0 , meaning it admits the same profile as the NBS of the previous section. Thus, similarly to the gravitational potential, in  \nthe matching region φ 0 is approximately constant and we can match the constant C of the inner solution with the central value of the scalar field, C = φ 0 (0).", 'III.3. Generic background for diluted configurations of matter around BHs': "The background spacetime obtained in the previous section for the case of parasitic BHs has similar properties to other diluted distributions of matter surrounding BHs. In particular, both for anisotropic fluids [4] and superradiant (scalar) clouds [17], the g tt metric function reduces to the Schwarzschild one with an exponential redshift pre-factor, as in Eq. (55).  \nWe now show this is a generic property when the compactness of the matter configuration is much smaller than the BH compactness. The tt and rr -component of Einstein's equations given, respectively, in Eq. (36) and Eq. (37) can be written for an arbitrary (sphericallysymmetric) energy-momentum tensor as [18]  \nm ' = 4 πr 2 ρ, (61)  \nwhere  \nlog ' ( A ) = 2 m r 2 ( 1 -2 m r ) + 4 πr 1 -2 m r p rad , (62)  \nρ = -T t t , p rad = T r r , (63)  \nare, respectively, a matter density and radial pressure term. For the case of scalar fields we identify  \nρ φ = 1 2 ( 1 -2 m r ) ( | φ ' 0 | 2 + µ 2 + ω 2 A 1 -2 m r | φ 0 | 2 ) ,  \n(64)  \np φ rad = 1 2 ( 1 -2 m r ) -1 ( | φ ' 0 | 2 -µ 2 -ω 2 A 1 -2 m r | φ 0 | 2 ) . (65)  \nWe now repeat the same steps as for the parasitic BH. In the inner region, we consider the matter density to be negligible, from which we directly obtain Eq. (52) and (53), and the BH solely controls the dynamics.  \nIn the outer region, we impose an ansatz for A as in Eq. (55) and substituting in Eq. (62) and expanding to first order in 1 /r  \nm ( r ) = M BH +4 π ∫ r r i dr ' r ' 2 ρ ( r ' ) , (66)  \nU ' = 4 π r 2 ∫ r r i dr ' r ' 2 ρ ( r ' ) + 4 πr p rad , (67)  \nwhere again imposed continuity with the mass function. For Newtonian configurations, the pressure term is negligible with respect to the matter density one. Differentiating with respect to r we arrive at  \nU '' + 2 r U ' = 4 π r 2 ρ, (68)  \nwhich once more is Poisson's equation for the Newtonian potential U of the generic matter distribtuion.", 'IV. FLUX FORMULAS AND ADDITIONAL RESULTS': 'As explained in the main text, the energy fluxes can be written in terms of the perturbations through their respective stress-energy tensors. In the case of gravitational perturbations, we use the Isaacson stress-energy tensor [7, 19]. The expressions are identical to the Schwarzschild case, as asymptotically the spacetimes are the same (with their respective masses). As such, we can just construct the Zerilli-Moncrief ψ Z function through H 1 and K [6] and compute the flux through  \n/lscript ˙ E g ∞ = 1 64 π ( /lscript +2)! ( /lscript -2)! | ˙ ψ Z | 2 , /lscript + m even , (69)  \nwhere the explicit form for ψ Z can be seen in Ref. [6]. Alternatively, we can also extract the amplitude through the asymptotic form of K , as  \nK ( r →∞ ) ≈ ∂ψ Z ∂r , (70)  \nand use this to compute the flux. For the scalar fluxes, we use the stress-energy tensor for the scalar field. The expression for this case can be read directly from Ref. [14, 17], being  \n/lscript ˙ E Φ ∞ =( ω + /lscript Ω p )Re [ √ ( ω + /lscript Ω p ) 2 -µ 2 ] × | φ + ( r →∞ ) | 2 , (71)  \nwhere φ + above is the Fourier amplitude of the scalar field perturbation defined in Eq. (5).  \nWe tested the above numerical relativistic expressions for BS by comparing the scalar fluxes with the analytical NBS presented in [14] (Eq. (149) there), accounting for relativistic modifications regarding the scalar field amplitude at the particle position and angular frequency  \nφ + ( r →∞ ) ≈ -8 π 3 2 m p φ 0 ( r p ) m m 2 -1 Y m m ( π 2 , 0 ) ( r p Ω) m 2 m 2 +2 (Ω /µ ) m 2 +1 Γ ( m + 3 2 ) . (72)  \nThe result can be seen in Fig. 1. The relativistic numerical results match the analytical one in the regime of validity predicted in the Newtonian approximation. Even going beyond it, the analytical prediction works remarkably well to reproduce the functional behavior of the flux.  \nIn the main text, we mention the energy flux going through the horizon in the resonant regime in the case of parasitic BHs. The same procedure used above can be transplanted to compute the fluxes going into the BH.  \nFIG. 1. Comparison between the NBS analytical expression (solid lines) with the full relativistic BS fluxes (dotted lines). In the regime of validity of the Newtonian approximation, the analytical expression matches very well our numerical the results. Even in more compact configurations, the analytical expression reproduces the correct behavior of the fluxes with the orbital radius.  \n<!-- image -->  \n- [1] Y. Kojima, Equations governing the nonradial oscillations of a slowly rotating relativistic star, Phys. Rev. D 46 , 4289 (1992).\n- [2] G. Allen, N. Andersson, K. D. Kokkotas, and B. F. Schutz, Gravitational waves from pulsating stars: Evolving the perturbation equations for a relativistic star, Phys. Rev. D 58 , 124012 (1998), arXiv:gr-qc/9704023.\n- [3] V. Cardoso, K. Destounis, F. Duque, R. P. Macedo, and A. Maselli, Black holes in galaxies: Environmental impact on gravitational-wave generation and propagation, Phys. Rev. D 105 , L061501 (2022), arXiv:2109.00005 [grqc].\n- [4] V. Cardoso, K. Destounis, F. Duque, R. Panosso Macedo, and A. Maselli, Gravitational Waves from Extreme-MassRatio Systems in Astrophysical Environments, Phys. Rev. Lett. 129 , 241103 (2022), arXiv:2210.01133 [gr-qc].\n- [5] F. J. Zerilli, Gravitational field of a particle falling in a schwarzschild geometry analyzed in tensor harmonics, Phys. Rev. D 2 , 2141 (1970).\n- [6] N. Sago, H. Nakano, and M. Sasaki, Gauge problem in the gravitational selfforce. 1. Harmonic gauge approach in the Schwarzschild background, Phys. Rev. D 67 , 104017 (2003), arXiv:gr-qc/0208060.\n- [7] K. Martel and E. Poisson, Gravitational perturbations of the Schwarzschild spacetime: A Practical covariant and gauge-invariant formalism, Phys. Rev. D 71 , 104003 (2005), arXiv:gr-qc/0502028.\n- [8] J. Barranco, A. Bernal, J. C. Degollado, A. DiezTejedor, M. Megevand, D. Nunez, and O. Sarbach, Selfgravitating black hole scalar wigs, Phys. Rev. D 96 , 024049 (2017), arXiv:1704.03450 [gr-qc].\n- [9] V. Cardoso, T. Ikeda, R. Vicente, and M. Zilh˜ao, Parasitic black holes: The swallowing of a fuzzy dark matter soliton, Phys. Rev. D 106 , L121302 (2022), arXiv:2207.09469 [gr-qc].\n- [10] DLMF, NIST Digital Library of Mathematical Func-  \nWith the amplitude of the scalar field at φ + ( r H ), the scalar flux can be obtained through  \n/lscript ˙ E Φ r H = ( ω + /lscript Ω p ) 2 | φ + ( r H ) | 2 . (73)  \nWith the above expression, we can compute the rate /lscript E Φ p , H explained in the main text. In Fig. 2 we show the quadrupole rate flux into the horizon as a function of the orbital radius. It is evident that for this configuration it is orders of magnitude smaller than the rate of change through the fluxes at infinity (see Fig. 1 in the main text). Nonetheless, it is interesting to see the resonant behavior of the fluxes going into the horizon, indicating that resonant configurations can enhance scalar accretion onto the central BH.  \nFIG. 2. Energy loss rate due to the quadrupolar scalar flux at the horizon for an EMRI evolving in the same NBS with M Φ µ = 0 . 23 hosting a parasitic BH of mass M BH /M Φ = 0 . 02 as in Fig. 1 of the main text. The peaks correspond to the resonant excitation of the proper modes of the scalar configuration, which produces a resonant accretion of scalar by the parasitic BH.  \n<!-- image -->  \ntions, http://dlmf.nist.gov/, Release 1.0.26 of 2020-0315, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.  \n- [11] T. Regge and J. A. Wheeler, Stability of a Schwarzschild singularity, Phys. Rev. 108 , 1063 (1957).\n- [12] S. L. Detweiler and E. Poisson, Low multipole contributions to the gravitational selfforce, Phys. Rev. D 69 , 084019 (2004), arXiv:gr-qc/0312010.\n- [13] V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Geodesic stability, Lyapunov exponents  \nand quasinormal modes, Phys. Rev. D 79 , 064016 (2009), arXiv:0812.1806 [hep-th].  \n- [14] L. Annulli, V. Cardoso, and R. Vicente, Response of ultralight dark matter to supermassive black holes and binaries, Phys. Rev. D 102 , 063022 (2020), arXiv:2009.00012 [gr-qc].\n- [15] R. Ruffini and J. A. Wheeler, Introducing the black hole, Phys. Today 24 , 30 (1971).\n- [16] J. D. Bekenstein, Transcendence of the law of baryonnumber conservation in black hole physics, Phys. Rev. Lett. 28 , 452 (1972).  \n- [17] 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- [18] C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Astrophysical signatures of boson stars: quasinormal modes and inspiral resonances, Phys. Rev. D 88 , 064046 (2013), arXiv:1307.4812 [gr-qc].\n- [19] R. A. Isaacson, Gravitational Radiation in the Limit of High Frequency. II. Nonlinear Terms and the Ef fective Stress Tensor, Phys. Rev. 166 , 1272 (1968).'}