Assessing the impact of transient orbital resonances

One of the primary sources for the future space-based gravitational wave detector, the Laser Interferometer Space Antenna, are the inspirals of small compact objects into massive black holes in the centres of galaxies. The gravitational waveforms from such Extreme Mass Ratio Inspiral (EMRI) systems will provide measurements of their parameters with unprecedented precision, but only if the waveforms are accurately modeled. Here we explore the impact of transient orbital resonances which occur when the ratio of radial and polar frequencies is a rational number. We introduce a new Effective Resonance Model, which is an extension of the numerical kludge EMRI waveform approximation to include the effect of resonances, and use it to explore the impact of resonances on EMRI parameter estimation. For one-year inspirals, we find that the few cycle dephasings induced by 3:2 resonances can lead to systematic errors in parameter estimates, that are up to several times the typical measurement precision of the parameters. The bias is greatest for 3:2 resonances that occur closer to the central black hole. By regarding them as unknown model parameters, we further show that observations will be able to constrain the size of the changes in the orbital parameters across the resonance to a relative precision of 10% for a typical one-year EMRI observation with signal-to-noise ratio of 20. Such measurements can be regarded as tests of fundamental physics, by comparing the measured changes to those predicted in general relativity.


I. INTRODUCTION
Gravity is the weakest of the fundamental interactions and is particularly challenging to study in a laboratory experiment. Fortunately, the Universe has plenty of gravitational phenomena which can be observed and used to improve our understanding. For hundreds of years we have been studying the cosmos through the observation of electromagnetic radiation, but in the last few years it finally became possible to also listen to the Universe through the observation of gravitational waves (GWs) [1].
In fact, these ripples in the fabric of spacetime allow us to study gravity and the Universe from a totally new perspective. Indirect evidence for gravitational wave emission was first observed in 1974 by Hulse and Taylor in a study of the orbital decay of a binary pulsar [2,3]. The first direct observation of a gravitational wave signal was made in September 2015 [1], by the ground based laser interferometer detectors LIGO [4]. This and subsequent detections with the LIGO and Virgo detectors [5] have opened up a new era of gravitational wave observations [6][7][8][9][10][11][12][13][14][15][16], which have permitted studies of the properties of the astrophysical population of compact objects, such as black holes and neutron stars, and have led to some of the most stringent tests of General Relativity.
The ground based gravitational wave detectors have a frequency range from about 10Hz to a few kHz, which permits observations of merging binaries of stellar mass compact objects at low redshift, z < 2. The European Space Agency plans a new space-based gravitational wave observatory, the Laser Interferometer Space Antenna (LISA), to open a gravitational window on the Uni- * lorenzo.speri@aei.mpg.de verse in the low-frequency range of 10 −4 −10 −1 Hz. LISA is expected to resolve thousands of overlapping gravitational wave signals from a wide variety of sources [17].
The inspiral of a stellar-origin compact object (CO) into a massive black hole (MBH) in the centre of a galaxy is one of the potential sources of gravitational waves for LISA. Observations of such Extreme Mass Ratio Inspiral (EMRI) systems have a huge scientific potential. The compact object typically completes 10 4 − 10 5 cycles in band, during which time it is orbiting in the strong field region close to the central black hole. Therefore, EMRI signals encode a detailed map of the background spacetime of the central MBH and offer a unique opportunity to measure the properties, evolution and environment of MBHs [18][19][20][21][22][23], to test for deviations from General Relativity (GR) [24,25] and to constrain cosmological parameters [26,27]. We expect to observe between a few and a few hundred EMRIs over the LISA mission duration [28][29][30][31], with the large uncertainty driven by uncertainties in the rates of the various dynamical processes that can lead to the capture of compact objects by a MBH [23,32].
In order to realize the scientific potential of these astrophysical sources we have to correctly detect and characterize EMRI signals. A sequence of LISA data challenges showed the feasibility of identifying individual EMRIs in data sets without other sources, and using narrow parameter priors [33][34][35]. However, any EMRI data analysis technique will rely on the existence of waveforms that are both accurate enough to faithfully identify astrophysical EMRI signals and that are quick enough to generate in the large number required to search the wide parameter space of potential signals (see for example recent work in [36,37]).
Accurate modeling of EMRI waveforms can be accom-arXiv:2103.06306v2 [gr-qc] 11 Jun 2021 plished within the framework of black-hole perturbation theory, regarding the mass ratio η ∼ 10 −4 − 10 −7 as a small expansion parameter. Within this framework, the evolution of the orbit of the compact object is governed by the Kerr geodesic equations with a forcing term, called the gravitational self-force, which takes into account the finite size and mass of the body and its back reaction on the background Kerr spacetime [38,39]. Generic EMRI orbits have three fundamental frequencies: the radial frequency ω r , the polar frequency ω θ and the azimuthal frequency ω φ . When the radial and polar frequencies become rational multiples of each other, a transient orbital resonance occurs and the motion of the small body is restricted to a subspace within the full orbital torus [40]. The standard adiabatic approximation to compute EMRI phase space trajectories, which involves averaging over the full torus of the orbit in phase space, is then no longer valid. Instead, additional terms in the evolution equations must be taken into account, which cause the system to exhibit a qualitatively different behavior, which has been proven to affect the detectability of EMRIs [41].
Our work aims to investigate two questions related to EMRI parameter estimation in the presence of resonances. The first is to explore the size of parameter biases that arise when the presence of resonances is ignored while building a waveform model. The second is to explore with what precision the effect of the resonance, i.e., the amount by which the constants describing the orbit change over the resonance, can be measured from an observation. EMRI resonances have been studied in a number of previous works [42][43][44][45][46][47][48][49][50][51], but this paper is the first attempt to explore resonant EMRI parameter estimation.
We achieve this goal using an Effective Resonance Model. This is a phenomenological model that adds transient orbital resonances into the Numerical Kludge (NK) model framework, which is an approximate waveform model that can quickly generate generic EMRI inspirals and compute their GW signals [52,53].
This work is organized as follows: in Section II we review the general approach to EMRI modeling, we briefly describe the Numerical Kludge model and we present the GW data analysis tools we employ. Subsequently, in Section III, we review the theoretical understanding of transient orbital resonances which will be used to construct the Effective Resonance Model in Section III B. Finally, in Section IV, we show the results of our investigation using the new resonance phenomenological models. We compute the dephasings induced by 3:2 resonances over the parameter space and assess the corresponding parameter biases that would be induced by ignoring resonances, and we explore the measurability of resonance effects over a wide range of resonance strengths.
Notation Throughout, we will use a spacelike signature (−, +, +, +) for the metric and work with a system of geometrized units in which G = c = 1. We use the Einstein summation convention for repeated indices. Greek letters will be used to indicate a sum over all spacetime indices, whereas latin letters will be used to indicate a sum over spatial only.

