The Fast and the Fiducial: Augmented kludge waveforms for detecting extreme-mass-ratio inspirals

TL;DR

The paper introduces an improved augmented analytic kludge (AAK) waveform model for EMRIs, achieving high accuracy and speed, facilitating detection and data analysis for future space-based gravitational wave observatories like LISA.

Contribution

The paper presents a new version of the AAK model with enhanced accuracy and efficiency, suitable for EMRI detection and data analysis in LISA mock data challenges.

Findings

Waveform overlaps exceed 0.97 with fiducial models

Waveform generation is 5-15 times faster than previous models

Potential viability for EMRI detection with semi-coherent template banks

Abstract

The extreme-mass-ratio inspirals (EMRIs) of stellar-mass compact objects into massive black holes are an important class of source for the future space-based gravitational-wave detector LISA. Detecting signals from EMRIs will require waveform models that are both accurate and computationally efficient. In this paper, we present the latest implementation of an augmented analytic kludge (AAK) model, publicly available at github.com/alvincjk/EMRI_Kludge_Suite as part of an EMRI waveform software suite. This version of the AAK model has improved accuracy compared to its predecessors, with two-month waveform overlaps against a more accurate fiducial model exceeding 0.97 for a generic range of sources; it also generates waveforms 5-15 times faster than the fiducial model. The AAK model is well suited for scoping out data analysis issues in the upcoming round of mock LISA data challenges. A simple analytic argument shows that it might even be viable for detecting EMRIs with LISA through a semi-coherent template bank method, while the use of the original analytic kludge in the same approach will result in around 90% fewer detections.

Tables

Table 1. Table 1: Number of EMRI events when using AAK and AK waveforms in a semi-coherent detection search, for the 12 astrophysical EMRI-population models considered in BEA2017 (M1–12) and two different plunge criteria (Schwarzschild and Kerr innermost stable circular orbits). Events are further divided into four black-hole mass bins, where M 10 := lg ⁡ ( M / M ⊙ ) assign subscript 𝑀 10 lg 𝑀 subscript 𝑀 direct-product M_{10}:=\lg{(M/M_{\odot})} . Event counts are rounded to the nearest 10 for all models except M11, where they are rounded to the nearest 1.

		Number of events in mass range
Plunge	Population							$M_{10}>6$		Total
criterion	model	AAK	AK	AAK	AK	AAK	AK	AAK	AK	AAK	AK
Schwarzschild	M1	20	0	240	10	110	10	10	0	380	20
	M2	30	0	190	10	70	10	0	0	290	10
	M3	20	0	310	10	510	40	40	10	880	50
	M4	70	0	280	20	80	20	0	0	440	40
	M5	0	0	10	0	20	0	0	0	30	0
	M6	20	0	270	10	210	10	20	0	520	30
	M7	230	0	2190	60	1040	100	60	10	3530	180
	M8	0	0	30	0	10	0	0	0	50	0
	M9	20	0	210	0	110	10	10	0	350	20
	M10	30	0	240	10	100	10	10	0	370	10
	M11	0	0	0	0	1	0	0	0	1	0
	M12	230	10	2420	70	1730	130	180	30	4560	230
Kerr	M1	20	0	260	10	230	10	80	10	590	30
	M2	20	0	210	0	160	10	50	10	440	20
	M3	10	0	360	10	1000	60	240	50	1620	120
	M4	50	0	300	20	140	30	30	10	520	70
	M5	0	0	10	0	40	0	40	10	90	10
	M6	20	0	300	10	430	30	200	50	960	80
	M7	190	0	2390	60	2110	150	730	120	5420	330
	M8	0	0	30	0	30	0	10	0	70	0
	M9	20	0	230	0	160	10	30	0	430	20
	M10	30	0	240	10	100	10	10	0	370	10
	M11	0	0	0	0	1	0	0	0	1	0
	M12	190	0	2700	60	3710	210	1830	410	8440	690

Equations

$$\dot{\mathbf{E}}=\mathbf{F}(\mu ,M,\mathbf{S},\mathbf{E});$$

$$\dot{\mathbf{X}}=\mathbf{G}(\mu ,M,\mathbf{S},\mathbf{E});$$

$$h(t)=H(\mathbf{X},\mathbf{R}).$$

$$(\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}}{)}_{\tilde{\mathbf{L}}}:=\left(\frac{(\hat{\mathbf{R}}\cdot \hat{\mathbf{L}})\hat{\mathbf{L}}-\hat{\mathbf{R}}}{{\mathcal{S}}_{\mathbf{L},\mathbf{R}}},\frac{\hat{\mathbf{R}}\times \hat{\mathbf{L}}}{{\mathcal{S}}_{\mathbf{L},\mathbf{R}}},\hat{\mathbf{L}}\right),$$

$$(\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}}{)}_{\mathrm{AK}}:=\left(\frac{\hat{\mathbf{R}}\times \hat{\mathbf{L}}}{{\mathcal{S}}_{\mathbf{L},\mathbf{R}}},\frac{\hat{\mathbf{L}}-(\hat{\mathbf{L}}\cdot \hat{\mathbf{R}})\hat{\mathbf{R}}}{{\mathcal{S}}_{\mathbf{L},\mathbf{R}}},-\hat{\mathbf{R}}\right),$$

$$h_{+}=\sum _{n=1}^{\infty }{h}_n^{+},h_{\times }=\sum _{n=1}^{\infty }{h}_n^{\times }$$

$$h_n^{+}=(1+(\hat{\mathbf{R}}\cdot \hat{\mathbf{L}}{)}^2)(b_n\sin 2\tilde{\gamma }-a_n\cos 2\tilde{\gamma })+(1-(\hat{\mathbf{R}}\cdot \hat{\mathbf{L}}{)}^2)c_n,$$

$$h_n^{\times }=2(\hat{\mathbf{R}}\cdot \hat{\mathbf{L}})(b_n\cos 2\tilde{\gamma }+a_n\sin 2\tilde{\gamma }),$$

$$a_n$$

$$b_n$$

$$c_n=2\mathcal{A}{J}_n(ne)\cos n\Phi ,$$

$$\hat{\mathbf{L}}=\hat{\mathbf{S}}\cos \iota +\left(\frac{\hat{\mathbf{z}}-(\hat{\mathbf{z}}\cdot \hat{\mathbf{S}})\hat{\mathbf{S}}}{|\hat{\mathbf{z}}-(\hat{\mathbf{z}}\cdot \hat{\mathbf{S}})\hat{\mathbf{S}}|}\cos \alpha +\frac{\hat{\mathbf{S}}\times \hat{\mathbf{z}}}{|\hat{\mathbf{S}}\times \hat{\mathbf{z}}|}\sin \alpha \right)\sin \iota ,$$

$$h_I=\frac{\sqrt{3}}{2}(F_I^{+}{h}_{+}+F_I^{\times }{h}_{\times }),$$

$$h_{II}=\frac{\sqrt{3}}{2}(F_{II}^{+}{h}_{+}+F_{II}^{\times }{h}_{\times }),$$

$$F_I^{+}$$

$$F_I^{\times }$$

$$F_{II}^{+}$$

$$F_{II}^{\times }$$

$$(\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}}{)}_{\mathbf{S}}:=\left(\frac{\hat{\mathbf{R}}\times \hat{\mathbf{S}}}{{\mathcal{S}}_{\mathbf{S},\mathbf{R}}},\frac{\hat{\mathbf{R}}-(\hat{\mathbf{R}}\cdot \hat{\mathbf{S}})\hat{\mathbf{S}}}{{\mathcal{S}}_{\mathbf{S},\mathbf{R}}},\hat{\mathbf{S}}\right),$$

$$(\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}}{)}_{\mathrm{NK}}:=\left(\frac{\hat{\mathbf{R}}\times \hat{\mathbf{S}}}{{\mathcal{S}}_{\mathbf{S},\mathbf{R}}},\frac{\hat{\mathbf{S}}-(\hat{\mathbf{S}}\cdot \hat{\mathbf{R}})\hat{\mathbf{R}}}{{\mathcal{S}}_{\mathbf{S},\mathbf{R}}},-\hat{\mathbf{R}}\right),$$

$$h_{+}=\frac{1}{2}{h}_{ij}{H}_{ij}^{+},h_{\times }=\frac{1}{2}{h}_{ij}{H}_{ij}^{\times }$$

$$h_{ij}=\frac{2}{|\mathbf{R}|}\left(P_{ik}{P}_{jl}-\frac{1}{2}{P}_{ij}{P}_{kl}\right){\ddot{I}}_{kl},$$

$$H_{ij}^{+}=({\hat{\mathbf{x}}}_i{\hat{\mathbf{x}}}_j-{\hat{\mathbf{y}}}_i{\hat{\mathbf{y}}}_j{)}_{\mathrm{NK}},$$

