Edinburgh Research Explorer Constraining the spin parameter of near-extremal black holes using LISA

We describe a model that generates ﬁrst order adiabatic EMRI waveforms for quasi-circular equatorial inspirals of compact objects into rapidly rotating (near-extremal) black holes. Using our model, we show that LISA could measure the spin parameter of near-extremal black holes (for a (cid:38) 0 . 9999 ) with extraordinary precision, ∼ 3-4 orders of magnitude better than for moderate spins, a ∼ 0 . 9 . Such spin measurements would be one of the tightest measurements of an astrophysical parameter within a gravitational wave context. Our results are primarily based oﬀ a Fisher matrix analysis, but are veriﬁed using both frequentist and Bayesian techniques. We present analytical arguments that explain these high spin precision measurements. The high precision arises from the spin dependence of the radial inspiral evolution, which is dominated by geodesic properties of the secondary orbit, rather than radiation reaction. High precision measurements are only possible if we observe the exponential damping of the signal that is characteristic of the near-horizon regime of near-extremal inspirals. Our results demonstrate that, if such black holes exist, LISA would be able to successfully identify rapidly rotating black holes up to a = 1 − 10 − 9 , far past the Thorne limit of a = 0 . 998 .


I. INTRODUCTION
Extreme mass ratio inspirals (EMRIs) are one of the most exciting possible sources of gravitational radiation for the space-based detector LISA [1], but also one of the most challenging to model and extract from the data stream.An EMRI involves the slow inspiral of a stellarorigin compact object (CO) of mass µ ∼ 10M into a massive black hole in the centre of a galaxy.For a central black hole with mass M ∼ 10 (5−7) M , EMRIs emit gravitational waves (GWs) in the mHz frequency band and so are prime sources for the LISA detector.EM-RIs begin when, as a result of scattering processes in the stellar cluster surrounding the massive black hole, the CO becomes gravitationally bound to the primary.The subsequent inspiral of the CO towards the horizon of the primary is driven by radiation reaction through the emission of gravitational waves.EMRI waveforms are very complicated and EMRIs can be present in the LISA frequency band for several years prior to plunge, so modelling the full observable signal is a complex task [2].EMRI orbits are expected to be both eccentric and inclined even up to the last few cycles before plunging into the primary black hole.For these reasons, EMRIs pose a challenging problem for both waveform modellers [2][3][4] and data analysts.
This same complexity also makes EMRIs one of the richest sources of gravitational waves.Typically an EMRI will be observable for 1/(mass ratio) ∼ 10 5−7 cycles before plunge and the emitted gravitational waves thus provide a very precise map of the spacetime geometry of the primary hole [5][6][7][8].Through accurate detection and parameter inference, one can conduct tests of general relativity to very high precision [6,9].
It is well known that the information about the source is carried through the time evolution of the phase in a gravitational wave [10,11].The slow evolution of EM-RIs means that a large number of cycles can be observed during the inspiral, which will provide constraints on the parameters of the source with remarkable precision [12,13].Previous work has indicated that LISA will be able to place constraints on the dimensionless Kerr spin parameter, a, of the primary black hole in an EMRI, at the level of 1 part in 10 4 for moderately spinning, a ∼ 0.9, primaries [3,4,14,15].In this paper, we explore how well LISA will be able to measure the spin parameter for very rapidly rotating black holes, i.e., systems in which the spin parameter is close to the maximum value allowed by general relativity.
Super massive black holes with a large spin parameter are abundant throughout our universe.Observations indicate that massive BHs reside in the centres of most galaxies, where these black holes are known to accrete matter and hence are predicted to have very high spins [11,[16][17][18][19][20].The dimensionless Kerr spin parameter of a Kerr black hole, cannot exceed 1, since the resulting spacetime contains a naked singularity no longer encased within a well defined horizon.Thorne [21] showed that a moderately spinning black hole cannot be spun up by thin-disc accretion above a spin of a ≈ 0.998.How-ever, in principle primordial black holes could be formed with spins exceeding that value [22]."Near-extremal" black holes with spins close to the limit of a = 1 have interesting properties and we focus our attention on these here.
This past decade, researchers [23][24][25][26][27][28][29][30][31][32] have explored the rich properties of near-extremal EMRIs.The gravitational radiation emitted from these systems is unique, and would prove a smoking gun for the existence of these near-extremal systems (see [24]).In this paper, we show qualitatively that the inspiraling dynamics of the compact object into an near-extremal massive black hole is very different from that into a moderately spinning black hole, and these differences are reflected in the emitted gravitational waves.As such, in order to detect and correctly perform parameter estimation on these near-extremal sources, it is essential to update our family of waveform models to include them.We will argue throughout this work that, if observed, nearextreme black holes offer significantly greater precision measurements on the Kerr spin parameter than moderately spinning systems.In particular, LISA will have the capability to successfully conclude whether the central object in an EMRI system is truly a near-extremal black hole.Thus, if near-extremal black holes exist, LISA observations of EMRIs may be one of the best ways to find them.
In this paper we will consider only EMRIs on circular and equatorial orbits around near-extremal primary black holes.This choice is made primarily for computational convenience, but there are also astrophysical scenarios that produce such systems.As discussed in [33], compact objects can form within accretion disks around massive black holes.When these objects fall into the central black hole, the resultant EMRI will be circular and equatorial.Super-Eddington accretion can provide a means to spin up a black hole past the Thorne limit [34], and so it is not unreasonable to expect that this EMRI formation channel would be more important for near-extremal systems.The standard EMRI formation channel, involving capture of a compact object via scattering interactions, tends to form EMRIs with moderate initial eccentricities.However, this eccentricity decreases during the inspiral due to the emission of gravitational radiation [35].This decrease in eccentricity continues until the orbit reaches a critical radius at which is starts to increase again [36,37].The critical radius moves closer to the last stable orbit as the spin parameter increases and for near-extremal systems is located within the regime where transition from inspiral to plunge occurs [38,39].Additionally, the increase in eccentricity is a subdominant effect throughout the transition regime [40].As the spin increases, we therefore expect that for an object captured at a fixed radius, the amount of eccentricity dissipated before the critical radius increases, and the eccentricity gained af-ter the critical radius decreases.Therefore, even in the standard capture picture it is reasonable to assume the eccentricity is small at the end of the inspiral.We will show in this paper that very precise measurements of spin for near-extremal systems are possible, but this precision comes from observation of features [24] in the final phase of the inspiral, which is where the near-circular assumption is most likely to be valid.
The main results of the paper are given in figures 9 and 10 in section VI.Readers who wish to understand why near-extremal systems offer greater precision spin measurements than moderate spin systems should direct their attention to section III.This paper is organised as follows.In section II, we set notation and discuss the trajectory of a compact object on a circular and equatorial orbit around a near-extremal Kerr BH.In section III, we show that the spin dependence of kinematical quantities appearing in the radial evolution rather than radiation-reactive effects dominate the spin precision measurements for near-extremal EMRI systems.Our Teukolsky based waveform generation schemes are outlined in section IV.We discuss prospects for detection in section V, arguing that LISA is more sensitive to heavier mass systems M ∼ 10 7 M than lighter systems M ∼ 10 6 M .Our Fisher Matrix results are presented in section VI.
Here we show that we can constrain the spin parameter ∆a ∼ 10 −10 , even when correlations amongst other parameters are taken into account.Finally, in section VII, we perform a Bayesian analysis to verify our Fisher matrix results, before finishing with conclusions and outlooks in section VIII.

II. BACKGROUND
We consider the inspiral of a secondary test particle of mass µ on a circular, equatorial orbit around a primary super massive Kerr black hole with mass M and Kerr spin parameter a 1 where the mass ratio is assumed small η = µ/M 1.The secondary is on a prograde orbit aligned with the rotation of the primary black hole with a > 0 and dimensionful angular momentum L > 0. Unless stated otherwise, throughout this paper any quantity with an over-tilde is dimensionless, e.g., r = r/M and t = t/M etc.The one exception is the dimensionless spin parameter, which we denote by a without a tilde.Quantities with an over-dot will denote coordinate time derivatives, e.g., ṙ = dr/dt.We use geometrised units such that G = c = 1.
In Boyer-Lindquist [41] coordinates, the metric of a Kerr black hole for θ = π/2 is given by where ∆ = r2 − 2r + a 2 and a is the dimensionless spin parameter introduced earlier.This is related to the mass, M , and angular momentum, J, of the Kerr black hole via a = J/M and lies in the range a ∈ [0, 1].
The event horizon is located on the surface defined by ∆ = 0, when Introducing an extremality parameter 1 the event horizon is at The trajectory of the secondary confined to the equatorial plane of a central Kerr hole is governed by the Kerr geodesic equations [42] r2 dr dτ in which τ denotes the proper-time coordinate for the inspiraling object.The dimensionless conserved quantities Ẽ = E/µ and L = L/(M µ) are related to the energy, E, and angular momentum, L, measured at infinity.For the circular and equatorial orbits considered here, the energy Ẽ and angular frequency Ω can be expressed analytically in which the dimensionless angular frequency Ω is defined through Ω = Ω/M = dφ/dt.Circular orbits exist only outside the innermost stable circular orbit (ISCO).For radii smaller than the ISCO, the secondary will start to plunge towards the horizon of the primary.The ISCO for equatorial orbits is at [43] For near extremal orbits, using (3) and ( 7), and expanding for 1, we obtain and deduce The radial coordinate separation between the ISCO and horizon is determined by the spin parameter.In the limit, → 0, then risco → r+ → 1 in Boyer-Lindquist coordinates.

A. Radiation Reaction
To compute circular and equatorial adiabatic inspirals, a detailed knowledge of the radial self force is required (see, for example, [44] for a detailed review).In this paper, we will work at leading order, including the radiative (dissipative) part of the radial self force at first order, but neglecting first order conservative effects and all second order in mass-ratio effects.The first order dissipative force can be computed by solving the Teukolsky equation [45].The rate of emission of energy is given by where ĖGW = −(u t ) −1 F t , • denotes coordinate time averaging over several periods of the orbit.The quantity u t is the t component of the four velocity and F t the t component of the gravitational self force at first order in the mass ratio η.This expression is valid only if η 1, that is, when orbits evolve adiabatically such that the timescale on which the orbital parameters evolve is much longer than the orbital period.
The fluxes Ė∞ and ĖH denote the (dimensionless and orbit averaged) dissipative fluxes of gravitational radiation emitted towards infinity and towards the horizon respectively.From here on, we shall drop the angular brackets Ė → Ė, to avoid cumbersome notation.The quantities |m| ≤ l are angular multiples which appear in the decomposition of the emitted radiation into a sum of spheroidal harmonics.The components of the fluxes Ė are obtained by numerical solution of the Teukolsky equation sourced by a point particle (the secondary).
There exists an open source code in the Black Hole Perturbation Toolkit (BHPT) [46] to do this for circular and equatorial orbits -specifically the Teukolsky package.
The enhancement of symmetry in the near horizon geometry of extreme Kerr [47] provides an additional tool to compute the fluxes Ė analytically from first principles.See [23][24][25][26][27][28][29][30][31][32] for a description of this work.For circular equatorial orbits near the horizon of a nearextremal black hole, there is a remarkably simple approximation for the total flux [25], which takes the form The quantities CH and C∞ are constants representing the emission towards the horizon and infinity respectively.These constants are given analytically in equations ( 76) and (77) of [25] and codes in the BHPT can be used to evaluate them.Numerically evaluating them and summing the contribution of the first |m| ≤ l = 30 modes gives CH ≈ 0.987 and C∞ ≈ −0.133.Eq.( 11) is useful when working within the near-horizon geometry of the rapidly rotating hole, but it breaks down far from the horizon and extra terms would be required to compute reliable fluxes.
All the numerical work presented in this paper, which is found in section V onwards, will use the exact fluxes obtained from BHPT.However, to understand our numerical results, we develop a set of new analytic tools in sections III A and III C.These will partially make use of the leading contribution to (11) This differs from (11) by O( ) contributions since r − 1 measures the BL radial distance to the extremal horizon and not the radial distance to the near-extremal horizon r+ .The approximation (12) can be derived from first principles by solving the Teukolsky equation in the NHEK region 1 .Our numerical analysis based on the BHPT, suggests the spin dependence of certain observables, to be discussed in section III B, is better captured by (12).Table I compares the flux at risco computed using BHPT to that obtained from the near-extremal approximations of Eq.( 11) and Eq.( 12).This table corroborates that ( 12) is a good approximation to the total energy flux, particularly in the limit as a → 1, where it outperforms the full expression, (11).
The radial evolution of the secondary can be found by taking a coordinate time derivative of the circular energy relation ( 5) where we defined P GW := ĖGW .As the ISCO is approached, the denominator ∂ r Ẽ tends to zero, marking a break down of the quasi-circular approximation.The ODE ( 13) is easily numerically integrated given an expression for the flux P GW .
The outgoing gravitational wave energy flux measured at infinity has a harmonic decomposition [11] where and Here Ė∞ m is the relativistic correction to the Newtonian expression for the flux in harmonic m.
In this work, we shall consider two different waveform models.For the analytic discussion in section III, we will use the waveform model in [11], whereas for the numerical analysis in later sections, we will use the full Teukolsky based waveform.
Let us first review the main features of the model discussed in [11] for the waveform observed by the detector in the source frame.This model is written Some remarks are in order.First, we ignore the m = 1 contribution since, as argued in [11], this is subleading to the m ≥ 2 contributions.Second, the amplitude ×m corresponds to the root mean square (RMS) amplitude of gravitational waves emitted towards infinity in harmonic m.These are averaged • over the viewing angle 2 and over the period of the waves.Third, the oscillatory phase depends on the initial phase φ 0 and the frequency fm of each waveform harmonic is given by Table I: NHEK fluxes at the ISCO computed using the approximations Eq. ( 11) (denoted Ė+ NHEK ) and Eq. ( 12) (denoted ĖNHEK ), and computed exactly using BHPT (denoted ĖExact and based on the first thirty m and l modes).
The relation between the RMS amplitude and the outgoing radiation flux in harmonic m is where D = D/M is the distance to the source from earth.Using Eq.( 14), we can rewrite h o,m as for m ≥ 2. We note that the effect of the averaging is that this waveform model does not represent the waveform measured by any physical observer.However, it captures the main physical features of the waveform which encode information about the source parameters.
Given the nature of our orbits, our parameter space will only be six dimensional θ = {r 0 , a, µ, M, φ 0 , D}, where r0 stands for the initial size of the circular orbit.We stress this waveform model does not include the LISA response functions, which affect the amplitude evolution of the signal and induce modulations, due to Doppler shifting, through the motion of the LISA spacecraft [48,49].Since these response functions do not depend on the intrinsic parameters of the system that we are most interested in, we omit these here and, consequently, they will also be omitted in our analytic discussion based on this waveform model.
Let us now review the full Teukolsky based waveform model that we will use in our numerical study.This is given by where depends on the radial Teukolsky amplitude at infinity, Z ∞ ml (r, a), and the viewing angle (θ, φ).The latter dependence is through the spin-weight minus 2 spheroidal harmonics −2 S am Ω ml (θ, φ) = −2 S am Ω ml (θ) exp(iφ).This work will consider two viewing angles: face on (θ, φ) = (0, 0) and edge on (θ, φ) = (π/2, 0).Using the identities where barred quantities are complex conjugates, we can write equation (20) as for the edge-on case, and as for the face-on case.Note we have neglected higher order l modes with m = 2 fixed in the last equation since the Teukolsky amplitudes Z ∞ l2 for l > 2 are negligible in comparison to the dominant quadrupolar l = m = 2 mode.Figure 1. in [25] further justifies our claim that higher order m modes when l = 2 can be ignored for face-on sources.Furthermore, the only spheroidal harmonics that are non-vanishing at θ = 0 are those with m = −s, or m = 2 [50,51].
To perform our numerics, the spheroidal harmonics are calculated using the SpinWeightedSpheroidalHarmonics mathematica package in the BHPT, whereas the Teukolsky amplitudes Z ∞ ml are calculated using the Teukolsky package from the same toolkit.For reasons discussed later, we generate both amplitudes and spheroidal harmonics for a fixed spin parameter a = 1 − 10 −9 .For the remainder of this study, we will only consider the plus polarised signal h(t; θ) ≡ h + (t; θ) for the face-on and edge-on observations.We finish this waveform discussion with a comment regarding the relation between the two models considered in this work.The (dimensionful) Teukolsky amplitudes are related to the energy flux for each (l, m) mode by Hence | Z∞ ml | ∼ M Ω η Ėlm .Averaging over the sky and ignoring the phase of the radial amplitude Z ∞ ml , the Teukolsky waveform (20) reduces to (16).Our numerical results indicate that the spin precision measurements are driven by the radial trajectory given by (13), which is common to both ( 16) and (20), while not being largely influenced by the spin dependence on the waveform amplitude.Given this fact and since it is analytically much easier to analyse the waveform model (16), this is the one being discussed in the analytics section III to explain the increase in the spin precision measurement for near-extremal primaries.

C. Gravitational Wave Data Analysis
The data stream of a gravitational wave detector, d(t) = h(t; θ) + n(t), is typically assumed to consist of probabilistic noise n(t) and (one or more) deterministic signals, h(t; θ), with parameters θ.Assuming that the noise is a weakly stationary Gaussian random process with zero mean, the likelihood is [52] with inner product Here b(f ) is the continuous time fourier transform (CTFT) of the signal b(t) and S n (f ) the power spectral density (PSD) of the noise.Here we use the analytical PSD given by Eq.( 1) in [53].We do not include the galactic foreground noise in the PSD to ensure all noise realisations generated through S n (f ) are stationary.This is not a serious restriction as for the sources we consider here, the majority of the GW emission is at higher frequencies where the galactic foreground lies below the level of instrumental noise in the detector.The optimal signal to noise ratio (SNR) of a source is given by This is the SNR that would be realised in a matched filtering search and is a measure of the brightness, or ease of detectability, of a gravitational wave signal.Measures of the similarity of two template waveforms h 1 := h(t; θ 1 ) and h 2 := h(t; θ 2 ) are the overlap If O(h 1 , h 2 ) = 1 then the shape of the two waveforms matches perfectly.Waveforms with O(h 1 , h 2 ) = 0 are orthogonal, being as much in phase as out of phase over the observation.Consider θ = θ 0 +∆θ for ∆θ a small deviation around the true parameters θ 0 .Assuming that the waveform h(t; θ) has a valid first order expansion 3 in ∆θ, we substitute into (27) and expand up to second order in ∆θ where ∆θ i bf = (Γ −1 ) ij (∂ j h|n) and Γ ij is the Fisher Matrix given by The Fisher Matrix Γ ∼ ρ 2 and therefore ∆θ scales like The linear signal approximation is therefore valid for high SNR, ρ 1.The Fisher Matrix Γ, evaluated at the true parameters θ 0 , provides an estimate of the width of the likelihood function (27).Hence, it can be used as a guide to how precisely you can measure parameters.The inverse of the Fisher matrix is an approximation to the variance-covariance matrix Σ on parameter precisions The square route of the diagonal elements of the inverse fisher matrix provide estimates on the precision of parameter measurements, accounting for correlations between the parameters.

III. ANALYTIC ESTIMATES OF SPIN PRECISION
Before discussing numerical results on the measurement precisions for the parameters θ of near-extremal 3 In the literature, this is called the linear signal approximation.
EMRIs, we would like to develop some analytic tools that will allow us to understand the precisions we find numerically.In particular the fact that spin measurements for near-extremal primaries are noticeably tighter than those obtained for more moderately rotating primaries.Throughout this section, we will use the waveform model ( 16) for analytical convenience.We will focus on the spin-spin component of the Fisher matrix in the following analytic discussion.Our numerical and statistical analysis will be more general and employ the Teukolsky based model (20).In future work, we will extend this analytic considerations to multiple parameter study.
If all other parameters were known perfectly, the estimated precision on the spin parameter would be Thus, to compare precisions between near-extremal (denoted ext) and moderately rotating (denoted mod ) primaries one is led to study the ratio Consider the (semi-analytic) gravitational wave amplitude ( 16) where we have chosen the initial phase φ 0 = 0 for simplicity.The Fisher Matrix depends on the PSD of the detector.In the numerical calculations presented later we will use the full frequency dependent PSD, but to derive our analytic results we will approximate The rationale for this is that EMRIs evolve quite slowly and so the total change in the PSD over the range of frequencies present in the signal is small.Between 1 mHz and 100 mHz, the (square root of the) LISA PSD changes by just one order of magnitude, which is much smaller than the three orders of magnitude improvement in spin measurement precision that we find numerically.Additionally, the difference in the ISCO frequencies across all combinations of mass and spin considered in our numerical analysis is less than a factor of 2.5.PSD variations can not therefore explain the numerical results, and so we can ignore these in deriving the analytic results which do explain the numerics.Under this approximation We additionally assume that the choice of f • does not depend on the spin, and therefore the ratio ( 37) is independent of S n (f • ).Again, this approximation could introduce at most an order of magnitude uncertainty, and most likely much less than that.Once the Fisher matrix is written in the form (39), we can use the semi-analytic waveform model (38) to evaluate it.In appendix A, we argue the dominant contribution can be approximated by Here t0 is the coordinate time at which the observation starts and tcut is the coordinate time at the end of the observation.For the results in this paper, we analyse ∼ 1 year long signals and fix tcut independently of spin, such that all inspirals terminate before risco is reached.
As seen in ( 40), a proper understanding of the precision in the spin measurement requires quantifying the spin dependence of the inspiral trajectory of the secondary, i.e. ∂ a r.

A. Spin dependence on the radial evolution
Our primary goal here is to understand the spin dependence on the radial trajectory of the secondary (∂ a r) for any spin parameter a of the primary.
The trajectory of the secondary is the integral of the inspiral equation This follows from energy conservation, where Ẽ(r, a) is the energy of a circular orbit (5) and is the energy rate carried away by gravitational waves (10).While Ẽ(r, a) is kinematic, that is, derived through geodesic properties, P GW is dynamic, that is, it is a radiation reactive term determined by solving Teukolsky's equation for a point particle source.The former is under analytic control, whereas the latter typically requires numerical treatment.
The quantity ∂ a r captures the change in the secondary's trajectory when the spin parameter a of the primary varies, keeping the remaining primary and secondary parameters fixed, including t.More explicitly, the integral r(r 0 , a) of ( 41) depends on the initial condition r( t0 ) = r0 and it depends on the spin parameter a both through (∂ r Ẽ) and (P GW ) information, but not through t, which is simply labelling the points in the trajectory.We will comment on the possible spin dependence on the initial condition r0 below.
One possibility to compute ∂ a r is to integrate (41) and to take the spin derivative explicitly afterwards.A second, equivalent, way is to observe r is a monotonic function of t at fixed spin and initial radius r0 .Hence, it can be used as the integration coordinate to study ∂ a r(r).To do this, notice that the total spin derivative of the kinematic and dynamic functions in (41), at fixed r0 and t, is To ease our notation, all spin partial derivatives in the rhs, and in the forthcoming discussion, should be understood as computed at fixed r0 and t.Defining u = ∂ a r (to ease notation) and computing the total spin derivative of equation ( 41), we obtain Plugging in the radial velocity using (41) one obtains This is a first order linear ODE, valid for any spin and for any location of the secondary, whose solution describes the desired spin dependence in the radial trajectory ∂ a r(r).
Its general solution is a sum of the homogeneous solution u h and a particular solution u p .It will depend on an initial condition u(r 0 ).The initial condition of the radial trajectory is from which we deduce u(t = 0) = 0 4 .The homogeneous version of equation ( 44) is equivalent to where k 0 is an arbitrary integration constant.We follow a standard approach and look for a particular solution of the form u p = k(r, a)u h .Plugging this into (44) gives Combining our results, we obtain This is valid for any spin, for any location of the secondary and for any flux P GW .This analytic result will allow us to determine what the dominant source of the spin dependence is in different regions of the trajectory.
In figure 1 we show the near perfect agreement between the solution to (49) and our numerical calculation of ∂ a r using finite difference method the method used to calculate year-long trajectories used for our Fisher matrix results in later sections, for both moderately and rapidly rotating primaries.Following [11], we express the energy flux as a relativistic correction factor, Ė, times the leading order Newtonian flux Plugging this into Eq.( 49) gives we see that the first and third terms are kinematic, i.e., driven by geodesic physics, whereas the second is dynamical, i.e., driven by the energy flux.Comparison between these terms at different stages of the inspiral, as a function of the spin, can help us to determine what the driving source of spin dependence is in each case.In the next subsection, we investigate the contribution of both the geodesic and radiation reactive terms to ∂ a r.49) with k 0 = 0 corresponding to ∂ a r( r0 ) = 0.In both plots, the solid colours (blue and violet) are ∂ a r calculated using a fifth order stencil method.In each plot, the intrinsic parameters given in the titles.