II. EMRI MODELING
In this section we give a general overview of EMRI waveform modeling, and highlight how this should be modified to include the effects of resonances. We will also summarise how the Numerical Kludge model is constructed. Lastly, we will review some tools from gravitational wave data analysis that we will use to assess the impact of resonances on gravitational waves from EMRIs.

A. From Kerr geodesics to adiabatic evolution
Kerr geodesics are a good approximation to the evolution of EMRIs on short time scales (for order of several orbital periods). Their properties are also a useful basis for understanding the motion of the compact object at higher orders in the mass ratio. From the symmetries of the Kerr space-time, it is possible to identify four integrals of the geodesic equations of motion: the energy, E, the axial angular momentum, L z , Carter's constant, Q and the mass of the compact object, µ. The first two are direct consequences of the stationarity and axial symmetry of Kerr space-time. Carter's constant is a conserved quantity associated with a tensorial Killing vector field [54] and, in the weak field or non-spinning limits, it corresponds to the squared angular momentum component parallel to the equatorial plane (Q ≈ L 2 x + L 2 y ). The mass is conserved because of the normalization of the four-momentum g νµ Kerr p ν p µ = −µ 2 , i.e., because of the conservation of the Hamiltonian H = g νµ Kerr p ν p µ /2. The existence of four constants of motion J β = (E, L z , Q, µ) allows us to reduce the geodesic equations from four second-order differential equations to four first-order equations. In Boyer-Lindquist coordinates, (t, r, θ, φ), the first-order equations are: where Σ(r, θ) = r 2 + a 2 M 2 cos 2 θ, ∆ = r 2 − 2M r + a 2 M 2 and M and a are the Kerr black hole mass and dimensionless spin parameter, respectively, and τ is the proper time. By using the Carter-Mino time parameter λ defined via dτ /dλ = Σ [54,55], the equations in r and θ decouple and the radial and polar motion can be determined if constants and initial conditions are specified. The trajectory of EMRIs can be modelled with Kerr geodesics only on a short time scale, of the order of the orbital time scale ∼ O(1), compared to the inspiral time scale ∼ O(1/η). If we want to describe the long inspiral of EMRIs, it is necessary to take into account the impact of the gravitational field of the compact object on the background Kerr space-time. During the inspiral, the compact object slowly deviates from geodesic orbits and this can be interpreted as an effective acceleration/force due to the so-called Gravitational Self-Force (GSF) acting on the compact object [56].
By using the Action-Angle formalism it is possible to rewrite the geodesic equations to include the effect of the first order GSF as [56,57]: where the evolution of the Boyer-Lindquist coordinates is described using the angle variables q α = (q t , q r , q θ , q φ ) [58] and the corrections due to the GSF are encoded in G and g. We note that here the action variables, J α , are not precisely the same as the standard constants of motion introduced earlier, J β , but are related by a transformation. This distinction is not important for the present discussion.
At "zeroth" order in η the motion is completely determined by the fundamental frequencies ω α , that define the fundamental modes of the evolution of the system and reduce to the usual Keplerian frequencies in the Newtonian limit. Expressions for these quantities can be found in [59,60]. The radial fundamental frequency ω r is associated with the radial motion and it is zero for circular orbits. The polar and azimuthal fundamental frequencies, ω θ and ω φ , are associated to precessional motion of the orbit.
The four constants of motion J α determine the "shape" of the orbit, whereas the phase variables give information on the (time dependent) location of the object within the orbit, and the orbit's orientation [56]. The Action-Angle formalism not only makes the periodicities of the system obvious, but it also facilitates the decomposition of any dynamical field as an expansion in the fundamental frequencies.
The equations (2a-2b) can be further simplified by using the two-time scale expansion developed by Flanagan and Hinderer [57], which consists of separating the evolution of the EMRI into two time scales, a long time scale J/J ∼ O(1/η) associated with the evolution of the constants of motion, and a short time scale q/ω ∼ O(1) associated with the evolution of the angle variables. By taking an average of the flux equation (2b) it is possible to simplify the evolution of the constants of motion: (3) where we have used an expansion of G α (q r , q θ , J α ) as a Fourier series and have introduced the 2-torus average over the phase variables, defined by For ergodic trajectories, i.e., phase-space filling trajectories, the phase variables evolve rapidly and the exponential term averages to zero e i(lqr+mq θ ) qr q θ ≈ 0 [57].
The adiabatic approximation consists of using the averaged equations for the evolution of the constants of motions. In this way, the evolution of averaged constants can be found independently of the evolution of the phases. The resulting J α (t) can then be used to solve equation (2a) for the phase variables by also dropping the oscillating terms g.
However, if the argument of the exponential does not rapidly oscillate, the average cannot be taken. As a consequence, the adiabatic approximation breaks down and it is necessary to take into account another secular term. This is precisely what happens when resonances occur, and we will return to this in Section III.

