Geodesics and gravitational waves in chaotic extreme-mass-ratio inspirals: the curious case of Zipoy-Voorhees black-hole mimickers

Due to the growing capacity of gravitational-wave astronomy and black-hole imaging, we will soon be able to emphatically decide if astrophysical dark objects lurking in galactic centers are black holes. Sgr A*, one of the most prolific astronomical radio sources in our galaxy, is the focal point for tests of general relativity. Current mass and spin constraints predict that the central object of the Milky Way is supermassive and slowly rotating, thus can be conservatively modeled as a Schwarzschild black hole. Nevertheless, the well-established presence of accretion disks and astrophysical environments around supermassive compact objects can significantly deform their geometry and complicate their observational scientific yield. Here, we study extreme-mass-ratio binaries comprised of a minuscule secondary object inspiraling onto a supermassive Zipoy-Voorhees compact object; the simplest exact solution of general relativity that describes a static, spheroidal deformation of Schwarzschild spacetime. We examine geodesics of prolate and oblate deformations for generic orbits and reevaluate the non-integrability of Zipoy-Voorhees spacetime through the existence of resonant islands in the orbital phase space. By including radiation loss with post-Newtonian techniques, we evolve stellar-mass secondary objects around a supermassive Zipoy-Voorhees primary and find clear imprints of non-integrability in these systems. The peculiar structure of the primary, allows for, not only typical single crossings of transient resonant islands, that are well-known for non-Kerr objects, but also inspirals that transverse through several islands, in a brief period of time, that lead to multiple glitches in the gravitational-wave frequency evolution of the binary. The detectability of glitches with future spaceborne detectors can, therefore, narrow down the parameter space of exotic solutions that, otherwise, can cast identical shadows with black holes.


