Hyperboloidal method for frequency-domain self-force calculations

Gravitational self-force theory is the leading approach for modeling gravitational wave emission from small mass-ratio compact binaries. This method perturbatively expands the metric of the binary in powers of the mass ratio. The source for the perturbations depends on the orbital configuration, calculational approach, and the order of the perturbative expansion. These sources fall into three broad classes: (i) distributional, (ii) worldtube, and (iii) unbounded support. The latter, in particular, is important for emerging second-order (in the mass ratio) calculations. Traditional frequency domain approaches employ the variation of parameters method and compute the perturbation on standard time slices with numerical boundary conditions supplied at finite radius from series expansions of the asymptotic behavior. This approach has been very successful, but the boundary conditions calculations are tedious, and the approach is not well suited to unbounded sources where homogeneous solutions must be computed at all radii. This work develops an alternative approach where hyperboloidal slices foliate the spacetime, and compactifying coordinates simplify the boundary treatment. We implement this approach with a multidomain spectral solver with analytic mesh refinement and use the scalar-field self-force on circular orbits around a Schwarzschild black hole as an example problem. The method works efficiently for all three source classes encountered in self-force calculations and has distinct advantages over the traditional approach. For example, our code efficiently computes the perturbation for orbits with extremely large orbital radii ($r_p>{10}^{5}M$) or modes with very high spherical harmonic mode index ($\ell \ge 100$). Our results indicate that hyperboloidal methods can play an essential role in self-force calculations.

Figure 1

Carter-Penrose diagram for the Schwarzschild exterior region. Thin, dashed lines depict standard Schwarzschild time surfaces $t=\text{constant}$ extending between the bifurcation sphere $\mathcal{B}$ at the horizon $r=2M$ and spacelike infinity $i^0$ as $r\to \infty$. The intersection of these time surfaces near $\mathcal{B}$ and $i^0$ imply a coordinate singularity. The domain must be truncated and boundary data must be imposed near $\mathcal{B}$ and $i^0$. Thick, solid lines depict hyperboloidal time surfaces $\tau =\text{constant}$ extending between the black hole horizon ${\mathcal{H}}^{+}$ at $\sigma =1$ and future null infinity ${\mathcal{I}}^{+}$ as $\sigma =0$ given by Eq. (25). These coordinates provide a smooth foliation on the full exterior domain which means that both the horizon and null infinity can be included in the computational domain. No external boundary conditions are needed to study perturbations.

Figure 2

AMR to better resolve functions with strong gradients. Bottom panel: a Chebyshev-Lobatto grid is considered for the coordinate $\chi \in [-1,1]$ according to Eq. (97). Middle panel: for $x_{\mathrm{B}}=-1$, the AMR map (96) populates the grid points around the left boundary. Top panel: for $x_{\mathrm{B}}=1$, the AMR map (96) populates the grid points around the right boundary. Examples with AMR parameter $\kappa =3$.

Figure 3

Left panel: real parts of hyperboloidal retarded field ${\overline{\varphi }}_{\ell m}^{\mathrm{ret}}$ for the angular mode $(\ell ,m)=(1,1)$ with $r_p=6M$ and $N=60$. The numerical domain extends from future null infinity, $\sigma =0$, to the future event horizon, $\sigma =1$. The inset demonstrates exponential decay of Chebyshev coefficients indicating spectral convergence. Right panel: same fields as in the left panel but for angular mode $(\ell ,m)=(100,0)$ and resolution $N=120$. We need higher resolution at high mode numbers because of the steep gradient around the particle. Insets demonstrate slower spectral convergence than in left panel.

Figure 4

Numerical convergence for the retarded field $|{\varphi }_{\ell m}^{\mathrm{ret}}|$ with $(\ell ,m)=(1,1)$ as function of particle’s orbit $r_p$. Despite the exponential decay of error against a reference solution with $N^{\mathrm{ref}}=150$, the decay rate is slower at large orbits $r_p$ due to steep gradients.

