Null infinity waveforms from extreme-mass-ratio inspirals in Kerr spacetime

We describe the hyperboloidal compactification for Teukolsky equations in Kerr spacetime. We include null infinity on the numerical grid by attaching a hyperboloidal layer to a compact domain surrounding the rotating black hole and the orbit of an inspiralling point particle. This technique allows us to study, for the first time, gravitational waveforms from large- and extreme-mass-ratio inspirals in Kerr spacetime extracted at null infinity. Tests and comparisons of our results with previous calculations establish the accuracy and efficiency of the hyperboloidal layer method.

We describe the hyperboloidal compactification for Teukolsky equations in Kerr spacetime. We include null infinity on the numerical grid by attaching a hyperboloidal layer to a compact domain surrounding the rotating black hole and the orbit of an inspiralling point particle. This technique allows us to study, for the first time, gravitational waveforms from large-and extreme-mass-ratio inspirals in Kerr spacetime extracted at null infinity. Tests and comparisons of our results with previous calculations establish the accuracy and efficiency of the hyperboloidal layer method.

I. INTRODUCTION
Black hole perturbation theory plays a prominent role in obtaining physical insight into the quantitative behavior of astrophysical systems [1][2][3][4][5][6][7][8]. In recent years there has been considerable interest in computing the gravitational waveforms emitted by the radiation-reaction driven inspiral of a small, compact object into a large supermassive black hole. The waveforms from such extreme-mass-ratio inspirals (EMRIs) are expected to be observed by future, space-based, gravitational wave detectors such as the Laser Interferometric Space Antenna.
In these methods, the waveforms can be calculated by numerically solving the equations describing black hole perturbations. When the central, supermassive black hole is assumed to be rotating, the most common approach is to solve the Teukolsky equations that describe curvature perturbations of a Kerr spacetime [4,5].
A difficulty in solving the Teukolsky equations numerically-especially in the time domain-is the construction of suitable outer boundary conditions. Typically, the computational domain is artificially truncated to calculate the solution in a finite domain. One needs to provide boundary data along the timelike boundary of this domain such that the artificial boundary is transparent to outgoing waves from the interior. The outer boundary problem is already nontrivial for the evolution of scalar waves in Minkowski spacetime [31][32][33]. There has been substantial progress in constructing good outer boundary conditions for wave equations in Schwarzschild spacetime [34][35][36][37], which have been recently generalized to Kerr spacetime [38]. Such sophisticated boundary con-ditions, however, can be difficult to implement. To our knowledge, there is no numerical implementation of the improved boundary conditions in Kerr spacetime.
Another difficulty in the numerical solution of the Teukolsky equation is the extraction of gravitational waves as measured by idealized observers at null infinity. Astrophysical sources of gravitational radiation are typically thousands of light years away, whereas the computational truncation is performed at a few thousand Schwarzschild radii. One argues that the dynamics in the asymptotic domain is negligible. However, there are certain effects, such as the backscattering off curvature [39,40], where asymptotic properties of the solution are essentially different from the corresponding properties at finite distances. Such differences can be relevant to gravitational wave astronomy. For example, it has been found that the asymptotic formula relating the curvature perturbation ψ 4 to the gravitational wave strain (Eq. (25)) is invalid for polynomially decaying solutions even at large distances used in numerical wave extraction [41]. Further, we know that the phase of the gravitational wave signal measured at null infinity and at finite distances may differ substantially in EMRIs [42,43].
A clean solution to both the outer boundary and the radiation extraction problems is to include null infinity in the computational domain. No artificial boundary conditions are needed when the solution is calculated globally in space, and the idealized gravitational wave signal can be read off directly at null infinity.
There are two methods to include null infinity in the computational domain: the characteristic method and the hyperboloidal method. These methods correspond to approaching null infinity along null and spacelike directions respectively.
The characteristic method is very well developed [44]. It has been applied recently to the unambiguous computation of gravitational waves from binary black hole mergers simulated via Einstein equations [45][46][47]. Concerning the Teukolsky equation, it has been used in Schwarzschild spacetime for black hole collisions in the close limit [48,49]. The construction of a null foliation in Kerr spacetime, however, is rather complicated [50][51][52][53]. As a result, there is no numerical computation of gravitational perturbations of Kerr spacetime at null infinity.
The second approach to including null infinity in the computational domain is the hyperboloidal approach [54]. The construction of suitable hyperboloidal foliations is simpler than characteristic ones because spacelike surfaces are more flexible than null surfaces [55]. Various studies have demonstrated the accuracy and efficiency of the hyperboloidal method in dealing with perturbations of black hole spacetimes, but most of them are restricted to spherical symmetry. For example, numerical studies have been performed on Yang-Mills equations in a Schwarzschild spacetime [56,57], spherically symmetric wormholes [58] and boson stars [59], self-force computations [60] and gravitational perturbations [41,61] including extreme-mass-ratio inspirals [42,43].
The only numerical computations of black hole perturbations using the hyperboloidal method without spherical symmetry deal with scalar perturbations [62][63][64][65]. One goal of this paper is to present the application of the hyperboloidal method to Teukolsky equations in Kerr spacetime for the calculation of gravitational waveforms at null infinity. The test problem we study to demonstrate the method is of astrophysical interest. The gravitational waveforms we compute are emitted by the perturbations due to a point particle inspiralling into the central Kerr black hole.
Another goal of this paper is to solve a difficulty of applying the hyperboloidal method in Kerr spacetime. The hyperboloidal foliation of Kerr spacetime constructed in [55] uses a transition zone in which a truncated Cauchy surface is matched to a hyperboloidal surface. This technique requires tuning a large number of foliation parameters and results in a sudden drop of characteristic speeds in the transition zone decreasing computational efficiency [41,60,62]. There are two solutions to this problem. One is to use a single smooth surface avoiding the transition zone of [55]. This technique has recently been applied by Rácz and Tóth in a detailed study of polynomial decay rates of a scalar field in Kerr spacetime [64]. They construct the first smooth, horizon-penetrating, hyperboloidal foliation of Kerr spacetime for their study. Modifying the coordinates everywhere smoothly leads to very efficient numerical computations. However, it may be favorable due to practical reasons to employ standard coordinates in the strong field region near the rotating black hole, especially when particles are present in the computational domain, so that the description of the particles' motion does not need to be changed.
In this paper we use the recently proposed hyperboloidal layer ( [66]) instead of a transition zone. The new layer technique has been applied in studies of EMRIs using the Regge-Wheeler-Zerilli formalism completed by EOB-resummed analytical radiation reaction (at linear order in the mass ratio) [42,43]. In this paper, we demonstrate the construction and application of a hyperboloidal layer for the Teukolsky formalism and the kludge approach in Kerr spacetime.
The paper is organized as follows. In Section II we present the standard Teukolsky formalism with the trans-formations typically used in numerical calculations. In Section III we present how the hyperboloidal layer technique is applied to Kerr spacetime and to the Teukolsky equations as two simple coordinate transformations. Numerical results obtained with this formalism are presented in Section IV. After describing briefly the numerical method, we present a comparison of waveforms, an improved agreement of the energy fluxes with frequencydomain computations for several circular-equatorial orbits, a discussion of recoil velocities, and an example simulation that demonstrates remarkable gains in efficiency with the new method. A discussion of our results and an outlook can be found in Section V.

