Environmentally-induced chaos: Extreme-mass-ratio systems of rotating black holes in astrophysical environments

TL;DR

This paper investigates how astrophysical environments around rotating black holes induce chaos in the orbital dynamics of extreme-mass-ratio inspirals, affecting gravitational-wave signals and challenging existing models.

Contribution

It introduces the first analysis of geodesics around rotating black holes with matter halos, revealing environmental effects cause non-integrability and chaos in orbital motion.

Findings

Environmental effects break Kerr spacetime symmetries.

Chaos manifests through resonant islands and chaotic layers.

Resonances can be significantly prolonged, impacting gravitational-wave signals.

Abstract

Extreme-mass-ratio inspirals, in which a stellar-mass object orbits a supermassive black hole, are prime sources of millihertz gravitational waves for upcoming space-based detectors. While most studies assume idealized vacuum backgrounds, realistic extreme-mass-ratio binaries are embedded in astrophysical environments such as accretion disks, stellar clusters, or dark matter spikes, disks, and halos, which can significantly alter the orbital dynamics. We explore bound geodesics around general-relativistic solutions describing rotating black holes surrounded by matter halos for the first time, mapping how environmental effects interfere with the spacetime symmetries of vacuum spinning (Kerr) black holes. In particular, we find that the loss of a Carter-like constant leads to geodesic non-integrability and the onset of chaos. This manifests through the formation of resonant islands and chaotic layers around transient orbital resonances in phase space--features that are otherwise completely absent in integrable Kerr geodesics. Resonant islands, which are extended, non-zero volume regions in phase space, encapsulate periodic orbit points. Non-integrability dictates that all geodesics inside the resonant island share the periodicity of the resonance. Thus, the lifespan of resonances around non-Kerr objects can be significantly enhanced beyond the predicted lifetime of Kerr resonances. Consequently, these effects can leave distinct imprints on gravitational-wave signals, with significant implications for gravitational-wave modeling and parameter inference of astrophysical extreme-mass-ratio inspirals.

Tables

Table 1. Table 1: Initial conditions used in the code to obtain all solutions used throughout this work. In all cases, M BH = 1 M_{\rm BH}=1 .

$J/M_{\mathrm{BH}}^2$	$r_h$	${\Omega }_h$	$ℳ$	$a_0$
0.17514520	0.49649779	0.030211615	9.9998521	99.998521
0.30021574	0.48609324	0.052252992	9.9992814	99.992814
0.44286783	0.46705390	0.078376538	9.9976902	99.976902
0.63465392	0.42787085	0.11637563	9.9924688	100.00646
0.79463800	0.34784164	0.19207557	9.9936337	201.33424
0.80351794	0.36333864	0.17138511	10.005515	130.07169
0.80479191	0.37707168	0.15483751	9.9837236	100.05776
0.99565384	0.29486837	0.20766210	9.9991003	99.991003

Equations

$$T_{\nu }^{\mu }=\text{diag}(-\epsilon ,0,P_t,P_t),$$

$$\epsilon =\frac{M_{\mathrm{halo}}{a}_0}{2\pi r(r+a_0{)}^3},$$

$$m(r)=M_{\text{BH}}+\frac{M_{\mathrm{halo}}{r}^2}{(a_0+r{)}^2}{\left(1-\frac{2M_{\text{BH}}}{r}\right)}^2,$$

$$d{s}^2=-f(r)d{t}^2+\frac{d{r}^2}{1-2m(r)/r}+r^{2}d{\Omega }^2,$$

$$f(r)$$

$${\rm Y}$$

$$+2\sqrt{M_{\mathrm{halo}}/\xi }\arctan \left[\frac{r+a_0-M_{\mathrm{halo}}}{\sqrt{M_{\mathrm{halo}}\xi }}\right],$$

$$\xi$$

$$T_{\mu \nu }=$$

$$+(p_r-p_1)k_{\mu }{k}_{\nu }+(p_2-p_1)s_{\mu }{s}_{\nu },$$

$$d{s}^2=-f\frac{N_{-}^2}{N_{+}^2}d{t}^2+$$

$$+r^2{\sin }^2\theta {\left(d\phi -\omega dt\right)}^2],$$

$$u^{\mu }=\frac{1}{\sqrt{-(g_{tt}+2\Omega g_{t\phi }+{\Omega }^2{g}_{\phi \phi })}}(1,0,0,\Omega ),$$

$$\kappa =\frac{f}{8r_h\sqrt{gh}},A_h=32\pi r_h^2{\int }_0^{\pi }\frac{g\sqrt{h}}{f}\sin \theta \text{differential}\theta ,$$

$$\epsilon =\frac{\mathcal{M}\left(a_0+r_h\right)}{2\pi r(r+a_0{)}^3}{\left(1-\frac{r_h}{r}\right)}^2{b}^{-5},$$

$$d{s}^2=g_{tt}d{t}^2+2g_{t\phi }dtd\phi +g_{rr}d{r}^2+g_{\theta \theta }d{\theta }^2+g_{\phi \phi }d{\phi }^2,$$

$${\ddot{x}}^{\kappa }+{\Gamma }_{\lambda \nu }^{\kappa }{\dot{x}}^{\lambda }{\dot{x}}^{\nu }=0,$$

$$-E/\mu =g_{tt}\dot{t}+g_{t\phi }\dot{\phi },L_z/\mu =g_{t\phi }\dot{t}+g_{\phi \phi }\dot{\phi }.$$

$$\dot{t}=\frac{E{g}_{\phi \phi }+L_z{g}_{t\phi }}{\mu \left(g_{t\phi }^2-g_{tt}{g}_{\phi \phi }\right)},\dot{\phi }=\frac{E{g}_{t\phi }+L_z{g}_{tt}}{\mu \left(g_{tt}{g}_{\phi \phi }-g_{t\phi }^2\right)}.$$

$${\dot{r}}^2+\frac{g_{\theta \theta }}{g_{rr}}{\dot{\theta }}^2+V_{\text{eff}}=0,$$

$$V_{\text{eff}}\equiv \frac{1}{g_{rr}}\left(1+\frac{g_{\phi \phi }{E}^2+g_{tt}{L}_z^2+2g_{t\phi }E{L}_z}{{\mu }^2\left(g_{tt}{g}_{\phi \phi }-g_{t\phi }^2\right)}\right).$$

$$g_{tt}(r,\theta )$$

$$+\frac{(1+r_h/r{)}^{4}g(r,\theta ){\sin }^2\theta {\omega }^2(r,\theta )}{r^{2}f(r,\theta )},$$