B. Comparison of radial evolution for moderate and near-extremal black holes
Despite the universality of ( 49) or ( 52), the dependence on the energy flux makes it not feasible to analytically integrate ∂ a r along the entire secondary trajectory.However, we can integrate (49) in specific regions of the secondary trajectory.
It is possible to prove that d∂ a r/dr < 0 and hence that ∂ a r grows monotonically over the inspiral.It is therefore natural to study the behaviour of ∂ a r close to ISCO, where its contribution to the Fisher matrix (40) will be maximal.We first compare the kinematic and dynamical contributions to (53).Using results from the BHPT, we have numerically calculated the spin derivative of Ė for two primaries with spin parameters a = 0.9 and a = 1 − 10 −6 .These are compared with the kinematic sources in (53) in figure 2. These figures show that for both spin parameters.This suggests it is the kinematic sources in (53) that drive the spin dependence of the secondary trajectory, particularly close to ISCO.Although we have only verified it for two choices of spin parameter, we will assume this approximation holds for any spin parameter a ≥ 0.9.We first consider moderately spinning black holes close to ISCO.Dropping the dynamical contribution to (53), we can compare the two remaining terms.The angular velocity piece is bounded and order one, but ∂ r Ẽ tends to zero at ISCO.This means that ∂ 2 ar Ẽ/∂ r Ẽ dominates close to ISCO, allowing us to use the approximation Since, for moderate spins, the variation of Ω and Ė with radius close to ISCO is negligible compared to the variation in ∂ 2 ar Ẽ, we will approximate them by their values at risco .This allows us to integrate (49) to give the spin dependence of the radial trajectory for moderately spinning black holes where k mod is an arbitrary constant.Since it follows from eq. (A5) in [40] that ∂ r Ẽ(r isco ) = 0.For moderately rotating primaries and near ISCO, we can expand  53) for a spin of a = 0.999999.The bottom plot is the same but for a spin parameter of a = 0.9.Notice that in these two cases the kinematical quantities dominate over the relativistic correction terms.
Using these expansions in Eq. ( 56) we deduce ∂ a r = kmod /(r − risco ) + ∂ a risco , where kmod = k mod Ω10/3 isco Ė0 (a, risco )/∂ 2 r Ẽ(r isco ).Assuming that r0 is sufficiently close to risco that this approximation holds throughout the range [r isco , r0 ], we can use the boundary condition (45) to determine kmod = ∂ a risco (r isco − r0 ) and hence We now repeat this analysis for near-extremal primaries.Near ISCO, the energy flux can be approximated by the NHEK flux (x ≡ r − 1 1) Using this approximation, there is no explicit spin dependence and so the ∂ a P GW term in Eq. ( 47) vanishes.Expanding (57) for x = r − 1 1 and 1, the denominator involves while the numerator has the expansion We conclude Using this approximation in Eq. ( 47) we find where we have used x isco ≈ 2 1/3 2/3 and P GW defined in (60).This is valid for x 1 and includes the corrections due to x ∼ x isco .Assuming r0 is close to ISCO, so that the initial condition (45) holds, we conclude that the spin dependence in the near-ISCO region of a nearextremal black hole is Figure 3 compares ( 59) and (65) to the full ∂ a r computed numerically without using the near-ISCO approximations.We see that the approximations are very accurate in the region close to the ISCO where they are valid.
Before continuing, we will comment further on the choice of flux (12) instead of (11).The latter has an explicit dependence on r+ = 1 + .Consequently, it carries an additional spin dependence.In particular, ∂ a r+ = −a/ .Thus, for near-extremal primaries this spin dependence can induce extra diverging sources for ∂ a r with a very specific sign.We can easily compute their effects by integrating the ODE with such an energy flux source.The result one finds does not agree with the numerical evaluation of ∂ a r generated from the BHPT, which computes the exact flux5 .We conclude that (12) appears to capture the spin dependence of our observable (the amplitude of the gravitational wave) more accurately than (11).This is in fact the reason we chose to work with (12).59) and (65).The purple and blue curves are the true solutions to ∂ a r obtained numerically without near-ISCO simplifications.We see both approximations capture the leading order behaviour of the spin derivative of the radial trajectory very well.
Let us close this discussion with a brief comparison between the analytic results for moderate and nearextremal spins.We write r − risco ∼ δ > η 2/5 , the latter inequality ensuring that we avoid entering the transition region [40,54].Expanding Eq. ( 61) we find for near-extremal black holes The first term is dominant unless |δ x 2 isco ∼ 4/3 .The constraint δ > η 2/5 ensures this is only violated if > η 3/10 .This will be satisfied for all the cases that we consider in this paper, but we emphasise this is not a physical constraint.When this constraint is violated, additional terms become important in the expansion which we have ignored, and these ensure that ∂ r Ẽ → 0 at risco .We conclude the scaling of ∂ r Ẽ is δ/ 2/3 for near-extremal black holes, compared to δ for moderate spins.
It follows using ( 59) and (65) that The spin dependence on the radial trajectory for nearextremal primaries is larger than for moderately rotating ones.