$$H_{ij}^{\times }=({\hat{\mathbf{x}}}_i{\hat{\mathbf{y}}}_j+{\hat{\mathbf{y}}}_i{\hat{\mathbf{x}}}_j{)}_{\mathrm{NK}},$$

$$P_{ij}=({\delta }_{ij}-{\hat{\mathbf{z}}}_i{\hat{\mathbf{z}}}_j{)}_{\mathrm{NK}},$$

$$I_{ij}=\mu x_i{x}_j,$$

$$\Sigma \frac{dr}{d\tau }=\pm \sqrt{V_r},$$

$$\Sigma \frac{d\theta }{d\tau }=\pm \sqrt{V_{\theta }},$$

$$\Sigma \frac{d\varphi }{d\tau }=V_{\varphi },$$

$$\Sigma \frac{dt}{d\tau }=V_t,$$

$$V_r(r)=P^2-(r^2+(L_z-aE{)}^2+Q)\Delta ,$$

Full text

I Introduction

The key sources for future space-based gravitational-wave (GW) detectors operating in the millihertz frequency band will include the inspirals of stellar-mass compact objects—typically stellar-origin black holes, but also potentially neutron stars or white dwarfs—into massive black holes (MBHs) at the centres of galaxies. Such systems are called extreme-mass-ratio inspirals (EMRIs), since the typical mass of the inspiralling object is $\sim 10M_{\odot}$ while that of the central MBH is $\sim 10^{4}$–$10^{7}M_{\odot}$ (to give emission at millihertz frequencies). EMRIs occur because an MBH is typically surrounded by clusters of stars, and various dynamical processes (including two-body scattering, tidal splitting of binaries and stripping of giant stars) can lead to compact objects being captured by and subsequently inspiralling into the MBH (see AEA2007 ; A2012 for comprehensive discussions of the astrophysical channels leading to EMRI formation).

In 2013, ESA identified GW detection from space AEA2013 as the science theme to be addressed by the third L-class mission (L3) in its Cosmic Vision scientific programme, with a provisional launch date of 2034. A call for mission proposals to address this science theme was issued in October 2016, and ESA has now selected a three-satellite interferometry mission—the Laser Interferometer Space Antenna (LISA) AEA2017 . EMRI detection formed a key part of both the L3 science case AEA2013 and the specification of the L3 mission requirements AEA2017 . Between a few and a few hundred EMRI observations are expected over the mission lifetime GEA2004 ; G2009 ; BCCG2016 ; BEA2017 , and these have tremendous potential for the purposes of astrophysics, cosmology and fundamental physics.

The observation of as few as ten EMRIs can provide a measurement of the slope of the black-hole mass function to better precision than is currently known GTV2010 . EMRIs can also be used as standard sirens to investigate the expansion history of the Universe, and hence to constrain cosmological parameters MH2008 . Finally, EMRI observations provide exquisite probes of gravitational physics; in particular, they can be used to map out the spacetime structure around the central MBH to high precision, allowing us to test if that geometry is described by general relativity or an alternative theory (see e.g. MCG2017 ). We refer the reader to GVLB2013 for a comprehensive review of testing relativity with EMRIs and other LISA sources.

The scientific objectives described above will only be achieved if EMRIs can be successfully identified and characterised in the noisy LISA data stream. This is difficult because EMRI orbits are generally expected to be both eccentric and inclined to the equatorial plane of the MBH, such that the emitted GWs contain a complex superposition of three fundamental frequencies. Since an EMRI waveform depends on 14 different parameters and $\sim 10^{5}$ waveform cycles can be observed in the LISA frequency band, there are a huge number of independent EMRI signals that must be searched for in the data.

Early work GEA2004 showed in theory that EMRI data analysis is possible using a semi-coherent approach, without actually demonstrating the effectiveness of such a method in practice. Between 2006 and 2011, the feasibility of LISA data analysis in general was explored through a sequence of mock LISA data challenges AEA2006 ; AEA2007a ; BEA2008 ; BEA2008a ; BEA2010 . These showed in practice that it is possible to correctly identify individual EMRIs in data sets without other sources, and using narrow parameter priors BEA2008 ; BEA2008a ; BEA2010 ; BGP2009 . The successful approaches typically used techniques such as Markov chain Monte Carlo to stochastically explore the EMRI parameter space and find the best-fit parameter values.

Much work is still needed to move from these simple initial demonstrations to a practical suite of tools for EMRI characterisation. However, one thing that all of these data analysis techniques have in common is the need for models of EMRI signals that might be present in the data. These models need to be sufficiently faithful to astrophysical EMRI signals in order to identify them in the data, but also computationally inexpensive enough that they can be generated in the large numbers required for grid-based or stochastic searches.

EMRI waveforms can be modelled accurately using black-hole perturbation theory, which exploits the extreme mass ratio to describe the inspiral through expansions in mass ratio. Perturbation theory can be used to construct the gravitational self-force acting on the inspiralling object; these calculations are challenging both theoretically and computationally (see B2009 ; PPV2011 for recent reviews). The self-force at leading order in mass ratio is well understood and has been successfully computed, but tracking the phase of an EMRI accurately over $\sim 10^{5}$ cycles requires second-order self-force calculations, which have yet to be performed (although the theory of such calculations has recently been worked out P2015 ). Additional complicating features such as transient resonances BCCG2016 ; FH2012 also require more work to fully understand.

While fully consistent self-force waveforms for EMRIs will probably be available by the time that LISA is taking data, these are also likely to be extremely computationally intensive. Hence there is a clear need for computationally efficient, but faithful, models to use in data analysis. Two such “kludge” models for generic inspirals exist in the literature; these were constructed to scope out LISA data analysis issues, e.g. EMRI event rates, parameter estimation precision, etc. The “analytic kludge” (AK) model BC2004 uses the quadrupole emission from a Keplerian orbit PM1963 as its base, imposing on top of this relativistic effects such as precession of the periapsis and the orbital plane, and radiation-reaction-driven evolution of the orbital parameters. The AK model is very fast to evaluate, and was used as the reference model in past mock data challenges for this reason. However, it does not provide a good match to true EMRI signals (as gauged by comparison to perturbative calculations) and so is not appropriate for LISA data analysis.

The “numerical kludge” (NK) model builds the trajectory of an EMRI from an exact Kerr geodesic, the parameters of which are then evolved using expressions derived from post-Newtonian (PN) expansions and fits to perturbative calculations GG2006 . A waveform is computed from the resulting trajectory by identifying the Boyer–Lindquist coordinates of the Kerr spacetime with flat-space spherical polar coordinates BEA2007 . NK waveforms show high fidelity with more accurate Teukolsky-based waveforms H2001 ; DH2006 and can readily be extended to include additional physical effects, such as conservative corrections to the evolution HG2009 or a prescription for the force acting on the orbit GEA2011 ; WEA2012 . These waveforms therefore satisfy the requirement of faithfulness to true EMRI signals, but they are more computationally expensive to generate than AK waveforms.

In this paper, we describe an “augmented analytic kludge” (AAK) model that uses information from the NK model to improve the faithfulness of AK waveforms without significantly increasing their computational cost. This is achieved by mapping the parameters of the AK model to match the frequencies of NK waveforms. The AAK model was first introduced in CG2015 ; here we further improve the model, give details on the released implementation at github.com/alvincjk/EMRI_Kludge_Suite, and discuss its implications for EMRI detection with LISA. AAK waveforms are much closer to astrophysical signals than AK waveforms, and so can detect/localise candidate signals in LISA data with sufficient sensitivity/accuracy such that follow-up parameter estimation with perturbative waveforms is feasible. They are also about an order of magnitude faster to evaluate than NK waveforms. Hence the AAK model is well suited for use in future mock data challenges, and possibly even in the final LISA data analysis infrastructure.

A brief overview of the AK and NK models is given in Sec. II, while the AAK model is presented in Sec. III. In Sec. III.1, we first introduce the Kerr fundamental frequencies and the parameter-space map they induce in the AK model, before describing the technical details of the AAK implementation in Sec. III.2. The performance of the augmented model is then compared to that of the original AK model in Sec. III.3, with the more accurate but slower NK model used as the benchmark for both. In Sec. IV, the application of the AAK model to EMRI detection with LISA is considered; we provide an analytic estimate for the threshold signal amplitude required in a semi-coherent search, and assess the viability of the AK and AAK models for real LISA data analysis.

II Kludge waveform models

A kludge in the context of EMRI modelling is any approximate model that uses a combination of formalisms to generate waveforms quickly and extensively for data analysis. Kludge waveforms capture many qualitative features of more accurate EMRI waveforms, and (owing to their modular construction) can be modified to incorporate self-force information as it becomes available.