B. Numerical Kludge Model for EMRI Waveforms
The Numerical Kludge (NK) [52] is one of a number of fast but approximate EMRI models that have been developed to substitute the computationally expensive GSF calculations when exploring EMRI data analysis. Other models include the Analytical Kludge (AK) [18], the Augmented Analytical Kludge (AAK) [36,37,61] and the Near Identity Transform (NIT) [62].
The Numerical Kludge (NK) generates a waveform for a one-year generic EMRI inspiral, sampled with a 10s time step, in about one minute, and its implementation is well suited to include resonances. Here we briefly summarise how the NK is constructed. The model uses a Keplerian parameterization of the orbit to simplify the treatment of turning points in the r and θ motion [59,63] Here r p and r a denote the turning points of the orbits at periapsis and apoapsis, and we have introduced the eccentricity, e, and semi-latus rectum, p, which are defined in terms of the turning points through r p = p/(1 + e), r a = p/(1 − e). We denote the relativistic anomaly by ψ and the "polar phase" by χ. The latter two are 2π periodic but they are not canonical phase variables. If we replace the Carter constant by an inclination angle 1 defined by tan ι = √ Q/L z , the motion can be described in terms of (p, e, ι) instead of (E, L z , Q) and the flux evolution becomes: dJ dt = f NK (a, M, µ, p, e, ι) .
The computation of the gravitational waveform proceeds via the following steps: (i) The parameters defining the waveform can be grouped in the following sets: intrinsic parameters extrinsic parameters: where the two sky-position angles (θ S , φ S ) are the co-latitude and azimuth in an ecliptic based coordinate system, the two angles (θ K , φ K ) define the direction of the spin of the massive black hole, with respect to the same coordinate system as the sky position 2 , φ 0 is the azimuthal initial phase and D L [Mpc] is the luminosity distance.
the length of the inspiral t max , and the sampling interval ∆t.
(ii) The evolution of the constants of motion J = (E, L z , Q) is computed using the adiabatic approximation, i.e., replacing the first order radiative selfforce G with its 2-torus averaged value f NK = G qr q θ , and neglecting the conservative piece of the GSF. The two torus average is approximated with the dissipative part of the first order self-force, obtained from second order post-Newtonian formulae and fits of Teukolsky-based inspirals [53,64]. We will have to modify this flux evolution equation to take into account resonances.
(iii) The inspiral trajectory is calculated by integrating the Kerr geodesic equations, evaluated for the instantaneous constants of the motion computed from the evolution equations and dropping all the forcing terms g α in Eq. (2a). Then, we identify the Boyer-Lindquist coordinates (r BL , θ BL , φ BL ) → (r sph , θ sph , φ sph ) with a set of flat-space spherical polar coordinates. By doing this, the particle is effectively forced to move along a curved path in flat spacetime.
(iv) The gravitational waveform is then calculated by applying the quadrupole formula in the TT gauge to the pseudo-flat-space particle orbit, obtaining h + , h × , which can then be projected into the observer direction.
(v) The low-frequency approximation to the LISA response function, described in [18,65], is then used to transform the waveform polarizations h +,× into the two LISA response functions h I and h II .
Despite the several approximations, the NK waveforms have been shown to be remarkably faithful, for periastron r p 5M , when compared to Teukolsky-based inspirals [52]. The NK procedure works well because it takes into account the most important physics, allowing phenomenological waveforms for generic EMRIs to be quickly generated without losing too much faithfulness [52].

C. Waveform analysis methods
Detection and parameter estimation of gravitational wave signals employs a wide range of techniques used to identify and characterize one or more signals, h(t), present in the output, s(t), of a gravitational wave detector. We assume that the output of a gravitational wave detector is composed of a signal h(t; λ), dependent on the source parameters, λ, and instrumental noise n(t): By assuming that the noise is weakly stationary, Gaussian and ergodic with zero mean, we can write the likelihood for the parameters λ as [66]: where we have introduced the inner product · | · The tilde indicates the continuous Fourier transform and S n (f ) is the one-sided spectral density of noise in the detector. From a practical point of view the spectral density represents our information on the detector sensitivity. In this work we adopt the LISA PSD described in [67]. Given a waveform h, the optimal matched filtering signal-to-noise ratio (SNR) is defined by This is a measure of the detectability of a particular signal. The precision with which observations will be able to determine the parameters of a system can be estimated using the Fisher Information Matrix, which is where l is the log-likelihood, E is the expected value, ∂ i ≡ ∂/∂λ i and to obtain the second line we we have used the likelihood given in eq. (7). Formally, the Fisher Matrix provides a lower bound on the variance of any unbiased estimatorλ of the parameters of the signal, but it is also provides a Gaussian approximation to the shape of the likelihood, valid in the high SNR limit. The squareroots of the diagonal elements of the inverse Fisher Matrix thus provide an estimate of the precision with which the corresponding parameter can be measured in an ob- If we use an approximate waveform model h m (λ) to estimate the parameters λ 0 of a signal actually described by a model h t (λ), the recovered parameters will be affected by systematic errors. In the linear signal approximation, the shift in the peak of the likelihood due to statistical and systematic errors is given by [68]: The statistical error is determined by the noise realisation and scales with the SNR as δλ stat ∼ 1/ρ. By contrast, the systematic error does not depend on the SNR and encodes the bias on the recovered parameters that arises from mis-modeling. If the bias from ignoring a resonance is smaller than the statistical error, then the effects of the resonance are negligible.
The inner product can be used to define an overlap O(a, b) in the usual way: This expresses how "similar" two signals a(t) and b(t) are. If two signals are identical then the overlap is 1. We define also the Mismatch M as M = |1 − O| We will often compare two kinds of waveforms: resonant waveforms which are signals produced with the NK including the Effective Resonance Model and nonresonant waveforms which are produced with the standard NK model. We will focus on a specific subset: of the full parameter space, where the resonance coefficients C = (E, L z , Q) will be defined in the next section. These are new parameters of the waveform, which are set to zero for non-resonant waveforms, i.e., if we want to ignore resonances. Note that we are not including the phase angles or extrinsic parameters in this list. This is for computational convenience, as the accurate evaluation of a Fisher Matrix on the full EMRI parameter space is very challenging. While the phase angles will correlate with other parameters, we expect our conclusions to be minimally impacted by ignoring them. Passing through a resonance leads to a change in the orbital frequencies relative to a non-resonant trajectory, and hence a growing phase discrepancy after the resonance. This growing phase difference is likely to dominate over any small phase differences that accumulate during the resonance itself, and will not be degenerate with phase offsets encoded in ψ 0 and χ 0 .
In all the studies reported here we have set the noise realisation, n = 0, since our focus is on the systematic biases rather than the statistical ones. We also set to zero the initial phases and choose the extrinsic parameters to be (θ K , φ K , θ S , φ S , D L ) = (π/8, 0, π/4, 0, 200 Mpc). After computing a waveform with these parameters, we then rescale the distance in order to fix the signal-tonoise ratio to 20. This is thought to be roughly the detection threshold for EMRIs and hence will be a fairly typical value for observed systems. For the remainder of this study all events will have SNR of 20. For further technical details and a report of tests used to validate our Fisher Matrix calculations, we refer the reader to the appendices V A and V B.