C. Precision of spin measurement
In the previous section we showed the effect that the spin parameter has on the radial trajectory.This was achieved by studying the general linear ODE for ∂ a r, Eq. (49).By arguing that the kinematic terms dominate the behaviour of ∂ a r, for both the near-extreme and moderately spinning black holes, analytic solutions were found near the ISCO.We were able to conclude that ∂ a r grows much more rapidly close to the ISCO for near-extreme black holes than for moderately spinning black holes.We also emphasise that Eq. (54) shows that corrections to ∂ a r of the form ∂ a Ė are subdominant.We now explore the consequences of these results for the precision of spin measurements, computed using the Fisher Matrix formalism.
Due to the large number of observable gravitational wave cycles that are generated while the secondary is within the strong field gravity region outside the primary Kerr black hole, extreme mass-ratio inspirals will provide measurements of the system parameters with unparalleled precision [1].In particular, it has been shown that our ability to constrain the spin parameter a is expected to be O(10 −6 ) [3,4,15,55].It has also been shown that the measurements are more precise for prograde inspirals into more rapidly spinning black holes, when the secondary spends more orbits closer to the event horizon of the primary (see Fig. (11) in [55]).In subsequent sections we show through numerical calculation that spin measurements are even more precise for EMRIs into near-extremal black holes.We now try to understand this result using Eq.(40).
Inspection of (40) suggests there are two main effects: the dependence on t2 and the dependence on (∂ a r) 2 .First, the fact that t ∼ O(η −1 ) follows from integrating (41), and therefore the contribution to the Fisher matrix due to t2 is large and scales like η −2 .Second, ∂ a r is monotonically increasing as the secondary spirals inwards.Thus, its maximal contribution comes from the region close to ISCO, which supports results in [55].Eq. (66) shows this contribution is largest in the last stages before entering into the transition regime.As changes in Ω close to ISCO are negligible, the factor ( Ωt ) 2 is, approximately, the square of the number of cycles, a proxy widely used in the literature in discussions of the precision of measurements.Our estimate (40) confirms this intuition and shows the spin precision will be further increased by large values of the radial spin derivative, ∂ a r.