Introduction
The Schwarzschild spacetime [1] is unequivocally the simplest and most remarkable black hole (BH) solution of general relativity. It describes a vacuum compact object with an event horizon and a static, spherically-symmetric exterior. The remarkable symmetry properties of Schwarzschild geometry places it at the top of Einstein field equation solutions in what regards its simplicity and singular externally-observable property; the gravitational mass. Its ultimate successor, the Kerr metric [2], is undoubtedly the most successful and astrophysically-relevant solution of the vacuum field equations that describes a spinning, stationary and axisymmetric BH with an oblate, spheroidal external geometry. Despite the fact that Kerr BHs possess less symmetries than Schwarzschild spacetimes, they compensate by including spin; a very crucial aspect of most astrophysical compact objects, and further form an integrable (separable) system of equations of motion for massless and massive particles due to the existence of the Carter constant [3]. The above statement of integrability translates to the absence of chaos in geodesics around Kerr BHs, and is traced trivially to Schwarzschild spacetime [4].
From an astrophysical perspective, a significant volume that surrounds BHs in the Universe consists of plasma, accretion disks, matter configurations and halos that extend significantly far away. Thus, it is hard for one to realize BHs in pure vacuum, especially those occupying the centers of galaxies, where a perplex and highly dynamical environment is present. The historic observations of shadows from supermassive compact objects in the center of M87* galaxy [5] and Sgr A* in our galaxy [6] has opened a new realm of observational yield with BH imaging. An accretion disk can, in principle, deform the surrounding geometry of a BH so that it continuously deviates from a Schwarzschild or Kerr description, causing degeneracies between a multitude of other exact, though more exotic, solutions that may mimic the observed shadow of supermassive compact objects.
To study potential degeneracies present in current electromagnetic observations, the literature usually operates in spacetime deviations from the Kerr description, known as bumpy/parameterized [7][8][9][10][11][12][13][14][15] or non-Kerr compact objects [16][17][18][19][20][21][22]. A distinctive feature of some of these objects is the absence of Carter symmetry, which leads to chaotic phenomena in particle-dynamics. In what regards Sgr A*, though, spin is not as crucial as the deviation from spherical symmetry itself, since current constraints on its spin conclude that it is less than 10% of its maximal allowed [23,24]. Therefore, a plethora of Sgr A*-like objects in the Universe can be sufficiently modeled as Schwarzschild (or slowly-rotating Kerr) BHs or exotic BH mimickers.
Gravitational-wave (GW) astrophysics has been proven to be an extraordinary tool to break the aforementioned degeneracies, thus synergies between GW and shadow observations are necessary in order to search for the ultimate spacetime description of known astrophysical compact objects. To that end, the LIGO/Virgo/Kagra collaboration has been flooding our databases with numerous GW detections from coalescing compact objects [25]. Such detectors, although extremely successful so far (at the level of groundbreaking), have the unfortunate attribute of being plagued by planetary noise, since they are placed on Earth, and present unavoidable limitations due to their scale, that play a crucial role in their target span and precision.
The Laser Interferometer Space Antenna (LISA) [26] is a space-borne GW detector that will open new realms in GW astrophysics, due to its unprecedented level of accuracy [27][28][29][30]. It will target, in particular, mHz sources of GWs which are undetectable with current ground-based detectors. One of the prime objectives of LISA (and other space-based programs [31][32][33]) is the detection of GWs from extreme-mass-ratio inspirals (EMRIs) [34,35], which involve a primary supermassive compact object, such as those lurking in galactic cores, and a secondary stellar-mass compact object. We expect that the detection of EMRIs will allow for exquisite measurements of EMRI parameters, the multipolar structure of the primary of the EMRI [36][37][38][39][40][41], and, in general, very accurate tests of GR [42]. Currently, the most significant tool of perturbations theory to study EMRI evolution is the gravitational self-force [43,44], though it has been proven to be a daunting task to undertake and has only been calculated up to second order till date [45,46]. However, other techniques, such as those that will be used in this work, has proven accurate enough to capture the dissipative effects of EMRI evolution, when compared to Teukolsky-based waveforms [47], which utilize weak field post-Newtonian techniques to find the fluxes together with general-relativistic kludge schemes to evolve the EMRI accurately enough [34,48,49] even when the primary spacetime of the EMRI is quasi-Kerr in nature [8,22,[50][51][52][53][54][55]. Environmental effects and spacetime deformations in EMRIs should, therefore, be taken into account in order to maximize the scientific yield from these sources [36][37][38]. Fortunately or not, it is atypical to find integrable spacetimes, especially when complex astrophysical environments are involved or the primary is an exotic compact object [78]. To that end, a significant spacetime degeneracy is present in Sgr A* and M87*. Current investigations have concluded that an abundance of exotic objects can cast shadows that are practically indistinguishable (in specified regions of their parameter space) from those of Schwarzschild (or Kerr) BHs [79][80][81][82][83][84][85][86], therefore assuming that Sgr A* and M87* are BHs, only from their shadow silhouette, can lead to significant misinterpretation. Indeed, we need to invoke geodesics [87], accretion disk analyses [88] and GW ringdown tests [89][90][91][92][93][94][95][96][97][98][99][100][101] in order to understand if supermassive compact objects are typical BHs or exotic in nature.
In this study, we will investigate the simplest deformation of Schwarzschild geometry, the Zipoy-Voorhees (ZV) metric [102,103] (also known as the γ -metric, that has also been generalized in [104,105] to include electric charge), that describes a vacuum, static, and spheroidal solution of Einstein equations which is continuously connected to Schwarzschild by a deformation parameter δ (in our convention). The ZV metric, presenting a spheroidal deformation of an otherwise static geometry can pose as a good model for a static compact object surrounded by a compact environment or an accretion disk, such as those residing in galactic centers. Recent shadow investigations [81] have shown that when the deformation parameter δ > 1 then very precise measurements will be needed in order to rule out an exotic compact object described by this geometry. A peculiarity of this solution is the appearance of a curvature singularity at its surface, thus characterizing it as a naked singularity. Another important feature for our analysis is that the ZV metric has a non-zero mass quadrupole and nonintegrable geodesics [106][107][108], despite earlier claims of integrability [109,110]. Since non-integrable EMRIs present very distinctive characteristics in phase space, such as prolonged resonant islands where geodesics share the same rational ratio of orbital frequencies [22,53,[111][112][113][114] and discontinuities in the GW frequencies during island crossings ('glitches') [54,55], here we examine if similar effects are present in ZV EMRIs, and in particular GW glitches from crossings of subdominant resonant islands.
We confirm that the conclusions of Refs. [106][107][108] are correct, that is the ZV metric is non-integrable due to the existence of chaotic layers of plunging geodesics and a series of resonant islands of stability. We further choose a primary supermassive compact object described by the ZV metric and evolve EMRIs with stellar-mass secondaries (with fixed mass ratio) to assess the effect of non-integrability at the orbital and waveform level. We find plateaus in the dissipative evolution of the ratio of radial and polar frequencies, that designate the crossing of resonant islands, and subsequently observe glitches in the GW frequency evolution of the EMRI. Due to the atypical structure of the ZV primary, a variety of successive resonant islands accumulate close to the edge of bound geodesics. By evolving EMRIs through successive resonances, with generic initial conditions, we find that consecutive glitches can appear in short timescales (of order of several days) in the GW frequency evolution of the inspiral. Since typical glitches experienced by non-Kerr EMRIs, with the same mass ratio as the ZV one, are usually separated by months or even years of dissipative evolution [54,55], the detection of multiple GW glitches in brief periods of time may demonstrate that very slowly-rotating supermassive compact objects are not Schwarzschild (or slowlyrotating Kerr) BHs. These findings will contribute in placing tighter constraints on exotic geometries, such as naked singularities, and narrow down the parameter space of viable BH mimicker primaries that can imitate the shadow of supermassive BHs.
In what follows we use geometrized units so that the gravitational constant and speed of light satisfy G = c = 1.

the Zipoy-Voorhees metric
The ZV spacetime [102,103] describes a two parameter family of exact vacuum solutions to the Einstein equations that are static, axisymmetric and asymptotically flat. The line element in Erez-Rosen coordinates [115,116] reads From Eqs. (2.1), (2.2), it is straightforward to obtain the Schwarzschild limit when δ = 1. Since at the Schwarzschild limit, k = M, δ can be interpreted as a measure for how much more (or less) mass M = kδ the ZV object has when compared to a Schwarzschild BH. Subsequently, the deformation δ captures the oblateness of the compact object. When δ > 1 the ZV geometry describes a spacetime around the central object that is more oblate than a Schwarzschild BH, while when 0 < δ < 1 the central object is more prolate. For δ = 0 (which identically means M = 0) we obtain the Minkowski spacetime. According to the no-hair theorem, when 0 < δ = 1 holds the event horizon is broken; a true curvature singularity appear at r = 2k [117], besides the typical one at r = 0, and the ZV metric describes a naked singularity [118]. Specifically, when writing the solution in Weyl coordinates (ρ, z), which are associated with the prolate spheroidal coordinates (x, y) that many pathological solution of GR are written in (see e.g. the Manko-Novikov solution [18]), then it has been shown that the curvature singularity naturally appears along the line segment ρ = 0, if |z| < k. For more information regarding the nature of the ZV naked singularity see [118][119][120][121]. Interestingly, when spherical symmetry is broken the geometry obtains a non-zero quadrupole moment M 2 = δ(1 − δ 2 )M 3 /3 [122,123] which eventually leads to the Carter constant (or any other higher order Killing tensor) being broken [106][107][108].

Orbital dynamics
Regardless of its peculiar causal structure, the ZV primary we will focus on should possess as many 'good' features as those of astrophysical compact objects. It has been shown that the line element (2.1) has an innermost stable circular orbit (ISCO) when δ > 1/ √ 5 at [88,107] Specifically, the typical ISCO of Schwarzschild BHs that is a special stable, circular and equatorial orbit of the separatrix manifold, that separates between plunge and bound motion, changes structure when δ = 1 and asymptotic manifolds emanate from an unstable orbit, called the Lyapunov orbit [124,125]. Orbits that cross the Lyapunov orbit will plunge unless certain E and L z [107] turn the Lyapunov orbit to the ISCO and the orbit becomes stable. For more information regarding the relation between the separatrix and last stable orbits, as well as other special last stable orbits, we refer the reader to [126]. Furthermore, if δ ≥ 1/2 then the geometry has both an ISCO and a photon sphere (PS), where unstable null geodesics accumulate [127,128], at [88] Notice that at the Schwarzschild limit where δ = 1, r PS = 3M and r ISCO = 6M as expected. Hereafter, we will focus on spacetime deformations that are larger than 1/2 in order to have an exotic central object with a PS and an ISCO that are fundamental hypersurfaces of astrophysically relevant (exotic) compact objects.

Geodesics
A first-order approximation to EMRI evolution can be accomplished through geodesics of a point-particle of mass μ which plays the role of the secondary orbiting around the primary supermassive compact object. The geodesic equations readẍ where κ λν are the Christoffel symbols of spacetime, x κ is the four-position vector, x κ is the four-velocity vector and the overdot denotes differentiation with respect to proper time τ .
Stationary and axisymmetric spacetimes, such as Eq. (2.1), possess metric tensor components that are t-and φ-independent. Therefore, they admit at least two conserved quantities (due to stationarity and axisymmetry) throughout geodesic evolution, namely the energy E and z-component of the orbital angular momentum L z The t-and φ-momenta can be expressed with respect to the conserved quantities and the non-zero metric tensor components. Together with the conservation of the rest mass μ of the secondary, (preservation of four-velocity) which leads to g λνẋ λẋ ν = −1, the geodesics of test particles present three constants of motion. Specifically, the conservation of the secondary's four-velocity leads to a constraint for bound orbitṡ where the effective potential V eff has the form The curve defined when V eff = 0. i.e. the curve of zero velocity (CZV), can be used in order to choose proper initial conditions that lead to bound orbits in the external vicinity of the primary. Bound geodesic motion in integrable systems can, generically, be characterized by three orbital frequencies. These frequencies are associated with the radial rate of transition between the periapsis and apoapsis of the geodesic (ω r ), the rate of longitudinal oscillations through the equatorial plane (ω θ ) and the revolution around the primary (ω φ ). Generic trajectories with irrational ratios of orbital frequencies span on two-dimensional tori and fill them completely. To the contrary, when the ratio of two orbital frequencies is a rational number then the geodesic is periodic (or resonant) and returns to its initial position after a number of oscillations defined by the multiplicity of the resonance [4]. Such orbits are special in the sense that they are not phase-space filling and therefore, can directly affect the evolution of EMRIs when encountered [22, 53-55, 111, 129-141].
Regardless of the fact that the ZV metric is non-integrable, when the deformation parameter δ is not too far from unity [4] then the majority of generic orbits are still characterized by orbital frequencies and elements in accord to the KAM theorem discussed below. Close to resonances, indirect and direct chaotic phenomena appear due to the non-integrability of spacetime [107], but we stay close enough to δ = 1 so that pure chaos never emanates in a direct manner (see Sec. 3.3 for more details on indirect and pure chaos). Therefore, we are still able to define orbital frequencies and only an almost zero-volume of the parameter space has regions of pure chaos (which we cannot spot in our analysis) where orbital frequencies are ill-defined. This will become more obvious in the following sections where all our imprints of chaos are indirect in nature.

Inspirals
To construct the inspiral trajectory we numerically integrate the coupled system of r , θ equations, after utilizing Eq. (3.4), augmented with post-Newtonian (PN) fluxes for the energy and angular momentum, respectively [48][49][50]142]. This treatment, though approximate, takes into account the dominant contribution of the secondary's radiative backreaction to the spacetime geometry, at second PN order, and results to an adiabatic evolution of the EMRI through time-dependent shifts onto the energy and z-component of angular momentum of the secondary. Since the inspiral evolves very slowly, the orbit is treated, in small timescales, as a geodesic, while for long timescales the trajectory is driven adiabatically through successively damped geodesics. This method, known as the hybrid kludge scheme, has been shown to perform very well when compared to Teukolsky-based Kerr waveforms for EMRIs [47].
Notice that the ZV metric has a non-trivial multipolar structure due to the deformation parameter δ, and in particular a non-zero mass quadrupole tensor, as opposed to that of Schwarzschild BHs. At second PN order, the kludge scheme [48] involves the mass quadrupole moment M 2 (for Kerr), thus to construct a more appropriate inspiral around ZV compact objects we apply a modification to the fluxes (see [51][52][53][54][55]111]) in order to include the quadrupole moment M 2 of the ZV metric, which represents the effect of δ on the evolution of E, L z , and set the spin parameter to zero. The adiabatic approximation, together with the flux augmentation, employed here has recently been found to provide results qualitatively equivalent to evolutions with instantaneous selfforce in non-Kerr electromagnetic analogues, which indicates that the methods we use can in principle describe resonance-crossings with sufficient accuracy in EMRIs [141]. Nevertheless, more accurate inspirals can be built by directly solving the wave equation resulting from metric perturbations and calculating the GW emission at the object and infinity, though this is a an almost impossible task and out of the scope of the phenomenology we are after in this article.
We assume linear variations of the momenta as in [54,55,143] where E 0 , L z,0 are the initial energy and z-component of the angular momentum, respectively, and d E/dt | 0 , d Lz/dt | 0 are the radiation fluxes calculated at the beginning of the inspiral, through the equations in [50,55]. T r is the time that the orbit takes to travel from the periapsis to apoapsis and back, while the fluxes (3.7) and (3.8) are updated every N r cycles for the whole EMRI evolution.

