Deformation of horizons during a merger

TL;DR

This paper models the deformation of black hole horizons during an extreme mass ratio merger by numerically evolving perturbations and analyzing the horizon geometry, providing insights into horizon dynamics during such events.

Contribution

It introduces a numerical approach using the Zerilli formalism to study horizon deformation in EMR mergers, including horizon geometry analysis.

Findings

Horizon deformation patterns during merger are characterized.

Apparent horizon geometry is analyzed in detail.

Numerical evolution of perturbations is successfully implemented.

Abstract

We model an extreme mass ratio merger (EMR) as a point particle radially plunging into a large Schwarzschild black hole. We assume that the mass of the point particle, $\mu$, is much smaller than the black hole mass M. Under this assumption we can employ the Zerilli formalism modified to include a source term which arises from the energy-momentum tensor of the small object. We solve the Zerilli equation by numerically evolving initial data. Then, we ray trace the null geodesics of the event horizon from after the merger backward in time to extract the geometry of the perturbed event horizon. Further, we take advantage of the axisymmetry of the setup to locate the apparent horizon and study its geometry.

Equations

$$p_{\mu \nu }^{l0}d{x}^{\mu }d{x}^{\nu }=Y^{l0}(\theta )\left(f{H}_0^{l0}d{t}^2+2H_1^{l0}dtdr+\frac{1}{f}{H}_2^{l0}d{r}^2+r^2{K}^{l0}d{\Omega }^2\right)$$

$$\left[-\frac{{\partial }^2}{\partial t^2}+\frac{{\partial }^2}{\partial r^{\ast 2}}-V(r)\right]\psi (r,t)=S(r,t)$$

$$\psi (r,t)=\frac{r}{\lambda +1}\left\{K+\frac{f}{\Lambda }\left[H_2-r\frac{\partial }{\partial r}K\right]\right\}$$

$$V(r)=\frac{2f}{r^2{\lambda }^2}\left[{\lambda }^2(\Lambda +1)+\frac{9M^2}{r^2}\left(\Lambda -\frac{2M}{r}\right)\right]$$

$$S(r,t)=$$

$$\displaystyle-\frac{f^{2}}{\Lambda r}\left[\lambda(\lambda-1)r^{2}+(4\lambda-9)Mr+15M^{2}\right]Q^{tt}+r(\Lambda-f)Q^{rr}\Bigg{\}}$$

$$Q^{ab}=8\pi \int T^{ab}{Y}^{\ast }(\theta ,\varphi )d\Omega$$

$$u^{\nu }=(\tilde{E}/f,-({\tilde{E}}^2-f{)}^{1/2},0,0)$$

$$Q^{ab}=\mu \frac{8\pi }{r^2}\frac{u^a{u}^b}{u^t}\delta (r-r_p(t))Y^{\ast }(\theta (t),\varphi (t)).$$

$$S(r,t)=16\pi \frac{\sqrt{(2l+1)/4\pi }}{l(l+1)}\frac{\mu }{\tilde{E}}\frac{f^3}{\Lambda }\left\{{\delta }^{\prime }(r-r_p(t))-\left[\frac{(\lambda +1)r-3M}{r^{2}f}-\frac{6M{\tilde{E}}^2}{\Lambda r^{2}f}\right]\delta (r-r_p(t))\right\}$$

$$\frac{t}{2M}=\tilde{E}\sqrt{1-\frac{r}{r_0}}\sqrt{\frac{r}{2M}}\frac{r_0}{2M}+2\text{arctanh}\left(\frac{\sqrt{\frac{2M}{r}-\frac{2M}{r_0}}}{\tilde{E}}\right)+\tilde{E}\left(1+\frac{4M}{r_0}\right){\left(\frac{r_0}{2M}\right)}^{3/2}\text{arctan}\left(\sqrt{\frac{r_0}{r}-1}\right)$$

$${\psi }_N=-{\psi }_S+\left({\psi }_W+{\psi }_E\right)\left(1-\frac{{\Delta }^2}{2}{V}_S\right)$$

$${\psi }_N=$$

$$+{\psi }_W\left[1-\frac{V_S}{4}(A_1+A_3)\right]-\frac{1}{4}\left(1-\frac{V_b}{4}{A}_3\right)\iint dAS(r,t)$$

$$\iint dAS(r,t)=$$

$$\pm \kappa \left\{\frac{f(t_b)}{\Lambda (t_b)}[1\mp {\dot{r}}_p^{\ast }(t_b){]}^{-1}+\frac{f(t_t)}{\Lambda (t_t)}[1\pm {\dot{r}}_p^{\ast }(t_t){]}^{-1}\right\}$$

$$H_2=K=2\mu \frac{\sqrt{4\pi /(2l+1)}}{(1+M/2{\bar{r}}_0)(1+M/2\bar{r})}\frac{{\bar{r}}_{<}^l}{{\bar{r}}_{>}^{l+1}}$$

$$U=-e^{-\kappa u}\text{and}V=e^{\kappa v}.$$

$$\frac{d^2{X}^{\mu }}{d{V}^2}=-{\hat{\Gamma }}_{\sigma \gamma }^{\mu }\frac{d{X}^{\sigma }}{dV}\frac{d{X}^{\gamma }}{dV}+g\frac{d{X}^{\mu }}{dV}$$