D. Comparison of spin measurement precision for moderate and near-extremal black holes
The Fisher matrix estimate (40) depends on the spin derivative of the radial evolution, on the duration of the inspiral and on the energy flux.Eq. (66) shows that, at a fixed distance to the corresponding ISCO, ∂ a r is larger for a near-extremal primary than for a moderately rotating primary.As a consequence of time dilation near the black hole horizon, Ė → 0 near the ISCO for near-extremal primaries, but remains finite for moderately rotating ones.This means that the energy flux for near-extremal inspirals is much smaller than that for moderate spins, but the duration of the inspiral is longer.However, we can write (40) as an integral over the BL radial coordinate r.In that case, the integrand is proportional to close to the relevant ISCO.While the energy fluxes are much smaller for near-extremal inspirals, the ratio of energy fluxes appearing above is an order one quantity for all spin parameters.The expression above is therefore a product of factors that have been argued to be either of comparable magnitude or much smaller for moderately rotating primaries.We therefore expect the precision of spin measurements to be much higher for near-extremal EMRIs.
A quantitative comparison between the near-extremal and the moderately spinning sources requires a precise calculation of the ratio (37) computed along the entire respective trajectories.In general, this is a hard analytic task since both energy fluxes ĖGW and Ė∞ m must be handled through numerical means and long observations of inspirals (starting in the weak field) require calculations performed in the frequency domain where S n (f ) shows non-trivial (non-constant) behaviour.This would be no more straightforward than direct numerical computation of the Fisher Matrix and so we do not pursue it here.
For any sources whose trajectory lies entirely lie in the near-ISCO region, these analytic approximations allow us to compute the ratio (37) reliably.This can be exploited to obtain an analytic approximation to the Fisher Matrix for such sources and this calculation will be pursued elsewhere.Additionally, earlier arguments tell us that it is the near-ISCO regime that dominates the spin precision and so these expressions are sufficient to understand the increase in spin precision seen for near-extremal inspirals.
a. Near-extremal source.From (65), it follows Since √ r∂ a r grows fast and the rate of change of r and Ω is small, near-ISCO, we can approximate (40) by Here, tcut is the time at the end of the integration, where x = x cut .Our approximations break down when the transition regime breaks down, so we can assume x cut ∼ η 2/5 + 2/3 , which is a small quantity.Using d Ẽ∞ m /d t = η C∞m x and assuming x 0 ≥ x x isco , so that the trajectory can be approximated by x( t) ≈ x 0 e −y with y = α η t 2 ( CH + C∞ ) η t, the Fisher matrix reduces to with e ycut = x 0 /x cut and where in the last step we used x 0 x cut .b. Moderately spinning source.Using the same kind of approximations as above, but taking into account the different energy flux and different trajectory one can approximate the Fisher matrix for moderate spins by where δ ≡ rcut−risco r0−risco < 1 and c. Ratio of Fisher matrices.Within these approximations, the ratio (37) now reduces to The most relevant feature for our current discussion is the quotient dependence The first two denominator factors increase the ratio, since x cut 1 and r0 − risco < 1.The last could in principle be large, due to the logarithmic term.However δ and x cut have similar scaling and therefore x cut F (δ) 1. We deduce that the spin component of the Fisher Matrix is much larger for near-extremal inspirals than for moderate spins.This is confirmed by the numerical results that will be reported in subsequent sections.
We finish by noting that the Fisher matrices increase in magnitude as the trajectory is cut off closer to risco .In the case of moderate spin, we already noted the logarithmic dependence of F (δ) as δ → 0. This has previously been observed in the literature, see for example Fig. (11) in [55].For near-extremal EMRIs, if x cut ∼ x isco ∼ 2/3 , then for fixed x 0 and as → 0 the spin Fisher matrix scales as Γ aa ∼ (log( )/ ) 2 .We deduce that observing the latter stages of inspiral is im-portant for precise parameter measurement, for any primary spin.
In summary, we have derived an analytic approximation, valid close to ISCO, for the spin component of the Fisher Matrix.This indicates that this component is much larger for near-extremal spins and therefore we expect much more precise measurements of the spin parameter in that case.The approximation depends sensitively on certain quantities, such as the cut-off radius, x cut , that are somewhat arbitrary.However, for any choice the near-extremal precision is a few orders of magnitude better.This provides support for the numerical results that we will obtain in Sec.(VI), which show a similar trend.

IV. WAVEFORM GENERATION
In this section we provide more details on how we construct the waveform model used to compute the Fisher Matrix in the next section.The waveform model was previously given in Eq. ( 16) and Eq.(20).Here, we de-scribe how the various terms entering these equations are evaluated.