Detecting chaos
To understand the phase space structure of orbits around the ZV primary we can employ well-known tools in order to gain further intuition regarding orbital phenomena and chaotic imprints. A typical example is the Poincaré surface of section which is constructed by successive intersections of geodesics, with varying initial conditions, on a surface of section (here, we choose the equatorial plane) with strictly positive (or strictly negative) direction of intersection. The structure of the Poincaré map can instantly reveal if pure chaos is present, through disorganized intersections, or indirect imprints of non-integrability with the appearance of resonant islands, that encapsulate stable periodic points [4], and exist due to non-integrability, in accord to the Kolmogorov-Arnold-Moser (KAM) and Poincaré-Birkhoff theorems [144][145][146].
Since the ZV metric does not have a Carter-like constant these features are present in its orbital phase space [107]. More specifically, indirect chaos, which we deal with in this work, is connected to the appearance of islands of stability that surround resonances in Poincaré maps. This is due to the fact that KAM curves at resonances disintegrate onto a set of stable and another set of unstable periodic points, instead of forming typical KAM curves. Stable periodic points are encapsulated by islands of stability (resonant KAM curves), where the rotation number is shared through all geodesics residing in the island, while unstable periodic points are sources of chaotic orbits that shield the islands with extremely thin chaotic layers of disordered intersections that are practically indistinguishable in most nonintegrable cases of non-Kerr BHs. Direct or pure chaotic orbits, on the other hand, occupy a much larger and clearly distinguishable volume of phase space in Poincaré maps and correspond again to disordered intersections on a Poincaré map (see e.g. Figure 7 in Ref. [53] or Figs. 9-12 in Ref. [107]). In any case, pure chaos, or even thin chaotic layers around resonant islands, are not discernible in the Poincare map of Fig. 1 due to the small deformation parameter we utilized.
Another tool to detect chaos is the rotation number. We calculate it by tracking the angle ϑ between two successive intersections on the Poincaré map relative to the fixed central point of the map which corresponds to a circular orbit that intersects the surface of section exactly at the same point. The rotation number is defined as the accumulation of many angles measured between consecutive intersections as [4] for which when N → ∞ (with N the number of angles measured), Eq. (3.9) converges to the radial and polar orbital frequency ratio ν ϑ = ω r /ω θ . Calculating consecutive rotation numbers for different geodesics, by smoothly varying one of the parameters or initial conditions of the system while keeping the rest fixed, leads to a rotation curve. Integrable systems demonstrate monotonous rotation curves. On the other hand, nonintegrable systems possess transient plateaus with a non-zero width when geodesics occupy resonant islands. This designates a crucial aspect of resonant islands, that is when a geodesic is inside the island it shares the same rational ratio ω r /ω θ with the stable periodic point which leads to the plateau formation; a phenomenon that does not appear in integrable systems whatsoever, even though resonances still exist but occupy only a single point in phase space. Inflection points also appear when trajectories pass through the intersection of resonant islands where unstable periodic points reside and chaotic layers, that surround resonant islands, emanate. Nevertheless, by changing the initial conditions the orbit can be driven through the island and give rise to a typical plateau. Rotation curves are not only a tool that is used for geodesics but can also be employed in dissipative scenarios, where the rotation number evolves with respect to time. In the case of an EMRI, one can use selected timesteps of the inspiral as initial conditions for a geodesic evolution. Through the non-dissipative trajectory, the Poincaré map and eventually the rotation number of each timestep can be calculated in order to plot a series of rotation numbers as the EMRI evolves with time. The same attributes hold here as well, namely monotonous dissipative rotation curves for integrable systems and appearance of plateaus for non-integrable EMRIs. For more information regarding integrability and chaos in EMRIs see the following series of works and references therein [22,[53][54][55].

Geodesic and inspiral evolution
By solving the coupled r , θ second-order ordinary differential equations, together with the first-order decoupled equations for t and φ from Eq. (3.4), we obtain bound orbits that reside inside CZVs and never plunge nor escape to infinity. To check the precision of our geodesics we evolve the constraint equation (3.5) for ∼ 10 4 revolutions and find that it is satisfied within one part in 10 8 − 10 10 depending on the initial conditions and deformation of spacetime.
To guarantee numerical accuracy for inspirals, we calculate the 4-velocity in each update of the fluxes and check its conservation along a geodesic evolution with initial conditions the energy, z-component of angular momentum, position and velocity at every update timestep. For all simulations presented hereafter, the constraint is satisfied to within a part in ∼ 10 8 for the first 10 4 crossings through the equatorial plane.
As a first step, we reproduce some qualitative features of the ZV metric, namely its non-integrability [107], which will be later used to choose initial conditions in  [22] for a zoom into the encapsulated structure of these islands). With a pedantic search on the available parameter space, we can easily spot four resonant islands with different multiplicity, which we designate in Fig. 1 with colored curves. Notice how those do not surround the central point, but rather the stable points in phase space where the true resonant orbits emanate, and that the multiplicity defined by the denominator of the resonance (the longitudinal oscillations through the equatorial plane) corresponds to the number of islands. It is noteworthy that the resonances appearing here, besides the 1/2-resonant island, are not those that strongly affect an inspiral, especially when assessing their impact in EMRI parameter estimation [129,137,138], yet we will later see that they can also significantly contribute to it if the spacetime is non-integrable, mainly because they accumulate before plunge.
The rotation curves presented in Fig. 2 promote the previous discussion perfectly. Various plateaus and inflections appear right where we expect the rotation number to have a rational ratio. The plateaus clearly designate that any geodesic residing in a resonant island shares the same rational ratio of orbital frequencies with the center of the island, where stable periodic orbits emanate. Even when an inflection point shows up at the rotation curve, meaning that the geodesic lies between the tips of two resonant islands where unstable periodic orbits exist, with a certain change in the initial velocityṙ (0), we can access the island as shown in Fig. 3. For   For the rest of the discussion we will focus on oblate deformations for the following reasons: (i) oblate deformations are usually the ones enabling the strongest effects of non-integrability and chaos (see [22, 52-55, 112, 140]), (ii) as seen in Fig. 2, an oblate deformation drives the resonances closer to the central object thus we expect amplified chaotic effects and (iii) δ > 1 is an optimal choice to imitate BHs since for these cases the ZV object possesses a PS, an ISCO and produces similar shadows [81].
So far, the discussion involved zeroth-order approximations of EMRIs since there was no radiation loss. Turning on the fluxes (rates of change of orbital energy and angular momentum) (3.7) and (3.8) leads to a numerical integration that is more intricate but the results are very interesting since the secondary inspirals adiabatically towards the primary due to fluxes through GWs. Figure 4 presents a typical inspiral of a secondary that transverses the 1/3-resonant island. The initial conditions are not fine tuned so what is presented here is a generic feature. For this inspiral, we have updated Since the revolution period in this case is ∼ 90M, the total evolution time of the EMRI is t total = 1.6 × 10 5 M ∼ 1800 revolutions (roughly 9 days for the mass of the primary considered). At a certain time, the fluxes of energy and angular momentum become more negative abruptly. This instant of time designates the entry of the secondary into the 1/3 plateau of the island (see right panel in Fig. 4). Right after the entry, the rate of change of orbital energy reaches a local maximum, at a special instant of time which is shown with a vertical dashed line, beyond which it continues to decline due to the inspiral. This behavior is not at all peculiar, but rather has a true physical meaning. Regardless of the fact that the secondary transverses the island while sharing the same rotation number throughout it, the secondary only reaches the stable periodic point at the center of the island at a certain time designated with the vertical dashed line. This is when the orbit becomes exactly periodic. Periodic orbits are the closest to circular ones. Since circular orbits emit monochromatic radiation at a single frequency (twice the revolution frequency) their energy emission is minimized (thus the rates of change of energy and angular momentum are also minimized (in absolute value)). Since a resonance emits quasi-monochromatic radiation, at the time that the inspiral passes through the perfect resonance, the fluxes are also maximized. The exit from the island takes place when the flux value at the entry is met again and the flux curve changes gradient (most obviously shown later in Fig. 6). The secondary can spend ∼ 100 cycles in perfect resonance which translates to roughly half a day or ∼ 5% of the whole EMRI evolution. The picture is qualitatively the same for the flux of angular momentum. Figure 5 presents an even more intriguing EMRI evolution around the ZV primary, that undergoes two consecutive island crossings during a single evolution, namely the 2/7 and 1/4 islands. Initial conditions are slightly (but not strongly) fine tuned so what is presented here is quasi-generic, in a sense that it can occur for a small but non-zero set of initial conditions. For this inspiral, we have updated the fluxes N r = 1100 times, every T r = 200M. Since the revolution period in this case is ∼ 80M, the total evolution time of the EMRI is t total = 2.2 × 10 5 M ∼ 2800 revolutions (roughly 12 days for the mass of the primary considered). The time that the secondary spends in the first island, in perfect resonance, is ∼ 200 cycles which translates to roughly a day. After exiting the first island, the orbit enters shortly after a subsequent resonant island and occupies it for another ∼ 50 cycles which translate to one fifth of a day. In total the inspirals experiences 250 cycles of perfect resonance (without taking into account preand post-resonant effects) which correspond to ∼ 10% of the whole evolution. Similar analyses for non-Kerr EMRIs with the same mass ratios used here (μ/M = 10 −6 ) experience 250 cycles in the most prominent 2/3-resonant island [54,55], though due to the much slower inspiral (resonances appear further from the central object) this only corresponds to ∼ 3% of the whole evolution spent in a single island. Practically, the fact that a multitude of islands gather close to the plunge gives us the ability to probe subdominant resonances consecutively, in a short period of time, and have an EMRI that experiences a significant fraction of its evolution in resonance. Zooming in on the left plot of Fig. 5 (see Fig. 6), we observe a complete agreement with what is demonstrated in Fig. 4. Especially for the 2/7-resonant island (left plot in Fig. 6), the initial condition chosen here allows for the secondary to remain in it for a substantial amount of time that leads to a significant maximization of the energy flux. This means that the orbit enters deep into the island and probes the stable periodic point for a significant amount of time till it exits. The initial condition may look special, since a typical crossing would not maximize fluxes considerably, nevertheless, plateaus on dissipative rotation curves are still generic for a wide range of initial conditions. In this particular subfigure, one can also locate not only the entry in the island (with the sudden drop of the energy flux discussed above) but also the exit where a change of the gradient occurs. The same practically happens in all fluxes when an EMRI transverses an island of stability but only very prominent and long plateaus can reveal evidence of not only the entry but also the exit from the resonance.