II. THE TEUKOLSKY FORMALISM
Curvature perturbations of Kerr spacetime are governed by a separable equation in the Newman-Penrose formalism [4,5]. The original calculations involve the Boyer-Lindquist representation of the Kerr metric. It is convenient for numerical applications to introduce certain transformations of the Teukolsky equation. We review the common transformations to set up the equations to which the hyperboloidal method will be applied.

A. Teukolsky equation in Boyer-Lindquist coordinates
The Kerr metric in Boyer-Lindquist coordinates (t, r, θ, ϕ) reads where Σ := r 2 + a 2 cos θ 2 , and △ := r 2 + a 2 − 2M r. We denote the mass of the Kerr spacetime by M , its angular momentum by aM . The Teukolsky equation describes curvature perturbations with spin weight s in an adapted Newman-Penrose orthonormal frame. The homogeneous Teukolsky equation in Boyer-Lindquist coordinates reads where D 2 = (r 2 + a 2 ) 2 /∆ − a 2 sin 2 θ. In the following, we review the transformations that are typically applied to the Teukolsky equation for numerical computations [13,67,68].

B. Transformations of the Teukolsky equation
Boyer-Lindquist coordinates have undesirable features for numerical computations. Their time slices intersect at the bifurcation sphere leading to a singular radial coordinate at the horizon, in a similar way as Schwarzschild time slices do. There are two ways to deal with the coordinate singularity at the horizon: transform to a horizonpenetrating time slicing, or push the horizon to infinite coordinate distance with the tortoise coordinate.
Horizon-penetrating coordinates are regular at the event horizon. They are typically used in combination with excision. The numerical treatment of the inner boundary is clean with this method because there are no characteristics coming out of the horizon. In addition, the radiation absorbed by the black hole can be calculated accurately.
The most common approach in perturbation theory, however, is to use the tortoise coordinate. The tortoise coordinate pushes the horizon to infinite coordinate distance. The infinite domain is then truncated at a finite distance, where artificial ingoing boundary conditions are imposed. These lead to a contamination of the physical solution due to numerical reflections from the inner boundary. But these reflections are small because the potential terms in the Teukolsky equation fall off exponentially fast in the tortoise coordinate, so the ingoing boundary conditions are quite accurate. In our study, we focus on the outgoing radiation and accept small reflections from the inner boundary. The truncation of the infinite domain at a finite coordinate implies that the efficiency of the code is not optimal because a large region in negative tortoise coordinate needs to be included in the computational domain.
The relation of the tortoise coordinate to horizonpenetrating coordinates near the event horizon is similar to the relation of standard Cauchy coordinates to hyperboloidal coordinates near null infinity [41,69]. A fundamental difference, however, is that the structure of the solution in the strong field region can be interesting, whereas the asymptotic solution is essentially simple in asymptotically flat spacetimes. The additional resolution that the tortoise coordinate provides near the black hole can be relevant in certain studies. This question should be analyzed in more detail. We choose to employ the tortoise coordinate, and leave the study of horizonpenetrating coordinates to future work.
The angular coordinate also needs to be transformed to cure the infinite twisting near the horizon [5,67]. Further, a rescaling respecting the fall-off behavior of the curvature perturbations is necessary. Finally, we separate the azimuthal dependence to arrive at a 2 + 1 dimensional system. These transformations are listed as follows: 1. Introduce the tortoise coordinate r * , 2. Introduce the azimuthal coordinateφ, 3. Rescale the curvature perturbation of spin weight s according to its asymptotic fall-off behavior, 4. Separate the azimuthal dependence using the mode number k, In the following, we drop the subscript k from ψ k for conciseness of notation. The Teukolsky equation becomes [68] These are standard transformations for numerical com-putations. Hyperboloidal compactification adds only two additional transformations to this list: a transformation of the time coordinate, and a compactification of the radial coordinate.

