{'fraction_non_alphanumeric': 0.03668676779215259, 'fraction_numerical': 0.010201099806072842, 'mean_word_length': 4.416402877697842, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Towards adiabatic waveforms for inspiral into Kerr black holes: II. Dynamical sources and generic orbits

This is the second in a series of papers whose aim is to generate ``adiabatic'' gravitational waveforms from the inspiral of stellar-mass compact objects into massive black holes. In earlier work, we presented an accurate (2+1)D finite-difference time-domain code to solve the Teukolsky equation, which evolves curvature perturbations near rotating (Kerr) black holes. The key new ingredient there was a simple but accurate model of the singular source term based on a discrete representation of the Dirac-delta function and its derivatives. Our earlier work was intended as a proof of concept, using simple circular, equatorial geodesic orbits as a testbed. Such a source is effectively static, in that the smaller body remains at the same coordinate radius and orbital inclination over an orbit. (It of course moves through axial angle, but we separate that degree of freedom from the problem. Our numerical grid has only radial, polar, and time coordinates.) We now extend the time-domain code so that it can accommodate dynamic sources that move on a variety of physically interesting world lines. We validate the code with extensive comparison to frequency-domain waveforms for cases in which the source moves along generic (inclined and eccentric) bound geodesic orbits. We also demonstrate the ability of the time-domain code to accommodate sources moving on interesting non-geodesic worldlines. We do this by computing the waveform produced by a test mass following a ``kludged'' inspiral trajectory, made of bound geodesic segments driven toward merger by an approximate radiation loss formula.


I. INTRODUCTION
A. Background The extreme mass ratio limit of general relativity's two-body problem has been a major focus of work in recent years. This limit corresponds to a stellar mass compact object that orbits and perturbs a massive black hole. The system generates gravitational waves (GWs) which drive the small body to inspiral into the large black hole. Measuring such "extreme mass ratio inspiral," or EMRI, events is a major goal for space-based GW antennae, particularly the LISA mission 1 . EMRIs should be measurable to a redshift z ∼ 0.5 − 1. The event rate at this range is estimated to be high enough that a multiyear LISA mission should measure dozens to hundreds of EMRI events [1]. Because the smaller body only slightly perturbs the larger black hole's spacetime, EMRI GWs are expected to provide an exceptionally clean probe of black hole properties. We expect to use EMRIs to measure black hole masses and spins with extremely good accuracy [2], and even to test how well the spacetime meets the rather stringent constraints that the "no-hair" theorems of general relativity impose on black holes [3,4,5,6]. Understanding EMRI sources will require us to compare measured waves with theoretical models that are as accurate as possible. This goal motivates much recent EMRI work. The waves are sufficiently complicated that simply detecting them in LISA's datastream will be a challenge. Techniques for finding these events are currently being developed and tested through the "Mock LISA Data Challenges", or MLDCs (see Refs. [7,8] for overviews of recent MLDCs). An important input to these challenges (and to the development of EMRI measurement techniques more generally) are waveform models that capture the true complexity of EMRI events (see [9,10] for discussion of recent work to include EMRI waves in the MLDCs).
This paper presents a further step in our program to construct accurate EMRI wave models. As discussed in the introduction to Ref. [11] (Paper I), our goal is to make "adiabatic" waveforms -waveforms built by separately treating the long-time dissipative evolution and the short-time conservative motion. In our present analysis, we take the short-time motion to be a geodesic orbit of the background spacetime; our approach thus amounts to approximating the inspiral trajectory as a sequence of geodesic orbits. As discussed by Pound and Poisson [12], this limit is more properly a "radiative" or "dissipative" approximation, since we do not include conservative self interactions. It may be possible to augment this analysis with at least some conservative effects [13], so we believe the program we are developing is capable of building truly adiabatic inspiral waveforms as described in [12].
We will describe our goal as "adiabatic" waveforms, but the reader should bear in mind that the approximation we are currently developing is more restricted than this.
Geodesic orbits are described (up to initial conditions) by three conserved constants: energy E, axial angular momentum L z , and "Carter constant" Q. Using black hole perturbation theory, we compute the rate at which these three constants evolve; fast and accurate frequencydomain codes make it possible to compute these rates of change fairly easily [14,15,16]. We then build the parameter-space trajectory [E(t), L z (t), Q(t)] followed by the small body; choosing initial conditions, it is simple to build the coordinate-space worldline [r(t), θ(t), φ(t)] of a particular inspiral. From this worldline, we build the source to a time-domain code. The output of this code is, at last, our model EMRI wave.

B. Time-domain black hole perturbation theory
Since the frequency-domain portion of this program is already well in hand, our current focus is on the timedomain code. In essence, our goal is to build a code which takes as input any physically reasonable worldline, and provides as output the waveform produced by a small body on this worldline. In Paper I, we demonstrated an accurate (2+1)D numerical code to solve, in the time domain, the wave equation for curvature perturbations to a black hole -the Teukolsky equation [17]. Our code evolves the Weyl curvature scalar Ψ 4 , constructed by projecting the vacuum curvature onto appropriate components of a null tetrad; see Paper I for details. The azimuthal dependence of Ψ 4 is separated out (due to the φ symmetry of black holes); the dependence on the Boyer-Lindquist coordinates r, θ, and t is found by evolving Ψ 4 on a numerical r-θ grid.
As is common in black hole perturbation theory, we treat the smaller body as a Dirac-delta point particle, leading to a singular source for the Teukolsky equation. In the frequency domain, the delta can be dealt with analytically, and presents no great challenge. By constrast, accurately computing the effect of a sharp source on the time-domain code's numerical grid can be extremely challenging. In Paper I, we presented a new technique for treating the singular source term. Our innovation was to model the delta as a series of finite impulses, with the largest impulse located close to the delta's argument, falling off rapidly as we move away from this "central" spike. Importantly, this approach allows us to accurately model the derivatives of the delta function. Since the Teukolsky equation source depends on first and second derivatives of the delta (as well as the delta itself), this appears to give us an accuracy boost relative to other finite-difference delta representations (such as a truncated Gaussian), which may accurately capture the delta's behavior, but not do so well with the derivatives.

C. This paper
Paper I focused on the properties of this new source representation. To clarify this focus, we studied very simple orbits: We only considered the (astrophysically unlikely) case of circular, equatorial black hole orbits. We now extend this to include inclined, eccentric and generic orbits, as well as non-geodesic inspiral sequences.
A particle in a circular, equatorial orbit has constant radial and angular coordinate, confining it to a fixed location on the r-θ grid. Eccentricity means that the orbit oscillates radially, crossing radial grid zones. Similarly, orbital inclination results in angular grid crossing. We quickly discovered that these new motions introduce high frequency numerical noise. This noise can be controlled by combining a low pass filter with a higher order discretization of the delta function; details are given in Sec. II. Aside from this mild extension of the basic formalism presented in Paper I, it was not terribly difficult to use our new source term to handle a broad class of astrophysically interesting orbits. To validate our results, we present in Sec. III extensive comparisions with waveform snapshots computed in the frequency domain [14], demonstrating graphically and quantitatively (with appropriate overlap integrals) excellent agreement between the two techniques.
As extensively discussed in the introduction to Paper I and here, our goal is to compute the waves from inspiral of a small body through a sequence of orbits. As a proof-of-concept demonstration of the feasibility of this idea, we present a simple example of inspiral in Sec. IV. In this example, we evolve through our geodesic sequence using a "kludged" approximation to the rates of change of orbital constants, using the code described in Refs. [18,19]. These waveforms are not reliable EMRI models, but they illustrate the ease with which we can handle the effect of radiation emission on the motion of the source. Computing waves from an inspiral is no more of a computational challenge than computing waves from a bound geodesic.
The next step will be to combine accurate radiative backreaction with our time-domain solver to compute "adiabatic" EMRI waveforms (albeit ones that still neglect conservative self interactions). Plans for this next step are described in our final summary, Sec. V.

II. DYNAMICALLY VARYING DISCRETE DELTA FUNCTIONS
In Paper I, we presented a method for representing a Dirac delta function and its derivatives on a discrete numerical grid. In that paper, we only considered a delta with fixed radial and angular position. Naive application of the discrete delta models presented in Paper I leads to instabilities when the particle moves in the numerical grid. The following argument outlines the root cause of these instabilities. Consider the function δ[x − α(t)], where x k ≤ α(t) ≤ x k+1 ; i.e., the delta's peak varies with time and lies between two discrete grid points. Let x i represent any discrete point on our grid, and let h = x k+1 − x k = x k − x k−1 be the grid resolution. Naive application of the results from Paper I might lead us to model the delta function with the impulse weights (This "two impulse" delta is in fact just the simplest representation we developed in Paper I, but is useful for the following discussion.) Each t n defines a time slice of our r-θ grid. As α varies from one time slice to another, so do the coefficients at x k and x k+1 . The frequency spectrum of δ k (t n ) and δ k+1 (t n ) will reflect the amount of variation in α. A large variation in α will produce a high frequency component in the Fourier transform of the time series of each weight. These variations couple to the time derivatives in the homogeneous part of the Teukolsky equation. Consequently, the solution contains spurious high frequency features of numerical origin. Consider the extreme limit of this effect: α changes so rapidly that the delta's peak moves across a grid zone in a single time slice: so that The weight of the delta function very suddenly becomes zero at x k+1 as we step from t = t 0 to t = t 1 ; likewise, the weight at x k−1 very suddenly jumps from non-zero to zero in this step. The coupling of this sudden change to numerical time derivatives drives instabilities in our code, in a manner reminiscent of the initial burst of radiation that occurs due to the sudden appearance of the particle at the start of our evolution; see Fig. 2 of Paper I. This problem is substantially mitigated by using a delta representation with a wider stencil; examples of this are described in Paper I. Wide stencils reduce the amount by which each weight changes from step to step, thereby reducing numerical noise. Another useful tool is to increase the order of the delta representation, thereby increasing the smoothness of the delta and its derivatives. This is particularly important since the Teukolsky equation is a second-order differential equation; some smoothness in the derivatives is necessary to prevent the differential operator from seeding excessive noise. Finally, residual high frequency noise can be removed by convolving the source with a low pass filter 2 . These three techniques are each described in the following subsections.
Each of these techniques smear out the delta function, pushing us away from the idealization of a zero width singularity. Choosing between stability (which tends to push us to a wider delta) and faithful representation of the singularity (which pushes us to a narrow delta) leads us to an optimization problem; we tune our delta representation in a way that (hopefully) minimizes numerical noise and maximizes accuracy. Note also that, in addition to high-frequency noise generated by abrupt movement of the delta across the grid, spurious excitations of the quasinormal modes of the black hole also appear due this motion. This source of "noise" appears to be controlled by grid resolution -wider grids lead to less pointlike deltas, which spuriously excite these modes. This spurious contribution to the EMRI waves can be mitigated with a form of Richardson extrapolation [20]. We discuss this further in Sec. III and the Appendix.

A. Higher order delta functions
Discrete delta representations based on linear and cubic interpolation were derived in Paper I. We now extend this process to arbitrary polynomial order, equipping us with an entire family of discrete delta functions.
As in Paper I, we start from the defining integral, Let x k+n−1 ≤ α ≤ x k+n ; the reason for our somewhat idiosyncratic choice of subscripts will become clear as we proceed. For clarity, we will not explicitly write out the time dependence of α; the reader should bear in mind that α = α(t) in all that follows. Rewriting Eq. (2.7) as a sum over a finite step size, we have The function f (α) can be approximated by the Lagrange interpolating polynomial, where 2n is the order of interpolation and Inserting this in Eq. (2.8) leaves us with comparing coefficients of f (x i ) allows us to read off δ i , The weights for derivatives of the delta function can be obtained similarly. Writing the identities as sums gives us We now insert centered finite difference formulae for the Expressions (2.18) and (2.19) are in a form that makes it simple to read off δ ′ i and δ ′′ i . For example, δ ′ j can be calculated by setting f (x j ) = 1 and f (x l ) = 0, l = j. It is straightforward to verify that setting n = 1 and n = 2 reproduces the weights given by the two-point linear hat and the cubic formulae (described in Paper I) respectively. We also note that the delta derivative coefficients are non-zero for i ∈ [k − 1, k + 2n].

B. Wider stencils at a given interpolation order
In Paper I, we generalized the two-point linear hat delta function such that it can be represented over a larger number of points. Similarly, we develop a procedure to widen the stencil of the generalized model obtained from Eqs. (2.13), (2.18), and (2.19).
Consider a model for δ i obtained from Eq. (2.13) for some n = m. Then, Our goal is to widen this representation by some integer factor w such that the coefficients are non-zero for a wider range of grid points. Let us label the weights of this wider representation by δ w i , with δ w i = 0 for i ∈ [k, ..., k + 2wm − 1]. It should be emphasized that this is different from simply using Eq. (2.13) with n = wm; we have not changed the polynomial order, it remains fixed at 2m.
For concreteness, let us choose w = 2, doubling the number of points in the delta representation. We infer the coefficients δ 2 i at gridpoints i = k, k + 2, k + 4, . . ., k + 4m − 2, by widening the grid by a factor of two: We evaluate δ i with h → 2h, x k+j → x k+2j to get Finally, we need δ 2 i at the intermediate points i = k + 1, k + 3, . . ., k + 4m − 1. We do this by exploiting the translational symmetry of the problem. Momentarily reinsert the time dependence of the δ's and α. Now consider the hypothetical situation where i.e, α(t) changes by a grid spacing from t 0 to t 1 . We must have We can turn this equation the other way around to read off the coefficient δ 2 k+2j+1 at t 0 : Simply replace α(t) with α(t)−h in the formula for δ 2 k+2j (t 0 ) to obtain δ 2 k+2j+1 (t 0 ). Since there was nothing special about our time slice, t 0 , we find for any moment t n . Though we chose w = 2 for concreteness, the above argument can be generalized to any integer w. Since our result holds for all time slices, we again suppress the time dependence to obtain expressions for any integer w: (2.28) These techniques carry over to the derivatives as well: (2.32) These should be used with Eqs. (2.13), (2.18) and (2.19) to widen the Teukolsky source term by any factor w.

C. Smoothing the source with a Gaussian filter
Further control of numerical noise can be achieved by filtering high frequency components in the source term. This requires a convolution of the source with a discrete low pass filter. We use a Gaussian filter because it maximizes the uncertainty principle -it can be localized in both position and frequency with greatest efficiency.
Consider a source of the form (2.33) Delta function identities allow us to rewrite this as On a discrete grid, this becomes If the delta function and its derivatives span 2n + 2 grid points, with x k+n−1 ≤ α ≤ x k+n , then s i = 0 for i ∈ [k − 1, . . . , k + 2n]. The source s i is zero everywhere else on the grid. The Gaussian filter is given by this normalization guarantees that the integrated value of any function convolved with the filter is unchanged. We now convolve the source with the filter to obtain where sg k is the smoothed source term. This indicates A wide filter spreads the source over a large domain on the numerical grid and thus increases errors, although it eliminates spurious harmonics. We have found that using a wide stencil followed by a narrow Gaussian smoother works very well to reduce numerical noise and minimize errors from an insufficiently pointlike source.

D. Order of convergence of the filtered delta
Paper I discussed in detail the convergence of a code that uses a discrete delta. Crucial background is given by Ref. [21] and summarized in Paper I. The key point is that the moment controls the delta's convergence properties. Clearly, M 0 = 1 (otherwise the delta is not properly normalized); in the continuum limit, M r = 0 for r > 0. For the discrete delta, the smallest non-zero value of r for which M r = 0 sets the order of convergence. In particular, if M r = 0, then a code which uses this delta will be no higher than rth-order convergent.
We now show that, if a delta representation is secondorder convergent before smoothing with the Gaussian filter (M 0 = 1, M 1 = 0, M 2 = 0), it will remain secondorder convergent after smoothing. Upon convolving the discrete delta with the Gaussian smoother, we find (2.41) Let us denote the moments of the smoothed delta by M g r . As discussed in Sec. II C, the convolution does not change the delta's normalization as long as the Gaussian filter is itself properly normalized; thus We now examine the next higher moment of the smoothed delta: The first term on the final line of (2.43) gives zero: Since The Gaussian filter's symmetry property c j = c −j has been applied in the last step. Hence, we find M g 1 = M 1 = 0.
Evaluating the second moment proceeds similarly, but we find in the end terms involving p j=−p j 2 c j which do not vanish. Thus, M g 2 is the first non-vanishing moment of the discrete delta, demonstrating that the Gaussianfiltered discrete delta function exhibits second-order convergence. The argument can be extended to the delta derivatives as well. The smoothed Teukolsky source term will thus be second-order convergent.

III. WAVEFORMS AND COMPARISONS FOR GENERIC GEODESIC KERR ORBITS
We now present the waveforms generated by a point particle in a geodesic orbit around a Kerr black hole. The code used to generate these waves is discussed in detail in Paper I; the only important change to that discussion is that the source term uses the techniques presented in Sec. II above. We begin by reviewing Kerr black hole geodesics, sketching the numerical scheme used to solve the equations of motion. We then examine differ-ent classes of eccentric and inclined orbits and compare the waveforms against those obtained from a frequencydomain code whose details are given in Ref. [14]. We compute the correlation between the two waveforms in order to measure our level of agreement with frequencydomain waveforms.
Our numerical grid is laid out in Boyer-Lindquist coordinats and uses (δr, δθ, δt) = (0.04M, π/60, 0.02M ) for the radial, angular and temporal resolutions. The source term is constructed using Eqs. (2.13), (2.18) and (2.19) with n θ in the range 3-9 (depending on the orbit) for the angular delta-function and n r = 2 for the radial delta. We use a Gaussian filter of width b = 1.5δθ to smooth higher harmonic noise.

A. Geodesics in Kerr spacetime
The source term for the time-domain code takes as input the worldline of the perturbation's source. Here, we neglect radiation reaction and assume that the point particle follows a bound geodesic trajectory around the central massive black hole. This bound trajectory can be computed by numerically integrating the geodesic equations. We now briefly review how we massage the geodesic equations to put them into a form that makes for accurate numerical calculation; this material is presented in greater depth in Sec. IIC of Ref. [19].
The normal "textbook" presentation of the equations governing Kerr black hole geodesics is [See, e.g., Ref. [22], Eqs. (33.32a-d).] Here, Σ = r 2 + a 2 cos 2 θ, ∆ = r 2 − 2M r + a 2 (where a = | S|/M is the black hole's spin angular momentum per unit mass). The constants of motion are orbital energy E, axial angular momentum L z , and Carter constant Q. This form of the equations of motion is not well suited to numerical studies; in particular, dr/dτ and dθ/dτ pass through zero and change sign when the orbiting body goes through turning points associated with those motions. A handy way to eliminate these problems is to eliminate the turning points by remapping the coordinates r and θ to parameters which accumulate secularly. The following parameterization, inspired by the Newtonian limit, has been found to work extremely well even deep in the strong field of rapidly rotating black holes: cos θ = cos θ min cos χ . (3.6) In the Newtonian limit, p is the orbit's semi-latus rectum, and e is its eccentricity; θ min is the minimum value of θ reached by the orbiting body, and is used to define the orbit's inclination θ inc Once E, L z , and Q are specified, p, e, and θ inc are fully determined. It is then a straightforward matter to turn Eqs. (3.1) and (3.2) into expressions for dψ/dτ and dχ/dτ ; see Ref. [19] for details. The resulting expressions behave extremely well for all bound orbits outside the black hole's event horizon. A numerical integrator for these variables allows us to compute the dynamics of our orbiting body's Teukolsky equation source term. Before moving on, we note that, within the context of the dissipative-only or radiative approximation to inspiral, it is simple to modify these equations to build the worldline of an inspiralling body: We simply allow the orbital "constants" (E, L z , and Q; or, p, e, and θ inc ) to evolve according to the inspiral law. Reference [19] uses approximate radiation reaction, based on fits to strongfield radiation reaction calculations in regimes where it is well understood, to compute the inspiral worldlines which underlie the "kludge" waveforms. We use this prescription for evolving the constants in Sec. IV to demonstrate this code's ability to compute inspiral waves.

B. Comparison with frequency-domain waveforms
To validate our waveforms, we compare with the "snapshots" generated using the frequency-domain code described in Ref. [14]. This code uses the fact that bound Kerr geodesics are fully described by three frequencies (radial Ω r , latitudinal Ω θ , and axial Ω φ ) to build the waveform from a geodesic orbit as a sum over harmonics of these frequencies [23]. Since both the time-domain and frequency-domain codes solve the same master equation, they should produce identical waveforms for identical orbits, so long as each code is sufficiently accurate.
To quantify the accuracy with which a time-domain waveform X = (x 1 , x 2 , . . . , x n ) agrees with a frequencydomain waveform Y = (y 1 , y 2 , . . . , y n ), we use the following correlation measure: The sums in all cases are from i = 1 to i = n.) This coefficient is identical to the match between two waveforms defined by Owen [24] in the white noise limit [noise spectral density S h (f ) = constant]. One might expect the waveforms' mean valuesx andȳ to equal zero. However, finite duration effects can make these quantities slightly non-zero, so it is useful to explicitly do this subtraction.
A useful reformulation of Eq. (3.8) is Note that r XY is always between −1 and 1; a value close to 1 indicates that the two waveforms are well correlated. Note also that the correlation depends on how many points n are used in comparing the two waveforms (or equivalently, the span of time over which we compare the waves). We have found that as long as n > ∼ several hundred, we get consistent results: Changing n for a given comparison only causes small variations in the fourth significant digit of r XY .
It is of course possible to concoct other measures of how well two waveforms agree. Ideally, disagreements between waveforms should be quantified in terms of their observational significance. For example, Cutler and Vallisneri have demonstrated that it is not unusual for waveforms with a match of 0.9999 to differ significantly in their estimates of the parameters which describe the source [25]. For our present purpose, r XY is sufficient to demonstrate that our time-domain code produces high quality waveforms; whether they are sufficiently high quality to be used for GW measurement purposes will need to be re-examined at a later time.
An important step in producing accurate waveforms is to perform runs at multiple resolutions, then estimate (and eliminate) the waveform error using a form of Richardson extrapolation [20]. This plays a crucial role in reducing "noise" from spurious excitation of the large black hole's quasinormal modes. The details of this extrapolation technique are described in Appendix A. Tables I, II, III, IV, V and VI list the correlation coefficients for the m = 2 and m = 3 azimuthal modes of different classes of orbits. The coefficient is greater than 0.99 for a large fraction of parameter space. Time domain runs corresponding to each column required about 125 CPU hours on an Apple MacPro processor. That code was compiled using the Intel C++ compiler. The frequency domain code's cost is about 3-4 CPU hours per waveform when attempting to get both asymptotic energy fluxes to accuracies of about 0.1% to 1% on a machine using a 3.2 GHz Pentium 4 Xeon processor. We also show (Figs. 1, 2, and 3) examples of the waves, computed with both time-and frequency-domain codes, to give the reader a visual sense of the overlap.   Table I.

IV. INSPIRAL WAVEFORMS
Having demonstrated that the finite-impulse source works well for astrophysically relevant generic black hole orbits, we now examine how well we do evolving through a sequence of such orbits. Since each orbit in the sequence is no different than the orbits that we validated against in Sec. III B, we anticipate no great difficulty here. Indeed, the biggest challenge is choosing a method to evolve   III: Correlation between time-and frequency-domain waveforms for the m = 2 mode for a range of generic orbits. All symbols have the same meaning as in Table I. To quickly produce results that are qualitatively correct, we presently make this tra- jectory using the "kludge" inspiral treatment described in Ref. [18], and used to make model waveforms in Ref. [19]. The "kludge" uses a somewhat idiosyncratic mix of post-Newtonian backreaction formulae combined with  numerical results from frequency-domain backreaction in the circular, inclined (e = 0, θ inc = 0) and eccentric, equatorial (e = 0, θ inc = 0) limits to estimate the properties of EMRI waves. By construction, the results agree very well with Teukolsky-based inspirals in those limits; for the generic case, they produce plausible inspirals. Figure 4 shows our waveform for a "kludge" inspiral. We took the large black hole to have spin a = 0.5M , and set the mass ratio to µ/M = 0.016. The orbit was initially chosen to have semi-latus rectum p = 10M , eccentricity e = 0.5, and inclination θ inc = 0.5 radians. This figure shows features reminiscent of the geodesic snapshots shown in Figs. 1, 2, and 3; in addition, one can clearly see evolution of the wave's properties. The increase in the wave's frequency, largely due to the decay of the orbit's semi-latus rectum, is quite clear. Perhaps less obvious is a signature of the eccentricity's decay. This is illustrated most clearly by comparing the lower left and lower right panels of Fig. 4, which zoom onto early and late portions of the inspiral. Early on, the waveform is dominated by a series of high-frequency bursts; these occur when the small body passes through periapsis and "whirls" most rapidly about the massive black hole. There is then a relatively quiet section while the body "zooms" out to apoapsis, and then comes in to "whirl" at periapsis again. As eccentricity shrinks, the difference between periapsis and apoapsis becomes smaller. The high-frequency bursts crowd closer and closer together, approaching a continuum sinusoid as the eccentricity approaches zero.
Although this inspiral model is somewhat unphysical, we expect that it shares many properties with true adiabatic inspiral waveforms. In particular, the spectral evolution of a wave like that in Fig. 4 should be quite similar to the evolution of real EMRI waveforms. It should be emphasized that computing the waveform shown in Fig. 4 required about as much computational effort as computing the geodesic snapshot waves, Figs. 1, 2, and 3 (modulo a factor ∼ 4-5 since the waveform in Fig. 4 lasts ∼ 4-5 times longer than the others). Given a robust code to generate the inspiral worldline of EMRI systems, the waveforms that our code produces should be a useful tool for examining issues in LISA measurement and data analysis.

V. SUMMARY AND FUTURE WORK
We have now shown that the finite impulse delta representation of the time-domain Teukolsky equation's source works very well for complicated and astrophysically relevant orbits. In our previous analysis ( [11], Paper I), we confined ourselves to the simplest circular, equatorial black hole orbits. The basic ideas from Paper I  work well even when the source arises from highly inclined and highly eccentric orbits, and when the source evolves through a sequence of those orbits. It is now a relatively straightforward matter to compute the waves arising from a body following any reasonably behaved worldline in the spacetime of a black hole.
The primary complication arising from these more generic orbit classes is that the orbiting body will cross zones within our numerical grid. The source thus becomes dynamical; the finite-impulse delta must likewise be dynamical to represent it. The evolution of the impulses that we use to represent the delta can seed numerical noise, reducing the calculation's accuracy. We have found that minor extensions of Paper I's basic techniques greatly mitigate the impact of this source of numerical noise. In particular, by using a higher-order representation (Sec. II A), the delta is smoothed enough that the coupling to the Teukolsky equation's second-order differential operators does not seed much error. Widening the delta's stencil (Sec. II B) also helps, since the fractional change in a given impulse will be less if the delta is represented by more impulses. Finally, residual high frequency noise not removed by these techniques can be taken out by convolving the Teukolsky source term with a low-pass (Gaussian) filter (Sec. II C). It's worth emphasizing that we smooth the entire source term, not just the delta function (which would arguably make our delta rather similar to the truncated Gaussian [26,27] which this technique was designed to improve upon).
Comparison with results from the frequency-domain [14] demonstrates that the waveforms generated with this source term are of very high quality (Sec. III). Visually, the waveforms lie on top of one another in every case that we have examined; a quantitative overlap integral demonstrates that waveforms from the two calculations are often more than 99% correlated. A key step in achieving such high quality results is to estimate the largest errors in our time-domain calculations, and then subtract that estimate from our result. We do this by performing these calculations at two different grid resolutions; under the assumption that our dominant error is quadratic in grid spacing, we then estimate the magnitude of our error (Appendix A). The excellent agreement we achieve with frequency-domain results validates this approach, at least for all the cases we have considered.
So far, our main physics accomplishment is excellent agreement between time-and frequency-domain approaches to waveform calculation. It should be emphasized, however, that for waveform calculations, there will be a large set of circumstances in which time-domain codes are more efficient. For generic orbits, a frequencydomain code may require the calculation and summation of many thousand multipoles and Fourier modes. A time-domain code "automatically" sums over all modes (except the m index), so that (in principle) it is no more difficult to compute the waves from a highly inclined, highly eccentric black hole orbit than from an orbit with modest inclination and eccentricity.
The real payoff of this tool will come when we allow the source to radiatively decay, evolving through a sequence of orbits. As a demonstration that this can be done, we use a "kludged" inspiral to compute a body's inspiral, and then use that inspiral as the source for our timedomain solver in Sec. IV. Though not a physically accurate inspiral, this scenario shares many properties with the actual adiabatic inspiral. In particular, it demonstrates the computational advantage of a robust timedomain code for computing inspiral waveforms, given the worldline the inspiraling body follows.
Future work will address our goal of complete waveforms for the EMRI problem, in the context of the dissipation-only approximation to EMRI dynamics. We have recently extended our frequency-domain code to include the evolution of Carter's constant in the radiative backreaction limit [15], and will use this code to produce the radiation reaction data describing an inspiraling body. With this step in hand, no issue of principle stands in the way of coupling the time-and frequency-domain approaches to make usefully accurate EMRI waveforms.  Here we describe the variation of Richardson extrapolation which we use to estimate and eliminate the largest errors arising from our finite difference scheme. In Ref. [11], we showed that our algorithm is second order convergent. This means that we can write the solution at any given resolution as Ψ c = Ψ t + a 1 δr 2 + a 2 δθ 2 + a 3 δrδθ + O(δ 3 ) , where Ψ c is the computed solution and Ψ t is the "true" solution. The final term O(δ 3 ) indicates that additional error terms will be third order in the grid spacing (and higher). The spatial and temporal dependences of Ψ c and Ψ t have been suppressed. We now perform runs at two different resolutions, (δr 1 , δθ 1 ) and (δr 2 , δθ 2 ), with all other parameters fixed. The resolutions are chosen such that Neglecting higher order terms, the two results can be written Ψ c1 ≃ Ψ t + a 1 δr 2 1 + a 2 δθ 2 1 + a 3 δr 1 δθ 1 , (A3) Ψ c2 ≃ Ψ t + a 1 δr 2 2 + a 2 δθ 2 2 + a 3 δr 2 δθ 2 .
To the extent that neglect of higher-order errors is warranted, this estimates the largest source of error. Using Eq. (A3) we can now estimate the "true" value: Ψ t ≃ Ψ c1 − (a 1 δr 2 1 + a 2 δθ 2 1 + a 3 δr 1 δθ 1 ) Figure 5 illustrates the improvement that this variant of Richardson extrapolation can yield. We plot h + at two different resolutions: (δr 1 , δθ 1 ) = (0.04, π/60) and (δr 2 , δθ 2 ) = (0.026667, π/90). We also show the extrapolated waveform, and the frequency-domain prediction. The particle is in a geodesic orbit with parameters p = 6M , θ inc = 45 • , e = 10 −4 and the black hole has a spin of a = 0.9M . The two time-domain calculations each differ noticeably from the frequency-domain result; the extrapolated waveform by contrast agrees very well. This excellent agreement can be regarded as a modified three-level convergence test, whose first two levels are the time domain waveforms and third level is the frequency domain waveform. If the code were not second order convergent, our assumption for the functional form of the errors in Eq. (A7) would be erroneous. This would lead to a substantial disagreement between the extrapolated and frequency domain waveforms.