Time-domain metric reconstruction for self-force applications

TL;DR

This paper introduces a new, computationally efficient time-domain method for calculating the gravitational self-force in Kerr spacetime, enabling more flexible and stable orbit evolution simulations.

Contribution

The paper presents a novel time-domain metric reconstruction approach from curvature scalars for GSF calculations, improving efficiency and stability over existing methods.

Findings

Method successfully applied to circular orbits in Schwarzschild geometry.

Reduces computational cost compared to traditional time-domain approaches.

Allows for flexible orbit types and self-consistent orbit evolution.

Abstract

We present a new method for calculation of the gravitational self-force (GSF) in Kerr geometry, based on a time-domain reconstruction of the metric perturbation from curvature scalars. In this approach, the GSF is computed directly from a certain scalar-like self-potential that satisfies the time-domain Teukolsky equation on the Kerr background. The approach is computationally much cheaper than existing time-domain methods, which rely on a direct integration of the linearized Einstein's equations and are impaired by mode instabilities. At the same time, it retains the utility and flexibility of a time-domain treatment, allowing calculations for any type of orbit (including highly eccentric or unbound ones) and the possibility of self-consistently evolving the orbit under the effect of the GSF. Here we formulate our method, and present a first numerical application, for circular geodesic orbits in Schwarzschild geometry. We discuss further applications.

Figures

Equations

$$h_{\alpha \beta }=h_{\alpha \beta }^{\mathrm{rec}}+h_{\alpha \beta }^{\mathrm{comp}}+h_{\alpha \beta }^{\mathrm{gauge}},$$

$$l^{\alpha }{h}_{\alpha \beta }^{\mathrm{rec}}=0\text{(IRG)},n^{\alpha }{h}_{\alpha \beta }^{\mathrm{rec}}=0\text{(ORG)},$$

$$g_K^{\alpha \beta }{h}_{\alpha \beta }^{\mathrm{rec}}=0,$$

$$e_1^{\alpha }=:l^{\alpha }$$

$$e_2^{\alpha }=:n^{\alpha }$$

$$e_3^{\alpha }=:m^{\alpha }$$

$$e_4^{\alpha }=:{\bar{m}}^{\alpha }$$

$$\Delta$$

$$\Sigma$$

$$\varrho$$

$$\left(\frac{(r^2+a^2{)}^2}{\Delta }-a^2{\sin }^2\theta \right)\frac{{\partial }^2\Phi }{\partial t^2}+\frac{4Mar}{\Delta }\frac{{\partial }^2\Phi }{\partial t\partial \phi }$$

$$+\left(\frac{a^2}{\Delta }-\frac{1}{{\sin }^2\theta }\right)\frac{{\partial }^2\Phi }{\partial {\phi }^2}-{\Delta }^{-s}\frac{\partial }{\partial r}\left({\Delta }^{s+1}\frac{\partial \Phi }{\partial r}\right)$$

$$-\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial \Phi }{\partial \theta }\right)+(s^2{\cot }^2\theta -s)\Phi$$

$$-2s\left(\frac{M(r^2-a^2)}{\Delta }-r-ia\cos \theta \right)\frac{\partial \Phi }{\partial t}$$

$$-2s\left(\frac{a(r-M)}{\Delta }+\frac{i\cos \theta }{{\sin }^2\theta }\right)\frac{\partial \Phi }{\partial \phi }=0,$$

$${\mathbf{D}}_l^4{\bar{\Phi }}^{\mathrm{IRG}}$$

$${\Delta }^2{\tilde{\mathbf{D}}}_n^4{\Delta }^2{\bar{\Phi }}^{\mathrm{ORG}}$$

$${\mathbf{D}}_l$$

$${\tilde{\mathbf{D}}}_n$$

$${\tilde{\mathcal{L}}}_{-1}{\tilde{\mathcal{L}}}_0{\tilde{\mathcal{L}}}_1{\tilde{\mathcal{L}}}_2{\bar{\Phi }}^{\mathrm{IRG}}-12M{\Phi }_{,t}^{\mathrm{IRG}}$$

$${\mathcal{L}}_1{\mathcal{L}}_0{\mathcal{L}}_{-1}{\mathcal{L}}_{-2}{\bar{\Phi }}^{\mathrm{ORG}}+12M{\Phi }_{,t}^{\mathrm{ORG}}$$

$${\mathcal{L}}_s$$

$${\tilde{\mathcal{L}}}_s$$

$$h_{\alpha \beta }^{\mathrm{rec}}=\mathrm{Re}\left(e_{\mathbf{a}(\alpha }{e}_{\mathbf{b}\beta )}{\mathcal{D}}^{\mathbf{ab}}\Phi \right),$$

$${\mathcal{D}}^{11}$$

$${\mathcal{D}}^{33}$$

$${\mathcal{D}}^{13}$$

$${\mathcal{D}}^{22}$$

$${\mathcal{D}}^{44}$$

$${\mathcal{D}}^{24}$$

$${\mathbf{D}}_n$$

$${\mathbf{D}}_m$$

$${\mathbf{D}}_{\bar{m}}$$

$$\beta$$

$$\mu$$

$$\tau$$

$$h_{\alpha \beta }^{\pm }=h_{\alpha \beta }^{\mathrm{rec}\pm }+h_{\alpha \beta }^{\mathrm{comp}\pm }+h_{\alpha \beta }^{\mathrm{gauge}\pm },$$

$$F_{\alpha }^{\pm }:=-\frac{1}{2}\mu ({\delta }_{\alpha }^{\beta }+u_{\alpha }{u}^{\beta })(2{\nabla }_{\nu }{h}_{\beta \mu }^{\pm }-{\nabla }_{\beta }{h}_{\mu \nu }^{\pm })u^{\mu }{u}^{\nu },$$

$$F_{\alpha }^{\ell \pm }:=\underset{r\to r_{\mathrm{p}}^{\pm }}{\lim }\sum _{m=-\ell }^{\ell }{Y}_{\ell m}({\theta }_{\mathrm{p}},{\phi }_{\mathrm{p}})\int F_{\alpha }^{\pm }{\bar{Y}}_{\ell m}(\theta ,\phi )d\Omega ,$$

$$F_{\alpha }=\sum _{\ell =0}^{\infty }\left[\frac{1}{2}(F_{\alpha }^{\ell +}+F_{\alpha }^{\ell -})-B_{\alpha }\right],$$

$$h_{\alpha \beta }^{\mathrm{rec}\pm }=\mathrm{Re}\left(e_{\mathbf{a}(\alpha }{e}_{\mathbf{b}\beta )}{\mathcal{D}}^{\mathbf{ab}}{\Phi }^{\pm }\right).$$

Full text

I Introduction

The extreme-mass-ratio regime of the binary black hole inspiral problem is most naturally tackled using the gravitational self-force (GSF) approach, which is based on a systematic expansion of the field equations in the small mass ratio (usually without any weak-field or slow-motion assumptions). At leading order, one has a pointlike particle moving in a geodesic orbit around the large black hole. At 1-GSF order, the particle experiences an effective self-force due to interaction with its own linear gravitational perturbation; nonlinear self-interaction effects are accounted for at subsequent orders. There is now a rigorously established formulation of self-forced motion through 2-GSF order, for general vacuum curved spacetimes (at least for nonspinning particles) Gralla and Wald (2008); Gralla (2012); Poisson et al. (2011); Pound (2012). There also exist a number of methods that recast the formal equation of motion into practical regularization schemes in the case of orbits around a Kerr black hole—cardinal among these are mode-sum regularization Barack and Ori (2000) and the puncture method Barack et al. (2007); Vega and Detweiler (2008); see Barack (2009); Wardell and Gopakumar (2015) for reviews. Finally, there now exist a variety of computational strategies and working codes that implement the regularization schemes numerically and calculate the GSF for orbits in Schwarzschild or Kerr geometries (see, e.g., van de Meent (2016) and references therein). Actual calculations have so far been restricted to 1-GSF order, but second-order results are expected soon Pound and Miller (2014); Pound (2014, 2015); Miller et al. (2016). This computational program is strongly motivated by the prospect of observing extreme-mass-ratio inspirals with planned gravitational-wave observatories in the millihertz band.

All GSF calculations require information about the local metric perturbation near the particle. Existing schemes differ on the precise type of information required and on how it is obtained in practice. More concretely, schemes may be classified according to the regularization method they rely on (e.g., mode sum vs. puncture), the version of field equations that are being solved (e.g., linearized Einstein’s equation vs. Teukolsky’s equation), or, relatedly, the gauge in which the perturbation is computed. They can also be categorized according to whether the perturbation is solved for in the frequency domain (FD) or in the time domain (TD). In FD methods one solves a single (or a coupled set of) ordinary differential equations (ODEs) for each frequency-harmonic mode of the perturbation. In TD methods one instead directly solves the partial differential equations (PDEs) that govern the time evolution of the perturbation, or of its individual multipole or azimuthal modes.

Two main computational approaches have been responsible for much of the progress in the field so far. The first, and more direct, is based on numerically solving the linearized Einstein’s equations in the Lorenz gauge. This can be (and has been) done in conjunction with either mode-sum or puncture regularizations, and is most naturally implemented in the TD. (A FD version of this approach has been developed for calculations in Schwarzschild spacetime Akcay et al. (2013), but extending it to Kerr is problematic due to the impossibility of separating the Lorenz-gauge perturbation equations into frequency-harmonic modes in that case.) The Lorenz-gauge/TD approach offers great versatility: codes can be readily implemented for any type of orbit, including highly eccentric or unbound ones, and provide a natural framework for studying the orbital evolution under the GSF effect. But the method has two significant disadvantages. First, it is computationally very expensive: It involves the numerical evolution of a set of 10 coupled PDEs for each multipole or azimuthal mode of the perturbation, over sufficiently long time to ensure that spurious radiation from imperfect initial conditions is sufficiently suppressed. Second, numerical evolutions have been shown Dolan and Barack (2013) to be contaminated by certain nonphysical Lorenz-gauge modes that grow linearly in time, the elimination of which remains an open problem.

The second, less direct computational approach is based on metric reconstruction from curvature scalars. Instead of tackling the complicated set of metric perturbation equations, one works within the elegant framework of the Newman-Penrose formalism, with the numerical task now reduced to solving the Teukolsky equation for either of the two Weyl scalars $\Psi_{0}$ or $\Psi_{4}$. Since this equation is fully separable into frequency-harmonic modes, even in the Kerr case, the problem further reduces to solving a set of ODEs in the FD. This approach has been gaining much momentum in the past few years, with the formulation of a metric reconstruction procedure suitable for GSF calculations Keidl et al. (2007, 2010); Merlin et al. (2016), and with the derivation of a 1-GSF equation of motion based on a reconstructed metric Pound et al. (2014). The Teukolsky/FD approach offers great computational efficiency, and ought to be the method of choice for high-precision calculations. Its main weaknesses are that its efficiency degrades quickly with increasing orbital eccentricity (see van de Meent (2016) for a detailed discussion of this point) and that its application to orbits that are not strictly periodic is more subtle. In particular, the incorporation of GSF backreaction effects is much less straightforward than in a TD treatment.

The purpose of this paper is to propose a new computational approach to the GSF, which combines the advantages of the Lorenz-gauge/TD and Teukolsky/FD methods, while avoiding some of their deficiencies. Essentially, our method is a TD application of the metric reconstruction approach. Instead of reconstructing the metric perturbation from a sum over frequency-harmonic modes, as is traditional, we obtain it directly from a certain time-dependent Hertz potential, which, in turn, is obtained by numerically solving the Teukolsky equation in the TD. Our approach thus retains the utility and flexibility of a TD treatment, but offers a much more computationally efficient platform compared to the Lorenz-gauge/TD approach: one solves a single scalarlike equation instead of a coupled set of ten, and the problem of unstable modes is altogether avoided.