III. A HYPERBOLOIDAL LAYER FOR THE TEUKOLSKY EQUATION
We present the hyperboloidal layer method in the language of standard transformations to make its application as straightforward as possible for future studies.

A. Hyperboloidal compactification
A hyperboloidal surface is an everywhere spacelike surface that approaches null infinity [54]. It has been known already in the early days of numerical relativity that such surfaces are beneficial for studying gravitational radiation [70][71][72][73]. A systematic study of the hyperboloidal initial value problem started with Friedrich's work showing that future null infinity is smooth for certain nontrivial classes of initial data [74]. There has been substantial effort to numerically simulate such data using Friedrich's conformally regular field equations [75,76].
Even though the benefits of hyperboloidal surfaces for radiation problems beyond the Einstein equations have been clear, it took more than two decades before the hyperboloidal method could be applied successfully to systems other than the conformally regular field equations. The first application of this kind is the study of magnetic monopoles by Fodor and Rácz [77,78]. Their hyperboloidal foliation is based on a scri-fixing gauge in Minkowski spacetime first explicitly constructed by Moncrief [79] (see also [76]). It was suggested that scri-fixing gauges should be beneficial also for studies on black hole spacetimes, however, the construction of a good hyperboloidal coordinates turned out to be difficult [80][81][82][83][84][85].
The general construction of suitable hyperboloidal, scri-fixing gauges has been presented in [55]. We follow [55] to present the hyperboloidal compactification. The technique consists of introducing a compactifying radial coordinate ρ and a suitable time coordinate τ .

Radial compactification
We map the infinite positive domain in the tortoise coordinate r * to a finite domain in the compactifying coordinate ρ to compute the solution in an unbounded physical domain using a finite numerical grid. After compactification, the outer boundary of the numerical grid corresponds to infinity with respect to the physical coordinate. Infinity has a rich structure in spacetime. What part of infinity is included in the computational domain with radial compactification depends on the nature of the time surfaces, as discussed in the next subsection.
The radial compactification can be performed conveniently with the following definition of a compactifying coordinate: The choice of Ω determines the properties of the compactification. Its zero set corresponds to infinity with respect to r * . For example, a common choice encountered in the literature is Ω = 1 − ρ. We need more freedom in the choice of compactification, therefore we will state the general properties that the function Ω needs to satisfy. These properties are those that are satisfied by a conformal factor in the Penrose conformal compactification picture [86,87]. For all finite r * we require Ω > 0. Denoting the zero set of Ω as S (the ρ-coordinate location of r * -infinity), we require where Ω ′ ≡ dΩ/dρ. It is known that compactifying the radial coordinate along Cauchy surfaces results in resolution problems [88]. These problems can be avoided by a hyperboloidal time transformation that changes the compactification at spatial infinity to a compactification at null infinity.

