The Fredholm Approach to Charactize Gravitational Wave Echoes

TL;DR

This paper introduces a novel Fredholm determinant-based method to compute and analyze gravitational wave echoes from compact objects, enabling efficient and accurate characterization of signals related to strong gravity physics.

Contribution

The paper presents a new Fredholm approach for perturbative and non-perturbative analysis of gravitational wave echoes, improving computational efficiency and clarity.

Findings

Method effectively characterizes individual echoes.

Perturbative and non-perturbative treatments are feasible.

Numerical implementation is resource-efficient and accurate.

Abstract

Gravitational waves are a sensitive probe into the structure of compact astronomical objects and the nature of gravity in the strong regime. Modifications of near-horizon physics can imprint on the late time ringdown waveform, leaving behind a train of echoes, from which useful information about new physics in the strong gravity regime can be extracted. We propose a novel approach to compute the ringdown waveform and characterize individual echoes perturbatively using the Fredholm determinants, which can be intuitively represented via a diagrammatical scheme. Direct non-perturbative treatments can also be easily implemented for some cases. Numerically, the method is also effective and accurate even for relatively low resources.

Equations

$$-\text{partialderivative}[2]\Phi (t,x)t+\text{partialderivative}[2]\Phi (t,x)x-V(x)\Phi (t,x)=0,$$

$$\Phi (t,x)\sim e^{i\omega (x-t)},x\to \infty .$$

$$\Phi (t,x)\sim e^{-i\omega (x+t)}+R(\omega )e^{i\omega (x-t)},x=-L,$$

$$\varphi (\omega ,x)={\int }_0^{\infty }\Phi (t,x)e^{i\omega t}\text{differential}t,$$

$${\varphi }^{\prime \prime }(\omega ,x)+({\omega }^2-V(x))\varphi (\omega ,x)=\mathcal{I}(\omega ,x),$$

$$\mathcal{I}(\omega ,x)=i\omega \Phi (0,x)-{\partial }_t\Phi (0,x).$$

$$\left(\text{derivative}[2]x+{\omega }^2\right)G(x,y)=\delta (x-y)$$

$$f(x)=\varphi (x)-\lambda \int \text{differential}yK(x,y)\varphi (y),$$

$$f(x)=\int \text{differential}yG(x,y)\mathcal{I}(y)$$

$$K(x,y)=G(x,y)V(y).$$

$$K(x,y)=\sum _{j=1}^{\infty }{\alpha }_j(x){\beta }_j(y).$$

$$Y_i=X_i-\lambda \sum _j{B}_{ij}{X}_j,$$

$$\begin{array}{rl}{X}_i& =\int \mathrm{d}x{\beta }_i(x)\varphi (x),Y_i=\int \mathrm{d}x{\beta }_i(x)f(x),\\ B_{ij}& =\int \mathrm{d}x{\beta }_i(x){\alpha }_j(x).\end{array}$$

$$X_i=\frac{{\Delta }_i(\lambda )}{\Delta (\lambda )},$$

$$\varphi (x)=f(x)+\frac{\lambda }{\Delta (\lambda )}\sum _j{\alpha }_j(x){\Delta }_j(\lambda ).$$

$$\sum _j{\alpha }_j(x){\Delta }_j(\lambda )$$

$$=\sum _j{\alpha }_j(x)\left|\begin{array}{ccccc}{c}_{11}& \cdots & Y_1& c_{1,j+1}& \cdots \\ c_{21}& \cdots & Y_2& c_{2,j+1}& \cdots \\ & & \cdots & & ⋮\end{array}\right|$$

$$=-\sum _j(-1{)}^j{\alpha }_j(x)\left|\begin{array}{ccccc}{Y}_1& c_{11}& \cdots & {\hat{c}}_{1j}& \cdots \\ Y_2& c_{21}& \cdots & {\hat{c}}_{2j}& \cdots \\ & & \cdots & & ⋮\end{array}\right|$$

$$=\left|\begin{array}{cccc}0& {\alpha }_1(x)& {\alpha }_2(x)& \cdots \\ Y_1& c_{11}& c_{12}& \cdots \\ Y_2& c_{21}& c_{22}& \cdots \\ & \cdots & & ⋮\end{array}\right|,$$

$$\varphi (x)=f(x)-\lambda \int \mathrm{d}yf(y)\frac{{\Delta }_K(x,y,\lambda )}{\Delta (\lambda )},$$

$${\Delta }_K(x,y,\lambda )=\left|\begin{array}{cccc}0& {\alpha }_1(x)& {\alpha }_2(x)& \cdots \\ {\beta }_1(y)& c_{11}& c_{12}& \cdots \\ {\beta }_2(y)& c_{21}& c_{22}& \cdots \\ & \cdots & & ⋮\end{array}\right|.$$

$$\Phi (t,x)={\int }_{-\infty +\beta i}^{\infty +\beta i}\varphi (\omega ,x)e^{-i\omega t}\frac{\text{differential}\omega }{2\pi },$$

$$\varphi (x)=f(x)+x{\int }_0^{1}y\varphi (y)\text{differential}y,$$

$$\varphi (x)$$

$$=f(x)+\frac{3}{2}x{\int }_0^{1}yf(y)\text{differential}y.$$

$$G(x,y)=\frac{e^{i\omega |x-y|}}{2i\omega }.$$