TD evolution of the Teukolsky equation has long been used in studies of black hole perturbations, notably by G. Khanna and collaborators Lopez-Aleman et al. (2003). Applications include the study of late-time behavior of vacuum perturbations outside a Kerr black hole Burko and Khanna (2014), investigations into the Cauchy horizon singularity inside black holes Burko et al. (2016), and the modelling of gravitational radiation from particle orbits Taracchini et al. (2014). In these studies, the Teukolsky equation is numerically solved for $\Psi_{0}$ (or $\Psi_{4}$) via time evolution in 2+1 dimensions (2+1D), and the relevant physics (e.g., asymptotic behavior at null or timelike infinity) is directly read off that curvature scalar. In our approach, a time evolution of the Teukolsky equation is used only as a first step in a procedure for reconstructing the local metric perturbation at the particle. Such a direct TD reconstruction has not been attempted so far, to the best of our knowledge.

Furthermore, in our method we do not solve the Teukolsky equation for the physical Weyl scalar $\Psi_{0}$ (or $\Psi_{4}$), as in existing codes. Rather, we solve it for a certain Hertz potential $\Phi$ that is not the Weyl scalar corresponding to the physical perturbation. The field $\Phi$ satisfies the TD Teukolsky equation in vacuum (with boundary conditions similar to those of the Weyl scalar), but the source term in that equation differs from that of $\Psi_{0}$ (or $\Psi_{4}$). A key ingredient of the formulation work to be presented in this paper is a derivation of the point-particle source function for the Hertz potential.111In related literature, the Hertz potential is usually denoted by $\Psi$, with a lowercase $\psi$ denoting the Weyl scalars: $\psi_{0}$ and $\psi_{4}$. Here, to avoid potential confusion, we use $\Phi$ for the Hertz potential; lowercase $\phi$ and $\psi_{0}/\psi_{4}$ will be reserved for the 1+1D projections of $\Phi$ and $\Psi_{0}/\Psi_{4}$, respectively.

Our method is designed with a 1+1D implementation in mind, that is a numerical evolution (either Cauchy or characteristic) on a 2D grid with one temporal dimension and one spatial dimension. In the case of a Schwarzschild background, the Teukolsky equation naturally separates into decoupled 1+1D evolution equations for each multipole (spin-weighted spherical harmonic) mode of $\Phi$. Not so in the Kerr case, where different multipoles ($\ell$ modes) remain coupled. Nonetheless, a 1+1D implementation is still a viable route even in the Kerr case, as demonstrated in Ref. Dolan (2013), at least in situations where mode coupling is relatively weak. In a 1+1D treatment, the particle’s worldline splits the numerical grid into two disjoint domains. The formulation of a source term for $\Phi$ then translates into the prescription of junction conditions that the field $\Phi$ must satisfy on the interface between the two domains—specifically, the jumps in the value of $\Phi$ and a sufficient number of its derivatives across the particle’s worldline. The bulk of our formulation work will be devoted to deriving these junction conditions.

To recap, our goal here is to formulate a TD evolution problem for a certain Hertz potential $\Phi$, from which the physical metric perturbation and the GSF can be derived directly by taking derivatives. We will formulate our method for arbitrary particle orbits in Kerr spacetime, and then, as a proof of principle, we shall present a full numerical implementation for the case of circular geodesic orbits in Schwarzschild spacetime. We leave the full numerical implementation of our method in Kerr to future work.

The structure of the paper is as follows. In Sec. II we review the derivation of the GSF from a reconstructed metric, as proposed by Pound et al. in Pound et al. (2014). We also review the reconstruction procedure starting from a Hertz potential $\Phi$. This will provide the necessary formal background for our work. Section III describes our new method and formulates a practical TD evolution scheme for $\Phi$. This is done for generic orbits in Kerr spacetime. Section IV specializes to circular geodesic orbits in Schwarzschild spacetime, giving all necessary implementation formulas in explicit form for that case. In Sec. V we then present an illustrative numerical implementation of the method, for circular geodesic orbits in Schwarzschild. We conclude in Sec. VI by laying out of a program for the full numerical implementation of our method in Kerr.

Throughout this work we use standard geometrized units, in which $G=1=c$. We adopt the metric signature $({-}{+}{+}{+})$. Greek letters are used for spacetime indices, and a comma denotes a partial derivative, as in $A_{,\alpha}:=\partial A/\partial x^{\alpha}$. Boldface roman indices, as in $e^{\alpha}_{\bf a}$, run over $1,\ldots,4$ and identify tetrad legs. Complex conjugation is denoted by an overbar, as in $\bar{m}^{\alpha}$. Parenthetical indices are symmetrized, as in $A_{(\alpha\beta)}=(A_{\alpha\beta}+A_{\beta\alpha})/2$. We will consider a Kerr background with metric $g^{K}_{\alpha\beta}$, mass parameter $M$ and angular-momentum parameter $a$. We will usually adopt standard Boyer-Lindquist coordinates $(t,r,\theta,\varphi)$. Our sign conventions for the Weyl scalars, the spin coefficients of the Newman-Penrose formalism, and the Hertz potential, are consistent with those of Ref. Merlin et al. (2016) (as summarized in Appendix A therein).

II Self-force from a reconstructed metric: a review

This section reviews the relevant background, while setting up notation and conventions. We first describe the (TD version of the) standard procedure for reconstruction of vacuum metric perturbations in Kerr geometry, then review the case where the perturbation is sourced by an orbiting point particle, and finally summarize the method by Pound et al.. Pound et al. (2014) for computing the GSF experienced by the particle, using the reconstructed metric as input.

II.1 Metric reconstruction in vacuum

A reconstruction procedure for vacuum perturbations was developed long ago by Chrzanowski (Chrzanowski, 1975) and Cohen and Kegeles (Kegeles and Cohen, 1979), with later contributions from Wald (Wald, 1978), Stewart (Stewart, 1979), Lousto and Whiting (Lousto and Whiting, 2002) and others. As is common, we shall refer to it here as the CCK procedure. Considering a vacuum perturbation $h_{\alpha\beta}$ of a Kerr black hole geometry, with corresponding Weyl scalars $\Psi_{0}$ and $\Psi_{4}$, the CCK method prescribes the reconstruction of $h_{\alpha\beta}$ from either $\Psi_{0}$ or $\Psi_{4}$. More precisely, the reconstruction procedure returns a perturbation $h^{\rm rec}_{\alpha\beta}$ that is equal to $h_{\alpha\beta}$ up to (i) some gauge perturbation $h^{\rm gauge}_{\alpha\beta}$, and (ii) a “completion” piece $h^{\rm comp}_{\alpha\beta}$ representing a four-parameter family of simple, stationary and axisymmetric vacuum solutions: mass and angular momentum perturbations of Kerr, and perturbations away from Kerr into Kerr-NUT or C-metric geometries Wald (1973). Thus, the original perturbation is given by

[TABLE]

where $h^{\rm comp}_{\alpha\beta}$ and $h^{\rm gauge}_{\alpha\beta}$ are not fixed within the CCK procedure, and must be determined separately.

There are two variants of the reconstruction procedure, returning $h^{\rm rec}_{\alpha\beta}$ in two different gauges, known as the “ingoing” and “outgoing” traceless radiation gauges (IRG and ORG, respectively). The corresponding gauge conditions are

[TABLE]

along with the trace-free condition

[TABLE]

where $g_{K}^{\alpha\beta}$ is the inverse of the background Kerr metric $g^{K}_{\alpha\beta}$. Here we have introduced Kinnersley’s null tetrad $e^{\alpha}_{\bf a}$ (${\bf a}=1,2,3,4$), whose legs are given, in Boyer-Lindquist coordinates, by

[TABLE]

where

[TABLE]

The legs $e^{\alpha}_{\bf n}$ are all null and mutually orthogonal, except $l^{\alpha}n_{\alpha}=-1$ and $m^{\alpha}\bar{m}_{\alpha}=1$.

The CCK procedure starts with the derivation of a suitable Hertz potential $\Phi$ (one such potential for ORG and another for IRG—call them $\Phi^{\rm ORG}$ and $\Phi^{\rm IRG}$, respectively). This Hertz potential is required to satisfy two differential equations. The first is the master (source-free) Teukolsky equation,

[TABLE]

with $s=+2$ for $\Phi^{\rm ORG}$ and $s=-2$ for $\Phi^{\rm IRG}$. The second equation links $\Phi$ to one of the given Weyl-scalar perturbations $\Psi_{0}$ or $\Psi_{4}$:

[TABLE]

where

[TABLE]

and $\bm{D}_{l}^{4}:=\bm{D}_{l}\bm{D}_{l}\bm{D}_{l}\bm{D}_{l}$, etc. Alternatively, instead of (7) or (8), one can require

[TABLE]

where

[TABLE]

In can be shown Ori (2003) that the combination of (II.1) and (7) [or (II.1) and (10)] determines $\Phi^{\rm IRG}$ uniquely, and similarly the combination of (II.1) and (8) [or (II.1) and (11)] determines $\Phi^{\rm ORG}$ uniquely.

Given $\Phi$, the vacuum perturbation $h^{\rm rec}_{\alpha\beta}$ is obtained via

[TABLE]

where the symmetrization is over the tensorial indices $\alpha\beta$, and ${\cal D}^{\bf ab}$ are certain second-order differential operators. Explicitly, the nonvanishing operators ${\cal D}^{\bf ab}$ are given by

[TABLE]

for the IRG, and by

[TABLE]

for the ORG. Here we have introduced

[TABLE]

and

[TABLE]

II.2 Metric reconstruction for a particle source

The above reconstruction procedure is guaranteed to return a valid vacuum solution $h_{\alpha\beta}^{\rm rec}$ whenever the Weyl scalar one starts with ($\Psi_{4}$ or $\Psi_{0}$) is a solution to the source-free Teukolsky equation. But when matter sources are present, the procedure can fail to return a valid solution even at vacuum points away from any sources. This is true even in the simplest example of a static particle in flat space, in either the IRG or the ORG. As shown first in Barack and Ori (2001), the reconstructed perturbation develops a stringlike singularity, which emanates from the particle in the radial null direction, either outward or inward, or in both directions (depending on the specific choice of gauge within each of the IRG or ORG classes). Reference Pound et al. (2014) introduced a categorization of reconstructed perturbations based on the form of singularity—the two “half-string” classes and the “full string” class—and showed that any reconstructed metric belongs to one of the classes (assuming continuity away from the string). A “no-string” gauge may be formed by joining together the two “regular sides” of two opposite half-string perturbations along a closed surface $\cal S$ through the particle. Such a construction was first introduced by Friedman and collaborators in Keidl et al. (2007, 2010). The no-string perturbation is free from stringlike singularities, but has a gauge discontinuity (and also delta-function distributions Pound et al. (2014)) on the interface $\cal S$.

Specializing to a particle in a bound orbit around a Kerr black hole, we let the particle’s worldline be represented (in Boyer-Lindquist coordinates) by $x_{\rm p}=\{t,r_{\rm p}(t),\theta_{\rm p}(t),\varphi_{\rm p}(t)\}$, where the radius $r_{\rm p}(t)$ is compactly supported outside the Black hole’s event horizon. It is convenient to choose the interface $\cal S$ to be the surface $r=r_{\rm p}(t)$, and we denote the interior of $\cal S$ [i.e., $r<r_{\rm p}(t)$] by ${\cal S}^{-}$ and its exterior [i.e., $r>r_{\rm p}(t)$] by ${\cal S}^{+}$. Then a half-string reconstructed metric has a string singularity in either ${\cal S}^{-}$ or ${\cal S}^{+}$ (and is smooth elsewhere off the particle), and a full-string metric has a string singularity in both. A no-string reconstructed metric is smooth anywhere in both vacuum domains ${\cal S}^{-}$ and ${\cal S}^{+}$ but has a discontinuity across the interface $\cal S$ between them.

Our time-domain metric reconstruction prescription in this work will yield the no-string vacuum perturbations $h_{\alpha\beta}^{\rm rec\pm}$, where henceforth in this paper $\pm$ denotes values in the corresponding vacuum domains ${\cal S}^{\pm}$. We denote the corresponding total, completed perturbations in ${\cal S}^{\pm}$ by

[TABLE]