Time transformation
A suitable time transformation shall leave the exterior timelike Killing field invariant so that no gauge dynamics is introduced into the stationary background. The new time coordinate τ satisfies then ∂ τ = ∂ t , which is achieved by a transformation of the form τ = t − h(x i ), where x i are the spatial coordinates. The function h is called the height function. Most useful foliations of stationary spacetimes are given via such a relation. We are interested in radially outgoing gravitational waves, therefore we choose the height function to depend only on the tortoise coordinate r * The height function is chosen such that the surfaces τ = const. are hyperboloidal [55]. The specific choice we use is presented in the next section (see Eq. (7)). In summary, hyperboloidal compactification consists of two simple transformations that can be written in the form (4) and (5) with free functions Ω and h. The specific choice of these free functions depends on the background spacetime and its coordinatization. There is a large freedom in their choice which allows us to use standard coordinates in an interior domain as described below.

B. Hyperboloidal layer for Kerr spacetime
An essential advantage of the hyperboloidal method is that it requires modifications only in the asymptotic domain, which allows us to use arbitrary coordinates in the strong field regime.
A transition zone is used in [55] for matching arbitrary coordinates near the central black hole to hyperboloidal scri-fixing coordinates in the asymptotic domain. Subsequent numerical experiments showed that the transition zone requires the adjustment of many arbitrary parameters and leads to a sudden drop of outgoing characteristic speeds decreasing the efficiency of the hyperboloidal method [41,60,62].
A new technique that avoids these deficiencies is the hyperboloidal layer presented in [66]. The idea is to attach a radial shell to the truncated computational domain in which the foliation is hyperboloidal and the radial coordinate is compactified. The layer needs to be attached in a sufficiently smooth way for the stability of the numerical computation.
We would like to have the new coordinate ρ agree with the tortoise coordinate in the strong field region. Denoting the location of the interface between the interior domain and the layer by R * , we set Ω = 1 for ρ < R * so that ρ ≡ r * in the interior. We also require that the coordinates agree up to a certain order at the interface. These conditions are fulfilled by the following expression [43,66] where Θ denotes the step function and S denotes the coordinate location of the outer boundary. We emphasize that various similar choices for the radial compactification are possible. The original construction of a height function for a hyperboloidal layer [66] requires unit characteristic speed in the layer with respect to the layer coordinates τ and ρ. An equivalent but more intuitive requirement is that outgoing null rays have the same representation in the layer coordinates as in the interior coordinates [43]. In Schwarzschild spacetime with the tortoise coordinate the requirement reads τ − ρ = t − r * . From this relation together with (5) we get In the domain where Ω = 1, we have h = 0 so that we obtain the standard coordinates. For the r * derivative of the height function we get We refer to H as the boost function. It vanishes in the interior domain where standard coordinates are used with Ω = 1; it is unity at infinity, where Ω = 0. The Refs. [43,66] use the prescription (8) in Minkowski spacetime with standard coordinates, and in Schwarzschild spacetime with the tortoise coordinate. In this paper we apply it for the Kerr spacetime in Boyer-Lindquist tortoise coordinates. This procedure leads to a regular hyperboloidal compactification in Kerr spacetime, even though the outgoing null rays do not have the simple form t − r * , because the Kerr metric in the tortoise Boyer-Lindquist coordinates has asymptotically the same form as the Schwarzschild metric in tortoise Schwarzschild coordinates, or the Minkowski metric in standard coordinates. This feature is another evidence for the flexibility of hyperboloidal coordinates. Only the asymptotic behavior of the null cone is relevant for the regularity of hyperboloidal compactification, whereas in the characteristic method the exact local form of null cone plays an essential role, thereby restricting the coordinate transformations [50][51][52][53].
The choices (4) and (5) with the free functions explicitly given in (6) and (7) 6. Introduce a time coordinate τ , The compactification function Ω(ρ) and the height function h(r * ) are determined via (6) and (7). With the coordinate transformations at hand, we can now proceed to transform the Teukolsky equation. It is a priori not clear that a simple coordinate transformation will lead to a regular hyperboloidal compactification of the Teukolsky equations. There are various possibilities to study the regularity of such compactification. One is to derive conformal Teukolsky equations with respect to a conformally extended, regular Kerr metric. It has been shown in a previous application of hyperboloidal compactification using the Teukolsky formalism in Schwarzschild spacetime (called the Bardeen-Press equation after [3]) that such an analysis is not necessary [61]. However, the method of [61] still requires the calculation of the Teukolsky equations ab initio in a general orthonormal Newman-Penrose frame which is then adapted to a hyperboloidal scri-fixing gauge to ensure regularity of hyperboloidal compactification. The method we use in this paper is easier to apply.
A general operator for an equation in the independent variables {t, r * , θ}, such as the Teukolsky equation (3), can be written as O[ψ] = 0, where O is the differential operator The coefficients depend on r * , and θ; they can be read off from (3). With the transformations (9) and (10) the derivative operators with respect to the old coordinates can be written in terms of the new coordinates as We used (8) in the expression for ∂ r * . We write the transformed operator as where and H ′ := dH/dρ. All coefficients are functions of ρ, and θ.
The transformed system is by construction equivalent to the original system (3) in the interior domain ρ < R * where we have H = 0. We need to make sure that the hyperboloidal compactification is regular at infinity where H = 1. Note that asymptotically We see by inspection that the coefficients of the Teukolsky equation (3) have the correct asymptotic behavior ensuring regularity of the hyperboloidal compactification. The transformed coefficiencts have explicitly finite values at future null infinity. We mention that this feature is not special to the Teukolsky equations. Similar regular hyperboloidal compactifications for other wave equations have been studied as mentioned in the Introduction.
We also need to make sure that the transformed system is pure outflow at the outer boundary so that no boundary conditions are required. This condition, along with the regularity, can be checked by evaluating the transformed equations at infinity, that is, at {ρ = S}. We obtain for s = −2 D∂ 2 τ ψ = −S(S − R * )∂ τ ∂ ρ ψ + 2∂ 2 θ ψ −4(4M + ia(k − 2 cos θ))∂ τ ψ +2 cot θ∂ θ ψ − 2(3 + k 2 − 4k cos θ + cos(2θ))ψ, whereD = S(S − R * ) − 2a 2 sin 2 θ. The principal part of the operator is of the form where C 1 and C 2 are coefficients. Therefore, the modes of the system propagate either along the boundary at {ρ = S}, or out of the domain along the level sets of characteristics τ − ρ. As expected, there are no incoming modes.