$$\alpha (x)=\frac{e^{i\omega x}}{2i\omega },\beta (y)=e^{-i\omega y}V(y).$$

$$\varphi (x)=f(x)+\int \text{differential}{x}_1\frac{e^{i\omega (x-x_1)}}{2i\omega }\frac{V(x_1)f(x_1)}{1-\frac{1}{2i\omega }\int \text{differential}yV(y)}.$$

$$\varphi (x)=f(x)+\int \frac{G(x,x_1)V(x_1)f(x_1)}{1-V_0/i\omega }\text{differential}{x}_1,$$

$$\varphi (x)$$

$$=f(x)+\lambda \int \text{differential}yf(y)\frac{{\sum }_{n=0}^{\infty }\frac{(-\lambda {)}^n}{n!}{A}_n(x,y)}{{\sum }_{n=0}^{\infty }\frac{(-\lambda {)}^n}{n!}{D}_n},$$

$$\begin{array}{rl}& A_n(x,y)=\int \prod _{i=1}^n\text{differential}{x}_i\left|\begin{array}{ccccc}K(x,y)& K(x,x_1)& K(x,x_2)& \cdots & K(x,x_n)\\ K(x_1,y)& K(x_1,x_1)& K(x_1,x_2)& \cdots & K(x_1,x_n)\\ & \cdots & & ⋮& \\ K(x_n,y)& K(x_n,x_1)& K(x_n,x_2)& \cdots & K(x_n,x_n)\end{array}\right|,\end{array}$$

$$\begin{array}{rl}& D_n=\int \prod _{i=1}^n\text{differential}{x}_i\left|\begin{array}{ccccc}K(x_1,x_1)& K(x_1,x_2)& K(x_1,x_3)& \cdots & K(x_1,x_n)\\ K(x_2,x_1)& K(x_2,x_2)& K(x_2,x_3)& \cdots & K(x_2,x_n)\\ & \cdots & & ⋮& \\ K(x_n,x_1)& K(x_n,x_2)& K(x_n,x_3)& \cdots & K(x_n,x_n)\end{array}\right|.\end{array}$$

Full text

1. I Introduction

2. II Fredholm approach of ringdown waveforms

3. II.1 The Fredholm integral equation

4. II.2 Examples

5. II.3 Perturbation scheme

6. II.4 Convergence and error estimate

7. III Echoes in the Fredholm approach

8. III.1 Delta potential with a mirror

9. III.2 Double delta potential

10. III.3 Perturbation scheme

11. IV Diagrammatical representation

12. V Numerical results in the Fredholm approach

13. VI Conclusions

I Introduction

The recent detections of gravitational waves (GWs) generated from some of the most extreme astronomical processes Abbott et al. (2016, 2017) have opened up a new window to probe the strong gravity regime and extreme astrophysical environments, providing fresh opportunities to test the fundamental laws of physics at the deepest levels Yunes et al. (2016).

Black holes are the most profound prediction of General Relativity (GR), which can be formed as the end states of gravitational collapses or other extreme processes. They are ideal playgrounds to theoretically test the principles and theories of quantum gravity, and also have played important roles in model astronomy even before the detections of GWs. The defining feature of a black hole is the existence of the event horizon, which is one-way membrane for all matter sources and thought to be of paramount importance to understand the fundamental properties of quantum gravity. However, we can not directly observe the event horizon itself due to its very nature, not even in the recent GW events Cardoso and Pani (2017a); Abramowicz et al. (2002). Indeed, the exact physical nature of the event horizon is still being heavily debated: When quantum corrections are included, it has been suggested that at presumably the Planck distance away from the event horizon there may exist a firewall that burns any infalling observer, contrary to the classical picture that the observer does not encounter anything extraordinary when passing through the horizon Almheiri et al. (2013). Whilst for a classical black hole the horizon in the tortoise coordinate is at negative infinity, in the firewall proposal the black hole has a surface structure immediately outside the would-be event horizon at some finite tortoise coordinate associated with the Planck length. Similarly, if there exist horizonless exotic compact objects (ECOs) whose spacetimes resemble those of black holes outside some small distances away from the would-be the event horizons, the GWs generated from the corresponding processes of those ECOs will also fit the current observational data Cardoso and Pani (2017b, a). Various ECO models have been proposed such as boson stars Schunck and Mielke (2003); Wheeler (1955), fuzzballs Mathur (2005), gravastars Mazur and Mottola (2004), and many others Holdom and Ren (2017); Barcelo et al. (2011); Barceló et al. (2015). A firewall or an ECO surface outside the would-be horizon can effectively introduces a barrier in the effective potential of the metric perturbation equation or a reflective boundary at a finite tortoise coordinate Cardoso and Pani (2017a). GWs can be trapped between the photosphere potential barrier and the reflective boundary, bounce back and forth within the trap and then transmit and propagate to the observer far away. This leads to a sequence of secondary pulses, dubbed echoes, after the initial ringdown pulse Cardoso et al. (2016a); Cardoso and Pani (2017a); Cardoso et al. (2016b).