where $h^{\rm comp\pm}_{\alpha\beta}$ and $h^{\rm gauge\pm}_{\alpha\beta}$ are the contributions from the “completion” and “gauge” perturbations in the corresponding domains [recall Eq. (1)]. The completion piece $h^{\rm comp\pm}_{\alpha\beta}$ has been derived recently in Refs. Merlin et al. (2016); van de Meent (2017a) (for arbitrary bound geodesic motion in Kerr spacetime), and is given there in explicit analytic form. A gauge adjustment $h^{\rm gauge\pm}_{\alpha\beta}$ may need to be included in certain applications (see Refs. Shah and Pound (2015) or van de Meent (2017b) for examples), but in this work we will set $h^{\rm gauge\pm}_{\alpha\beta}\equiv 0$ for simplicity. Thus, our attention in the rest of this paper will be focused on the construction of the perturbation $h^{\rm rec\pm}_{\alpha\beta}$.

II.3 Self-force from a reconstructed metric

In Ref. Pound et al. (2014), Pound et al. obtained a formulation of the GSF starting from a reconstructed metric (in either the IRG or the ORG), complete with practical mode-sum formulas. Two different schemes were described, one based on (either of the two) half-string gauges, and another based on the no-string gauge. The prescription for calculating the GSF in the no-string scheme is as follows. First, given the reconstructed (and completed) perturbations $h_{\alpha\beta}^{\pm}$, introduce the one-sided “force” fields

[TABLE]

defined on ${\cal S}^{\pm}$, respectively, where $\mu$ is the mass of the particle, $\nabla_{\alpha}$ is a covariant derivative compatible with the background (Kerr) geometry, and $u^{\alpha}$ is any smooth extension of the particle’s four-velocity to form a vector field in spacetime. Next, expand each coordinate component of the fields $F_{\alpha}^{\pm}$ in spherical harmonics $Y_{\ell m}(\theta,\varphi)$ on spheres of constant $r,t$, and for each multipole contribution (summed over azimuthal number $m$) take the one-sided radial limits to the particle:

[TABLE]

where $d\Omega:=\sin\theta d\theta d\varphi$. These one-sided $\ell$-mode force contributions are each finite (bounded), and generally grow in amplitude as $\sim\ell$ for large $\ell$. According to Pound et al. (2014), given $F_{\alpha}^{\ell\pm}$, the physical GSF is calculated via

[TABLE]

where $B_{\alpha}$ is the standard, Lorenz-gauge regularization parameter (which depends on the details of the orbit, but not on $\ell$). An analytical expression for $B_{\alpha}$, for generic orbits in Kerr, was first derived in Barack and Ori (2003) and it is given explicitly in Barack (2009).

The mode-sum formula (21) can be implemented in either the IRG or the ORG. While $F^{\ell\pm}\sim\ell$ at large $\ell$, the average $\frac{1}{2}(F^{\ell+}_{\alpha}+F^{\ell-}_{\alpha})$ approaches a constant ($\ell$-independent) value, and the entire summand in (21) is guaranteed to fall off at least as $\sim\ell^{-2}$. Hence, the mode sum converges at least as $\sim 1/\ell$.

In conclusion, we see that knowledge of the “no-string” perturbations $h^{\rm rec\pm}_{\alpha\beta}$ (and their derivatives at the particle) provides sufficient input for calculating the physical GSF. Since $h^{\rm rec\pm}_{\alpha\beta}$ are each a vacuum solution in its corresponding vacuum domain (${\cal S}^{\pm}$), they can be reconstructed from suitable Hertz potentials $\Phi^{\pm}$ by applying the vacuum reconstruction formula (13) in each of the domains:

[TABLE]

In the next section we will lay out our proposed method for deriving the appropriate potentials $\Phi^{\pm}$, via a direct time-domain evolution of the Teukolsky equation. Once $\Phi^{\pm}$ are at hand, the GSF is obtained (in a completed no-string gauge) by substituting for $\Phi^{\pm}$ in Eq. (22), then consequently for $h^{\rm rec\pm}_{\alpha\beta}$ in Eq. (18), for $h^{\pm}_{\alpha\beta}$ in Eq. (19), for $F_{\alpha}^{\pm}$ in Eq. (20), and finally for $F_{\alpha}^{\ell\pm}$ in the mode-sum formula (21).

III 1+1D evolution scheme for the Hertz potentials $\Phi^{\pm}$

Our goal here is to formulate a practical evolution scheme for $\Phi^{\pm}$, suitable for numerical integration in 1+1D. Such a formulation requires three components. First, starting from the master Teukolsky equation (II.1), we need to introduce a suitable decomposition of $\Phi^{\pm}$ into multipole modes and obtain a time-evolution equation for the “time-radial” piece of each of the modes. Second, we need to derive the physical boundary conditions for the fields $\Phi^{\pm}$ and for their multipole modes. Third, we need to formulate junction conditions for the multipole modes on the surface $\cal S$. These three components of the problem are dealt with, each in turn, in the next three subsections.

We aim here to give a full formulation for both $\Phi^{\rm IRG}$ and $\Phi^{\rm ORG}$. Since the former is a solution of the master equation (II.1) with $s=-2$, and the latter is a solution of that equation with $s=+2$, we shall use the notation

[TABLE]

which will allow us to unify the treatment of both gauges. Correspondingly, the vacuum potentials in ${\cal S}^{\pm}$ will be denoted by $\Phi_{s}^{\pm}$, with $s=-2$ for IRG and $s=+2$ for ORG.

III.1 The Teukolsky equation in a 1+1D form

As far as we know, in general, the master Teukolsky equation (II.1) cannot be separated in the time domain, i.e. without first decomposing $\Phi_{s}$ into frequency modes $\sim e^{i\omega t}$; only the azimuthal dependence can be separated, using $\Phi_{s}=\sum_{m}\Phi_{sm}e^{im\varphi}$. A full separation of the angular dependence becomes possible in the special case of $a=0$ (Schwarzschild background), or for $t$-independent perturbations. In both special cases, separation is achieved using a basis of spin-weighted spherical harmonic functions ${}_{s}\!Y_{\ell m}(\theta,\varphi)$ Goldberg et al. (1967). For spins $s=\pm 2$, relevant here, these functions are derived from the standard spherical harmonics $Y_{\ell m}(\theta,\varphi)$ via

[TABLE]

and they satisfy the differential equation

[TABLE]

Since the spherical harmonics satisfy $\bar{Y}_{\ell m}\equiv(-1)^{m}Y_{\ell,-m}$ (where the sign factor is conventional), it is easy to see that we have the symmetry relation

[TABLE]

Our strategy here will be to expand $\Phi_{s}^{\pm}$ in ${}_{s}\!Y_{\ell m}(\theta,\varphi)$ even in the Kerr case. The resulting field equations for the time-radial modes will exhibit coupling between modes of different $\ell$ (though modes of different $m$ will remain decoupled). We will then formulate our evolution problem for that mode-coupled set.

We thus expand the fields $\Phi_{s}^{\pm}$ in the form

[TABLE]

Here, the radial factor $(r\Delta^{s})^{-1}$ is introduced to regulate the behavior of $\phi^{\pm}_{s\ell m}$ at infinity and on the event horizon; its particular form will be explained in Sec. III.2 below. The azimuthal coordinate $\tilde{\varphi}$ is a horizon-regularized version of the Boyer-Lindquist $\varphi$, defined through

[TABLE]

where

[TABLE]

it satisfies $\partial\tilde{\varphi}/\partial r=a/\Delta$. Had we instead used the standard $\varphi$ in Eq. (27), the field $\phi^{-}_{s\ell m}$ would exhibit irregular oscillations $\sim e^{\pm im\Omega_{\rm H}r_{*}}$ near the event horizon Teukolsky (1972). Here $\Omega_{\rm H}=a/(2Mr_{+})$ is the horizon’s angular velocity, and $r_{*}$ is the standard “tortoise” coordinate, satisfying

[TABLE]

Note that $\tilde{\varphi}=\varphi+O(a/r)$, and hence $\tilde{\varphi}$ approaches the standard $\varphi$ at large $r$.

Inserting the expansion (27) into the master equation (II.1), and using Eq. (III.1), we obtain

[TABLE]

where $\tilde{D}$ is a certain partial differential operator independent of $\theta,\varphi$. Note how the two terms $\propto\sin^{2}\theta$ and $\propto\cos\theta$ prevent a full separation of variables when $a\neq 0$ and the field is $t$ dependent. Following the treatment (and notation) of Ref. Barack (2000), we proceed by re-expanding the angular functions ${}_{s}\!Y_{\ell m}\cos\theta$ and ${}_{s}\!Y_{\ell m}\sin^{2}\theta$ in spin-weighted harmonics:

[TABLE]

[TABLE]

Here we have dropped the index $m$ for clarity, and the coefficients are

[TABLE]

and

[TABLE]

Substituting back into (III.1) and using the orthogonality property of the functions ${}_{s}\!Y_{lm}$, we finally separate out the angular dependence, arriving at a mode-coupled set of equations for the time-radial part. For each $s,\ell,m$ (with $s=\pm 2$, $\ell\geq 2$ and $|m|\leq\ell$), it has the form

[TABLE]

where, to avoid confusion, we have omitted the label $\pm$ associated with the domains ${\cal S}^{\pm}$, and used superscript for the multipole label $\ell$. In this equation, $\hat{D}$ is yet another time-radial differential operator (independent of $\theta,\varphi$), and $\cal I$ is a functional describing coupling between each of the $\ell$-modes and their nearest and next-to-nearest neighbours:

[TABLE]

The coupling disappears when $a=0$ or the perturbation is $t$ independent.

To write Eq. (36) explicitly in a convenient form, we introduce the advanced and retarded time coordinates,

[TABLE]

respectively (which reduce to the standard double-null Eddington-Finkelstein coordinates in the Schwarzschild case). Then the 1+1D modal Teukolsky equation (36) takes the explicit form

[TABLE]

where $\partial_{u}$, $\partial_{v}$ and $\partial_{t}$ are taken with fixed $v$, $u$ and $r$, respectively, and we have dropped the indices $s,m$ for improved readability. The radial functions in this equation read

[TABLE]

[TABLE]

[TABLE]

[TABLE]

It is important to reiterate that, even in the Kerr case, the 1+1D Teukolsky equation (III.1) exhibits only a finite coupling between spherical-harmonic modes: each $\ell$ mode couples only to its nearest and next-to-nearest neighbors. This is a remarkable property that can bring much simplification in practice (see below). We also note that the coupling terms contained in $\cal I$ in Eq. (III.1) are expected to be subdominant (compared to, e.g., the term $\propto W$ in that equation) in the problem of interest to us here: these terms are relatively suppressed by factors of orders $(a\omega)^{2}$ or $a\omega$, where $\omega$ is a characteristic frequency of the perturbation, which, in most relevant scenarios (and for relevant values of $m$) is considerably smaller than $1/M$, giving suppression factors considerably smaller than unity.222The situation is less clear for strongly bound orbits around a near-extremal black hole, where one may expect $\omega\sim 1/M$. This case will require further investigation. We suspect that coupling between modes remains subdominant even in that case, thanks to the time-delay effect noted in Gralla et al. (2016) (recalling that coupling terms all involve $t$ derivatives). This can be used to one’s advantage in numerical implementations, as we further discuss in Sec. III.3 below.

The term $\propto\phi^{\ell}_{,tt}$ may be eliminated from Eq. (III.1) by introducing a modified ($\ell,m,s$-dependent) radial coordinate $\tilde{r}_{*}$, satisfying

[TABLE]

The expression under the square root here is positive definite, and smaller than $1$, since $0<4Ka^{2}<1$ and also (it can be shown) $0<C_{0}^{\ell}<1$ for all relevant $\ell,m,s$. We note that $\beta\to 1$ in both limits $r\to\infty$ and $r\to r_{+}$, meaning the modified coordinate $\tilde{r}_{*}$ coincides with the usual $r_{*}$ in both limits. Introducing also the modified advanced and retarded times $\tilde{v}:=t+\tilde{r}_{*}$ and $\tilde{u}:=t-\tilde{r}_{*}$, Eq. (III.1) becomes

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