IV. NUMERICAL RESULTS
In this section, we present tests using the hyperboloidal compactification of the Teukolsky equation. Our main result is the efficiency and accuracy of the hyperboloidal layer method providing unambiguous waveforms at future null infinity from large-and extreme-mass-ratio inspirals in Kerr spacetime.

A. The numerical method
For the discretization of the continuous equation we use a two-step Lax-Wendroff algorithm as in [67]. After performing the transformations presented in the previous sections, we rewrite the vacuum equation in the form where the coefficients with a tilde are obtained by dividing the coefficients of the operator in (11) by −A τ τ . We put the equation (12) in first-order form (in ρ and τ ) by defining a new field variable We chose to use this first-order form of the equation for our numerical evolutions, because it has been discovered through extensive experimentation [67] that such a form is ideally suited for long stable evolutions. Then equation (12) with source terms T takes the form where is the solution vector. The subscripts R and I refer to the real and imaginary parts respectively. The matrices M , A and L are obtained from (12) as in [13]. Here it will suffice to simply indicate the final form taken by these matrices: The angular derivatives are encoded in L.
The main advantage of casting the equation (12) into the form (15) is that the system has advantageous properties in the variable ρ. The matrix M has a complete set of linearly independent eigenvectors with real eigenvalues. This is not a rigorous statement on the hyperbolicity of the system because the matrix L contains second-order angular derivatives. Nevertheless, experiments show that the system is numerically well-behaved. A time-explicit evolution scheme is developed for this system of linear partial differential equations using the two-step, secondorder Lax-Wendroff finite-difference method. We write (15) as where Each iteration consists of two steps. In the first step, the solution vector between grid points is obtained from This is used to compute the solution vector at the next time step, The angular subscripts are dropped in the above equation for conciseness. All angular derivatives are computed using second-order, centered finite difference expressions. Symmetries of the spheroidal harmonics are used to determine the angular boundary conditions: For even |m| modes we have ∂ θ ψ = 0 at θ = 0, π; for odd |m| modes we have ψ = 0 at θ = 0, π. We impose "ingoing" boundary conditions at the inner radial boundary. We do not need to apply any boundary condition at the outer boundary.
The small compact object inspiraling into the central supermassive black hole is modeled as a point particle in the large-or extreme-mass-ratio limit. The representation of point particle in mathematics is performed typically by a Dirac delta distribution. More specifically, the particle energy-momentum tensor takes the form where u α and µ refer to the 4-velocity and the rest mass of the particle, respectively. Note thatτ ≡ dτ /ds appears in the denominator of this expression, which tends to ∞ as the particle approaches the horizon (we denote by s the proper time along the particle's trajectory). Thus, the particle source-term on the right-hand-side of (15) smoothly decays away as the particle approaches the horizon, thereby allowing the evolution equation to gradually transition into its homogeneous form. This source-term is constructed using the energy-momentum tensor describing a point particle moving in Kerr spacetime depicted above. The explicit expression for T is very lengthy and not particularly illuminating. Here, we simply point out that the final expression is built using Dirac delta functions in r and θ, as well as first and second derivatives of the delta functions in these variables. These terms have coefficients that are complex functions of the black hole's physical parameters and also the trajectory of the point particle. Details on the particle source-term, along with the representation of the delta function and its derivatives on a numerical grid, are given in [13][14][15]89]. The trajectory of the particle in a decaying orbit around the central black hole is computed separately and is then used to calculate the source-term mentioned above, which in turn is fed into the Teukolsky equation solver code. The computation of the decaying trajectory can be broken into three distinct pieces: An early time adiabatic inspiral, in which the particle is approximated as evolving through a sequence of bound orbits of the central black hole that are calculated using an "energy balance" approach; a late time geodesic plunge, in which the particle falls into the black hole and radiation reaction is ignored; and an intermediate transition regime that smoothly connects these two pieces. Details on how these steps are handled and used by the Teukolsky equation solver code can be found in [15]. It is straightforward to make use of decaying orbital trajectories from other approaches, such as the EOB or the self-consistent selfforce approaches.