$$\frac{d^2\delta \theta }{d{V}^2}=-\delta {\hat{\Gamma }}_{VV}^{\theta }$$

$$\frac{d\delta U}{dV}=e{\kappa }^2\delta g_{VV}$$

$$\left(\frac{\partial }{\partial v}-\kappa \right)\frac{\partial \delta {\theta }^l}{\partial v}=\kappa \left(\frac{\partial }{\partial v}-\kappa \right)\frac{\partial {\psi }^l}{\partial v}{\partial }_{\theta }{Y}^l$$

$$\left(\frac{\partial }{\partial v}-\kappa \right)\delta r^l=-\frac{1}{2\kappa }\left(\frac{\partial }{\partial v}-\kappa \right)\frac{\partial {\psi }^l}{\partial v}{Y}^l$$

$$\delta {\theta }^l=\kappa {\psi }^l{\partial }_{\theta }{Y}^{l0}+\frac{1}{\kappa }{\theta }_0^l{e}^{\kappa v}+{\theta }_1^l$$

$$\delta r^l=-\frac{1}{4\kappa }{\partial }_v{\psi }^l{Y}^{l0}+R_0^l{e}^{kv}$$

$$R_0^l$$

$${\theta }_0^l$$

$${\theta }_1^l$$

$$\Delta {\psi }^l={\psi }^l{|}_{v_{+}}-{\psi }^l{|}_{v_{-}}\Delta {\psi }_{,v}^l={\psi }_{,v}^l{|}_{v_{+}}-{\psi }_{,v}^l{|}_{v_{-}}$$

$$\left(\frac{\partial }{\partial v}-\kappa \right)\delta r^0=-\frac{\mu }{2M}$$

$${\tilde{n}}_{\mu }=\left(-\sqrt{f}+\frac{1}{2\sqrt{f}}{p}_{tt}\right){\delta }_{\mu }^t$$

$${\tilde{k}}^{\mu }=\frac{1}{\sqrt{2}}\left({\tilde{n}}^{\mu }+s^{\mu }\right)\text{and}{\tilde{l}}^{\mu }=\frac{1}{\sqrt{2}}\left({\tilde{n}}^{\mu }-s^{\mu }\right)$$

$$\tilde{\tau }(r,\theta )=r-\tilde{h}(\theta )$$

Full text

I Introduction

From the gravitational radiation observed by LIGO’s interferometers we now have strong evidence pointing towards the existence of binary black hole systems and their inevitable mergers Abbott:2016 . To predict the waveforms of the radiation emitted by these mergers complex simulations are executed on super computers using state-of-art algorithms. These simulations can accommodate a wide range of the possible parameters of the binary system, for e.g., initial spins, spin orientations, mass ratios, initial separation and velocities, etc.

In a very specific type of merger the calculations can be handled by much simpler methods. This happens when we consider the head-on merger of two non-spinning black holes which have an extreme mass ratio (EMR): that is the mass, $\mu$, of one of the black holes is much smaller than the mass, $M$, of the other, i.e. $\mu/M\rightarrow 0$. Then focusing on the large black hole, one can model the smaller black hole as a point particle with mass $\mu$. Hence, the problem is reduced to the radial in-fall of a point particle into a Schwarzschild black hole.

This approach to the problem was pioneered by Zerilli Zerilli:1971 and has since been studied extensively Martel:2001 ; Spallicci:2010 ; Lousto:1997 . At the heart of the solution is the assumption that the gravitational field produced by the particle drops off rapidly and so only perturbs the surrounding geometry. Further, there is no back-reaction on the particle and thus it follows a geodesic in the spacetime of the large Schwarzschild black hole. Under these assumptions Einstein’s equations can be linearized.

Here we will study how, in this perturbed spacetime, the event and apparent horizons of the large black hole are deformed by the gravitational effects of the stress-energy tensor of the point particle. For the deformation of the event horizon, a prominent role is played by caustics: points where null rays in the vicinity of the large black hole’s horizon cross each other and join the horizon. Recently, similar results have been shown in Emparan:2016 . There the strategy was to take the limit $M\rightarrow\infty$ and model the large black hole as a Rindler-type planar horizon accelerating towards the small black hole, which is assumed to be Schwarzschild. Then as the small black hole approaches, the null generators of the planar horizon are deflected towards the small black hole due to its Schwarzschild geometry and as such the deformation of the planar event horizon and the smaller hole’s event horizon can be computed. Another study that shows a similar deformation of the event horizon is Hamerly:2010 , where the approach is the one we take: approximating the small black hole as a point particle and solving for the perturbed geometry. In that paper the influence of the particle is treated as an impulse in the frequency domain of the perturbations.

By contrast, in this study we will do a direct numerical evolution, in the time domain, of Brill-Linquist type initial data for binary black holes. As in previous studies, we find that the event horizon of the large black hole “reaches out” in a cusp towards infalling particle. By contrast, for the apparent horizon we find that it appears to recede away from the approaching particle. While this is at first counter-intuitive, a little consideration shows that this isn’t so strange. Intuitively, as the particle gets closer there are null geodesics inside the black hole, close to the horizon, which get pulled outwards by the particle’s gravitational field. Hence the zero expansion surface moves inwards. Similar behaviour has been observed examining the marginally outer trapped surfaces (MOTS) in slices of initial data for arbitrary mass ratio head-on collisions Mosta:2015 .

