Gravitational waves (GWs) provide a unique probe of compact-object mergers. While the LIGO–Virgo–KAGRA collaboration has now observed hundreds of such events [1]–[4], upcoming detectors will greatly enhance the precision of GW astronomy. Future space-based detectors such as the Laser Interferometer Space Antenna (LISA) [5], TianQin [6], and Taiji [7] will open the milliHertz window where extreme-mass-ratio inspirals (EMRIs) are among the most promising sources. EMRIs consist of a stellar-mass compact object inspiraling into a massive black hole (MBH), accumulating tens to hundreds of thousands of cycles in the strong-field regime. This makes them unparalleled probes of MBH spacetimes and their astrophysical environments [8]–[11].

A compelling scenario is the presence of an ultralight bosonic cloud around the primary black hole (BH), formed via superradiant instabilities [12]. Such environments can leave a detectable imprint on the GW signal emitted by an EMRI [13], [14], potentially revealing the existence of new fundamental bosons. Environmental effects of boson clouds in EMRIs have been extensively studied using Newtonian approximations [15]–[20], which are insufficient to capture the dynamics in the late inspiral phase, while relativistic calculations have only considered circular, equatorial EMRIs in scalar clouds [21]–[25]. However, realistic EMRIs are expected to possess non-vanishing eccentricity [26]. In this Letter, we reveal a new class of signatures for EMRIs embedded in scalar clouds, within the fully relativistic framework developed in our companion paper [27]. We show that eccentricity triggers a dense series of resonant transitions between cloud states in the strong-field regime, amplifying the imprint of the cloud’s presence in EMRI signals. These resonances arise from the split between azimuthal and radial frequencies in the strong-field regime, a feature that is absent in Newtonian approximations.

Related resonant phenomena were first identified as a mechanism for the floating or sinking of circular inspirals in scalar-tensor theories [28], and were subsequently shown to be excited across a discrete spectrum of harmonics in eccentric systems [29]. In those works, the resonances were sourced by an intrinsic scalar charge. Here we show that similar resonances can arise solely from the gravitational interaction between a standard relativistic perturber and a scalar cloud environment. While previous studies [15], [16], [19], [20], [30]–[32] also identified resonances in scalar clouds and discussed the associated orbital backreaction, their analyses were primarily Newtonian and focused on the early stages of the inspiral.
We consider a free complex scalar field \(\mathbf{\Phi}\) minimally coupled to gravity in the presence of a point-like perturber with mass \(m_p\), described by the action (we use geometrized units \(G=c=1\)) \[\require{physics} S=\int{\dd}^{4}x\sqrt{-\mathbf{g}}\left(\frac{\mathbf{R}}{16\pi}-\mathfrak{L}[\mathbf{\Phi}]\right)-m_{p}\int_{\gamma}{\dd}\tau,\] where \(\mathfrak{L}[\mathbf{\Phi}]=\partial_{\mu} \mathbf{\Phi}\partial^{\mu} \mathbf{\Phi}^*+\mu^2|\mathbf{\Phi}|^2\) is the scalar field’s Lagrangian density and \(\hbar \mu\) its mass; and \(\gamma\) is the point-particle’s worldline parametrized by its proper time \(\tau\). Varying the action with respect to the metric and the scalar field, one gets the Einstein-Klein-Gordon (EKG) field equations, \[\begin{align}\label{eq:32EKG32field32equation} \mathbf{G}_{\mu\nu} = 8\pi(T^{\mathbf{\Phi}}_{\mu\nu}+T^p_{\mu\nu}),\qquad \Box \mathbf{\Phi} = \mu^2 \mathbf{\Phi}, \end{align}\tag{1}\] where \(\mathbf{G}_{\mu\nu}:= \mathbf{R}_{\mu\nu}-\mathbf{g}_{\mu\nu}\mathbf{R}/2\) and \(\Box:=\mathbf{g}_{\mu\nu}\nabla^\mu\nabla^{\nu}\). The stress-energy tensor of the scalar field, \(T_{\mu\nu}^\mathbf{\Phi}\), reads \[\label{eq:32stress-energy32tensor32definition} T_{\mu\nu}^\mathbf{\Phi}=2\partial_{(\mu}\mathbf{\Phi}\partial_{\nu)}\mathbf{\Phi}^*-\mathbf{g}_{\mu\nu}\left(\partial_\alpha\mathbf{\Phi}\partial^\alpha\mathbf{\Phi}^*+\mu^2|\mathbf{\Phi}|^2\right),\tag{2}\] whereas the stress-energy tensor of the point-particle (describing a small compact object), \(T_{\mu\nu}^p\), is given by \[\require{physics} \label{eq:32stress-energy32tensor32of32the32particle} T^{p}_{\mu\nu}=m_p\int_{\gamma} u_\mu u_\nu\frac{\delta^{(4)}[x^\mu-\mathbf{x}_p^\mu(\tau)]}{\sqrt{-\mathbf{g}}}\dd\tau,\tag{3}\] with \(u^{\mu}=d\mathbf{x}_p^{\mu}/d\tau\) the particle’s four-velocity. This theory has a global \(U(1)\) symmetry with associated conserved Noether current and charge, respectively, \[\require{physics} \label{eq:32current} j^\mu=-i\left(\mathbf{\Phi}^*\partial^\mu\mathbf{\Phi}-\mathbf{\Phi}\partial^\mu\mathbf{\Phi}^*\right),\quad Q_\Sigma = \int_\Sigma \dd\Sigma_\mu \,j^\mu,\tag{4}\] where \(\Sigma\) is a spacelike hypersurface.