$$g_{rr}(r,\theta )$$

$$g_{\theta \theta }(r,\theta )$$

$$g_{\phi \phi }(r,\theta )$$

$$g_{t\phi }(r,\theta )$$

$${\nu }_{\vartheta }=\frac{1}{2\pi \mathcal{N}}\sum _{i=1}^{\mathcal{N}}{\vartheta }_i.$$

Full text

I Introduction

As gravitational-wave (GW) detections from ground-based interferometers are increasing exponentially [1, 2], our understanding of gravitation becomes more precise [3]. In the verge of an era where GW detections are systematically becoming far more sensitive and accurate than the initial, though still ground-breaking, detection of 2015 [4], novel space-borne detectors, like the Laser Interferometer Space Antenna (LISA) [5], will unlock a new window of the Cosmos through the detection of gravitational radiation emitted by enticing theorized sources, thus transforming our understanding of GW astrophysics and astronomy [6, 7].

The main categories of LISA targets that regard fundamental physics [8, 9, 10] and astrophysics [11, 12] alone are enthralling; supermassive black-hole (BH) mergers, such as those occurring during galactic collisions that include such cosmic giants in their very core [13, 14], compact binary systems in the Milky Way, and infalling stellar-mass compact objects into supermassive BHs, known as extreme-mass-ratio inspirals (EMRIs) [15, 16, 17]. LISA will also be sensitive to the measurement of the stochastic GW background which is very important for gravitation, cosmology and particle physics. LISA will contribute to the understanding of the early Universe, the nature of dark matter and dark energy, as well as potential deviations from standard cosmological models and “beyond the Standard Model” physics [18].

EMRIs are one of the most promising targets of LISA due to their expected formation rates in galactic centers [19, 20, 21, 22], and the fact that the secondary stellar-mass companion is anticipated to perform thousands of revolutions in the strong-field regime before plunging into the supermassive primary. Their extremely large lifetime instantly renders numerical relativity impractical. New techniques are being developed, in order to speed up the process of waveform modeling in EMRIs [23, 24, 25], though there is still an abundance of phenomenology that has to be taken into account. Their unique cosmic design is tailored for precision tests of General Relativity (GR) [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40], probing the existence of new fundamental fields [41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51], and even testing the existence of astrophysical environments around EMRIs [52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79].

Even though, significant work has been done under the assumption of a vacuum spacetime, the majority of binary merger events should, in principle, take place inside astrophysical environments that occupy our Universe [80, 81, 82]. In fact, astrophysical environments pose a key feature in various contexts of gravitational radiation; from the inspiral till the ringdown [83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100]. Realistic EMRIs, are, thus, expected to interact with surrounding matter. These environmental components can alter the orbital structure of EMRIs. We now have more than enough evidence that matter environments, luminous or otherwise, are present in the surroundings of compact objects [101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114]. Therefore, a huge effort is being spent into employing the addition of environmental effects in GW generation and propagation from EMRIs and systematically streamlining the process.

The bulk of EMRI analyses are treating environments within perturbation theory or with Newtonian and post-Newtonian schemes that approximate the self-gravity of matter, estimate dynamical friction, and simplify gravitational redshift or peculiar motion [115, 116, 117, 118, 119, 120]. Even though successful at the time, due to the demands of LISA parameter inference a growing need for first-principle calculations of BHs in astrophysical environments is inevitable. To go beyond crude estimates, a general-relativistic, exact solution of the field equations has been obtained [121], that describes a BH immersed in a Hernquist-type dark matter halo [122]. This technique helped the community to extend further the construction of general-relativistic spacetimes to other exact and numerical BH solutions with different dark matter distributions around them [123, 124, 125, 126]. The axial and polar GW fluxes of circular equatorial EMRIs including general-relativistic, static and spherically-symmetric BHs surrounded by generic environments have been investigated in a relativistic setup in Refs. [127, 128, 129, 130]. The results strongly suggest the need for a paradigm-shifting strategy in order to accurately capture environmental effects in EMRIs.

A milestone regarding BHs in environments was achieved recently; a fully-relativistic numerical configuration of a rotating BH embedded in a Hernquist-type matter halo [131]. These solutions, and the analysis of EMRIs that will be performed with them, will lead us one step closer to the ultimate goal; a generic perturbation scheme for general-relativistic BHs in astrophysical environments that calculates the GW fluxes of generic EMRIs. This technique will provide unlimited access to precise EMRI waveform models, and to date it only exists in spherical symmetry [127, 100]. EMRIs evolve adiabatically due to the large mass disparity between the primary and secondary. This allows one to perceive the evolution of the secondary, at leading-order, as a geodesic of a massive test-particle in a fixed background. At next-to-leading order, radiation reaction effects take place and drive the secondary through a drift of successively-damped geodesics.

Formally, a general, stationary and axisymmetric spacetime possesses two Killing vector fields associated with stationarity and axisymmetry. These Killing vector fields are associated with two constants (integrals) of motion; the conserved energy and azimuthal angular momentum. Vacuum Kerr spacetimes have the privilege of possessing yet another constant of motion, namely the Carter constant [132], which is associated to a rank-two Killing tensor field. The conservation of the energy, azimuthal angular momentum and Carter constant of the test particle, are enough to deem the geodesics separable, i.e., integrable. This results to four, decoupled, first-order ordinary differential equations for the degrees of freedom of the metric [133].

On the other hand, if the, very fragile, Carter symmetry is broken, then even though the temporal and azimuthal geodesics still form decoupled and first-order differential equations, the radial and polar motion of geodesics consist of a coupled second-order dynamical system. In this case, the geodesics of spacetime are non-integrable, i.e. they lack a separation constant, and chaotic phenomenology in the orbital phase space emanates [133]. These phenomena appear in many distinctive ways in various dynamical systems. The most relevant indicator of chaos in EMRIs is the formation of a Birkhoff chain around resonant points in a Poincaré surface of section. The chain includes nested island formations around the periodic stable points, also known as resonant islands or islands of stability.

Resonances in Kerr EMRIs are by themselves challenging for GW modeling and data analysis [134]. This is due to the fact that when the secondary crosses a Kerr resonance, it undergoes a transient phase from a quasi-periodic to periodic orbit and back that can last up to dozens of cycles [135, 136, 137, 138]. If the geodesics are non-integrable, the resonant islands that surround stable periodic points become regions of shared periodicity for all geodesics that occupy the island. This inevitably leads to resonances that last much longer [139, 140, 141, 142] ($\sim$ hundreds of orbits of purely resonant motion) and designates the existence of non-integrability and indirect chaos.