III. TRANSIENT ORBITAL RESONANCES
Resonances are a common phenomenon in nature in which there is a change in the evolution of some physical quantity due to coherent superposition of a forcing term. When two frequencies of a system become commensurate, i.e., their ratio is a rational number, a resonance can occur and induce a distinctive change in the evolution. Flanagan and Hinderer found that transient orbital resonances occur in EMRIs when the radial frequency ω r and the polar frequency ω θ become commensurate [40]. During a resonance, the phase-space trajectory is not space-filling and therefore the averaging procedure of the adiabatic approximation breaks down. This leads to a deviation of the evolution of the orbital constants from the adiabatic trajectory. The phase evolution also shifts away from the standard evolution leading to an overall dephasing in the emitted gravitational waves. To study the impact of these resonances, it is necessary to build a model that takes them into account and that can be quickly used to study their impact. To this end, we review in the next sections the main features of transient orbital resonances in EMRIs and use them as a basis to construct an Effective Resonance Model. This phenomenological model for resonances is then implemented in the NK model framework and used to investigate the impact of these phenomena on parameter estimation.

A. Properties
The evolution of the constants of motion (eq. (2b)) can be decomposed in a Fourier series in the radial and polar angle variables, q r and q θ . Using this expansion, the change in the fluxes is manifestly expressed by two terms: the secular term G 0,0 and the rapidly oscillating term. The latter one can be averaged out using the adiabatic approximation for ergodic trajectories. However, if the phase Ψ = lq r + mq θ of the oscillatory term is stationary, the adiabatic approximation is no longer valid and additional contributions add up to the secular term G 0,0 . By expanding the phase variables in terms of the fundamental frequencies, we can understand when this occurs: The phase Ψ becomes stationary in the proximity of a resonance, i.e., when the resonance condition: l * ω r0 + m * ω θ0 = 0 is satisfied around τ = τ 0 for a given Fourier mode l * , m * ∈ Z. Therefore, when the ratio of the radial and polar fundamental frequencies, m * /l * , is a rational number, the respective term G l * ,m * and all the higher harmonic terms s l * , s m * , for s ∈ Z, (i.e., integer multiples of l * , m * ) are also approximately secular and the adiabatic approximation breaks down. In fact, during a resonance, the phases q r , q θ cover only a specific non-ergodic trajectory and the phase-space average of the normally oscillating term is non zero. The flux equation (3) during a resonance becomes: The resonance ends when the phase Ψ starts oscillating again, therefore, when the third term in equation (10) is significantly different from zero. Thus we can define the duration of the resonance as: We want to stress that there is not a sharp transition between the resonance regime and adiabatic regimes but a smooth one and, therefore, the resonance duration is somehow an arbitrary definition up to overall factors. Since the resonance condition could be satisfied for every pair of integers (l, m), one might think that an EMRI is basically always on resonance. However, only the low order resonances are significant. In fact, if the resonance order is high it means that the number of cycles in r with respect to the cycles in θ is high enough that the average over the 2-torus is effectively done. Furthermore, the duration of the resonance scales as the inverse of the order of the resonance |l| + |m| and only low order resonances are long enough to affect the evolution.
The fundamental frequencies and the constants of motion change typically on a time scale of J/J ∼ ω/ω ∼ 1/η. When an EMRI system passes through a resonance, the evolution of the constants of motion deviates from the standard adiabatic evolution. This deviation is caused by the additional contributions to the fluxes in Eq. (11) and its size is of order ∆J ∼Jτ res ∼ √ η, and analogously for the fundamental frequencies ∆ω ∼ √ η. Therefore, the evolution of the constants of motion receives a "kick", which can be either positive or negative depending on the phase at resonance Ψ 0 and on G l * ,m * . This deviation is bigger than any higher order corrections in the mass ratio. The exact expression for the total change in the constants of motion after resonance l * :m * can be found in [40,41]. The jumps in the constants cause the phases to shift away from the standard adiabatic evolution. This leads to a cumulative dephasing of order ∆q = ∆ω T ∼ 1/ √ η by the end of the inspiral, where the inspiral timescale T ∼ 1/η. These shifts of the phases can be particularly large in small mass ratio systems, which motivates investigating their effect on EMRI parameter estimation.
To do so we need to build a model for resonances.

B. Effective Resonance Model
Here we develop a phenomenological model of resonances based on the main features of transient orbital resonances discussed in the previous section. This is done without reducing the computational efficiency of the Numerical Kludge. We build our model with two main ingredients: a criterion to specify when a resonance starts, and a modification of the flux equation during resonance. The Numerical Kludge evolves firstly the phase-space trajectory of the constants of motion J = (E, L z , Q) and, since the fundamental frequencies are only functions of J, the condition for resonance can be checked at each time step. If it is satisfied, we modify the flux evolution to include the effects of the resonance.
The 3:2 resonances are considered to have the strongest impact on the inspiral [41], and are encountered by EMRI systems long enough before plunge that a significant phase shift can accumulate [48]. We will therefore explain how to construct the Effective Resonance Model for the case of a 3:2 transient orbital resonance, but the extension to other resonances is straightforward. All frequencies and times will refer to coordinate time, unless clearly stated otherwise.
The resonance condition is satisfied when the ratio of the polar and radial frequency is 3/2: ω θ0 /ω r0 = 1.5 at a given τ = τ 0 . However, the adiabatic approximation does not break down at the specific instant of time τ = τ 0 where ω θ0 /ω r0 = 3/2, but in the neighbourhood of the resonance. The question is where?
Due to the symmetry of the phase around τ 0 , we assume that the start of the resonance, τ start , is half the resonance duration before the resonance condition is satisfied at τ 0 , i.e., τ start = τ 0 − τ res /2. Since we do not know τ 0 beforehand, we use the following approach to trigger the start of the resonance in our model. If we define a threshold function as ξ = ω θ /ω r − 1.5, which is zero at τ 0 , it is possible to find the value of this function at τ start = τ 0 − τ res /2 such that the total duration of the resonance is τ res . To do this, we expand the threshold function around τ 0 up to the first order Thus, the threshold function at τ start = τ 0 − τ res /2 is ξ * = ξ(τ 0 − τ res /2) = −π/(τ res ω r0 ). Assuming that τ res and ω r0 are approximately constants in the neighborhood of the resonance, we can compute ξ = ω θ /ω r − 1.5 and ξ * = −π/(τ res ω r0 ) and will trigger the start of the resonance when ξ > ξ * . This is a universal choice which is independent of the specific EMRI system since the threshold is dimension-free. We verified that by using this criterion to turn on the resonance correction, the ratio 3/2 was reached exactly at τ 0 , as expected. So, if |ω θ0 /ω r0 − 1.5| < |ξ * | and if the orbit is noncircular and non-equatorial, we modify the flux evolution of the Numerical Kludge by implementing an Effective Resonance in coordinate time: dJ where we define the resonance coefficients C = (E, L z , Q) as the fractional change in flux on resonance, and the "impulse" function w as Evolution of the constants of motion through a 3:2 resonance to represent the evolution of the second term in the flux equation (11) due to the resonance. In Figure (1) it is shown how the effective model changes the evolution of J and implements a smooth deviation in the constants of motion during the resonance. The size of the changes are different for each of the constants of motion and they are of order ∼ 0.5% by the end of the resonance.