A. Energy Flux
Both waveform models (16) and (20) depend on the radial trajectory r( t, a, η, Ė).The amplitude evolution using the Teukolsky formalism depends on the spheroidal harmonics −2 S am Ω ml (θ, φ) and Teukolsky amplitudes at infinity Z ∞ ml (r, a).The energy flux at infinity Ė∞ ml (r, a) is related to the Teukolsky amplitudes Z ∞ ml through equation (26).Thus, to accurately generate the waveforms ( 16) and ( 20) far from the horizon where near-extremal simplifications can not be made, the various radiation reactive terms Z ∞ ml , Ė( Ė), Ė∞ m ( Ė∞ m ) have to be handled numerically.This section outlines our numerical routines to do so.
We use the BHPT to calculate the first order dissipative radial fluxes ĖGW for a = 1 − {10 −i } i=9 i=3 from which Ė in Eq. ( 51) can be computed.We used the Teukolsky mathematica script in the toolkit and tuned the numerical precision to ∼ 240 decimal digits to avoid numerical instabilities when computing ĖGW in the near-horizon regime for rapidly rotating holes.For moderately spinning holes a 0.999, we used the tabulated data in Table II of [11].
Each coefficient appearing in Eq. ( 19) is itself a sum over l modes, Ė∞ m = ∞ |l|=m Ė∞ ml .Both the sum over l and the sum over m in Eq. ( 19) can be truncated without appreciable loss of accuracy.As discussed in [32], near-extremal EMRIs require a significant number of harmonics to be included to obtain an accurate representation of the gravitational wave signal.To illustrate for a high spin of a = 1 − 10 −9 , we used the BHPT to compute Ė∞ ml for harmonics |m| ≤ l ∈ {2, . . ., 15}. Figure 4 illustrates the convergence as the number of harmonics is increased.
Based on these results, we go further by including harmonics with l ≤ l max = 20 to calculate the total energy flux using the Teukolsky package in the BHPT.In the same numerical routine, we compute Ė∞ m = lmax |l|=m Ė∞ ml using l max = 20 for m ≤ 20.These formulas are rearranged to obtain Ė and Ė∞ m using ( 51) and ( 14).Finally for our Teukolsky based waveforms used in numerics section VI, we use the BHPT to extract the Teukolsky amplitudes Z ∞ ml (r, a) and build an interpolant for each viewing angle (θ, φ) = (π/2, 0) and (θ, φ) = (0, 0).To summarise, we use (77) in (20) to compute Fisher matrices numerically in section VI.To aid our analytic study, we use the computed Ė∞,m in the waveform model ( 16) when evaluating the ratio (74).

B. Radial trajectory & Waveform
The radial trajectory can be constructed by numerically integrating the ODE (13) using an interpolant for Ė(r) and suitable initial conditions.As before, we use the spin independent initial condition r( t0 = 0) = r0 .Fig. (5) shows some example radial trajectories for various spin parameters, computed using flux data from the BHPT.In the high spin regime, the exponential decay of the radial coordinate is prominent as discussed in [24,56].Throughout our simulations, the observation ends after a fixed amount of time, chosen such that this is before the transition to plunge for all parameter values used to compute the Fisher Matrix.This is important to avoid introducing artifacts from the termination of the waveform, given that the transition to plunge is not properly included in this waveform model.It is clear from Figure (5) that larger the spin parameter, the longer the secondary spends in the dampening regime.See equation (22) of [24] for further details.
The spin dependence of the radial evolution can be calculated by integrating (13) and then taking numerical derivatives.We consider two reference cases, both with component masses µ = 10M and M = 2 × 10 6 M , but  The higher the spin parameter, the more time the secondary spends in the throat before plunge.The lower panel shows the corresponding inspiral trajectory.The dampening is clearly shown when the primary is near maximal spin, as seen in [24].
with different spin parameters a = 0.9 and a = 1−10 −6 .We compute one year long trajectories, with r(0) = 5.08 in the first case and r(0) = 4.315 in the second.The spin derivative of the radial evolution can be calculated by perturbing the spin and using the symmetric difference formula for δ 1 Figure 6 plots the quantity |∂ a r| 2 appearing in the Fisher matrix estimation (40).By inspection, it is clear that |∂ a r| 2 is largest when the spin parameter is close to unity and when the radius is close to risco , matching our analytical conclusions using approximations (65) and (56).
Using the semi-analytic model ( 16) we now evaluate the estimate (74), for the same two systems, but different r 0 to ensure that the assumptions made in deriving Eq. (74) still hold (r 0 = 2.85 for a = 0.9 and r0 = 1.2 for a = 1 − 10 −6 ).We choose termination points rcut = risco + λ with λ ∼ {λ ext = 10 −4 , λ mod = 10 −2 }, just Spin dependence on radial evolution a = 0.999 999 a = 0.9 Figure 6: The blue curve is ∂ a r for a = 0.999999.The orange curve is ∂ a r for a = 0.9.Notice that the spin dependence on r grows rapidly in the near-ISCO region of the rapidly rotating hole.
outside the transition region.Finally, the expression C ∞,m was calculated using the high_spin_fluxes.nbmathematica notebook in the BHPT, including harmonics up to m = 10.We find the ratio to be giving a rough estimate that the spin precision increases by at least two orders of magnitude for these two sources.This verifies claims made in section III C. When correlations with other parameters and the shape of the PSD are ignored, we predict a precision on the spin parameter roughly two orders of magnitude higher than for moderately spinning black holes.
To generate gravitational waveforms for the numerical study we use the Teukolsky waveform model (20).The waveform depends on parameters θ = {a, r0 , µ, M, φ 0 , D}.We will consider two classes of near-extremal source, differentiated by the magnitude of their component masses and mass ratio.The first "heavier" source has parameters and the second "lighter" source has where D edge and D face refer to the distance if each source is viewed edge-on/face-on respectively.The distances are fine tuned6 so that we achieve a signal to noise ratio of ρ ∼ 20.This is discussed later in section V.The lighter source is sampled with sampling interval ∆t s ≈ 4 seconds and the heavier one with ∆t s ≈ 25 seconds.We note here that ∆t s = M ∆ t where ∆ t is the dimensionless sampling interval used to integrate (13).The sampling interval is chosen from Shannon's sampling theorem such that ∆t s < 1/(2f max ), where are the highest frequencies present in the waveform for the edge-on and face-on cases respectively.To illustrate, near-extremal waveforms with parameters θ light for both edge-on and face-on viewing angles are plotted in Fig. (7).The lighter source is interesting because it exhibits both an "inspiral" regime and a exponentially decaying regime that we will refer to as "dampening".The heavier source is interesting because the dampening regime lasts more than one year and so the signal is in the dampening region for the entire duration of the observation.In the next section, we discuss detectability of these two types of sources by LISA.

V. DETECTABILITY
The LISA PSD reaches a minimum around 3mHz, and is fairly flat within the band from 1 to 100mHz.For an edge-on near-extremal inspiral with primary mass of ∼ 10 7 M , the dominant harmonic has a frequency of ∼ 3.2mHz at plunge, while the m = 20 harmonic has frequency of 64 mHz.Such heavy sources are thus ideal systems for observing the near-ISCO dynamics.For the lighter mass considered, 2 × 10 6 M , the near-ISCO dynamics are at frequencies a factor of 5 higher, where the LISA PSD starts to rise.While the near-ISCO radiation will still be observable for these systems, its relative contribution to the signal will be relatively reduced.We therefore expect to obtain more precise spin measurements for the heavier of the two reference systems.
The discrete analogue of the optimal matched filtering SNR defined in Eq. ( 29) Here N is the length of the time series, ∆t s the sampling interval (in seconds) and f i = i/N ∆t s are the Fourier frequencies.In Eq.( 83), the discrete time Fourier transform (DTFT) h(f j ) is related to the CTFT through ĥ(f ) = ∆t s • h(f ).To avoid problems with spectral leakage, prior to computing the Fourier transform, we smoothly taper the end points of our signals using the Tukey window here n is defined through tn = n∆ t.The tunable parameter α defines the width of the cosine lobes on either side of the Tukey window.If α = 0 then our window is a rectangular window offering excellent frequency resolution but is subject to high leakage (high resolution).If α = 1 then this defines a Hann window, which has poor frequency resolution but has significantly reduced leakage (high dynamic range).For the heavier source, we use α = 0.25 to reduce leakage effects significantly and frequency resolution is not a problem since the fre- Edge on: ( , ) = ( /2, 0) Figure 7: A near-extremal waveform with parameters θ light viewed face-on (left) and edge-on (right).The dampening region lasts ∼ 55 days.The edge-on case is asymmetric due to the large number of l = 20 modes and shows prominent relativistic beaming near the ISCO as observed in figure 3.b) of [24].
form.Computing the SNR in this way gives ρ ∼ 20 for the light and heavy sources respectively when viewed both edge-on and face-on under the configuration of parameters θ light and θ heavy .In all cases we marginally exceed the threshold of ρ ≈ 20 which is typically assumed to be required for EMRI detection in the literature [4,55].
As mentioned above, the lighter source exhibits two regimes of interest -the initial gradually chirping phase, where the waveform resembles those for moderately spinning primaries, and then the exponentially damped phase while the secondary is in the near-horizon regime.It is natural to ask what proportion of the SNR, and later what proportion of the spin measurement precision, is contributed by each regime.For both edge-on and face-on systems, we separate the two parts of the waveform using Tukey windows and compute the SNR contributed by each part to find For the face-on source, there is just a single dominant harmonic, and the frequency of this harmonic is such that it lies in the most sensitive part of the LISA frequency range.This helps to enhance the relative SNR contributed by the dampening region.The edge-on source, by contrast, has multiple contributing harmonics, which are spread over a range of frequencies, and the proportional contribution of the dampening region to the overall SNR is therefore diminished.
For a non-evolving signal the SNR accumulates like √ T obs , where T obs is the total observation time.The predampening regime lasts 308 days, and so from duration alone we would expect a fraction 308/365 ≈ 93% of SNR to be accumulated there.The difference to what we find above is explained by differences in amplitudes of the individual harmonic(s).The heavier system is within the dampening regime throughout the last year of inspiral and so all of the SNR of ρ ∼ 20 is accumulated there.This may seem counter-intuitive given the exponential decay of the signal during the dampening regime.However, the exponential decay rate is relatively slow, a large number of harmonics contribute to the SNR and the emission is all within the most sensitive range of the LISA detector.This is clear from looking at the time-frequency spectrogram of the heavier signal shown in Fig. (8).What we learn from this figure is that there are a significant number of harmonics that have comparable power to the dominant m = 2 harmonic.We see also that the angular velocity at each harmonic, and thus f m , shows little rate of change for M ∼ 10 7 and η ∼ 10 −6 .This is consistent with [24,32], where it was shown that a large number of m harmonics is required to produce an accurate representation of the gravitational wave signal for a near-extremal EMRI, particularly for near edge-on viewing angles.For moderately spinning black holes a ∼ 0.9 there are not as many dominant harmonics, so those waveforms are cheaper to evaluate.
We are now ready to move on to compute Fisher Matrix estimates of parameter measurement precisions.This will be the focus of the next section.