The paper is organized as follows. Section II reviews the perturbation formalism for the radial in-fall of a point particle, Section III outlines the numerical method used to evolve the initial data, Section IV shows how the event and apparent horizons are deformed, and finally we conclude with Section V. Appendix A examines gauge choices for the $\l=1$ modes, Appendix B examines the physical effects of one of the approximations that we make in the course of our calculations and Appendix C tests numerical convergence of our schemes.

II Perturbations by radially in-falling point particle

We assume the particle creates a small disturbance in the spacetime which can be modelled by standard perturbation theory methods. We use the approach of Zerilli Zerilli:1971 with the Moncrief wavefunction Moncrief:1974 and express the metric as $\tilde{g}_{\mu\nu}=g_{\mu\nu}+p_{\mu\nu}$, where $g_{\mu\nu}$ is the Schwarzschild metric and $p_{\mu\nu}$ is the perturbation. Given the spherical symmetry of the Schwarzchild spacetime, $p_{\mu\nu}$ can be expanded in terms of odd and even modes with indices $l$ and $m$ arising from the spherical harmonics. If we assume the particle falls along the $z$-axis, then we can see that the problem is axisymmetric and only even modes are excited Martel:2001 . We write the perturbation as,

[TABLE]

where we have taken $m=0$ again because of axisymmetry, $f(r)=1-2M/r$ and the $H_{0},H_{1},H_{2}$ and $K$ are functions of $r$ and $t$, we omit the $l0$ where there is no chance of confusion. We work in the Regge-Wheeler gauge and Schwarzschild coordinates.

Then invoking Einstein’s equations up to linear order in $p_{\mu\nu}$, we find that the equations can be decoupled into one inhomogeneous wave equation:

[TABLE]

where $r^{*}=r+2M\text{ln}(r/2M-1)$ is the tortoise coordinate,

[TABLE]

is the Zerilli-Moncrief function,

[TABLE]

is the even potential with $\lambda=(l+2)(l-1)/2$, $\Lambda=\lambda+\frac{3M}{r}$, and

[TABLE]

is the source term for radial infall. The tensor $Q^{ab}$ where $a,b$ run over $r,t$ is given by,

[TABLE]

where $Y^{*}(\theta,\phi)$ is the complex conjugate of a spherical harmonic (notice the omitted $l$), $T^{ab}=\mu\int u^{a}u^{b}\delta^{4}(x^{\nu}-x_{p}^{\nu}(\tau))\sqrt{-g}d\tau$, is the stress tensor of the particle, $\mu$ is the mass of the particle, $\tau$ is the proper time along the particle’s trajectory, $x_{p}(\tau)$, and the four-velocity of the particle is,

[TABLE]

where $\tilde{E}$ is the conserved energy. Thus

[TABLE]

Here $Y^{*}(\theta(t),\phi(t))$ means that the spherical harmonic is evaluated along the trajectory of the particle, since the particle moves radially and its angular coordinates do not change, this is a constant. Inserting (9) into (5)

[TABLE]

where $r_{p}(t)$ is the trajectory of the particle and $Y^{*l0}(0)=\sqrt{(2l+1)/4\pi}$. The equation for the trajectory starting from rest is,

[TABLE]

where $r_{0}$ is the initial position of the particle and $\tilde{E}=\sqrt{1-\frac{2M}{r_{0}}}$.

III Numerical method

We use the numerical method of Lousto and Price Lousto:1997 with the modification by Martel and Poisson Martel:2001 . The method is a finite difference algorithm on a staggered grid with a step size $\Delta=\Delta r^{*}/4M=\Delta t/2M\approx 0.009$ over a domain bounded by two null hypersurfaces, $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, and one spacelike surface $\Sigma_{0}$. $\mathcal{H}_{1}$ approximates the event horizon and is given by, $u=t_{f}-r^{*}_{p}(t_{f})$ where $t_{f}$ is determined by (11) with $r_{p}(t_{f})/2M=1.00001$. $\mathcal{H}_{2}$ is a null hypersurface which approximates future null infinity and is chosen to be, $v/2M=(t+r^{*})/2M\approx 1500$. Finally $\Sigma_{0}$ is the $t=0$ surface which is the moment of time symmetry for the initial data.

After discretizing the domain we encounter two kinds of cells, as shown in FIG. 1; ones that do not contain the particle trajectory (Type I) and ones that that do (Type II). For Type I cells we do not have a contribution from the source term, $S(r,t)$. So we may discretize (2) with the following formula,

[TABLE]

where $S=(r^{*}_{S},t_{S})$ is the base point, $W=(r^{*}_{S}-\Delta,t_{S}+\Delta)$, $E=(r^{*}_{S}+\Delta,t_{S}+\Delta)$ and $N=(r^{*}_{S},t_{S}+2\Delta)$. This scheme is accurate to $\mathcal{O}(\Delta^{4})$.

For Type II cells there is a contribution from the source term and (2) now takes the form111There is a typo in Martel:2001 for equations (III) and (III), see Spallicci:2010 .

