# Fourth order indirect integration method for black hole perturbations: even modes

Patxi Ritter, Alessandro D. A. M. Spallicci, Sofiane Aoudia, Stéphane Cordier

► To cite this version:

Patxi Ritter, Alessandro D. A. M. Spallicci, Sofiane Aoudia, Stéphane Cordier. Fourth order indirect integration method for black hole perturbations: even modes. Conference on Theory Meets Data Analysis at Comparable and Extreme Mass Ratios, Jun 2010, Waterloo, Canada. hal-00565872

HAL Id: hal-00565872

https://hal.science/hal-00565872v1

Submitted on 5 Aug 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Fourth order indirect integration method for black hole perturbations: even modes

Patxi Ritter¹,² Alessandro D.A.M. Spallicci¹‡
Sofiane Aoudia³ Stéphane Cordier²

¹ Université d'Orléans, Observatoire des Sciences de l'Univers en région Centre, LPC2E Campus CNRS, 3A Av. Recherche Scientifique, 45071 Orléans, France

² Université d'Orléans, Laboratoire de Mathématiques - Analyse, Probabilités, Modélisation - Orléans, MAPMO, Rue de Chartres, 45067 Orléans, France

³ Max Planck Institut für Gravitationphysik, A. Einstein, Am Mühlenberg 1, 14476 Potsdam, Deutschland

**Abstract.** On the basis of a recently proposed strategy of finite element integration in time domain for partial differential equations with a singular source term, we present a fourth order algorithm for non-rotating black hole perturbations in the Regge-Wheeler gauge. Herein, we address even perturbations induced by a particle plunging in. The forward time value at the upper node of the (r*, t) grid cell is obtained by an algebraic sum of i) the preceding node values of the same cell, ii) analytic expressions, related to the jump conditions on the wave function and its derivatives, iii) the values of the wave function at adjacent cells. In this approach, the numerical integration does not deal with the source and potential terms directly, for cells crossed by the particle world line. This scheme has also been applied to circular and eccentric orbits and it will be object of a forthcoming publication.

PACS numbers: 04.25.Nx, 04.30.Db, 04.30.Nk, 04.70.Bw, 95.30.Sf

# 1. Introduction

In the scenario of the capture of compact objects by a supermassive black hole of mass *M*, the seized object is compared to a small mass *m* (henceforth the particle or the source) perturbing the background spacetime curvature and generating gravitational radiation. A comprehensive introduction to the general relativistic issues related to EMRI (Extreme Mass Ratio Inspiral) sources is contained in a topical volume [1].

Schwarzschild-Droste (SD) [2–5] (see Rothman [6] for a justification of this terminology), black hole perturbations have been hugely developed in the Regge-Wheeler (RW) gauge, before in vacuum [7] and after in the presence of a particle by Zerilli [8–11]. The first finite difference scheme in time domain was proposed by Lousto and Price [12]. The initial conditions, reflecting the past motion of the particle and the initial amount of gravitational waves, were parametrised by Martel and Poisson [13].

If the gravitational radiation emitted and the mass of the captured object are to be taken into account for the determination of the motion of the latter, it is necessary to compute the derivatives of the perturbations that implies the third derivative of the wave function $\Psi(r^*, t)$, see e.g. [14]. For a given accuracy $O(h)$ of the third derivative

‡ Corresponding author: spallicci@cnrs-orleans.fr
of $\Psi$, the error on $\Psi$ itself should be $O(h^4)$. Effectively, the reminder ought to be $O(h^5)$ due to the presence in the mesh of the particle that lowers by one more degree the convergence order of the code for geometrical effects [15]. We have therefore developed a fourth order scheme.

The complexity in assessing the continuity of the perturbations at the position of the particle and the compatibility of the self-force to the harmonic (Lorenz-de Donder§) gauge [16, 17] has led researchers to convey their efforts to this gauge, as commenced by Barack and Lousto [19]. Conversely, work in harmonic gauge is made cumbersome by the presence of a system of ten coupled equations which replace the single wave equation of the RW gauge.

We have proposed [20, 21] a finite element method of integration, in RW gauge, based on the jump conditions that the wave function and its derivatives have to satisfy for the SD black hole perturbations to be continuous at the position of the particle. We first deal with the radial trajectory and the associated even parity perturbations, while in a forthcoming paper we shall present the circular and eccentric orbital cases, referring thus to both odd and even parity perturbations.

The main feature of this method consists in avoiding the direct and explicit integration of the wave equation (the potential and the source term with the associated singularities) whenever the grid cells are crossed by the particle. Indeed, the information on the wave equation is implicitly given by the jump conditions on the wave function and its derivatives. Conversely, for cells not crossed by the particle world line, the integrating method might retain the previous approach by Lousto [22] and Haas [15]. Among the efforts using jump discontinuities, although in a different context, it is worthwhile to mention those of Haas [15], Sopuerta and coworkers [23–25] getting the self-force in a scalar case. For the geodesic gravitational case, like Sopuerta and coworkers, Jung et al. [26], Chakraborty et al. [27] rely on spectral methods; Zumbusch [28], Field et al. [29] use a discontinuous Galerkin method; Hopper and Evans [30] work partially in frequency domain. Among recent results not based on jump discontinuities but concerning fourth order time domain codes, the one proposed by Thornburg [31] deals with and adaptive mesh refinement, while Nagar and coworkers replace the delta distribution with a narrow Gaussian [32, 33].

For the computation of the back-action, this method ensures a well behaved wave function at the particle position, since the approach is governed by the analytical values of the jump conditions at the particle position.

In [21] we have provided waveforms at infinity and the wave function at the position of the particle at first order. Herein, we focus instead on the improvement of the algorithm at fourth order and refer to [21] for all complementary information. The features of this method can be summarised as follows:

* Avoidance of direct and explicit integration of the wave equation (the potential and the source term with the associated singularities) for the grid cells crossed by the particle.

* Improvement of the reliability, since analytic expressions partly replace numerical ones (the replacement is total at first order [20, 21]).

* Applicability of the method to generic orbits, assuming that the even and odd wave equations are satisfied by $\Psi$, respectively $R$, being $C^{-1}$.

<footnote>§ FitzGerald is considered to have also identified the harmonic gauge [18].</footnote>

<footnote>‖ A $C^{-1}$ continuity class element, like a Heaviside step distribution, may be seen as an element which after integration transforms into an element belonging to the $C^0$ class of functions.</footnote>
Geometric units ($G = c = 1$) are used, unless stated otherwise. The metric signature is $(-, +, +, +)$.

## 2. The wave equation

The wave function (its dimension is such that the energy is proportional to $\int_0^\infty \dot{\Psi}^2 dt$), in the Moncrief form [34] and RW gauge [7], is defined by

$$ \Psi_l(t, r) = \frac{r}{\lambda + 1} \left[ K^l + \frac{r - 2M}{\lambda r + 3M} \left( H_2^l - r \frac{\partial K^l}{\partial r} \right) \right], \quad (1) $$

where $K(t, r)$ and $H_2(t, r)$ are the perturbations, and the Zerilli [9] normalisation is used for $\Psi_l$. The wave equation is given by the operator $\mathcal{Z}$ acting on the wave function

$$ \mathcal{Z}\Psi^l(t, r) = \partial_{r^*}^2 \Psi^l(t, r) - \partial_t^2 \Psi^l(t, r) - V^l(r)\Psi^l(t, r) = S^l(t, r), \quad (2) $$

where $r^* = r + 2M \ln(r/2M - 1)$ is the tortoise coordinate and the potential $V^l(r)$ is

$$ V^l(r) = \left(1 - \frac{2M}{r}\right) \frac{2\lambda^2(\lambda+1)r^3 + 6\lambda^2Mr^2 + 18\lambda M^2r + 18M^3}{r^3(\lambda r + 3M)^2}, \quad (3) $$

being $\lambda = 1/2(l - 1)(l + 2)$. The source $S^l(t, r)$ includes the derivative of the Dirac distribution (the latter appear in the process of forming the wave equation out of the ten linearised Einstein equations)