The form (III.1) (without a $\phi^{\ell}_{,tt}$ term) is more conveniently amenable to a finite-difference representation, especially in double-null-type coordinates.

We envisage a numerical implementation of Eq. (III.1) using a finite-difference scheme based on $\tilde{u},\tilde{v}$ coordinates. For stationary modes of the perturbation, and in the Schwarzschild case, the entire second line of (III.1) vanishes, and the equation takes a simple form, ready for numerical implementation $\ell$ by $\ell$. In the general Kerr case, the equation can first be recast in a matrix form, introducing the vector variable $\vec{\phi}:=\{\phi^{\ell=1},\phi^{\ell=2},\ldots,\phi^{\ell_{\rm max}}\}$, where $\ell_{\rm max}$ is a suitable cutoff; in many problems, including the GSF problem, the large-$\ell$ truncation error may be controlled and made sufficiently small. The matrix equation can then be discretized and solved as in the Schwarzschild case, this time obtaining all $\ell$ modes at once. Since the coupling between modes is finite and “weak” (modes couple only to nearest and next-to-nearest neighbors), the matrices involved are band diagonal and hence comfortably amenable to numerical manipulation.

III.2 Boundary conditions for $\phi^{\pm}_{s\ell m}$

Typically, and in the GSF problem in particular, we require the physical, “retarded” solution for the metric perturbation, i.e., the one for which there is no radiation coming in from past null infinity, and no radiation coming out of the past event horizon. In frequency-domain implementations, and also in Cauchy-type time evolutions, this requirement is imposed via a choice of boundary conditions on suitable timelike boundaries. In a characteristic-type 1+1D evolution of the kind we have in mind here, the numerical domain has no timelike boundaries, and the solution is completely determined once initial data are specified on two initial characteristic rays. In principle, one should be able to select the retarded solution via a suitable choice of characteristic initial data (though in practice we shall adopt a much simpler approach, described in Sec. V below). But even though boundary conditions are not actively imposed at each time step, knowledge of the form of physical boundary conditions is still important, for a number of reasons. First, as we shall see, such knowledge informs our choice of regulator functions [specifically, the factor $(r\Delta^{s})^{-1}$ in Eq. (27)] that control the behavior of the numerical field at large retarded and advanced times. Our choice of regulator will be such that the physical solution is globally bounded in magnitude. Second, once a numerical solution is obtained, it is important to check that it is indeed the desired physical solution. Third, there are cases where pieces of the physical field may be determined analytically (see Sec. IV.1 for an example), and in such cases the choice of a particular solution requires knowledge of the physical boundary conditions.

Our goal now, therefore, is to prescribe physical boundary conditions for the time-radial fields $\phi^{\pm}_{s\ell m}(t,r)$. In Sec. III.2.1 we will consider the behavior at null infinity, and in Sec. III.2.2 the behavior at the event horizon. Section III.2.3 will discuss the special case of stationary modes. In all cases we will base our analysis on a study of the asymptotic form of solutions to the 1+1D time-domain Teukolsky equation (III.1). We recall the coordinates $(\tilde{u},\tilde{v})$ are interchangeable with the standard $(u,v)$ in both asymptotic limits.

III.2.1 Behavior at null infinity

Consider solutions of (III.1) that for $r\gg M$ have the asymptotic forms $\sim r^{\alpha}e^{-i\omega u}$ (“outgoing waves”) or $\sim r^{\beta}e^{-i\omega v}$ (“incoming waves”), for some frequency $\omega\neq 0$ (the case $\omega=0$ will be considered separately below). We determine the powers $\alpha$ and $\beta$ by substituting each of these asymptotic-form Ansätze in Eq. (III.1), then expanding in powers of $1/r$ (at fixed $t$), noting $\tilde{V}(r)=s/r+O(1/r^{2})$ and $\tilde{K}(r),\tilde{U}(r),\tilde{W}(r)=O(1/r^{2})$. From the leading-order term of each expansion we readily find $\alpha=0$ and $\beta=2s$. Thus we have the two asymptotic solutions

[TABLE]

It is easy to check, using (27) and the reconstruction formula (22), that the first solution yields a reconstructed metric with a large-$r$ asymptotic form $h_{\alpha\beta}^{{\rm rec}+}\sim e^{-i\omega u}/r$ (in suitable Cartesian coordinates), representing a purely outgoing wave—hence the designation “physical”. The second solution yields a perturbation $h_{\alpha\beta}^{{\rm rec}+}\sim e^{-i\omega v}$ (multiplied by some factor of $r$), which does not represent a purely outgoing wave—hence “nonphysical”.

We note that, for the physical solution, the magnitude of $\phi_{s\ell m}^{+}$ approaches a constant, generally nonzero value at future null infinity ($v\to\infty$ for any fixed $u$). This behavior, which is computationally desirable, was achieved by introducing the factor $(r\Delta^{s})^{-1}\sim r^{-2s-1}$ in Eq. (27). Without this factor, the physical solution would blow up as $\sim r^{3}$ for $s=-2$ (IRG), or would fall off rapidly, as $\sim r^{-5}$, for $s=+2$ (ORG), both types of behavior being computationally problematic.

Note also that, for $s=+2$, any nonphysical solution blows up rapidly (as $\sim r^{4}$) at infinity. This means that, in an ORG implementation, a numerical solution $\phi_{s\ell m}^{+}$ that can be demonstrated to remain bounded at infinity (even as the finite-difference step size tends to zero) is automatically guaranteed to be the physical solution, i.e., the one satisfying the correct, outgoing boundary conditions. An IRG calculation does not share this convenient feature: for $s=-2$, nonphysical modes decay rapidly (as $\sim r^{-4}$) at infinity and would be hard to identify in the data. In this sense, it is computationally advantageous to calculate the external field $\phi_{s\ell m}^{+}$ in the ORG.

III.2.2 Behavior at the event horizon

Moving on to consider the behavior along the horizon, examine the form of (III.1) for small $\Delta$. Noting $\tilde{K}(r),\tilde{V}(r),\tilde{W}(r)=O(\Delta)$ while $\tilde{U}(r)=O(\Delta^{0})$, we find that, at leading order in $\Delta$, the equation reduces to

[TABLE]

This equation admits the two pure-mode asymptotic solutions

[TABLE]

It can be checked that the first solution yields a reconstructed perturbation that (in suitable, horizon-regular coordinates such as $\{v,r,\theta,\tilde{\varphi}\}$) has the asymptotic form $h_{\alpha\beta}^{\rm rec-}\sim e^{-i\omega v}$, representing a purely ingoing wave at the future event horizon. This is therefore the “physical” solution.333From Eq. (27) we see that, for the physical solution, the Hertz potential $\Phi_{s}^{-}$ itself is $\propto\Delta^{-s}$ at the horizon. This apparent irregular behavior owes itself simply to the irregularity of the Kinnersley tetrad at the horizon; see, e.g., Sec. V of Barack and Ori (1999). For the physical solution, $\phi^{-}_{s\ell m\omega}$ approaches a finite, generally nonzero value at the horizon ($u\to\infty$ for any fixed $v$), which is computationally desirable. Indeed, to achieve this was the purpose of our introduction of a regulator factor $\Delta^{-s}$ in Eq. (27).

It can also be checked that the second solution in (III.2.2) produces a reconstructed perturbation that, in coordinates regular on the past event horizon, has the asymptotic behavior $h_{\alpha\beta}^{\rm rec-}\sim e^{-i\omega u}$ there.444Note that our coordinate $\tilde{\varphi}$ is not regular on the past horizon, and it is this coordinate irregularity that gives rise to the singular factor $e^{-2im\Omega_{\rm H}r_{*}}$ in the “nonphysical” solution in (III.2.2). When we say that this solution is nonphysical we do not refer to this coordinate irregularity but to the fact that the solution represents outgoing waves at the past horizon. This solution thus represents nonphysical outgoing waves at the past horizon. Note that for $s=-2$ the nonphysical modes blow up as $\sim\Delta^{-2}$ at the horizon. Hence, in an IRG reconstruction it is sufficient to demonstrate the boundedness of the internal numerical solution $\phi^{-}_{s\ell m\omega}$ at the horizon in order to establish that it represents the true, physical solution (up to numerical error). In the ORG, on the other hand, nonphysical modes decay rapidly (as $\sim\Delta^{2}$) near the horizon and would be hard to identify there if they existed in the numerical data. Thus, there is a computational advantage in calculating the internal field $\phi_{s\ell m}^{-}$ in the IRG.

III.2.3 Stationary modes

Finally, let us consider the asymptotic behavior of stationary, $t$-independent (or, equivalently, $\omega=0$) modes. Substituting the Ansatz $\phi\sim r^{\alpha}$ in Eq. (III.1) and considering the leading-order beavior at $r\gg M$, one obtains $\alpha=s-\ell$ or $\ell+s+1$, and hence the two asymptotic solutions

[TABLE]

The designations ‘physical’ and ‘nonphysical’ here come from examining the behavior of the physical perturbation associated with each solution. To eliminate gauge dependence, it is instructive to consider, for example, the corresponding Weyl scalars $\Psi_{0}$ or $\Psi_{4}$. Using Eq. (27) with (7) and (8) [or with (10) and 11)] shows $\Psi_{0},\Psi_{4}\sim r^{-\ell-3}$ for the first solution in (III.2.3), and $\Psi_{0},\Psi_{4}\sim r^{\ell-2}$ for the second, in both the IRG and ORG cases. Thus, the first solution corresponds to a perturbation whose curvature decays at infinity (hence ‘physical’), while the second solution (‘nonphysical’) corresponds to a perturbation whose curvature does not decay at infinity (and, for $\ell>2$, it actually blows up there).

We move on to consider the behavior on the event horizon. Two stationary asymptotic solutions there are

[TABLE]

where, importantly, the first solution has a regular Taylor expansion on the horizon, while the second solution also contains a high-order contribution of the form $(\log\Delta)\times$ a Taylor series in $\Delta$. It can be checked that the corresponding Weyl scalars have the behavior $\Psi_{s}\sim\Delta^{-s}$ for the first solution, and $\Psi_{s}\sim\Delta^{0}$ (plus higher-order log terms) for the second solution. As explained (e.g.) in Sec. V of Ref. Barack and Ori (1999), carefully taking into account the irregularity of the Kinnersley tetrad at the horizon, for a smooth physical perturbation it is not the Weyl scalars themselves that are regular (smooth) at the horizon, but rather the product $\Delta^{s}\Psi_{s}$. Applying this criterion to our solutions, we have that $\Delta^{s}\Psi_{s}$ is perfectly smooth for the first solution (“physical”), but non-smooth for the second solution (“nonphysical”). In the latter case, $\Delta^{s}\Psi_{s}$ blows up like $\Delta^{-2}$ for $s=-2$, while for $s=+2$ the differentiability is spoiled by the $\log\Delta$ term.

We note that, in terms of our variable $\phi^{\pm}(t,r)$, the stationary piece of the physical perturbation is bounded everywhere, just like the rest of the perturbation. For nonphysical stationary perturbations, $\phi^{-}$ blows up on the horizon in the IRG case, and $\phi^{+}$ blows up at infinity in both the IRG and ORG cases. This is true for all relevant values of $\ell$, i.e. $\ell\geq 2$.

III.2.4 Summary

In summary, our choice of time-radial fields $\phi^{\pm}_{s\ell m}(t,r)$ is such that, in terms of these variables, the physical, retarded solution is bounded both at infinity and on the horizon—and, in fact, anywhere else in the computation domain. This is a convenient feature, computationally. Furthermore, we have noted that all nonphysical ORG solutions $\phi^{+}_{s\ell m}$ blow up at infinity, while all nonphysical IRG solutions $\phi^{-}_{s\ell m}$ blow up at the horizon. Thus, boundedness of an ORG solution at infinity implies that the correct outgoing boundary conditions are satisfied at infinity, and boundedness of an IRG solution at the horizon implies that the correct ingoing boundary conditions are satisfied on the horizon.