Two widely used kludges are introduced briefly in this section: the AK waveform of Barack & Cutler BC2004 , which is very fast to compute and provides the basis for our new model, and the NK waveform of Babak et al. BEA2007 , which we take as a fiducial model for calibration and benchmarking purposes. Other approximate EMRI models exist but at varying levels of implementation (e.g. ST2003 ; FEH2016 ), and we do not consider them in this work (apart from the PN fluxes of Sago & Fujita SF2015 , which are used as part of the AAK model).

Assuming the spin of the compact object is negligible, an EMRI can be described by 14 parameters: the two masses $(\mu,M\gg\mu)$ of the system, the three components of the central black hole’s spin vector $\mathbf{S}$, three constants $\mathbf{E}$ describing the compact object’s (instantaneous) orbit, the three components of the compact object’s position vector $\mathbf{X}$ with respect to the black hole, and the three components of the system’s position vector $\mathbf{R}$ with respect to the Solar System.

Of these 14 degrees of freedom, seven are extrinsic to the source: two in $\mathbf{S}$ and one in $\mathbf{X}$ (corresponding to spatial rotation of the source), three in $\mathbf{R}$ (corresponding to spatial translation), and one in $\mathbf{E}$ (corresponding to temporal translation). The parameters of an EMRI model are often chosen to decouple the intrinsic degrees of freedom from the extrinsic ones, which are generally cheaper to search over during data analysis BCV2003 .

Schematically, the main ingredients of a kludge waveform model are then (i) the evolution of the orbital constants along the inspiral (i.e. the “phase-space” trajectory), using PN or fitted fluxes $\mathbf{F}$:

[TABLE]

(ii) the construction of the compact object’s worldline (i.e. the “configuration-space” trajectory), using geodesic or flux-derived expressions $\mathbf{G}$:

[TABLE]

and (iii) the generation of the waveform field $h$ at the detector, using some weak-field multipole formula $H$:

[TABLE]

II.1 Analytic kludge

In the AK model BC2004 , both the orbital trajectory and the waveform are computed in a flat-space approximation, with relativistic effects such as inspiralling and precession added separately. The trajectory is built out of rotating Keplerian ellipses. Radiation reaction is introduced in phase space, where the orbital constants describing a Keplerian ellipse are evolved with PN equations. In configuration space, the orientation of this ellipse is also evolved with PN equations to simulate relativistic precession. The waveform is then generated using the Peters–Mathews mode-sum approximation for Keplerian orbits PM1963 , in which the mass quadrupole moment is decomposed into harmonics of the Keplerian orbital frequency.

Since a Keplerian orbit is confined within the plane normal to its angular momentum vector $\mathbf{L}$, the AK waveform is constructed in an $\mathbf{L}$-based coordinate frame

[TABLE]

and projected transverse to the wave frame

[TABLE]

where the normalisation factor $\mathcal{S}_{\mathbf{L},\mathbf{R}}:=(1-(\hat{\mathbf{L}}\cdot\hat{\mathbf{R}})^{2})^{1/2}$. These two frames are made time-varying (with respect to a fixed heliocentric and ecliptic-based frame C1998 ) through the forced precession of $\mathbf{L}$.

The two waveform polarisations in the transverse–traceless gauge (with the usual $(+,\times)$ convention for $(\hat{\mathbf{x}},\hat{\mathbf{y}})_{\mathrm{AK}}$) are given by the $n$-mode sums

[TABLE]

with

[TABLE]

[TABLE]

where $\tilde{\gamma}$ is an azimuthal angle in the orbital plane measuring the direction of periapsis with respect to $(\hat{\mathbf{R}}\cdot\hat{\mathbf{L}})\hat{\mathbf{L}}-\hat{\mathbf{R}}$ (i.e. the orthogonal projection of $\hat{\mathbf{z}}_{\mathrm{AK}}$ onto the plane normal to $\hat{\mathbf{L}}$).111We have changed some of the notation in BC2004 for consistency, since we are constructing a hybrid model using different formalisms. For example: the notation for the angles $\gamma$ and $\tilde{\gamma}$ has been swapped; the notation $\nu$ for the orbital frequency $\dot{\Phi}/(2\pi)$ is unused; the notation for the inclination $\lambda$ is now $\iota$. The functions $(a_{n},b_{n},c_{n})$ describe the changing mass quadrupole moment of a Keplerian orbit with mean anomaly $\Phi(t)$, eccentricity $e$ and orbital angular frequency $\dot{\Phi}$, and are given by PM1963

[TABLE]

[TABLE]

[TABLE]

where the $J_{n}$ are Bessel functions of the first kind, and $\mathcal{A}=(\dot{\Phi}M)^{2/3}\mu/|\mathbf{R}|$ in the extreme-mass-ratio limit.

In the ecliptic-based coordinate system, the sky position $\hat{\mathbf{R}}\equiv(\theta_{S},\phi_{S})$ of the source and the black-hole spin orientation $\hat{\mathbf{S}}\equiv(\theta_{K},\phi_{K})$ are effectively constant. It is convenient to represent $\hat{\mathbf{L}}$ in ecliptic coordinates with respect to $\hat{\mathbf{S}}$. We have

[TABLE]

where $\hat{\mathbf{z}}=[0,0,1]^{T}$ is normal to the ecliptic plane, $\iota$ is the inclination angle between $\hat{\mathbf{L}}$ and $\hat{\mathbf{S}}$, and $\alpha$ is an azimuthal angle in the spin-equatorial plane measuring the direction of $\hat{\mathbf{L}}-(\hat{\mathbf{L}}\cdot\hat{\mathbf{S}})\hat{\mathbf{S}}$ with respect to $\hat{\mathbf{z}}-(\hat{\mathbf{z}}\cdot\hat{\mathbf{S}})\hat{\mathbf{S}}$ (i.e. the angle between the orthogonal projections of $\hat{\mathbf{L}}$ and $\hat{\mathbf{z}}$ onto the plane normal to $\hat{\mathbf{S}}$). Furthermore, since $\tilde{\gamma}$ is neither intrinsic nor extrinsic, it is useful to define the purely intrinsic parameter $\gamma:=\tilde{\gamma}-\beta$, where $\beta=\beta(\hat{\mathbf{R}},\hat{\mathbf{S}},\hat{\mathbf{L}})=\beta(\theta_{S},\phi_{S},\theta_{K},\phi_{K},\iota,\alpha)$ is an azimuthal angle in the orbital plane measuring the direction of $\hat{\mathbf{L}}\times\hat{\mathbf{S}}$ with respect to $(\hat{\mathbf{R}}\cdot\hat{\mathbf{L}})\hat{\mathbf{L}}-\hat{\mathbf{R}}$.

Only one of the six parameters comprising $\mathbf{E}=(e,\iota,\dot{\Phi})$ and $\mathbf{X}=(\Phi(t),\gamma,\alpha)$ in the above Keplerian setup changes with time. In the AK model, $(e,\dot{\Phi},\gamma,\alpha)$ are promoted to functions of time and evolved with mixed-order PN expressions that depend on $(\mu,M,a=|\mathbf{S}|/M,\mathbf{E})$ BC1975 ; B1991 ; JS1992 ; R1996 , while $\iota$ is approximated as constant (since the inclination angle of a typical EMRI varies extremely slowly H2000 ). The Keplerian orbit shrinks and circularises as $\dot{\Phi}(t)$ and $e(t)$ increase and decrease respectively. From (12), the time dependence of the orbital orientation $\hat{\mathbf{L}}(t)$ is confined to $\alpha(t)$, where $\dot{\alpha}$ is precisely the angular rate of Lense–Thirring precession. Finally, the angular rate of periapsis precession is given by $\dot{\gamma}+\dot{\alpha}$ since $\gamma(t)$ is measured with respect to $\hat{\mathbf{L}}(t)\times\hat{\mathbf{S}}$.

While the waveform field is effectively planar at the Solar System and may be calculated in the fixed heliocentric frame (as opposed to a detector-centric one), the rotational and orbital motion of LISA in the ecliptic plane must be factored into the detector’s response to the GW. In the standard LISA framework, the waveform polarisations $h_{+,\times}$ are transformed into the response functions $h_{I,II}$ via

[TABLE]

where the antenna pattern functions ACST1994

[TABLE]

[TABLE]

[TABLE]

[TABLE]

depend on the sky location $(\theta_{D},\phi_{D})$ and polarisation angle $\psi_{D}$ of the source in a detector-based coordinate system, and hence rotate with respect to $\hat{\mathbf{R}}$ as the plane of the detector along its orbit precesses around the ecliptic plane. Doppler modulation of the waveform phase (through $\Phi(t)$) is also included to correct for the orbital motion of the detector itself.

The AK model was the first waveform model used to investigate the precision of LISA parameter estimation over the full (modulo compact-object spin) EMRI parameter space BC2004 . Due to its computational efficiency, the model has also been employed in past mock LISA data challenges to generate injected signals in simulated data and parametrised templates for search algorithms BEA2008 ; BEA2008a ; BEA2010 ; BGP2009 . However, the approximate waveforms it produces are demonstrably inaccurate, and will result in reduced detection and parameter estimation performance if used to analyse data sets containing realistic EMRI signals.