C. Resonance Coefficients
Determining the resonance coefficients C as a function of the EMRI parameters is possible, but is very challenging. The value of C for a specific trajectory can be estimated by using Eq. (11) following the procedure described in [47]. However, this procedure is currently too computationally expensive for our purposes. As far as we are aware, the only available resonance coefficient estimates are given in [47] and their largest values reach ∼ 0.01.
We can check that the Effective Resonance Model reproduces an evolution of the orbital constants and phases similar to that one given in [40]. There, Flanagan and The sign of C can be positive or negative according to the phase Ψ 0 at which the compact object reaches the resonance. Determining the correct phase Ψ 0 at resonance can be problematic since it would require the waveform model to be perfectly in phase from the beginning of the inspiral up to the start of the resonance. Therefore the resonance coefficients C are considered as parameters of our Effective Resonance Model. Once we have predictions from the GSF, the resonance coefficients could be used as a new way to test GR. For the time being we take as fiducial values C = (−0.01, −0.01, −0.01).

D. Impulse function
The "impulse" function w(t; τ res , τ start ) expresses the rising and the disappearing of the additional term during the resonance, and its functional form has been chosen such that: (i) it is symmetric with respect to τ 0 = τ start + τ res /2 ; (ii) it is smooth, positive and normalized to τ res ; (iii) it rises from zero with zero slope w (t = τ start ) = 0, and symmetrically it decays to zero with zero slope w (t = τ start + τ res ) = 0 ; (iv) it resembles the functional form of the correction in Eq. (11): exp ix 2 ∼ cos x 2 ; The first condition (i) is imposed to preserve the symmetry with respect to τ 0 as expressed in equation (11). The impulse function is normalized to τ res because the total change in the constants of motion is proportional to the resonance duration. Condition (iii) and (ii) avoid abrupt behaviors in the evolution of the constants. Conditions (i, ii, iii) are chosen according to the study of the semianalytic solution of resonances given in [69]. Condition (iv) is suggested by the functional form of the equations regulating the resonances. Since condition (iv) is guided by our intuition, we test that the final waveform is not affected by this choice. Therefore, we compute the mismatch M between oneyear EMRI waveforms with different impulse functions in the Effective Resonance Model (initial conditions given in [40]). In Table (I) we report the mismatch between the waveforms generated with our standard (S) impulse function w of eq.(14) and four different impulse functions: Blackman-Nuttall window (B-N), Blackman-Harris window (B-H), Tukey window (T) and Nuttall window(N) [70][71][72]. We also report the mismatches for different resonance coefficients C. As shown from Table I, the mismatches are particularly low for all the considered cases, ∼ 10 −4 − 10 −6 . Interestingly, the value of the mismatch is quadrupled if the resonance coefficients is doubled. An explanation for this behavior will be provided in Sec. IV B. Since the largest expected resonance coefficients C = (−0.01, −0.01, −0.01) lead to a mismatch of order 10 −4 , there is little difference between using different impulse functions.
We conclude that the functional form of the impulse function does not matter too much, as long as the phase immediately after resonance is correctly modelled. In  fact, we do not aim to exactly describe the evolution of the inspiral during the resonance, but to correctly account for their effects afterward, which is where large dephasings can accumulate. As long as we are back in phase after the resonance, and the resonance duration is relatively short, we can use the Effective Resonance Model to match astrophysical resonant EMRIs without significant loss of signal-to-noise ratio. We will use the impulse function of Eq. (14) for the rest of this work Our aim is to study how resonances influence the gravitational wave signals from EMRIs and how parameter estimation could be affected. Even though we are using a phenomenological model and we do not use the exact value of the resonance coefficients, this model is a reasonable starting point to understand at which order of magnitude resonances are going to affect the gravitational waves and EMRI parameter estimation. In addition, it provides a new tool for observing and measuring resonances. The Numerical Kludge combined with this phenomenological model has the advantage that it is fast enough to scan a large EMRI parameter space and to capture the crucial physics of resonances and EMRI evolution.