Astrophysical chaos has proven to be a prominent effect in non-Kerr EMRIs, that can significantly affect the GW frequency evolution when a resonant island is crossed. The studies conducted so far have employed theory-agnostic primaries that are tailored to break integrability [143, 144, 145], pathological solutions of GR [146, 147, 148], exotic compact object primaries [149, 150], and even EMRI analogs [151, 152]. Even though the above classes of EMRIs discussed might lead to chaotic phenomena both at the orbital and GW level, they are either exotic, arise from modifications of gravity or include secondary effects that modify the geodesic equations, such as the secondary spin [153] and other relativistic effects [154, 155].

In this work, we embark into a tentative geodesic treatment of general-relativistic solutions of the field equations, describing rotating BHs surrounded by a Hernquist-type dark matter distribution [131], in an attempt to examine if there are imprints of non-integrability when, relativistically-constructed, rotating BHs in matter environments are present111Very early attempts for unraveling chaotic phenomena in Schwarzschild BHs perturbed with matter halos, that are introduced in the form of dipolar, quadrupolar and octupolar perturbations to background have been successful, though the contribution of the matter halos was added by hand to vacuum GR backgrounds and introduced various pathologies [156, 157].. This analysis results to a novel class of chaos that is environmentally-driven and purely relativistic. Thus here, we lay the groundwork for future EMRI studies, that will include radiation reaction, and provide powerful insights into how realistic astrophysical environments influence spacetime symmetries.

II Black holes in astrophysical environments

In what follows, we review the first exact, fully-relativistic, solution to the field equations that describes a static and spherically-symmetric BH in the center of a matter halo environment222Of course, other general-relativistic solutions exist but they do not include a halo of matter but rather a thin disk around the BH, see e.g. [158, 159, 160]. [121]. Afterwards, we revisit the generalization of this exact solution to a general-relativistic and fully-numerical spacetime that describes a stationary and axially-symmetric rotating BH surrounded by a matter environment [131].

II.1 Exact, static and spherically-symmetric BHs in matter halos

The first, exact solution of Einstein’s equations that describes a static and spherically-symmetric BH placed in the core of a dark matter halo was obtained in [121]. The basic principle of the construction of the environment is the assumption of many gravitating masses, that follow all possible spherical geodesics around a central point where the BH exists, thus building an Einstein cluster [161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 125, 176]. By integrating all geodesics over a finite radius range leads to introducing an anisotropic fluid with zero radial pressure (since the Einstein cluster assumes geodesics) and non-zero tangential pressure $P_{t}$ that drives the fluid, such that

[TABLE]

where $\varepsilon$ is the energy density of the matter profile. Even though in [121] a density profile was used that describes dark matter halos around galaxies [177, 178, 179, 180], namely the Hernquist density profile [122]

[TABLE]

where $M_{\rm halo}$ is the mass of the halo, $a_{0}$ its length scale and $M_{\rm halo}/a_{0}$ defines the halo compactness, there have been other configurations that have been constructed in an equivalent manner though they do not lead to an exact, but rather a numerical solution. The assumption of spherical symmetry, together with a Hernquist-inspired matter distribution

[TABLE]

with $M_{\text{BH}}$ the mass of the primary central BH, leads to the geometry333Recently, it was shown in Ref. [176] that this, and other geometries following from the Einstein cluster construction, are also solutions to a particular vector-tensor theory, where the vector encodes the astrophysical environment.

[TABLE]

with

[TABLE]

At small scales, Eq. (4) describes a BH of mass $M_{\text{BH}}$, while at large scales the Newtonian potential corresponds to that of the Hernquist profile (2), dominated by the halo’s mass $M_{\rm halo}$. The spacetime consists of a BH event horizon at $r=r_{h}=2M_{\text{BH}}$, a curvature singularity at $r=0$, while the configuration has an Arnowitt-Deser-Misner (ADM) mass equal to $M_{\textrm{ADM}}=M_{\rm halo}+M_{\text{BH}}$. For galactic dark matter halos, the inequality $M_{\text{BH}}\ll M_{\rm halo}\ll a_{0}$ must hold while the compactness should satisfy the inequality $M_{\rm halo}/a_{0}\lesssim 10^{-4}$ [179], in order to describe the dark matter halo of a galaxy. Even so, in the context of BH environments, the compactness, in principle, can be treated as a parameter of spacetime as long as $M_{\rm halo}<2(a_{0}+2M_{\text{BH}})$, in order to avoid spacetime ambiguities and pathologies.

II.2 Including spin: general-relativistic, rotating BHs in matter halos

Astrophysical BHs are generally expected to rotate and to reside in galaxies that themselves possess intrinsic angular momentum. Incorporating rotation into the models discussed in the previous section is therefore crucial for a complete understanding of how environments influence astrophysical BHs. However, the absence of spherical symmetry appears to rule out the possibility of obtaining closed-form analytic solutions.

Using numerical methods, Ref. [131] has recently obtained stationary and axisymmetric solutions to the Einstein equations, describing rotating BHs immersed in an astrophysical environment. In the following, we describe briefly the formalism and the numerical method used in Ref. [131], to which the reader is referred to for a more detailed discussion.

We consider solutions to the Einstein field equations, $G_{\mu\nu}=8\pi T_{\mu\nu}$, sourced by an anisotropic fluid, with stress-energy tensor

[TABLE]

where $\varepsilon$, $p_{r}$, $p_{1}$ and $p_{2}$ are the energy density, radial pressure and transverse pressures in the co-moving frame of the fluid, respectively. The vectors $u^{\mu}$, $k^{\mu}$ and $s^{\mu}$ are, respectively, the four-velocity of the fluid and spacelike vectors that define the directions of anisotropy such that $u^{\mu}u_{\mu}=-1$, $k^{\mu}k_{\mu}=1$, $s^{\mu}s_{\mu}=1$, $u^{\mu}k_{\mu}=0$, $u^{\mu}s_{\mu}=0$, and $k^{\mu}s_{\mu}=0$.

We use a metric that is stationary, axially-symmetric and circular, in quasi-isotropic coordinates, described by four functions, $f$, $g$, $h$, $\omega$, of $r$ and $\theta$, namely

[TABLE]

where $N_{\pm}=(1\pm r_{h}/r)$. In these coordinates, the four-velocity of the fluid is

[TABLE]