We model the system perturbatively through a two-parameter expansion controlled by the small mass-ratio \(q \equiv m_p/M \ll 1\) (with \(M\) the MBH mass) and the background scalar field amplitude \(\epsilon \ll 1\). The metric and scalar field are expanded as [21], [23] \[\mathbf{g}_{\mu\nu} = g_{\mu\nu} + \sum_{i,j} \epsilon^i q^j h_{\mu\nu}^{(i,j)},\quad \mathbf{\Phi} = \sum_{i,j} \epsilon^i q^j \phi^{(i,j)},\] where here \(g_{\mu\nu}\) is the Schwarzschild background \(\require{physics} \dd s^2 = -f \dd t^2 + f^{-1} \dd r^2 + r^2 \dd \Omega^2\) with \(f=1-2M/r\). Substituting these expansions into the EKG equations yields a hierarchy of equations, which are solved order-by-order.

At order \(\mathcal{O}(\epsilon^1 q^0)\) the metric is just Schwarzschild and the scalar field satisfies the homogeneous Klein-Gordon (KG) equation \[(\square^{(0)} - \mu^2)\phi^{(1,0)} = 0\,,\] which admits quasi-bound state solutions. Exploiting spherical symmetry, we decompose the field as \(\phi^{(1,0)} = R^{b}_{n_i\ell_i}(r)Y_{\ell_im_i}(\theta, \varphi)e^{-i\omega t}\) with indices \(\{n_i, \ell_i, m_i\}\) characterizing the solutions. Here, \(Y_{\ell_im_i}\) are scalar spherical harmonics. Quasi-bound state solutions with radial function \(R^{b}_{n_i\ell_i}\) and their corresponding complex eigenfrequencies \(\omega:=\omega_{n_i\ell_i}\) are found by imposing exponential decay at spatial infinity and purely ingoing waves at the event horizon and using the continued-fraction method developed in Refs. [33], [34]. Their total mass is \(M_b=\omega\, Q[\phi^{(1,0)},g]\); for concreteness, we will take \(\epsilon\equiv \sqrt{M_b/M}\).

While superradiant growth of quasi-bound states requires a spinning BH [34]–[38], for light scalars (\(M\mu \ll 1\)) superradiant extraction drives the BH toward a small saturation spin, and the scalar cloud is localized far from the horizon, where spin effects are subleading. Therefore, the Schwarzschild metric provides a good approximation for the cloud profile and associated radiation fluxes [21], [23].

At order \(\mathcal{O}(\epsilon^0 q^1)\) the metric perturbation \(h_{\mu\nu}^{(0,1)}\) is governed by the linearized Einstein equations \[\delta \mathbf{G}_{\mu\nu}[h^{(0,1)}] = 8\pi T_{\mu\nu}^{p}[g].\] We solve these equations within the Regge-Wheeler-Zerilli formalism, with metric reconstruction performed according to the procedure in Refs. [39], [40]. Metric perturbations are sourced by a small compact object following stable bound timelike geodesics in the Schwarzschild geometry.