VI. NUMERICS: FISHER MATRIX
We now compute (33) numerically without making the simplifying assumptions used in Sections III B and III C. We will use one simplification, which is to ignore the spin dependence in Ė, Z ∞ lm (r, a) and −2 S am Ω lm (θ, φ) and fix these at the values computed for a = 1 − 10 −9 using the BHPT.We argued in Eq. ( 54) that the spin dependence of the flux correction is a subdominant contribution in the near-ISCO regime, and this is further justified in Appendix B (see Fig. 16 in particular).While ∂ a Ė does grow as the ISCO is approached, it remains sub-dominant to the spin dependence of the kinematic terms.This approximation is probably conservative in the sense that we are removing information about the spin from the waveform model and so the true measurement precision is most likely higher.Nonetheless we expect this to be a small effect, and have verified that relaxing this assumption does not significantly change the result for the heavier reference source (see Figure 9).We note that we make this assumption only for computational convenience.Waveform models used for parameter estimation on actual LISA data should use the most complete results available to ensure maximum sensitivity and minimal parameter biases.
To compute the waveform derivatives required to evaluate (33), we use the fifth order stencil method for δx 1 and f i = f (x + iδx).To avoid numerical instability of ∂ a h for the near-extremal spin values of a ≤ 1 − 10 −9 , we ensure that δx < 1 − a so the perturbed waveform does not have spin exceeding a = 1.We further assume that ∂ a tend is zero so there is no spin dependence on the total observation time.
Fisher matrix estimates of parameter measurement precisions for all three sources viewed edge-on are shown in Figure 9.We do not present the results for a face-on observation as they are near-equivalent to the measurements presented in figure 9 for equivalent SNR.
We see from this figure that we should be able to constrain the spin parameter of near-extremal EMRI sources to a precision as high as ∆a ∼ 10 −10 , even when accounting for correlations amongst the waveform parameters.This is true for both the lighter and the heavier sources viewed edge-on and face-on, with a constraint a factor of a few better for the heavier source.The right panel of the figure compares the contribution to the measurement precision for the lighter source from the two different phases of the signal.We see that the high spin precision comes almost entirely from the observation of the dampening regime and this phase of the signal contributes much more information than we would expect based on its contribution to the total SNR 7 .
The spin measurement precision for the near-extremal systems is three orders of magnitude better than for the system with moderate spin, while all other parameter measurements are comparable.
Comparing to the exact Fisher matrix result with spin dependence included in all the various terms, we see that the two precisions are almost identical : the exact result offers precisions that are marginally better in comparison to our approximate result (removing spin dependence from the corrections).This Figure thus justifies ignoring the spin dependence of Ė, since relaxing that assumption makes almost no difference to the results.This numerically confirms our belief that the spin dependence in the corrections to the fluxes are subdominant in the analysis leading to (40).In the same plot 9, we also compare results of near-extremal black holes to moderately spinning holes.A direct comparison shows an increase in the spin precision by ∼ 3 orders of magnitude, which agrees with the intuition given by the earlier analytic analysis, Eq. (79).
To our knowledge, these are the first circular and equatorial parameter precision studies for EMRIs that have employed Teukolsky-based adiabatic waveforms, rather than approximate waveform models (or "kludges"), which have been used for many studies [3,4,15].Comparing our results for the moderately spinning system to these previous studies, we find that our results are very comparable, but a factor of a few tighter.This could be because we are including only a subset of parameters and ignoring the details of the LISA response, or because we have a more complete treatment of relativistic effects.A more in depth study addressing both of these limitations would be needed to understand the origin of the differences.However, the agreement between our results and previous studies is sufficiently close, and considerably less than the difference we find between the moderate and near-extremal spin cases, to give us confidence that our results are not being unduly influenced by these simplifications.
In Figure (10) we show how the parameter estimation precision for the source with parameters θ light changes as we vary the spin parameter, while keeping all other parameters unchanged.We present results for both face-on and edge-on viewing angles.This shows that while the measurement precision for most of the parameters is largely independent of spin in the near-extremal regime, the spin precision steadily increases as a → 1.We note that even at a spin of 1 − 10 −9 , the measurement precision satisfies the constraint ∆a < |1 − a| and The black diamonds show the precisions obtained when including the spin-dependence of the relativistic corrections, Ė in the waveform model for the heavy source.(Right plot) Parameter measurement precisions for the source with parameters θ light , computed using the full waveform (blue asterix), only the inspiral phase (blue dot) and only the dampening phase (green diamond).
therefore a LISA EMRI observation would be able to resolve that the system was not maximally extremal, i.e., that a < 1.We stop at 1 − 10 −9 since the derivative using Eq.( 87) begins to misbehave.Due to large condition numbers, inverting Fisher matrices for EMRI sources is a highly non-trivial task.In appendix C, we provide multiple diagnostic tests of our Fisher matrix algorithm and verify that, in the single parameter case, the spin parameter precision is a suitable representation of the 1σ width of the Gaussian likelihood as shown in figure 18.These single parameter tests of the Fisher matrix are useful tests to verify that a single parameter algorithm yields sensible results.However, real instabilities of the numerical procedure are prominent the moment the inverse of the Fisher matrix is performed when correlations are present.Hence, it is both necessary and sufficient to verify our Fisher matrix calculations using an independent procedure.The next section is dedicated to performing a parameter estimation study on both near-extremal EMRIs with parameters θ light and θ heavy .

VII. NUMERICS: MARKOV CHAIN MONTE CARLO
The Fisher Matrix is a local approximation to the likelihood, valid in the limit of sufficiently high signalto-noise ratio.We can verify that this local approximation is correctly representing the parameter measurement uncertainties by numerically evaluating the likelihood using Markov Chain Monte Carlo.To reduce the computational cost of these simulations we use a faceon viewing profile and thus only consider the m = 2 harmonic.We have shown in figure (10) that parameter precision measurements are not largely dependent on the choice of viewing angle for the lighter source.We have further verified this claim for the heavier source.

A. Markov Chain Monte Carlo
Markov Chain Monte Carlo (MCMC) methods were developed for Bayesian inference to sample from the posterior probability distribution, p(θ|d), which is given by Bayes' theorem as log p(θ|d) ∝ log p(d|θ) + log p(θ) where p(d|θ) is the likelihood function, and p(θ) is the prior probability distribution on the parameters.In our context the likelihood is given by Eq. ( 27) and we will assume independent priors such that (90) We generate a data set d(t) = h(t; θ tr )+n(t) by specifying the waveform parameters, θ tr , of the injected signal and generating noise in the frequency domain We use MCMC to sample from the posterior distribution (90), employing a standard Metropolis algorithm with proposal distribution q(θ |θ i−1 ) equal to a multivariate normal distribution, centred at the current point and with a fixed covariance.We take flat priors on all of the waveform parameters, since the goal is to check the validity of the Fisher matrix approximation to the likelihood.The algorithm proceeds as follows 1.We start the algorithm close to the true values θ 0 = θ tr + δ for ||δ|| 1.For iteration i = 1, 2, . . ., N 2. Draw new candidate parameters θ ∼ q and generate the corresponding signal template h(t; θ ).
3. Using (90), compute the log acceptance probability We note that we are using a symmetric proposal distribution and so the usual proposal ratio is not required.Since we know the true parameters we can start the algorithm in the vicinity of the true parameters and do not need to discard the initial samples as burn-in, allowing us to generate useful samples more quickly.In principle, the MCMC algorithm should converge for any choice of proposal distribution, but proposals that more closely match the shape of the posterior should lead to more rapid convergence.As we expect that the proposal should be approximated by the Fisher Matrix, we set the covariance matrix of the normal proposal distribution to be equal to the inverse Fisher matrix, evaluated at the known injection parameters.

