Worldtube excision method for intermediate-mass-ratio inspirals: Scalar-field model in $3+1$ dimensions

Binary black hole simulations become increasingly more computationally expensive with smaller mass ratios, partly because of the longer evolution time, and partly because the lengthscale disparity dictates smaller time steps. The program initiated by Dhesi et al. [Phys. Rev. D 104, 124002 (2021)] explores a method for alleviating the scale disparity in simulations with mass ratios in the intermediate astrophysical range (${10}^{-4}\lesssim q\lesssim {10}^{-2}$), where purely perturbative methods may not be adequate. A region (“worldtube”) much larger than the small black hole is excised from the numerical domain, and replaced with an analytical model approximating a tidally deformed black hole. Here we apply this idea to a toy model of a scalar charge in a fixed circular geodesic orbit around a Schwarzschild black hole, solving for the massless Klein-Gordon field. This is a first implementation of the worldtube excision method in full $3+1$ dimensions. We demonstrate the accuracy and efficiency of the method, and discuss the steps toward applying it for evolving orbits and, ultimately, in the binary black-hole scenario. Our implementation is publicly accessible in the spectre numerical relativity code.

Figure 1

Illustration of the computational domain: Shown is the equatorial plane, with height-deformation proportional to the value of the scalar field. The grid lines correspond to the DG-element boundaries of the 3-D numerical evolution. The central blue/green peak represents the region inside the worldtube, where the approximate solution is dominated by the singularity of the scalar field at the point-charge. Left of the peak an excision region is cut out within the horizon of the central black hole. A zoomed-out view of the entire domain is shown in Fig. 2.

Figure 2

The equatorial plane of the domain, depicting the steady-state solution of the scalar field $\mathrm{\Psi }$. The scalar charge creates an outward propagating spiral as it orbits the central black hole. Figure 1 shows a tilted perspective zoomed into the center of the same plane with the spiral arms visible in the background.

Figure 3

The steady-state solution of the scalar field ${\mathrm{\Psi }}^{\mathcal{N}}$ along the comoving $x$-axis and $z$-axis of the domain.

Figure 4

Top panel: the relative error of the regular field at the position of the charge compared to the value computed in [59]. Each cross represents the settled error at the final simulation time. The dashed lines are best fits for the relation $ϵ\propto R^{\alpha }$. Bottom panel: the local convergence order between simulations of neighboring worldtube radii.

Figure 5

Top panel: the absolute difference between the radial and angular derivatives of the regular field at the position of the particle and the reference values of Eqs. (51) and (52). Each cross represents the final value of a simulation. The dashed lines are best fits to the power-law relation $\propto R^{\alpha }$. Bottom panel: the local convergence order between the simulations of adjacent worldtube radii as defined in Eq. (48). The simulations at order $n=1$ show a convergence order consistent with the predicted rate $\alpha =1$. The $n=2$ simulations show a higher convergence rate likely due to dominant higher order terms. Some anomalies are not visible.

Figure 6

The relative error of the scalar field $\mathrm{\Psi }$ along the $z$-axis compared to the analytical solution ${\mathrm{\Psi }}_z$ given by Eq. (54). We show two simulations with worldtube radii $1.6M$ and $0.6M$ for each order, 0, 1 and 2. The error decreases with higher order or smaller worldtube radius, as expected. A small region is cut out around $z=r_p=5M$, where Eq. (54) converges too slowly to be calculated to sufficient accuracy in practice.

Figure 7

Top panel: the relative difference between the numerical, retarded field ${\mathrm{\Psi }}^{\mathcal{N}}$ along the $z$-axis and the analytical solution ${\mathrm{\Psi }}_z$ given in Eq. (54), integrated between $z=10M$ and $z=100M$ using the $L_1$-norm of Eq. (55). Each cross represents the final value of a simulation. The straight lines are a best fit to the power-law relation $ϵ\propto R^{\alpha }$. Bottom panel: the local convergence order between the simulations of adjacent worldtube radii as defined in Eq. (48). The simulations with $n=1$ and $n=2$ show a constant convergence order consistent with the predicted rate $\alpha =n+2$. The $n=0$ simulations show a higher convergence rate at larger worldtube radii but approach the expected value at smaller radii.

Figure 8

The relative error along the $x$-axis between $x=10M$ and $x=100M$ for two sample simulations of each order. The worldtube is centered at $x=5M$ and cut out from the plots.

Figure 9

Top panel: the relative error $ϵ$ integrated along the comoving x-axis compared to a reference solution ${\mathrm{\Psi }}_{\mathrm{ref}}$ with $n=2$ and $R=0.4M$. Each cross represents a simulation, and the straight lines are a best fit of the relation $ϵ\propto R^{\alpha }$. Bottom panel: the local convergence order as defined in Eq. (48). At first and second order, the simulations reproduce the expected convergence order of $\alpha =n+2$. At zeroth order, the local convergence approaches the expected order for smaller worldtube radii, likely indicating that the error still has contributions from higher-order terms in this regime.

Figure 10

Time steps $\mathrm{\Delta }t$ of the worldtube simulations presented here.