In the equatorial plane, geodesics are determined by the dimensionless semi-latus rectum \(p\) and eccentricity \(e\), which are related to the conserved specific energy \(\mathcal{E}\) and angular momentum \(\mathcal{L}\) via [41], [42] \[\mathcal{E}^2=\frac{(p-2-2e)(p-2+2e)}{p(p-3-e^2)},\quad\mathcal{L}^2=\frac{p^2M^2}{p-3-e^2}.\] The requirement of a stable bound orbit implies \(p \geq p_{\rm LSO}\equiv 6+2e\), with the equality at the last stable orbit (LSO); for smaller values the orbit is unstable and the particle plunges. The source’s eccentric motion restricts the metric perturbation to a discrete frequency spectrum \[\sigma_{mn} = m\,\Omega_\varphi + n\,\Omega_r, \qquad m,n\in\mathbb{Z},\] with azimuthal and radial orbital frequencies, respectively, \(\Omega_\varphi(p,e)\) and \(\Omega_r(p,e)\) [41].

The influence of the point particle on the scalar field first appears at \(\mathcal{O}(\epsilon^1 q^1)\) through the inhomogeneous KG equation, \[\label{eq:1195KG} (\square^{(0)} - \mu^2)\phi^{(1,1)} = \nabla_{\mu}^{(0)}\bar{h}_{(0,1)}^{\mu\nu}\partial_{\nu}\phi^{(1,0)} + h_{(0,1)}^{\mu\nu}\nabla_{\mu}^{(0)}\partial_{\nu}\phi^{(1,0)},\tag{5}\] where \(\bar{h}_{\mu\nu}^{(0,1)} = h_{\mu\nu}^{(0,1)} - \frac{1}{2}g_{\mu\nu}h^{(0,1)}\) is the trace-reversed metric perturbation. We solve for the scalar perturbation using a Green’s function method, imposing outgoing boundary conditions at spatial infinity and ingoing ones at the MBH horizon; further details are given in [27].
The perturbative framework just described implies that, in addition to the usual GW fluxes, the EMRI exchanges energy and angular momentum with the ambient scalar field through gravitational scattering. In a steady-state regime, this exchange is encoded in the fluxes of energy and angular momentum carried by the scalar perturbations to spatial infinity and through the MBH horizon [43]–[45].

The orbit-averaged scalar fluxes are nonvanishing at order \(\mathcal{O}(\epsilon^{2}q^{2})\), since they arise from quadratic terms in \(\phi^{(1,1)}\). The energy fluxes are defined with respect to the stationary Killing vector \(\xi^\mu_{(t)}=(\partial_t)^\mu\) as [21], [23] \[\require{physics} \label{eq:dotE95Phi} \begin{align} \dot{E}^{\Phi,\infty} &= -\lim_{r\to +\infty} r^2\int \dd\Omega\, T^{\phi^{(1,1)}}_{\mu r}\xi_{(t)}^{\mu}\,,\\ \dot{E}^{\Phi,H} &= \lim_{r\to 2M} 4M^2 \int \dd\Omega\, T^{\phi^{(1,1)}}_{\mu t}\xi_{(t)}^{\mu}\,. \end{align}\tag{6}\] The angular momentum fluxes \(\dot{L}^{\Phi,\infty/H}\) are obtained analogously by replacing \(\xi_{(t)}^\mu\) with the axial Killing vector \(\xi_{(\varphi)}^\mu=(\partial_\varphi)^\mu\). The charge evolves due to analogous fluxes to infinity and through the horizon, \(\dot{Q}^{\infty/H}\), which are also order \(\mathcal{O}(\epsilon^{2}q^{2})\). These induce changes in the cloud mass and angular momentum as \(\dot{M}_b = \omega\, \dot{Q}\) and \(\dot{J}_b = m_i\, \dot{Q}\).