B. Waveforms
After solving the Teukolsky equation with a particle source in the time-domain, we compute the gravitational waveforms h + and h × using the simple relationship shown below that is valid in the far field, ψ → r 2 Note that in our study we have direct access to the far field because null infinity is part of the computational domain. This allows us to use the above relation as an equality and to extract the waveform cleanly. Figures 1 and 2 depict the gravitational waveform h ℓm measured at future null infinity, with ℓ = m = 2. The mass-ratio for this evolved binary system is µ = 10 −3 and the spin of the central Kerr black hole is a = 0.5. The system undergoes a decay due to gravitational wave emission over approximately 123 full cycles. The computational domain is ρ ∈ [−50, 50] and θ ∈ [0, π] with grid sizes of 3125 and 32 respectively. The hyperboloidal layer starts at ρ = 14, therefore the particle's orbit never crosses into the layer. There are a total number of 562, 500 time-steps in this computation.
The figures depict the waveforms from both the original Cauchy code [13,89] and the new code with the hyperboloidal layer in Kerr spacetime developed in this work. The Cauchy code approaches spatial infinity, and therefore cannot provide direct access to waveforms at null infinity. To obtain data from the Cauchy code to compare with that from the hyperboloidal layer code, we extract the waveforms at multiple radii on the spatial grid and extrapolate them along fixed values of retarded time, to infinity, using a simple second order polynomial expansion. We do not employ higher order extrapolation for the Cauchy code because the truncation error becomes larger than extrapolation error with the given resolution and second order finite differencing. Nevertheless, the waveforms from the two codes agree so well that it is difficult to make out any difference between the two in these figures. In Figure 3 we show the differences in the waveform amplitude and phase for the same data. The relative difference between the amplitudes is at the level of 10 −3 and the maximum difference in phase is 0.06 rad. A better agreement between extrapolated and null infinity waveforms can be obtained with more accurate extrapolation algorithms and higher-order methods in spherical symmetry [43].

C. Fluxes
Gravitational waves carry energy away from a binary system thus causing it to inspiral. The calculation of the energy carried away by these waves uses the standard expression, in a postprocessing step after time evolution. In standard computations, the limit to infinite radius in the above equation is approximated by application of the formula at finite but large radii, and a subsequent extrapolation of the corresponding values to infinity. One major advantage of the hyperboloidal method is that the limit in Eq. (26) can be replaced by local evaluation of the integral along future null infinity, which is the boundary of the computational domain. Therefore, the extrapolation of fluxes is not necessary, and the energy can be computed cleanly. There are more subtleties in the computation of the fluxes related to the integration in time from the infinite past [90]. We ignore these and set the integration constant to zero. We observe empirically that this procedure does not lead to large errors.
In Table I we show numerical values of the energy fluxes computed at null infinity by the new hyperboloidal layer code for several circular-equatorial orbits. Comparing these values with those from very high-accuracy frequency-domain approaches [91], we note an agreement at the 99.95% level. We interpret this high level of agree- ment with the frequency-domain fluxes as strong evidence for the accuracy of the code, especially as compared with the original time-domain Cauchy code that achieves agreement with frequency-domain fluxes at the level of 99% [13,89].
In Table I, we also show the standard three level convergence rates. We obtain clear second-order convergence, as expected. The fluxes depicted in the table were obtained via Richardson extrapolation using the data from the different grid resolutions. The computational domain is ρ ∈ [−50, 50] and θ ∈ [0, π] with grid sizes 1250 × 32, 2500 × 64 and 5000 × 128. The hyperboloidal layer starts at ρ = 14, therefore these orbits never cross into that layer throughout the simulations.