II.2 Numerical kludge

In the NK model BEA2007 , the orbital trajectory is computed in curved space with a treatment that is fully relativistic up to the evolution of orbital constants GHK2002 ; GG2006 , i.e. it is built out of Kerr geodesics. The three constants of motion for a geodesic are evolved with Teukolsky-fitted PN equations, which introduces radiation reaction. In configuration space, precession effects are obtained for free by integrating the geodesic equations along the phase-space trajectory. The curved-space coordinates of the compact object’s worldline are then associated artificially with coordinates in flat space, and the waveform is generated using the standard quadrupole formula (or variants that include additional contributions from higher-order moments of mass B1973 ; P1977 ).

The NK waveform is constructed in an $\mathbf{S}$-based coordinate frame

[TABLE]

and projected transverse to the wave frame

[TABLE]

where the normalisation factor $\mathcal{S}_{\mathbf{S},\mathbf{R}}:=(1-(\hat{\mathbf{S}}\cdot\hat{\mathbf{R}})^{2})^{1/2}$. Aligning the $z$-axis with $\mathbf{S}$ is a more natural choice for the NK model, since the compact object’s worldline is computed in Boyer–Lindquist coordinates. The two wave frames (5) and (19) are related by a (time-varying) rotation about $\mathbf{R}$.

Using the standard quadrupole formalism, the two waveform polarisations in the transverse–traceless gauge (with the $(+,\times)$ convention for $(\hat{\mathbf{x}},\hat{\mathbf{y}})_{\mathrm{NK}}$) are given by

[TABLE]

with

[TABLE]

where $\ddot{I}_{ij}(t)$ is the second time derivative of the source’s mass quadrupole moment $I_{ij}(t)$. The polarisation tensors $H^{+,\times}_{ij}$ and transverse projection tensor $P_{ij}$ are given by

[TABLE]

[TABLE]

[TABLE]

where $\delta_{ij}$ is the Kronecker delta.

In the extreme-mass-ratio limit, the mass quadrupole moment is simply

[TABLE]

where the $x_{i}(t)$ are Cartesian components of the compact object’s position vector $\mathbf{X}$ with respect to the frame (18) centred on the black hole. Although (21) (with (25)) is a weak-field equation in flat-space coordinates, the NK model specifies and calculates $(x_{1},x_{2},x_{3})=(r\sin{\theta}\cos{\phi},r\sin{\theta}\sin{\phi},r\cos{\theta})$ in Boyer–Lindquist coordinates. The self-consistency of this approach clearly degrades further into the strong field, but does not severely impact the effectiveness of the NK waveforms as an approximation to Teukolsky-based ones BEA2007 .

A timelike Kerr geodesic is described fully by three first integrals of motion: the orbital energy $E$, the projection $L_{z}$ of the orbital angular momentum $\mathbf{L}$ onto $\mathbf{S}$, and the quadratic Carter constant $Q$. Along such an orbit, $(r(t),\theta(t),\phi(t))$ are obtained by integrating the geodesic equations for a test particle in the Kerr spacetime; these are written in canonical form as MTW1973

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\tau$ is proper time along the worldline and $\Sigma=r^{2}+a^{2}\cos^{2}\theta$. The potential functions $V_{r,\theta,\phi,t}$ are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with $P=E(r^{2}+a^{2})-aL_{z}$ and $\Delta=r^{2}-2Mr+a^{2}$.

In practice, it is convenient to work with alternative parametrisations of $(E,L_{z},Q)$. For a bound orbit, the geodesic may be specified by the parameters $(r_{p},r_{a},\theta_{\mathrm{min}})$ (the values of $r$ at periapsis and apoapsis, and the minimal value of $\theta$ respectively), which fully describe the range of motion in the radial and polar coordinates. The roots of $V_{r}$ determine $r_{p}$ and $r_{a}$, while the roots of $V_{\theta}$ determine $\cos{\theta_{\mathrm{min}}}$ (the maximal value of $\cos{\theta}$). Another parametrisation is $(e,\iota,p)$ (the quasi-Keplerian eccentricity, inclination and semi-latus rectum); these are defined in terms of $(r_{p},r_{a},\theta_{\mathrm{min}})$ as

[TABLE]

Finally, since the configuration-space parameters $(r,\theta)$ oscillate between the bounds $r_{p}\leq r\leq r_{a}$ and $\theta_{\mathrm{min}}\leq\theta\leq\pi-\theta_{\mathrm{min}}$, it is useful to define

[TABLE]

such that $\psi$ (the quasi-Keplerian true anomaly) and $\chi$ are the phases of radial and polar motion respectively.

The orbital constants $\mathbf{E}=(E,L_{z},Q)$ in the above geodesic setup do not vary with time. Radiation reaction is added to the NK model by evolving $\mathbf{E}$ with fluxes that depend on $(\mu,M,a,\mathbf{E})$ (note that the inclination $\iota(t)=\tan^{-1}(\sqrt{Q}/L_{z})$ is correctly time-dependent in this model). These fluxes are mixed-order PN expressions that have been fitted to the results of Teukolsky-based computations for circular inclined orbits GG2006 . Integrating the geodesic equations along the phase-space trajectory then gives $\mathbf{X}=(\psi(t),\chi(t),\phi(t))$, complete with relativistic precession. Once the waveform polarisations $h_{+,\times}$ have been calculated via (20)–(25), the LISA response functions $h_{I,II}$ may be obtained through the method outlined in Sec. II.1.

Waveforms from the NK model display excellent agreement with Teukolsky-based waveforms in the strong-field regime; they are reliable up to a closest approach of $r_{p}\approx 5M$, with typical matches of over 0.95 BEA2007 . NK waveforms might even be accurate enough to serve as templates in actual LISA detection algorithms. However, they are still slightly expensive to generate in large numbers due to the relatively elaborate construction of the phase- and configuration-space trajectories, while added computational cost also arises in the parameter conversion $(E,L_{z},Q)\leftrightarrow(e,\iota,p)$, the handling of plunge, etc.

III Augmented analytic kludge

The AK model is 5–15 times faster than the NK model at generating year-long waveforms sampled at 0.2 Hz for a generic $(10^{1},10^{6})M_{\odot}$ EMRI with low initial eccentricity ($e_{0}\lesssim 0.3$); this speed-up is enhanced for longer waveform durations, but diminished for higher initial eccentricity (since more modes must be summed in the Peters–Mathews approximation).222The sums in (6) must be truncated at some arbitrary number of modes $N$, which directly affects both the speed and accuracy of the AK model. This number may be specified by setting a threshold for the relative power radiated into the $N$-th harmonic, and has been experimentally determined to scale linearly with eccentricity BC2004 . We use $N=\lfloor 30e_{0}\rfloor$ as the default value for both the AK and AAK models. However, AK waveforms suffer from severe dephasing with respect to NK waveforms, even at the early-inspiral stage. In Fig. 1, the AK waveform for a $(10^{1},10^{6})M_{\odot}$ EMRI with initial semi-latus rectum $p_{0}=15M$ matches the qualitative features of the corresponding NK waveform, but is a full cycle out of phase within three hours. This is due to the mismatched frequencies in the two models.

In Secs III.1 and III.2, we describe the construction of a hybrid model that capitalises on the benefits of both kludges. The AK model is augmented with an initial map to the fundamental frequencies of Kerr geodesic motion, which corrects the instantaneous phasing as shown in Fig. 1. Over longer timescales, the mapped orbital trajectory is further improved through self-consistent PN evolution and a local polynomial fit to the phase-space trajectory of the NK model. Fast algorithms for higher-order fits and plunge handling have been incorporated in the latest implementation of the AAK model, which has been made publicly available at github.com/alvincjk/EMRI_Kludge_Suite as part of a software suite for generating kludges.

The initial version of the AAK model yields waveforms that can remain phase-coherent with NK waveforms for over two months, but with overlap values lower than 0.97 CG2015 . This is the commonly chosen minimal match for a waveform template bank that corresponds to a $90\%$-ideal observed event rate O1996 , and thus ensures the equivalent localisation of any signal detected with such banks of AAK and NK templates. In Sec. III.3, we report further improved results for the present AAK implementation. Two-month overlaps higher than 0.97 are found for EMRIs with varying spin and eccentricity; however, the overlaps still degrade with proximity to plunge, due to the divergence of the AAK and NK trajectories deep within the strong field.

III.1 The fundamental-frequency map

The geodesic equations (26)–(29) take a simple form with the choice of a timelike parameter $\lambda=\int d\tau/\Sigma$ C1983 ; M2003 ; this decouples (26) and (27), and for a bound orbit makes the radial and polar components of motion manifestly periodic with respect to $\lambda$. For the azimuthal and temporal components (whose potentials depend only on $(r,\theta)$), overall rates of evolution may be obtained by averaging (32) and (33) over many periods of radial and polar motion.