On the other hand, if gravity is modified in the strong field regime, which, despite precise GR tests in the solar system and other weak field environments Will (2014), has not been well constrained, extra features may arise in the effective potential of the perturbations around the black hole. For example, in mass-varying massive gravity Huang et al. (2012), because of the large enhancement of the graviton mass near the black hole horizon, the graviton potential can induce an extra potential barrier to the left of the photosphere potential peak, and thus GWs can also be trapped in between, leading to echoes in the late time ringdown waveform Zhang and Zhou (2018). Wormhole scenarios can also produce effective potentials with extra barriers Cardoso et al. (2016a); Wang et al. (2018).

Therefore, GW echoes might be a sensitive probe of near-horizon new physics, and recently there have been a lot of interests in studying the phenomenological and theoretical aspects of GW echoes Barceló et al. (2017); Price and Khanna (2017); Nakano et al. (2017); Maselli et al. (2017); Völkel and Kokkotas (2017a); Bueno et al. (2018); Völkel and Kokkotas (2017b); Mark et al. (2017); Wang and Afshordi (2018); Pani and Ferrari (2018); Oshita and Afshordi (2019); Mannarelli and Tonelli (2018); Tsang et al. (2018); Testa and Pani (2018); Correia and Cardoso (2018); Burgess et al. (2018); Urbano and Veermäe (2019); Nielsen et al. (2019); Lo et al. (2019); Wang et al. (2018, 2019a); Konoplya et al. (2019); Conklin et al. (2018); Wang et al. (2019b); Li and Piao (2019); Ren (2019); Gálvez Ghersi et al. (2019); Cardoso et al. (2019); Uchikata et al. (2019); Saraswat and Afshordi (2019); Tsang et al. (2019); Maggio et al. (2019). Based some simple templates, it has been claimed that evidence of the existence of echoes is already contained in the current GW data Abedi et al. (2017a, b). However, for the current data, the statistical significance of the evidence is not sufficiently high to be a convincing detection Ashton et al. (2016); Westerweck et al. (2018). In Price and Khanna (2017), it is pointed out that for echoes coming from spacetime reflection conditions, the later ones are not simple repetition of the former in the sequence, which suggests that more sophisticated templates are generally desirable to accurately extract the echo signals. In Wang et al. (2018, 2019b); Li and Piao (2019), generic templates of non-equal interval echoes and possible scenarios that lead to un-equal intervals are initiated so as to accommodate more generic classes of ECOs or exotic scenarios to maximize possible information extraction in the GW data analysis.

In Mark et al. (2017), near-horizon exotic features are parametrized by a frequency dependent reflectivity at a fixed radius, and the full Green’s function is then split into the Green’s function of the corresponding black hole plus an extra term that can be expanded as a power series of the product of the reflectivities at the photosphere and the boundary and the phase delay between them. Consequently, the waveform observed far way can be re-processed as a power series of this product convoluted with the frequency domain wave amplitude at the would-be horizon. In this way, different terms in this series generates different individual echoes, which formalizes the physical picture of waves bouncing back and forth in the cavity formed by the photosphere barrier and the reflective boundary. Ref Correia and Cardoso (2018) proposed another another way to re-process the individual echoes perturbatively using the Dyson series. In this approach, the Green’s function for the free wave equation (that is, without any potential) with the reflectively boundary is utilized instead, and the wave equation is solved with the Dyson series, a series of multi-dimensional integrals of the product of the Green’s functions and the potentials. The individual echoes can be obtained via further expansions in terms of the reflectivity of the boundary.

In this paper, we propose yet another method to efficiently compute the ringdown waveforms and intuitively re-process the individual echoes. Our method is based the Fredholm approach to solve the wave equation with the Fredholm determinants. An intuitive diagrammatical scheme can be utilized to evaluate two kinds of the Fredholm determinants. One kind of them are obtained by constructing all possible loop diagrams connecting all the vertices, and the other kind are obtained by cutting one of the propagators in each of the loop diagrams. Introducing a new kind of vertex, individual echoes can be re-processed as diagrams containing different numbers of the new vertices. For many cases, non-perturbative solutions can also be straightforwardly obtained without the need of re-summations of the perturbative series. Numerically, waveforms can be obtained in a simple yet efficient manner.

The paper is organized as follows. In Section II, we will introduce the Fredholm formalism to compute the waveforms, perturbatively expand the solution with the Fredholm determinants, prove the convergence of the perturbative expansion and put estimates on the errors. In Section III, we will show how to non-perturbatively obtain the re-summed echo waveforms for a couple of toy models and how to extract individual echoes perturbatively. In Section IV, we will introduce our diagrammatical scheme which can be used to compute generic waveform or isolate individual echoes. In Section V, we will numerically compute the waveforms in this formalism, check its reliability and show its efficiency.

II Fredholm approach of ringdown waveforms

In the ringdown phase, the metric perturbations for a spherical compact star can be expanded with spherical harmonics $Y_{lm}(\theta,\phi)$, and for a given $l$ mode the linear perturbation $\Phi(t,x)$ satisfies a wave equation of the form

[TABLE]

where the effective potential $V(x)$ depends on the background metric, the spherical harmonic number $l$ and the spin of the perturbation mode. $x$ is the tortoise (radial) coordinate that makes the background metric conformally flat for the temporal and radial coordinates. In order to obtain the ringdown waveform, the following asymptotic behavior is imposed at spatial infinity

[TABLE]

On the left hand side, to capture various boundary conditions of potential ECOs, we impose the following boundary condition

[TABLE]