B. Results
We compute MCMC posteriors for the two signals h({θ heavy , θ light }; t) for the waveform model (20) for the face-on case only.As before, we construct waveforms ignoring the spin dependence in Ė, the Teukolsky amplitudes Z ∞ lm and the spheroidal harmonics −2 S a Ω lm (θ, φ).We evaluate these for a fixed spin parameter of a = 1 − 10 −9 .
We remind the reader the main drive for the tight constraints on the spin parameter is due to the spin dependence induced through the kinematic terms present in (20), as discussed in section VI, and not its dynamical terms, justifying our approximation.
The priors on a, φ 0 and D for both sources were The priors on µ, M and r0 were chosen differently for the heavy and light source as The prior on a ensures that we do not move outside the range in which our approximations are valid, a 0.9999.The tight priors on the individual component masses helped to improve the computational efficiency of our algorithm.However, there was no evidence of the MCMC chains reaching the edges of the priors in our simulations, so we are confident these restrictions are not influencing the results.
Evaluating the likelihood for EMRI waveforms is an expensive procedure.In order to obtain a sufficient number of samples from the posterior, we used high performance computing facilities and ran 20 unique chains for N = 40, 000 iterations.All chains analysed the same input data set, but with different initial random seeds.This ensures that the dynamics of the chains are different but the noise realisations are the same for each MCMC procedure.
The marginal posterior distributions and twodimensional contour plots for the two sources are shown in figures (11) and (12).These plots confirm the high precisions of parameter measurements that were seen with the Fisher Matrix.The relative uncertainties ∆θ/θ are similar for the two sources, although we can measure the spin parameter more precisely for the heavier source.For the most part the posteriors are unimodal, apart from the spin posterior of the lighter source.We have verified that the secondary modes are real features of the likelihood, and correspond to the waveform phase shifting by one cycle within the late dampening regime.We also note that shifts in the peak of the posterior away form the true value are larger for the heavier source than for the lighter source.This appears to be due to the particular noise realisation.For other noise realisations the noise-induced biases for the heavier source are smaller.For noise-free data sets, we find posterior distributions peaked at the true parameters, as expected.
The primary reason for doing the MCMC simulations was to verify the Fisher Matrix results found earlier.In figures 13 and 14, we plot the marginalised posteriors on the parameters {θ light , θ heavy } alongside a Gaussian distribution with variance given by the Fisher matrix and centred at the mean value of the posterior distributions p(θ|d).
These results nicely confirm the accuracy of the Fisher matrix results for these sources.In each case, the 1σ precision predicted by the Fisher matrix is slightly smaller than the width of the numerically computed posterior.This is to be expected as the Fisher matrix also provides the Cramer-Rao lower bound on parameter uncertainties.However, the difference is very small.We are thus confident that all of our Fisher matrix predictions are accurate, including the exploration of parameter space shown in Figure 10.We conclude that even at the near-threshold SNR of this source, ρ ∼ 20, the Fisher matrix can be used to confidently estimate the precision of parameters for near-extremal sources.

VIII. CONCLUSION
This work has shown that we expect tighter spin constraints on the Kerr spin parameter for EMRIs when the primary is rapidly rotating rather than moderately rotating.We argued that the spin precision is sourced through the rapid growth of the spin dependence on the radial trajectory throughout the near-ISCO regime.It was shown that radiation-reaction effects are subdominant to geodesic motion, which is what ultimately drives the excellent constraint on the spin parameter.
This work has tried to answer the question of how well LISA can constrain the parameter of an EMRI source, conditioned on the primary rapidly rotating.In the earlier sections of this paper, we were able to derive a general linear ODE whose solution describes the spin dependence on the radial trajectory.Using this ODE, we were able to argue that this spin dependence is much more significant for near-extremal primaries rather than moderately spinning primaries a ∼ 0.9.It was also shown that the spin dependence on the radial evolution is dominated more by geodesics, rather than the radiation reactive terms.
Using reasonable assumptions, analytic Fisher scalars were derived for both rapidly rotating and moderately spinning black holes.It was then shown, explicitly, that the precision on spin parameters for near-extreme black holes should exceed that of moderately spinning ones by a few orders of magnitude.The exceptional precision is governed by the growth of the spin dependence of the trajectory near the ISCO.
Given that the kinematical terms dominate the precision, we were able to build a near-extremal waveform model by removing the spin dependence in corrections to the radial fluxes.This meant we could build an interpolant, valid for quasi-circular and equatorial inspirals, into rapidly spinning black holes with a 0.999.Using this approximation, we were able to build a waveform model valid for near-extremal spins that could be used for parameter precision and estimation studies.
Our analysis showed that the LISA optimum masses for parameter estimation and precision studies are heavier mass systems M ≈ 10 7 M .For heavier sources, the emitted frequencies will lie within the minimum of the LISA PSD during the dampening regime; when the signal provides most accurate information about the source.Analysing different parts of the signal, we were able to conclude that the end part of the signal will be "loud" in the data stream for spins approaching near-extremality a → 1.For lighter mass sources M ∼ 10 6 M , the accumulated SNR of the dampened part of the signal is weak.Fully numerical Fisher matrix analysis revealed that we can constrain the spin parameter of a near-extremal EMRI ∼ 3 orders of magnitude higher than moderately spinning black holes a ∼ 0.9.For very near-extremal primaryes a ∼ 1 − 10 −9 , one is able to constrain the spin parameter with precision ∆a ∼ 10 −10 .This is, as far as we know, the tightest constraint on astrophysical orbital parameter found in the literature.
In the final section, with the view of verifying our Fisher matrix calculations, we performed a parameter estimation study on these near-extremal sources.We showed that the signal parameters are able to be extracted with precision comparable to our Fisher matrix estimates.
To conclude, we can safely say that the spin parameter of near-extremal EMRIs can be measured with excellent precision.We are also confident that the increase in precision is governed by the spin dependence on the radial trajectory.If the dampened part of the signal is not observed (where ∂ a r 1), then we cannot make such precise statements on the spin of the primary black hole.
There are a few extensions to this work.The first would be to include the results of the transition from in-spiral to plunge in [40] to our analysis here.We suspect that, although the SNR accumulation would be weak, the extra contribution from ∂ a r could only improve the precision measurement.Another interesting direction to explore would be the effect of systematic errors (through innacurate waveform modelling) in our EMRI PE studies.Through the Cutler and Vallisneri formalism [57], one could investigate whether parameter estimation/precision studies on near-extremal waveforms could be dominated by systematic errors rather than statistical ones induced purely by noise realisations from the detector.Given that we have shown precisions on the spin parameter on the order ∆a ∼ 10 −10 , it would be of utmost importance to have immensely accurate waveform models in order to keep the systematic error to a mini- mum.This would be to ensure that parameter estimates are not affected by (deterministic) biases due to waveform modelling errors rather than (probabilistic) biases due to noise fluctuations of the LISA detector.
Another obvious extension would be to introduce more parameters in the parameter space, such as eccentricity and inclination.In this case, we are unsure whether the kinematic terms would dominate over the flux corrections.As such, we would require a more complex and general waveform model.As seen in [24], eccentric trajectories of secondaries in the vicinity of near-extremal black holes exhibit inverse zoom-whirl behaviour.That is, the "zoom" part of the inspiral (near aperiastron) outputs larger amplitude radiation than the "whirl" phase.Perhaps this extra information, unique to near-extreme EMRIs, could provide even tighter constraints on the spin parameter a.It would be interesting to see whether the larger number of parameters will aid or hinder parameter precision/estimation studies.Finally, due to the precision on both ∆a and ∆M , near-extremal EMRIs would provide stringent tests of GR theories.One of which being the no-hair theorem which states that the Kerr black hole can be uniquely parametrized in terms of two charges; mass and spin.We believe that it should be possible to measure the multipolar moments more accurately when the primary is near-extremal than if it were only moderately spinning.A paper exploring potential tests of general relativity using the near-extremal Kerr spacetime would be useful for the literature.Figure 13: For the source with parameters θ heavy ,we compare the one-dimensional marginalised posterior distributions (orange histograms) to a Gaussian distribution (blue solid line), centred at the posterior mean, and with standard deviation set to the prediction of the Fisher matrix.