where $\Omega\equiv\Omega(r,\theta)$ is the angular velocity of the fluid. Following Ref. [131] we focus on the simplest case, where $\Omega=\omega$, such that the angular velocity of the dark matter halo is simply the frame-dragging angular velocity.

The total mass of the spacetime obeys $M_{\rm ADM}=M_{\rm BH}+M_{\rm halo}$, where $M_{\rm BH}=\frac{\kappa}{4\pi}A_{h}+2\Omega_{h}J_{h}$, is the mass of the BH444The quantities $\kappa$, $A_{h}$, $\Omega_{h}$ and $J_{h}$ are, respectively, the surface gravity, horizon area of the BH, angular velocity of the horizon, and the angular momentum of the BH., and $M_{\rm halo}$ is the mass of the halo. The total mass $M_{\rm ADM}$ and angular momentum $J_{h}$555Since $\Omega=\omega$, the total angular momentum is equal to the angular momentum of the BH. can also be computed from the asymptotic decay of the metric components $g_{tt}=-1+2M_{\rm ADM}/r+\mathcal{O}(r^{-2})$, and $g_{t\varphi}=-2J_{h}\sin^{2}\theta/r+\mathcal{O}(r^{-2})$, while the surface gravity and the area of the event horizon for the metric (9) are computed as

[TABLE]

at $r=r_{h}$.

To solve the field equations and obtain stationary, axially-symmetric BH solutions, we use the six independent combinations of equations that follow from the Einstein equations (see Ref. [131]). The system, however, contains eight undetermined functions of $r$ and $\theta$, the four metric functions and the four eigenvalues of the stress-energy tensor. The field equations, can only determine six of these functions, that we choose as the four metric functions and the two tangential pressures $p_{1}$ and $p_{2}$. To close the system we choose, from physical considerations, a profile for the energy density $\varepsilon$ of the dark matter halo, and provide an equation of state for $p_{r}$, namely, $p_{r}=0$, in accordance with the Einstein cluster. Motivated by the Hernquist profile, and by Ref. [121], we use the following profile for the energy density

[TABLE]

where $\mathcal{M}$ and $a_{0}$ are parameters with length dimensionality, and $b\equiv b(r)$ is a complicated function of $r$ presented in Ref. [131]. The parameter $a_{0}$ represents a typical length scale associated with the halo, i.e. the limiting radius that includes $95\%$ of the dark matter content of Eq. (12) beyond which its mass falls of as $r^{-4}$. Finally, the parameter $\mathcal{M}$ represents, to leading order in an expansion of powers of $\mathcal{M}/a_{0}$, the mass of the halo $M_{\rm halo}$, see Ref. [131]. The function $b(r)$ guarantees that the parameters $\mathcal{M}$ and $a_{0}$ have suitable physical meaning. Importantly, $b>0$, and therefore $\varepsilon\geq 0$. In the large $r$ limit, this profile for $\varepsilon$ agrees with the Hernquist profile, up to terms proportional to $r_{h}$.

The functions obey the following boundary conditions. At the horizon ($r=r_{h}$) we have $\partial_{r}f=\partial_{r}g=\partial_{r}h=\partial_{r}p_{1}=\partial_{r}p_{2}=0,\quad\omega=\Omega_{h}$, while asymptotically they obey $f=g=h=1,\quad\omega=p_{1}=p_{2}=0$. Regularity, axial symmetry and parity considerations imply $\partial_{\theta}f=\partial_{\theta}g=\partial_{\theta}h=\partial_{\theta}\omega=\partial_{\theta}p_{1}=\partial_{\theta}p_{2}=0$, at $\theta=0$ and $\theta=\pi/2$.

The input parameters are those describing the halo, i.e., $a_{0}$ and $\mathcal{M}$, and the ones associated with the spinning BH, namely $r=r_{h}$, $\omega=\Omega_{h}$, that altogether fully describe the spacetime. Once a solution is obtained with the numerical method, where we keep 20 digits of precision in every calculation, the ADM mass, BH mass, total angular momentum and mass of the halo can be computed with the expressions provided in the previous discussion. We find that there is a one-to-one correspondence between these quantities and the input parameters; solutions can be fully specified either by $(a_{0},\mathcal{M},r_{h},\Omega_{h})$ or by the quantities $(a_{0},M_{\rm halo},M_{\rm BH},J_{h})$. In the following, we choose to present our results in terms of dimensionless values constructed from the second set of quantities, namely the normalized halo compactness $a_{0}/M_{\rm halo}$, the ratio between the halo and BH mass $M_{\rm halo}/M_{\rm BH}$, and the dimensionless spin parameter $J_{h}/M_{\rm BH}^{2}$. These values are presented with an approximate sign ($\simeq$), in what will follow, due to the numerical nature of the solutions. Nonetheless, we present in Table 1 the precise values of $(a_{0},\mathcal{M},r_{h},\Omega_{h})$ used as initial conditions for each solution used throughout this work, as to ensure reproducibility of our results.

To solve the system of partial differential equations, we use the same code as in Ref. [131], which is publicly available at [181]. This code is a high-precision version of the code presented in Ref. [182], originally developed by one of the authors. This code combines a pseudospectral method with the Newton-Raphson root-finding algorithm. Details on the code, its validation and accuracy are provided in the appendix of Ref. [131]. We observe exponential convergence of the code as we increased its resolution, and only accepted a solution when the absolute estimated error was below $\mathcal{O}\left(10^{-8}\right)$.

III Geodesic motion

III.1 Equations of motion and conserved quantities

We operate on a generic metric, i.e. a stationary and axisymmetric spacetime written as

[TABLE]

where the metric tensor components in (13) are functions of $r$ and $\theta$, and the coordinate system $\left(t,r,\theta,\varphi\right)$ can be chosen to be of Boyer-Lindquist [183] or quasi-isotropic type [184]. A zero-order assumption for the evolution of an EMRI is to think of the motion of the secondary as a test-particle following the geodesics of the spacetime geometry of the primary. Under this simplistic assumption, which serves as a good proxy of a short-timescale EMRI, that corresponds to some dozens of revolutions around the primary, can be described by the geodesic equations

[TABLE]

where $\Gamma^{\kappa}_{\lambda\nu}$ are the Christoffel symbols of the primary, $x^{\kappa}$ is the four-position of the secondary, and an overdot denotes differentiation with respect to proper time $\tau$. The system of second-order, coupled differential equations resulting from Eq. (14) can be considerably simplified using constants of motion. The spacetime metric (13) bears two Killing vector fields resulting from stationarity and axisymmetry, giving rise to two constants of motion throughout geodesic motion, i.e. the specific energy $E$ and azimuthal angular momentum $L_{z}$ of the test particle of mass $\mu$