$$ S^l = \frac{2(r - 2M)\kappa}{r^2(\lambda + 1)(\lambda r + 3M)} \times \left\{ \frac{r(r - 2M)}{2U^0} \delta'[r - r_u(t)] - \left[ \frac{r(\lambda + 1) - 3M}{2U^0} - \frac{3MU^0(r - 2M)^2}{r(\lambda r + 3M)} \right] \delta[r - r_u(t)] \right\}, \quad (4) $$

$U^0 = E/(1 - 2M/r_u)$ being the time component of the 4-velocity, $E = \sqrt{1 - 2M/r_u}$ the conserved energy per unit mass, and $\kappa = 4m\sqrt{(2l + 1)\pi}$. The geodesic in the unperturbed SD metric $z_u(\tau) = \{t_u(\tau), r_u(\tau), \theta_u(\tau), \phi_u(\tau)\}$ assumes different forms according to the initial conditions. For radial infall of a particle starting from rest at finite distance $r_{u0}$, $r_u(t)$ is the inverse function in coordinate time $t$ of the trajectory in the background field, corresponding to the geodesic in proper time $\tau$ (u stands for unperturbed)

$$ \frac{t(r_u)}{2M} = \sqrt{1 - \frac{2M}{r_{u0}}} \sqrt{1 - \frac{r_u}{r_{u0}}} \left(\frac{r_{u0}}{2M}\right) \left(\frac{r_u}{2M}\right)^{1/2} + 2\operatorname{arctanh}\left(\frac{\sqrt{\frac{2M}{r_u} - \frac{2M}{r_{u0}}}}{\sqrt{1 - \frac{2M}{r_{u0}}}}\right) + \\ \sqrt{1 - \frac{2M}{r_{u0}}} \left(1 + \frac{4M}{r_{u0}}\right) \left(\frac{r_{u0}}{2M}\right)^{3/2} \operatorname{arctan}\left(\sqrt{\frac{r_{u0}}{r} - 1}\right). \quad (5) $$

The expressions above correspond to those in [14], where some of the errors of previously published literature on radial fall are indicated.
### 3. Jump conditions

From the visual inspection of the Zerilli wave equation (2), it is evinced that the wave function $\Psi$ is of $C^{-1}$ continuity class (the second derivative of the wave function is proportional to the first derivative of the Dirac distribution, in itself a $C^{-3}$ class element). Thus, the wave function and its derivatives can be written as (the $l$ index is dropped henceforth for simplicity of notation)

$$ \Psi = \Psi^{+} \Theta_{1} + \Psi^{-} \Theta_{2}, \qquad (6) $$

$$ \Psi_{,r} = \Psi_{,r}^{+} \Theta_{1} + \Psi_{,r}^{-} \Theta_{2} + (\Psi^{+} - \Psi^{-}) \delta, \qquad (7) $$

$$ \Psi_{,t} = \Psi_{,t}^{+} \Theta_{1} + \Psi_{,t}^{-} \Theta_{2} - \dot{r}_{u} (\Psi^{+} - \Psi^{-}) \delta, \qquad (8) $$

$$ \Psi_{,rr} = \Psi_{,rr}^{+} \Theta_{1} + \Psi_{,rr}^{-} \Theta_{2} + 2 (\Psi_{,r}^{+} - \Psi_{,r}^{-}) \delta + (\Psi^{+} - \Psi^{-}) \delta', \qquad (9) $$

$$ \Psi_{,tt} = \Psi_{,tt}^{+} \Theta_{1} + \Psi_{,tt}^{-} \Theta_{2} - 2 \dot{r}_{u} (\Psi_{,t}^{+} - \Psi_{,t}^{-}) \delta - \ddot{r}_{u} (\Psi^{+} - \Psi^{-}) \delta + \dot{r}_{u}^{2} (\Psi^{+} - \Psi^{-}) \delta', \quad (10) $$

$$ \Psi_{,tr} = \Psi_{,tr}^{+} \Theta_{1} + \Psi_{,tr}^{-} \Theta_{2} + (\Psi_{,t}^{+} - \Psi_{,t}^{-}) \delta - \dot{r}_{u} (\Psi_{,r}^{+} - \Psi_{,r}^{-}) \delta - \dot{r}_{u} (\Psi^{+} - \Psi^{-}) \delta', \quad (11) $$

where in shortened notation $\Theta_1 = \Theta[r - r_u(t)]$, and $\Theta_2 = \Theta[r_u(t) - r]$ are two Heaviside step distributions, while $\delta = \delta[r - r_u(t)]$ and $\delta' = \delta'[r - r_u(t)]$ are the Dirac delta - and its derivative - distributions. The dot and the prime indicate time and space derivatives, respectively.

#### 3.1. Jump conditions from the wave equation

For the computation of back-action effects, we need first order derivatives of the perturbations and thus third order wave function derivatives. To this end, we operate directly on the wave equation, Eq. 2. The source term is cast in the following form

$$ S(t, r) = G(t, r)\delta + F(t, r)\delta' = \tilde{G}_{r_u(t)}\delta + F_{r_u(t)}\delta', \qquad (12) $$

where $\tilde{G}_{r_u(t)} = G_{r_u(t)} - F'_{r_u(t)}$ and one of the properties of the Dirac delta distribution, namely $\phi(r)\delta'[r - r_u(t)] = \phi_{r_u(t)}\delta'[r - r_u(t)] - \phi'_{r_u(t)}\delta[r - r_u(t)]$, has been used at the position of the particle. The determination of the jump conditions imposes the transformation of Eq. 2 into the corresponding equation in (r,t) domain (the tortoise coordinate can't be inverted). Turning to the $r$ variable, we get ($f = 1 - 2M/r$)

$$ \begin{aligned} \partial_{r^*}^2 \Psi &= f f' \partial_r \Psi + f^2 \partial_r^2 \Psi \\ &= [f f' \Psi_{,r}^+ + f^2 \Psi_{,rr}^+] \Theta_1 + [f f' \Psi_{,r}^- + f^2 \Psi_{,rr}^-] \Theta_2 + f f' (\Psi^+ - \Psi^-) \delta \\ &\quad + 2 f^2 (\Psi_{,r}^+ - \Psi_{,r}^-) \delta + f^2 (\Psi^+ - \Psi^-) \delta', \end{aligned} \qquad (13) $$

$$ \begin{aligned} \partial_t^2 \Psi = & \Psi_{,tt}^{+} \Theta_1 + \Psi_{,tt}^{+} \Theta_2 - 2 \dot{r}_u \partial_t (\Psi^+ - \Psi^-) \delta - \ddot{r}_u (\Psi^+ - \Psi^-) \delta \\ & + \dot{r}_u^2 (\Psi^+ - \Psi^-) \delta', \end{aligned} \qquad (14) $$

$$ V \Psi = V \Psi^{+} \Theta_{1} + V \Psi^{-} \Theta_{2}. \qquad (15) $$

The notation $[\Psi]$ stands for the difference $(\Psi^+ - \Psi^-)_{r_u}$ and a likewise notation is used for the derivatives at the point $r_u$. Equating the coefficients of $\delta'$, and owing
to the above mentioned property of the delta derivative for which $(\Psi^+ - \Psi^-) \delta' = [\Psi] \delta' - [\Psi, r] \delta$, we get the jump condition for $\Psi$

$$[\Psi] = \frac{1}{f_{r_u}^2 - \dot{r}_u^2} F_{r_u}. \qquad (16)$$

Equating the coefficients of $\delta$, we get the jump condition on the space derivative

$$[\Psi, r] = \frac{1}{f_{r_u}^2 - \dot{r}_u^2} \left[ \tilde{G}_{r_u} + (f_{r_u} f'_{r_u} - \ddot{r}_u) [\Psi] - 2 \dot{r}_u \frac{d}{dr_u} [\Psi] \right], \quad (17)$$

and therefore the jump condition on the first time derivative

$$[\Psi, t] = \dot{r}_u \frac{d}{dr_u} [\Psi] - \dot{r}_u [\Psi, r]. \qquad (18)$$

Since $\mathcal{Z}\Psi^\pm = 0$, the coefficients of $\Theta_1$ and $\Theta_2$ ought to be equal. We thus obtain

$$[\Psi, tt] - f_{r_u} f'_{r_u} [\Psi, r] - f_{r_u}^2 [\Psi, rr] + V_{r_u} [\Psi] = 0, \qquad (19)$$

which is an equation with two unknowns. We circumvent the difficulty by using i) the commutativity of the derivatives, $[\Psi, tr] = [\Psi, rt]$, ii) the transformation $d/dt = \dot{r}_u d/dr_u$, and write

$$[\Psi, tt] = \frac{d}{dt} [\Psi, t] - \dot{r}_u [\Psi, tr] = \frac{d}{dt} [\Psi, t] - \dot{r}_u \left\{ \frac{d}{dt} [\Psi, r] - \dot{r}_u [\Psi, rr] \right\} = \\ = \dot{r}_u \frac{d}{dr_u} [\Psi, t] - \dot{r}_u^2 \frac{d}{dr_u} [\Psi, r] + \dot{r}_u^2 [\Psi, rr]. \qquad (20)$$

The jump condition on the second space derivative can now be expressed by

$$[\Psi, rr] = \frac{1}{f_{r_u}^2 - \dot{r}_u^2} \left\{ \dot{r}_u \frac{d}{dr_u} [\Psi, t] - \dot{r}_u^2 \frac{d}{dr_u} [\Psi, r] - f_{r_u} f'_{r_u} [\Psi, r] + V_{r_u} [\Psi] \right\}. \quad (21)$$

The other second derivatives are obtained by

$$[\Psi, tr] = [\Psi, rt] = \frac{d}{dt} [\Psi, r] - \dot{r}_u [\Psi, rr], \qquad (22)$$

$$[\Psi, tt] = \frac{d}{dt} [\Psi, t] - \dot{r}_u [\Psi, tr]. \qquad (23)$$

For the third order derivatives, we derive the wave equation with respect to $r$ and obtain

$$[\Psi, rrr] = \frac{1}{\dot{r}_u^2 - f_{r_u}^2} \left\{ \dot{r}_u^2 \frac{d}{dr_u} [\Psi, rr] - \dot{r}_u \frac{d}{dr_u} [\Psi, rt] + (f_{r_u}^{\prime 2} + f_{r_u} f_{r_u}^{\prime \prime} - V_{r_u}) [\Psi, r] + 3 f_{r_u} f_{r_u}^{\prime} [\Psi, rr] - V_{r_u}^{\prime} [\Psi] \right\}, \quad (24)$$

while deriving with respect to $t$, we obtain

$$[\Psi, ttt] = \frac{\dot{r}_u^2}{\dot{r}_u^2 - f_{r_u}^2} \left\{ f_{r_u}^2 \frac{d}{dr_u} [\Psi, rt] - \dot{r}_u^{-1} f_{r_u}^2 \frac{d}{dr_u} [\Psi, tt] + f_{r_u} f_{r_u}^{\prime} [\Psi, rt] - V_{r_u} [\Psi, t] \right\}, (25)$$
$$[\Psi_{,ttr}] = [\Psi_{,trt}] = [\Psi_{,rtt}] = \frac{d}{dr_u} [\Psi_{,tt}] - \dot{r}_u^{-1} [\Psi_{,ttt}] , \quad (26)$$

$$[\Psi_{,trr}] = [\Psi_{,rtr}] = [\Psi_{,rrt}] = \frac{d}{dr_u} [\Psi_{,tr}] - \dot{r}_u^{-1} [\Psi_{,ttr}] , \quad (27)$$

$$[\Psi_{,rrr}] = \frac{d}{dr_u} [\Psi_{,rr}] - \dot{r}_u^{-1} [\Psi_{,trr}] . \quad (28)$$

Finally, we similarly proceed for the fourth derivatives

$$[\Psi_{,tttt}] = \frac{\dot{r}_u^2}{\dot{r}_u^2 - f_{r_u}^2} \times \left\{ f_{r_u}^2 \frac{d}{dr_u} [\Psi_{,ttr}] - \dot{r}_u^{-1} f_{r_u}^2 \frac{d}{dr_u} [\Psi_{,ttt}] + f_{r_u} f'_{r_u} [\Psi_{,ttr}] - V_{r_u} [\Psi_{,tt}] \right\}, \quad (29)$$

$$[\Psi_{,tttr}] = [\Psi_{,ttrt}] = [\Psi_{,trtt}] = [\Psi_{,rttt}] \frac{d}{dr_u} [\Psi_{,ttt}] - \dot{r}_u^{-1} [\Psi_{,tttt}] , \quad (30)$$

$$[\Psi_{,ttrr}] = [\Psi_{,trtr}] = [\Psi_{,trrt}] = [\Psi_{,rttr}] [\Psi_{,rtrt}] = [\Psi_{,rrtt}] = \\ = \frac{d}{dr_u} [\Psi_{,ttr}] - \dot{r}_u^{-1} [\Psi_{,tttr}] , \quad (31)$$

$$[\Psi_{,trrr}] = [\Psi_{,rtrr}] = [\Psi_{,rrtr}] [\Psi_{,rrrt}] \frac{d}{dr_u} [\Psi_{,trr}] - \dot{r}_u^{-1} [\Psi_{,ttrr}] , \quad (32)$$

$$[\Psi_{,rrrr}] = \frac{d}{dr_u} [\Psi_{,rrr}] - \dot{r}_u^{-1} [\Psi_{,rrrt}] . \quad (33)$$

### 3.1.1. Jump conditions in explicit form.
We list hereafter the jump conditions in explicit form.

#### Jump conditions

$$[\Psi] = \frac{\kappa E r_u}{(\lambda + 1)(3M + \lambda r_u)} \quad (34)$$

#### First derivative jump conditions

$$[\Psi_{,t}] = - \frac{\kappa E r_u \dot{r}_u}{(2M - r_u)(3M + \lambda r_u)} \quad (35)$$

$$[\Psi_{,r}] = \frac{\kappa E [6M^2 + 3M\lambda r_u + \lambda(\lambda + 1)r_u^2]}{(\lambda + 1)(2M - r_u)(3M + \lambda r_u)^2} \quad (36)$$

#### Second derivative jump conditions

$$[\Psi_{,rr}] = - \frac{\kappa E [3M^3(5\lambda - 3) + 6M^2\lambda(\lambda - 3)r_u + 3M\lambda^2(\lambda - 1)r_u^2 - 2\lambda^2(\lambda + 1)r_u^3]}{(\lambda + 1)(2M - r_u)^2(3M + \lambda r_u)^3} \quad (37)$$

$$[\Psi_{,tr}] = \frac{\kappa E (3M^2 + 3M\lambda r_u - \lambda r_u^2) \dot{r}_u}{(2M - r_u)^2 (3M + \lambda r_u)^2} \quad (38)$$

$$[\Psi_{,tt}] = - \frac{\kappa E M}{r_u^2 (3M + r_u \lambda)} \quad (39)$$
Third derivative jump conditions

$$[\Psi_{rrr}] = \frac{\kappa E}{r_u (\lambda+1) (2M-r_u)^3 (3M+r_u \lambda)^4} \left[ 81 (\lambda+1) M^5 + 9r_u (19\lambda^2 + 18E^2\lambda + 3\lambda + 18E^2) M^4 + 9r_u^2 \lambda (7\lambda^2 + 24E^2\lambda - 14\lambda + 24E^2 + 3) M^3 + 3r_u^3 \lambda^2 (7\lambda^2 + 36E^2\lambda - 11\lambda + 36E^2 + 18) M^2 + 3r_u^4 \lambda^3 (8E^2\lambda - 7\lambda + 8E^2 - 1) M + 2r_u^5 \lambda^3 (\lambda+1) (E^2\lambda + 3) \right] \quad (40)$$

$$[\Psi_{trr}] = \frac{-\kappa E \dot{r}_u}{r_u (2M - r_u)^3 (3M + r_u \lambda)^3} \left[ 27M^4 + 6r_u (5\lambda + 9E^2 - 3) M^3 + 3r_u^2 \lambda (5\lambda + 18E^2 - 6) M^2 + 6r_u^3 \lambda^2 (3E^2 - 2) M + 2r_u^4 \lambda^2 (E^2\lambda + 1) \right] \quad (41)$$

$$[\Psi_{ttr}] = \frac{\kappa E}{r_u^3 (2M - r_u) (3M + r_u \lambda)^2} \left[ 39M^3 + 9r_u (3\lambda + 2E^2 - 2) M^2 + r_u^2 \lambda (4\lambda + 12E^2 - 13) M + 2r_u^3 \lambda^2 (E^2 - 1) \right] \quad (42)$$

$$[\Psi_{ttt}] = \frac{-\kappa E \dot{r}_u}{r_u^3 (2M - r_u) (3M + r_u \lambda)} \left[ 9M^2 + 2r_u (2\lambda + 3E^2 - 2) M + 2r_u^2 \lambda (E^2 - 1) \right] \quad (43)$$

Fourth derivative jump conditions

$$[\Psi_{rrrr}] = \frac{-3\kappa E}{r_u^2 (\lambda+1) (2M-r_u)^4 (3M+r_u\lambda)^5} \left[ 567 (\lambda+1) M^7 + 162 r_u (\lambda+1) (6\lambda + 16E^2 - 5) M^6 + 6r_u^2 (139\lambda^3 + 738E^2\lambda^2 - 123\lambda^2 + 162E^4\lambda + 441E^2\lambda - 171\lambda + 162E^4 - 297E^2 + 27) M^5 + 12r_u^3 \lambda (21\lambda^3 + 252E^2\lambda^2 - 85\lambda^2 + 135E^4\lambda - 24\lambda + 135E^4 - 252E^2 + 18) M^4 + 3r_u^4 \lambda^2 (21\lambda^3 + 344E^2\lambda^2 - 95\lambda^2 + 360E^4\lambda - 340E^2\lambda + 100\lambda + 360E^4 - 684E^2 + 24) M^3 + 2r_u^5 \lambda^3 (88E^2\lambda^2 - 47\lambda^2 + 180E^4\lambda - 260E^2\lambda + 25\lambda + 180E^4 - 348E^2 - 24) M^2 + 2r_u^6 \lambda^4 (6E^2\lambda^2 + 30E^4\lambda - 53E^2\lambda + 23\lambda + 30E^4 - 59E^2 + 11) M + 4r_u^7 \lambda^4 (\lambda+1) (E^4\lambda - 2E^2\lambda - 2) \right] \quad (44)$$

$$[\Psi_{trrr}] = \frac{3\kappa E \dot{r}_u}{r_u^2 (2M - r_u)^4 (3M + r_u \lambda)^4} \left[ 135M^6 + 27r_u (7\lambda + 32E^2 - 6) M^5 + 3r_u^2 (35\lambda^2 + 396E^2\lambda - 75\lambda + 108E^4 - 144E^2 + 18) M^4 + r_u^3 \lambda (35\lambda^2 + 612E^2\lambda - 120\lambda + 432E^4 - 594E^2 + 72) M^3 + r_u^4 \lambda^2 (140E^2\lambda - 45\lambda + 216E^4 - 306E^2 + 36) M^2 + 2r_u^5 \lambda^3 (6E^2\lambda + 24E^4 - 35E^2 + 9) M + 2r_u^6 \lambda^3 (2E^4\lambda - 3E^2\lambda - 1) \right] \quad (45)$$
$$[\Psi_{,ttrr}] = \frac{-\kappa E}{r_u^4 (2M - r_u)^2 (3M + r_u \lambda)^3} \left[ 1431 M^5 + 6r_u (251\lambda + 234E^2 - 210) M^4 + 9r_u^2 (59\lambda^2 + 160E^2\lambda - 148\lambda + 36E^4 - 66E^2 + 30) M^3 + 6r_u^3 \lambda (10\lambda^2 + 82E^2\lambda - 79\lambda + 54E^4 - 102E^2 + 48) M^2 + 2r_u^4 \lambda^2 (28E^2\lambda - 27\lambda + 54E^4 - 105E^2 + 52) M + 12r_u^5 \lambda^3 (E^2 - 1)^2 \right] \quad (46)$$

$$[\Psi_{,tttr}] = \frac{\kappa E \dot{r}_u}{r_u^4 (2M - r_u)^2 (3M + r_u \lambda)^2} \left[ 243 M^4 + 3r_u (61\lambda + 132E^2 - 64) M^3 + 3r_u^2 (12\lambda^2 + 92E^2\lambda - 49\lambda + 36E^4 - 48E^2 + 12) M^2 + 2r_u^3 \lambda (24E^2\lambda - 15\lambda + 36E^4 - 51E^2 + 14) M + 6r_u^4 \lambda^2 (E^2 - 1) (2E^2 - 1) \right] \quad (47)$$

$$[\Psi_{,tttt}] = \frac{-\kappa E}{r_u^6 (3M + r_u \lambda)} \left[ 189 M^3 + 2r_u (36\lambda + 84E^2 - 77) M^2 + 6r_u^2 (E^2 - 1) (10\lambda + 6E^2 - 5) M + 12r_u^3 \lambda (E^2 - 1)^2 \right] \quad (48)$$

While heuristic arguments [35,36] have been put forward to show that, for radial fall in the RW gauge, even metric perturbations belong to the $C^0$ continuity class at the position of the th particle, in [20,21] we have provided an analysis vis à vis the jump conditions that the wave function and its (first and second) derivatives have to satisfy for guaranteeing the continuity of the perturbations at the position of the particle. Therein, we have derived the same jump conditions (34 - 38) from the inverse relations (expressions giving the perturbations as function of the wave function and its derivatives) by fulfilment of the continuity conditions (equal coefficients for the two Heaviside distributions, and null coefficients for the Dirac distribution and its derivative).

## 4. The algorithm

The integration method considers cells belonging to two groups for cells never crossed by the world line, the integrating method may be drawn by previous approaches explored by Lousto [22] and Haas [15], whereas for cells crossed by a particle, we propose a new algorithm. The grid is in the $r^*, t$ domain.

Initial conditions require knowledge of the situation prior to $t = 0$. At fourth order, the wave function may be Taylor-expanded around $t = 0$. For the boundary conditions, simplicity suggests a sufficiently huge grid to avoid unwanted reflections.

### 4.1. Empty cells

Empty cells are those cells which are not crossed by the particle. In this case, the cell upper point is obtained by performing an integration of the wave equation over the entire surface A of the cell, identified by the nodes $\alpha, \beta, \gamma, \delta$. We briefly recall the algorithm used by Haas [15]. Therein, the sole numerical computation to be carried
out is represented by the product of the potential term and the wave function $V\Psi=g$. It is performed via a double Simpson integral, using points of the past light cone of the upper node $\alpha$, Fig. 1. We set $g_q = g(r_q^*, t_q) = V(r_q)\Psi(r_q^*, t_q)$, $V_q = V(r_q)$ and $\Psi_q = \Psi(r_q^*, t_q)$, where $q$ is one of the points shown in Fig. 1. The increment $h$ is defined as $h = \frac{1}{2}\Delta r^* = \frac{1}{2}\Delta t$ where $\Delta r^*$ is the spatial step and $\Delta t$ is the time step.

<figure><img src="image_2.png" /><figcaption>Figure 1: Set of points (circles and crosses) used for the integration of the $V\Psi=g$ term in the vacuum case. The crosses don't overlap with grid nodes; thus the field $g$ at these points, Eqs. (50,51), is approximated by the field at the nodes on the past light cone of the grid node $\alpha$.</figcaption></figure>

We have

$$ \iint_{Cell} g dA = \left(\frac{h}{3}\right)^2 \left[ g_{\alpha} + g_{\beta} + g_{\gamma} + g_{\delta} + 4(g_{\beta\gamma} + g_{\alpha\beta} + g_{\delta\gamma} + g_{\alpha\delta}) + 16g_{\sigma} \right] + \mathcal{O}(h^6), (49) $$

where the sum of the intermediate terms between nodes is given by

$$ \begin{aligned} g_{\beta\gamma} + g_{\alpha\beta} + g_{\delta\gamma} + g_{\alpha\delta} = & 2V_{\sigma}\Psi_{\sigma}\left[1 - \frac{1}{2}\left(\frac{h}{2}\right)^2 V_{\sigma}\right] + V_{\beta\gamma}\Psi_{\beta}\left[1 - \frac{1}{2}\left(\frac{h}{2}\right)^2 V_{\beta\gamma}\right] + \\ & V_{\delta\gamma}\Psi_{\delta}\left[1 - \frac{1}{2}\left(\frac{h}{2}\right)^2 V_{\delta\gamma}\right] + \frac{1}{2}\left[V_{\beta\gamma} - 2V_{\sigma} + V_{\delta\gamma}\right]\left(\Psi_{\beta} + \Psi_{\delta}\right) \\ & + \mathcal{O}(h^4). \end{aligned} \quad (50) $$

The last intermediate term $g_{\sigma}$ in Eq. 49 is evaluated using given nodes in the past light cone of $\alpha$, Fig. 1

$$ g_{\sigma} = \frac{1}{16} \left[ 8g_{\beta} + 8g_{\gamma} + 8g_{\delta} - 4g_{\gamma_1} - 4g_{\gamma_2} + g_{\mu_1} - g_{\mu_2} - g_{\mu_3} + g_{\mu_4} \right] + \mathcal{O}(h^4). \quad (51) $$

For the differential operators, an exact integration simply leads to

$$ \iint_{Cell} (\partial_{r^*}^2 - \partial_t^2) \Psi(r^*, t) dA = -4 \left[ \Psi_{\alpha} - \Psi_{\beta} + \Psi_{\gamma} - \Psi_{\delta} \right]. \quad (52) $$
Finally, we get

$$ \Psi_{\alpha} = -\Psi_{\gamma} + \Psi_{\beta} \left[ 1 - \frac{1}{4} \left(\frac{h}{2}\right)^2 (V_{\sigma} + V_{\beta}) + \frac{1}{16} \left(\frac{h}{2}\right)^4 V_{\sigma} (V_{\sigma} + V_{\beta}) \right] \\ + \Psi_{\delta} \left[ 1 - \frac{1}{4} \left(\frac{h}{2}\right)^2 (V_{\sigma} + V_{\delta}) + \frac{1}{16} \left(\frac{h}{2}\right)^4 V_{\sigma} (V_{\sigma} + V_{\delta}) \right] \\ - \left(\frac{h}{2}\right)^2 \left[ 1 - \frac{1}{4} \left(\frac{h}{2}\right)^2 V_{\sigma} \right] \left[ g_{\beta\gamma} + g_{\alpha\beta} + g_{\delta\gamma} + g_{\alpha\delta} + 4g_{\sigma} \right] . \quad (53) $$

For cells adjacent to cells crossed by the particle, the requirement of good accuracy suggests a different dealing for the computation of $g_{\sigma}$, since the past light cone of an adjacent cell can cross the path of the particle. In such a case, $g_{\sigma}$ is approximated by non-centred spatial finite difference expressions [15].

### 4.2. Cells crossed by the world line

For a given cell, our aim remains the determination of the wave function value at the upper node, now rebaptised $\alpha_0$. As in the previous section, we consider (fifteen) points both located in the past light cone of the $\alpha_0$ point and lying around a chosen point on the discontinuity $r_u(t)$, with the intent of determining $\Psi_{\alpha_0}$ by their linear combination. The non-regularity of the wave function due to the discontinuity, obviously entails a different value according to whether the discontinuity is approached from below ($\Psi^-$, left of the trajectory, Figs. 2 - 4) or above ($\Psi^+$, right of the trajectory, Figs. 2 - 4) the particle in radial fall. The same stands for the wave function derivatives. The addition of the jump condition to the value of the e.g. $\Psi^-$ ($\Psi^+$) wave function (or derivative of) allows to equate this sum to the value $\Psi^+$ ($\Psi^-$) of the wave function (or derivative of). This straightforward property turns being helpful for the achievement of the just mentioned linear combination of fifteen points. Incidentally, other linear combinations may be envisaged, though combinations of points located solely on one side of the discontinuity are to be avoided.

With reference to Figs. 2 - 4, there are three different cases depending upon how the trajectory of the particle crosses the cell wherein $\alpha_0$ lies. These three cases are further subdivided into three sub-cases, for a total of nine. In the following, we label by R the points on the right of the $[\alpha_0\alpha_6]$ line and by L the points on the left. Dealing with radial fall, and thereby with a 2D code, the up and down labels might be proper; nevertheless, we stick to right and left labels, given the orientation of the $r^*$ axis in the Figs. 2 - 4. For the first group of three, the trajectory crosses the $[\alpha_2\beta_1^R]$ and $[\alpha_0\beta_1^L]$ lines, Fig. 2; for the second group, the $[\alpha_2\beta_1^L]$ and $[\alpha_0\beta_1^L]$ lines, Fig. 3; finally for the third group, the $[\alpha_2\beta_1^R]$ and $[\alpha_0\beta_1^R]$ lines, Fig. 4.

We start considering the sub-case (1a) shown by Fig. 2, for which the trajectory crosses the line $[\alpha_0\alpha_2]$ at the point b. For compactness of the presentation of the final results, while we still adopt the same notation for the jump conditions, namely $[\Psi]_q$ for the difference $(\Psi^+ - \Psi^-)_{r_u=r_u(t_q)}$, for the jump derivatives instead, we rely henceforth on the notation $[\partial_{r^*}^n \partial_l^m \Psi]_q = (\partial_{r^*}^n \partial_l^m \Psi^+ - \partial_{r^*}^n \partial_l^m \Psi^-)_{r_u=r_u(t_q)}$, where $t_q$ is the coordinate time at the point $q = a, b$. We also define the lapse $\epsilon_b = t_{\alpha_0} - t_b$.

We recall that our aim is the determination of the value of $\Psi_{\alpha_0}^+$, knowing: i) $\epsilon_b$, ii) the jump (analytical) conditions on $\Psi$ and its derivatives at the point b; iii) the values of $\Psi$ on a set of fifteen points $\{\alpha, \beta, \gamma, \mu, \nu\}$ at the left and right sides of the world
line. A Taylor series is applied at each point around b up to fourth order, thereby obtaining

$$ \Psi_{\alpha_0}^{+} = \Psi^{+} (t_b + \epsilon_b, r_b^*) = \sum_{n=0}^{4} \frac{\epsilon_b^n}{n!} \partial_t^n \Psi_b^{+} + \mathcal{O}(\epsilon_b^5), \quad (54) $$

$$ \Psi_{\alpha_i}^{-} = \Psi^{-} (t_b - (ih - \epsilon_b), r_b^*) = \sum_{n=0}^{4} (-1)^n \frac{(ih - \epsilon_b)^n}{n!} \partial_t^n \Psi_b^{-} + \mathcal{O}(h^5), \quad (55) $$

$$ \Psi_{\beta_j^{R,L}}^{\pm} = \Psi^{\pm} (t_b - (jh - \epsilon_b), r_b^{\pm} \pm h) = \sum_{n+m \le 4} (-1)^m (\pm 1)^n \frac{h^n}{n!} \frac{(jh - \epsilon_b)^m}{m!} \partial_{r^*}^n \partial_t^m \Psi_b^{\pm} + \mathcal{O}(h^5), \quad (56) $$

$$ \Psi_{\gamma_k^{R,L}}^{\pm} = \Psi^{\pm} (t_b - (kh - \epsilon_b), r_b^{\pm} \pm 2h) = \sum_{n+m \le 4} (-1)^m (\pm 1)^n \frac{(2h)^n}{n!} \frac{(kh - \epsilon_b)^m}{m!} \partial_{r^*}^n \partial_t^m \Psi_b^{\pm} + \mathcal{O}(h^5), \quad (57) $$

$$ \Psi_{\mu_3^{R,L}}^{\pm} = \Psi^{\pm} (t_b - (3h - \epsilon_b), r_b^{\pm} \pm 3h) = \sum_{n+m \le 4} (-1)^m (\pm 1)^n \frac{(3h)^n}{n!} \frac{(3h - \epsilon_b)^m}{m!} \partial_{r^*}^n \partial_t^m \Psi_b^{\pm} + \mathcal{O}(h^5), \quad (58) $$

$$ \Psi_{\nu_4^{R,L}}^{\pm} = \Psi^{\pm} (t_b - (4h - \epsilon_b), r_b^{\pm} \pm 4h) = \sum_{n+m \le 4} (-1)^m (\pm 1)^n \frac{(4h)^n}{n!} \frac{(4h - \epsilon_b)^m}{m!} \partial_{r^*}^n \partial_t^m \Psi_b^{\pm} + \mathcal{O}(h^5), \quad (59) $$

for the indexes running as $i = 2, 4, 6, j = 1, 3$ and $k = 2, 4$ and concerning the $\alpha, \beta$ and $\gamma$ nodes, respectively. Our notation implies that the subscript R, L stands for R when the superscript $\pm$ corresponds to $+$, whereas R, L stands for L when $\pm$ corresponds to $-$.

$$ \begin{aligned} \Psi_{\alpha_0}^{+} &= \sum_{n=0}^{4} c_n \partial_t^n \Psi_b^{+} + \mathcal{O}(h^5) = \sum_{n=0}^{4} c_n (\partial_t^n \Psi_b^{-} + [\partial_t^n \Psi]_b) + \mathcal{O}(h^5) = \\ &= c_0 \Psi_b^{-} + c_1 \partial_t \Psi_b^{-} + c_2 \partial_t^2 \Psi_b^{-} + c_3 \partial_t^3 \Psi_b^{-} + c_4 \partial_t^4 \Psi_b^{-} + \sum_{n=0}^{4} c_n [\partial_t^n \Psi]_b + \mathcal{O}(h^5) \end{aligned} \quad (60) $$

For an accuracy at fourth order, all quantities $\mathcal{O}(h^5)$ are disregarded. The sum $\hat{S} = c_0 \Psi_b^{-} + c_1 \partial_t \Psi_b^{-} + c_2 \partial_t^2 \Psi_b^{-} + c_3 \partial_t^3 \Psi_b^{-} + c_4 \partial_t^4 \Psi_b^{-}$, Eq. 60, is composed by numerical derivatives of lower order than $\mathcal{O}(h^5)$, and therefore they can't be neglected. However, the computation of high order derivatives is often accompanied by numerical noise. Therefore, we replace this sum by a combination of wave function values in the $\alpha_0$ light cone. This is attained in two steps. The former involves taking fifteen wave function values on the two sides of the trajectory, that is $\{\Psi_{\alpha_i}^{-}, \Psi_{\beta_j}^{-}, \Psi_{\beta_j^{R}}^{+}, \Psi_{\gamma_k}^{-}, \Psi_{\gamma_k^{R}}^{+}, \Psi_{\mu_3^{L}}^{-}, \Psi_{\mu_3^{R}}^{+}, \Psi_{\nu_4^{-}}^{+}, \Psi_{\nu_4^{R}}^{+}\}$, Fig. 2. The latter employs the jump conditions to relate the fifteen mentioned points with $\{\Psi_{\alpha_i}^{-}, \Psi_{\beta_j}^{-}, \Psi_{\beta_j^{R}}^{-}, \Psi_{\gamma_k}^{-}, \Psi_{\gamma_k^{R}}^{-}, \Psi_{\mu_3^{-}}^{+}, \Psi_{\mu_3^{R}}^{-}, \Psi_{\nu_4^{-}}^{+}, \Psi_{\nu_4^{R}}^{-}\}$. For the former step, we define the sum $S$
$$S = \sum_{i} (A_{i} \Psi_{\alpha_{i}}^{-}) + \sum_{j} (B_{j}^{L} \Psi_{\beta_{j}^{L}}^{-} + B_{j}^{R} \Psi_{\beta_{j}^{R}}^{+}) + \sum_{k} (G_{k}^{L} \Psi_{\gamma_{k}^{L}}^{-} + G_{k}^{R} \Psi_{\gamma_{k}^{R}}^{+}) \\ + M_{3}^{L} \Psi_{\mu_{3}^{L}}^{-} + M_{3}^{R} \Psi_{\mu_{3}^{R}}^{+} + N_{4}^{L} \Psi_{\nu_{4}^{L}}^{-} + N_{4}^{R} \Psi_{\nu_{4}^{R}}^{+}. \qquad (61)$$

where {$A_i, B_j^L, B_j^R, G_k^L, G_k^R, M_3^L, M_3^R, N_4^L, N_4^R$} are constants.

We observe that the $\hat{S}$ sum entails only wave function values at the left of the b point on the trajectory. The jump conditions are once more exploited to relate the two domains $r^* < r_u^*(t)$ and $r^* > r_u^*(t)$. This specifically concerns six points {$\beta_j^R, \gamma_k^R, \mu_3^R, \nu_4^R$}. For instance, at the $\beta_j^R$ point, we can write

<figure><img src="image_3.png" /><figcaption>Figure 2: The three sub-cases for which the particle enters through the [$\alpha_2\beta_1^R$] side and leaves through the [$\alpha_0\beta_1^L$] side. The elimination of the $\Psi_b^-$ derivatives demands, Eq. 60, the utilisation of fifteen points, represented by circles, in the light cone of $\alpha_0$. Numerical efficiency suggests that the points are taken at both left and right sides of the $r_u^*(t)$ trajectory. In the three cases, the particle crosses the line [$\alpha_0\alpha_2$] at the point b. The background distinguishes two zones: one where $\Psi(r^* < r_u^*(t), t) = \Psi^-(r^*, t)$, the other where $\Psi(r^* > r_u^*(t), t) = \Psi^+(r^*, t)$, the path $r_u^*(t)$ representing the separation between the two zones.</figcaption></figure>
$$ \begin{aligned} \Psi_{\beta_j}^{+} &= \sum_{n+m \le 4} (-1)^m \frac{h^n}{n!} \frac{(jh-\epsilon_b)^m}{m!} (\partial_{r_*}^n \partial_t^m \Psi_b^- + [\partial_{r_*}^n \partial_t^m \Psi]_b) + \mathcal{O}(h^5) \\ &= \Psi_{\beta_j}^{-} + \sum_{n+m \le 4} (-1)^m \frac{h^n}{n!} \frac{(jh-\epsilon_b)^m}{m!} [\partial_{r_*}^n \partial_t^m \Psi]_b, \end{aligned} \quad (62) $$

where

$$ \Psi_{\beta_j}^{-} = \sum_{n+m \le 4} (-1)^m \frac{h^n}{n!} \frac{(jh-\epsilon_b)^m}{m!} (\partial_{r_*}^n \partial_t^m \Psi_b^-) + \mathcal{O}(h^5) . \quad (63) $$

<figure><img src="image_4.png" /><figcaption>Figure 3: The three sub-cases for which the particle enters through the [$\alpha_2\beta_1^L$] side and leaves through the [$\alpha_0\beta_1^L$] side. The elimination of the $\Psi_a^-$ derivatives demands the utilisation of fifteen points, represented by circles, in the light cone of $\alpha_0$. Numerical efficiency suggests that the points are taken at both left and right sides of the $r_u^*(t)$ trajectory. In the three cases, the particle crosses the line [$\beta_1^L\beta_1^R$] at the point a. The background distinguishes two zones: one where $\Psi(r^* < r_u^*(t), t) = \Psi^-(r^*, t)$, the other where $\Psi(r^* > r_u^*(t), t) = \Psi^+(r^*, t)$, the path $r_u^*(t)$ representing the separation between the two zones.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: The three sub-cases for which the particle enters through the $[\alpha_2\beta_1^R]$ side and leaves through the $[\alpha_0\beta_1^R]$ side. The elimination of the $\Psi_a^-$ derivatives demands the utilisation of fifteen points, represented by circles, in the light cone of $\alpha_0$. Numerical efficiency suggests that the points are taken at both left and right sides of the $r_u^*(t)$ trajectory. In the three cases, the particle crosses the line $[\beta_1^L\beta_1^R]$ at the point a. The background distinguishes two zones: one where $\Psi(r^* < r_u^*(t), t) = \Psi^-(r^*, t)$, the other where $\Psi(r^* > r_u^*(t), t) = \Psi^+(r^*, t)$, the path $r_u^*(t)$ representing the separation between the two zones.</figcaption></figure>

By application of the same transformation to the quantities $\Psi_{\gamma_k^R}^+, \Psi_{\mu_3^R}^+, \Psi_{\nu_3^R}^+$, Eq. 61 becomes

$$S - \Phi_{r_u^*}^{jump} = \sum_i (A_i \Psi_{\alpha_i}^-) + \sum_j (B_j^L \Psi_{\beta_j^L}^- + B_j^R \Psi_{\beta_j^R}^-) + \sum_k (G_k^L \Psi_{\gamma_k^L}^- + G_k^R \Psi_{\gamma_k^R}^-) \\ + M_3^L \Psi_{\mu_3^L}^- + M_3^R \Psi_{\mu_3^R}^- + N_4^L \Psi_{\nu_4^L}^- + N_4^R \Psi_{\nu_4^R}^-, \qquad (64)$$

where $\Phi_{r_u^*}^{jump}$ is an analytic function, composed by the jump conditions at the b point, weighted by coefficients issued by Eq. 62 or similar equations.

Having only $\Psi^-$ terms on the right hand side of Eq. 64, we can finally search the coefficients $\{A_i, B_j^L, B_j^R, G_k^L, G_k^R, M_3^L, M_3^R, N_4^L, N_4^R\}$ that satisfy the equation
$\hat{S} = S - \Phi_{r_u^*}^{jump}$, that is

$$c_0 \Psi_b^- + c_1 \partial_t \Psi_b^- + c_2 \partial_t^2 \Psi_b^- + c_3 \partial_t^3 \Psi_b^- + c_4 \partial_t^4 \Psi_b^- = \\ \sum_i (A_i \Psi_{\alpha_i}^-) + \sum_j (B_j^L \Psi_{\beta_j^L}^- + B_j^R \Psi_{\beta_j^R}^-) + \sum_k (G_k^L \Psi_{\gamma_k^L}^- + G_k^R \Psi_{\gamma_k^R}^-) \\ + M_3^L \Psi_{\mu_3^L}^- + M_3^R \Psi_{\mu_3^R}^- + N_4^L \Psi_{\nu_4^L}^- + N_4^R \Psi_{\nu_4^R}^- . \qquad (65)$$

Using the notation of Eqs. 62, 63, and by injection of Eqs. 55-59, a Taylor expansion of fourth order at the b point is applied to the right-hand side of Eq. 65. The system can be cast in a matrix form

$$T \cdot P = C, \qquad (66)$$

where $P$ is the unknown 15-vector formed by the coefficients $\{A_i, B_j^L, B_j^R, G_k^L, G_k^R, M_3^L, M_3^R, N_4^L, N_4^R\}$

$$P = (A_2, A_4, A_6, B_1^L, B_3^L, B_1^R, B_3^R, G_2^L, G_4^L, G_2^R, G_4^R, M_3^L, M_3^R, N_4^L, N_4^R)^t, \quad (67)$$

and $C$ is given by the 15-vector

$$C = (c_0, c_1, c_2, c_3, c_4, 0, \dots, 0)^t, \qquad (68)$$

while $T$ is the $(15 \times 15)$ matrix constructed from the Taylor coefficients in Eqs. 55-59 (see appendix). By inversion of $T$, we get $P$ and specifically

$$\begin{aligned} A_2 &= \frac{-27}{5}, & A_4 &= \frac{-9}{5}, & A_6 &= \frac{1}{5}, \\ B_1^L &= B_1^R = \frac{12}{5}, & B_3^L &= B_3^R = \frac{18}{5}, \\ G_2^L &= G_2^R = \frac{-9}{5}, & G_4^L &= G_4^R = \frac{-3}{5}, \\ M_3^L &= M_3^R = \frac{2}{5}, & N_4^L &= N_4^R = 0. \end{aligned}$$

The following equivalences path the last stretch of the way

$$\Psi_{\alpha_0}^{+} = S - \Phi_{r_u^*}^{jump} + \sum_{n=0}^{4} c_n [\partial_t^n \Psi]_b = S + \Phi_{r_u^*}^{(1)}(t_b), \qquad (69)$$

and explicitly, we get

$$\Psi_{\alpha_0}^{+} = -\frac{27}{5}\Psi_{\alpha_2}^{-} - \frac{9}{5}\Psi_{\alpha_4}^{-} + \frac{1}{5}\Psi_{\alpha_6}^{-} + \frac{12}{5}\left(\Psi_{\beta_1^{-}}^{L} + \Psi_{\beta_1^{+}}^{R}\right) + \frac{18}{5}\left(\Psi_{\beta_3^{-}}^{L} + \Psi_{\beta_3^{+}}^{R}\right) \\ -\frac{9}{5}\left(\Psi_{\gamma_2^{-}}^{L} + \Psi_{\gamma_2^{+}}^{R}\right) + \frac{3}{5}\left(\Psi_{\gamma_4^{-}}^{L} + \Psi_{\gamma_4^{+}}^{R}\right) - \frac{2}{5}\left(\Psi_{\mu_3^{-}}^{L} + \Psi_{\mu_3^{+}}^{R}\right) + \Phi_{r_u^*}^{(1)}(t_b), \qquad (70)$$

where $\Psi_{\beta_3^R}^{\pm} = \Psi_{\beta_3^L}^{\pm}$ for sub-case (1a), and $\Psi_{\beta_3^R}^{\pm} = \Psi_{\beta_3^L}^{-}$ for sub-cases (1b,1c); $\Psi_{\gamma_4^R}^{\pm} = \Psi_{\gamma_4^L}^{\pm}$ for sub-cases (1a,1b), and $\Psi_{\gamma_4^R}^{\pm} = \Psi_{\gamma_4^L}^{-}$ for sub-case (1c); and $\Phi_{r_u^*}(t_b)$ is an analytic function, that for the (1a) sub-case, it takes the value
$$
\begin{equation}
\begin{aligned}
\Phi_{r_u^*(t_b)}^{(1a)} = & -3[\Psi]_b - \frac{3(5\epsilon_b - 14h)}{5}[\partial_t\Psi]_b - \frac{3(\epsilon_b - 2h)(5\epsilon_b - 18h)}{10}[\partial_t^2\Psi]_b \\
& - \frac{5\epsilon_b^3 - 42h\epsilon_b^2 + 108h^2\epsilon_b - 96h^3}{10}[\partial_t^3\Psi]_b \\
& - \frac{5\epsilon_b^4 - 56h\epsilon_b^3 + 216h^2\epsilon_b^2 - 384h^3\epsilon_b + 240h^4}{40}[\partial_t^4\Psi]_b - \frac{12h}{5}[\partial_{r^*}\Psi]_b \\
& + \frac{2h^3}{5}[\partial_{r^*}^3\Psi]_b - \frac{12h(\epsilon_b - 2h)}{5}[\partial_{r^*}\partial_t\Psi]_b - \frac{6h(\epsilon_b^2 - 4h\epsilon_b + 5h^2)}{5}[\partial_{r^*}^3\partial_t\Psi]_b \\
& + \frac{2h^3(\epsilon_b - h)}{5}[\partial_{r^*}^2\partial_t^2\Psi]_b - \frac{2h(\epsilon_b^3 - 6h\epsilon_b^2 + 15h^2\epsilon_b - 11h^3)}{5}[\partial_{r^*}\partial_t^3\Psi]_b.
\end{aligned}
\tag{71}
\end{equation}
$$

The quantity $\Phi_{r_u^*}^{(1)}$ varies according to the different sub-cases: for the case (1b) of Fig. 2, the point $\beta_3^R$, whereas for the case (1c) the points $\beta_3^R$ and $\gamma_4^R$ are in the $r^* < r_u^*$ domain. Therefore

$$
\begin{equation}
\begin{aligned}
\Phi_{r_u^*(t_b)}^{(1b)} = & \frac{6}{10}[\Psi]_b + \frac{3(\epsilon_b - 4h)}{5}[\partial_t\Psi]_b + \frac{3(\epsilon_b^2 - 8h\epsilon_b + 18h^2)}{10}[\partial_t^2\Psi]_b \\
& + \frac{\epsilon_b^3 - 12h\epsilon_b^2 + 54h^2\epsilon_b - 66h^3}{10}[\partial_t^3\Psi]_b \\
& + \frac{\epsilon_b^4 - 16h\epsilon_b^3 + 108h^2\epsilon_b^2 - 264h^3\epsilon_b + 246h^4}{40}[\partial_t^4\Psi]_b \\
& + \frac{6h}{5}[\partial_{r^*}\Psi]_b + \frac{9h^2}{5}[\partial_{r^*}^2\Psi]_b + h^3[\partial_{r^*}^3\Psi]_b + \frac{3h^4}{20}[\partial_{r^*}^4\Psi]_b \\
& + \frac{6h(\epsilon_b - 5h)}{5}[\partial_{r^*}\partial_t\Psi]_b + \frac{9h^2(\epsilon_b - 3h)}{5}[\partial_{r^*}^2\partial_t\Psi]_b \\
& + \frac{3h(\epsilon_b^2 - 10h\epsilon_b + 17h^2)}{5}[\partial_{r^*}^3\partial_t\Psi]_b + \frac{9h^2(\epsilon_b - 3h)^2}{10}[\partial_{r^*}^2\partial_t^2\Psi]_b \\
& + \frac{h^3(5\epsilon_b - 11h)}{5}[\partial_{r^*}^2\partial_t^2\Psi]_b + \frac{h(\epsilon_b^3 - 15h\epsilon_b^2 + 51h^2\epsilon_b - 59h^3)}{5}[\partial_{r^*}\partial_t^3\Psi]_b,
\end{aligned}
\tag{72}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
\Phi_{r_u^*(t_b)}^{(1c)} = & \frac{3h^2}{5}[\partial_t^2\Psi]_b + \frac{h^2(3\epsilon_b - h)}{5}[\partial_t^3\Psi]_b + \frac{h^2(6\epsilon_b^2 - 4h\epsilon_b - 5h^2)}{20}[\partial_t^4\Psi]_b \\
& + \frac{3h^2}{5}[\partial_{r^*}^2\Psi]_b + \frac{h^3}{5}[\partial_{r^*}^3\Psi]_b - \frac{h^4}{4}[\partial_{r^*}^4\Psi]_b - \frac{6h^2}{5}[\partial_{r^*}\partial_t\Psi]_b \\
& + \frac{3h^2(\epsilon_b - h)}{5}[\partial_{r^*}^2\partial_t\Psi]_b - \frac{3h^2(2\epsilon_b - h)}{5}[\partial_{r^*}^3\partial_t\Psi]_b \\
& + \frac{3h^2(\epsilon_b^2 - 2h\epsilon_b - 5h^2)}{10}[\partial_{r^*}\partial_t^2\Psi]_b + \frac{h^3(\epsilon_b + 5h)}{5}[\partial_{r^*}^2\partial_t^2\Psi]_b \\
& - \frac{h^2(3\epsilon_b^2 - 3h\epsilon_b - 5h^2)}{5}[\partial_{r^*}\partial_t^3\Psi]_b.
\end{aligned}
\tag{73}
\end{equation}
$$

We thus have obtained, without direct integration of the singular source and the potential term, the value of the upper node. The equations shows three types of
terms: the preceding node values of the same cell, the jump conditions which are fully analytical quantities, and the wave function values at adjacent cells. Incidentally, at first order [21], the latter type of terms disappears and a simpler expression is obtained.

Similar relations are found for the other two remaining cases. For case 2, Fig. 3, we obtain (having defined the shift $\epsilon_a = t_{\beta_1^R} - r_a^*$)

$$ \Psi_{\alpha_0}^{+} = -\frac{27}{5}\Psi_{\alpha_2}^{+} - \frac{9}{5}\Psi_{\alpha_4}^{\pm} + \frac{1}{5}\Psi_{\alpha_6}^{\pm} + \frac{12}{5}(\Psi_{\beta_1^L}^{-} + \Psi_{\beta_1^R}^{+}) + \frac{18}{5}(\Psi_{\beta_3^L}^{-} + \Psi_{\beta_3^R}^{+}) \\ - \frac{9}{5}(\Psi_{\gamma_2^L}^{-} + \Psi_{\gamma_2^R}^{+}) + \frac{3}{5}(\Psi_{\gamma_4^L}^{-} + \Psi_{\gamma_4^R}^{+}) - \frac{2}{5}(\Psi_{\mu_3^L}^{-} + \Psi_{\mu_3^R}^{+}) + \Phi_{r_u^*(t_a)}^{(2)} \quad (74) $$

where $\Psi_{\alpha_4}^{\pm} = \Psi_{\alpha_4}^{+}$ for sub-cases (2a,2b), and $\Psi_{\alpha_4}^{\pm} = \Psi_{\alpha_4}^{-}$ for sub-case (2c); $\Psi_{\alpha_6}^{\pm} = \Psi_{\alpha_6}^{+}$ for sub-case (2a), and $\Psi_{\alpha_6}^{\pm} = \Psi_{\alpha_6}^{-}$ for sub-cases (2b,2c). For the (2a) sub-case, $\Phi_{r_u^*(t_a)}^{(2)}$ takes the following value

$$ \Phi_{r_u^*(t_a)}^{(2a)} = 4[\Psi]_a - \frac{22h}{5}[\partial_t\Psi]_a + \frac{22h^2}{5}[\partial_t^2\Psi]_a - \frac{7h^3}{3}[\partial_t^3\Psi]_a + \frac{17h^4}{30}[\partial_t^4\Psi]_a \\ + \frac{4(5\epsilon_a - 8h)}{5}[\partial_{r^*}\Psi]_a + \frac{2(\epsilon_a - h)(5\epsilon_a - 11h)}{5}[\partial_{r^*}^2\Psi]_a \\ + \frac{2(5\epsilon_a^3 - 24h\epsilon_a^2 + 33h^2\epsilon_a - 11h^3)}{15}[\partial_{r^*}^3\Psi]_a \\ + \frac{(\epsilon_a - h)(5\epsilon_a^3 - 27h\epsilon_a^2 + 39h^2\epsilon_a - 5h^3)}{30}[\partial_{r^*}^4\Psi]_a \\ - \frac{2h(11\epsilon_a - 17h)}{5}[\partial_{r^*}\partial_t\Psi]_a - \frac{h(\epsilon_a - h)(11\epsilon_a - 23h)}{5}[\partial_{r^*}^2\partial_t\Psi]_a \\ + \frac{2h^2(11\epsilon_a - 17h)}{5}[\partial_{r^*}\partial_t^2\Psi]_a - \frac{h(\epsilon_a - h)^2(11\epsilon_a - 29h)}{15}[\partial_{r^*}^3\partial_t\Psi]_a \\ + \frac{h^3(\epsilon_b + 5h)}{5}[\partial_{r^*}^2\partial_t^2\Psi]_b - \frac{h^3(35\epsilon_a - 41h)}{15}[\partial_{r^*}\partial_t^3\Psi]_a \quad (75) $$

For the same preceding reason, the sub-cases (2b, 2c) differ as the points $\alpha_4$ and $\alpha_6$ are or aren't in the $r^* > r_u^*$ domain. Therefore, we have

$$ \Phi_{r_u^*(t_a)}^{(2b)} = \frac{42}{10}[\Psi]_a - \frac{27h}{5}[\partial_t\Psi]_a + \frac{69h^2}{10}[\partial_t^2\Psi]_a - \frac{13h^3}{2}[\partial_t^3\Psi]_a + \frac{231h^4}{40}[\partial_t^4\Psi]_a \\ + \frac{3(7\epsilon_a - 11h)}{5}[\partial_{r^*}\Psi]_a + \frac{3(\epsilon_a - h)(7\epsilon_a - 15h)}{10}[\partial_{r^*}^2\Psi]_a \\ + \frac{7\epsilon_a^3 - 33h\epsilon_a^2 + 45h^2\epsilon_a - 15h^3}{10}[\partial_{r^*}^3\Psi]_a \\ + \frac{(\epsilon_a - h)(7\epsilon_a^3 - 37h\epsilon_a^2 + 53h^2\epsilon_a - 7h^3)}{40}[\partial_{r^*}^4\Psi]_a \\ - \frac{3h(9\epsilon_a - 13h)}{5}[\partial_{r^*}\partial_t\Psi]_a - \frac{3h(\epsilon_a - h)(9\epsilon_a - 17h)}{10}[\partial_{r^*}^2\partial_t\Psi]_a \\ + \frac{3h^2(23\epsilon_a - 31h)}{10}[\partial_{r^*}\partial_t^2\Psi]_a - \frac{3h(\epsilon_a - h)^2(3\epsilon_a - 7h)}{10}[\partial_{r^*}^3\partial_t\Psi]_a \\ + \frac{h^2(\epsilon_a - h)(157\epsilon_a - 109h)}{20}[\partial_{r^*}^2\partial_t^2\Psi]_a - \frac{h^3(65\epsilon_a - 69h)}{10}[\partial_{r^*}\partial_t^3\Psi]_a, $$
(76)

$$
\begin{equation}
\begin{aligned}
\Phi_{r_u^*(t_a)}^{(2c)} = & \frac{24}{10} [\Psi]_a - \frac{6h^2}{5} [\partial_t^2 \Psi]_a + \frac{8h^3}{5} [\partial_t^3 \Psi]_a - \frac{3h^4}{10} [\partial_t^4 \Psi]_a \\
& + \frac{12(\epsilon_a - 2h)}{5} [\partial_{r^*} \Psi]_a + \frac{6(\epsilon_a - 3h)(\epsilon_a - h)}{5} [\partial_{r^*}^2 \Psi]_a \\
& + \frac{2(\epsilon_a^3 - 6h\epsilon_a^2 + 9h^2\epsilon_a - 3h^3)}{5} [\partial_{r^*}^3 \Psi]_a \\
& + \frac{(\epsilon_a - h)(\epsilon_a^3 - 7h\epsilon_a^2 + 11h^2\epsilon_a - h^3)}{10} [\partial_{r^*}^4 \Psi]_a + \frac{12h^2}{5} [\partial_{r^*} \partial_t \Psi]_a \\
& + \frac{12h^2(\epsilon_a - h)}{5} [\partial_{r^*}^2 \partial_t \Psi]_a - \frac{6h^2(\epsilon_a + h)}{5} [\partial_{r^*} \partial_t^2 \Psi]_a \\
& + \frac{6h^2(\epsilon_a - h)^2}{5} [\partial_{r^*}^3 \partial_t \Psi]_a + \frac{h^2(\epsilon_a - h)(19\epsilon_a - 7h)}{5} [\partial_{r^*} \partial_t^2 \Psi]_a \\
& + \frac{2h^3(4\epsilon_a - 3h)}{5} [\partial_{r^*} \partial_t^3 \Psi]_a.
\end{aligned}
\tag{77}
\end{equation}
$$

Finally for case 3, Fig. 4, we have

$$
\begin{equation}
\begin{aligned}
\Psi_{\alpha_0}^{-} = & -\frac{27}{5}\Psi_{\alpha_2}^{-} - \frac{9}{5}\Psi_{\alpha_4}^{-} + \frac{1}{5}\Psi_{\alpha_6}^{-} + \frac{12}{5}(\Psi_{\beta_1^{-}}^{L} + \Psi_{\beta_1^{+}}^{R}) + \frac{18}{5}(\Psi_{\beta_3^{-}}^{L} + \Psi_{\beta_3^{+}}^{R}) \\
& -\frac{9}{5}(\Psi_{\gamma_2^{-}}^{L} + \Psi_{\gamma_2^{+}}^{R}) + \frac{3}{5}(\Psi_{\gamma_4^{-}}^{L} + \Psi_{\gamma_4^{+}}^{R}) - \frac{2}{5}(\Psi_{\mu_3^{-}}^{L} + \Psi_{\mu_3^{+}}^{R}) + \Phi_{r_u(t_a)}^{(3)}
\end{aligned}
\tag{78}
\end{equation}
$$

where $\Psi_{\beta_3^R}^{\pm} = \Psi_{\beta_3^L}^{+}$ for sub-case (3a), and $\Psi_{\beta_3^R}^{\pm} = \Psi_{\beta_3^L}^{-}$ for sub-cases (3b,3c); $\Psi_{\gamma_4^R}^{\pm} = \Psi_{\gamma_4^L}^{+}$ for sub-cases (3a,3b), and $\Psi_{\gamma_4^R}^{\pm} = \Psi_{\gamma_4^L}^{-}$ for sub-case (3c); and $\Phi_{r_u(t_a)}^{(3)}$ takes the values

$$
\begin{equation}
\begin{aligned}
\Phi_{r_u^*(t_a)}^{(3a)} = & -4[\Psi]_a + \frac{22h}{5}[\partial_t\Psi]_a - \frac{22h^2}{5}[\partial_t^2\Psi]_a + \frac{7h^3}{3}[\partial_t^3\Psi]_a - \frac{17h^4}{30}[\partial_t^4\Psi]_a \\
& -\frac{4(5\epsilon_a - 2h)}{5}[\partial_{r^*}\Psi]_a - \frac{2(\epsilon_a - h)(5\epsilon_a + h)}{5}[\partial_{r^*}^2\Psi]_a \\
& -\frac{2(5\epsilon_a^3 - 6h\epsilon_a^2 - 3h^2\epsilon_a + h^3)}{15}[\partial_{r^*}^3\Psi]_a \\
& -\frac{(\epsilon_a - h)(5\epsilon_a^3 - 3h\epsilon_a^2 - 9h^2\epsilon_a - 5h^3)}{30}[\partial_{r^*}^4\Psi]_a \\
& +\frac{2h(11\epsilon_a - 5h)}{5}[\partial_{r^*}\partial_t\Psi]_a + \frac{h(\epsilon_a - h)(11\epsilon_a + h)}{5}[\partial_{r^*}^2\partial_t\Psi]_a \\
& -\frac{2h^2(11\epsilon_a - 5h)}{5}[\partial_{r^*}^3\partial_t^2\Psi]_a + \frac{h(\epsilon_a - h)^2(11\epsilon_a + 7h)}{15}[\partial_{r^*}^4\partial_t\Psi]_a \\
& +\frac{h^2(\epsilon_a - h)(11\epsilon_a + h)}{5}[\partial_{r^*}^2\partial_t^2\Psi]_a + \frac{h^3(35\epsilon_a - 29h)}{15}[\partial_{r^*}^3\partial_t^3\Psi]_a,
\end{aligned}
\tag{79}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
\Phi_{r_u^*(t_a)}^{(3b)} = & -\frac{2}{5}[\Psi]_a - \frac{14h}{5}[\partial_t\Psi]_a + \frac{14h^2}{5}[\partial_t^2\Psi]_a - \frac{37h^3}{15}[\partial_t^3\Psi]_a + \frac{11h^4}{6}[\partial_t^4\Psi]_a \\
& -\frac{2(\epsilon_a - 4h)}{5}[\partial_{r^*}\Psi]_a - \frac{\epsilon_a^2 - 8h\epsilon_a - 2h^2}{5}[\partial_{r^*}^2\Psi]_a
\end{aligned}
\tag{79}
\end{equation}
$$
$$
\begin{aligned}
& - \frac{\epsilon_a^3 - 12 h \epsilon_a^2 - 6 h^2 \epsilon_a + 2 h^3}{15} [\partial_{r^*}^3 \Psi]_a \\
& - \frac{\epsilon_a^4 - 16 h \epsilon_a^3 - 12 h^2 \epsilon_a^2 + 8 h^3 \epsilon_a + 10 h^4}{60} [\partial_{r^*}^4 \Psi]_a \\
& - \frac{2 h (7 \epsilon_a + 5 h)}{5} [\partial_{r^*} \partial_t \Psi]_a - \frac{h (7 \epsilon_a^2 + 10 h \epsilon_a + h^2)}{5} [\partial_{r^*}^2 \partial_t \Psi]_a \\
& + \frac{2 h^2 (7 \epsilon_a + 5 h)}{5} [\partial_{r^*} \partial_t^2 \Psi]_a - \frac{h (7 \epsilon_a^3 + 15 h \epsilon_a^2 + 3 h^2 \epsilon_a - 7 h^3)}{15} [\partial_{r^*}^3 \partial_t \Psi]_a \\
& - \frac{h^2 (7 \epsilon_a^2 + 10 h \epsilon_a + h^2)}{5} [\partial_{r^*}^2 \partial_t^2 \Psi]_a - \frac{h^3 (37 \epsilon_a + 29 h)}{15} [\partial_{r^*} \partial_t^3 \Psi]_a,
\end{aligned}
\quad (80)
$$

$$
\begin{aligned}
\Phi_{r_u^*(t_a)}^{(3c)} = & -[\Psi]_a - h[\partial_t \Psi]_a + \frac{h^2}{10}[\partial_t^2 \Psi]_a + \frac{7h^3}{30}[\partial_t^3 \Psi]_a - \frac{23h^4}{120}[\partial_t^4 \Psi]_a \\
& + (h - \epsilon_a)[\partial_{r^*} \Psi]_a - \frac{5\epsilon_a^2 - 10h\epsilon_a - h^2}{10}[\partial_{r^*}^2 \Psi]_a \\
& - \frac{5\epsilon_a^3 - 15h\epsilon_a^2 - 3h^2\epsilon_a + 7h^3}{30}[\partial_{r^*}^3 \Psi]_a \\
& - \frac{5\epsilon_a^4 - 20h\epsilon_a^3 - 6h^2\epsilon_a^2 + 28h^3\epsilon_a + 23h^4}{120}[\partial_{r^*}^4 \Psi]_a \\
& - \frac{h(5\epsilon_a + h)}{5}[\partial_{r^*} \partial_t \Psi]_a - \frac{h(\epsilon_a - h)(5\epsilon_a + 7h)}{10}[\partial_{r^*}^2 \partial_t \Psi]_a \\
& + \frac{h^2(\epsilon_a - 7h)}{10}[\partial_{r^*} \partial_t^2 \Psi]_a - \frac{h(5\epsilon_a^3 + 3h\epsilon_a^2 - 21h^2\epsilon_a - 23h^3)}{30}[\partial_{r^*}^3 \partial_t \Psi]_a \\
& - \frac{h^2(\epsilon_a^2 - 14h\epsilon_a - 23h^2)}{20}[\partial_{r^*}^2 \partial_t^2 \Psi]_a + \frac{h^3(7\epsilon_a + 23h)}{30}[\partial_{r^*} \partial_t^3 \Psi]_a.
\end{aligned}
\quad (81)
$$

The jump conditions in the tortoise $r^*$ relate to those previously computed in the $r$ variable (the relations for mixed derivatives $(r^*, t)$ are easily inferred)

$$ [\Psi, r^*] = f_{r_u} [\Psi, r], \quad (82) $$

$$ [\Psi, r^* r^*] = f_{r_u} f'_{r_u} [\Psi, r] + f_{r_u}^2 [\Psi, rr], \quad (83) $$

$$ [\Psi, r^* r^* r^*] = f_{r_u} (f'^2 + ff'')_{r_u} [\Psi, r] + 3f_{r_u}^2 f'_{r_u} [\Psi, rr] + f_{r_u}^3 [\Psi, rrr], \quad (84) $$

$$
\begin{aligned}
[\Psi, r^* r^* r^* r^*] = f_{r_u} (f'^3 + 4ff'f'' + f^2f''' )_{r_u} [\Psi, r] &+ f_{r_u}^2 (7f'^2 + 4ff'')_{r_u} [\Psi, rr] \\
&+ 6f_{r_u}^3 f'_{r_u} [\Psi, rrr] + f_{r_u}^4 [\Psi, rrrr] .
\end{aligned}
\quad (85)
$$

## 5. Numerical implementation

Waveforms at infinity and at the particle position at first order are to be found in [21], as well as comparisons with other methods. Herein we are concerned on the numerical improvement. To this end, we have considered a distant observer, located at $r^* = 400(2M)$. The observer is reached by a pulse produced by a Gaussian, time-symmetric perturbation

$$ \Psi(r^*, t)_{t=0} = \exp[-(r^* - r_0^*)^2], \quad (86) $$

$$ \partial_t \Psi(r^*, t)_{t=0} = 0. \quad (87) $$
Fig. 5, obtained for $r_{u0} = 5(2M)$, shows the waveform produced in the homogeneous case. The convergence rate is computed as $(\epsilon^{(n)}(\xi))$ is the unknown error function of order $\approx 1$)

$$n = \log \left| \frac{\Psi(4h) - \Psi(2h)}{\Psi(2h) - \Psi(h)} \right| / \log(2) + \log \left| \epsilon^{(n)}(\xi) \right| / \log(2). \quad (88)$$

<figure><img src="image_6.png" /><figcaption>Figure 5: The waveform, $r_{u0} = 5(2M)$, of a Gaussian, time-symmetric initial pulse. The observer is located at $r^* = 400(2M)$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Convergence rates of the fourth and second order algorithms, $r_{u0} = 5(2M)$.</figcaption></figure>

Fig. 6, obtained for $r_{u0} = 5(2M)$, shows the fourth and second order convergence rates (we remind that the first order code [21] includes empty cells dealt at second order).
## 6. Conclusions

We have presented a fourth order novel integration method in time domain for the Zerilli wave equation. We have focused our attention to the even perturbations produced by a particle plunging in a non-rotating black hole. For cells crossed by the particle world line, the forward time wave function value at the upper node of the $(t, r^*)$ grid cell is obtained by the combination of the preceding node values of the same cell, analytic expressions related to the jump conditions, and the values of the wave function at adjacent cells. In this manner, the numerical integration does not deal directly nor with the source term and the associated singularities, nor with the potential term. In short, the direct integration of the wave equation is avoided. For other cells, we refer instead to already published approaches [15].

The scheme has also been applied to circular and eccentric orbits and it will be object of a forthcoming publication.

## Acknowledgements

The referees are thanked for careful reading and suggestions. The authors wish to acknowledge the FNAK (Fondation Nationale Alfred Kastler), the CJC (Confédération des Jeunes Chercheurs) and all organisations which stand against discrimination of foreign researchers.

## Appendix

Through Eq. 60, we have determined the value of $\Psi$ at the upper node of the cell as function of the analytic jump conditions and of the time derivatives of the wave function up to fourth order. The derivatives are evaluated at the point $b$ and weighted by five coefficients $c_0, c_1, c_2, c_3$ and $c_4$. Afterwards, the derivatives are converted into a linear combination of the wave function values taken on points at the left and right sides of the trajectory. Indeed, Eq. 65 represents such a system of linear equations. By injection of Eqs. 56-59 into Eq. 65, we get

$$ 
\begin{aligned}
& A_2 T_{\alpha_2}^{(0,0)} \Psi_b^- + A_2 T_{\alpha_2}^{(0,1)} \partial_t \Psi_b^- + A_2 T_{\alpha_2}^{(0,2)} \partial_t^2 \Psi_b^- + \dots + A_2 T_{\alpha_2}^{(1,3)} \partial_{r*} \partial_t^3 \Psi_b^- && c_0 \Psi_b^- \\
& + && + \\
& A_4 T_{\alpha_4}^{(0,0)} \Psi_b^- + A_4 T_{\alpha_4}^{(0,1)} \partial_t \Psi_b^- + A_4 T_{\alpha_4}^{(0,2)} \partial_t^2 \Psi_b^- + \dots + A_4 T_{\alpha_4}^{(1,3)} \partial_{r*} \partial_t^3 \Psi_b^- && c_1 \partial_t \Psi_b^- \\
& + && + \\
& \vdots && + \\
& + && c_3 \partial_t^3 \Psi_b^- \\
& \mathcal{N}_4^R T_{\nu_4^R}^{(0,0)} \Psi_b^- + \mathcal{N}_4^R T_{\nu_4^R}^{(0,1)} \partial_t \Psi_b^- + \mathcal{N}_4^R T_{\nu_4^R}^{(0,2)} \partial_t^2 \Psi_b^- + \dots + \mathcal{N}_4^R T_{\nu_4^R}^{(1,3)} \partial_{r*} \partial_t^3 \Psi_b^- && c_4 \partial_t^4 \Psi_b^- \\
& = && c_2 \partial_t^2 \Psi_b^- , (89)
\end{aligned}
$$

where $T_p^{(n,m)}$ represent the Taylor series coefficients at $p$ in the neighbourhood of $b$ and the indexes correspond to $n^{th}$ space and $m^{th}$ time derivatives. The wave function at $p$ is thus given by

$$ \Psi_p^{\pm} = \sum_{n+m \le 4} T_p^{(n,m)} \partial_{r^*}^n \partial_t^m \Psi_b^{\pm} + \mathcal{O}(h^5) . \qquad (90) $$
An example shows the procedure which is applicable to all cases. We pick the node $\alpha_2$, Eq. 55, where $T_{\alpha_2}^{(n,m)} = (-1)^n \frac{(2h-\epsilon_b)^n}{n!}$ and remind that $T_p^{(0,0)} = 1 \forall p$. By grouping the derivatives, we get

$$ 
\begin{aligned}
& \left( A_2 T_{\alpha_2}^{(0,0)} + A_4 T_{\alpha_4}^{(0,0)} + \cdots + N_4^R T_{\nu_4^R}^{(0,0)} \right) \Psi_b^- \\
& + \qquad c_0 \Psi_b^- \\
& \left( A_2 T_{\alpha_2}^{(0,1)} + A_4 T_{\alpha_4}^{(0,1)} + \cdots + N_4^R T_{\nu_4^R}^{(0,1)} \right) \partial_t \Psi_b^- \\
& + \qquad c_1 \partial_t \Psi_b^- \\
& \left( A_2 T_{\alpha_2}^{(0,2)} + A_4 T_{\alpha_4}^{(0,2)} + \cdots + N_4^R T_{\nu_4^R}^{(0,2)} \right) \partial_t^2 \Psi_b^- \\
& + \qquad c_2 \partial_t^2 \Psi_b^- \quad . \quad (91) \\
& + \qquad c_3 \partial_t^3 \Psi_b^- \\
& \vdots \\
& + \qquad c_4 \partial_t^4 \Psi_b^- \\
& \left( A_2 T_{\alpha_2}^{(1,3)} + A_4 T_{\alpha_4}^{(1,3)} + \cdots + N_4^R T_{\nu_4^R}^{(1,3)} \right) \partial_{r*} \partial_t^3 \Psi_b^- \\
\end{aligned}
$$

By identification, we obtain a linear system, that is cast in the form

$$ \underbrace{\begin{pmatrix} 1 & \cdots & 1 & \cdots & 1 & \cdots & 1 & \cdots & 1 \\ T_{\alpha_1}^{(0,1)} & \cdots & T_{\beta_1^L}^{(0,1)} & \cdots & T_{\gamma_2^L}^{(0,1)} & \cdots & T_{\mu_3^L}^{(0,1)} & \cdots & T_{\nu_4^R}^{(0,1)} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ T_{\alpha_1}^{(0,4)} & \cdots & T_{\beta_1^L}^{(0,4)} & \cdots & T_{\gamma_2^L}^{(0,4)} & \cdots & T_{\mu_3^L}^{(0,4)} & \cdots & T_{\nu_4^R}^{(0,4)} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ T_{\alpha_1}^{(1,0)} & \cdots & T_{\beta_1^L}^{(1,0)} & \cdots & T_{\gamma_2^L}^{(1,0)} & \cdots & T_{\mu_3^L}^{(1,0)} & \cdots & T_{\nu_4^R}^{(1,0)} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ T_{\alpha_1}^{(1,3)} & \cdots & T_{\beta_1^L}^{(1,3)} & \cdots & T_{\gamma_2^L}^{(1,3)} & \cdots & T_{\mu_3^L}^{(1,3)} & \cdots & T_{\nu_4^R}^{(1,3)} \end{pmatrix}}_{\mathbb{T}} \underbrace{\begin{pmatrix} A_2 \\ \vdots \\ B_1^L \\ \vdots \\ G_2^L \\ \vdots \\ M_3^L \\ \vdots \\ N_4^R \end{pmatrix}}_{\mathbb{P}} = \underbrace{\begin{pmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \\ c_4 \\ 0 \\ \vdots \\ 0 \\ \vdots \\ 0 \end{pmatrix}}_{\mathbb{C}} , (92) $$

where the upper indexes $(n,m)$ cover all combinations such that $n+m \le 4$. Finally, by inversion of the $\mathbb{T}$ matrix, the unknown terms of the $\mathbb{P}$ vector are identified.

## References

[1] Blanchet L Spallicci A and Whiting B 2011 *Mass and motion in general relativity* (Berlin: Springer)

[2] Droste J 1915 *Kon. Ak. Wetensch. Amsterdam* **23** 968. English translation: 1915 *Proc. Acad. Sci. Amsterdam* **17** 998

[3] Droste J 1916 *Het zwaartekrachtsveld van een of meer lichamen volgens de theorie van Einstein* Doctorate thesis Dir. Lorentz H A (Rijksuniversiteit van Leiden)

[4] Droste J 1916 *Kon. Ak. Wetensch. Amsterdam* **25** 163. English translation: 1917 *Proc. Acad. Sci. Amsterdam* **19** 197

[5] Schwarzschild K 1916 *Sitzungsber. Preuss. Akad. Wiss., Phys. Math. Kl.* 189. English translation with foreword by Antoci S and Loinger A 1999 arXiv:physics/9905030v1 [physics.hist-ph]

[6] Rothman T 2002 *Gen. Rel. Grav.* **34** 1541

[7] Regge T and Wheeler J A 1957 *Phys. Rev.* **108** 1063

[8] Zerilli F J 1969 *The gravitational field of a particle falling in a Schwarzschild geometry analyzed in tensor harmonics*, Doctorate thesis Dir. Wheeler J A (Princeton University)

[9] Zerilli F J 1970 *Phys. Rev. Lett.* **24** 737

[10] Zerilli F J 1970 *J. Math. Phys.* **11** 2203

[11] Zerilli F J 1970 *Phys. Rev. D* **2** 2141. Erratum, 1973 in *Black holes* Les Houches 30 July 31 August 1972 ed DeWitt C and DeWitt B (New York: Gordon and Breach Science Publ.)
[12] Lousto C O and Price R H 1997 *Phys. Rev. D* **56** 6439

[13] Martel K and Poisson E 2002 *Phys. Rev. D* **66** 084001

[14] Spallicci A 2011 in *Mass and motion in general relativity*, eds. Blanchet L Spallicci A Whiting B (Berlin: Springer) p 561 arXiv:1005.0611v1 [physics.hist-ph]

[15] Haas R 2007 *Phys. Rev. D* **75** 124011

[16] Lorenz L 1867 *Philos. Mag.* **34** 287

[17] de Donder T 1921 *La gravifique Einsteinienne* (Paris: Gauthier-Villars)

[18] Hunt B J 1991 *The Maxwellians* (New York: Cornell University Press)

[19] Barack L and Lousto C O 2005 *Phys. Rev. D* **72** 104026

[20] Aoudia S and Spallicci A to appear 12th *Marcel Grossmann Meeting* Paris 12-18 July 2009 ed. Damour T Jantzen R T and Ruffini R (Singapore: World Scientific) 2010 arXiv:1003.3107v3 [gr-qc]

[21] Aoudia S and Spallicci A 2011 *Phys. Rev. D* **83** 064029 arXiv:1008.2507v3 [gr-qc]

[22] Lousto C O 2005 *Class. Q. Grav.* **22** S543

[23] Sopuerta C F and Laguna P 2006 *Phys. Rev. D* **73** 044028

[24] Cañizares P and Sopuerta C F 2009 *Phys. Rev. D* **79** 064020

[25] Cañizares P and Sopuerta C F 2011 arXiv:1101.2526v1 [gr-qc]

[26] Jung J-H Khanna G and Nagle I 2007 Int. J. Mod. Phys. C **20** 1827

[27] Chakraborty D Jung J-H Khanna G 2011 arXiv:1103.1551v1 [physics.comp-ph]

[28] Zumbusch G 2009 *Class. Q. Grav.* **26** 175011

[29] Field S E Hesthaven J S and Lau S R 2009 *Class. Q. Grav.* **26** 165010

[30] Hopper S and Evans C R 2010 *Phys. Rev. D* **82** 084010

[31] Thornburg J 2010 arXiv:1006.3788v1 [gr-qc]

[32] Nagar A Damour T and Tartaglia A 2007 *Class. Q. Grav.* **24** S109

[33] Bernuzzi S and Nagar A 2010 *Phys. Rev. D* **81** 084056

[34] Moncrief V 1974 *Ann. Phys. (N.Y.)* **88** 323

[35] Lousto C O 2000 *Phys. Rev. Lett.* **84** 5251

[36] Lousto C O and Nakano H 2009 *Class. Q. Grav.* **26** 015007