D. Kick velocities
The asymmetric emission of gravitational waves during the inspiral of the small black hole into the supermassive central black hole leads to a recoil due to conservation of momentum. We compute this recoil by first computing the linear momentum flux carried away by the gravitation waves, and then performing a simple time-integral of this flux [15]. To compute the momentum flux, we use the expression where the n i is the unit direction vector in spherical coordinates. As for Eq. (26), the limit to infinity is replaced by local evaluation of the integral at future null infinity, which allows us to read off the momentum flux without an approximation in the radial direction. The integration constant resulting from starting the integration at τ = 0 is set arbitrarily to zero. In Figures 4, 5 and 6 we depict the results from the computation of linear momentum flux and its timeintegral. The mass-ratio for this computation is µ = 10 −3 with the central black hole spin a = 0.7. The computational domain is ρ ∈ [−50, 50] and θ ∈ [0, π] with grid size 3125 × 32. The hyperboloidal layer starts at ρ = 14, therefore the particle's orbit never crosses into that layer throughout the simulation. There are 2, 250, 000 timesteps in this computation. Once again we superimpose results from the Cauchy code and the new hyperboloidal code which agree to a very high level. The final values of the recoil speed from the two different codes, 1.645×10 −8 and 1.641 × 10 −8 , agree to three significant digits, consistent with our previous results.

E. Efficiency
Previous sections provided evidence for the accuracy of the hyperboloidal method. The gain in accuracy in observables, such as the waveforms or energy fluxes, is mainly due to the direct numerical access to asymptotic quantities at future null infinity, which is part of the computational domain. The difference in those quantities, however, is rather small, essentially because extrapolation in given background spacetimes already gives reliable results. The basic advantage of the hyperboloidal method in this context is that it simplifies the accurate extraction of physical observables from numerical simulations: It eliminates the extrapolations.
The remarkable feature of the hyperboloidal method is that it provides improved accuracy at negative cost. The numerical computation is much more efficient along hyperboloidal surfaces than along Cauchy surfaces. One may think that the calculation of waveforms at null infinity introduces additional computational expense for numerical simulations. In fact, the opposite is true: One gains more accuracy for less computational cost.
The explanation of this property of hyperboloidal evolutions is that the standard method of using truncated Cauchy foliations with artificial outer boundary conditions is unnatural to study radiative solutions. An efficient computation of the asymptotic gravitational wave signal can only be achieved if the surfaces, along which the solution is computed, follow the outgoing signal closely. Gravitational waves propagate to infinity along null rays. Null rays and Cauchy surfaces diverge infinitely from each other in the asymptotic region. Cauchy surfaces are therefore inadequate for computing outgoing asymptotic radiation. Further, their truncation at a finite distance introduces errors through artificial outer boundary data contaminating the computation.
To avoid such contamination in a standard Cauchy code the outer boundary is usually chosen to be causally disconnected from the region of interest. This implies that the duration of the simulation is limited by the size of the numerical grid. When null infinity is part of the computational domain, however, no such restriction applies. Therefore, in principle, we can perform arbitrary long evolutions, only restricted by the accumulation of the numerical truncation and the physical approximation errors. The efficiency of a hyperboloidal code as compared to a Cauchy-type code increases with the duration of the simulation.
To demonstrate the efficiency of our hyperboloidal Teukolsky code, we performed an evolution of an EMRI that lasts for one million M . The spin of the central black hole is 0.9, the mass ratio is 10 −4 . There are about 10, 111 orbital cycles in this simulation. The computational domain is ρ ∈ [−50, 50] and θ ∈ [0, π] with the grid size 5000 × 32. The hyperboloidal layer starts at ρ = 14, therefore the particle's orbit never crosses into that layer throughout the simulation. There are over a 100 million time-steps involved in this computation. In Fig. 7 we show the amplitude and frequency against time in a halflogarithmic scale. The amplitude and frequency of the gravitational waveform at null infinity is fairly constant over most of the simulation, as expected. The main content of this plot is to demonstrate that such simulations in time-domain are now possible with the new method. Using the Cauchy code for such a simulation, would have required us to place the outer boundary at least as far as 500, 000M . So the Cauchy code would take 5000 times longer to perform this simulation with equivalent local grid resolution. This number can be increased even further using horizon-penetrating coordinates.
A million M in geometric units for a supermassive black hole of, say, 6 million solar masses corresponds to about one Earth year in real time. The simulation for Fig. 7 has been performed on a computer cluster using on the order of 1000 processor-cores in just under a day. This shows that we can perform long-time simulations of EMRIs much faster than the expected real time duration of the events, and that large parameter studies in perturbation theory are computationally feasible with time-domain codes.