Gravitational waves
In this section we approximate the dominant GW emission of an inspiraling stellarmass secondary around an oblate ZV supermassive compact object and search for imprints of non-integrability in the waveforms produced by such EMRIs when detected by a LISA-like interferometer. For the GW modeling, we use the quadrupole approximation described below.

Quadrupole formula
The radiative component of the metric perturbation introduced by the secondary at luminosity distance d from the source T can be read at the transverse and traceless gauge as where Q i j is the symmetric and trace-free (STF) quadrupole tensor with t being the coordinate time measured at very large distances from the detector. The source term of the secondary (which is treated as a point particle through a delta function) is where Z(t) = (x(t), y(t), z(t)) the position vector in pseudo-Cartesian coordinates and the trajectory components with respect to flat spherical coordinates, under the assumption that our space-borne detector is positioned at infinity. Even though we have identified the Schwarzschild-like coordinates (r , θ, φ) of the secondary's trajectory with flat-space coordinates, known in the literature as the "particle-on-a-string" approximation, and we assume a finite luminosity distance d from the source, such prescription is not strictly valid. Nevertheless, it has been found to work very well when generating EMRI waveforms in GR [47]. GWs can be projected onto two polarizations, + and ×, with the introduction of two unit vectors, p and q, which are defined with respect to a third unit vector n that points from the source to the detector. The triplet of unit vectors ( p, q, n) is chosen so that they form an orthonormal basis. The polarization tensor components read i j and allow for the metric perturbation to be written as The GW polarization components can then be described in terms of the position, Z i (t), velocity, v i (t) = d Z i /dt, and acceleration a i (t) = d 2 Z i /dt 2 vectors as [143] h +, LISA's response to an incident GW is rather complicated and depends on the antennae response patterns (see [142,147] for the full equations). Here we assume a detector that lies at a luminosity distance d with orientation n = (0, 0, 1) with respect to the source, and utilize where α = (I , I I ) is an index representing the different antenna pattern functions F + α , F × α which can be found in Refs. [142,147,148]. For phenomenological purposes we will use a single-channel approximation, that is set α = I in Eq. (5.11), since it is enough to accommodate the fundamental parts of gravitational radiation emitted by the EMRI, i.e. its phasing.