Assuming that the inspiral proceeds on a timescale much longer than that to reach a steady-state energy-momentum exchange (i.e., using the adiabatic approximation), the point-particle’s specific energy \(\mathcal{E}\) and angular momentum \(\mathcal{L}\) evolve according to: \[\begin{align}\label{eq:balance95laws} &\dot{\mathcal{E}}=-m_p^{-1} \bigl( \dot{E}^{g,\infty}+\dot{E}^{g,H}+\dot{E}^{s,\infty}+\dot{E}^{s,H}\bigr),\\ &\dot{\mathcal{L}}=-m_p^{-1} \bigl(\dot{L}^{g,\infty}+\dot{L}^{g,H}+\dot{L}^{s,\infty}+\dot{L}^{s,H}\bigr), \end{align}\tag{7}\] where \(\dot{E}^{g,\infty/H}\) and \(\dot{L}^{g,\infty/H}\) are the GW energy and angular momentum fluxes at infinity and at the horizon, and we defined \(\dot{E}^{s,\infty/H}\:=\:\dot{E}^{\Phi,\infty/H}+\omega\:\dot{Q}^{\infty/H}\) and \(\dot{L}^{s,\infty/H}\:=\:\dot{L}^{\Phi,\infty/H}+m_i\:\dot{Q}^{\infty/H}\). These dissipation rates are computed in full generality (including eccentricity and orbital inclination) for a Schwarzschild background in the companion paper [27]. Here we focus on equatorial, prograde orbits. We neglect the conservative effects of the cloud’s self-gravity on \(\mathcal{E}\) and \(\mathcal{L}\), which are non-vanishing at \(\mathcal{O}(\epsilon^2)\), and \(\mathcal{O}(\epsilon^2q^2)\) corrections to the GW fluxes [21], [25]. While these terms have not been computed yet for dipolar clouds, they are expected to be approximately degenerate with a redshift at small \(M\mu\) (low compactness) or with a mass shift at large \(M\mu\) (high compactness) [22], [25]. Their relative importance at intermediate \(M\mu\) for dipolar clouds is unclear, but they should not impact the existence of the resonances studied here. Eq. 7 also assumes that the orbit can be evolved only from radiative fluxes in the steady-state regime. An extension to this picture is to include non-radiative (local) torques using the torque-balance law of Ref. [46].

Here, our emphasis is on the role of resonances as characteristic signatures of scalar clouds in the signals of EMRIs. These resonances occur at an infinite set of orbital frequencies (that we label with \(n_f\)) given by \[\label{eq:resonances} \sigma^{n_f}_{mn}=\Re(\omega_{n_f \ell_f}) - \Re(\omega_{n_i\ell_i})\,,\tag{8}\] where we recall that \(\omega_{n \ell}\) corresponds to an eigenfrequency of the KG equation with quantum numbers \(\{n,\ell\}\) (note that in a Kerr background the eigenfrequencies would also depend on \(m\)). Eccentric EMRIs can generically excite such resonances in the strong-field regime, where space-based detectors will probe these systems.

Fig. 1 compares the steady-state dissipation rates from gravitational radiation and from scalar-field scattering as a function of the semi-latus rectum, for an EMRI with \(e=0.6\) in a prograde dipolar cloud \(\{n_i,\ell_i,m_i\}=\{0,1,1\}\) with mass \(M_b/M=0.05\) and coupling \(M\mu=0.2\) (for an MBH of mass \(M=10^6\,M_\odot\), this corresponds to bosons of mass \(\hbar\mu\sim 3\times 10^{-17}\,\mathrm{eV}\)). For comparison, we also show the horizon loss rates for circular geodesics (\(e=0\)) and for \(e=0.8\). The total GW energy and angular-momentum fluxes are computed with the FastEMRIWaveforms package [47]–[49]. The scalar loss rates are dominated by the horizon terms \(\dot{L}^{s,H}\) and \(\dot{E}^{s,H}\), as in the circular case [21], [23], [24]. The step-like behavior seen in \(\dot{L}^{s,\infty}\) is also expected and was first explained within Newtonian approximations [18], [50].

Most strikingly, a series of resonances that grow in importance with increasing eccentricity and are absent in the circular case, can be clearly seen in the scalar-field torque. All these resonances arise from transitions to cloud states with \(\ell_f=m_f=0\), occurring at \(\sigma_{-1, n}^{n_f}\). The most prominent broad peaks, especially noticeable for \(e=0.8\), correspond to resonances with the \(n_f=0\) state, in the \(n =1, 2, 3, 4\) harmonics (ordered by decreasing \(p\)). These broad resonances substantially enhance the torques and orbital loss rates, as can be seen by the nearly two orders of magnitude difference between the on-resonance \(\dot{L}^{s,H}\) for \(e=0.8\) and for the circular case.