From the radial and polar periods $\Lambda_{r,\theta}$, the average azimuthal rate $\langle d\phi/d\lambda\rangle$ and the average temporal rate $\langle dt/d\lambda\rangle$ (denoted $\Gamma$ by analogy with the Lorentz factor), we may define three angular and dimensionless fundamental frequencies $\Omega_{r,\theta,\phi}$ for the test particle’s motion with respect to coordinate time. In terms of $(r_{p},r_{a},\theta_{\mathrm{min}})$, these frequencies are written as DH2004 ; SF2015

[TABLE]

[TABLE]

where $\Lambda_{r,\theta}$ and $\Gamma$ are given by

[TABLE]

[TABLE]

Expressions for $\Omega_{r,\theta,\phi}$ in terms of $(e,\iota,p)$ have been derived by Schmidt S2002 ; these are less compact, but have more utility in practical implementations.

In terms of the fundamental frequencies, the periapsis and Lense–Thirring precession rates are given by $\Omega_{\phi}-\Omega_{r}$ and $\Omega_{\phi}-\Omega_{\theta}$ respectively. These vanish in the Newtonian limit, where $\Omega_{r,\theta,\phi}$ approach a single orbital frequency $\Omega$ from below, i.e. $\Omega_{r}\nearrow\Omega_{\theta}\nearrow\Omega_{\phi}\nearrow\Omega$. The frequency $\Omega$ is then related to $(e,p)$ by Kepler’s third law:

[TABLE]

In the AK model, $\Omega$ is associated with the quantity $\dot{\Phi}M$; however, periapsis and Lense–Thirring precession are added on top of $\dot{\Phi}$ via $\dot{\gamma}+\dot{\alpha}$ and $\dot{\alpha}$ respectively, and so $\Omega$ is the lowest frequency by construction. This inconsistency with the relativistic case leads to mismatched frequencies when supplying identical parameters to the two models, since the same value of $p$ in (40) specifies the radial AK frequency while approximating the azimuthal NK frequency. In other words, the frequencies in the AK model are generally too high.

A three-dimensional endomorphism over the AK space of orbits is induced by requiring that the radial, polar and azimuthal frequencies $(\dot{\Phi},\dot{\Phi}+\dot{\gamma},\dot{\Phi}+\dot{\gamma}+\dot{\alpha})$ for any $(e,\iota,p)$ have the same values as the (dimensionful) relativistic frequencies $\omega_{r,\theta,\phi}:=\Omega_{r,\theta,\phi}/M$. We map the parameters $(M,a,p)$ rather than $(e,\iota,p)$ to unphysical values; this gives better results since periapsis and Lense–Thirring precession are more directly determined by the central mass and its rotation respectively. The map $(M,a,p)\mapsto(\tilde{M},\tilde{a},\tilde{p})$ is given implicitly by solving the algebraic system of equations

[TABLE]

[TABLE]

[TABLE]

for the unphysical set $(\tilde{M},\tilde{a},\tilde{p})$, which is defined as the root closest to the physical set $(M,a,p)$ with a Euclidean metric on parameter space.

Substituting $(\tilde{M},\tilde{a},\tilde{p})$ for $(M,a,p)$ in the AK model provides an instantaneous correction of its frequencies at any point along the inspiral trajectory $(e(t),\iota(t),p(t))$. In principle, applying the map along the entirety of a fiducial inspiral will keep the AK waveform phase-coherent with relativistic waveforms (generated from that trajectory) until plunge. However, such an inspiral is usually more expensive to compute (as in the case of the NK model), and the additional cost from evaluating the map itself also scales linearly with the number of points sampled along the trajectory. Complications also arise as the compact object approaches the point of plunge, where the fundamental frequencies diverge and the map (41)–(43) is no longer well-defined.

III.2 Implementation

In order to retain the main advantage of the AK model, computational costs are kept as low as possible by evaluating the map at a small number of points and relying on independent evolution of the orbital constants over long timescales. Firstly, the NK phase-space trajectory is generated at and around the specified initial time $t_{0}$ over a user-defined timescale $T_{\mathrm{fit}}$, which depends on the radiation-reaction timescale $T_{\mathrm{RR}}:=M^{2}/\mu$ and specifies the duration over which the AAK model is calibrated. The timescale $T_{\mathrm{fit}}$ and number of sample points $N_{\mathrm{fit}}$ may be adjusted adaptively based on the proximity of the initial point $(e_{0},\iota_{0},p_{0}):=(e(t_{0}),\iota(t_{0}),p(t_{0}))$ to plunge; they typically satisfy $0.1T_{\mathrm{RR}}\lesssim T_{\mathrm{fit}}\lesssim 10T_{\mathrm{RR}}$ (over which $\iota$ is approximately constant) and $N_{\mathrm{fit}}\lesssim 10$ (to ensure an added computational cost of $\lesssim 1\%$). Evaluation of the map at each of the $N_{\mathrm{fit}}$ points gives a local “best-fit” trajectory $(\tilde{M}(t),\tilde{a}(t),e(t),\tilde{p}(t))_{\mathrm{fit}}$ in the AAK phase space, where the unphysical $(\tilde{M}(t),\tilde{a}(t))_{\mathrm{fit}}$ change with time.

The global trajectory in the AAK phase space is generated independently from the NK model, and hence more rapidly. From the mapped initial point $(\tilde{M}_{0},\tilde{a}_{0},e_{0},\tilde{p}_{0})$ (which also lies on the best-fit trajectory by construction), a global PN trajectory $(\tilde{M},\tilde{a},e(t),\tilde{p}(t))_{\mathrm{PN}}$ is obtained by evolving $(e(t),\tilde{p}(t))_{\mathrm{PN}}$ with 3PN $O(e^{6})$ expressions given by Sago & Fujita SF2015 , while $(\tilde{M},\tilde{a})_{\mathrm{PN}}$ are left constant. Higher-order 4PN $O(e^{6})$ expressions SF2015 have also been tested, but these seem to result in poorer agreement with NK waveforms (which use fluxes of up to 3PN), possibly due to the known divergence of certain expansions beyond 3PN order YB2008 .

Sampling the PN trajectory at each of the $N_{\mathrm{fit}}$ points allows the calculation of $(\tilde{M}(t),\tilde{a}(t),e(t),\tilde{p}(t))_{\mathrm{fit}}-(\tilde{M},\tilde{a},e(t),\tilde{p}(t))_{\mathrm{PN}}$ over the duration $T_{\mathrm{fit}}$. This difference trajectory is fitted to polynomials in time and extrapolated over the lifetime of the inspiral; it is then added to the PN trajectory, giving the final global trajectory $(\tilde{M}(t),\tilde{a}(t),e(t),\tilde{p}(t))$. In the initial AAK implementation, the coefficients of the quadratic fit are given by second-order finite-difference quotients (i.e. $N_{\mathrm{fit}}=3$), which works well but only for small values of $T_{\mathrm{fit}}$. The present version uses a quartic least-squares fit, which allows the choice of longer $T_{\mathrm{fit}}$ and consequently gives better long-term phase agreement with NK waveforms.

After the phase-space trajectory $(\tilde{M}(t),\tilde{a}(t),e(t),\tilde{p}(t))$ has been constructed, the configuration-space evolution of $(\Phi,\gamma,\alpha)$ is performed with the appropriate combinations (41)–(43) of the fundamental frequencies, given by Sago–Fujita expressions SF2015 that are consistent at 3PN $O(e^{6})$ with the phase-space evolution. The waveform field is then generated as in the AK model. For illustrative purposes, a flowchart summary of the entire AAK algorithm (as presently implemented) is shown in Fig. 2.

The augmentations to the AK framework are focused on improving the phase information of its waveforms, since the amplitude of a GW signal is measured far less precisely than its phase. However, the calculation of amplitude in the original AK model ($\mathcal{A}$ in (9)–(11)) is a decent approximation since it is based on $\dot{\Phi}M$, which is assigned a value $\approx\Omega_{\phi}$ that turns out to be correct for this purpose (see discussion around (40)). Hence it is the AAK amplitude that is shifted away from the fiducial NK value through the mapping of frequencies and the unphysical evolution of $\tilde{M}$. A simple adjustment is made to reverse this shift; in (9)–(11) for the AAK model, the amplitude is now given by

[TABLE]

where the azimuthal frequency $\omega_{\phi}$ and the physical black-hole mass $M$ are used in place of $\dot{\Phi}=\omega_{r}$ and $\tilde{M}$ (i.e. $M$ in the original AK model) respectively.

Finally, a fast method of plunge handling has been added to the present AAK implementation; this feature is useful in general, but especially when generating large numbers of AAK waveform templates for search algorithms. The compact object plunges when its instantaneous orbit along the phase-space trajectory $\mathbf{E}(t)$ becomes unstable, i.e.