Gravitational-wave frequency evolution and cumulative glitches
To comprehend how a ZV EMRI imprints non-integrable effects, such as plateaus, onto its emitted waveform, we calculate the GW in the Einstein-quadrupole approximation for the particular inspiral outlined in Sec. 4 that crosses two consecutive islands, namely the 2/7 and 1/4-resonances. We focus solely on this example in order to demonstrate that even subdominant resonances affect EMRI evolution, when the spacetime presents chaotic features.
We perform a Fourier transform on the extracted waveform from the inspiral and plot a density spectrogram that displays the evolution of one of the harmonics with respect to time. We use the same methodology for the spectrogram as in [54,55], i.e. we cut the waveform in time segments with a particular window size and offset and perform consecutive Fourier transforms, in order to overcome the uncertainty between frequency and time resolution. Here, since the inspiral is evolving fast and the island crossings have different timescales, we use distinct window size and offsets for the chunks of evolution around the two resonant islands. This is necessary for demonstration purposes since the first resonance is crossed much slower than the second one.
In Fig. 7 we show the spectrogram of the GW extracted from the EMRI that successively transverses two resonant islands of stability. In both cases, when the island is met the waveform frequency evolution loses monotonicity. The resonant crossings  Fig. 7) or a rapid glitch 2 (right plot in Fig. 7). Both instances designate that the system is non-integrable, especially the encounter of the first island that shares a striking similarity with plateaus appearing in rotation curves. However, this is an outcome of short-term occupancy in the islands and altogether a smaller deformation of spacetime, with respect to those presented in non-Kerr EMRIs [54,55]. If the inspiral is given more time in the island, the manifestations should resemble more those of discontinuous GW glitches. Yet, the fact that resonant island crossings are imprinted in the GW of a ZV EMRI, in a short period of time, can serve as a 'smoking gun' for a non-integrable BH mimicker in the center of our galaxy, since these phenomena do not appear in Schwarzschild EMRIs that have monotonous rotation curves and spectrograms.