[TABLE]

where the $A_{i}$’s are the areas illustrated in right panel of FIG. 1 and $\iint dAS(r,t)$ is the integration of the source term over the cell,

[TABLE]

where, $l(l+1)\tilde{E}\kappa=16\pi\mu\sqrt{(2l+1/(4\pi)}$, $\dot{r}^{*}_{p}(t)=-\sqrt{\tilde{E}^{2}-f(t)}/\tilde{E}$, $f(t)=f(r_{p}(t))$, $\Lambda(t)=\Lambda(r_{p}(t))$, $t_{b}$ $(t_{t})$ is the time when the particle enters (leaves) the cell. Further, in the expression on the second line, the upper (lower) sign for the first term (function of $t_{b}$) is used when the particle enters the cell on the right (left) of $r^{*}_{S}$. For the second term (function of $t_{t}$) the upper (lower) sign is used when the particle leaves the cell on the right (left) of $r^{*}_{S}$.

The initial data at $t=0$ is constructed from the Brill-Lindquist Brill:1963 solution and is therefore conformally flat. It has been analysed in detail in Martel:2001 ; Lousto:1997 ; Lousto:1996 , and is given by,

[TABLE]

where $\bar{r}=r(1+\sqrt{f})^{2}/4$ is the isotropic radius, and $\bar{r}_{<}(\bar{r}_{>})$ is the smaller (greater) of $\bar{r}$ and $\bar{r}_{0}$. Note that with this choice of initial data, when we evolve from the $\Sigma_{0}$ surface, we have not specified what the $\psi_{S}$ points will be. We circumvent this problem by doing a Taylor expansion, $\psi(-\Delta,r^{*})=\psi(\Delta,r^{*})-2\Delta\partial_{t}\psi(t_{0},r^{*})+\mathcal{O}(\Delta^{3})$, since the data is time symmetric we can take $\psi(-\Delta,r^{*})=\psi(\Delta,r^{*})$. Since we will be interested in the horizon deformations, we will extract $\psi$ on the $\mathcal{H}_{1}$ surface. FIG. 2 shows the waveforms, $\psi(v)$, for $l=2$ to $l=25$ extracted on the surface $\mathcal{H}_{1}$ as a function of the ingoing coordinate $v$ with $r_{0}/2M=30$. This choice of $r_{0}$ makes $\tilde{E}\approx 1$. Notice that for large $l$ we get a smaller amplitude. Hence, we introduce a cutoff, $l_{\text{cut}}$, at $l_{\text{cut}}=25$. The (small) physical effect of this approximation is discussed in Appendix B.

In Appendix C we have shown how the code converges, where the method used to show this is the same as in Martel:2001 . Further, as a small test of our code we present the results of a simulation with a set of parameters chosen from Martel:2001 . Specifically, we extract the $l=2$ mode for $\psi$ at the surface that approximates null infinity, $\mathcal{H}_{2}$. We do this for a particle which starts from rest with a starting position of $r_{0}/2M=40$. In Martel:2001 the result with these parameters are presented in FIG. 3 (first plot, with $\alpha=1$), here it is presented in FIG. 8. By visual inspection these curves are virtually identical.

IV Horizons deformation

IV.1 Event Horizon

To find the deformation of the event horizon we use a method similar to the one employed in Hamerly:2010 , which calculates how each null generator of the horizon gets perturbed from its usual Schwarzschild position. For another method which locates the whole event horizon as a null surface, see Anninos:1994 . Consider the null generators of the event horizon of the static Schwarzschild spacetime in Kruskal-Szekeres coordinates, in terms of $u=t-r*$ and $v=t+r*$, these coordinates are defined as,

[TABLE]

The geodesic equation for the generators is given by,

[TABLE]

where $X^{0}=V(V),X^{1}=U(V)=0,X^{3}=\theta(V),X^{4}=\phi(V)$, $\hat{\Gamma}^{\mu}_{\sigma\gamma}$ denotes the Christoffel symbols in the Kruskal-Szekeres coordinates and we have kept the non-affine parameter $g$; for the static case we have $g=0$. In our case there is a perturbation present and so we add a perturbative term to each of the $X^{i}$ where $i$ runs over $U$ and $\theta$. We have then, $X^{i}\rightarrow X^{i}+\delta X^{i}$. Also, the Christoffel symbols will be perturbed, $\hat{\Gamma}^{\mu}_{\sigma\gamma}\rightarrow\hat{\Gamma}^{\mu}_{\sigma\gamma}+\delta\hat{\Gamma}^{\mu}_{\sigma\gamma}$. Inserting this ansatz into (17) and appropriately evaluating on the horizon we get for $i=\theta$;

[TABLE]

for $\delta U$ we may use the null condition for the generators,

[TABLE]

Switching to ingoing coordinates and expressing the perturbed quantities in terms of $\psi$, we can write (18) and (19) as,

[TABLE]

[TABLE]

where $l$ is the mode number for the spherical harmonics. The general solution is given by,

[TABLE]

[TABLE]

where $R_{0}^{l},\theta_{0}$ and $\theta_{1}$ are constants.

To find the horizon deformation we need to evolve the equations ‘backwards in time’. We start from the null generator of a static black hole at $U=0$ and $V\rightarrow\infty$ and integrate backwards to extract the deformation of the horizon. Hence, we impose the requirement that $\delta\theta$ and $\delta r$ vanish for $v\gg v_{m}$ where $v_{m}$ is time when the particle and black hole meet. Since at late times the perturbation vanishes we choose $R_{0}^{l}=\theta_{0}=\theta_{1}=0$. For $v=v_{m}$, note that there is a discontinuity in $\psi$ (cf. FIG. 2). To get a continuous deformation we choose:

[TABLE]

where,

[TABLE]

where $v_{+(-)}$ are the limits from the right (left) of $v_{m}$.

The effect of the monopole term, $l=0$, is calculated in the following manner: since the monopole term is isotropic, we use (19) for $v>v_{m}$ and arrive at,

[TABLE]

which has a general solution $\delta r^{0}=2\mu+C_{0}e^{\kappa v}$, where $C_{0}$ is a constant. For $v<v_{m}$ the right hand side of (28) will be zero and as such we have $\delta r^{0}=C_{0}e^{\kappa v}$. We fix $C_{0}=0$ for $v>v_{m}$ and for continuity at $v_{m}$ we fix $C_{0}=2\mu e^{\kappa(v-v_{m})}$ for $v<v_{m}$.

The $l=1$ mode is a little more complicated. Intuitively it corresponds to linear motions of the black hole and translations of the coordinate system. Ideally one would like to pick the “centre-of-mass” gauge for which the system is stationary at infinity. Unfortunately that gauge turns out to be unsuitable, with key components of the metric diverging on the horizon222While this can be seen algebraically in the discussion of Appendix A as of now we don’t have a good intuitive explanation for this result. We will return to this point in future work.. Hence we instead perform an infinitesimal gauge transformation to go into the singular gauge detweiler2004low in which the $l=1$ contribution to horizon deformations vanishes. These points are discussed in more detail in Appendix A.

Event horizons generated by these simulations are shown in FIG. 3–6. In each of these the mass of the large black hole is $M=1/2$ while the mass of the particle is $\mu=0.015$ (FIG. 3), $\mu=0.025$ (FIG. 4,6) and $\mu=0.075$ (FIG. 5). In the first three figures we see how the horizon deforms prior to $v_{m}$. The blue line is the event horizon and black dot illustrates the location of the particle. The green dots are null geodesics that are in the vicinity of the black hole. Given that we are working in perturbation theory the smaller the ratio of $\mu/M$ the more accurate our results and so in particular FIG. 3 ($\mu/M=0.03$) should be the most accurate, FIG. 4 ($\mu/M=0.05$) should still be pretty good but in FIG. 5 ($\mu/M=0.15$) is leaving the regime of applicability for the approximation. This last case is included mainly to show the effects of the approximation breaking down.

Then several interesting effects can be observed:

1. (i)

Before the merger there is a “bubble” of geodesics initially hugging close to the horizon that gradually moves away before joining the event horizon through caustics Siino:1997ix . It appears that the last geodesic to join the horizon lies on the axis of symmetry in the top half of the bubble and the caustics end at the point when the particle crosses the event horizon. This is most clearly seen for $\mu/M=0.03$. The disappearance of the caustics appears to be slightly post-merger for $\mu/M=0.05$ however this is probably an effect of the perturbation theory beginning to fail as the lag is even more dramatic for $\mu/M=0.15$.

Similar behaviours for the “bubble” of geodesics destined to join the event horizon can be seen in full non-EMR numerical relativity Bohn:2016 as well as the exact Schwarzschild calculations of Emparan and Martínez emparan_email . 2. (ii)

As a consequence of these caustic points, the event horizon appears to be stretched towards the infalling particle. Similar effects can be seen in full numerical relativity Bohn:2016 , the exact Schwarzschild calculations Emparan:2016 and the Hamerly-Chen impulse approximation Hamerly:2010 . 3. (iii)

Post-merger (FIG. 6) the sharp caustic point rounds out and the event horizon returns to sphericity. Again similar effects can be seen in the aforementioned studies.

Physically these effects are all a result of the focussing effects of the infalling particle. Because of its presence, some near-horizon null geodesics that would have otherwise escaped are captured on the event horizon. The associated caustics occur away from the original $r=2M$ surface in the direction of the particle. After it crosses the event horizon, the particle can no longer focus external geodesics.

Before moving on it is worth reiterating the limitations of the approximations that we have made. First we are working in first order perturbation theory and so close to the particle, where fields become strong we cannot expect high accuracy. Then within the perturbation theory the point particle is further approximated by the $l=25$ cut-off. The effect of this is to replace the point particle with a shell of stress-energy of radius $r=r_{p}(t)$. As is discussed in Appendix B, this is strongly peaked at $\theta=0$ but is also non-zero (and sometimes even energy-condition violating) at other $\theta$.

Combined these approximations model the physical spacetime well away from the particle but poorly close to $r=r_{p}(t)$ and $\theta=0$. In particular a real point mass would be masked by its own horizon and this is something that we cannot see in our model. Given that we do not have a horizon around the point mass we also cannot see any merger of horizons like that seen in Emparan:2016 .

That said we can have a high degree of confidence in features seen after the particle crosses the event horizon. Keep in mind that the event horizon is anchored in the future not the past: null rays are traced back from the future. As such once the particle is safely inside the black hole the problems associated with the approximations similarly disappear. For small $\mu/M$ it is only in the approach where the various approximations can manifest themselves.

IV.2 Apparent Horizon

Next we turn to the apparent horizon333Here we use apparent horizon in its colloquial sense as the outermost MOTS for the large black hole. and see how it is deformed by the infalling particle. The method we use is a standard one Dennison:2006 ; baumgarte2010numerical modified for small perturbations from the background geometry. We foliate the spacetime into spacelike slices, $\Sigma$, with the normal vector,

[TABLE]

which is the Schwarzschild time with a normalization $\tilde{n}^{\mu}\tilde{n}_{\mu}=-1$. Let $S$ be a closed two-dimensional surface in $\Sigma$ and let $\tilde{s}^{\mu}$ be the spacelike normal to $S$, hence, $\tilde{s}^{\mu}\tilde{s}_{\mu}=1$. For all points on $S$ we have outward and inward pointing null geodesics whose tangents, denoted $\tilde{k}^{\mu}$ and $\tilde{l}^{\mu}$ respectively, can be written as,

[TABLE]

Then the metric, $\tilde{m}_{\mu\nu}$, on $S$ can be written as, $\tilde{m}_{\mu\nu}=\tilde{g}_{\mu\nu}+\tilde{k}_{\mu}\tilde{l}_{\nu}+\tilde{l}_{\mu}\tilde{k}_{\nu}$. Further, the expansion of the outgoing null rays is given by, $\Theta=\tilde{m}^{\mu\nu}\tilde{\nabla}_{\mu}\tilde{k_{\nu}}$. The apparent horizon is then defined to be a marginally outer trapped surface (MOTS) on which the expansion of the outward null geodesics vanishes, $\Theta=0$.

We can define $S$ as a level surface of some scalar function $\tilde{\tau}(x^{i})=0$, where the $x^{i}$ are coordinates on $\Sigma$. In the axisymmetric case, using Schwarzschild coordinates, we write,

[TABLE]

and a normal and tangent to $S$ are

[TABLE]

Then the condition for vanishing expansion can be written as

[TABLE]

where,

[TABLE]

$\tilde{\gamma}_{\mu\nu}$ is the induced metric on $\Sigma$, the indices $A,B$ run over $(r,\theta)$, and $\tilde{\gamma}^{(2)}$ is the determinant of $\tilde{\gamma}_{AB}$. Since our approach is perturbative to linear order we will use the ansatz $\tilde{h}(\theta)=2M+\delta h(\theta)$ and keep terms linear in $\delta$. Further evaluating $\psi(v)$ on $\mathcal{H}_{1}$, which allows us to switch to ingoing coordinates, (33) becomes,

[TABLE]

Notice that $\left(\partial_{\theta}^{2}+\cot\theta\partial_{\theta}\right)$ is the Laplacian on the unit sphere. To solve (37) we assume that $\delta h$ can be expanded in ($m=0$) spherical harmonics,

[TABLE]

inserting this into (37) we find,

[TABLE]

Inserting (39) into (38) we can construct the deformation of the apparent horizon. This is shown as the magenta curve in FIG.  3-5. For early times, when the particle is far from the black hole the apparent horizon coincides with $r=2M$. As for the event horizon the effects of the approaching particle are concentrated on the piece of the horizon facing the particle. However the apparent horizon behaves quite differently from the event horizon. While the event horizon puckered out towards the particle the apparent horizon recedes (first seen about $v=270.30$ in all simulations) and appears to become concave outwards as the particle comes closer.

Intuitively this behaviour can be understood in the following way. As the particle approaches the black hole, it bends light rays towards it. In particular this increases the outward expansion of congruences of null geodesics inside the $r=2M$ surface and thereby increases their (initially negative) expansion. Thus the new $\Theta=0$ surface is inside the old $r=2M$ surface. This recession does not violate the rule that apparent horizons cannot decrease in area Hayward:1993wb ; Ashtekar:2004cn ; Andersson:2005gq ; Bousso:2015qqa : our approximation is only first order and area change is a second order effect. To the order of accuracy of our approximation the area can’t change.

There is a gauge issue to consider here. Apparent horizons are creatures of their particular foliation and in general the time-evolved apparent horizon of one foliation will not match that from a different slicing of spacetime. Hence a proper physical understanding of the recession would include an investigation of its foliation dependence. That is beyond the scope of the current paper but we note that this is not the first time that such an effect has been observed. A similar effect was seen in the analysis of the MOTSs in initial data for head-on collisions in a non-perturbative setting Mosta:2015 444Referring back to the concerns of the previous paragraph, in that paper the total horizon area was seen to increase in spite of the recession from the smaller black hole.. Hence we expect this recession to be robust for reasonable choices of foliation.

Finally we note that $v=270.35$ is probably as far as we can expect our linear perturbation theory to return reasonable results for the apparent horizon. For example with ${\mu}/{M}=0.05$, by $v=270.35$ the apparent horizon of the small black hole (that we are modelling as a particle) is approaching that of the large black hole. Hence we are now in a region where the gravitational field that we have assumed to be small is becoming large. Indeed just after this point, our horizon finding methods fail with the apparent horizon diverging through and away from the event horizon.

This undesirable (and physically impossible!) behaviour results from trying to push the horizon through an area where our approximations are no longer good: the perturbations are no longer small and we are probing at a resolution where the $l=25$ cut-off becomes problematic. As was emphasized in Booth:2017fob , identifying the location of an apparent horizon strongly depends on local features of the geometry: there are many surfaces of vanishing null expansion running through all points of a spacetime and the only thing distinguishing a MOTS is that it is a closed surface. As such, once we lose control of those details in a region through which the horizon should pass, we can no longer locate it.

In our case, at $v=270.35$ the apparent horizon should close near the infalling particle. However our description of the geometry is not good in that region and so the closure does not happen in the expected way. By contrast the highly non-local event horizon calculations are more robust. Problems in one small region only affect a few geodesics and so the global evolution is not seriously perturbed.

V Discussion

In this paper we modelled an EMR merger by using perturbation techniques for Schwarzschild black holes. Specifically, we used the stress tensor of a point particle to model the gravitational field of a small black hole/point particle and used Zerilli’s formalism to solve the linearized Einstein’s equations. We then studied the deformation of the event and apparent horizons as the particle approaches the black hole (FIGS. 3–6).

Given the teleological nature of event horizons we used final post-merger boundary conditions and integrated backwards in time to get the full horizon. To ensure a smooth horizon evolution we fixed constants in (24) to (26), whereby we see that the discontinuity in the initial data plays a major role in the deformation of the event horizon. This discontinuity arises from doing a multipole expansion of the Brill-Lindquist type initial data Lousto:1996 . Such a dependence on initial data for perturbations is lost when analysing solutions of (2) in the frequency domain. We also saw how there is a “bubble” of null rays that emerges between the particle and the black hole just before the merger, shown in green in FIG. 3–4. These are null rays that join the event horizon at the top of the black hole by crossing each other and forming caustics. We noted that this bubble of null rays and associated caustics are a standard feature of black hole mergers and commonly seen in full non-linear treatments (for example Bohn:2016 ).

For locating the apparent horizon we used a standard approach (baumgarte2010numerical, ) simplified to the case of linear perturbations. Perhaps the most interesting observation was that as the particle makes its final approach to the black hole, the region of the apparent horizon facing the particle appears to recede away from it, ultimately becoming concave outwards in the diagrams. This apparently non-intuitive behaviour can be understood as resulting from null rays in the immediate vicinity of the unperturbed apparent horizon being attracted towards particle. Hence relative to the unperturbed metric the apparent horizon moves inwards. At the same time we noted that the recession doesn’t result in a decrease in the area of the horizon: in our first order expansion the area is invariant.

It is interesting to compare the quite different evolutions of the event versus apparent horizons. The event horizon is “attracted” towards the particle as one might naively expect, while the apparent horizon initially recedes. Further during the final approach the event horizon puckers to form a caustic cusp as new null generators join the horizon while the apparent horizon remains smooth as long as it can be found

The appearance of caustics is not a consequence of our point source. For one we have noted that our cut-off of higher $l$ modes “smears out” the particle and so despite our initial set-up we are not really dealing with a point source. However more generally caustics are a result of geodesic focussing and this does not require a point source. For example they are also seen in Emparan:2016 where the smaller black hole has a Schwarzschild geometry and is not a point source, as well as in the full non-linear treatment of mergers with arbitrary initial parameters Bohn:2016 ; Matzner941 .

The event horizon grows during the approach (thanks to its final boundary conditions) while the apparent horizon does not (thanks to the linear approximations). However this is realistic. During this physical process most of the mass change in the large black hole comes directly from the infalling particle. The event horizon grows “in anticipation” of the absorption however the apparent horizon should not substantially change in size until the particle actually is absorbed.

Overall we have found that for an EMR merger many key features can be captured by using standard perturbation methods. However some features escape us. Specifically any features related to the event and apparent horizon of the smaller black hole are lost, simply because we use a point particle approximation. Perhaps the matched asymptotic methods suggested in Emparan:2016 could be used to focus on the event horizon geometry of the smaller black hole. Further the appearance of a common apparent horizon Bishop1982 ; Cadez1973 ; Anninos:1994ay ; Cook:1992 that is usually seen in arbitrary mass ratio mergers is also lost. This feature in principle could also be recovered by upgrading the particle to a Schwarzschild geometry rather than a point stress tensor. In future papers we intend to return to this issue as well as considering more general interactions than the head-on collisions studied here.

Acknowledgements.

IB thanks Abhay Ashtekar and Aaron Zimmerman for a conversation at the 15th Canadian Conference on General Relativity and Relativistic Astrophysics which inspired this paper. UH would like to thank Hari Kunduri, Harald Pfeiffer and Aaron Zimmerman for their insights and suggestions. We would also like to thank Roberto Emparan for an enlightening correspondence following the posting of the first version of this paper to the arXiv which significantly influenced later versions and Eric Poisson for his valuable comments contributed as part of UH’s PhD defence. IB is supported by NSERC Discovery Grant 261429-2013. During the course of this research UH was partially funded by a fellowship from the School of Graduate Studies at Memorial University as well as stipends from NSERC Discovery Grants 261429-2013 and 418537-2012.

Appendix A Choices for the $l=1$ gauge

The $l=1$ perturbation is given by,

[TABLE]

where $r_{p}(t)$ is the trajectory, $\Theta$ is the Heaviside function and,

[TABLE]

Clearly, this is a gauge where the perturbation is non-zero outside a sphere of radius $r=r_{p}(t)$. As mentioned above, when computing the event horizon deformation we need to evaluate the right hand side of Eqs. (18) and (19) at $r=2M$. We find,

[TABLE]

where $f_{0}(v)\sim r_{p}(v)$ where $r_{p}(v)$ is the trajectory of the particle in Eddington-Finkelstein coordinates (see for example muller2008falling ). The term in the square brackets is finite and hence the whole quantity is undefined at $r=2M$. We can utilize the gauge transformation used by Zerilli in Zerilli:1971 to switch to a gauge where the perturbation is non-zero inside a sphere of radius $r=r_{p}(t)$. We find in this gauge,

[TABLE]

where the dots denote derivatives w.r.t. $t$. Again, we see a similar situation as above where the term in the square brackets is finite and the whole quantity is undefined at $r=2M$.

It is possible to go into a gauge where the perturbation is concentrated at the location of the particle in the form of a delta function. This gauge can be achieved by appropriately adding a Heaviside function to Zerilli’s gauge transformation. In this gauge we have,

[TABLE]

We get the contribution from the delta function when we integrate in Eq. (19) but notice that the quantity in square brackets in front of the delta function vanishes when $r_{p}(t)\rightarrow 2M$ and $r\rightarrow 2M$. This nulls the effect of the delta function. Choosing the trivial solution for Eq. (19) we see that in this gauge the $l=1$ perturbation does not effect the horizon deformation. The same reasoning applies to Eq (18) as $\delta\Gamma^{\theta}_{VV}\sim\partial_{\theta}\delta g_{VV}$. Further, since we compute the apparent horizon before the particle touches $r=2M$, the apparent horizon is not affected.

Appendix B Energy condition violations and $l$ cut-off

In this appendix we demonstrate that the introduction of an $l=25$ cut-off in the spherical harmonic expansion of the wave-form is physically equivalent to the introduction of a spherical matter field at $r=r_{p}(t)$. For $l=25$ the energy density is strongly peaked at $\theta=0$ however it is not a perfect point source. While we do not believe that this approximation significantly affects our simulation it is useful to keep in mind when thinking about the results of the simulations.

The stress tensor associated with a point particle is given by,

[TABLE]

where $\{a,b\}=\{r,t\}$, $u^{a}$ is the four velocity of the particle and $x_{p}^{\nu}$ is the position of the particle. This can be expanded in ($m=0$) spherical harmonics,

[TABLE]

where

[TABLE]

Further, it can be shown that,

[TABLE]

Now, since the particle falls along the $z$-axis, its $\theta$ coordinate does not change: $\theta_{p}=0$. This implies that $Y_{l}(0)=\sqrt{\frac{2l+1}{4\pi}}P_{l}(\cos(0))=\sqrt{\frac{2l+1}{4\pi}}$ for all $t$. We have then,

[TABLE]

If a finite cut-off is introduced in the $l$ series in Eq. 51 at $l_{\text{cut}}=25$, we get dependence on $\theta$ as shown in Fig. 7. Notice that there are regions where the sum is negative, this implies that there are regions where $T^{tt}_{\text{cut}}<0$. This means that there are small energy condition violations associated with introducing a finite $l$ cut off.

Such energy condition violations can, for example, cause the apparent horizon to be outside the event horizon in a dynamical setting (see, for example, Scheel:1994yn ). However they are very weak in our example and so relatively innocuous.

Appendix C Numerical convergence and waves at infinity

To test our implementation of the numerical algorithm we make comparisons with the results of Martel:2001 . Specifically, we reproduce the result presented in FIG. 3 of Martel:2001 where the particle falls in from rest with a starting position of $r_{0}/2M=40$ and the $l=2$ master function is evaluated at null infinity, $\mathcal{H}_{2}$. This is shown in FIG. 8, where it can be seen that we get the same result (for $\alpha=1$), we get identical results for other values $l$ and $r_{0}$ also.

We also demonstrate quadratic convergence of the code in an identical manner to the Appendix of Martel:2001 by computing $\delta\psi$, as defined in the Appendix of Martel:2001 . For completeness the definition of $\delta\psi$ is reviewed here. The numerical algorithm is constructed to be convergent at second order, hence we expect,

[TABLE]

where $\psi_{N}$ is the numerical solution, $\psi_{\text{exact}}$ is the exact solution, $\Delta$ is step size and $\rho$ is an error function that is independent of $\Delta$. We may then define,

[TABLE]

Notice that we have $\delta\psi(n\Delta)=n^{2}\delta\psi(\Delta)$. We calculate $\delta\psi(\Delta)$ at a resolution of $\Delta=0.02$ and $\Delta=0.04$. In FIG. 9 the black curve is $\delta\psi$ with a resolution $\Delta=0.04$ and the grey curve is $4\delta\psi$ with a resolution of $\Delta=0.02$, the agreement of these curves indicates quadratic convergence.