Unfortunately, we cannot make a stronger statement: In the ORG case, boundedness of the solution on the horizon does not necessarily mean that boundary conditions are satisfied there, because nonphysical ORG modes are subdominant on the horizon. Similarly, in the IRG case, boundedness of the solution at infinity does not necessarily imply, in general, that boundary conditions are satisfied there, because nonphysical IRG modes are subdominant at infinity (with the exception of stationary modes, which blow up there). This situation suggests that, at least from a computational point of view, it would be convenient to work with a mixed-gauge field composed of the ORG $\phi^{+}_{s\ell m}$ and the IRG $\phi^{-}_{s\ell m}$. However, the formulation of jump conditions on the orbit (see below) would then be harder, and one would also need to generalize the self-force formulation to accommodate the possibility of such a mixed-gauge perturbation. It would be worth exploring the mixed-gauge idea in future work, but here we shall stick with the more straightforward single-gauge approach.

III.3 Jump conditions for $\phi^{\pm}$ across $\cal S$

Finally, we need a set of rules that relate the fields $\phi^{+}_{s\ell m}$ and $\phi^{-}_{s\ell m}$ along the particle’s timelike trajectory in the 1+1D domain. Specifically, we need the “jumps”

[TABLE]

as well as the jumps in the first derivatives—say, $\left[\phi_{,\tilde{u}}\right]$ and $\left[\phi_{,\tilde{v}}\right]$. (In this subsection we occasionally, where possible, omit the indices $s\ell m$ for brevity.) The jumps in higher-order derivatives may also be required, depending on the particular numerical method implemented and on the order of numerical convergence sought. Our goal now is to describe a method for determining these jumps, as functions along the particle’s worldline. We will assume that the jumps $\left[\psi_{\pm 2}\right]$ in the physical Weyl scalars, and in their derivatives, are already known. These jumps may be deduced directly, in analytic form, from the source term of the (1+1D version of the) Teukolsky equation, without needing to solve the equation for the Weyl scalars themselves—a specific example will be worked out explicitly in Appendix A. Thus, we will be seeking to determine the jumps in the Hertz potential and its derivatives in terms of the jumps in the physical Weyl scalars and their derivatives.

Our starting point is Eqs. (7), (8), (10) and (11), as applied to the vacuum solutions $\Phi_{s}^{\pm}$. Recall these are relations that the vacuum Hertz potential must satisfy given the physical Weyl scalars $\Psi_{0}$ or $\Psi_{4}$. First, we need to obtain the 1+1D version of these equations. To this end, we substitute the expansion (27) for $\Phi^{\pm}_{s}$ on the left-hand side of each of these relations, and on the right-hand side we substitute for $\Psi_{4}=:\varrho^{4}\Psi_{s=-2}$ and $\Psi_{0}=:\Psi_{s=+2}$ using the analogous expansions

[TABLE]

for $s=\pm 2$.

Considering first the relations (7) and (8) (whose frequency-domain versions are often referred to as “radial inversion formulas”), we observe that the differential operators on the left-hand side do not couple between different $\ell$ modes: the equations relate each (spin-weighted spherical-harmonic) $\ell$ mode of the Hertz potentials $\bar{\Phi}^{\pm}_{s}$ and their derivatives to the same $\ell$ mode of the Weyl scalars $\Psi^{\pm}_{-s}$ (with $m$ opposite in sign). Using (26) and the orthogonality of ${}_{s}\!Y_{\ell m}$ we obtain, for each $\ell,m$,

[TABLE]

[TABLE]

where the differential operators are

[TABLE]

here $\partial_{v}$ is taken with fixed $u$, and $\partial_{u}$ is taken with fixed $v$. We observe that (56) and (57) are effectively ordinary differential equations for $\bar{\phi}^{{\rm IRG}\pm}_{\ell m}$ and $\bar{\phi}^{{\rm ORG}\pm}_{\ell m}$, with independent variables $v$ and $u$, respectively.

In general, the alternative relations (10) and (11) (“angular inversion formulas”) are less useful here, because they feature operators that do mix between different $\ell$ modes, giving rise to mode-coupled relations; even worse, the equations relate each $\ell$ mode of $\Phi^{\pm}$ to an infinite number of $\ell$ modes of $\Psi_{s}^{\pm}$. The coupling disappears only in the Schwarzschild case, $a=0$, where, in fact, the angular inversion formulas become extremely simple (since the operators $\tilde{\cal L}_{s}$ and ${\cal L}_{s}$ reduce to, respectively, spin-lowering and spin-raising operators when they act on spin-weighted spherical harmonics). This simplicity was noted previously by Lousto and Whiting in Lousto and Whiting (2002), where they considered 1+1D metric reconstruction in Schwarzschild. Here, however, our ambition is to treat the more general Kerr case, so we will utilize the “radial” inversion equations—even as (in the next two sections) we consider a Schwarzschild example.

Focusing thus on the 1+1D inversion formulas (56) and (57), we first note that the general solution for each of these two ODEs can be written down in a simple closed form (involving four nested integrals with respect to $v$ or $u$, respectively). Based on the form of these general solutions it is straightforward to show (in analogy with Ori’s analysis in Ori (2003)) that $\phi^{{\rm IRG}\pm}_{\ell m}$ and $\phi^{{\rm ORG}\pm}_{\ell m}$ each admits a unique particular solution that also satisfies the vacuum Teukolsky equation (III.1), as required. This confirms the uniqueness of the Hertz potential in the reconstruction procedure.

However, the above closed-form particular solutions for $\phi^{{\rm IRG}\pm}_{\ell m}$ and $\phi^{{\rm ORG}\pm}_{\ell m}$ involve integrals of the Weyl scalars ($\psi_{2}^{\pm}$ and $\psi_{-2}^{\pm}$, respectively) along rays extending to infinity and down to the event horizon, and to evaluate them in practice would require solving the appropriate sourced Teukolsky equations as a preliminary step. In our method we wish to bypass this preliminary step, and work directly with the Hertz potential; we wish to have no recourse to knowledge of the Weyl scalars themselves (except the values of their jumps across the particle, which, as mentioned, are accessible directly from the source of the Teukolsky equation). Our goal, therefore, is to express the jumps in $\phi$ (and its derivatives) in terms of the jumps in $\psi$ (and its derivatives) alone. In what follows we describe a procedure that achieves that.

Let us start with the IRG case (the ORG case will follow analogously). Evaluating the difference between the ‘$+$’ and ‘$-$’ versions of Eq. (56) at $r=r_{\rm p}(t)$ yields a relation of the form

[TABLE]

in which $f_{n}(r_{\rm p})$ are some smooth functions along the orbit, and the jumps on the right-hand side are assumed known. For a reason that will become clear shortly, we also need the $\tilde{v}$ derivative of (56), which yields

[TABLE]

where $\tilde{f}_{n}(r_{\rm p})$ are some other smooth functions, and the jumps on the right-hand side are also known. Our goal is to determine the jumps $\left[\phi\right]$ and $\left[\phi_{,\tilde{v}}\right]$, as well as $\left[\phi_{,\tilde{u}}\right]$, as functions along the orbit (here and in the following discussion we omit the labels $\ell,m$ and IRG for brevity).

Using an overdot to denote $d/d\tau$, where $\tau$ is proper time along the orbit, we write

[TABLE]

where $\tilde{u}_{\rm p}(\tau)$ and $\tilde{v}_{\rm p}(\tau)$ are the values of $\tilde{u}$ and $\tilde{v}$ at a worldline point with proper time $\tau$. Equation (61) gives $\left[\phi_{,\tilde{u}}\right]$ in terms of $\dot{\left[\phi\right]}$ and $\left[\phi_{,\tilde{v}}\right]$. In Eq. (62) we replace $\left[\phi_{,\tilde{v}\tilde{u}}\right]$ in favor of $\left[\phi\right]$, $\left[\phi_{,\tilde{v}}\right]$ and $\left[\phi_{,\tilde{u}}\right]$ using the Teukolsky equation (III.1), and hence express $\left[\phi_{,\tilde{v}\tilde{v}}\right]$ in terms of $\left[\phi\right]$ and $\left[\phi_{,\tilde{v}}\right]$ (and their $\tau$ derivatives) alone. (For $a\neq 0$ this relation will involve also the jump in the coupling terms $\cal I$. Let us ignore these terms for a moment to simplify the discussion; we shall return to them momentarily.) Next, we write

[TABLE]

and use (III.1) again to replace all mixed-derivative jumps with lower-order-derivative jumps. We thus express $\left[\phi_{,\tilde{u}\tilde{u}}\right]$ and (in turn) $\left[\phi_{,\tilde{v}\tilde{v}\tilde{v}}\right]$ in terms of $\left[\phi\right]$ and $\left[\phi_{,\tilde{v}}\right]$ (and their first and second $\tau$ derivatives) alone.

Proceeding in a similar way, we can determine $\left[\phi_{,\tilde{u}\tilde{u}\tilde{u}}\right]$ and $\left[\phi_{,\tilde{v}\tilde{v}\tilde{v}\tilde{v}}\right]$ in terms of $\left[\phi\right]$, $\left[\phi_{,\tilde{v}}\right]$ and their first, second and third $\tau$ derivatives; and finally we can determine $\left[\phi_{,\tilde{u}\tilde{u}\tilde{u}\tilde{u}}\right]$ and $\left[\phi_{,\tilde{v}\tilde{v}\tilde{v}\tilde{v}\tilde{v}}\right]$ in terms of $\left[\phi\right]$, $\left[\phi_{,\tilde{v}}\right]$ and their first, second, third and fourth $\tau$ derivatives. Equations (59) and (60) can thus be written as a coupled set of ODEs for $\left[\phi\right]$ and $\left[\phi_{,\tilde{v}}\right]$:

[TABLE]

[TABLE]

where $a_{n}(\tau),\ldots,d_{n}(\tau)$ are certain smooth functions along the worldline (depending only on $r_{\rm p},\dot{r}_{\rm p},\ldots,d^{4}r_{\rm p}/d\tau^{4}$, as well as on $\ell$). The terms collected under ‘${\cal I}$ terms’ are certain linear combinations of the coupling term ${\cal I}(\phi^{\ell\pm 1},\phi^{\ell\pm 2})$ and its $v$ and $u$ derivatives (up to third derivatives), which have entered the relations via our use of the Teukolsky equation (III.1). Unfortunately, the general explicit form of Eqs. (III.3) and (III.3) is too unwieldy to be presented here in any meaningful way, but it can be straightforwardly obtained using computer algebra tools, following the procedure described above. In Sec. IV we will present explicit expressions for the special case of circular geodesic orbits in Schwarzschild spacetime [where all ${\cal I}$ terms drop, and the coefficients $a_{n},\ldots,d_{n}$ depend only on $r_{\rm p}$(=const)].

Let us for the moment ignore the coupling terms in Eqs. (III.3) and (III.3). Then these equations constitute a coupled set of fourth-order ODEs for the jumps $\left[\phi_{\ell m}\right]$ and $\left[\phi_{\ell m,\tilde{v}}\right]$ as functions along the orbit. How to solve these equations depends on the particular problem under consideration. If the perturbation is sourced by a particle on a fixed bound geodesic orbit, then the jumps may be assumed to exhibit the same periodicity as the orbit (two fundamental frequencies, in general), and the ODEs (III.3) and (III.3) can then be conveniently recast as a set of algebraic equations, one for each frequency (the assumption of periodicity then effectively selects a particular solution of the ODEs).555By solving the jump equations (III.3) and (III.3) frequency by frequency we would not be reverting to the standard frequency-domain approach to metric reconstruction: the field equation (III.1) would still be solved in the time domain. If the perturbation is sourced by a slowly evolving orbit (e.g., under the effect of the self-force), then one should still be able to obtain a frequency-by-frequency algebraic solution at some initial point along the orbit, then solve the set (III.3) and (III.3) as ODEs, starting from these initial conditions. For nonperiodic (parabolic- or hyperbolic-type) orbits, initial conditions may be formulated at infinity, using an asymptotic analysis. In Sec. IV we give the explicit physical solution of Eqs. (III.3) and (III.3) for the special case of circular geodesic orbits in Schwarzschild spacetime.