Concluding remarks
Probing the spacetime around supermassive compact objects with EMRIs is one of the prime targets for space-borne detectors like LISA. Compact objects are either spinning and/or surrounded by astrophysical environments, especially those residing in active galactic nuclei. Therefore, spherical symmetry is rather fragile and is usually broken. When there are not enough spacetime symmetries to guarantee the integrability of geodesics around these objects, then LISA may be able to detect particular phenomenological imprints of such manifestations like GW glitches around expected transient orbital resonances [54,55] that differ significantly from instrumental glitches [149].
We showed that deformations from spherical symmetry which keep the spacetime static, described by the ZV geometry, have the potential to exhibit GW glitches when involved as primary objects in EMRIs, due to the non-integrability of test-particle dynamics. These phenomena not only appear for single resonant island crossings, such as those occurring in non-Kerr EMRIs, but also are present in a cumulative and short-timescale manner, where the secondary transverses two (and possibly more) subdominant resonant islands. Note though that this can also occur in non-Kerr EMRIs if long-lasting evolutions are to be performed, though the effect of subdominant resonances are usually suppressed [54,55].
When supermassive compact objects, such as those residing in galactic centers, have multiple interpretations due to the plethora of objects that can cast similar shadows, GW astronomy with LISA can, in principle, distinguish between models that are integrable or not through the detection of glitches in gravitational waveforms. Synergies between shadow observations and space-borne detectors can, therefore, narrow down the parameter space of solutions describing supermassive objects in galactic centers, and in particular can emphatically decide if M87* and Sgr A* are described by Kerr and Schwarzschild geometries, respectively, unless the EMRI includes multiple [57,150] or spinning [139,151] secondaries.
Our analysis only deals with the phenomenological imprints of non-integrability. Nevertheless, if one wants to utilize the aforementioned phenomenology in practice, a consistent glitch modeling analysis for non-integrable EMRIs should be carried out, in a systematic way in order to understand, first, if these phenomena are clearly detectable with space interferometers, second, to which extent they affect parameter estimation, and third, to which degree these effects differ from standard transient resonances experienced by integrable EMRIs which are already sufficiently modeled with PN techniques [137,138] and gravitational self-force [129,130,136,152,153].
Author Contributions All authors have contributed equally for the execution of calculations and the presentation of results, as well as the creation of this manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL. This work was supported by the DAAD program for the "promotion of the exchange and scientific cooperation between Greece and Germany IKYDAAD 2022" (57628320). K.D. and K.D.K. are grateful for hospitality provided by the Section of Astrophysics, Astronomy, and Mechanics, Department of Physics of the National and Kapodistrian University of Athens, Panepistimiopolis Zografos GR15783, Athens, Greece. K.D. acknowledges financial support provided under the European Union's H2020 ERC, Starting Grant Agreement No. DarkGRA-757480 and the MIUR PRIN and FARE programmes (GW-NEXT, CUP: B84I20000100001).

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflict of interest
We the authors hereby declare that there are no competing interests of financial or personal nature.
Ethical approval Ethical approval for this work is not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.