[TABLE]

where $V_{r}$ is given in (30) with $\mathbf{E}=(E,L_{z},Q)$. This point is termed the last stable orbit $\mathbf{E}_{\mathrm{LSO}}$, and is precisely the point at which the discriminant $\mathcal{D}(a,\mathbf{E})$ of the quartic polynomial $V_{r}(r)$ changes sign from positive (four real roots) to negative (two real roots) R1922 . Since $\mathcal{D}$ is a simple analytic function of the quartic coefficients, it is computationally trivial to check for stability at every integration step for the phase-space trajectory, provided the evolution is done in terms of $(E,L_{z},Q)$.333The discriminant method is applicable to the NK model, and may speed it up slightly. Currently, the NK implementation precomputes (for the specified value of $a$) an interpolated $p_{\mathrm{LSO}}$ surface over the relevant region of $(e,\iota)_{\mathrm{LSO}}$ space, by finding and examining the roots of $V_{r}(r)$ numerically. It then checks for $p<p_{\mathrm{LSO}}(e,\iota)$ when generating the phase-space trajectory.

Plunge detection in the AAK model is far less straightforward than in the other two kludges, partly because the evolution is performed in terms of the quasi-Keplerian orbital parameters, and the computational benefits of the discriminant method are nullified by having to convert $(e,\iota,\tilde{p})\to(E,L_{z},Q)$. Furthermore, $(\tilde{M},\tilde{a},\tilde{p})$ are unphysical; the inverse of the map (41)–(43) is computationally expensive and (more crucially) ill-defined at plunge, and so cannot be used to obtain the physical parameters for stability calculations.

To circumvent these issues, the AAK model uses (40) with $\Omega\approx\Omega_{\phi}$ to obtain an approximation for the physical parameter $p$. While generating the phase-space trajectory, it checks (between the least- and most-bound orbits H2001 ) the stability of $(e,\iota,p)$ at every radiation-reaction interval $T_{\mathrm{RR}}$. Once the stability changes across an interval, it then bisects that interval to find $p_{\mathrm{LSO}}$, and smoothly zeroes the waveform over ten additional orbits with a one-sided Planck-taper window MRS2010 . The added computational cost associated with this algorithm is $\lesssim 1\%$. Although the approximation for $p$ is crude, the phase-space trajectories in the AAK and NK models are generally divergent to begin with, and the plunge points for both models may differ significantly even if a more accurate expression is used.

III.3 Benchmarking results

The specified initial state of an EMRI in the AAK model is described by the intrinsic parameters $(\mu,M,a,e_{0},\iota_{0},\gamma_{0},\psi_{0})$ and the extrinsic parameters $(p_{0},\theta_{S},\phi_{S},\theta_{K},\phi_{K},\alpha_{0},D)$, where $D:=|\mathbf{R}|$.444For a source at cosmological redshift $z$, the values of $D$ and $(\mu,M)$ are replaced with the luminosity distance $D(1+z)$ and the redshifted masses $(\mu(1+z),M(1+z))$ respectively. In configuration space, transformations from $\mathbf{X}_{\mathrm{AAK}}=(\psi,\gamma,\alpha)$ to $\mathbf{X}_{\mathrm{AK}}=(\Phi,\gamma,\alpha)$ and $\mathbf{X}_{\mathrm{NK}}=(\psi,\chi,\phi)$ are required for a comparison of the three waveform models. We have chosen to specify the quasi-Keplerian true anomaly $\psi$ in the (shared) AAK parameter space, since there is no closed-form expression for $\psi$ in terms of the mean anomaly $\Phi$. The conversion $\psi\to\Phi$ is given by the Keplerian expressions C2005

[TABLE]

where $E$ is known as the eccentric anomaly.

On the other hand, the AAK model retains the AK parameters $(\gamma,\alpha)$; these have explicit meanings in the (intrinsic) $\mathbf{L}$-based coordinate frame

[TABLE]

where the normalisation factor $\mathcal{S}_{\mathbf{L},\mathbf{S}}:=(1-(\hat{\mathbf{L}}\cdot\hat{\mathbf{S}})^{2})^{1/2}$ and $\hat{\mathbf{L}}(\alpha)$ is given in ecliptic coordinates by (12). The unit position vector of the compact object with respect to (47) is $\hat{\mathbf{r}}_{\mathbf{L}}=[\cos{(\psi+\gamma)},\sin{(\psi+\gamma)},0]^{T}$, and a change of basis to the $\mathbf{S}$-based coordinate frame (18) gives

[TABLE]

where the orthogonal matrices $\mathbf{Q}_{\mathbf{S}}:=[\hat{\mathbf{x}}|\hat{\mathbf{y}}|\hat{\mathbf{z}}]_{\mathbf{S}}$ and $\mathbf{Q}_{\mathbf{L}}:=[\hat{\mathbf{x}}|\hat{\mathbf{y}}|\hat{\mathbf{z}}]_{\mathbf{L}}$ are formed from the triads in (18) and (47) respectively. It is then straightforward to obtain $(\chi,\phi)$ from $\hat{\mathbf{r}}_{\mathbf{S}}=[\sin{\theta}\cos{\phi},\sin{\theta}\sin{\phi},\cos{\theta}]^{T}$, via (35).

As a generic example, we consider a prograde EMRI with redshifted component masses $(\mu,M)=(10^{1},10^{6})M_{\odot}$, spin $a=0.5M$, and initial orbital parameters $(e_{0},\iota_{0},p_{0})=(0.1,\pi/6,8.25M)$. The initial semi-latus rectum is chosen such that the compact object plunges approximately one year after entering the LISA band at a representative frequency $f_{\mathrm{GW}}=2.7\,\mathrm{mHz}$, where $f_{\mathrm{GW}}$ is defined as twice the azimuthal orbital frequency (i.e. the dominant GW harmonic at low eccentricity). In the AAK model, the fitting timescale and number of sample points are set to $T_{\mathrm{fit}}\leq 10T_{\mathrm{RR}}$ and $N_{\mathrm{fit}}\leq 10$ respectively, with inequality in the case of adaptive adjustments.

One important result from our comparison studies is that the AK model can lead to an overestimation of SNR if used without modification. This is due to the fact that the frequencies in the AK model are generally too high for any given $(e,\iota,p)$, as mentioned in Sec. II.1. The SNR of a signal $h=h_{I}+ih_{II}$ is defined as $\rho:=\sqrt{\langle h|h\rangle}$, where $\langle\cdot|\cdot\rangle$ is the standard matched-filtering inner product (on the space of finite-length time series) between two waveforms, i.e. CF1994

[TABLE]

Throughout this paper, the noise power spectral density $S_{n}(f)$ is taken to be a LISA noise model for the L6A5 (six links, five-million-kilometre arms) configuration known as classic LISA, assuming the original mission requirements and including confusion noise from the foreground of Galactic white-dwarf binaries KEA2016 .

At a luminosity distance of 5 Gpc, an NK signal (sampled at 0.2 Hz) from the example EMRI described above has a one-year SNR of $\rho=30.8$. When using the AAK model to generate the signal, a comparable value of $\rho=32.1$ is obtained. However, an AK signal from the same EMRI has $\rho=57.8$. To illustrate this, we consider the characteristic strain $h_{c}$ of a signal and the noise amplitude $h_{n}$; these are given respectively by MCB2015

[TABLE]

such that

[TABLE]

With these definitions, the area between $h_{c}$ and $h_{f}$ on a log–log plot gives an indication (but not an approximation) of SNR, and allows the relative detectability of signals to be estimated. The characteristic strain for the three signals and the LISA noise amplitude are shown in Fig. 3, where the excess power in the AK signal at higher frequencies is evident, along with the consequent boost to SNR. This error is likely to persist for $M\gtrsim 10^{6}M_{\odot}$, but may be mitigated for less massive central black holes as the maxima of the three $h_{c}$ curves are blueshifted past the minimum of $h_{n}$.

For the purposes of this paper (where the NK model is taken as fiducial), the phase accuracy of the AK and AAK models is assessed by how well their waveforms overlap with NK waveforms. The overlap $\mathcal{O}$ between two waveforms $a$ and $b$ is defined as

[TABLE]

which takes the value of one for identical waveforms and zero for orthogonal waveforms. In CG2015 , the overlaps $\mathcal{O}(h_{\mathrm{AK}}|h_{\mathrm{NK}})$ and $\mathcal{O}(h_{\mathrm{AAK}}|h_{\mathrm{NK}})$ over two and six months are computed for the example EMRI, as well as for the same source with (i) $\mu=10^{0}M_{\odot}$, (ii) $a=0.8M$ and (iii) $e_{0}=0.5$. The AK and initial AAK models have virtually identical computation times $\tau$, and are both quicker than the NK model with typical speed-up factors of $\sigma:=1-\tau/\tau_{\mathrm{NK}}\approx 0.9$ (except in the case of $e_{0}=0.5$, where $\sigma\approx 0.4$). However, the AAK model yields overlaps that are consistently higher, and by 2–3 orders of magnitude in most cases.