where $R(\omega)$ is the frequency dependent reflectivity and $L$ is to the left of the photosphere potential barrier. When the left hand side is also open ($x\to-\infty$), similar to the case of the Schwarzschild black hole, we may set $R(\omega)=0$ and $L\to\infty$. Equivalently, we may represent the effects of the extra potential barriers with some specific $R(\omega)$ at $x=-L$.

After Laplacian transform

[TABLE]

the wave equation in frequency domain is given by

[TABLE]

where ′ stands for a derivative with respect to $x$ and

[TABLE]

In the following we will often suppress the $\omega$ dependence $\phi(x)\equiv\phi(\omega,x)$. To get the solution in the time domain, we can solve Eq. (5) for different frequencies and perform the inverse Laplacian transform.

II.1 The Fredholm integral equation

In this section, we shall solve Eq. (5) by the Fredholm method. Suppose that $G(x,y)$ is the Green’s function of

[TABLE]

with appropriate boundary conditions. Then Eq. (5) can be cast as the Fredholm integral equation,

[TABLE]

where

[TABLE]

and the Kernel function is

[TABLE]

For later convenience, an explicit order parameter $\lambda$ has been introduced to count the expansion order in terms of interaction $V(y)$, which may be set to 1 in the end results. Assuming that $\{\alpha_{i}(x)\}$ and $\{\beta_{i}(y)\}$ are appropriate sets of base functions for radial variables $x$ and $y$ respectively, the kernel function can be expanded as

[TABLE]

With this expansion, Eq. (8) can be then written as a system of linear algebraic equations

[TABLE]

where we have defined

[TABLE]

The solution of Eq. (12) is given by

[TABLE]

where $\Delta(\lambda)=\det(I-\lambda B)$, $B=(B_{ij})$ and $\Delta_{i}(\lambda)$ is almost ${\Delta}(\lambda)$ but with the $i$-th column of $I-\lambda B$ replaced by $Y=(Y_{1},Y_{2},\dots)^{T}$. Plugging the solution back into Eq. (8), we find that

[TABLE]

Notice that $\sum_{j}\alpha_{j}(x)\Delta_{j}(\lambda)$ can be rewritten as

[TABLE]

where $c_{ij}$ are the matrix elements of $I-\lambda B$ and $\hat{c}_{1j}$, $\hat{c}_{2j}$, etc. denote missing elements in the corresponding places of the determinant. Therefore, the solution in the frequency domain is given by

[TABLE]

where

[TABLE]

Then the solution in the time domain can be written as

[TABLE]

where $\beta$ is chosen such that all the singularities of $\phi(\omega,x)$ are below the integration path in the complex $\omega$ plane.

II.2 Examples

Let us do a couple of examples to demonstrate how this formalism works. A trivial example would be when $K(x,y)=xy$ within the interval $y\in[0,1]$ and $K(x,y)=0$ otherwise.

[TABLE]

where we have set $\lambda=1$. For this case, the solution is simply

[TABLE]

Next, let us consider a system with a Schwarzschild-like boundary condition, for which case the free Green’s function is

[TABLE]

We seek the asymptotical wave solution observed far away, that is, we want to evaluate the waveforms at large and positive $x$, where the observer is located, so we have $|x-y|=x-y$. Then the kernel in this case is manifestly separable, we can set $K(x,y)=\alpha(x)\beta(y)$ with

[TABLE]

Thus the far-field solution is not much different from the trivial case above, and we can get

[TABLE]

To compare with the example used in Correia and Cardoso (2018), let us restrict to the Dirac delta potential $V(x_{1})=2V_{0}\delta(x_{1}-y)$, with $V_{0}$ being a positive constant, and we have

[TABLE]

and the pole at $\omega=-iV_{0}$ gives the quasi-normal mode of the system. This is the same as the re-summed Dyson series (see Correia and Cardoso (2018))

Therefore, for cases where the kernel function $K(x,y)$ can be separated with a sum of finite terms, the Fredholm approach can directly give the nonperturbative solution. We want to emphasize that Eq. (11) with an infinite summation can always be achieved, but an exact separation with a finite summation is not always possible. However, numerically, the kernel can be always approximated with a finite separable sum for a given precision, as we will see shortly.

II.3 Perturbation scheme

We can also obtain perturbative solutions in the Fredholm approach by expanding both $\Delta_{K}(x,y,\lambda)$ and ${\Delta}(\lambda)$ in terms of $\lambda$, and we have

[TABLE]

where the Fredholm determiants $A_{n}(x,y)$ and $D_{n}$ are given by

[TABLE]

[TABLE]

We will show why the expansion of $\Delta_{K}(x,y,\lambda)$ and ${\Delta}(\lambda)$ can be achieved in terms of $A_{n}(x,y)$ and $D_{n}$ respectively in Section IV via a diagrammatical scheme. In doing so, one obtains a Pade series in terms of $\lambda$, which is a double Taylor expansion both in the numerator and the denominator and whose coefficients are related to the eigen-polynomials of transition matrix $B$. As we shall later, if the kernel contains some exotic features near the would-be horizon such as a “mirror” with a frequency dependent reflectivity $R(\omega)$, one can also get a power series in $R(\omega)$ by re-processing the $\lambda$ Pade series. Thus, we have reduced the problem of solving the differential equation Eq. (5) to evaluating the determinants and integrations of the kernel $K(x,y)$. As we shall see later, these determinants can be conveniently computed with the diagrammatical method outlined in Section IV.