V. DISCUSSION AND OUTLOOK
We presented the first numerical evolution scheme for computing gravitational perturbations of a rotating black hole spacetime including future null infinity. We tested our method on the inspiral of a point particle into a central rotating black hole. The main advantages of the method are: (i) increase in accuracy by direct radiation extraction at future null infinity; (ii) a clean solution to the outer boundary problem in Kerr spacetime; and (iii) a great gain in efficiency for long time simulations.
Our method is based on the hyperboloidal compactification of Kerr spacetime presented in [55]. We also showed, for the first time, how the hyperboloidal layer developed in [66] can be employed in a Kerr spacetime. The layer technique allowes us to use arbitrary coordinates in an interior compact domain including the central black hole and the orbit of the particle, while compactifying the asymptotic domain along a hyperboloidal foliation. Because the hyperboloidal layer is attached to the interior domain sufficiently smoothly, only minor modifications of existing computational infrastructure is needed (see [43] for a related study in Schwarzschild spacetime) We showed that the hyperboloidal compactification of the Teukolsky equations can be performed by a simple coordinate transformation. There is no need for studying a conformally regular Kerr spacetime, or for calculating the equations ab initio in a general orthonormal Newman-Penrose frame as has been done in [61]. The coordinate transformation of the equations as given in Boyer-Lindquist coordinates leads naturally to a regular hyperboloidal compactification. Only the asymptotic behavior of the null cone is relevant for the hyperboloidal method, which implies that the simple prescription devised in Minkowski and Schwarzschild spacetimes can be used in Kerr spacetime with respect to Boyer-Lindquist tortoise coordinates. This provides evidence for the flexibility of the hyperboloidal approach, as opposed to the characteristic approach which restricts the coordinates locally for compactification at future null infinity [50][51][52][53].
We compared our results to a previous study of EM-RIs in Kerr spacetime [15]. The emphasis in this paper is not on new physics but on the numerical accuracy and efficiency of the hyperboloidal method as applied in Kerr spacetime. The waveforms at null infinity differ only slightly from their previous computation obtained by extrapolation. We showed that the hyperboloidal method gives better accuracy in the energy fluxes at infinity. Our strongest case for efficiency is the simulation of an inspiral that lasts one million M in geometric units, which would have taken about 5000 times longer using the standard method.
In the future this method should be applied further to physically interesting problems. One could, for example, look for the new ringdown frequency mode for perturbations of a Kerr black hole discovered by Mino and Brink [92] (see also [93]) or for transient resonances expected for certain orbits in Kerr spacetime [94]. Such studies require high accuracy, so the description of the particle's motion needs to be improved, for example by including conservative effects in the kludge approach, by reading off the source terms from an EOB description, or by computing the self-consistent evolution driven by the high-order in mass ratio self-force of the particle.
It is clear that there is much room for development in the accurate computations of EMRIs. It would be interesting to compare the different approaches to the EMRI problem in the same setting, similar to comparative studies performed in nonlinear numerical relativity [95]. Comparison of physical invariants on given background spacetimes is simpler, therefore such a study should be easier to perform in the perturbative setting than in the nonlinear one.
On a numerical level, aside from the technical improvements such as the use of high-order finite differencing or multidomain pseudospectral methods, implementation of horizon-penetrating coordinates near the black hole should lead to a further improvement in efficiency. Horizon-penetrating, hyperbolodial surfaces have certain advantages in perturbation theory over the standard Boyer-Lindquist or Schwarzschild time surfaces that intersect at the bifurcation sphere and approach spatial infinity [69]. In our context, horizon-penetrating coordinates should allow us to calculate accurately the absorbed fluxes at the horizon via the ingoing radiation represented by the curvature component ψ 0 .
Our calculation of an inspiral that corresponds to roughly a year in real time (depending on the mass of the central supermassive black hole) but takes only about one day of computation is a strong evidence for the efficiency of our numerical infrastructure (see Section IV E). Note that the situation is reversed when fully nonlinear Einstein equations are solved for equal-mass, stellar, binary black-hole mergers. There, a typical computation takes about a month, but a stellar black-hole merger takes only milliseconds in real time. Therefore, it is imperative that the efficiency of nonlinear Einstein codes is increased. Given the physical simplicity of an equal-mass binary black-hole configuration, a combination of technical improvements including the hyperboloidal method, along with analytical insight into the problem, should allow us to simulate such systems much more efficiently in the near future.
The application of the hyperboloidal method to generic computations with the fully nonlinear Einstein equations is an outstanding problem. Recently there has been studies on various aspects of the hyperboloidal approach such as initial data [96], or gauge conditions [97]. There are also suggestions on the evolution problem that do not require explicit regularity of the equations at future null infinity [98][99][100]. The only successful numerical implementation of such a formalism is by Rinne in axisymmetry [101], but no generic computation could be presented so far. We hope that the insight provided by studies in perturbation theory will help the application of the hy-perboloidal method to nonlinear Einstein equations.