[TABLE]

Equations (15) can be written as two, first-order, decoupled, ordinary differential equations for the $t$ and $\varphi$ momenta, so that

[TABLE]

There are two remaining geodesics for $r$ and $\theta$ which, in general, form a second-order dynamical system of coupled differential equations. Nevertheless, test particles in geodesic motion provide a third constant of motion, i.e., the conservation of their rest mass, or equivalently their four-velocity $g_{\lambda\nu}\dot{x}^{\lambda}\dot{x}^{\nu}=-1$. This leads to a constraint equation

[TABLE]

with $V_{\text{eff}}$ a potential of the form

[TABLE]

Equation (17) characterizes bound geodesic motion, through the potential $V_{\text{eff}}$. In turn, when $V_{\text{eff}}=0$ a curve appears (under a choice of $\theta\in[0,\pi]$); the curve of zero velocity (CZV), where the $r$ and $\theta$ momenta are bound to be null there, i.e., $\dot{r}=\dot{\theta}=0$.

Under the hypothesis of another constant of motion, the geodesics for all degrees of freedom can be separated into decoupled, first-order, differential equations. If though the primary spacetime does not allow for further spacetime symmetries then the $r$ and $\theta$ motion remains coupled, and of second differential order. A Kerr BH, though, possesses a rank-two Killing tensor field, that associates the remaining components of the angular momentum with the Carter constant [132]. This successfully decouples the motion of $r(\tau)$ and $\theta(\tau)$. Since the primary configurations we will use here do not necessarily have a Carter constant (or any constant related to a higher-rank Killing tensor field), we will work under the assumption of the absence of a fourth constant of motion, and the phase space of orbits will inform us if a separation constant exists in our numerical spacetime. To evolve orbits we will utilize the coupled second-order differential equation system for $r$ and $\theta$, together with Eqs. (16) and (17), without any further symmetry assumptions.

III.2 Evolving geodesics on numerical backgrounds

Adopting quasi-isotropic coordinates, the line element (13) can be described through four independent functions $(f,g,h,\omega)$ of $(r,\theta)$, where the metric tensor components result directly from Eq. (9), i.e.,

[TABLE]

The solutions have been found numerically following the procedure described in the previous section. The details of the numerical implementation are provided also in [182, 131]. The geodesic evolution of massive particles around rotating BHs in a halo of Hernquist-type matter, uses the numerically-constructed metric tensor components (19), as in [149], and geodesics are evolved with the use of high-order interpolation functions that are constructed from the data of the metric tensor. For all simulations performed herein, we find that the constraint equation (17) is satisfied to within one part in $10^{8}$ for the first $\sim 20$ thousand revolutions (and $\sim 20$ thousand intersections through the equatorial plane). In all cases we fix the mass ratio of the EMRI to $\mu/M_{\textrm{BH}}=10^{-6}$ (even though at the geodesic level it just leads to a change of scales) and choose the constant energy and azimuthal angular momentum as $E/\mu=0.9$ and $L_{z}/\mu=3M_{\textrm{BH}}$, respectively. This choice stems from the fact that the most dominant resonances are located in the strong-field regime, and the aforementioned values for $E$ and $L_{z}$ gives us access to these resonances. Furthermore, due to the effect of circularization of orbits that takes place during GW emission, the eccentricity should be significantly decreased, with respect to the initial eccentricity at the astrophysical capture of the secondary, that can be very high. The particular choice of $E$ and $L_{z}$ facilitate eccentricities approximately smaller than $0.4$, thus are good quantities to initialize a somewhat more realistic inspiral. Nevertheless, there is a large choice of parameter space values where small eccentricities can be achieved, thus our choice is neither fine-tuned, nor random.

III.3 Non-integrability and signatures of chaos

The geodesics of massive test-particles around Kerr BHs possess four constants of motion, namely the energy, the azimuthal angular momentum, the Carter constant and the four-velocity of the particle. Thus, the system of equations is (Liouville-)integrable and does not present chaotic features [133]. Since the Carter symmetry is very fragile it may not exist when the multipolar structure of spacetime has been significantly deformed (see however Ref. [185]). The absence of a Carter symmetry, or in general any other (higher-rank) Killing tensor field, leads to the non-integrability of geodesics; the radial and polar sectors of the geodesics are not separable, but rather form a second-order, coupled dynamical system. Some cases, where the Carter constant ceases to exist, have been discussed in the Introduction. The absence of a fourth constant of motion inevitably leads to non-integrability, which allows for chaotic effects [133]. Full-blown chaos is extremely rare, if not completely absent, at least in astrophysical EMRIs. The underlying effects of non-integrability, though, can manifest themselves indirectly in the phase space of orbits and in particular close to orbital resonances. Transient orbital resonances already affect EMRI evolution, even in Kerr, where geodesics are integrable [136, 137, 186, 187], and introduce significant dephasing effects. Non-integrable EMRIs are very likely to amplify these dephasing effects and can introduce even stronger, clear-cut phenomenology of indirect chaos [139, 140, 141].

In integrable systems, geodesics are orbits that occupy the surface of a two-dimensional torus in the available four-dimensional phase space. The torus of each geodesic is called an invariant torus and is described by the minimum and maximum radii of the orbit, associated with the eccentricity, and by the minimum and maximum polar angles, associated with the inclination. Together with the azimuthal orbital period, any trajectory can be described by three fundamental frequencies; the frequency of oscillation from the periapsis (minimum radius) to the apoapsis (maximum radius) and back, $\omega_{r}$, the oscillation frequency through the equatorial plane, between the maximum and minimum polar angle, $\omega_{\theta}$, and the revolution frequency $\omega_{\varphi}$.

If one considers a two-dimensional surface that cuts through a foliage of invariant tori, named Poincaré surface of section, then each torus defines a closed curve on this surface [188]. This curve is called an invariant curve. Each torus corresponds to a characteristic pair of frequencies, as discussed above. Not only the two frequencies, but also the ratio between them varies continuously over successive torii. If the ratio of the frequencies is an irrational number, an orbit continuously winds around its corresponding torus covering its surface densely. This kind of orbit is called quasi-periodic, or generic. A generic orbit goes repeatedly through a surface of section defining a succession of points which eventually cover densely the corresponding invariant curve on the surface of section. In the case where the ratio of frequencies is a rational number $n/m$, where $n,m\in\mathbb{N}$, the orbit repeats itself after $m$ windings. The orbit is periodic and the corresponding torus is called resonant. In this case, the resonant invariant curve consists of an infinite number of $m$-multiplets of periodic points. Each $m$-multiplet represents an $m$-multiple periodic orbit.