ACKNOWLEDGMENTS
We give thanks to Maarten van de Meent for both providing high spin Teukolsky amplitudes and fluxes used throughout this work and for advice on the waveform generation.We also give thinks to Niels Warburton for advice on how to use the BHPTand general comments on the results.Finally, we wish to acknowledge Lorenzo Speri and Andrea Antonelli for useful comments on the manuscript and for many fruitful discussions/support leading to the completion of this work.This work made use of data hosted as part of the Black Hole Perturba-  IV in Ref. [11].The figure on the right shows the GRCs evaluated at a fixed radial coordinate r whilst varying the spin parameter a.The GRCs for a > 0.999 were computed using the BHPT [46].from this is that the behaviour of the relativistic corrections becomes somewhat universal as a → 1.In other words, as the near extremal parameter approaches unity, the corrections approach a constant value.This means that the correction values at some r will not differ much as the spin parameter is changed providing the spin parameter is close to one.That is, spin dependence on Ė is weak "far" from ISCO.
We emphasise that this does not mean that ∂ a Ė ≈ 0 throughout the entire inspiral.Here we have shown three radial coordinates which are considerably far away from the ISCO in the near-extreme spin parameter case.The closer in radial coordinate to the ISCO, the larger the spin dependence in the relativistic corrections.This can be observed in Fig. (16).For a weak gravitational Spin dependence on corrections a = 0.9 a = 0.999 999 Figure 16: Here we compute ∂ a Ė for the two spin parameters a = 0.999999 (orange) and a = 0.9 (blue).The growth in ∂ a Ė for rapidly rotating holes is larger closer to ISCO than for moderately spinning holes.
field, the spin dependence on the flux corrections for near-extreme inspirals is small, approaching Ė → 0 for a → 1 as the ISCO is approached.Much smaller than when compared to moderately spinning holes, where the total energy flux is approximately constant.Only when very close to the ISCO does the spin dependence on Ė grow significantly for near-extremal holes.
A comparison of our interpolant, Ė for a = 1 − 10 −9 with high spinning fluxes computed from the BHPT are given in Fig. (17).A useful approximant for the relativistic corrections Ė(r) is given by Ė(x) = ax arctan(bx c ), x = r − 1 (B1) with constants (a, b, c) = (0.6897912, 1.084803, −0.936685).This approximates the high spin corrections Ė(r, 1 − 10 −9 ) for r 6 with largest fractional error ∼ 1% at the ISCO.It is well known throughout the literature [6,[58][59][60][61][62][63]] that a Fisher matrix analysis, when not handled carefully, has the potential to under-estimate (or, worse, over-estimate) precision measurements.One first has to be sure that the regime of SNR is high enough so that the linear signal approximation applies when truncating the perturbed waveform at first order.Numerical derivatives also exhibit convergence problems and are very problematic if your waveforms are not smooth.One must also make sure to taper each derivative so  that "hard cut-offs" do not feature in the waveform [64].
Finally, one has to invert the Fisher matrix and, if the Fisher matrix is ill conditioned (as normally the case for EMRIs), it could lead to catastrophic errors in in the elements of Γ −1 .The condition number of a matrix is defined through where λ i is an i th eigenvalue of the matrix Γ.In the case of EMRIs, the computation of the fisher matrix is affected by numerical instabilities.A small perturbation to the systems (intrinsic) parameters lead to a large overall change in the waves phase evolution.It is for this reason EMRI observations permit parameter constraints with such high precision.Since waveforms are sensitive to small perturbations of the source parameters, the numerical derivatives ∂ θ h(θ; t) are large and so the elements of the Fisher matrix can be enormous in magnitude.In other words, there is a significant amount of "information" about the parameters θ encoded in the waveforms, which are then reflected by the large elements of the Fisher matrix.However, not all derivatives are large and the differences between the size of the elements for different parameters leads to significantly varying eigenvalues.In our case, the measurements of spin and distance are ∼ 8 orders of magnitude apart.A consequence of this is that the condition number of the Fisher matrix is Cond(Γ) ∼ 10 21 .In light of these known instabilities, we spend the next few subsections providing suitable tests to verify our Fisher matrix in the single parameter case, assuming all other parameters are known perfectly.In order to verify Fisher matrices over multiple parameters, Bayesian techniques like MCMC are required (see sec.VII).
We used three different methods to verify the Fisher matrix calculations without performing a Bayesian analysis: 1. Verification of the linear signal approximation h(a + ∆a) ≈ h(a) + ∆a ∂h ∂a .
2. Overlaps defined through (30) -perturbing the spin parameter by ∆a given by the Fisher matrix should return overlaps close to one.
3. Likelihood -The log-likelihood is maximised at the true parameters θ 0 .If a parameter is perturbed by the Fisher matrix estimate then it should be a measure of the 1σ width of this loglikelihood.
As a proof of principle of these methods, we will compute the precision on the spin parameter for a heavy source θ heavy (see (80)), viewed face-on with ρ ∼ 20.
For a source with this configuration of parameters, we found ∆a NHEK ∼ 2 × 10 −10 numerically.
a. Linear-Signal approximation In the derivation of the Fisher matrix, we used the linear-signal approximation so a first test would be to test whether it is valid in our analysis.To test whether the expansion is valid in the regime of SNR we are considering, we compute the overlap O h(a + ∆a)|h(a) + ∆a ∂h ∂a a=atrue ≈ 1−10 −5 (C2) Hence we conclude that our waveform model at ρ ∼ 20 does not violate the linear-signal approximation.
b. Overlaps In the limit of high SNR, the inner product (28)  where (h|h) = ρ 2 as in (29).Since the SNR is fixed, it's easy to show that (∂ a h|h) = 0 and Γ aa = −(∂ 2 a h|h).It can then be shown for ρ 1 Substituting ∆a = ∆a NHEK into the left hand side of (C3), we numerically find an agreement of ∼ 0.01%.
c. Likelihood In the high SNR limit, the precision pridicted by the fisher matrix should approximate the 1σ width of the likelihood function.Using Eq.( 27), we can write the log-likelihood as log p(d|a) ∝ (d|h) − 1 2 (h|h).
Since the noise realisation in the data stream d(t) = h(t; a) + n(t) induces a bias to the maximum likelihood estimate, and does not affect the likelihood width, we shall ignore the noise in this case.As such, we will consider a zero noise approximation and use d = h(a) with signal templates h := h(a + ∆a; t).Substituting this into the likelihood above we find that log p(d|a) ≈ 1 2 (ρ 2 − 1). (C4) Here we have assumed that Γ −1 aa = ∆a 2 .Calculating log p(d|a) for d = h(a + ∆a) with ∆a our Fisher matrix estimate we found that log p(d|a) ≈ 199.46 which agrees with the above formula, to a precision of 0.05%.
Since we are interested in only one parameter, it is evaluate the likelihood function on a grid of spin parameter values.In doing so, we find Fig.(18).The area between the yellow line and red line is approximately 31.51%, which is a reasonable approximation to the true 1 − σ ≈ 34% width of likelihood.
To conclude these subsections, we are confident that our Fisher matrix approximations in the single parameter study give a good guide to the spin parameter uncertainty.

Figure 1 :
Figure 1: The dashed curves (black dashed and yellow dashed) on each figure is the solution to (49) with k 0 = 0 corresponding to ∂ a r( r0 ) = 0.In both plots, the solid colours (blue and violet) are ∂ a r calculated using a fifth order stencil method.In each plot, the intrinsic parameters given in the titles.

Figure 2 :
Figure 2: The top plot compares the kinematic and radiation reaction quantities given in (53) for a spin of a = 0.999999.The bottom plot is the same but for a spin parameter of a = 0.9.Notice that in these two cases the kinematical quantities dominate over the relativistic correction terms.

Figure 3 :
Figure 3: The yellow dashed and black dashed curves are solutions to (59) and (65).The purple and blue curves are the true solutions to ∂ a r obtained numerically without near-ISCO simplifications.We see both approximations capture the leading order behaviour of the spin derivative of the radial trajectory very well.

= 2 Figure 4 :
Figure 4: Comparison of the total energy flux at infinity (black curve) including different harmonic Ė∞ lm contributions.Note that at r ≈ 1.3, the l = 2 harmonic energy flux Ė∞ 2 contributes ∼ 32% of the total energy flux, whereas including the first l max = 11 harmonics (violet curve) contributes more than ∼ 98%.

Figure 5 :
Figure5: The top panel shows how dr/d t varies with r.The higher the spin parameter, the more time the secondary spends in the throat before plunge.The lower panel shows the corresponding inspiral trajectory.The dampening is clearly shown when the primary is near maximal spin, as seen in[24].

Figure 8 :
Figure 8: Here we plot the spectrogram of h(θ heavy ; t) viewed edge on.We see 20 tracks in the time-frequency plane corresponding to the m ∈ {1, . . ., 20} harmonics.The colorbar shows that the m = 2 harmonic (second lowest track in frequency) is dominant, but that there are several other harmonics which contribute significantly to the radiated power.

Figure 9 :
Figure9: (Left plot) Parameter measurement precision, as estimated using the Fisher Matrix formalism, for the three reference sources, with parameters θ light (green diamonds), θ heavy (purple crosses) and θ mod (blue asterisks).The black diamonds show the precisions obtained when including the spin-dependence of the relativistic corrections, Ė in the waveform model for the heavy source.(Right plot) Parameter measurement precisions for the source with parameters θ light , computed using the full waveform (blue asterix), only the inspiral phase (blue dot) and only the dampening phase (green diamond).

Figure 10 :
Figure 10: We keep θ light \{a} fixed and vary a = 1 − 10 −i for i ∈ {4, . . ., 9} while computing estimates on the precision of the measured parameters using the Fisher Matrix.Results are shown for sources viewed face-on (left) and edge-on (right).

4 . 5 .
Draw u ∼ U [0, 1].(a) If log u < log α we accept the proposed point and set θ i = θ .(b) Else we reject the proposed point and set θ i = θ i−1 .Increment i → i + 1 and go back to step 2 until N iterations have been completed.

Figure 11 :
Figure 11: The diagonal plots represent the marginalised posterior distributions on the parameters θ heavy .The plots below the diagonal are the joint two-dimensional posterior distributions.The red lines indicate the true values of the injected signal.

Figure 12 :
Figure 12: As Figure 11, but now for the source with parameters θ light .

Figure 14 :
Figure 14: As Figure 13, but now for the lighter source with parameters θ light .

Figure 15 :
Figure 15: The figure on the left shows the values of the GRCs Ė(r) evaluated at a fixed spin parameter.The purple dashed line (bottom left of leftmost plot) uses the near-extremal approximation to the flux(11).The other colours are the tabulated values GRCs presented in Table IV in Ref.[11].The figure on the right shows the GRCs evaluated at a fixed radial coordinate r whilst varying the spin parameter a.The GRCs for a > 0.999 were computed using the BHPT[46].

Figure 17 :
Figure17: On the left panel, we have used the BHPT to compute the total energy flux for near-extremal spin parameters with our interpolant (black dashed line) overlaid.The right panel is a zoom in of the left plot giving visual aid as to why our interpolant can be used to approximate a larger regime of spin parameters.

Figure 18 :
Figure 18: blue curve is the likelihood (C5) evaluated on a grid of points.The red the true value a = 1 − 10 −6 and orange the precision measurement predicted by the 1 parameter Fisher matrix ∆a NHEK .