How should the $\ell$-mode coupling terms in Eqs. (III.3) and (III.3) be dealt with, in the Kerr case? As already mentioned, we expect the $\cal I$ term in the field equation (III.1) to be small, in general, in a certain relative sense. The coupling terms in the jump equations (III.3) and (III.3) will be small in the same sense. One could then incorporate these terms perturbatively, using an iterative scheme: In the first iteration, Eqs. (III.3) and (III.3) are solved for each $\ell,m$ with the coupling terms dropped. The solutions are then used to calculate the $\cal I$ terms in Eqs. (III.3) and (III.3), and the equations are solved again, with these $\cal I$ terms as sources. One keeps iterating in this manner until a sufficiently convergent solution is achieved; the smaller the coupling terms are in relative magnitude, the less iterations should be required. How computationally tasking this procedure may prove to be would depend on the number of iterations required and on whether the jump equations are solved as ODEs or via a mode decomposition. Yet we expect the computational cost of calculating the jumps to be negligible compared with the cost of solving the field equation in the time domain.

At any rate, once the jumps $\left[\phi\right]$ and $\left[\phi_{,\tilde{v}}\right]$ have been computed as functions along the orbit, the jumps in any higher $\tilde{v}$ and $\tilde{u}$ derivatives of $\phi$ are calculable algebraically using the order-reduction procedure described a couple of paragraphs above: $\left[\phi_{,\tilde{u}}\right]$ is obtained from (61), then $\left[\phi_{,\tilde{v}\tilde{v}}\right]$ is obtained from (62), and so on.

Finally, the ORG version of the problem is dealt with in a completely analogous manner, this time starting with Eq. (57). The procedure again yields ODEs of the form (III.3)-(III.3), but now with $\left[\phi_{,\tilde{v}}\right]$ replaced with $\left[\phi_{,\tilde{u}}\right]$ and, on the right-hand side, $\bar{\psi}_{2,\ell,-m}$ replaced with $\bar{\psi}_{-2,\ell,-m}$. The explicit form of the coefficients $a_{n}(\tau),\ldots,d_{n}(\tau)$ and of the $\cal I$ terms will also differ.

III.4 Summary of proposed evolution scheme

In summary, our evolution scheme for the Hertz potential consists of the 1+1D evolution equation (III.1), the asymptotic boundary conditions (III.2.1) and (III.2.2) (“physical”), and jump conditions across the particle’s orbits, given as solutions to Eqs. (III.3) and (III.3). We have sketched how the actual jump conditions for the modes of the Hertz potential may be obtained in practice, depending on the orbital configuration.

In the next two sections we will illustrate the application of our method with the concrete example of a particle moving in a fixed circular geodesic orbit around a Schwarzschild black hole. Section IV will formulate the evolution problem as applied to this case, including explicit expressions for the jumps across the particle. Section V will present a numerical implementation. For concreteness and brevity, we will concentrate on the IRG problem.

IV circular orbits in Schwarzschild spacetime

Specialized to $a=0$ (Schwarzschild case) and $s=-2$ (IRG), the 1+1D Teukolsky equation (III.1) becomes

[TABLE]

with

[TABLE]

[TABLE]

where we have introduced

[TABLE]

Note that, in the Schwarzschild case, (i) the 1+1D Teukolsky equation does not couple between $\ell$ modes, so individual modes evolve independently of each other; (ii) the equation has no reference to the azimuthal number $m$ (as expected, by virtue of the background’s spherical symmetry); and (iii) the modified coordinates $\tilde{u},\tilde{v}$ reduce to the standard coordinates $u,v$.

Next consider the jump equations (III.3)-(III.3). For $a=0$, the coupling terms drop. Specializing further to a circular geodesic orbit with radius $r_{\rm p}=r_{0}(={\rm const}>2M)$, the coefficients on the left-hand side of (III.3)-(III.3) work out to be

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

We have used here the fact that, for circular geodesics, $\dot{u}_{\rm p}=\dot{v}_{\rm p}(=\gamma_{0})$ and all higher-order $\tau$ derivatives of $u_{\rm p}$ and $v_{\rm p}$ vanish.

Furthermore, for a circular geodesic orbit we may assume that $[\phi_{\ell m}]$ depends on time solely via $e^{-im\Omega t}$, where

[TABLE]

is the orbital angular velocity. Hence, in Eqs.  (III.3)-(III.3) $d/d\tau$ may be replaced with $-im\Omega(dt_{\rm p}/d\tau)=-im\Omega\gamma_{0}$. These equations then become algebraic, taking the simple form

[TABLE]

where $a_{\Sigma}=\sum_{n=0}^{3}(-im\Omega\gamma_{0})^{n}a_{n}$ and similarly for $b_{\Sigma}$, $c_{\Sigma}$ and $d_{\Sigma}$. The solutions are

[TABLE]

where $\tilde{\Delta}:=a_{\Sigma}d_{\Sigma}-b_{\Sigma}c_{\Sigma}$, and we have restored all labels. We find

[TABLE]

which, we note, is positive definite.

Equations (78) and (79) give the desired jumps in the Hertz potential and its $v$ derivative in terms of the jumps in the Weyl scalar ($\Psi_{0}$) corresponding to the physical perturbation. The latter jumps are easily obtained from the distributional source of the Teukolsky equation satisfied by the Weyl scalar. In Appendix A we show how the jumps $\left[\partial_{v}\bar{\psi}_{2,\ell,-m}\right]$ and $\left[\bar{\psi}_{2,\ell,-m}\right]$ are derived, and give them explicitly for circular geodesic orbits in Schwarzschild. Once $\left[\phi\right]$ and $\left[\phi_{,v}\right]$ are known, the jumps in other derivatives of $\phi$ can be obtained iteratively, by means of the procedure described above. For instance,

[TABLE]

and so on.

IV.1 Analytical solutions for $m=0$

Axisymmetric modes (those with $m=0$) are also stationary, and admit simple analytic solutions. We write these solutions here explicitly, as they will be useful in testing our numerical implementation in the next section.

For $m=0$ we have $\phi^{\ell}_{,t}=0$, and the homogeneous Teukolsky equation (67) reduces to

[TABLE]

where a prime denotes $d/dr$. Two linearly independent solutions are

[TABLE]

where ${\sf P}_{\ell}^{m}$ and ${\sf Q}_{\ell}^{m}$ are associated Legendre functions of the first and second kinds, respectively, $x:=(r-M)/M$, and the solutions have been normalized so as to render the Wronskian $\ell$ independent:

[TABLE]

The solution $\phi^{\ell}_{P}$ blows up as $\sim r^{\ell-1}$ at infinity but has a regular Taylor expansion at the event horizon. The solution $\phi^{\ell}_{Q}$ is regular at infinity (where it falls off as $\sim r^{-\ell-2}$) but blows up as $\sim f^{-2}$ at the horizon.

Recalling Eqs. (III.2.3) and (III.2.3), we see that $\phi^{\ell}_{P}$ satisfies physical boundary conditions on the horizon (but not at infinity), while $\phi^{\ell}_{Q}$ satisfies physical boundary conditions at infinity (but not on the horizon). Therefore, a unique physical solution is given by

[TABLE]

The coefficients $C^{\pm}_{\ell}(r_{0})$ are determined from the two jump conditions $\phi^{+}_{\ell 0}(r_{0})-\phi^{-}_{\ell 0}(r_{0})=\left[\phi_{\ell 0}\right]$ and $(\phi^{+}_{\ell 0})^{\prime}(r_{0})-(\phi^{-}_{\ell 0})^{\prime}(r_{0})=\left[\phi^{\prime}_{\ell 0}\right]$, giving

[TABLE]

where we have substituted for the Wronskian from Eq. (84). The jumps $\left[\phi_{\ell 0}\right]$ and $\left[\phi^{\prime}_{\ell 0}\right]=(2/f_{0})\left[\phi_{\ell 0,v}\right]$ are calculated using Eq. (78) with Eqs. (104) and (105). For $m=0$ we obtain the simple expressions

[TABLE]

where ${\cal Y}^{\ell}:={}_{2}\!\bar{Y}_{\ell 0}\left(\frac{\pi}{2},0\right)$, and ${\cal Y}^{\ell}_{\theta}$ and ${\cal Y}^{\ell}_{\theta\theta}$ are, respectively, the first and second derivatives of ${}_{2}\!\bar{Y}_{\ell m}(\theta,0)$ with respect to $\theta$, evaluated at $\theta=\pi/2$. We note that ${\cal Y}^{\ell}$ and ${\cal Y}^{\ell}_{\theta\theta}$ vanish for all odd values of $\ell$, while ${\cal Y}^{\ell}_{\theta}$ vanishes for all even values of $\ell$. Inspecting Eqs. (87) and (88) we consequently find that the jumps $\left[\phi_{\ell 0}\right]$ and $\left[\phi^{\prime}_{\ell 0}\right]$ are real for even $\ell$ and imaginary for odd $\ell$.

In summary, the axially symmetric piece of the IRG $\ell$-mode Hertz potential in and out of a circular geodesic orbit in Schwarzschild spacetime ($\phi^{-}_{{-2}\ell 0}$ and $\phi^{+}_{{-2}\ell 0}$, respectively) is given analytically by Eq. (IV.1), with the coefficients $C^{\pm}_{\ell}(r_{0})$ given in Eq. (IV.1) [with Eqs. (87) and (88)]. The $m=0$ solution is purely real for even $\ell$ and purely imaginary for odd $\ell$. We have checked that our analytic solution (IV.1) agrees with that obtained in Ref. Merlin et al. (2016) using a different method, namely starting with $\Psi_{4}$ and using the “angular” inversion formula (10).

V Numerical Implementation

In this section we present an illustrative numerical implementation of our method, specializing to circular geodesic orbits in Schwarzschild spacetime. Our goal here is twofold: First, we will be able to provide some test for our formulation by comparing numerical results with results obtained analytically or using other methods. Second, we aim to illustrate the kind of numerical implementation strategy we have in mind, which, we believe, can be used to tackle more general cases. Here we do not seek sophistication in our numerical technique, and we do not concern ourselves with questions of computational performance or code optimization. We leave such matters to future work.

V.1 Method

Our code solves for the IRG Hertz-potential modes $\phi^{\rm IRG}_{\ell m}(r,t)$ by evolving a finite-difference version of the hyperbolic equation (67) on a fixed 1+1D grid in double-null ($u,v$) coordinates, subject to the jump conditions (78) and (79). This double-null approach follows the general strategy set out in Refs. Barack and Lousto (2005); Barack and Sago (2007, 2010). Our finite-difference method, detailed in Appendix B, is a straightforward second-order convergent scheme. By this we mean that its global (time-accumulated) error scales as $h^{2}$, where $h\times h$ are the null-coordinate dimensions of a single grid cell. (To achieve this requires that the local finite-difference error at each grid cell scales with a higher power of $h$; see Appendix B for details.) Higher-order convergence may be desirable in future applications. It can be achieved in a straightforward manner, following, e.g., the methods of Lousto (2005) or Barack and Sago (2010).

The code, implemented in Mathematica, takes as input the orbital radius $r_{0}$ and the mode numbers $\ell,m$, and returns the field $\phi^{\rm IRG}_{\ell m}(t,r)$. The evolution starts from initial data on two characteristic rays $v=r_{*}(r_{0})=:v_{0}$ and $u=-r_{*}(r_{0})=:u_{0}$, which intersect on the particle’s orbit at $t=0$. The numerical integration then proceeds along successive $v={\rm const}>v_{0}$ rays, with the jump conditions (78) and (79) imposed whenever the particle’s worldline [represented by the line $v=u+2r_{*}(r_{0})$] is crossed. The future boundaries of the grid are taken at some large values of $v$ and $u$ (approximating null infinity and the event horizon, respectively), so the numerical domain has no timelike causal boundaries where boundary conditions would have been required. For initial data we take, for simplicity, $\phi^{\rm IRG}_{\ell m}(v=v_{0})=0=\phi^{\rm IRG}_{\ell m}(u=u_{0})$. This choice, which violates the jump conditions at $r=r_{0}$, produces an initial burst of nonphysical (“junk”) radiation, which, however, dies out over time (with some $\ell$-dependent inverse-power law). The evolution proceeds until the level of junk radiation (as determined from the residual nonstationarity of the numerical solution) drops below a set threshold, and one then discards the early, junk-dominated part of the data.