In the torque \(\dot{L}^{s,H}\), a clear substructure consisting of a series of narrower peaks can also be seen. These secondary peaks arise from resonances in the subsequent harmonic \(n+1\) with \(n_f \geq 1\) states, so that close to the peak associated with the \((n,n_f=0)\) resonance, there is always a sequence of narrower peaks arising from \((n+1,n_f \geq 1)\) resonances. Interestingly, the accumulation of resonances occurs to the left of the position where \(\dot{L}^{s,\infty}\) has a sharp jump, which is the opposite behavior of what happens for resonances that occur at larger \(p\) [32].

The eccentricity-induced resonances discussed here are absent in a Newtonian treatment, as they arise from \(\Omega_{\varphi} \neq \Omega_{r}\) in the strong-field regime. In the weak-field regime, \(\Omega_{\varphi} \approx \Omega_{r}\), and increasing \(n\) instead produces resonances at larger \(p\), opposite to the behavior described here.
To assess the impact of eccentricity and resonances in the EMRI dynamics, we model the orbital evolution under the adiabatic approximation. The trajectory in the \(\{p, e\}\) parameter space is driven by the losses in specific energy \(\mathcal{E}\) and angular momentum \(\mathcal{L}\), including both gravitational and scalar channels. This evolution is governed by \[\require{physics} \begin{align} &\frac{\dd p}{\dd t} = \frac{1}{H} \left( \frac{\partial \mathcal{L}}{\partial e} \frac{\dd\mathcal{E}}{\dd t} - \frac{\partial \mathcal{E}}{\partial e} \frac{\dd\mathcal{L}}{\dd t} \right), \\ &\frac{\dd e}{\dd t} = \frac{1}{H} \left( \frac{\partial \mathcal{E}}{\partial p} \frac{\dd\mathcal{L}}{\dd t} - \frac{\partial \mathcal{L}}{\partial p} \frac{\dd\mathcal{E}}{\dd t} \right), \end{align}\] where \(H\) is the Jacobian of the transformation from \((\mathcal{E}, \mathcal{L})\) to \((p, e)\), given by \(H = \frac{\partial \mathcal{E}}{\partial p} \frac{\partial \mathcal{L}}{\partial e} - \frac{\partial \mathcal{L}}{\partial p} \frac{\partial \mathcal{E}}{\partial e}\).