When non-integrable perturbations are presented in an integrable system, two theorems, namely the Kolmogorov-Arnold-Moser (KAM) [189, 190] and Poincaré-Birkhoff [191] theorems, dictate the modification of orbital dynamics. The KAM theorem ensures that orbits sufficiently away from resonances, are smoothly shifted but, in general, most of the tori are deformed but not destroyed. These tori are called KAM tori. Thus, the corresponding surface of section looks very much like the surface of section of the corresponding integrable system. The generic invariant curves in non-integrable systems are called KAM curves. Successive KAM curves organize around a common central point, which corresponds to a spherical orbit (zero eccentricity).

Close to resonances, the KAM curves disintegrate into two sets of periodic points, in accord with the Poincaré-Birkhoff theorem. The stable periodic points, from which stable resonant orbits emanate, are encapsulated by KAM curves that form a nested island, i.e., a resonant island. A phase-orbit of such a resonant case visits all the $m$ islands of the $n/m$-resonance, moving successively to the next $n$-th island along the aforementioned closed curve at every winding, forming eventually the KAM curves inside every island. On the other hand, the unstable periodic points generate chaotic orbits that surround the resonant islands with thin layers. The structure around resonances of non-integrable systems is called a Birkhoff chain. The crucial aspect of resonant islands, and their ultimate significance in EMRI dynamics, is the fact that the rational ratio of the central, stable periodic point $\omega_{r}/\omega_{\theta}=n/m$ is shared throughout all KAM curves residing inside the island, no matter on which KAM curve inside the chain of islands it belongs. This property is not shared by the non-resonant KAM curves, since the ratio of frequencies on them is irrational and it varies smoothly from one KAM curve to another.

In a sense, integrable EMRIs experience resonances that occupy a “zero-volume point” in phase space, while non-integrable EMRIs exhibit prolonged “non-zero-volume” resonances where the secondary is locked in perfect resonance for a significant amount of revolutions around the primary. For example, in non-Kerr EMRIs, these prolonged resonances may last up to $\sim 200-300$ cycles depending on the primary [141, 142], without taking into account pre- and post-resonant effects [137, 136, 135] or the conservative part of gravitational self-force [192].

The existence of resonant islands serves as an indirect proof of non-integrability in EMRIs, therefore sketching a Poincaré map can inform us regarding the system’s symmetries. Furthermore, through the Poincaré map, the rotation number can be calculated in order to scan for resonant islands. It practically tracks the angle $\vartheta$ between successive intersections on KAM curves, relative to the fixed central point of the Poincaré map. The rotation number is defined as the summation of $\mathcal{N}$ angles $\vartheta$ measured between successive intersections of each geodesic in the surface of section, i.e.,

[TABLE]

When $\mathcal{N}\rightarrow\infty$, Eq. (20) asymptotes to the ratio $\nu_{\vartheta}=\omega_{r}/\omega_{\theta}$, thus the rotation number can serve as a resonance browser. Consecutive rotation numbers found from different initial conditions of geodesics, by smoothly varying one of the system’s parameters while keeping the rest fixed, forms a rotation curve.

Non-integrable systems display discontinuities in the monotonicity of a rotation curve through the formation of plateaus with non-zero widths, when geodesics transverse resonant islands. Therefore, the width of the plateau corresponds to the width of the resonant island and is a general measure of the timescales involved in resonant islands of non-integrable EMRIs. Inflection points can also appear in the rotation curve when trajectories pass through unstable periodic points. In what follows, we will examine the characteristic features of rotating BHs in halos of matter with geodesic evolutions, in order to search for signatures of non-integrability driven by astrophysical environments. We note that Poincaré maps and rotation curves are not the only diagnostics of chaotic phenomena. For more diagnostics that can be applied in curved spacetimes see [193].

IV Chaos in relativistic, rotating, non-vacuum black holes

In what follows, we will use relativistic, rotating BHs in a Hernquist-type halo, with a variety of spin parameters and halo compactness, while fixing the central BH mass to unity, since all our numerical evolutions lie at the geodesic level, where all observables can be appropriately rescaled for any BH mass to simulate an EMRI, thus $M_{\textrm{BH}}=1$ is merely chosen for simplicity. In turn, the test particle that will be used has the fixed constants of motion $E/\mu=0.9$, $L_{z}/\mu=3M_{\textrm{BH}}$. We assume a test-particle mass so that the mass-ratio is kept to $\mu/M_{\textrm{BH}}=10^{-6}$, even though this is merely a scale factor in the geodesic evolution scheme. We utilize the numerical configurations in quasi-isotropic coordinates as described in Eq. (13) and normalize the radial quasi-isotropic coordinate $r$ with respect to the event horizon radius $r=r_{h}$ of each configuration.

As a proof-of-principle, we numerically construct a rapidly-rotating BH with spin $J/M^{2}_{\rm BH}\simeq 0.996$ and a Hernquist-type halo with compactness $M_{\textrm{halo}}/a_{0}\simeq 10^{-1}$, $M_{\textrm{halo}}\simeq 10M_{\rm BH}$. The choice of parameters is exaggerated in order to sketch a clear-cut phase-space structure of various bound geodesics that we evolved, which might show chaotic phenomenology, as discussed in Subsection III.3. The initial position of geodesics is chosen such that $r(0)$ is varying inside the CZV, while the rest of the position vector components are fixed to $\dot{r}(0)=0,\,\theta(0)=\pi/2$ and $\dot{\theta}(0)$ is chosen in accord with the constraint in Eq. (17) so that the resulting orbits remain bounded inside the CZV.

As each geodesic evolves, we dynamically capture the position and velocity $(r(t),\,\dot{r}(t))$ of the orbit, each time it crossed the equatorial plane, with $\dot{r}(t)>0$. This process defines a curve of successive intersections for each geodesic. Putting together the various curves that correspond to different geodesics we sketch the Poincaré map in Fig. 1. It is evident that we have managed to capture, at least, three resonances, i.e., two resonant islands namely the $\omega_{r}/\omega_{\theta}=2/3$, $2/5$, shown with red and blue, respectively, and a thin chaotic layer, that surrounds the $1/2$-resonant island, shown in green.