How can one be sure that the late-time solution thus obtained is the physical, “retarded” solution sought for? Our artificial choice of characteristic initial data means that the early part of the evolution likely contains a nonphysical component that violates the retarded boundary condition (it would contain, in particular, radiation that comes in from past null infinity and out of the past horizon). Since the numerical solution satisfies the correct jump conditions on the particle (up to numerical error), this nonphysical component may be thought of as a homogeneous (vacuum) perturbation superposed on the true, inhomogeneous physical solution. With our choice of initial conditions, we also know that this vacuum perturbation has an initial compact support. It is a well known feature of the vacuum Teukolsky equation (see, e.g., Barack (2000)) that solutions of an initial compact support die off at late time and have no manifestation at timelike infinity. Thus, we expect the nonphysical component of the numerical solution to die off at late time, and the full solution to relax to its physical value.666In order to obtain a late-time solution other than the retarded one, one would need to adjust the form of characteristic initial data at late retarded and advanced times, so as to represent radiation coming out of past null infinity and/or out of the past horizon. Our argument is that whenever the initial data are compactly supported, as in our implementation, the solution will relax to the retarded one at late time.

As discussed in Sec. III.2, one can also check a posteriori whether the numerical solution is the physical one. For stationary ($m=0$) modes, in our IRG case, all nonphysical solutions blow up either at infinity or at the horizon (or at both), so demonstrating boundedness of our solution (even in the limit of diminishing step size) should suffice to establish that the solution is indeed the physical one. For nonstationary ($m\neq 0$) modes, boundedness on the horizon implies that the internal solution $\phi^{-}$ is the physical one, though, in our IRG implementation, we have no such direct test for the external solution $\phi^{+}$.

Our numerical experiments indeed suggest that the Hertz-potential field always settles down to the true, physical solution at late time (to within our controlled numerical accuracy). This is shown most convincingly by comparing with analytic results (for $m=0$) and frequency-domain numerical calculations (for $m\neq 0$), in which the correct boundary conditions were explicitly imposed. We present some of this evidence below.

V.2 Sample results

Figure 1 shows sample results for $(\ell,m)=(2,0)$, demonstrating the relaxation of the numerical solution towards the analytically known solution (IV.1) at late time. The figure displays the behavior of the fields $\phi_{20}^{+}$ and $\phi_{20}^{-}$ on the particle (fixed on a circular geodesic orbit at $r_{0}=7M$), as a function of time. The inset illustrates how the asymptotic agreement with the analytical solution improves with decreasing grid size. Our highest-resolution run in this case (with $h=M/8$) took about 3 minutes on a modest laptop, outputting the correct result with a mere $\sim 5\times 10^{-5}$ relative error. If we use as a rough measure of numerical error the difference between the highest-resolution value and the value obtained via a Richardson extrapolation to $h\to 0$ (assuming quadratic convergence), we see that the extrapolated value is consistent with the analytical result to within the estimated error.

Figures 2 and 3 further test the $m=0$ solution for a variety of $\ell$ values, again comparing with the analytical solution (IV.1). Higher multipoles exhibit smaller-scale features, the resolution of which demands smaller step sizes and hence extra computational resources. On the other hand, higher multipoles also relax faster at late time, allowing a shorter evolution, which somewhat alleviates the computational burden.

Figures 4 and 5 display numerical results for a nonaxisymmetric mode of the Hertz potential ($m=1$). In this case we do not have analytical results, but we can compare with numerical solutions obtained using the frequency-domain approach of Merlin and Shah Merlin and Shah (2015), who kindly provided us with numerical data generated by their code. As demonstrated in Fig. 4, our results agree with theirs to within our (small) estimated numerical error. Figure 5 illustrates the behavior of the $m=1$ solution on a $t$=const slice, showing waves away from the particle and a discontinuity on it.

Once we have at hand the Hertz potential, reconstruction of the IRG metric perturbation becomes straightforward, through a mode-by-mode application of Eq. (22). Figure 6 shows, as an illustration, some of the components of the $(\ell,m)=(2,0)$ mode of reconstructed perturbation in the vicinity of the particle. The analytical solution is also shown, for comparison. We observe $h_{\alpha\beta}^{\rm rec+}\neq h_{\alpha\beta}^{\rm rec-}$ on the particle, as expected.

VI Summary and future applications

A recent body of work Keidl et al. (2007, 2010); Pound et al. (2014); Merlin and Shah (2015); Shah and Pound (2015); van de Meent and Shah (2015); van de Meent (2016); Merlin et al. (2016); van de Meent (2017a), reviewed in Sec. II, showed how the gravitational self-force experienced by a point particle in orbit around a Kerr black hole can be calculated from a reconstructed metric perturbation, starting from frequency-domain solutions of the fully separated Teukolsky equations. In this approach, a suitable version of the metric perturbation is reconstructed as a sum over frequency-harmonic modes (subject to a certain regularization procedure). In the current paper we have developed a time-domain version of that approach, which circumvents the need to calculate individual frequency modes of the perturbation. Instead, each multipole ($\ell,m$ mode) of the perturbation is reconstructed (numerically) directly as a function of time. This should enable tackling some self-force problems that are not directly amenable to a frequency-domain treatment, including the very important problem of self-consistently evolving the orbit under backreaction from the self-force.

In our new method, developed in Sec. III, the metric perturbation is reconstructed multipole by multipole from scalarlike potentials $\phi_{\ell m}^{+}(t,r)$ and $\phi_{\ell m}^{-}(t,r)$ defined in the vacuum regions $r>r_{\rm p}(t)$ and $r<r_{\rm p}(t)$, respectively. Each of these fields is a solution of a 1+1D version of the vacuum Teukolsky equation [Eq. (III.1), or its alternative form (III.1)], which, for $a\neq 0$, exhibits coupling between different $\ell$ modes. The crucial and most difficult step in our method is the determination of jump conditions that relate between $\phi_{\ell m}^{+}$ and $\phi_{\ell m}^{-}$ across the particle’s trajectory $r=r_{\rm p}(t)$ in the 1+1D domain. A general method for doing so, for arbitrary motion in Kerr, was developed in Sec. III.3. Section IV derived the jump conditions explicitly for the example of circular geodesic orbits in Schwarzschild, and Sec. V presented a full numerical implementation for that case.

Looking ahead, let us first consider the extension of our treatment to an arbitrary motion (for now still in Schwarzschild spacetime). The most nontrivial aspect of the extension concerns the explicit determination of jump conditions for $\phi_{\ell m}^{\pm}$. For this, one needs to refer back to the general jump equations (III.3) and (III.3), which are a coupled set of fourth-order ODEs along the orbit. Strategies for solving this set were discussed in Sec. III.3. For strictly periodic orbits (including any bound geodesic orbit, even in the Kerr case), one can decompose the jumps by means of a discrete Fourier series, and solve the resulting ODEs frequency mode by frequency mode. But we see the main potential of our method in its ability to tackle nonperiodic setups, including the problem of parabolic/hyperbolic flybys or captures, and the problem of self-consistent evolution. In such cases, the jump equations may have to be solved as ODEs, with suitable initial conditions. For orbits that start or end at infinity, such conditions may most conveniently be imposed there. In the case of self-consistent evolution, initial values for the jump conditions could be well approximated by assuming exact orbital periodicity at an initial moment where the evolution is strongly adiabatic. In both cases, it remains to explore the numerical stability of the set of ODEs and develop a robust and efficient method for solving them along the orbit.

Moving on to implementation in Kerr, the obvious additional complication comes from $\ell$-mode coupling. The 1+1D Teukolsky equation (III.1) [or its alternative version (III.1)] couples between each $\ell$ mode and its nearest and next-to-nearest neighbors, so the numerical evolution problem would need to be recast in a matrix form, with “all” $\ell$ modes solved for simultaneously. In practice, a large-$\ell$ cutoff would need to be introduced, and its error controlled. As already mentioned, we expect the narrow band-diagonal form of the matrix operator to be an advantage, computationally. We have preformed preliminary experiments evolving our 1+1D Teukolsky equation in vacuum, for Kerr, with very encouraging results.

An additional complication, in the Kerr case, is the occurrence of coupling terms also in the jump equations (III.3) and (III.3). Here, too, the equation would need to be recast in a suitable matrix form and solved for all relevant $\ell$ modes simultaneously. We do not expect this additional hurdle to present a serious difficulty in practice.

The method developed in this paper is, to the best of our knowledge, the only method thus far proposed capable of fully tackling the problem of time-domain calculations of metric perturbations from a point particle in Kerr spacetime. (The Lorenz-gauge approach of Dolan and Barack Dolan and Barack (2013) is hampered by gauge-instability problems that are yet to be resolved.) We are therefore keen to see it further developed, and have embarked on a program to study further applications.

Acknowledgements

We are grateful to Cesar Merlin and Maarten van de Meent for providing us comparison data generated, upon our request, using their respective frequency-domain codes. We thank Charalampos Markakis for useful inputs during early discussions leading to this project. We gratefully acknowledge support from the European Research Council under the European Union’s Seventh Framework Programme FP7/2007-2013/ERC, Grant No. 304978.

Appendix A Calculating the discontinuity in the Weyl scalar and its derivative

In Sec. IV we described the derivation of the (IRG) Hertz potential for circular geodesic orbits around a Schwarzschild black hole. Equations (78) and (79) therein gave necessary jump conditions for the $\ell$-mode Hertz potential across the particle’s orbit in the $1+1$D domain. These jumps were expressed in terms of the jumps in the corresponding $\ell$ mode of the Weyl scalar $\Psi_{0}$ associated with the physical perturbation—more precisely, in terms of $\left[\psi_{s=2,\ell,m}\right]$ and $\left[\partial_{v}\psi_{s=2,\ell,m}\right]$. In this appendix we show how the latter jumps are derived directly from the source of the $s=2$ Teukolsky equation, and provide explicit expressions for them.

Let $\Psi_{0}\equiv\Psi_{s=2}$ be the Weyl scalar associated with the physical metric perturbation sourced by the point particle. We expand it in $s=2$ spherical harmonics, as in Eq. (55):

[TABLE]

Thus $\psi_{2\ell m}\equiv\psi^{+}_{2\ell m}$ for $r>r_{0}$ and $\psi_{2\ell m}\equiv\psi^{-}_{2\ell m}$ for $r<r_{0}$, where $\psi^{\pm}_{2\ell m}$ are the “homogeneous” $\ell$-mode fields featured in Eq. (55). The field $\psi_{2\ell m}$ satisfies the $1+1$D $s=2$ inhomogeneous Teukolsky equation (here specialized to $a=0$ and $s=2$) with a suitable distributional source term $T$ corresponding to the energy-momentum of our geodesic pointlike particle:

[TABLE]

where hereafter indices $s\ell m$ are dropped for brevity, and the potentials are

[TABLE]

[TABLE]

We now need the explicit form of the $\ell$-mode source $T$, for the case of a point particle in a circular geodesic orbit. This can be worked out as explained (e.g.) in Sec. III.A of Ref. Merlin et al. (2016) (where $T$ was derived for $s=-2$ and $m=0$), by using Eq. (A10) therein with the energy-momentum tensor (21) therein, and then decomposing into spin-2 spherical harmonics. We obtain the distributional form

[TABLE]

where the coefficients are given by

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

We recall $y=M/r_{0}$ and $f_{0}=1-2M/r_{0}$, and we have further introduced here

[TABLE]

and

[TABLE]