II.4 Convergence and error estimate

To numerically calculate the waveforms to a desired accuracy, we can choose a sufficiently large interval $[a,b]$ for the integration in the Fredholm integral equation and truncate the Pade expansion to a sufficiently large order $N$. In the following, we shall prove the convergence of the Pade series and then estimate of the errors arising from the different sources of approximation.

First, note that Eq. (31) and Eq. (32) allow us to use the Hadamard inequality to find the upper bounds:

[TABLE]

where $\tilde{M}$ is the upper bound of $|K_{N}(x,y)|$. So the numerator of Eq.(30) is bounded by

[TABLE]

Since $\frac{z^{n}}{n!}<e^{z}$ for positive number $z$, we get

[TABLE]

Since $\lim_{n\rightarrow\infty}(n!)^{1/n}/(n+1)^{1/2}=\infty$, we see that when $\Delta_{K}(x,y,\lambda)$ is expanded in terms of $A_{n}(x,y)$, the convergence radius for $\lambda$ is infinite. Of course, a similar argument holds for $\Delta(\lambda)$. This means that the solution we have constructed in Eq. (30) is convergent, and the limit of $N\to\infty$ leads to the exact solution.

Setting $\lambda=1$, our approximate solution is given by

[TABLE]

and the numerical errors can be bounded the following three parts:

[TABLE]

where $\Delta^{(N)}_{K}(x,y)$ and $\Delta^{(N)}$ denote the $N$-th terms of $\Delta_{K}(x,y)$ and $\Delta$ respectively.

To estimate the errors introduced by truncating the ${\lambda}$ expansion to order $N$, we note that Eq. (II.4) implies that $|\Delta_{K}(x,y)-\Delta^{(N)}_{K}(x,y)|$ is bounded by

[TABLE]

Also, we can show that $|\Delta-\Delta^{(N)}|$ is bounded $\epsilon_{N}$ in a similar manner, particularly suppressed by $1/\sqrt{N!}$. Therefore, one finds

[TABLE]

which is suppressed by $1/\sqrt{N!}$.

To estimate the errors arising from numerical integrations in a finite interval, we need to get a handle of how integrations in the infinite intervals $(-\infty,a)$ and $(b,\infty)$ will contribute to $A_{n}$ and $D_{n}$. For concreteness, we shall do the estimate on the case of the Schwarzschild effective potential, but similar estimates apply for the cases of ECOs as long as the effective potential still falls off sufficiently fast away from the peak:

[TABLE]

where $l$ is the $l$ from the spherical harmonics, $r(y)$ is implicitly determined by $y=r(y)+\ln({r(y)-1})$, and in this subsection we have restricted to the scalar perturbation and chosen the units such that $2GM=1$ ($G$ being the Newton’s gravitational constant and $M$ being the mass of the black hole).

Obviously, we need to choose the interval $[a,b]$ such that $V(y=a,b)$ are sufficiently small. When $r(y)$ is close to 1, which corresponds to $y$ being large and negative, we have

[TABLE]

which leads to

[TABLE]

On the other hand, when $r$ is large, which corresponds to the case when $y\to\infty$, we have

[TABLE]

which leads to

[TABLE]

Note also that the maximum value of $G(x,y)=\frac{e^{i\omega|x-y|}}{2i\omega}$ is just $w=\frac{1}{2|\omega|}$ and

[TABLE]

which can be bounded by the Hadamard inequality. Denoting $A^{[a,b]}_{n}(x,y)$ as $A_{n}(x,y)$ but with the integration limiting to the interval $[a,b]$, we have

[TABLE]

where we have defined

[TABLE]

and chosen $a$ and $b$ such that $\tilde{V}-\epsilon(a,b)>\frac{\tilde{V}}{2}$. So the finite interval error in the numerator is bounded by

[TABLE]

where $\tilde{f}(\tilde{m})\equiv\tilde{m}+\sqrt{2}\tilde{m}^{2}+\frac{\sqrt{6}}{2}\tilde{m}^{3}+\frac{\sqrt{5}}{2}\tilde{m}^{4}e^{\tilde{m}^{2}/2}$. For the denominator, similarly, we have

[TABLE]

Therefore, we see that the errors introduced by using the finite interval $[a,b]$ to compute $A_{n}(x,y)$ and $D_{n}$ is bounced by $\epsilon(a,b)$, which decreases exponentially with $a$ ($a$ being negative) but only inversely proportionally with $b$. This suggests that while a relatively small $|a|$ is sufficient, it is imperative to choose relatively large $b$.

Combining Eq. (II.4) and Eq. (54), we can get

[TABLE]

We see that the bound is independent of $a$ and $b$. Therefore, the first two terms in Eq. (II.4) are bounded by

[TABLE]

where $m_{2}$ is the maximum value of $f(y)$ in $(-\infty,a)\cup(b,\infty)$. The bound is also proportional to $e^{a}/{b}$.

These are of course the upper bounds, and in actual calculations there may be cancellations, which will lead to better accuracies.

III Echoes in the Fredholm approach