Figure 1, left panel, shows the structure of the map. The black curves encircle the central point of the map. These correspond to the KAM curves, away from resonances, that smoothly depart from the invariant curves of integrable systems. In fact, most non-integrable systems associated with EMRIs possess KAM curves that are similar to invariant ones, due to the fact that the perturbation inserted in the Hamiltonian is sufficiently small, in order to not completely disintegrate the whole phase space of orbits. Therefore, the structure of phase space does not entirely break down. Into the strong-field regime, though, numerous islands appear that do not encircle the central point of the Poincaré map, but rather encapsulate the stable periodic points of resonances. The island structure around the resonance $2/5$ and $2/3$, shown with blue and red colors, respectively, serve as imprints of chaos and non-integrability of the system under consideration. Figure 1, right panel, shows a clear, zoomed-in, depiction of the islands’ nested structure. Interestingly, we do not only observe resonant islands, but we also show a chaotic layer, that surrounds densely the $1/2$-resonant islands. This occurs due to the fact that we have not crossed through the resonant island, but rather traversed within the vicinity of the chaotic zone forming around them. This discussion is very relevant when one sketches the corresponding rotation curve that results from the structure of successive intersections in the map.

In Fig. 2, left panel, we demonstrate the rotation curve resulting from successive bound geodesics in the BH configuration of our study. Here, we used a rather dense arrangement of successive orbits to achieve a smooth curve. We observe that in most of its parts, the curve increases monotonously as $r(0)$ is increased. A more thorough look in the curve reveals the outcome of the emergence of resonant islands in the Poincaré surface. When $\nu_{\vartheta}$ encounters the resonance $2/3$, we can clearly identify a plateau formation (see red box in Fig. 2, that corresponds to the inset in the left panel) that breaks the curve’s monotonicity and effectively renders the system in study non-integrable. The inset in the left panel of Fig. 2 shows the corresponding plateau encountered when $\nu_{\vartheta}=2/3=\omega_{r}/\omega_{\theta}$, while the right panel shows the plateau formed due to the $\nu_{\vartheta}=2/5$ resonant island. We note that we do not show how the crossing through the unstable $1/2$-periodic point is imprinted in the rotation curve. This is due the chaoticity of successive intersections in the Poincaré map, that leads to rotation numbers that are not well-defined. The result is an erratic oscillation (or overall an inflection point) in the rotation curve (see e.g. Fig. 9, bottom left panel, in [149]). The appearance of plateaus, and inflection points, instantly corroborates the Poincaré-Birkhoff theorem, which states that geodesics residing inside resonant islands share the same rational ratio $\omega_{r}/\omega_{\theta}$, or $\nu_{\vartheta}$ when the number of intersections measured ($\mathcal{N}$) is sufficiently large. In this analysis, we have recorded approximately 5 to 10 thousand intersections per geodesic, that corresponds to an accuracy of 3 to 4 digits in the resulting rotation numbers. A thorough check of geodesics in further non-vacuum, rotating BH configurations, with smaller BH spin, or smaller compactness, reveals a qualitatively similar structure in the Poincaré map and the rotation curve. Hence, adding rotation to the static solution of Ref. [121], breaks integrability and leads to chaotic phenomenology, that extends the lifetime of an non-integrable EMRI in resonance [139, 140, 141].

From Fig. 1 and 2, we can conclude that the resonance $2/3$ is the strongest, i.e., having the largest island width (and corresponding plateau in the rotation curve) in non-Kerr EMRIs. This is due to the multiplicity and location of the resonance $2/3$ in the strong-field regime [133]. This is the underlying reason why the literature is usually studying the effects of the $2/3$ resonance in Kerr [134, 136, 135, 137] and non-Kerr EMRIs [139, 140, 141, 142].

In what follows, we present results regarding the behavior of the $2/3$-resonant island, and in particular its width, when the spacetime’s BH spin and the halo compactness are varied. The following is an order-of-magnitude analysis and will help us qualitatively understand the role of the BH spin and the halo in EMRIs, which combined break the Carter constant and the eventual integrability of geodesics around these objects.

IV.1 $2/3$ island-width dependence on the BH spin and halo parameters

Besides providing a proof-of-principle for the non-integrability of general-relativistic BHs surrounded by a matter halo, it is interesting to examine how the different parameters of spacetime affect the chaotic signatures, i.e., the width of the resonant islands. For simplicity, we focus on the $2/3$-resonant island since it provides the widest islands in phase-space, and plateaus in the corresponding rotation curves. Nevertheless, we have checked that the effect observed for the $2/3$-resonance holds for all other resonant islands that we could probe.

Figure 3 shows how the $2/3$-island width changes when we vary the BH spin and the halo compactness. In both cases, increasing the respective parameter leads to an enhancement of the island’s width. Therefore, rapidly rotating BHs surrounded by compact matter halos are expected to provide the most pronounced imprints of non-integrability and chaos666See also Ref. [194].. Thus, although EMRIs of supermassive BHs in galaxies, surrounded by a dark matter halo may not provide detectable chaotic signatures, stellar-mass BHs surrounded by a compact environment, such as BH hair or boson clouds, can produce much more prominent imprints of chaotic phenomenology. Interestingly, the latter configuration can be supported in lower mass-gap EMRIs [195].

Other studies regarding geodesics and EMRIs around non-Kerr objects [140, 196] have also demonstrated that the increment of the BH’s spin leads to an increase in the resonant-island widths, therefore this result was expected and simultaneously corroborates the validity of our analysis. On the other hand, the island widths are also increasing with the halo compactness. Hence, we may extrapolate that the mass quadrupole moment, and higher multipoles are deformed in a prolate manner due to the increment of the halo’s compactness (see e.g., the similarities in the multipole moments of Hartle-Thorne [196] and Manko-Novikov objects [140] that produce broader islands for larger prolate deformation parameters; exactly as we observe in Fig. 3, right panel).

V Conclusions

Resonances are choreographically well-defined periodic orbits, that return to their initial position after a number of oscillations. When a dynamical system is non-integrable, resonances are extended to resonant islands, which occupy an atypical, non-zero volume in phase space, where the periodicity of the central resonance is shared throughout the geodesics in the island. This leads to orbital and, most importantly, GW observables [141] that may be detected with LISA through abnormally-large cumulative dephasing events, resulting in imprecise parameter inference [134].