with ${\cal Y}_{\theta}$ and ${\cal Y}_{\theta\theta}$ being, respectively, the first and second derivatives of ${}_{2}\!\bar{Y}_{\ell m}\left(\theta,\Omega t\right)$ with respect to $\theta$, evaluated at $\theta=\pi/2$. We have taken $\varphi_{\rm p}(t)=\Omega t$, with $\Omega$ being the orbital angular velocity, and used $d{\cal Y}/dt=-im\Omega{\cal Y}$ and $d^{2}{\cal Y}/dt^{2}=-m^{2}\Omega^{2}{\cal Y}=-m^{2}(M/r_{0}^{3}){\cal Y}$.

Next, we substitute the Ansatz

[TABLE]

where $\Theta(\cdot)$ is the standard Heaviside step function, into the field equation (90). We note all resulting terms that are proportional to $\Theta(r-r_{0})$ or to $\Theta(r_{0}-r)$ vanish, by virtue of $\psi^{\pm}(t,r)$ being homogeneous solutions. The remaining terms are supported on $r=r_{0}$ only, and are each proportional to $\delta$, $\delta^{\prime}$ or $\delta^{\prime\prime}$. We use the distributional identities

[TABLE]

[valid for any smooth function $F(r)$] to reexpress the coefficients of $\delta$, $\delta^{\prime}$ and $\delta^{\prime\prime}$ in terms of $r_{0}$ (and $t$) only, and then compare the values of these coefficients across both sides of Eq. (90), recalling the form of $T$ in Eq. (A). From the coefficient of $\delta^{\prime\prime}$ one immediately obtains

[TABLE]

and subsequently comparing the coefficients of $\delta^{\prime}$ and $\delta$ uniquely determines the jumps $\left[\psi\right]=\psi^{+}(r_{0})-\psi^{-}(r_{0})$ and $\left[\psi^{\prime}\right]=(\psi^{+})^{\prime}(r_{0})-(\psi^{-})^{\prime}(r_{0})$. We obtain

[TABLE]

[TABLE]

We recall $\lambda=(\ell+2)(\ell-1)$.

In summary, and restoring all suppressed indices, Eqs. (104) and (105) give the desired jumps $\left[\psi_{s=2,\ell,m}\right]$ and $\left[\psi^{\prime}_{s=2,\ell,m}\right]$ in explicit form. The jump in the $v$ derivative, needed as input for Eqs. (78) and (79), is obtained via

[TABLE]

Appendix B Finite-difference scheme

In this appendix we describe the finite-difference scheme used in Sec. V for numerically computing the Hertz-potential modes in the time domain. Our treatment is based on a characteristic evolution of the vacuum field equation $(\ref{Teukolsky1+1Sch})$ on a fixed 1+1D uniform mesh in double-null $(v,u)$ coordinates. The numerical domain is depicted in Fig. 7. We start with initial conditions on the rays $v=v_{0}$ and $u=u_{0}$, and evolve along successive rows (“rays”) of constant $v$. At each step of the integration, we approximate the value of the field at a grid point using the already-known values at a few grid points in its recent “causal past”. The form of the finite-difference formula applied depends on the position of the grid point with respect to the particle’s worldline (represented by the vertical line $v-v_{0}=u-u_{0}$): For points that are sufficiently far from the worldline (those marked ‘V’ in Fig. 7) we apply a certain “vacuum” formula, whereas for points in the vicinity of the worldline (‘L’, ‘R’ or ‘RR’) or on it (‘P’) we apply modified formulas that involve the known jumps in the value of the field and its derivatives across the worldline. Below we shall describe the finite-difference schemes applied for each type of point.

But first we introduce some notation. We let $h$ denote the constant stepping interval (in each of $v$ and $u$), so that the $u,v$-coordinate dimensions of a single grid cell are $h\times h$. Consider an arbitrary grid point C at the coordinate position $(u,v)=(u_{c},v_{c})$, for some $u_{c}$ and $v_{c}$. With reference to the point C, we denote by $\phi(n,k)$ (where $n,k\in\mathbb{N}$) our finite-difference approximant for the value of the field at the grid point with coordinates $(u,v)=\left(u_{c}-nh,v_{c}-kh\right)$. Our goal is to obtain, for each such point C, the value $\phi(0,0)$ in terms of the values $\phi(n\geq 0,k\geq 0)$ (excluding $n=0=k$) that are computed in previous steps of the integration. For our second-order-convergence scheme it will prove sufficient to use the six values $\phi(1,0)$, $\phi(0,1)$, $\phi(1,1)$, $\phi(2,0)$, $\phi(0,2)$ and $\phi(2,2)$. This requires storing field values at three successive $v$=const rays at a time.

We note that the field value at a P-type point is ambiguous: the field there admits, in general, two different one-sided values $\phi^{-}$ and $\phi^{+}$. We have found it convenient to assign to each P point a concrete one-sided value, and in practice we take it to be the left-hand value $\phi^{-}$. (The choice of side here is arbitrary; assigning the value $\phi^{+}$ would work just the same.) This amounts to taking the orbit to pass “just to the right” of the line of P points. Henceforth, whenever we refer to the value $\phi(n,k)$ at a P-type point, we mean $\phi^{-}(n,k)$. Of course, given $\phi^{-}$ at a P-type point, the right-hand value may be immediately recovered using

[TABLE]

where $[\phi]$ is the known jump in $\phi$ at that point.

B.1 V-type (vacuum) grid points

V-type points are those with coordinates $(u,v)$ satisfying $\Delta v-\Delta u\geq 3h$ or $\Delta v-\Delta u\leq-2h$, where henceforth $\Delta v:=v-v_{0}$ and $\Delta u:=u-u_{0}$. Let us consider an arbitrary V-type grid point at $(u,v)=(u_{c},v_{c})$. To obtain a finite-difference expression for the field value at that point, it is convenient to consider the integral of both sides of the field equation (67) over the grid cell with vertices $(u_{c},v_{c})$, $(u_{c}-h,v_{c})$, $(u_{c},v_{c}-h)$ and $(u_{c}-h,v_{c}-h)$ (i.e., the grid cell with our V-type point as its upper vertex). We represent that integral with the symbol $\int_{\Diamond}$.

Let us consider Eq. (67) term by term. First, we have

[TABLE]

which is exact, and does not involve any finite-difference approximation. Next, for the fourth term in (67) we obtain

[TABLE]

where $r_{c}$ is the value of $r$ at the point $(u_{c},v_{c})$. This finite-difference approximation suffices for our purpose: since the total number of V-type grid points scales as $h^{-2}$, the accumulated error from the omitted $O(h^{4})$ terms should be of $O(h^{2})$ at most, consistent with the sought-for quadratic convergence.

To obtain finite-difference approximations for the second and third terms in (67), we apply the following, more systematic procedure (following Barack and Sago (2010)). We first formally write $\phi$ as a double Taylor expansion in $u$ and $v$ about the point $(u_{c},v_{c})$, keeping terms up to $O(h^{2})$:

[TABLE]

This should be valid in a vacuum neighborhood of $(u_{c},v_{c})$. Taking $(u_{c},v_{c})$ as our “$(0,0)$” point, we then apply Eq. (110) at each of the vacuum points $(1,0)$, $(0,1)$, $(1,1)$, $(2,0)$, $(0,2)$ and $(2,2)$. This yields six algebraic equations for the six coefficients $c_{ij}$ ($0\leq i,j\leq 2$) in terms of the values $\phi(1,0)$, $\phi(0,1)$, $\phi(1,1)$, $\phi(2,0)$, $\phi(0,2)$ and $\phi(2,2)$, assumed known. We solve these equations and substitute the coefficient values back into Eq. (110), to obtain an approximation formula for $\phi(u,v)$ valid through $O(h^{2})$ in the vicinity of $(u_{c},v_{c})$. Using this formula with a Taylor expansion of the (smooth) potentials $U(r)$ and $V(r)$ about $r=r_{c}$, we obtain approximation formulas for the two terms $U(r)\phi_{,u}$ and $V(r)\phi_{,v}$, valid through $O(h)$ near $(u_{c},v_{c})$. Finally, integrating over the grid cell, we obtain

[TABLE]

[TABLE]

in which the radial functions $U$, $U_{,r_{*}}$, $V$ and $V_{,r_{*}}$ are all evaluated at $r=r_{c}$.

Adding together the expressions (108)–(B.1) and equating to zero [the cell integral of the right-hand side of (67)], we obtain the desired finite-difference formula for $\phi(0,0)$ in the vacuum case:

[TABLE]

where the only nonvanishing coefficients are

[TABLE]

Here, all radial functions are evaluated at $r=r_{c}$.

B.2 Near-particle grid points

RR, R, L and P-type points are those with coordinates $(v,u)$ satisfying $\Delta v-\Delta u=+2h$, $+h$, $-h$ and [math], respectively. For such points, the above vacuum scheme (113) does not quite work as it stands, because it involves field values at points on both sides of the particle’s worldline, where the field has a jump discontinuity.

To account for the discontinuity we apply the following procedure (again following Barack and Sago (2010)). For a given near-particle point (RR, R, L or P) with coordinates $(u_{c},v_{c})$, we write down two separate Taylor expansions,

[TABLE]

each suitable for points on the corresponding side of the particle: $\phi^{-}$ for points on the “left” ($r\leq r_{0}$) and $\phi^{+}$ for points on the “right” ($r>r_{0}$). We apply the formula to the six points $(1,0)$, $(0,1)$, $(1,1)$, $(2,0)$, $(0,2)$ and $(2,2)$ as before. This now yields 6 equations for the 12 coefficients $c_{ij}^{\pm}$ ($0\leq i,j\leq 2$), so the system is underdetermined. However, the known jump conditions across the particle’s worldline provide additional constraints: we get 6 additional relations between the $c_{ij}^{\pm}$’s by imposing the known jumps $[\phi]$, $[\phi_{,u}]$, $[\phi_{,v}]$, $[\phi_{,uu}]$, $[\phi_{,uv}]$ and $[\phi_{,vv}]$ at some point along the worldline. We choose that point to be $(0,2)$ for the RR case, $(0,1)$ for the R case, $(1,0)$ for the L case, and $(0,0)$ for the P case. Solving the resulting set of 12 equations, we obtain the 12 coefficients $c^{\pm}_{ij}$ in terms of the field values at the above six points and the above 6 jumps. Substituting back into (115) again gives an approximation formula for $\phi^{\pm}$ in the vicinity of $(u_{c},v_{c})$, which can be used to construct approximate expressions for each of the terms in the field equation (67) near $(u_{c},v_{c})$. Finally, we integrate these expressions over the grid cell, using either ‘$+$’ or ‘$-$’ values, as appropriate (the former for cases RR and R and the latter for case L; for P points we use $\phi^{+}$ on the right-hand half of the cell and $\phi^{-}$ on its left).

Following this procedure, and collecting all terms, we arrive at the near-particle finite-difference formula

[TABLE]

in which $X\in\{RR,R,L,P\}$, $\phi_{\rm Vac}(0,0)$ is the vacuum expression given in Eq. (113), and the form of the jump terms $J_{\rm X}$ depends on the point type in question:

[TABLE]

Here, all radial functions are evaluated at the point $(0,0)$ under consideration, while the jumps are evaluated at the appropriate worldline point, indicated in subscript. Note that we have truncated the expressions at $O(h^{2})$, leaving an $O(h^{3})$ error in Eq. (116). We allow ourselves to do so because the total number of near-particle (RR, R, L and P-type) points scales only like $1/h$, so the accumulated error from the omitted $O(h^{3})$ local term should scale as $h^{2}$ at most, consistent with quadratic convergence. We also note that $J_{\rm P}$ may be written, through the relevant $O(h^{2})$, in the alternative form

[TABLE]

which involves jumps in the field and its first derivatives only [at the expense of requiring the jumps at the point $(1,1)$ as well]. This form may be simpler to use in practice.

Equations (113) and (116) constitute our finite-difference scheme. We see that, at the required order, the implementation of our scheme requires knowledge of the jumps in the field and its first derivatives only. Jumps in higher derivatives may be required for scheme with convergence faster than quadratic.