Echoes in the ringdown waveform are signatures of new physics in the strong gravity regime. They can arise if the effective potential is modified with a reflective mirror or an extra potential barrier to the left of the photosphere peak. In the following, we shall first demonstrate how to re-process the echoes in the Fredholm approach with two separable examples, which allows for direct nonpeturbative treatments, and then show how to extract echo signals systematically by perturbation theory.

III.1 Delta potential with a mirror

First, we consider a Dirac delta potential

[TABLE]

with a reflective mirror at $x=-L$. In this case, using Eq. (77), the integral kernel can be written as

[TABLE]

Using the property of the delta function, the kernel can be written as

[TABLE]

with

[TABLE]

Calculating the transition matrix and setting $\lambda=1$, it is straightforward to get

[TABLE]

where $f(x)$ is simply the first term in the Born approximation

[TABLE]

The denominator of Eq. (62) can be expanded as

[TABLE]

where $\binom{n}{k}$ is the binomial coefficients “$n$ choose $k$” and we have defined reflectivity of the delta potential

[TABLE]

Substituting Eq. (66) into Eq. (62), we get

[TABLE]

which is a power series in terms of the reflectivity $R(\omega)$ of the mirror. The the $n$-th echo in the frequency domain is just the wave packet that has bounced $n$ time by the mirror, thus the term with the $n$-th power of $R(\omega)$ Mark et al. (2017), and therefore reads

[TABLE]

for $n\geq 1$; $\phi_{0}$ is given by Eq. (III.1) with $R_{\delta}^{-1}$ setting to 0. This is in complete agreement with Correia and Cardoso (2018), which uses the Dyson series approach.

III.2 Double delta potential

Now, we want to compute the echoes for the double delta potential

[TABLE]

with an open boundary on the left hand side, that is, the left hand boundary is at $x=-L$ with $L\to\infty$. As we shall see momentarily, this is very much the same as the previous case in the Fredholm approach, while in the Dyson series approach it is more involved if one needs to get the full non-perturbative solution.

In this case, the kernel function is

[TABLE]

which can be written as

[TABLE]

with

[TABLE]

Going through steps similar to the case of the delta potential with a mirror, we can get

[TABLE]

which is the full non-perturbative solution in terms of $R$. The $n$-th echo signal is simply the $n$-th term in the expansion of Eq. (76) in terms of $R$. This is different from the Dyson series approach where one gets the perturbative solutions order by order, which seems to be difficult to be re-summed to a full non-perturbative form.

Notice that the phase delay factor $Re^{2i\omega|P|}$ appears in denominators of the terms in Eq. (76) and there are $R$ factors in numerators, so when expanded in terms of $R$, the $R^{n}$ term will come with either a phase delay factor $e^{2ni\omega|P|}$ or $e^{(2n+1)i\omega|P|}$. This basically reflects the physical intuition that the $n$-th echo comes from waves that are reflected $n$ or $n+1$ times between the two potential barriers.

III.3 Perturbation scheme

In the Fredholm formalism, we can also decompose the waveforms into different echo signals via a perturbative expansion. For the case where the boundary $x=-L$ is a reflective mirror. The Green’s function in Eq. 7 is modified to

[TABLE]

On the other hand, for the case where the potential can be split into two localized parts

[TABLE]

the kernel can be written as

[TABLE]

where often $R(\omega)$ just has a trivial $\omega$ dependence. For either case, the kernel function can be split into two parts

[TABLE]

Substituting this into the definition of $D_{n}$, we find

[TABLE]

where $D_{n}^{j}$ is the sum of all possible ways to change $j$ of the $n$ columns of $D_{n}$ from $\bar{K}$ to $Q$ and $A_{n}^{j}(x,y)$ is the sum of all possible ways to change $j$ of the $n+1$ columns of $A_{n}(x,y)$ from $\bar{K}$ to $Q$, namely,

[TABLE]

[TABLE]

Substituting all these into Eq. (30) and expanding perturbatively in terms of $R$, we get the waveform of the echoes. Explicitly, if we define

[TABLE]

($A_{-1}^{0}(x,y)$ being defined to be 0 for convenience), we can re-write

[TABLE]

where $c_{j}$ can be expressed by $a_{j}$ and $b_{j}$ as:

[TABLE]

and $c_{0}$ is defined as $c_{0}=a_{0}/b_{0}$.

Again the integrals of the determinants can be systematically computed by the diagrammatical method outlined in Section IV.

In summary, the $n$-th echo can be written as

[TABLE]

where $c_{n}$ is defined in Eq. (86).

IV Diagrammatical representation

In this section, we will show how $\Delta_{K}(x,y;\lambda)$ and $\Delta(\lambda)$ can be expanded with terms defined in Eq. (31) and Eq. (32) respectively, which can be elegantly done with a diagrammatic scheme. We will see that echoes can also be represented within this scheme.

This scheme is based on the realization that $K(x_{i},x_{j})$ can be viewed as a 2-point “vertex” while an integration $\int\mathrm{d}x$ can viewed as a “propagator”. In addition, for each diagram with an odd number of loops, we need to assign a minus sign, which is simply due to the antisymmetric nature of the determinant. Furthermore, by Eq. (11) and Eq. (13), we can also equate the integration of $K$ to the trace of matrix $B$, so we can directly translate the diagrams to traces of matrix $B$. It might be instructive to give a few examples:

[TABLE]