In order to study astrophysical EMRIs, a vacuum Kerr background is typically utilized with manual additions that express the environment’s contribution to the dynamics of orbits. These approximates are not solutions to GR, but rather include the environmental effects through, either Newtonian/post-Newtonian schemes, or with fully-relativistic treatments that mimic the dynamical interaction of the secondary with the environment, known as dynamical friction [76]. A reasonable understanding of the main astrophysical features in GW generation and propagation has, thus, been achieved [57, 62, 61, 80]. The results so far are quite straightforward; the dynamics of EMRIs in astrophysical environments need to be treated in a fully-relativistic manner, rather than with post-Newtonian methods, in order to capture the correct environmental dephasing of GWs from EMRIs [77, 76, 73]. Nevertheless, at the level of the metric, the spacetime used in contemporary studies still assumes a vacuum (integrable) Kerr BH perturbed by a non-backreacting matter environment.

Here, we surpass crude estimates, at the level of the metric, by constructing and employing general-relativistic solutions of non-vacuum GR, that describe rotating BHs at the center of matter halo environments [131]. We have evolved geodesics in order to capture the leading-order effects in EMRI dynamics. Our analysis demonstrates that when the environment backreacts on the rotating Kerr BH, it produces a non-vacuum solution of GR that significantly deforms the multipolar spacetime structure to the point that the Carter symmetry is broken. The absence of a Carter constant alters the dynamics of geodesics around transient resonances, in a manner that designates the emergence of chaotic phenomena and non-integrability, dubbed environmental chaos. We have explicitly found resonant islands in the phase space of orbits, and thin chaotic layers that form around them, thus astrophysical environments can indeed lead to non-integrability when considered as a general-relativistic part of spacetime, and not just as a perturbation around a vacuum Kerr system. The rotation curves corroborate the aforementioned results, due to the formation of plateaus when the geodesics are inside the respective resonant island. We, finally, performed a tentative analysis of the widths of various $2/3$-resonant islands resulting by varying the BH spin and the compactness of the halo. We observe that the increment of both spin and compactness, typically, leads to larger, in span, resonant islands.

Even though our work lies at the geodesic level, it has been shown that when radiation reaction is taken into account, the resulting non-Kerr EMRIs can potentially spend hundreds of revolutions inside a resonant island [139, 140], where the periodicity of the central resonance is shared throughout it. The crossing of the secondary through successively-damped, “resonantly-locked”, geodesics leads to a GW observable, namely a GW glitch in the frequency evolution of the non-Kerr EMRI [141]. These glitches have a clear astrophysical interpretation [142, 148, 149], can last for weeks [142], instead of minutes or hours as in Kerr EMRIs [135], and are not associated with the typical instrumental glitches that occur due to detector noise [197, 198, 199, 200, 201, 202, 203]. Thus, we can extrapolate that if our analysis is extended at the inspiral level, then similar glitches will appear in the GW frequency evolution of general-relativistic, matter-enriched rotating EMRIs.

In the dawn of GW astronomy, the astrophysical environments around realistic, spinning BHs, such as those residing in the core of galaxies, have encountered a significant interest from the GW community, especially for sources like EMRIs that are prime targets of the LISA mission. Our geodesic analysis demonstrates that the assumption of vacuum Kerr BHs with a perturbative, but not backreacting, addition of the astrophysical environment is not enough to describe the general-relativistic nature of the full system. By constructing a general-relativistic BH within a matter halo from first principles, using the Einstein cluster formalism [121], a complete understanding of environmental effects around BHs can be achieved at the level of the metric. Nevertheless, EMRIs evolve under radiation reaction, hence they need to be treated with fully-relativistic schemes at the dynamical level as well [76, 77, 73] (instead of using Newtonian or post-Newtonian treatments).

The relativistic treatment of inspiral dynamics has already been achieved for exact GR solutions that describe matter-enriched, static and spherically-symmetric BHs [127, 204, 100]. The final missing piece to the puzzle is to generalize the fully-relativistic treatment of EMRIs of Refs. [127, 100] to include the rotation of the primary. Nevertheless, it is easier said than done because the rotating astrophysical BHs used in this work do not have a closed analytic form, but rather are constructed numerically. On top of that, the generalization of the relativistic techniques prescribed in [127, 204, 76, 100] needs to be accomplished, i.e., to transform from the Regge-Wheeler-Zerilli-matter equations (zero BH spin) to Teukolsky-matter relativistic equations (non-zero BH spin).

Even if the above generalization proves to be a thorny and time-consuming process, we can still infer some rough detectability aspects, that are an outcome of the current geodesic work. A typical resonance in a vacuum Kerr EMRI produces dephasing of dozens of radians [134, 205, 135], without taking into account the dephasing that an astrophysical environments can further introduce. Imagining that a zero-volume single point of geodesic periodicity in phase space creates such an effect in Kerr EMRIs [136], we could cautiously extrapolate that a finite volume of geodesic periodicity in phase-space, which is introduced by the formation of resonant islands, will produce even bigger resonant enhancement through unexpectedly extended resonance timescales. This could be in the form of beyond-Kerr-resonance dephasing of hundreds of radians, or even abrupt dephasing kicks when the island is crossed (see Refs. [141, 142] for the resulting GW imprint of a resonant-island crossing). These effects alone might serve as conclusive evidence of the existence of astrophysical environments, and in general non-integrability.

Thus, it is important that resonances, and in particular possible crossings through resonant islands, to be taken into consideration with novel state-of-the-art methods (see e.g. Ref. [206, 73]). Since GW modeling of generic Kerr EMRIs is a highly active subject of research, so does the GW modeling of resonant-island crossings should be explored for non-vacuum, general-relativistic spinning EMRIs. The combination of environmental effects, resonances, and astrophysical chaos, could form, in the future, the amalgamation of a well-explored branch of precise GW modeling for realistic EMRIs within astrophysical environments. One, though, needs to stay cautious for potential modifications of gravity, supermassive exotic compact object primaries, or triple systems [207, 208, 209, 210, 211, 212, 213, 214], in order to avoid numerous degeneracies.

Acknowledgements.

The authors would like to warmly thank Vitor Cardoso for fruitful discussions. K.D. acknowledges financial support provided by FCT–Fundação para a Ciência e a Tecnologia, I.P., under the Scientific Employment Stimulus – Individual Call – Grant No. 2023.07417.CEECIND/CP2830/CT0008. K.D. would also like to 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. P.F. is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1 - 390900948 (the Heidelberg STRUCTURES Excellence Cluster). This project has received funding from the European Union’s Horizon MSCA-2022 research and innovation programme “Einstein Waves” under grant agreement No. 101131233.