Figure 5

The $\ell$ mode contributions to the $t$ component of SSF, $F_{\ell t}$, for a particle on a circular orbit of radius $r_p=10M$. The modes of $F_t^{\ell }$ converge exponentially until machine precision roundoff is encountered near $\ell =17$.

Figure 6

The $\ell$ mode contribution to the $r$ component of the SSF, $F_{\ell r}$, for a particle on a circular orbit of radius $r_p=10M$. For large $\ell$ the unregularized modes (blue squares) grow linearly. After subtracting the leading regularization parameter, the modes (light blue triangles) tend to a constant for large $\ell$. Further subtracting the next regularization parameter, the regular modes (orange circles) falls off as ${\ell }^{-2}$ for large $\ell$. The convergence of the $\ell$-mode sum is then accelerated using higher-order regularization parameters with each additional parameter changing the large $\ell$ behavior by ${\ell }^{-2}$. After all the known regularization parameters are subtracted, the modes quickly reach machine roundoff. Note that the agreement with the expected large-$\ell$ behavior is excellent out to $\ell =100$ (when the contributions are above machine precision).

Figure 7

Top panel: effective source solution using four domains with angular mode $(\ell ,m)=(1,1)$ and particle location $r_p=6M$. Bottom panel: decay of Chebyshev coefficients in all domains demonstrates spectral convergence. The inset displays the real part of residual field ${\overline{\varphi }}_{\ell m}^{\mathcal{R}}$ with a separate color for each domain.

Figure 8

Same setup as in Fig. 7 but for a high angular mode $(\ell ,m)=(100,0)$.

Figure 9

Convergence tests for the residual field $|{\varphi }_{\ell m}^{\mathcal{R}}|$ displaying the error against a reference solution with $N^{\mathrm{ref}}=150$ according to Eq. (95). Top panel: angular modes $(\ell ,m)=(1,1),(50,50)$, and (100,0) with the particle at $r_p=6M$. Higher angular modes require slightly higher resolution. Bottom panel: various particle locations $r_p/M=\{3.01,10,100,1000\}$ with a fixed angular mode $(\ell ,m)=(1,1)$. The exponential decay rate is lower for larger orbits.

Figure 10

The $\ell$-mode contribution to the $r$ component of self-force, $F_{\ell r}$, computed using the effective-source method. The direct output of the hyperboloidal calculation with the effective source are shown as the (orange) circles which falls off as ${\ell }^{-2}$. We then use higher-order regularization parameters to accelerate the convergence of the series. The more rapidly convergent series quickly reaches machine precision roundoff. The results presented here for the regularized force are, as expected, the same as the results from the mode-sum approach—see Fig. 6. With our setup the effective-source method is more efficient than the distributional source and mode-sum approach. This is because with the distributional source large gradients of the field occur near the particle which necessitates $N=150$ Chebyshev nodes in each domain where the effective-source only requires $N=50$.

Figure 11

The $r_p$ derivative of the scalar field, ${\overline{\psi }}_{\ell m}^{\mathrm{ret}}$, computed for $r_p=6M$ and $(l,m)=(1,1)$. The source for ${\overline{\psi }}_{\ell m}^{\mathrm{ret}}$ is unbounded but our approach handles it with ease. The inset shows the exponential convergence (until machine roundoff is reached) for the Chebyshev coefficients in each domain.

Figure 12

Convergence of ${\psi }_{\ell m}^{\mathrm{ret}}$ with increasing number of Chebyshev nodes, $N$, for different orbital radii. In all cases the convergence is exponential but for large radius orbits the convergence can be quite slow. The rate of convergence can be improved with analytic mesh refinement—see Fig. 15.

Figure 13

In self-force calculations for large orbits, strong gradients form in the compact coordinate $\sigma$ around the particle.

Figure 14