In this language, the $n$-dimensional determinant $D_{n}$ can be denoted as a sum of all possible allowed diagrams, each of which contains $n$ vertices connected $n$ propagators without loose ends (ie, each diagram being a multiplication of loops). For example, for $D_{4}$, we have

[TABLE]

This of course can be verified by the well-known result for the determinant. Since $\Delta(\lambda)$ is just $\det(I-\lambda B)$, we have

[TABLE]

and then $D_{4}$ is simply $\mathrm{d}^{4}\Delta(\lambda)/\mathrm{d}\lambda^{4}|_{\lambda=0}$, which agrees with Eq. (94). Note that a simple consistent check is that for each $D_{n}$ expansion the sum of coefficients should add up to zero.

To compute the waveform in a finite interval, we use the approximation $K_{N}(x,y)=\sum_{i=1}^{N}\alpha_{i}(x)\beta_{i}(y)$. In this case, we can truncate the $\lambda$ expansion of $\Delta(\lambda)$ at order $N$, as higher order $D_{n}$’s vanish identically.

For $\Delta_{K}(x,y;\lambda)$, the diagrams are not closed loops and extra “cutting rules” are needed, as it contains “external line” $x$ and $y$. To find these rules, we note that

[TABLE]

Conversely, $A_{n}(x,y)$ can be obtained by cutting out a propagator of $D_{n+1}$. (Of course, $A_{0}(x,y)$ can not be obtained by cutting, but it is simply $A_{0}(x,y)=K(x,y)$.) When doing the cutting, it is important to remember that we need to apply the Leibniz rule to cut each part of the diagram (eg, a diagram may have 3 parts, ie, 3 loops multiplied together), and for each diagram we cut exactly once. For example, by cutting a line in the diagrams of $D_{4}$, we can get

[TABLE]

where we have used the shorthand notation $\alpha\cdot B^{n}\cdot\beta=\sum_{i,j}\alpha_{i}(x)(B^{n})_{ij}\beta_{j}(y)$. The coefficients in front of the diagrams are obtained by

[TABLE]

where $a_{i}$ denotes a loop with $i$ vertices and $b_{i}$ denotes a line with $i$ vertices. Again, for each $A_{n}$ expansion, the sum of coefficients should add up to zero, as the cutting does not change the total number of diagrams.

Therefore, the waveform can be represented as

[TABLE]

Note that the cutting operation does not commute with the exponential map, so we need to compute the exponential first and then perform the cutting.

Having established the general diagrammatical formalism to compute waveforms, it is straightforward to extract the echoes form it. Note that, for scenarios with echoes, the kernel function may be split into two parts

[TABLE]

This suggests that we can add an extra 2-point vertex representing $Q(x,y)$:

[TABLE]

(In this notation, of course, the 2-point circle vertex represent $\bar{K}(x,y)$ or not $K(x,y)$.) Therefore, all the loop diagrams in the general formalism should be replaced with

[TABLE]

$\Delta(\lambda)$ then is expanded as

[TABLE]

with

[TABLE]

where we have used the short note $\bar{K}_{ij}$ and $Q_{ij}$ for $\bar{K}(x_{i},x_{j})$ and $Q(x_{i},x_{j})$. The cutting rules also apply, so we have

[TABLE]

with

[TABLE]

Therefore, to calculate the $n$-th echo signal, we should collect all the diagrams with $n$ 2-point solid vertices.

Let us summarize the diagrammatical rules that are needed in order to compute the waveform and echoes:

• $\bullet$

$\bar{K}(x_{i},x_{j})$ is represented as a (2-point) circle vertex;

• $\bullet$

$Q(x_{i},x_{j})$ is represented as a (2-point) solid vertex;

• $\bullet$

Diagrams with $n$ solid vertices contribute to the $n$-th echo waveform;

• $\bullet$

$\int\mathrm{d}x$ is represented as a “propagator”;

• $\bullet$

For each diagram with an odd number of loops, assign a minus sign;

• $\bullet$

$D_{n}$ includes all loop diagrams with $n$ vertices;

• $\bullet$

$A_{n}$ is obtained by cutting open one of the loops in $D_{n+1}$.

V Numerical results in the Fredholm approach

We have shown that the Fredholm approach can be used perturbatively to process individual echo waveforms and for separable cases to conveniently compute the exact waveform analytically. In this section, we shall demonstrate how to compute the waveform for generic (non-separable) cases.

To numerically compute the full waveform, the simplest way is simply to discretize the Fredholm integral equation Eq. (8) in the finite interval $[a,b]$ as

[TABLE]

where $N$ is an even number and we have defined

[TABLE]

Here we use the Simpson’s rule for numerical integration, as coded in the numerical factor $\eta_{j}$. Thus the waveform computation reduces to a linear algebraic problem, for which there exist a vast collection of fast numerical libraries to facilitate efficient and precision computations. To get the waveform in the time domain, we evaluate $\phi(\omega,x)$ for a large number of $\omega$’s and then we perform the inverse Laplace transformation (see Eq. (22)).

As an aside, we want to point out a simple application of the Fredholm formalism to the anti-scattering problem — extract the form of the potential via the waveform. With the method we introduced here, it is straightforward to re-construct the potential from the waveform in the frequency domain. A simple reformulation of Eq. (109) gives

[TABLE]