The speed-up factors for the AAK model are preserved by the present implementation, while its overlaps are increased across the board due to the enhanced fitting algorithm. Most notably, there is substantial improvement for EMRIs with higher initial eccentricities, although this is partly attributable to the use of a higher sample rate than that in CG2015 . The overlap and timing performance of the AAK model across 2–6 months with varying compact-object mass, black-hole spin and initial eccentricity is summarised by the plots in Figs 4–6.

In Fig. 4, the overlap $\mathcal{O}(h_{\mathrm{AAK}}|h_{\mathrm{NK}})$ is computed for the example EMRI, as well as for the same source with $0.5\leq\lg{(\mu/M_{\odot})}\leq 2$. The minimum value of $\mu\approx 3M_{\odot}$ is chosen since the upper bound of $10T_{\mathrm{RR}}$ for $T_{\mathrm{fit}}$ exceeds six months for a less massive compact object, and so the overlaps in that regime will not show significant improvement (even the six-month overlap is already $>0.97$). We also do not consider intermediate-mass black holes with $\mu>100M_{\odot}$. Instead of choosing $p_{0}$ such that the EMRI plunges after one year (as done in CG2015 ), we consider fixed $p_{0}=8.25M$ in this analysis; this leads the overlaps to degrade at larger rather than smaller mass ratios $\mu/M$. The two-month overlaps are $>0.97$ up to $\mu\approx 20M_{\odot}$, for which plunge occurs at around 5.6 months. There is greater speed-up over the NK model for longer waveform durations as expected, but all values of $\sigma$ are $\gtrsim 0.8$.

Overlaps for the example EMRI with varying spin $0.1\leq a/M\leq 0.9$ are shown in Fig. 5. Again, all values of $\sigma$ are $\gtrsim 0.8$, with greater speed-up for longer waveform durations. For fixed $p_{0}$, prograde EMRIs with lower spin start closer to plunge, and so the overlap values generally increase along with $a$. There appears to be an opposing effect at higher spin (possibly due to the additional degrees of freedom for error from fitting $\tilde{M}$ and $\tilde{a}$ in the AAK model) that causes a fall-off in the four- and six-month overlaps for $a\gtrsim 0.6$. The two-month overlaps are $>0.97$ across the full range of considered spins, which is perhaps unsurprising since the upper bound of $10T_{\mathrm{RR}}$ for $T_{\mathrm{fit}}$ is also two months for a $(10^{1},10^{6})M_{\odot}$ EMRI. However, we note here that while the phase accuracy of the AAK model may be arbitrarily increased in principle by taking $T_{\mathrm{fit}}>10T_{\mathrm{RR}}$, the computational requirement that $N_{\mathrm{fit}}\lesssim 10$ will likely reduce the quality of the trajectory fit at early times.

The effect of varying eccentricity for the example EMRI is illustrated in Fig. 6. We consider initial eccentricities $0.05\leq e_{0}\leq 0.5$, since there turns out to be no speed-up over the NK model ($\sigma\approx 0$) when generating a two-month AAK waveform with $e_{0}=0.5$. This is an important limitation of the Peters–Mathews approximation for waveform generation (see discussion in footnote 2), and will have to be addressed if the AAK model is to be useful in searches for high-eccentricity EMRIs. Nevertheless, the two- and even four-month overlaps are $>0.97$ for $e_{0}\lesssim 0.3$, with $\sigma\gtrsim 0.5$. As in the case of Fig. 5, there is a peak in the six-month overlap; this is probably due to variance in the fit for $e$, and is unlikely to carry any fundamental significance.

Finally, although parameter estimation with AAK waveforms is not the focus of this paper, we summarise here the results of a Fisher matrix calculation for a year-long AAK signal from the example EMRI considered throughout Sec. III.3. The Fisher information matrix $\boldsymbol{\Gamma}$ for a GW signal $h$ parametrised by $\boldsymbol{\lambda}$ is given by CF1994

[TABLE]

where $\boldsymbol{\lambda}=(\mu,M,\mathbf{S},\mathbf{E},\mathbf{X},\mathbf{R})$ for EMRIs. The parameter-estimation errors $\Delta\boldsymbol{\lambda}$ due to Gaussian noise have the normal distribution $\mathcal{N}(\mathbf{0},\boldsymbol{\Gamma}^{-1})$ in the case of high SNR, and so the root-mean-square errors in the general case can be approximated as

[TABLE]

For our AAK signal normalised to an SNR of $\rho=30$, we find that the masses and spin can be measured to within the fractional errors

[TABLE]

[TABLE]

[TABLE]

These errors are roughly an order of magnitude better than the corresponding values for the AK model BC2004 , and are more comparable to those cited for the NK model HG2009 (although the latter consider a circular, equatorial EMRI with $a=0.9M$). They are also consistent with (the lower end of) the values reported in the L2/L3 mission proposal for ESA’s Cosmic Vision programme AEA2013 , where the AK model was used but with a modified plunge criterion.

IV Data analysis application: Semi-coherent detection searches

GW detection is usually achieved via a template bank search, in which a large set of signal templates $h_{i}$ is compared against the noisy detector data $s$. A single coherent integral of template and data is calculated for each template in the bank and used as a detection statistic, i.e. a detection is claimed if any of these values exceeds a predetermined threshold. For EMRI detection, such a procedure is hampered by the extremely large number of templates needed to cover the parameter space; a previous estimate GEA2004 put this number at $N_{\mathrm{bank}}\sim 10^{40}$. In that same work, a computationally viable alternative was suggested: a semi-coherent search that involves splitting the time series data into $N$ segments, searching each segment separately with a smaller template bank, then combining the results to obtain a new detection statistic. Inevitably, such a search is less sensitive than the (computationally prohibitive) fully coherent search.

The minimum feasible number of segments for a semi-coherent EMRI search was estimated in GEA2004 by considering the computational resources anticipated to be available when LISA flies; it was found that $N\gtrsim 100$ for a mission lasting $\sim 10^{8}\,\mathrm{s}$, giving segments of length $\Delta T\lesssim 10^{6}\,\mathrm{s}$. In this work, we assume that computational resources are not a limiting factor and instead consider the loss in performance of the semi-coherent search compared to the fully coherent search, estimated as a function of $N$. For a waveform model to safely be used in a semi-coherent search, it must remain phase-accurate over the duration of each segment, and so the different dephasing times in the AK and AAK models will determine the maximum $\Delta T$ for each model and the corresponding loss in performance. This section will discuss whether either model is sufficiently accurate for real LISA data analysis.

For an EMRI template bank $\{h_{i}\,|\,i=1,2,\ldots,N_{\mathrm{bank}}\}$, we may define the fully coherent detection statistic for each template as

[TABLE]

i.e. the template SNR. If the measured data consists solely of Gaussian noise ($s=n$), it is straightforward to show that $\rho_{i}\sim\mathcal{N}(0,1)$ (each template SNR is distributed as a zero-mean normal random variable). In the presence of an EMRI signal ($s=h_{i}+n$), we have $\rho_{i}\sim\mathcal{N}(A,1)$ for the corresponding template, where $A:=\langle h_{i}|h_{i}\rangle^{1/2}$ is the signal amplitude. A detection is claimed if any $\rho_{i}$ exceeds a predetermined threshold $\rho_{*}$, which may be set by fixing a desired false-alarm probability $P_{F}$:

[TABLE]

In practice, $\rho_{i}$ is not used as a detection statistic because it is computationally cheaper to analytically search over several extrinsic parameters (rather than having to generate templates that vary in those parameters). Each template may be written as

[TABLE]

where the normalised templates $u_{i}$ satisfy $\langle u_{i}|u_{i}\rangle=1$. The waveform amplitude $A$, time-of-arrival $t_{c}$ and phase offset $\phi_{c}$ are extrinsic parameters that may be searched over for each $u_{i}$ at negligible additional cost BCV2003 .

We may now define a fully coherent phase-maximised detection statistic as

[TABLE]

i.e. $\rho_{i}$ maximised over $\phi_{c}$. Unlike the template SNR, $\tilde{\rho}_{i}$ is not normally distributed; in the no-signal case $s=n$, it follows a Rayleigh distribution with unit scale parameter, and has the probability density function BCV2003

[TABLE]

If $s=h_{i}+n$, then $\tilde{\rho}_{i}$ for the corresponding template follows a Rice distribution with unit scale parameter and offset parameter $A$; its probability density function is BCV2003

[TABLE]

where $I_{0}$ denotes the order-zero modified Bessel function of the first kind.