Top panel: Chebyshev coefficients of $\mathrm{Re}({\overline{\varphi }}_{\ell m}^{\mathrm{ret}})$ for $r_p=1000M$ and $(\ell ,m)=(1,1)$ in domain 2 between the particle and the horizon. The optimal value for the AMR parameter is $\kappa =3.85$. Bottom panel: optimal AMR parameters plotted against $r_p/M$ with the fit Eq. (100) for different angular modes.

Figure 15

Left panel: exponential error decay for large orbits with optimal AMR. Comparison with Fig. 12 demonstrates the power of AMR in calculating accurate self-force results for large radius orbits. Right panel: convergence for different spherical harmonic modes with $r_p/M={10}^6$ including a high-$\ell$ mode where machine precision is reached around $N=80$. This computation would be prohibitively resource intensive in standard self-force calculations.

Figure 16

Left panel: comparison of the Newtonian-normalized $t$ component of the self-force, $\widehat{F_t}$, with its 4.5PN expansion. The (dark blue) squares show our numerical results for $\widehat{F_t}$. This data approaches the leading (normalized) PN result for large radius orbits. When we subtract the leading PN term from the numerical data we get the (light blue) triangles. For large $r_p$ this data approaches a $V^2$ reference curve as expected from the PN series in Eq. (101). When we subtract the first subdominant term in the PN series and see that the residual (orange circles) falls off as $V^3$, as expected. We repeat this procedure with the remaining terms in the PN series to compute the other data and find agreement with the PN series (until machine roundoff is reached). This shows our numerical results are accurate even for extreme large radius orbits with $r_p\sim {10}^6$. Right panel: comparison of the $r_p$ derivative of the Newtonian-normalized $t$ component of the self-force obtained from flux-balance laws with the 4.5PN expression. This figure constructed in the same way as the left panel except we subtract terms from the PN series in Eq. (102). Again we see that our numerical results are accurate even for extreme large radius orbits with $r_p\sim {10}^6$. Accurate results for large radius orbits for self-force problem with unbounded sources are very difficult to achieve with the standard variation of parameters approach. Our hyperboloidal method can thus be instrumental in future precision comparisons with PN theory for self-force problems with unbounded sources, e.g., second-order self-force calculations.

Figure 17

AMR effect over ${\overline{\varphi }}_{\ell m}^{\mathrm{ret}}$ on domain 1 (between future null infinity and particle’s orbit). Top panel: Chebyshev coefficients of $\mathrm{Re}({\overline{\varphi }}_{\ell m}^{\mathrm{ret}})$ with $r_p=1000M$ and $(\ell ,m)=(1,1)$ in domain 1 between the particle and future null infinity. The improvement with AnMR is not as compelling as in domain 2 presented in Fig. 14. Bottom panel: optimal AMR parameters plotted against $r_p/M$ with the fit Eq. (c1).

Figure 18

AMR for scenarios with worldtube sources, requiring a code with 4 domains. Chebyshev coefficients of $\mathrm{Re}({\overline{\varphi }}_{\ell m}^{\mathcal{R}})$ for $r_p=1000M$ and $(\ell ,m)=(1,1)$ in domain 3 (between the particle ${\sigma }_p$ and the worldtube boundary ${\sigma }_{+}$). Here, the optimal value for the AMR parameter is $\kappa =3.5$ and the overall tendency on $r_p$ and $(\ell ,m)$ follows the tendency observed in Eq. (100).

Figure 19

Exponential error decay for large orbits with optimal AMR for worldtube sources. It demonstrates the challenging computation of $r_p/M={10}^6$ and $(\ell ,m)=(1,1)$ where machine precision is reached around $N=50$. In this particular example, the error saturates at $\sim {10}^{-6}$ because the residual field ${\overline{\mathrm{\Phi }}}_{\ell ,m}^{\mathcal{R}}$ assumes values $|{\mathrm{\Phi }}_{1,1}^{\mathcal{R}}|\sim {10}^{-13}$, i.e., the around roundoff limits of double float operations.