A. Dephasing caused by 3:2 resonances
The gravitational wave of an EMRI system consists of a superposition of multiple modes. Each mode is characterized by an amplitude and an oscillating sinusoid. The phase evolution of this oscillating term depends on a particular combination of the three Keplerian phases (ψ, χ, φ) associated with the EMRI. The overall gravitational wave phase of an EMRI system thus consists of a complex combination of the evolution of the phases (ψ, χ, φ). Resonances cause these phases to depart from the standard adiabatic evolution. In this section we will quantify the size of these dephasings in terms of how many cycles the resonant evolution deviates from from the non-resonant evolution, (∆ψ, ∆χ, ∆φ)/(2π), by the end of the one year observation. This analysis will serve as a verification of the Effective Resonance Model, and allow comparison with prior Dephasing over the parameter space work, but we will also try to quantify the dephasings in a different way. Previous investigations have considered a limited set of EMRI configurations and analysed the dephasings induced by resonances using a Post-Newtonian model for the resonant force [40,41]. The Post-Newtonian force model predicts a particular relationship between the parameters of the EMRI at resonance and the size of the flux changes on resonance, which may or may not be a good approximation to the true impact of the GSF. In the Effective Resonance Model, by contrast, the size of the flux changes is set by hand, and so we can compare the dephasing with these fixed and thus decouple the effect of the size of the flux change from the timing of the resonance within the inspiral. These results can thus be more easily reused once GSF calculations have accurately computed the on-resonance flux changes. In the following, we will use our fiducial resonance coefficients C = (−0.01, −0.01, −0.01) to estimate the dephasings. This can be regarded as a worst case scenario, but these dephasings can be readily rescaled to other choices.
By studying which regions of the EMRI parameter space lead to the largest dephasings, we can understand which EMRI configurations are mostly likely to be affected by 3:2 resonances. Thus, we evolve several EMRI systems with initial conditions e ∈ [0.1, 0.8], ι ∈ [0.1, 3.1], a ∈ [0.8, 0.95], with steps, ∆e = 0.1, ∆ι = 0.3 and ∆a = 0.05. The initial semi-latus rectum p/M is determined using the Black Hole Perturbation Toolkit [73] to start the inspiral near the resonance, more precisely at |ξ| = 3/2 − 1.498 (see Appendix V C for the parameter space location of the 3:2 resonance). We limit our analysis to configurations that have an initial periastron r p > 5 M for consistency with the limitations of the NK model. These parameter ranges are chosen such that they span the range of eccentricities e, inclination angle parameters ι and spins a expected for astrophysical EMRIs [30,[74][75][76]. We fix the mass ratio and the mass of the central black hole to η = 10 −5 , M = 10 6 M , respectively, and we set the initial phases to zero. For most of these configurations, the dephasing was computed over an observation time of one year. For systems with ι < 1.3 the location of the resonance in parameter space is such that the systems plunge in less than one year. In those cases the dephasing was calculated at plunge.
The dephasing for fixed inspiral length scales like ∆ω T 2 , where ∆ω is the change in the frequency derivative that accumulates over the resonance, and T is the time between the resonance and the end of the inspiral. The change in frequency derivative depends on the difference in the energy flux at resonance, and on the duration of the resonance. So, the largest phase shifts will be for systems that have large changes in the flux on resonance and long resonance durations. The longest resonance durations are for systems with high eccentricities and high ι, but the biggest phase shifts are for systems with smaller ι. These latter systems have resonances which occur closer to the central black hole, where the absolute magnitude of the dissipation rate is higher. For fixed T , this leads to a higher ∆J(τ res + τ start ) compared to other configurations and, as a consequence, higher phase shifts (∆ψ, ∆χ, ∆φ)/(2π), as we can see from Figure (2). The turnover in the dephasing for configurations with ι 1.0 is due to the fact that these systems plunge within one year. The total time over which the dephasing can accumulate decreases as ι decreases from ∼ 1.0 to zero, and this decrease compensates for the increased strength of the resonance effect to give a total dephasing that is decreasing.
Overall we conclude that, for fixed inspiral length, we expect the strongest impact on the waveform from the 3:2 resonances located in the strong field regime, i.e., closer to the central black hole.

B. Mismatch as a function of resonance strength
The dephasing due to a resonance accumulates over the inspiral, which means the overlap between a resonant and non-resonant waveform drops over time (or equivalently the mismatch increases). This has been shown to affect the detection of EMRI systems [41]. The size of the accumulated mismatch will depend also on the strength of the resonance. Here, we study the evolution of the mismatch, M, as a function of the resonance coefficients, C, for different EMRI initial conditions. By doing this we can assess which values of C lead to a significant dephasing.
We consider a one-year inspiral of an EMRI system with η = 10 −5 and M = 10 6 M . For simplicity we set the three resonance coefficients to be equal E = L z = Q, but have verified that the results are not significantly influenced by this choice. As we increase the size of the resonance coefficients, we calculate the mismatch M(h, h res ) between non-resonant and resonant waveforms for the two LISA channels. Since the results of the two channels are basically identical we show the results in Figure (3) for channel I only.
These results reflect some of the previous findings. The region of the parameter space where the dephasings are larger have a larger mismatch. In fact, prograde orbits (dash-dot-dotted green and dotted blue lines) show a larger mismatch than retrograde (dashed orange and, solid and dashdotted red lines). The spin parameter mildly affects the overall behavior. For higher eccentricity systems, the mismatch begins to increase significantly for smaller absolute values of the resonance coefficients. We find a direct proportionality between the mismatch log |M| and the resonance coefficients log C for C < 3 · 10 −4 . This is due to the fact that the mismatch Mismatch as a function of the resonance strength scales like the dephasing squared, and, as previously discussed, the dephasing at the end of the inspiral scales approximately like ∼ ∆ω ∼ C. Therefore the scaling log |M| ∝ log C is not unexpected.
A common criterion for assessing if two waveforms are indistinguishable is that the norm of the waveform difference, (δh|δh) < 1. This criterion was first introduced in [77] and [78], but was popularised by [79]. The mismatch is approximately (δh|δh)/(2ρ 2 ) and so the indistinguishability criterion is satisfied for M 1/(2ρ 2 ), which is approximately 10 −3 for the ρ ∼ 20 systems considered here. From Figure (3) we conclude that for resonance coefficients smaller than ∼ 10 −5 the waveforms are indistinguishable, regardless of the choice of system parameters, and so the effects of the resonance can be neglected. Furthermore, we infer that the mismatch strongly depends on the parameter space location of the EMRI system and that the largest mismatches, for a fixed size of the resonance coefficients, occur for high eccentricity, prograde orbits in the small resonance strength regime |C| < 10 −3 .
We have considered a worst case scenario for the mismatch because we did not maximize over the initial time, phases or other model parameters in order to find the best-fit template, as is normally done in GW data analysis. However, we do not expect the above results to drastically change with this maximization, as already demonstrated in [41]. For δλ bias /∆λ > 1 the bias induced by inaccurate modeling is larger than the bias induced by the typical noise fluctuations. This figure shows how the bias is distributed over many resonance strength realizations (resonance coefficient realizations). For the fiducial EMRI configuration (solid blue line), more than 95% of the biases are larger than two sigma, i.e. the cumulative distribution reaches ≈ 0.05 at δλ bias /∆λ = 2.

C. Systematic errors
The use of an approximate waveform model causes biases in parameter estimation. Therefore, the scientific potential of EMRI systems can be undermined if the systematic error δλ bias induced by inaccurate waveform modeling is larger than the size of the typical statistical error (Eq. 9) [68]. Statistical errors arise due to noise fluctuations and are therefore unavoidable, but are accounted for in the width of posterior distributions on parameters. Requiring the systematic error to be smaller than the statistical error ensures that the true parameters remain consistent with the posterior for single events, although errors can accumulate for populations [80]. The statistical error can be estimated as the square root of the diagonal elements of the inverse of the Fisher matrix ∆λ and so we are interested in the ratio where there is no sum over the index i, h m and h t refer to the approximate and true waveforms respectively, and each i refers to one of the considered parameters λ = (p/M, log η, log M, e, ι). The subscript I, II on the Fisher matrices mean that they have been calculated for the respective detector channels. Here we investigate the size of the errors that arise from ignoring the presence of resonances in the waveform model. The approximate waveform, h m , will be the NK waveform without resonance and the true waveform, h t , will be the waveform produced with the Effective Resonance Model.
We consider a fiducial EMRI configuration η = 10 −5 , M = 10 6 M , e = 0.1, ι = 1.3, p/M = 9.0447, a = 0.9. The resonance coefficients, C, are in principle a function of the waveform parameters and can be computed using the GSF. However, such calculations are beyond the scope of this paper and so we draw random values for the resonance coefficients from a uniform distribution (E ∼ L z ∼ Q ∼ U [−0.01, 0.01], where each resonance coefficient is independently sampled). Then, we look at the distribution of biases over repeated random draws of this kind. The domain of the uniform distribution is chosen to cover all the expected resonance coefficient values. To explore the impact of the other waveform parameters we repeat this procedure for other EMRI systems which differ from the fiducial system by the change of one parameter at a time. For each system, the semilatus rectum was adjusted to ensure the 3:2 resonance was encountered close to the start of the inspiral, i.e., at the start of the inspiral we had |ξ| = 3/2 − 1.498. In all cases the EMRI did not plunge within the one year observation time used.
We show the cumulative distribution of the biases for the various EMRI systems in Figure (4). The cumulative distribution of δλ bias /∆λ for the fiducial configuration (blue solid lines) shows that more than 95% of the biases are more than twice the statistical uncertainty, ∆λ, ("two sigma away"). Reducing the spin (red dashdotted lines) or the mass ratio η (violet dashed lines) does not significantly change the cumulative distribution. For these configurations we conclude that EMRI waveform models need to account for resonances. For the EMRI systems with higher eccentricity (orange dotted lines) and larger ι (green dash-dot-dotted lines), the proportion of systems with biases above 2 sigma is lower, but these proportions still exceed ∼ 85% and ∼ 25% respectively. These results are consistent with the mismatch analysis presented earlier, where we found typically lower mismatches for more retrograde systems. Overall, the conclusion from Figure (4), is that it is important to use waveform models that include resonances when performing EMRI parameter estimation. The biases will be less severe for retrograde and higher eccentricity systems, although even there the biases are not negligible and will continue to accumulate if the observation time is increased beyond one year.
The bias estimates here were calculated using the approximate formalism described in [68] and its validity will degrade for large biases, δλ bias /∆λ 10. This does not invalidate our findings but it is likely to affect the reliability of our estimate of the tails of the cumulative distribution of biases. The tails could be more accurately resolved by computing posterior distributions for each configuration. This is computationally expensive and so is beyond the scope of the current work, but we do not expect the conclusions to be significantly different in such an analysis.

D. Measurability of Resonances
We now investigate two things: the measurability of resonances and how the inclusion of the resonance coefficients C as unknown model parameters change the measurement precision of other intrinsic parameters. This is assessed using the Fisher information matrix. If ∆λ/λ < 1 then the parameter is measurable. For details about the validation and checks of our Fisher matrix calculations we refer the reader to Appendix V B.
As before, we consider a fiducial EMRI configuration and we calculate the relative measurement precision ∆λ/λ for uniformly randomly drawn resonance coefficients. We repeat this for the same set of reference systems considered in the previous section, which vary one parameter at a time away from the fiducial model. All the waveforms have been normalized to SNR = 20 and we compute measurement precisions for the parameters (p/M, η, M, e, ι, a, E, L z , Q). We show the results for the distribution of parameter measurement precision in rameters (p/M, η, M, e, ι, a) are at the sub-percent level in all cases, and vary by a factor of a few across the distribution of resonance coefficient values. This shows that the inclusion of the resonance coefficients as additional parameters does not significantly degrade the measurement precision for the intrinsic parameters. In addition, we find that these precision estimates are in agreement with the typical values found in [30]. It is clear from Figure (5) that the resonance coefficients cannot be constrained for the high ι configuration (green dashdotted). In fact as previously noticed, the parameter space location of the 3:2 resonance for high ι is particularly far (p/M ∼ 12) from the strong field regime and this leads to smaller dephasings with respect to the other configurations. Since the time to plunge for this configuration (green dashdotted) is three years, we expect that if we had evolved such an EMRI system for a longer time, the measurement precision of the resonance coefficients for this system would improve.
On the other hand, the most precise measurement estimates are obtained for the smaller mass ratio system (violet solid). This is probably due to the choice of increasing the observation time for this system. While that was done to ensure the amount of inspiral was comparable, it means that twice as many waveform cycles are also observed so we expect to measure parameters somewhat better. We find that for this EMRI configuration the resonance coefficients E, L z , Q are determined with relative precision better than 0.1, respectively, 54%, 35%, 81% of the time. Here, we report the median relative precision of the three resonance coefficients, (

V. DISCUSSION AND FUTURE OUTLOOK
Extreme Mass Ratio Inspirals are modeled using perturbation theory in the small mass ratio η. The effects of transient orbital resonances on the EMRI phase evolution scale as ∼ η −1/2 , therefore contributing more than postadiabatic corrections which scale as ∼ η 0 . It has been shown how resonances affect detection [41]. However, the impact of these phenomena on parameter estimation had not previously been investigated.
In this work we have explored this question by implementing a phenomenological model for transient orbital resonances: the Effective Resonance Model. This model allows us to efficiently study a wide variety of EMRI systems avoiding expensive gravitational self-force calculations, but at the cost of introducing three additional model parameters: the resonance coefficients C = (E, L z , Q). These encode the relative flux changes for each of the constants of motion that occur dur-ing resonance, and, therefore represent the resonance strength. Since the largest estimated resonance coefficients in the literature are of order ∼ 0.01, we adopt the fiducial values C = (−0.01, −0.01, −0.01) to study the dephasings induced by 3:2 resonances over the parameter space. For a one-year EMRI inspiral with parameters η = 10 −5 , M = 10 6 M , we find that the maximum dephasings amount to (14.8, 12.2, 2.0) cycles for (ψ, χ, φ) respectively, and occur for prograde orbits. This is due to the fact that the 3:2 resonances of these configurations are located closer to the central MBH where the absolute magnitude of the fluxes is higher.
We have also investigated how the mismatch between a resonant and non-resonant waveform depends on the resonance coefficients, finding a proportionality between the logarithm of the mismatch and the logarithm of the resonance coefficients. The study of the mismatch was consistent with the behavior of the dephasings over the parameter space. In addition, it revealed that more eccentric systems might be more affected by resonances in the small resonance strength regime |C| < 10 −3 . We also inferred that resonance coefficients with absolute values less than ∼ 10 −5 lead to waveforms that are indistinguishable from the waveforms without resonances and hence there is no measurable effect in the phase evolution.
We conducted a study of the systematic errors that would arise in estimates of the intrinsic parameters from neglecting resonances. We found that for the considered EMRI configurations the biases can reach values up to a dozen times larger than the statistical errors arising from noise fluctuations. Over all the resonant coefficient realizations, we find that more than 95% of the biases are larger than 2 sigma, for our fiducial EMRI configuration. This suggests that parameter estimates of resonant EM-RIs are likely to be biased if resonances are not taken into account in models used for parameter estimation. We also expect an increase in the parameter bias for longer signals.
Finally, we used the Effective Resonance Model to assess the measurability of parameters with resonant EM-RIs. The EMRI systems which are mostly affected by resonances provide measurements of the resonance coefficients with a relative precision ranging from 0.07 − 0.68 (median values).
The Effective Resonance Model presented here has a number of potential uses, beyond being used to scope out the impact of resonances as done in the current paper. A model similar to this could be used in analysis of LISA data to match EMRI signals over resonances, mitigating the induced biases without requiring expensive GSF calculations of the resonance coefficients. Furthermore, due to its general implementation, the Effective Resonance Model can also be adapted to represent tidal resonances [44], which have a similar form to transient orbital resonances but are caused by the tidal perturbation of a third object. Finally, such a model can be used to measure resonance coefficients in EMRI observations. Comparing the measured values to the self-force prediction of the resonance coefficients would provide a test of general relativity.
There are a number of ways in which this work could be extended. We have explored only 3:2 resonances here, so it would be interesting also to analyse higher order resonances, or systems that pass through multiple strong resonances. In addition, our parameter estimation and bias results were based on the linear signal approximation, and it would be useful to verify these using full Bayesian posterior calculations. Finally, once full GSF waveforms with resonances are available, the ability of the Effective Resonance Model to detect and characterise such systems should be assessed.
While the development of an accurate self-force waveform model will be necessary to have more precise and quantitative results, the Effective Resonance Model has provided a first step towards an understanding of the impact of transient orbital resonances on parameter estimation.

ACKNOWLEDGMENTS
We thank Matthias Bartelmann, Ollie Burke, Maarten van de Meent, Michael Katz, Vojtech Witzany and Béatrice Bonga for their useful suggestions and discussions. This work makes use of the Black Hole Perturbation Toolkit [73].

A. Data Analysis Caveats
The Numerical Kludge produces approximate numerical waveforms expressed as an array of data points spaced by a constant time sampling interval, ∆t. Since we are not dealing with continuous and analytical signals, we are going to briefly discuss different caveats that can occur when calculating Fourier transforms and derivatives.
The choice of a too large sampling interval can lead to aliasing and it imposes a limit on the maximum resolvable frequency f max = 1/(2∆t) as prescribed by the Nyquist-Shannon sampling theorem. On the other hand, a too small sampling interval requires a higher computational cost. Throughout this work we fixed ∆t = 10s, which is sufficient to resolve the frequencies accessible to the LISA sensitivity band without unduly increasing the computational time.
The waveform derivatives ∂ j h are calculated numerically using the five-point stencil formula to ensure a high stability of the derivatives and a numerical error that scales at fourth order in the derivative spacing. In order to avoid spectral leakage and to sample over the minimum possible frequency range we apply a Tukey window function with shape parameter 0.05 to every waveform and then zero-pad them to exploit the efficiency of the Fast Fourier Transform [81].

B. Fisher matrix validation
We make use of the following result, valid for high SNR in the linear signal approximation regime, to check the correctness of the Fisher matrix results and numerical derivative approximations: The choice of the perturbation γ must be small enough to respect the linear signal approximations but also big enough to test the 1σ region. We follow a similar approach to those described in [82,83], taking γ to be a linear combination of the eigenvectors, v i A , and corresponding eigenvalues, w A , of the Fisher matrix, Γ, where the index A runs over the N different eigenvectors. We draw coefficients c A from an N dimensional unit sphere and we define γ as follows For every computed Fisher matrix, we calculate separately the left hand side and right hand side, neglecting terms of order O(ρ −4 ) and higher in Eq. (16). By taking the ratio between the two sides, we find a maximum deviation of order 10 −5 . This is a way of checking both the validity of the Fisher matrix and the linear signal approximation. The Fisher matrices encountered in EMRI data analysis are known for having a high condition number, κ = max(w A )/min(w A ), which is the ratio between the maximum and minimum eigenvalue. A small perturbation in a Fisher matrix with high condition number κ can be amplified by a factor of κ in the inversion. Therefore, if we have inaccuracies in our Fisher matrix, Γ, for instance due to our numerical derivative approximation, they can lead to a wrong value of ∆λ i = (Γ −1 ) ii , invalidating our measurability conclusions. We perturb each element of every Fisher matrix with a deviation matrix F ij , where each element is drawn from a uniform distribution U [−10 −3 , 10 −3 ], and we calculate max all configurations This confirmed the stability of our Fisher matrix results.

C. Parameter space location of 3:2 resonances
In Figure (6) we show the parameter space location of 3:2 resonances, calculated using the Black Hole Perturbation Toolkit [73]. The values of p/M, e, ι, a just before resonance can be calculated by searching the roots of ω θ /ω r = 1.498. These values are used as initial conditions to evolve EMRI system through resonances.