From \(p(t)\) and \(e(t)\) we can compute the orbital frequencies \(\Omega_{\varphi,r}(t)\) and define the accumulated phase \(\require{physics} \Psi_{mn}(t) = \int^{t}_0 \sigma_{mn}(t')\dd t'\). To quantify the environmental impact of the scalar cloud, we compute the dephasing relative to a vacuum inspiral \[\require{physics} \Delta\Psi_{mn} = m\Delta\Psi_{\varphi}+n\Delta\Psi_{r},\quad \Delta\Psi_{\varphi,r} = \int \Delta\Omega_{\varphi,r}\dd t\, .\] We find that \(\Delta\Psi_{r}\) is typically slightly larger but comparable to \(\Delta\Psi_{\varphi}\) in magnitude. Moreover the mode with the largest amplitude in the gravitational waveform is typically the \(\{m,n\}=\{2,0\}\) mode [51]. We therefore focus on the azimuthal dephasing \(\Delta\Psi_{\varphi}\), since it provides a conservative measure of the environmental impact on the observed signal.

We consider systems with initial dimensionless semi-latus rectum \(p_i=10\) and eccentricities \(e_i \in \{0.2,~0.4,~0.6,~0.8\}\) immersed in a dipolar scalar cloud with the same properties as in Fig. 1. For these four systems, we show the orbital evolution trajectories in the \((p,e)\) plane as well as the corresponding dephasing in Fig. 2. The infinity loss rates (\(\dot{E}^{s\,,\infty}\), \(\dot{L}^{s\,,\infty}\)) have only a modest effect on the orbital evolution, since they are smaller than the horizon contributions across the parameter space considered here (see Fig. 1). We compute the infinity contributions using our fully relativistic framework for all evolutions except the \(e_i=0.8\) case, where the large number of radial harmonics required for convergence makes the computation costly; for this case we instead use the Newtonian approximation described in the companion paper [27]. We checked that for \(e_i=0.6\) this approximation changes the accumulated dephasing by only \(\sim4\%\), so we expect its use at \(e_i=0.8\) to have a similarly small impact.

For \(e_i=0.8\), there is a clear departure in eccentricity evolution from the vacuum case, with the eccentricity decreasing more slowly in the presence of the scalar cloud, as shown in the left panel of Fig. 2. This holds for all eccentricities, but the effect becomes more pronounced as \(e_i\) increases, owing to the enhancement of the scalar loss rates (mainly induced by the resonances). For a typical EMRI with \(q=10^{-5}\) and \(M=10^6\,M_\odot\), the total signal duration from \(p_i=10\) to the LSO is approximately one and a half years. Over this period, the accumulated dephasing \(\Delta\Psi_\varphi\) relative to vacuum reaches \(\sim10^4\) rad for \(e_i=0.8\) and \(\sim10^3\) rad for \(e_i=0.6\), while for smaller eccentricities (\(e_i=0.2\)) we obtain \(\sim 90\) rad, compared to \(\sim 70\) rad in the circular limit. Even after half a year of inspiral, the accumulated dephasing for \(e_i=0.8\) already exceeds the \(e_i=0.2\) case by nearly two orders of magnitude. This indicates that, while eccentricity does little to improve the detectability of the cloud for eccentricities as small as \(e_i=0.2\), the improvement becomes substantial at higher eccentricities.

The accumulated dephasing scales approximately linearly with the cloud mass \(M_b\). For a lighter cloud with \(M_b/M=5\times10^{-3}\) and the same initial semi-latus rectum \(p_i=10\), the total dephasing over the entire inspiral is \(\sim1.5\times10^3\) rad for \(e_i=0.8\) and \(\sim 9\) rad for \(e_i=0.2\), compared to \(\sim 7\) rad in the circular limit. This suggests that the effect should remain detectable for typical LISA signal-to-noise ratios [14], [52], even for small cloud masses.
In this Letter, we demonstrated that orbital eccentricity plays a key role in the relativistic evolution of EMRIs in scalar-cloud environments. Our most striking finding is the emergence of resonances in the strong-field regime, absent in both circular-orbit models and Newtonian approximations. Resonances substantially enhance the energy and angular-momentum exchange between the EMRI and the scalar cloud, an effect that becomes increasingly pronounced at larger eccentricities. As a result, eccentricity can greatly improve the detectability of scalar clouds, underscoring the importance of eccentricity in both the modeling and interpretation of EMRI signals in boson-cloud environments.

The existence of eccentric EMRIs is well-supported by standard EMRI formation mechanisms via two-body relaxation [26], [53], [54]. These mechanisms predict that EMRIs are injected directly into the strong-field region with high eccentricity and a small initial semi-latus rectum (\(p\sim 20\) [55]), which could help to bypass resonances encountered at larger separations and avoid long inspirals that could dissipate the cloud before the system enters the LISA band [19], [20], [31].

A natural extension of our work would be to consider a Kerr MBH. Although superradiant growth may spin the BH down to small values, the tri-periodic nature of generic Kerr orbits could still enrich further the resonant signatures identified here.
We are grateful to Seth Hopper for kindly sharing his code developed in Ref. [39] with us. We are grateful to Conor Dyson, Robrecht Keijzer and Thomas Spieksma for insightful discussions and comments on the draft of this manuscript. We acknowledge financial support provided by FCT – Fundação para a Ciência e a Tecnologia, I.P., through the ERC-Portugal program Project “GravNewFields”. We also thank the Fundação para a Ciência e Tecnologia (FCT), Portugal, for the financial support to the Center for Astrophysics and Gravitation (CENTRA/IST/ULisboa) through grant No. UID/PRR/00099/2025 and grant No. UID/00099/2025. QX gratefully acknowledges support from FCT grant 2025.01546.BD. RV gratefully acknowledges the support of the Dutch Research Council (NWO) through an Open Competition Domain Science-M grant, project number OCENW.M.21.375.

References↩︎