A detection threshold $\tilde{\rho}_{*}$ for the new statistic may be chosen in similar fashion to (59):

[TABLE]

In the following analysis, we consider a desired false-alarm probability of $10^{-3}$ for the entire template bank. By approximating $\{\tilde{\rho}_{i}\,|\,i=1,2,\ldots,N_{\mathrm{bank}}\}$ as a set of independent random variables, we may simply reduce this value by a trials factor of $N_{\mathrm{bank}}$ to obtain $P_{F}$ (the desired false-alarm probability for a single template). We may also assume that the time-of-arrival has been maximised over (e.g. using fast Fourier transforms BCCS1998 ), which incurs an additional trials factor of $N_{\mathrm{time}}\sim 10^{8}$ (for $t_{c}$ offsets of 1 s). Hence we set

[TABLE]

With $\tilde{\rho}_{*}$ fixed by (65), we now consider the detection probability

[TABLE]

i.e. the probability in the presence of a signal that $\tilde{\rho}_{i}>\tilde{\rho}_{*}$ for the corresponding template. The detection probability as a function of $A$ is shown as the blue curve in the top panel of Fig. 7, and the critical signal amplitude needed to achieve $P_{D}=0.9$ is $A_{c}\approx 17$. This threshold is the SNR required for an EMRI to be detected by an idealised fully coherent search; for a semi-coherent search, $A_{c}$ will increase with $N$ as detection sensitivity is lost.

In a semi-coherent search where the data is split into $N$ segments of equal length, the phase-maximised detection statistic $\tilde{\rho}_{i}$ is calculated as before but for each segment $j$, and the semi-coherent detection statistic is then taken to be the quadrature sum

[TABLE]

Although each of the $\tilde{\rho}_{i,j}$ follows a Rayleigh/Rice distribution, the distribution of $\Upsilon_{i}$ is analytically intractable. However, in the limit $N\to\infty$, the central limit theorem guarantees that $\Upsilon_{i}$ will be normally distributed. We find empirically that for $N\gtrsim 70$, the distribution of $\Upsilon_{i}$ is well approximated as Gaussian; the mean and variance in the respective absence/presence of a signal are given by

[TABLE]

[TABLE]

with $\mu_{k}\approx 2.00$ and $\sigma_{k}\approx 1.45$. The functional forms of $(\mu_{1},\sigma^{2}_{1})$ are motivated by the corresponding expressions for a noncentral chi-squared distribution.

Using this normal approximation for the distribution of $\Upsilon_{i}$, the detection threshold $\Upsilon_{*}$ is set by the fixed false-alarm probability (65) as

[TABLE]

while the detection probability is given by

[TABLE]

Along with the fully coherent expression (66) in the top panel of Fig. 7, the semi-coherent detection probability (71) is also plotted for several values of $N$; the threshold signal amplitude $A_{c}$ for which $P_{D}=0.9$ in each case is indicated by the corresponding vertical line.

In the bottom panel of Fig. 7, $A_{c}$ is shown as a function of $N$; through a reduced chi-squared fit, we find that this relationship is well approximated as the power law

[TABLE]

in the range $1\leq N\lesssim 10^{4}$ (well beyond the range plotted in Fig. 7). As the number of segments is increased, the threshold SNR at which an EMRI can be detected increases, i.e. the search becomes less sensitive and a greater number of events will be missed. However, (72) taken at face value suggests that for $N\lesssim 55$, the semi-coherent search is more sensitive than the fully coherent search ($A_{c}<17$). This conclusion is obviously incorrect, and is due to the breakdown of the normal approximation for the distribution of $\Upsilon_{i}$ at $N\lesssim 70$ (indicated on the plot by the transition to a thinner black curve).

Since the AK model typically dephases over a few hours ($\sim 10^{4}\,\mathrm{s}$), employing AK waveforms in a semi-coherent detection search would necessitate the use of $N\sim 10^{4}$ segments to cover the full mission lifetime of $\sim 10^{8}\,\mathrm{s}$. This incurs an intolerable loss of performance; the threshold SNR at which an EMRI can be detected would be raised from 17 to 59, where this number is estimated by using the power law (72) and dividing by an assumed model accuracy of 0.97 over the dephasing time. On the other hand, only $N\sim 10^{2}$ segments are required for the AAK model with its typical dephasing time of two months ($\sim 10^{6}\,\mathrm{s}$), which would raise the threshold SNR from 17 to 20. It was also found previously in GEA2004 that $N\sim 10^{2}$ is the minimum number of segments needed for a computationally feasible semi-coherent search.

In moving from the AK model to the AAK model, the lowering of the threshold SNR from 59 to 20 has a profound effect on the likely number of detected EMRIs. Tab. 1 shows how many events LISA would expect to observe, using the thresholds appropriate for both kludges. Numbers are reported for the 12 astrophysical EMRI-population models that were recently considered in BEA2017 . These populations make different assumptions about the characteristics of the MBHs that play hosts to EMRIs, the typical masses of the compact objects involved, and the amount of mass that the central black holes can accrete; we refer the reader to BEA2017 for full details. Event SNRs are computed with two different plunge criteria: either the innermost stable circular orbit for a Schwarzschild black hole (as in the original AK paper BC2004 ), or that for a Kerr black hole. As argued in BEA2017 , these two assumptions should give values that bracket the true SNR. Detected events are further divided into four mass ranges for the central black hole.

We see that some detections would still be expected if the AK model is used to analyse LISA data, although significantly fewer (by at least an order of magnitude) than if the more faithful AAK model is employed. Moreover, we would lose the events at the edges of the distribution (i.e. in the lowest and highest mass ranges), and so the astrophysical information provided by the LISA EMRI population would be reduced. The AAK model, on the other hand, would find a significant number of events under most of the population models, with its required threshold SNR of 20 for a semi-coherent search being compatible with the thresholds that have traditionally been assumed when assessing LISA’s capability for EMRI detection. Results for a fully coherent search are not shown in Tab. 1, but reducing the threshold SNR from 20 to 17 would increase the event rate by only a modest amount (50% or less); furthermore, it is unlikely that the computing power needed for such a search will be available even in the mid-2030s.

V Conclusion

We have developed an augmented variant of Barack & Cutler’s widely used EMRI kludge waveform model BC2004 ; the new AAK model retains the speed of its predecessor, while matching the phase evolution of more accurate but slower kludges over a significant fraction of the inspiral. With the latest implementation of the model released online at github.com/alvincjk/EMRI_Kludge_Suite as part of a kludge software suite, AAK waveforms will hopefully see widespread use in the next round of mock LISA data challenges.

One existing deficiency in the AAK model is the ill-defined nature of the fundamental frequency map (41)–(43) at the last stable orbit, due to the divergent Kerr frequencies; this complicates both plunge detection and the specification of orbital parameters at plunge. Another limitation is that the mode-sum approximation (6) becomes more expensive than the quadrupole formula itself at high eccentricities ($e_{0}\gtrsim 0.5$). Work is ongoing to resolve these issues, and to make the model as streamlined and robust as possible for the mock data challenges.

Regardless, the present implementation of the AAK model shows significantly improved accuracy over larger fractions of the inspiral as compared to the initial AAK implementation CG2015 (which itself extends the dephasing time of the original AK model from hours to months). The two-month overlaps of AAK waveforms with NK waveforms for a variety of EMRIs with different compact-object mass, black-hole spin and initial eccentricity are typically increased from $\lesssim 0.95$ (in the initial implementation) to $>0.97$ by the techniques presented in this paper. Computational efficiency for the present implementation is also retained to within 1% of the AK and initial AAK models, with all analytic kludges able to generate generic waveforms 5–15 times more quickly than the NK model (except when $e_{0}\gtrsim 0.3$).

We have also considered the performance of the AAK model in a data analysis application: the semi-coherent EMRI detection search proposed in GEA2004 . An analytic estimate is provided for the threshold SNR needed to detect an EMRI with a false-alarm probability of $\sim 10^{-3}$ when using a bank of $\sim 10^{40}$ templates. This threshold is 17 for a fully coherent search with templates that are 97% accurate over the inspiral lifetime, although such a search is likely to be out of reach computationally. If the AK model is used in a semi-coherent search, the threshold rises to 59 and around 95% fewer EMRI events will be detected (under various population models); however, a semi-coherent AAK search requires a lower threshold of 20 due to the model’s longer dephasing times, and will yield at least an order of magnitude more events than the AK search while remaining computationally feasible. This suggests that unlike the AK model, the AAK model might realistically be employed in actual LISA data analysis without much loss in detection sensitivity.

Acknowledgements.

We thank Stanislav Babak, Leor Barack, Christopher Berry and Anthony Lasenby for helpful discussions and/or comments on the manuscript. AJKC’s work was supported by the Cambridge Commonwealth, European and International Trust. CJM has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 690904, and from the STFC Consolidator Grant No. ST/L000636/1.