Defining matrix $C_{ij}=\eta_{j}G(x_{i},x_{j})\phi_{j}$ and vector $\tilde{\phi}_{i}=\frac{3N}{b-a}(\phi_{i}-f_{i})$, the potential vector is then simply given by

[TABLE]

where $C^{-1}$ is the inverse of matrix $C$. This is of course easy to solve. In the gravitational wave setting, this is not useful since we can only observe the waveform far away from the source. However, it may find its application in a laboratory setting for a system satisfying the Schroedinger equation.

In the following, we shall compute the waveforms of a few idealized models to demonstrate the efficiency and accuracy of our method. For all these simulations, the initial conditions we use are

[TABLE]

The models we consider are tabulated as follows,

[TABLE]

with

[TABLE]

where we have set $2GM=1$ and again $r(x)$ is obtained by inverting the relation $x=r+\ln({r-1})$ and in this section we shall set $l=2$ and $s=2$.

Fig. 1 is just the ringdown waveform generated by the Schwarzschild effective potential $V_{\rm gr}(l=2,s=2)$, observed at $x=30$. As we shall see shortly, even for a relatively narrow interval and small $N$, our method already generates very accurate waveforms.

Fig. 2 is a comparison between the two waveforms generated by Model (b) ($V(x)=V_{\rm gr}(x)+V_{\rm rec1}(x)$) and Model (e) ($V(x)=V_{\rm gr}(x)+V_{\rm rec2}(x)$), observed at $x=40$. Model (e) has an extra potential barrier to the left of the potential peak of Model (b). We see that Model (e) contains all the echoes of Model (b) and also has extra smaller echoes in between the bigger echoes of Model (b), as anticipated from the intuition of a propagating wave package transmitting and reflected between the potential barriers before reaching the observer to the right of the Schwarzschild potential barrier. In particular, the bigger echoes are due to the waves bouncing back and forward between the Schwarzschild potential peak and the rectangular potential at $x=-50$, while the smaller echoes come from waves bouncing back and forward between the rectangular potentials at $x=-50$ and $x=-80$. The latter ones are relatively smaller because they have to transmit through an extra potential barrier.

When there is a reflecting “mirror” at $x=-L$, we should use the Green’s function Eq. (77) with nonzero $R(w)$. For an illustration, we consider a simple frequency independent reflectivity

[TABLE]

where $r_{f}=1$ for the Dirichlet boundary condition and $r_{f}=-1$ for the Neumann boundary condition. The difference of these two boundary conditions is of course to flip the value of waveform $\phi(t,x)$ for echoes associated with odd powers of $R(\omega)$ (see Fig. 3).

In our numerical scheme, we have truncated the potentials for $x<a$ and $x>b$. We previously have shown that this is a consistency scheme by proving the convergence and estimating the errors. Numerically, the effectiveness of this approximation can also be verified by convergence studies, in which we do computations with different ranges of the interval $[a,b]$ and check whether the waveforms converge to a waveform with enlarged intervals. In Fig. 4, we plot the differences between the 3 waveforms generated by 3 different intervals in Model (a), ie, $V_{\rm gr}$ $(l=2$, $s=2)$. We see that the waveforms converge to the real values very quickly even for relatively small intervals.

Because of the good convergence properties, our scheme is insensitive to the boundary conditions. We also want to point out that the computation costs do not sensitively rely on the complexity of the effective potential, as we only need to solve a system of linear equations.

Finally, let us note that the Fredholm approach also presents a way to compute the quasi-normal modes, which are simply the solutions to $\det(B)=0$. For example, the fundamental quasi-normal mode for a rectangular potential

[TABLE]

can be analytically calculated and is $4.660-0.710i$, see Barausse et al. (2014). In our approach, we only need $N=26$ to reach the accuracy of 1%, ie, $\det(B)=10^{-2}$.

VI Conclusions

Apart from the confirmation of the most spectacular natural phenomena themselves, the newly arrival of the gravitational wave astronomy has allowed us to probe the astonishing properties of compact astronomical objects and gravitation itself, particularly in the strong field regime. While initially considered as rather exotic objects, black holes have long established their roles in modern astrophysics. The gravitational wave astronomy ushers in the exciting possibility of detecting potentially more surprising attributes of the black holes or even more exotic compact objects, leading to new discoveries in fundamental physics.

Echoes in the ringdown waveform are messengers of many near-horizon exotic new physics, and thus it is justified to carefully characterize and accurately describe the essential features of the waveforms containing them so as to help build reasonable templates to extract useful information in the observational data. We have proposed a new method to represent the individual echoes as terms of different orders in a Pade type perturbative expansion. Also, we have devised an intuitive diagrammatical scheme involving vertices, propagators and cutting rules to facilitate the perturbative expansion. We have also proven the convergence of the Pade expansion, estimated the errors and numerically demonstrated efficiency of the method. For some separable cases, our approach also allows for an easy computation of the exact non-perturbative wave solution.

Acknowledgments: We would like to thank Dong-Ming He, Zhi-Peng Li, Yun-Song Piao, Hao-Yu Shang, Ying-Jie Wang and Jun Zhang for helpful discussions. We are especially grateful to Wen-Bo Li for discussions on some mathematical aspects in the paper. SYZ acknowledges support from the starting grant from University of Science and Technology of China under grant No. KY2030000089 and is also supported by National Natural Science Foundation of China under grant No. GG2030040375.