Quantum Black Holes in the Sky

Black Holes are possibly the most enigmatic objects in our Universe. From their detection in gravitational waves upon their mergers, to their snapshot eating at the centres of galaxies, black hole astrophysics has undergone an observational renaissance in the past 4 years. Nevertheless, they remain active playgrounds for strong gravity and quantum effects, where novel aspects of the elusive theory of quantum gravity may be hard at work. In this review article, we provide an overview of the strong motivations for why"Quantum Black Holes"may be radically different from their classical counterparts in Einstein's General Relativity. We then discuss the observational signatures of quantum black holes, focusing on gravitational wave echoes as smoking guns for quantum horizons (or exotic compact objects), which have led to significant recent excitement and activity. We review the theoretical underpinning of gravitational wave echoes and critically examine the seemingly contradictory observational claims regarding their (non-)existence. Finally, we discuss the future theoretical and observational landscape for unraveling the"Quantum Black Holes in the Sky".


Introduction
Black holes (BHs) are very interesting "stars" in the Universe where both strong gravity and macroscopic quantum behavior are expected to coexist. Classical BHs in General Relativity (GR) have been thought to have only three hairs, i.e., mass, angular momentum, and charge, making observational predictions for BHs relatively easy [11,12] (compared to other astrophysical compact objects). For astrophysical BHs, due to the effect of ambient plasma, this charge is vanishingly small, leaving us with effectively two hairs for isolated black holes, with small accretion rates. In other words, finding conclusive deviations from standard predictions of these 2-parameter models, may be interpreted as fingerprints of a quantum theory of gravity or other possible deviations from GR. For example, the quasinormal modes (QNMs) of spinning BHs, which have been widely-studied over the past few decades (a subject often referred to as BH spectroscopy), only depend on the mass and spin of the Kerr BH (e.g., [13]). The ringdown of the perturbations of the BH is regarded as a superposition of these QNMs, and thus can be used to test the accuracy of GR predictions and no-hair theorem (e.g., see [14]). As a result, precise detection of QNMs from the ringdown phase (from BH mergers or formation) in gravitational wave (GW) observations may enable us to test the classical and quantum modifications to GR (e.g., [15]).
A concrete path towards this goal is paved through the study of "GW echoes", a smoking gun for near-horizon modifications of GR which are motivated from the resolutions of the proposed resolutions to the BH information paradox and dark energy problems [16,17]. The list of these models include wormholes [18], gravastars [19], fuzzballs [20], 2-2 holes [21], Aether Holes [17], Firewalls [16] and the Planckian correction in the dispersion relation of gravitational field [22,23].
Given their uncertain theoretical and observational status, GW echoes are gathering much attention from those who are interested in the observational signatures of quantum gravity, and the field remains full of excitement, controversy and confusion. In this review article, we aim to bring some clarity to this situation, from its background, to its current status, and into its future outlook.
The review article is organized as follows: In the next section, we provide builds the motivation to investigate the quantum signatures from BHs. In Sec. 3, we discuss theoretical models of quantum BHs, starting from the BH information loss paradox, and then its proposed physical resolutions that lead to observable signatures. In Sec. 4, we review how to predict the GW echoes from spinning BHs based on the Chandrasekhar-Detweiler (CD) equation, and also review the Boltzmann reflectivity model [23,29] for quantum black holes. Sec. 5 is devoted to the echo searches, where we summarize positive, negative, and mixed reported outcomes, and attempt to provide a balanced and unified census. In Sec. 6, we discuss the future prospects for advancement in theoretical and observational studies of quantum black holes, while Sec. 7 concludes the review article.
Throughout the article, we use the following notations: Furthermore, unless noted otherwise, we use the natural Planck units withh = c = 1 = G = 1.

Classical BHs
Schwarzschild and Kerr spacetimes are solutions of general relativity giving the spacetime configurations of BHs, which are the most dense objects in our universe. In some regions, matter accumulates and attracts more matter with its gravity which is classically always attractive. At the end, the force is so strong that even the light cannot escape from those regions, where then BHs form. The first and most important feature (the definition of BHs) is the formation of horizon. Inside the (event) horizon, all the light cones are directed into the singularity, and nothing can escape, unless it could travel faster than the speed of light. Therefore, horizons stand as the causal boundaries of BHs in Einstein's theory of Relativity.
Realistic BHs in the sky have different hairs (mass, spin and charge), and their dynamics share more complicated structure, thus, have different kinds of horizons. To list some of them, event horizons are defined as the boundaries where no light can escape to the infinite future. However, for a dynamically evolving BH, event horizons are teleological, i.e. we cannot predict them until we have the entire history of the spacetime. Apparent horizons, however, are predictable at a specific time without knowing the 5 of 84 future. Any surface has two null normal vectors and if expansion of both of them are negative, the surface is called "trapped". Apparent horizons are the outermost of all the trapped surfaces, which is why they are also known as the "marginally outer trapped surface".
Here is a simple example to distinguish these two horizons -we start with a Schwarzschild BH at time t 1 , now the event and apparent horizons coincide at the Schwarzschild radius. We throw a spherical null shell into the BH and let it cool down at t 2 . This process is perfectly described by Vaidya metric [30]. The apparent horizon changes immediately when the shell falls into the BH but the event horizon starts to expand earlier, even before the shell reaches it. It is because that after throwing the shell, the gravity of BH increases. Thus it is harder for light to escape from the BH to infinity. In other words, particles might be doomed to fall into a singularity, even before they had a chance to meet the infalling gravitating matter that is responsible for their fate. Therefore, the event horizon is modified earlier than the apparent horizon. While this result is counter-intuitive, it is a result of the formal definition of the event horizons, which requires the information about the entire history of spacetime, in particular, the future! Beyond the horizons, another intriguing trait of BHs is the curvature singularity, which sits at the centre of the BHs. Horizons can also be singular, but usually only coordinate singularities and (in classical General Relativity) removable by changing to a proper coordinate system. However, the singularities inside the BHs are where the general relativity breaks down and so far we do not have any good physics to describe them. We cannot chase the information lost into these singularities (using standard physics), which leads to the information paradox (more on this later).
Back in November 1784, John Michell, an English clergyman, advanced the idea that light might not be able to escape from a very massive object (at a fixed density). For example, light cannot escape from the surface of a star with the density of the sun, if it was 500 times bigger than the sun. Albert Einstein, later in 1915, developed general relativity. Soon after this, Karl Schwarzschild solved the Einstein vacuum field equation under spherical symmetry with a singular mass at the center, which was the first solution for BHs, the Schwarzschild metric.
While 20th century saw a golden age of general relativity with blooming of dozens of different BH solutions, the existence of BHs was not directly confirmed until one century later in 2015. LIGO-Virgo collaboration reported unprecedented detection of GWs from the binary BH merger events [31][32][33][34][35][36][37][38]. Numerical relativity is consistent with LIGO data at least up to quite near the horizon range. But the detection has not confirmed the existence of the horizons. We will discuss in this article how the detection opens a window for searching for quantum nature of the BHs beyond the general relativity.

Schwarzschild spacetime
The Schwarzschild spacetime was the first exact solution in the Einstein theory of general relativity. It models a static gravitational field outside a mass which has spherical symmetry, zero charge and rotation. Karl Schwarzschild found this solution in 1915, and four months later, Johannes Droste published a more concrete study on this independently. The metric in the Schwarzschild coordinate is: where M is the mass of the centre object, 2M is Schwarzschild radius and dΩ 2 = dθ 2 + sin 2 θdφ 2 is the metric on a 2-sphere. The metric describes gravitational field outside any spherical object without charges.
If the radius of the central object is smaller than the Schwarzschild radius, the object is then too dense to be stable, and will go through a gravitational collapse and form the Schwarzschild BH. Later in 1923, G.D.Birkhoff proved that any spherically symmetric solution of the vacuum Einstein field equation must be static and asymptotically flat. Hence, Schwarzschild metric is the only solution in 6 of 84 that case. For any static solution, the event horizon always coincides with the apparent horizon at r = 2M. In general relativity, Schwarzschild metric is singular at the horizon but , as stated above, this is only a coordinate artifact. That is to say, a free falling observer feels no drama going through the horizon. It takes the observer a finite amount of proper time but infinite coordinate time. Particularly, we can remove the singularity by a proper coordinate transformation. In contrast, the origin r = 0 is the intrinsic curvature singularity. The scalar curvature is infinite and the general relativity is no longer valid at this point.

Kerr spacetime
The Kerr spacetime [39], discovered by Roy Kerr, is a realistic generalization of the Schwarzschild spacetime. It describes the gravitational field of an empty spacetime outside a rotating object. The spacetime is stationary and has axial symmetry. The metric in the Boyer-Lindquist coordinate is: ∆ dr 2 + ρ 2 dθ 2 + r 2 + a 2 + 2Mra 2 ρ 2 sin 2 θ sin 2 θdφ 2 − 4Mra sin 2 θ where a = J/M, ρ 2 = r 2 + a 2 cos 2 θ, ∆ = r 2 − 2Mr + a 2 , Σ = (r 2 + a 2 ) 2 − a 2 ∆ sin 2 θ and ω = − g tφ g φφ = 2Mar Σ . The Cartesian coordinates can be defined as x = r 2 + a 2 sin θ cos φ, y = r 2 + a 2 sin θ sin φ, z = r cos θ. (4) There are two singularities easily reading from the coordinate where the g rr and g tt vanish. The first one gives r ± = M ± √ M 2 − a 2 corresponding to the horizon analog to the Schwarzschild metric. The larger root r + = M + √ M 2 − a 2 is the event horizon, while the other root is inner apparent horizon.
The second singularity is related to an interesting effect in the Kerr spacetime called frame-dragging effect: When reaching close to the Kerr BHs, the observers even with zero angular momentum (ZAMOs) will co-rotate with the BHs because of the swirling of spacetime from the rotating body. We assume that u α is the four-velocity of ZAMOs, and from the conservation of angular momentum g φtṫ + g φφφ = 0, where an overdot is differentiation with respect to the proper time of the observers τ. Thus, dφ dt = − g tφ g φφ . Because of this frame-dragging effect, there is a region of spacetime where static observers cannot exist, no matter how much external force is applied. This region is known as the "ergosphere" r ≤ M + √ M 2 − a 2 cos 2 θ.
The rotation also leads to another interesting feature, called "superradiance". That is, we can extract energy from scattering waves off the Kerr BHs. The exact formalism of superradiance is defined and discussed in Sec. 4.2. Finally, the Kerr spacetime also possesses a curvature singularity at the origin ρ 2 = r 2 + a 2 cos 2 θ. However, in contrast to Schwarzschild case, this singularity can be avoided since it is a ring at r=0 and θ = π/2, where z=0 and x 2 + y 2 = a 2 . In principle, observers can go through the ring without hitting the singularity. However, it is widely believed that the inner horizon, r − in Kerr spacetime is subject to an instability which would dim the analytic extension of Kerr metric beyond r − unphysical [40].

Blue-shift near horizon
As shown in the metric, different observers have different proper time. Hence, in the general relativity, the clocks at a gravitational field tick in a different speed in a different spacetime point. This is the blue(red)-shift effect, and it is extremely strong close to the dense object, especially near horizon. 7 of 84 Assuming static clocks in the Schwarzschild spacetime ds 2 = −dτ 2 = −(1 − 2M/r o )dt 2 , where τ is the proper (clock) time of an observer at distance r o . Hence, t is the proper time of an observer at infinity. The shifted wavelength λ o measured by observers at r o compared to observers at infinite is 2.1.4. Thermodynamics of Semi-classical BH Jacob Bekenstein and Stephen Hawking first proposed that the entropy of BHs is related to the area of their event horizons divided by the Planck area [41][42][43][44][45]. Furthermore, in 1974, Stephen Hawking showed that rather than being totally black, BHs emit thermal radiation at the Hawking temperature, T H = κ 2π , where κ is the surface gravity at the horizon [46][47][48]. This then lead to the celebrated Bekenstein-Hawking entropy formula S BH = A 4 [44,45], where A is the area of the event horizon. However, the nature of microstates of BHs that are enumerated by this entropy remains so far unknown. String theory associates it with higher dimensional fuzzball solutions, as discussed later in Sec. 3.5. Loop quantum gravity relates the quantum geometries of the horizon to the microstates [49]. Both these approaches can give the right Bekenstein-Hawking entropy, given specific assumptions and idealizations.
Interestingly, not only the entropy exists for the BHs, but also Brandon Carter, Stephen Hawking and James Bardeen [50] discovered the four laws of BH thermality analogous to the four laws of thermodynamics. The latter is presented in the parentheses.
• The zeroth law: A stationary BH has constant surface gravity κ. (A thermal equilibrium system has a constant temperature T H .) • The first law: A small change of mass for a stationary BH is related to the changes in the horizon area A, the angular momentum J, and the electric charge Q: dM = κ 8π dA + ΩdJ + ΦdQ, where Ω is the angular velocity and Φ is the electrostatic potential (Energy conservation: dE = TdS − PdV − µdN).
• The second law: The area of event horizon A never decreases in general relativity (The entropy of isolated systems never decreases). • The third law: BHs with a zero surface gravity cannot be achieved (Matter in a zero temperature cannot be reached).

Membrane Paradigm
As mentioned above, in classical general relativity, freely falling observers experience no drama as they cross the event BH horizons, at least not until they reach the singularity inside the BH. However, to a distant and static observer outside a BH, any infalling objects are frozen at the horizons due to the blue-shift effect. Hence, the BH interior can be regarded as an irrelevant region for the static observers. Based on this complementary picture near horizon, in 1986, Kip S. Thorne, Richard H. Price and Douglas A. Macdonald published the idea of membrane paradigm [51]. They use a classically radiating membrane to model the BHs, which is motivated as a useful tool to study physics outside BHs without involving any obscure behavior within BH interior.
Let us introduce a spherical membrane located infinitesimally outside the Schwarzschild radius, a.k.a. the stretched horizon. When the membrane is sufficiently thin, one can use the Israel junction condition to nicely embed the membrane in the Schwarzschild spacetime. The condition is 8 of 84 where f ab is the induced metric of the membrane, K ab is the extrinsic curvatures on its two sides, and T ab is its stress tensor. The infalling observer will cross the horizon and enter the BH interior without possibility of seeing the membrane. However the static observer outside the BH can remove irrelevant interior region from the remaining spacetime with a membrane. Assuming reflection symmetry K + ab = −K − ab , the Israel junction condition on the membrane becomes where the stress tensor T ab is no longer zero but has contribution from the extrinsic curvature on the membrane. Rewriting the left hand side of (7), one can obtain the following relation where σ a b is the shear, θ is the expansion, and κ is the surface gravity at the horizon. From the analogy with the energy momentum tensor of the 2-dimensional compressible fluid, one can read that the shear viscosity η and bulk viscosity ζ are given by η = 1/(16π) and ζ = −1/(16π), respectively. The negativity of the bulk viscosity implies gravitational instability for the expansion or compression at the horizon. The viscosity at the membrane lead to the thermal dissipation of infalling gravitational waves into the horizon and in this sense, the causality of a classical BH is dual to the viscosity of a 2-dimensional fluid at the stretched horizon. This was the earliest example of fluid-gravity correspondence, which is now active area of research.
Modifying Einstein gravity which revises the structure of BHs can provide a modified structure of the thin-shell membrane. For example, by adapting the transport properties of the membrane fluid, we can investigate various models of quantum BHs. As an example, we provide a simple idea [23] which relates the reflectivity of the "horizon" of quantum BHs to viscosity in the context of membrane paradigm. We start by perturbing the Schwarzschild spacetime, whose metric is g Sch µν . Within Regge-Wheeler formalism [52], the axial axisymmetric perturbation g µν = g Sch µν (r) + δg µν (r, θ, t) take teh form: where other δg µν components vanish, and 1 controls the order of perturbation. The membrane stands at r = r 0 + R(t, θ), where r 0 is its unperturbed position. We apply the Israel junction conditions K ab − K f ab = −4πT ab to Brown-York stress tensor as defined in [53]. The indexes µ, ν run over (t, r, θ, φ) in the 4d spacetime, while a, b run over (t, θ, φ) on the 3d membrane. We further assume that T ab is the energy stress tensor of a viscous fluid: where ρ 0 and p 0 (ρ 1 and p 1 ) are background (perturbation on) membrane density and pressure, and u a , η and ζ are fluid velocity, shear viscosity, and bulk viscosity, respectively. Plugging Eqs. (9)(10)(11)(12)(13) into the the Israel junction condition, we find in the zeroth order in : where g(r 0 ) = (1 − 2M/r 0 ) 1/2 and f (r 0 ) = 1 − 2M/r 0 . Assuming u φ = 0, equation of θφ component gives in next order of : We can further use ψ ω = 1 r 1 − 2M r h 1 (r) and the tortoise coordinate x = r + 2M log[r/(2M) − 1] to rewrite Eq. (16) as For the classical BHs with a purely ingoing boundary condition ψ ω ∝ e −iωx at the horizon, Eq. (17) gives η = 1 16π , which is consistent with the standard membrane paradigm. If instead we assume there is no longer horizon but a reflective surface with ψ ω = A out e iωx + A in e −iωx , Eq. (17) gives: Which relates the reflectivity of the membrane to the viscosity of the surface fluid.

Dawn of Gravitational Wave Astronomy
From 2015 onwards, the LIGO/Virgo collaboration reported unprecedented GW observations from binary BH merger events [31][32][33][34][35][36][37][38]. It is the first time that humankind can detect GWs after one century of the Einstein's general theory of gravity. In 2017, Rainer Weiss, Kip Thorne and Barry C. Barish won the Nobel Prize in Physics "for decisive contributions to the LIGO detector and the observation of Gravitational Waves".
The first and most prominent binary BH merger signal seen by LIGO, GW150914, matches well with predictions of numerical relativity simulations that settle into Kerr metric, but contrary to original claims, it could not confirm the existence of the event horizons [18]. However, it opened a new front to test general relativity in strong gravity regime and Kerr-like spacetimes (e.g., Quantum BHs) from modified gravity, which is the main topic of this review article. This is the dawn of GW astronomy, and we stand at the threshold of a new age. We are detecting even more compact binary merger events with a better sensitivity from the O3 run of LIGO/Virgo. Future experiments such as Einstein Telescope, Cosmic Explorer, and LISA are expected to improve this by orders of magnitude. More studies on the echo-emission mechanism as well as observational strategies will be crucial for taking advantage of these new observations, to shed light on the nature of quantum BHs. It is our point of view that the best bet is on a sustained synergy between theory and observation, relying on well-motivated theoretical models (such as the Boltzmann reflectivity, aether holes, 2-2 holes, or fuzzballs, discussed in this review) to provide concrete templates for data analysis, which in turn could be used to pin down the correct theory underlying quantum BHs. With some luck, this has the potential to revolutionize our understanding of fundamental physics and quantum gravity.

Quantum Gravity and Equivalence Principle
The Einstein's general theory of relativity is classical. However, in the Einstein field equation G µν = 8πGT µν , the classical spacetime geometry is related to stress energy tensor of quantum matter. For decades, scientist have tried to reconcile this inconsistency by embedding general relativity (or its generalizations) within some quantum mechanical framework, i.e. quantum gravity.
Conventional approach to quantizing Einstein gravity fails because it is not renormalizable. This implies that making predictions for observables, such as scattering cross-sections, requires knowledge of infinitely many parameters at high energies, leading to loss of predictivity. In the modern language, general relativity could at best be an effective field theory, and requires UV-completion beyond a cutoff near (or below) Planck energy (e.g., [54]).
Most proposals for this UV-completion involve replacing spacetime geometry with a more fundamental degree of freedom, such as strings (string theory) [55], discrete spins (loop quantum gravity) [56], spacetime atoms (causal sets) [57], or tetra-hydra (causal dynamical triangulation) [58]. More exotic possibilities include Asymptotic Safety [59], Quadratic Gravity [60], and Fakeon approach [61] that introduce a non-perturbative or non-traditional quantization schemes for 4d geometry. Yet another possibility is to modify the symmetry structure of General Relativity in the UV, as is proposed in Lorentz-violating (or Horava-Lifshitz) quantum gravity [62].
While proponents of these various proposals (with varying degrees of popularity) have claimed limited success in empirical explanations of some natural phenomena, it should be fair to say that none can objectively pass muster of predicitivity. As such, for now, the greatest successes of these proposals remain in the realm of Mathematics.
Due to this lack of concrete predictivity, the EFT estimates (discussed above) are instead commonly used to argue that the quantum gravitational effects should only show up at Planck scale ∼ 10 −35 meter or 10 28 eV, which is far from anything accessible by current experiments. However, such arguments miss the possibilities of non-perturbative effects (such as phase transitions) which depend on a more comprehensive understanding of the full phase space of the specific quantum gravity proposal.
For example, it has been shown that the non-perturbative quantum gravitational effects may lead to Planck-scale modifications of the classical BH horizons [63]. Proposed models like gravastars [19], fuzzballs [20,[64][65][66][67], aether BHs [17], and firewalls [68,69] amongst others [70][71][72] all drastically alter the standard structure of the BH stretched horizons with a non-classical surface. Soon after the first reported detection of gravitational waves, [18] discerned that Planck-length structure modification around horizons leads to a similar waveform as in classical GR, but followed by later repeating signalechoes -in the ringdown from the reflective surface that replaces the classical horizon. This discovery equals a new road leading to Rome -quantum nature of gravity -and has sparked off a novel area of modeling and searching for signatures from Quantum BHs. The next section will discuss the quantum theories of BH models and possible road maps to probe them, inspired by the detection of binary BH merger events in gravitational waves.

Evaporation of BHs and the Information Paradox
It was already recognized by Stephen Hawking in the 1970s that the evaporation of a BH leads to an apparent breakdown of the unitarity of quantum mechanics. Here, we will briefly review this problem, which is known as the BH information loss paradox [73]. In the context of quantum field theory in curved spacetime, the energy flux out of a BH horizon is obtained by specifying a proper vacuum state and fixing the (classical) background spacetime. However, a radiating BH must lose its mass in time, and so fixing 11 of 84 the background is valid only for a much shorter timescale than the evaporation timescale. One can roughly estimate the lifetime of a BH as follows: The energy expectation value of a Hawking particle is of the order of the Hawking temperature T H ≡ (8πM) −1 , which would be emitted over the timescale of t ∼ M. Then we can estimate the luminosity of the BH as and this gives its lifetime To be consistent with the result of a more rigorous calculation (see e.g. [74]), we need a factor of about 10 5 in (20) which is much much longer than the cosmic age of ∼ 4 × 10 17 [sec] for astrophysical BHs whose mass are M . It may be true that BHs evaporate due to the Hawking radiation, at least, until reaching the Planck mass. However, the gravitational curvature near the horizon eventually reaches the Planckian scale and the classical picture of background gravitational field would break down. As such, the possibility of leaving a "remnant" after the evaporation has been discussed (see e.g. [75][76][77][78]), but the most natural possibility would be that only Hawking radiation is left after the completion of the BH evaporation. If the Hawking evaporation just leaves the "thermal" radiation afterwards, one can immediately understand why the evaporation process is paradoxical. Let us suppose that a pure quantum state collapses into a BH and it radiates Hawking quanta until the BH evaporates. If the final state is a thermal mixed state, the evaporation is a process which transforms a pure to mixed state. Therefore, if the final state of any BH is a completely thermal state, one can say that the evaporation process is a non-unitary process. The information loss paradox can be also explained from the geometric aspect using the Penrose diagram. In quantum mechanics, the time-evolution of a quantum state is described by a unitary operator,Û, that maps an initial quantum state |in on a past Cauchy surface Σ i into a final quantum state |f on a future Cauchy surface Σ f . Since the unitary operator gives a reversible process, one can also obtain the initial state from the final state as |in =Û † |f .
Although this is true in a flat space, the argument is very controversial in the existence of an evaporating BH. Assuming a gravitational collapse forms a horizon and singularity, then it eventually evaporates, leaving behind a thermal radiation, the Penrose diagram describing the whole process is given by Fig. 1 The resulting density matrix, (25), is independent of the interior orthogonal basis |ψ j int due to the tracing operation. Therefore, the loss of the interior information results in a non-unitary evolution and an initial quantum state evolves to a mixed state after the BH evaporation.

BH complementarity
The BH complementarity has been one of the leading proposals for the retrieval of BH information, which was first put forth by by Susskind, Thorlacius, and Uglum [79]. According to a distant observer, due to the infinite redshift at a BH horizon, the Hawking radiation involves modes of transplanckian frequency whose energy can be arbitrarily large in the vicinity of the horizon. In the BH complementarity proposal, the energetic modes form the membrane, which can absorb, thermalize, and reemit information, on the BH horizon. They argue that such a picture regarding the retrieval of BH information by the stretched horizon is consistent with the following three plausible postulates: Postulate 1 (unitarity)-According to a distant observer, the formation of a BH and the evaporation process can be described by the standard quantum theory. There exists a unitary S-matrix which describes a process from infalling matter to outgoing non-thermal radiation.

of 84
Postulate 2 (semi-classical equations)-Outside the stretched horizon of a massive BH, physics can be approximately described by a set of semi-classical field equations. On the other hand, it has been presumed that a freely infalling observer would not observe anything special when passing through the horizon due to the equivalence principle. In this sense, there are two totally different and seemingly inconsistent scenarios that co-exist in the BH complementarity. However, the contradiction arises only when attempting to compare the experiments performed inside and outside horizon, which might be impossible due to a backreaction of the high-energy modes near the stretched horizon [80].

Firewalls
In 2012, Almheiri, Marolf, Polchinski and Sully (AMPS) argued [69] that the Postulates 1-3 in the BH complementarity and the Equivalence principle of GR are mutually inconsistent for an old BH [81][82][83], provided that the monogamy of entanglement is satisfied. Then they argued that the "most conservative" resolution is a violation of the equivalence principle near the BH and its horizon should be replaced by high-energetic quanta, so called "firewall", to avoid the inconsistency. Before introducing the original firewall argument in more detail, let us review a theorem in quantum information theory, the monogamy of entanglement. Let us consider three independent quantum systems, A, B, and C. The strong subadditivity relation of entropy is given by If A and B is fully entangled, we have S AB = 0 and S ABC = S C .
Then the strong subadditivity relation reduces to Since the left hand side in (28) is the mutual information of B and C, denoted by I BC , and it is a non-negative quantity, (28) reduces to which means that the quantum system B cannot fully correlate with C when B and A are fully entangled mutually. Therefore, any quantum system cannot fully entangle with other two quantum systems simultaneously. This is the monogamy of entanglement that is an essential theorem in the firewall argument. Let us consider an old BH, whose origin is a gravitational collapse of a pure state, with early Hawking particles A, late Hawking particle B, and infalling particle inside the horizon C. In order for the final state of the BH to be pure state, A and B should be fully entangled mutually, that is a necessary condition for the Postulate 1. On the other hand, created pair particles , B and C, are also fully entangled according to the quantum field theory in classical background (Postulate 2). That is, imposing the Postulate 1 and 2 inevitably results in that B is fully and simultaneously entangled with both A and C, which obviously 14 of 84 contradicts with the monogamy of entanglement. In order to avoid this contradiction, AMPS argued that there is no interior of BHs and the horizons should be replaced by energetic boundaries that the entanglement of Hawking pairs are broken. They called these boundaries "firewalls". According to this proposal, any object falling into a BH would burn up at the firewall, which contradicts the equivalence principle (in vacuum) and replaces the BH complementarity proposal. Although there are some updates of this proposal, based on ER=EPR conjecture [16,[84][85][86][87], backreaction due to gravitational schockwaves [88], and quantum decohenrence of Hawking pair due to the interior tidal force [89]), they do remain speculative, and at the level of toy models. However, on general grounds, if quantum effects lead to such an energetic wall at the stretched horizon, it could contribute to the reflectivity of BH which may be observable by merger events leading to the formation of BHs.

Gravastars
The gravitational vacuum condensate star (gravastar) was proposed as a final state of gravitational collapse by Mazur and Mottola [19]. According to the proposal, the resulting state of gravitational collapse is a cold compact object whose interior is a de Sitter condensate, which is separated from the outside black hole spacetime by a null surface. In this state, there is no singularity (with the exception of the null boundary) and no event horizon, which avoids the BH information loss paradox. Such gravitational condensation could be caused by quantum backreaction at the Schwarzschild horizon r = r g even for an arbitrarily large-mass collapsing object. One might wonder why the backreaction can lead to such a drastic effect for any mass since the tidal force which acts on an infalling test body can be arbitrarily weak for an arbitrarily large mass at the Schwarzschild radius. The argument is that considering a photon with asymptotic frequency ω near the Schwarzschild radius, the (infinite) blue-shift effect by which the local energy is enhanced ashω/ 1 − r g /r, could lead to a drastic effect at the Schwarzschild radius. This is unavoidable since any object is immersed in quantum vacuum fluctuations and virtual particles always exist around them. From this argument, the gravitational condensation has been expected to take place at the final stage of gravitational collapse. The authors in [19] also estimate the entropy on the surface of gravastar by starting with a simplified vacuum condenstate model which consists of three different equations of state where r 1 is the radius of interior region and δr is the thickness of the thin-shell of the gravastar. Then the obtained entropy of the shell was found out to be S ∼ 10 57 gk B (M/M ) 3/2 , where g is a dimensionless constant. Recently, the derivation of gravastar-like configuration was performed by Carballo-Rubio [90]. He derived the semi-classical Tolman-Oppenheimer-Volkoff (TOV) equation by taking into account the polarization of quantum vacuum and solved it to obtain the exact solution of an equilibrium stellar configuration. It also has its de Sitter interior and thin-shell near the Schwarzschild radius, which is consistent with the original gravastar proposal [19]. From the observational point of view, the shadows of a gravastar was investigated in [91] where they argue the shadows of a BH and gravastar could be distinguishable. In addition, tests of gravastar with GW observations have been discussed in e.g. [4,24,92].

Fuzzballs
Samir D. Mathur has proposed fuzzballs [93] as description of true microstates of the quantum BHs from string theory. A fuzzball state has the BH mass inside a horizon-sized region and a smooth (but higher-dimensional) geometry. Here are some crucial features of the conjecture: 1. Different fuzzball geometries represent different microstates of the quantum BH -fuzzball.
Application the AdS/CFT duality [94] suggests that the counting of the microstates is consistent with the Bekenstein-Hawking entropy. 2. Fuzzballs do not possess horizons. Instead, they end with smooth "caps" near where the horizons would have been. Every microstate has almost the same geometry outside the would-be horizon matching the classical BH picture for the outside observers. But the microstates differ from each other near the would-be horizons. 3. Fuzzball solves the information paradox by removing the horizon and singularity. The horizon is replaced by fuzzy matter and no longer vacuum. The particles created near the would-be horizon now have access to the information of fuzzball interior. Moreover, the higher-dimensional spacetime ends smoothly around the would-be horizon and is singularity-free. The infalling particles at the low frequencies interact with the "fuzz" for a relatively long time scale, while high frequency ones excite the microstates and lose their energy the same as in the classical BHs case. Hence, the traditional horizons only show up effectively from the point of view of an outside observers, over relatively short time scale M log(M).
How do these higher dimensional "microstates" with the smooth and horizonless geometries looks like? We, for the first time, show a specific reduced 4D fuzzball solution has an associated 4D effective fluid near the would-be horizon. The anisotropic pressure of the fluid is crucial to the horizonless geometry.
Applying Kaluza-Klein reduction of non-supersymmetric microstates of the D1-D5-KK system [95]. the metric in 4D is where parameters c 1 , c 5 , s 1 , s 5 , r 0 , n and p are related to the mass, angular momentum and charges of the solution.

Asymptotic behavior
Here, we study the asymptotic behavior of metric. As shown in Table 1, it behaves exactly like Schwarzschild metric when setting K=1, M=r g . They have different g tϕ compared with Kerr metric. As stated, it resembles the Schwarzschild BH far away, but present different geometry close to the horizon.

Matter field
We can now study the effective 4d matter stress tensor from the Einstein tensor of the 4d fuzzball geometry (33). For a sample choice of parameters, outlined in Table 2 where f (r, θ) and g(r, θ) are functions of coordinate r and θ. The concrete expressions depend on parameter setting. ρ, P 1 , P 2 and P 3 are the energy density and anisotropic pressure of the matter field. Their behavior near horizon is shown in Fig. 2. Pressure P 1 = −P 2 is an analytic result and true for any parameter setting, while relationship between the energy and pressure: −ρ = − 5981 1461 f (r, θ) − g(r, θ) = − 5981 1461 P 1 − P 3 is a numerical approximation. The approximation is exact far away from the fuzzball with parameters in Table  2 and for any given θ except 0 and π. The relationship changes near the horizon shown as in Fig 3 after averaging over θ. At around r ∼ 1000, we have radial pressure equals tangential pressure P 1 = P 3 . Similar to the metric, matter fields are singular at r = r 0 , θ = 0 and π.
The matter field has an anisotropic pressure. It is not traceless so it cannot be a simple electromagnetic field. We also checked that it cannot be a single scalar field. However, most fuzzball microstates still remain intractable, with no clear dimensional reduction or 4d geometry. To circumvent this obstacle, we have proposed a "mock fuzzball" spacetime [96] which captures the horizonless feature of the model with an anisotropic fluid. This conjecture leads to an interesting application to dynamic binary quantum BH merger simulations, as discussed later in Sec. 6.2.

Mock Fuzzballs
Here, we introduce a "mock" fuzzball geometry, based on the motivation to build a generic and macroscopic metric which captures important, coarse-grained properties of the fuzzball, i.e, the metric has neither horizon nor singularity and spacetime ends around the stretched horizon. Fig. 4 shows how the causal diagram of a Schwarzschild BH is changed to remove the horizon. We study different mock fuzzballs, and check the corresponding matter fields, a swell as potential observable effect. The simplest case of the mock fuzzball 1 is given by: where parameter 0 < b 1 is only important around 2M. It ensures that g tt doesn't vanish and thus remove the horizon. The geometry resembles a traditional BH far away. Besides, this metric is only valid where r > 2M to imitate a fuzzball metric which ends around the stretched horizon. The corresponding matter field is other components vanish (47) The anisotropic behavior of the pressure near the stretched horizon is shown as in Fig. 5 with b = 0.01 and M = 1 2 . The absolute value of two pressures are equal at r = 3M 1+b . Tangential pressure is much larger than 1 It turns out that this toy model spacetime coincides with the one proposed earlier by [97].

Figure 3.
Near-horizon relationship between density and pressure of the fuzzball solution from paper [95]. The dashed line is the asymptotic behavior, and the solid line is the real behavior of P 1 P 3 −ρ radial pressure near horizon, and opposite far away. It doesn't have any singularity at r=2M since all fields reach an extremum there.

Other mock fuzzballs
Besides the simplest case introduced in the last section, we can modify other terms in Schwarzschild metric to recover energy density which fuzzballs in Sec. 3.5 actually have. We also study charged and rotating BHs.

• Schwarzschild metric
The simplest mock fuzzball studied above has no energy density. However, we can recover energy density by changing g θθ and g ϕϕ within spherical symmetry: other components vanish (52) Small b and d ensures that ρ P r P t . In addition, b > 0, d < 1 ensure finite Ricci scalar without curvature singularity at r > 2M. Another possible modification to recover energy density is to assume that mass has a time dependence: • Extremal BH metric 21 of 84 For an extremal BH, the fuzzball has another interesting property: Proper length from somewhere near stretched horizon to "horizon" is finite, in contrast to the infinite throat in the traditional picture. Parameter c here captures the finite throat.
Assuming Einstein field equations, mock fuzzball geometries can only be sourced by matter fields with exotic (and anisotropic) equations of state. Considering simplest Schwarzschild mock fuzzball: There is no longer vacuum outside the horizon, which may potentially lead to observable effects, depending on how strongly the "fuzz" matter can interact with the detectors. To visualize the signal, we assume a geodesic observer radially falling towards the stretched horizon with zero velocity at infinity. We calculate (see Appendix A for details) two observable scalars: energy density U and energy flux F , as seen by the observer: where T µν is from Einstein field equation of the mock fuzzball, a ν is the unit detector area vector and u µ is the four-velocity of the observer with a µ u µ = 0. Both energy density and flux are finite and vanish when parameter b vanishes.

of 84
To visualize observable energy density and flux, we choose parameters M = 1 and b = 0.1 and compare it with the signal of Hawking radiation as Fig. 6 and Fig. 7. Both energy density and flux of the mock fuzzball are larger than those of Hawking radiation near the horizon, but drops faster away from the horizon. Energy density of the mock fuzzball is negative while Hawking radiation is positive. Flux of the mock fuzzball flows into BH while Hawking radiation flows into BH near horizon and changes direction away from the horizon.
Specially, energy density and flux of fuzzball are finite near horizon while Hawking radiation diverges. This is because we use Page's approximation of the Hawking radiation [98], which assumes the Hartle-Hawking state for expectation value of stress-energy tensor: The Hartle-Hawking state is a thermal equilibrium states of particles while the Hawking radiation is not in equilibrium. The inconsistency leads to infinite flux and energy density of Hawking radiation. More precise correction can be found in [99].

Aether Holes and Dark Energy
In 2009, Prescod-Weinstein, Afshordi, and Balogh [100] studied the spherically symmetric solutions of the Gravitational Aether proposal for solving the old cosmological constant problem [101,102]. Surprisingly, they showed that if one sets Planck-scale boundary conditions for aether near the horizons of stellar mass BHs, its pressure will match the observed pressure of dark energy at infinity.
In the Gravitational Aether proposal [101,102], the modified Einstein field equation is given by where G = 4 3 G N , and then energy-momentum tensor of aether is assumed to be a perfect fluid with stress-energy tensor T µν without energy density. Here, quantum vacuum energy decouples from the gravity, as only the traceless part of the matter energy-momentum tensor appears on the right-hand side of the field equations. It can be shown that the Bianchi identity and energy-momentum conservation completely fix the dynamics, and thus the theory has no additional free parameters, or dynamical degrees of freedom, compared to General Relativity.
The modified Schwarzschild metric is the vacuum solution with spherical symmetry in modified equations, and identical to a traditional equations sourced by the aether perfect fluid. Far away from the would-be horizon but close enough to the origin (2M r |p 0 | −1/2 ), the solution has the form Figure 6. The energy density of mock fuzzball, compared with that of Hawking radiation as seen by a radially infalling observer along a geodesic. Here the orange curve is positive and the black curve is negative. The energy density of mock fuzzball is larger than the Hawking radiation near horizon and drops faster than the Hawking radiation away from the horizon.
which can be compared to the de Sitter metric We see that assuming p 0 = − 2 3 ρ Λ , the g tt 's agree with each other. Therefore, the Newtonian observers (for 2M r |p 0 | −1/2 ) will experience the same acceleration as in the de-Sitter metric with the cosmological constant. However, on larger scales, one has to take into account the effects of multiple black holes and other matter in the Universe. The Planckian boundary conditions at the (would-be) horizon relates the pressure of the aether to the mass of the astrophysical BHs, −p 0 ∼ M −3 [100]. In particular, the BH masses within the range 10 M − 100 M , which correspond to the most astrophysical BHs in galaxies, yield aether pressures comparable to the pressure of Dark Energy, inferred from cosmic acceleration. Moreover, Ricci scalar is inversely proportional to g tt , so the event horizon where g tt = 0 has a curvature singularity, which is reminiscent of the firewall and fuzzball proposals discussed above.
In particular, the fuzzball paradigm is a good approach to remove the singularity. On the one hand, fuzzball gives an extra anisotropic matter field similar to the aether theory, which stands as a good evidence that quantum effects can modify the Einstein field equation with extra sources of 4d energy-momentum like aether. Furthermore, fuzzball is a regular and horizonless geometry, which might indicate the singularity is removable in the full quantum picture of BHs.

2-2 holes
In general relativity, gravitational collapse of ordinary matter will always leads to singularities behind trapping horizons [103]. In [21], Holdom and Ren revisited this problem with the asymptotically free quadratic gravity, which could be regarded as a UV completion of general relativity [21]. The quantum quadratic gravity (QQG), whose action is given by is famously known to be not only asymptotically free, but also perturbatively renormalizable [104][105][106][107]. However, it suffers from a spin-2 ghost due to the higher derivative terms, which is commonly regarded as a pathology of the theory. In [21], it is proposed that the ghost may not be problematic when M is sufficiently small, so that the poles in the perturbative propagators fall into the non-perturbative regime, and the perturbative analysis of ghosts is not reliable. Then it is conjectured [60] that the full graviton propagator in the IR, when M Λ QQG , the spin-2 ghost pole is absent in an analogy with the quantum chromodynamics (QCD) where the gluon propagator, describing off-shell gluons, also does not have a pole. Here Λ QQG is a certain critical value in QQG, analogous to confinment scale Λ QCD in QCD. Based 25 of 84 on this conjecture, the asymptotically free quadratic action in (81) may involve small quadratic corrections at super-Planckian scale, and so the super-Planckian gravity might be governed by the classical action Since gravitational collapse would involve the super-Planckian energy scale, applying the classical action (82) to such a situation is interesting from a point of view of the quantum gravitational phenomenology.
Then the authors in [21] found a solution of horizonless compact object, so-called 2-2 hole, in the classical quadratic gravity. 2-2 holes have an interior with a shrinking volume and a timelike curvature singularity at the origin. It also has a thin-shell configuration, leading to non-zero reflectivity at the would-be horizon, which may cause the emission of GW echoes [4]. Recently, 2-2 holes sourced by thermal gases were also investigated in [108,109].

Non-violent Unitarization
A separate class of possible approaches to the BH information paradox involves a violation Postulate 2 in BH complementarity, i.e. non-locality of field equations well outside the stretched horizon, which is dubbed as "nonviolent unitarization" by Steve Giddings [110]. Such a possibility would allow for transfer of information outside horizon around the Page time (e.g., [111,112]), but could also lead to large scale observable deviations from general relativistic predictions in GW and electromagnetic signals [113]. However, it is not clear whether this non-locality is only limited to BH neighborhoods, and if not, how it could affect precision experimental/observational tests in other contexts. Moreover, in contrast to GW echoes that we shall discuss next, it is hard to provide concrete predictions for astrophysical observations in the nonviolent unitarization scenarios.

Gravitational Wave Echoes: Predictions
GW echoes may be one of the observable astrophysical signals, a smoking gun, so to speak, for the quantum gravitational processes near BH horizons. A number of models of Exotic Compact Objects (ECOs) that we discussed above are expected to emit GW echoes. Some examples are wormholes [18], gravastars [24], and 2-2 holes [21]. Moreover, even Planckian correction in the dispersion relation of gravitational field [22,23,29] and the BH area quantization [114] may also lead to echo signals. Not only the specific models to reproduce GW echoes but also comprehensive modeling of echo spectra in non-spinning case [115], in spinning case [29,116], and in a semi-analytical way [117,118] have been investigated, which enable us to easily obtain echo spectra. In this section, we review the details of GW echoes by starting with the Chandrasekhar-Detweiler (CD) equation [119,120] that is a wave equation with a purely real angular momentum barrier in the Kerr spacetime. We also provide a short review of the GW ringdown signal, that is followed by the GW echo, and the superradiance of spinning BHs. The superradiance with a high reflectivity at the would-be horizon may cause the ergoregion instability, which we shall also discuss separately.

On the equations governing the gravitational perturbation of spinning BHs
The GW ringdown is one of the most important signals to probe the structure of BH since it mainly consists of discrete QNMs of BH characterized by mass and spin. In this subsection, we review that the QNMs can be obtained by looking for specific complex frequencies such that the mode functions of GWs 26 of 84 satisfy the outgoing boundary condition. Let us start with the CD equation [119,120] that is the wave equation for a spin-s field and has a purely real angular momentum barrier: where T is the source term and the potential V ij with i, j = ±1 is given by The functions in (84) are defined by Equation (84) gives four potentials, (i, j) = (−1, +1), (+1, −1), (+1, +1), and (−1, −1). One has to use the different potentials in order to cover the whole frequency space with the CD potentials because The CD equation is obtained as the generalized Darboux transformation of the Teukolsky equation [121] 2 . In the asymptotic regions, r * → ±∞, the CD equation reduces to the following wave equation whereω ≡ ω − mΩ H , in terms horizon angular frequency Ω H ≡ a/(2Mr + ), and horizon outre radius In the following, we will omit the subscripts of l, m, and s for brevity. One can read that the homogeneous solutions of (94) are given by the superposition of ingoing and outgoing modes where A, B, C, and D are arbitrary constants. The QNMs can be found by looking for the complex frequencies at which the homogeneous solution satisfies the outgoing boundary condition of A = D = 0. This is equivalent to looking for the zero-points of the Wronskian between the two homogeneous solutions where the two homogeneous solutions satisfy the following boundary conditions Recently, the QNMs of Kerr spacetime were precicely investigated in [122] by using the method developed by Mano, Suzuki, and Takasugi [123][124][125][126] that enables us to obtain the solution of the Teukolsky equation in an analytic way.

Transmission and reflection coefficients of the angular momentum barrier
The GW echoes are results of multiple reflections in the cavity between the would-be horizon (e.g., fuzzball/firewall) and angular momentum barrier (see Fig. 9) and so the amplitude of echoes is mainly determined by the reflectivities of the would-be horizon and angular momentum barrier. In this subsection, we review the calculation of the reflectivity of angular momentum barrier.
From the mode functions (97,98), one can obtain the energy conservation law for the incident, reflected, and transmitted waves by using another Wronskian relatioñ  which is constant for real frequency due to the reality of the angular momentum barrier in the CD equation.
Then we obtain the following relations by usingW BH From the above relations, one can read that the energy reflectivity and transmissivity for inward incident waves respectively, and those for outward incident waves are given by From (102,103), one can calculate the reflectivity/transmissivity of the angular momentum barrier by numerically solving the homogeneous CD equation (83). In the spinning caseā > 0, the energy reflectivity is greater than 1 for −mΩ H <ω < 0, a phenomenon that is often referred to as BH superradiance. The superradiance can be characterized by the amplification factor Z ≡ (I ←/→ ref ) 2 − 1, and when only interested in the low frequency region, one can use the analytic expression [127] Z 4Qβ sl where r − is the radius of inner horizon, β sl ≡ (l−s)!(l+s)! (2l)!(2l+1)!! and Q ≡ −  10. As can be seen from FIG. 10, the energy flux reflectivity exceeds 1, which means that the energy of a spinning BH is extracted by reflected radiation, within the superradiance regime. We will discuss the ergoregion instability caused by the superradiance and the reflectivity of would-be horizon in Sec. 4.4.

Transfer function of echo spectra and geometric optics approximation
When the GW echo is caused by an incident wave packet repeatedly reflected between the cavity, one can use the geometric optics approximation to predict the GW echo signal, which was first pioneered in [115]. Let us start with the calculation of the Green's function of GW ringdwon signal G BH by using the CD equation. It satisfies where we omit the subscripts of V ij . Once imposing the outgoing boundary condition, the Green's function G BH is uniquely determined as where r * < ≡ min(r * , r * ) and r * > ≡ max(r * , r * ). Therefore, when there is no reflectivity at the horizon, the Fourier mode of GWs at infinity and at the horizon can be obtained as where T(r * ) is the source for the inhomogeneous CD equation. If there is no reflection near the horizon, the relevant observable spectrum is only Z ∞ , and Z BH is irrelevant for observation. On the other hand, if reflection at the would-be horizon is caused by a certain mechanism, Z BH is also observable in addition to Z ∞ .
One can obtain echo spectra by using the geometric optics approximation, which should be reliable as long as the would-be horizon and angular momentum barrier are well separated in tortoise coordinates, r * . The amplitude of the first echo, Z (1) echo , can be estimated by and the second echo may have the amplitude of As such, one can obtain the amplitude of n-th echo as where R → BH ≡ B in /B out and T → BH ≡ ω/|ω|B −1 out . Since only the reflectivity and transmissivity of outgoing waves are involved in the echoes, we will not use R ← BH and T ← BH in the following, and so omit the symbol →. Summing up all contributions from n = 1 to n = ∞, one obtains

of 84
Note that here we assume that |RR BH | < 1, as otherwise the infinite sum of the geometric series does not converge. Finally, the spectrum taking into account the reflection at the would-be horizon is obtained When R = 0 we have K = 0 and so it reduces to (107). Once we specify a specific form of the source term T, one can obtain Z BH /Z ∞ . For example, let us assume the source term located at where S(ω) is a non-singular function in terms of frequency. Substituting this source term in (107) and (108), one obtains Note that this is independent of the function S(ω). Therefore, we finally obtain the following transfer function As discussed in [128], actually K + echo and K − echo represent two different trajectories of GWs in the cavity. Here we are interested in outgoing incident waves that is related to K + echo and so in the following we discard K − echo from the transfer function, which does not change the qualitative feature of resulting echo signals.
Once we determine the spectrum of injected GWs, Z ∞ , one can obtain a template of GW echoes. Using the spectrum of GW ringdown may be a good approximation to obtain a realistic template. In this case, Z ∞ is given by [129] where D o is the distance between the GW source and observer, ω lm0 is the most long-lived QNM,Ã lm0 is the initial ringdown amplitude, φ lm0 is the phase of ringdown GWs, and θ is the observation angle.
The amplitudeÃ lm0 is proportional to √ rd , where rd ≡ E GW /M and E GW is the total energy of GW ringdown [129]. To give a few examples, the ringdown + echo spectra are shown in FIG. 13.

Ergoregion Instability and the QNMs of quantum BH
As pointed out in the previous subsection, one should check if |RR BH | < 1 is satisfied when calculating the transfer function in the geometric optics picture. This is physically important to understand the ergoregion instability caused by the reflection at the would-be horizon. Since the common ratio of the geometric series is RR BH , the echo amplitude may be amplified and diverges when |RR BH | > 1. This is nothing but the ergoregion instability that prevents BHs from having high spins. One can also derive 32 of 84 Figure 11. QNMs with |R| = 0.7 (red filled circles) and 0.9 (blue filled circles), r * 0 = −40M,ā = 0.9, and = m = 2.
the criterion from the QNMs of ECOs. The echo QNMs ω n can be obtain by looking for the poles of the Green's function of echo GWs. That is, one can look for the poles from the zero points of the denominator of the transfer function and we obtain where ∆t echo ≡ 2|r * 0 |, δ ≡ arg[R] and δ ≡ arg[R BH ]. Then we obtain the real and imaginary parts of the echo QNMs The positivity of the imaginary part of QNMs, which leads to the instability, is equivalent to having |RR BH | > 1. Furthermore, we can see that the real parts of the QNM frequencies depend on the phases of R and R BH , while their imaginary part depends on their absolute values. We can also rewrite the imaginary part in terms of the amplification factor Then we obtain the analytic form of the imaginary part of QNMs in the low-frequency regime (Mω 1) where we used (104). This is the generalization of the analytic form of QNMs [128]. This analytic form is well consistent with numerically obtained QNMs in the low frequency region as is shown in FIG. 11. 33 of 84

Echoes from Planckian correction to dispersion relation and Boltzmann Reflectivity
The Boltzmann reflection of a BH horizon has been discussed in the context of (stimulated) Hawking radiation from the path integral approach [48], quantum tunneling approach [130,131], and Feynman propagator approach [132]. Recently, two of us studied the reflectivity of a BH for incident GWs from a Lorentz violating dispersion relation and argued that it can be approximated by a Boltzmann-like reflectivity [22]. More recently, three of us used general arguments from thermodynamic detailed balance, fluctuation-dissipation theorem, and CP-symmetry to show that the reflectivity of quantum BH horizons should be universally given by a Boltzmann factor [23,29]: The reflection of quantum BH might be understood as Hawking radiation stimulated by enormous number of incoming gravitons, and if that is so, having the dependence of the reflectivity on the Hawking temperature T H is natural. Furthermore, one can also avoid the ergoregion instability in this model [23,29]. In this subsection, we briefly review the Boltzmann reflectivity model from both theoretical and phenomenological aspects.

Boltzmann reflectivity from dissipation
The dissipative effects at the apparent horizon have been discussed from the point of view of the membrane paradigm [51,53], the fluctuating geometry around a BH [133,134], and the minimal length uncertainty principle [135]. Our approach to derive the Boltzmann reflectivity starts with a heuristic assumption to model the dissipative effects, which are expected in any thermodynamic system from fluctuation-dissipation theorem. Let us assume that the wave equation governing the perturbation of BH is given by [23]: whereγ is a dimensionless dissipation parameter, Ω(x) ≡ |ω|/ |g 00 (x)| is the blueshifted (or proper) frequency, and V is the angular momentum barrier. The form of the dissipation term is expected from the fluctuation-dissipation theorem near the horizon, where the Hawking radiation (quantum fluctuation/dissipation) and the incoming GWs (stimulation) are blue shifted. This dissipative modification to the dispersion relation becomes dominant only when the blueshift effect is so intense that the proper frequency is comparable to the Planck energy, Ω ∼ E Pl . Furthermore, from a phenomenological point of view, the dissipative term in (128) is similar to the viscous correction to sound wave propagation in terms of shear viscosity, ν, in Navier-Stokes equation, −i(4/3)νΩ∇ 2 (e.g., [136]). Let us solve the modified wave equation by imposing a physically reasonable boundary condition (see FIG. 12): The constant boundary condition in the limit of r * → −∞ means that the energy flux carried by the ingoing GWs cannot go through the horizon, and is either absorbed or reflected. That is consistent with the BH complementarity [79] or the membrane paradigm [51,137]. Although there is no unique choice of wave where δ(r) ≡ 1 − r g /r + (a/r) 2 is the blue shift factor in terms of the co-rotating frame [30,138]. In the near horizon limit (r * → −∞, see below for details), the CD equation reduces to the following form in the limit of r * → −∞: where κ + is the surface acceleration at the outer horizon, Q is defined as and r ± ≡ M(1 ± √ 1 −ā 2 ). The solution of (131) which satisfies the aforementioned boundary condition is and one can read that in the intermediate region, −κ −1 + log [QE Pl /(γ|ω|)] r * ±κ −1 + , ψω can be expressed as the superposition of outgoing and ingoing modes ψω = e πω/(2κ + ) A + e −iωr * + e −πω/(2κ + ) A * + e iωr * forω > 0, e −πω/(2κ + ) A − e −iωr * + e πω/(2κ + ) A * − e iωr * forω < 0, where A ± has the form of and Therefore, the energy reflectivity is given by and finally we obtain where δ wall is the phase shift at the would-be horizon and it is determined by A + or A − . Equation (138) then reproduces the Boltzmann energy flux reflectivity in (127). As we noted earlier, the same result can be independently derived using thermodynamic detailed balance or CP symmetry near BH horizons. When we further modify the dispersion relation by adding a quartic correction term where C d is a constant parameter andΩ andK are the proper frequency and proper wavenumber, respectively, the exponent of the Boltzmann factor is modified, and the analytic form can be obtain for C d γ by using the WKB approximation [22] |R| exp One of the essential differences between the modified dispersion relation model, that could give the Boltzmann reflectivity at the would-be horizon, and the Exotic Compact Object (ECO) model is the reflection radius r * 0 . In the former case, r * 0 depends on the frequency of incoming GWs,ω, and so the reflection surface is not uniquely determined. This is because the reflection takes place when the frequency |ω| reaches the Planckian frequency at which the modification in the dispersion relation becomes dominant. Therefore, the reflection radius depends on the initial (asymptotic) frequency of incoming GWs (see Equation 136). On the other hand, in the ECO scenario, the reflection radius would be fixed and it would stand at ∼ a Planck proper length outside the horizon. In this case, the reflection radius is given by which depends only on the mass of BH. For a detailed discussion of how one can observationally distinguish these scenarios, we refer the reader to [128].

Phenomenology of Boltzmann reflectivity
Here we summarize some interesting phenomenological aspects of the Boltzmann reflectivity model, and refer the reader to [29] for more details.
Nearly all the previous studies of GW echoes assume a constant reflectivity model, in which the ratio of outgoing to ingoing flux at the horizon is assumed to be independent of frequency. In contrast, the echo spectrum in the Boltzmann reflectivity model can be significantly different. This difference can be seen in the sample echo spectra for both models, shown in FIG. 13. As can be seen from the spectra, the echo amplitude is highly excited near m× the horizon frequency, since the Boltzmann reflectivity is sharply peaked around ω mΩ H ± T H and is exponentially suppressed outside this range. In the extremal limit a → 1, the Hawking temperature becomes zero and so the frequency range in which |R| ∼ 1 vanishes. Therefore, the peaks in echo spectrum is highly suppressed for a highly spinning BH (see FIG. 14) 3 .
This nature of the Boltzmann reflectivity suppresses the ergoregion instability at least up to the Thorne limitā ≤ 0.998. In [128], a more general case is investigated, where the Hawking temperature in the Boltzmann factor is replaced by the quantum horizon temperature, T H → T QH (e.g., as in Equation 140 above) and the ratio T H /T QH is constrained from the ergoregion instability by using |RR BH | < 1. The constraint is T H /T QH 0.5 up to the Thorne limit [128] and so the Boltzmann reflectivity (T H /T QH = 1) is safe up toā 0.998. As an example, we show a time domain function of ringdown and echo phases with T H /T QH = 0.6 in FIG. 15 by implementing the inverse Fourier transform of X = Z ∞ (1 + K + echo ), where we choose Z ∞ so that it reproduces the ringdown phase [128].
Other notable phenomenological properties of quantum BHs with Boltzmann echoes are [29]: • The QNMs of the quantum BH are approximately those of a cavity with a complex length |r * 0 | + i(4T QH ) −1 . • Forγ ∼ 1 (i.e. Planck-scale modifications), the first ∼ 20 echo amplitudes decay as inverse time 1/t, and then exponentially. • Each QNM of the classical BH can be written as a superposition QNMs of the quantum BH for t < ∆t echo . The superposition can be approximated as a geometric series, leading to a closed-form expression for echo waveforms. In particular, the first 20 echoes have approximate temporal Lorentzian envelopes around their peaks, whose width grows linearly width echo number.
The parameter dependence of echo spectrum in the Boltzmann reflectivity model and its consistency with the tentative detection of echo in GW170817 are also investigated in [128] in more detail.

Gravitational Wave Echoes: Observations
For the first time in modern science history we are able to probe the smallest possible theoretical scales or highest possible theoretical energies through GW echoes.
The direct observation of GWs [33] was a scientific breakthrough that has opened a vast new frontier in astronomy, providing us with possible tests of General relativity in the extreme physical conditions near the BH horizons. Motivated by the resolutions of BH information paradox that propose alternatives to BH horizons (see Section 3 above), several groups have searched the LIGO/Virgo public data for GW echoes [18,24] (see Section 4 and [25] for a review), which has led to claims (and counter-claims) of tentative evidence and/or detection [1][2][3][4][5][6][7][8]. While the origins of these tentative signals remain controversial [6][7][8][26][27][28] they motivate further investigation using improved statistical and theoretical tools, and well as new observations.
In astrophysics, GW echoes from quantum BHs can be seen as a transient signal, coming from the post-coalescence phase of the binary BH merger (Fig. 9) or formation of a BH (e.g., via collapse of a hypermassive neutron star; Fig. 16). This section will summarize the current status of observational searches for echoes [139].
In order to properly model echoes, we need a full knowledge of quantum BH nonlinear dynamics, which is so far nonexistent. Therefore, any strategy to search for echoes requires parametrizing one's ignorance, which has so far taken many shapes and form. Indeed, we need to keep a balance between having a simple tractable model (which may simply miss the real signal), or an exhaustive complex model (which may dilute a weak signal with look-elsewhere effects). Current search methods can be generally split into two: Parametrized template-based methods [1,2,6,7,9,140], and "model-agnostic" coherent methods [3][4][5]8].

Echoes from the Abyss: Echoes from binary BH mergers O1 by Abedi, Dykaar, and Afshordi (ADA) [1]
The first search for echoes from Planck-scale modifications of general relativity near BH event horizons using the public data release by the Advanced LIGO GW observatory was developed by Abedi, Dykaar, and Afshordi (ADA) [1]. In this search, a naive phenomenological template for echoes was introduced, leading to tentative evidence at false detection probability of 1% (or 2.5σ significance level 4 shown in Figs. 17, 18 and 19) for the presence of echoes [1]. This work was followed by comments, discussion, and controversy about the origin of this signal [6][7][8][26][27][28]. The ADA model was also later tested for LIGO/Virgo O2 independent events [2], which interestingly, yielded a similar percent-level p-value as O1 (see Section 5.1.3 below). The ADA search was the first phenomenological time-domain echo template search applied to real GW observations [1]. Using a standard GR inspiral-merger-ringdown template M(t), a naive model including five free parameters was proposed: 4 In [1], 2-tailed gaussian probability assigned to significance, e.g., p-value = 68% and 95% correspond to 1σ and 2σ respectively. For a Kerr BH with final mass M BH and dimensionless spin parameter a, the time delay of echoes from Planck-scale modifications of general relativity is: [1,3]: For the final BH (redshifted) masses and spins reported by the LIGO collaboration for each merger event, echo time delays ∆t echo and their errors constrained inside 1σ error are as follows [1]:

Search
In this analysis, ADA devised an echo waveform using theoretical best-fit waveform of Hanford M H,I (t) and Livingston M L,I (t) detectors (in real time series) for the BBH events, provided by the LIGO and Virgo collaborations. The search used the observed data release for the two detectors, h H,I (t) and h L,I (t) respectively, at 4096 Hz and for 32 sec duration. The devised phenomenological echo waveform which was then constructed using five free parameters: 1. ∆t echo : Time-interval between successive echoes, within their 1σ range (Eq. 144).   Figure 18. Same as Fig. 17, but over an extended range of x = t echo −t merger ∆t echo . The SNR peaks at the predicted value of 1 − 0.01 < x < 1 + 0.01 within gray rectangle have false detection probability of 0.11 (0.011) and significance of 1.6σ (2.5σ), for GW150914 (combined events) [1] (See also [27]  2. t echo : Time of arrival of the first echo, which is related to ∆t echo with corrections ∼ ±O(1%) × ∆t echo due to the non-linear dynamics of the merger. 3. t 0 : Truncation time for GR template with a smooth cut-off function, where ω I (t) is frequency of GR template as a function of time [142] and t merger is the time of maximum amplitude of the template. It is assumed that t 0 vary within the range t 0 ∈ (−0.1, 0)∆t echo .
Having this definition, a truncated template is introduced: 4. γ: Damping factor of successive echoes, varying between 0.1 and 0.9. 5. A: Over-all amplitude of the echo waveform with respect to the main event. This free parameter is fitted assuming a flat prior.
The search model of echoes having all the free parameters, assuming a (−1) n+1 factor due to the phase flip at each reflection, is given by: Fig. (20) shows this template using the best fit parameters within the range given above along with the main merger event GW150914.
Once this analysis has been completed for GW150914 (loudest event of O1), it has been repeated for rest of the events combined via SNR maximization:  Table 3. Values of best fit echo parameters of the model Eq. 147 of the highest SNR peak near the predicted ∆t echo (gray rectangle in Fig. 17 The proposed combination takes same γ and t 0 /∆t echo for all events, keeping ∆t echo and A's as free. The results are shown in Fig's (17-19) and Tables 3-4.
• Energy estimation [1]: Given a best-fit template for the echoes, one can provide an estimate for their total GW energy: Uchikata et al. [2] have examined GW echo signals for nine binary BH merger events observed by Advanced LIGO and Virgo during the first and second observation runs (O1 and O2 respectively). They have used several models for a number of searches leading to positive and negative results. In this part we bring their positive results and discuss the rest in part 5.3.3. In this search the critical p-value as 0.05 corrsponding to 2σ significance, (p-value below/above this value) indicates echo signals (are likely/unlikely) to be present in the data.
SNR is evaluated using a matched filter analysis defined using wherex(f) andh(f) are observed data and template in frequency domain respectively, and S n (f) is the noise power spectrum of detector. In this analysis they assume frequency band of f max = 2048 Hz, f min = 40 44 of 84 ent in the data. Our results events and the combined pical value; that is, echo sigamework do not exist in the the signals are too small to ent detector sensitivity. We ion of t 0 weakly affects SNR; t echo is a reasonable assumpcosts. onsider the best fit of the iniθ ini , which is different from ], so it might be inappropriirectly. However, we also ane same template as in Abedi h the same condition for the mparison to those given by hown in Appendix A 1. We 2 events with this template, as that of O1 events. Results . behavior of SNR in Fig. 2 for eters of GW150914 (C01) as and dotted lines correspond d and Livingston), Hanford, y. We can see a peak for the cases near T ∼ 1; however, se is located slightly outside he figure shows that ρ 2 osmpared to Fig. 7 in Ref. [16] st fit initial phase of the temng a frequency-independent s as a parameter, we also ane phase inversion is considlts are given in Appendix B. wer if the phase shift is not .

AND DISCUSSIONS
vitational wave echo signals merger events observed by during the first and second me that the spacetime is ent that a reflective membrane   is located near the event horizon radius. We use the template waveform given by Nakano et al. [23], in which the reflection rate and the phase shift at the potential barrier due to the angular momentum are calculated from Teukolsky equations. We assume a perfect reflection at the membrane; however, the phase shift at the membrane due to reflection is model dependent, so we assume the frequency-independent phase shift at both the membrane and the potential barrier as a parameter. The transmission rate given from the reflection rate strongly suppresses the lower frequencies contained in the template waveform. In addition to the echo parameters, we maximized the signal-to-noise ratio against the initial phase of the template. We used adjacent 4096-second data from open LIGO data for the background estimation, and evaluated the significance by p-values. We found no significant echo signals within our analysis. Since the method of analysis is slightly different from  Hz and normalization condition of (h|h) = 1 [2]. They then perform an echo search by maximizing SNR, following [1]. Fig. 21 presents SNR 2 for the best fit parameters of the event GW150914 (C01) 5 with additional consideration of the best fit initial phase 6 to the ADA [1] template.
In this part, we show the results using the same template given by ADA [1] except the cut-off parameter t 0 as described in 5.1.2 has been fixed to its best fit value, and set the search region of ∆t techo to its 90% (rather than 68%) credible regions in (a, M) space. Similar to ADA, the initial phase of the template is also fixed to zero. Since variation of t 0 weakly affects SNR and has an advantage in saving computational costs, they fixed t 0 = −0.1∆t echo .
Here are the results:  6 Additional consideration of the best fit initial phase has given negative results and is discussed in negative result part below 5.3.3. Therefore, this plot only justifies Uchikata et al. [2] reevaluation of ADA search [1] in Fig. 17 Analysis of Uchikata et al. [2] show that the six independent BBH O2 events in Table 6 have similarly small p-values for ADA echoes as O1. As shown in this table, the total p-value for the six O2 events is 0.039. Combining O2 with O1 events shown in Table 5, leads to the total p-value of 0.047.

Echoes from the Abyss: Binary neutron star merger GW170817 [3]
A binary neutron star merger event collapsing into a black hole (Fig. 16) can also enable us to test general relativity through GW echoes. Although, the current LIGO/Virgo/KAGRA detector sensitivity is blind to post-merger ringdown frequency of GWs, they can be sensitive to low frequency echoes harmonics where the ringdown frequency is suppressed by ln(M BH /M planck ) in Eq. 143 [24]. In other words, since the final mass of BNS merger (2-3 M ) is much smaller than that of the binary BH mergers [1], the lowest harmonics n/∆t echo ( n × 80 Hz) of echo chamber are shifted to the regime of LIGO sensitivity, for small n. Therefore, as first suggested by [4], an optimal model-agnostic search strategy could consist of looking for periodically spaced-harmonics in the frequency space.
Using this model-agnostic search applied to the cross-power spectrum of the two LIGO detectors, Abedi and Afshordi [3] found a tentative detection of echoes around 1.0 sec after the BNS merger, at f echo 72 Hz (see Figs. 22, 23, 24, and 25) using GW event data GW170817 provided by LIGO/Virgo collaboration [37,144]. As it is shown in Figs. 22 and 23, the main signal is also accompanied by secondary lower significance resonances at 73 Hz and t − t merger = 32.9 sec. It is worth noting that after this detection, Gill et al. [145] used independent Astrophysical considerations, based GW170817 electromagnetic follow-ups, to determine that the remnant of GW170817 must have collapsed into a BH after t coll = 0.98 +0. 31 −0.26 sec, which coincides with the detected GW echo signal at 1.0 second (see Fig. 25). This fining of Abedi and Afshordi is consistent with a 2.6 − 2.7 M "BH" remnant with dimensionless spin 0.84 − 0.87. For this signal, considering all the "look-elsewhere" effects, a significance of 4.2σ 7 (see Fig. 26), or a false alarm probability of 1.6 × 10 −5 has been reported, i.e. a similar cross-correlation within the expected frequency/time window after the merger cannot be found more than 4 times in 3 days of GW data. Total energy of detected GW echoes signal using simple assumptions is around ∼ 10 −2 M c 2 .

MODEL-AGNOSTIC SEARCH FOR ECHOES
In this part, we describe the method that leads to the detection of [3], in some detail. The final mass of GW170817 is within ∼ 2 − 3M which could form either a BH or a neutron star (NS). A BNS merger can end up in four possible ways: [144]: ; combining all harmonics with frequencies n × f , with n ∈ N) around the merger for the BNS gravitational-wave event GW170817, observed through cross-correlating the two LIGO detectors [3]. The possible resonance peaks of echoes found in this plot are marked with a green squares. The color scale shows the peak at f peak = 72 (±0.5) Hz and t − t merger 1.0 sec, is the highest peak in this diagram, from before and after the BNS merger (see Figs. 23, 24 and 25). A secondary peak at the same frequency but t − t merger 32.9 sec is also highlighted in this plot.  ; combining all harmonics with frequencies n × f , with n ∈ N) for the first peak (red) at 1.0 sec after the merger for the BNS merger gravitational-wave event GW170817, observed by x-correlating the LIGO detectors [3]. The same amplitude-frequency plot for a random time in data (yellow), is also shown for comparison. Solid area between 63 Hz and 92 Hz was the search frequency prior range for Planckian echoes.  [3]. After this detection Gill et al. [145] with independent Astrophysical considerations have also determined that the remnant of GW170817 must have collapsed to a BH after t coll = 0.98 +0. 31 −0.26 sec. Error-bar (in blue) is the time of collapse considering this independent observation in [145] compared to the detected signal of echoes which is also as a consequence of BH collapse. The shaded region is 0-1 sec prior range after the merger.   Abedi and Afshordi [3] have considered the possibility of first and second scenario.
As it was pointed out in part 5.1.1 the time delay for Planck scale echoes is given by, [3], There are two natural frequencies for the waveform of the echoes: The resonance frequencies (natural harmonics) of the echo chamber (formed by the angular momentum barrier and the near-horizon quantum structure), and the BH ringdown (or classical QNM) frequencies. The high frequency harmonics that are initially excited by the merger event decay quickly, while the low frequency harmonics live for longer time [146][147][148]. While the former captures the repeat period of the echoes, the latter describes the echo internal structure. Given that the ringdown frequencies are not resolved by LIGO detector, we can roughly approximate the observable signal as a sum of Dirac delta functions, repeating with the period ∆t echo : . (152) Therefore, the method searches for coherent periodic peaks of equal amplitude in cross-power spectrum of the two detectors at integer multiples of f echo ≡ ∆t −1 echo , in following steps: 1. Fundamental frequency of echoes f echo = ∆t −1 echo within the 90% credible region range for final BH mass and spin is given by: 63 ≤ f echo (Hz) ≤ 92.
2. The prior range for echoes search is 0 < t − t merger ≤ 1 sec. 3. Using amplitude spectral density (ASD) Wiener filter (rather than whiten) the data by dividing by noise variance PSD=ASD 2 (rather than ASD): where δt is the time shift between detectors. 4. Cross-correlating the obtained spectrograms and sum over all the resonance frequencies of n × f , (154) Since the polarizations of the LIGO detectors are opposite for GW170817, the real GW signals appears as peaks in −X(t, f ) (see Figs. 22 and 24).
The simplicity of the method (not having any arbitrary or ad-hoc cuts or parameters) along with its high significance are reasons for making this finding more interesting and reliable. (For GW170608 we use only 512 s of data, which is all the noise-subtracted data available.) A signal peak tends to persist over various changes of the window parameters more so than a noise peak. Figure 11 shows an example of the persistence of the signal peak as a function of N E for GW170104, which makes clear that an averaging of the correlations over N E will improve the signal.
The window parameters used are summarized in Table I along with the best-fit value of t d , the p-value and the frequency bandpass for each analysis. The bandpass turns out to be around the most sensitive region for the detectors. For smaller (larger) mass events, the upper (lower) end starts to sample higher noise levels, but it is still away from where the noise gets significantly larger. As we have mentioned earlier, it is convenient to express the bandpass as a Figure 27. Correlation vs. echo time delay ∆t echo vs. N E for GW170104 using method II [4]. Here N E is the number of frequency steps between spikes. 5.1.6. "GW echoes through new windows" by Conklin et al. [4,5] In this search, three methods (named as I, II, III) that are based on general properties of echoes has been suggested. Conklin et al. [4] have mostly focused on Method II, which is based on frequency windows, while the other two methods use time windows. Window functions in these methods allow us to find quasiperiodic structures in time and/or frequency domains. Method II turns out to be the most successful one. Accordingly, in this paper we just review this method. The search methods become more optimal using correlations of data in multiple detectors. Using the suggested method and search [4] find significant evidence for GW echoes, which is shown in Table 7 (and see Fig. 27 Table 7. The best-fit ∆t echo , p-value, bandpass and window parameters for the six signals [4]. Here N E is the number of frequency steps between spikes which has been chosen using injection and echoes model properties in [4]. 51 of 84 size of the error bars has nothing to do with the strength of the resonance signal for each event. The value η = 1.72 ± .06 means that δr ≈ 10 −28 Pl ≈ 10 12 × (proper Planck length). These results provide a good test of the spin and redshift dependence in (1). This is shown more clearly in Fig. 6, where the left plot is the same as  The time delay of echoes, ∆t echo which is similar to the Planckian echoes in Eq. 143, is paramaterized in [4] using where r 0 is the location of the quantum structure outside r + . Here η = 2 corresponds to proper Planck length (Eq. 143).
Taking into account of errors in final mass, spin, and redshift it is realised that (see Fig 28) the echoes found are consistent with η = 1.7 [4,5]. Best fit properties of the peaks for O1 and O2 events also shown in Table 8 [5]. Alternatively, this result can be interpreted as the energy scale for reflection from quantum horizons to be 6 ± 2 orders of magnitude below Planck energy [89].

Comment on: "Gravitational wave echoes through new windows by Conklin et al. [4,5]"
Conklin et al. [4] remark that "We have not found signals for the two earlier events, GW150914 and GW151012, which play a significant role in ADA results in [1]". While in their updated search [5] they indicate existence of signals for both GW150914 and GW151012, although with no p-value estimation which is crucial in this search. Therefore, it is hard to truly evaluate the significance of echoes reported in [5]. Moreover, it is not clear how much the choices made in Method II in [4] might have been affected by a posteriori statistics.

Results of Westerweck et al. and
Nielsen et al. [6,7] Westerweck et al. [6] re-analysed the same model proposed by ADA [1], using more background data and a modified procedure. They focused on the data analysis methods of ADA [1] and their significance 52 Table 9. Comparison of p-values obtained in [1] and using larger portion of data (4096 seconds of LOSC data) [6]. This data is divided into segments of 16  estimation, namely the concerns presented in [26] suggesting a different significance estimate using 4096 seconds of LOSC data. The results of p-value estimation and comparison with ADA results are shown in Table 9. In addition, Nielsen et al. [7] have searched for echoes signals in GW data via Bayesian model selection probabilities, comparing signal and no-signal hypotheses using ADA model [1]. Accordingly, calculation of Bayes factors for the ADA model in O1 events presented in Table 10.
In the following we explain the results given in different plots: 1. In Fig. 29 it is shown that depending on the overall amplitude of the injection, the signal either can be recovered or it would be difficult to recover. 2. Fig. 30 shows injected and recovered values for γ having different overall amplitude A used in [1].
Although this plot shows a preference for γ = 1, having high value of γ can be recovered easily. 3. Fig. 31 shows injection of echoes signals into the Gaussian noise having different overall amplitudes A.  Table 10. Results of Bayes factor [7] using ADA model. Gaussian noise hypothesis is preferred for negative values of Log Bayes factor. Echoes hypothesis is preferred for positive values of Log Bayes factor. Log Bayes values of < 1 are "not worth more than a bare mention" [7].

3
III. GENERAL REMARKS mediate problem arises regarding how strong signal should be for the three events. The two hole events GW150914 and GW151226 were the Advanced LIGO detectors with signifi-> 5.3σ and signal-to-noise ratios of 23.7 and ively [3]. The other event, LVT151012, had ignificance of only 1.7σ and a signal-to-noise combined between the two Advanced LIGO owever, in Table II of [21] we see that the ise ratio of the claimed echo signal is actually VT151012. er SNR of LVT151012 cannot be due to t projected number of echoes between the different ∆t echo leads to differing numbers of given duration: the 32 seconds of data used in (∼ 180) for LVT151012 and (∼ 110) for Although the number of echoes is larger for late echoes are strongly damped. They defactor of 10 over ∼ 22 echoes for the claimed litude γ ∼ 0.9. Thus in order for the echoes 2 to have a higher SNR than the echoes of their amplitude must be very high. In fact for the reported SNRs, the initial amplitude echo of LVT151012 would have to be about than that of GW150914 [28], while the origieak is about 2-3 times lower for LVT151012 on to GW150914's. This would require their to be about 2-3 times larger for LVT151012 150914. This seems to be confirmed by the h results in Table II of the updated work [21], A GW150914 = 0.091 and A LVT151012 = 0.34. e that far in the wave zone the gravitational of the echoes decays similarly to the signal of self, i.e. linearly with the distance from the s explicit astrophysical assumption, in addiin [20][21][22], is the basis for the above concern. ignificance of LVT151012 is rooted in its disean estimated distance being more than twice that of GW150914 and GW151226, we execho signals. While particular combinations arameters and signal morphologies may have ffects on the generation of echoes and their litudes, changing their relative significance, no extensive model to justify abandoning this e. red amplitude parameters suggest that a lot nal wave energy was emitted in the echoes: a alculation implies that the amount of energy he echoes was approximately 0.1 solar masses 914) and 0.2 solar masses (for LVT151012). be compared to the total estimated energy the original signal of 3 solar masses (for Best fit SNR 2

Maximised SNR 2 for injections in Gaussian noise
The matched filtering technique is able to recover signals with a variety of amplitudes. As shown here, the SNR depends on the amplitude of the signal. The amplitude found by ADA (A = 0.1) is close to the level that is found in pure Gaussian noise. An amplitude twice as large as this would be clearly identifiable in the data.
LOSC [23]. The parameters of the echo templates, in particular ∆t echo , depend on the mass and spin parameters of the final black hole. Instead of using only one initial waveform and generating all echo templates with this, one should use an initial waveform that corresponds to each set of echo parameters to be varied over. Using the single LOSC waveform is a simplification, restricting to only one choice of final mass and spin parameters for the echoed original event, while simultaneously varying over the final mass and spin values through ∆t echo .

IV. VALIDATION OF THE MATCHED-FILTER ANALYSIS
We wrote a separate implementation of the ADAsearch procedure, that we refer to as ADA AEI -search. No changes were made to the algorithm as described before, while the implementation itself is independent. The SNR 2 -results obtained with our implementation are similar to those shown in [20].
As a first check, we verify that the ADA AEI -search procedure can distinguish between pure noise and simulated echo-signals. For this, a known signal is injected into simulated noise. We simulate Gaussian noise with a Power Spectral Density (PSD) similar to that found for the detector data around each event (calculated from Figure 29. This plot shows whether it is possible to recover the potential signals with a variety of amplitudes in [6]. Here it can be seen that amplitudes less than A=0.1 in ADA [1] are difficult to be identified in data, while amplitude twice this value would be clearly identifiable. 4 the best-fit results of [20,21]. Fig. 2 shows ence of the SNR 2 peak on the injection amj . The effectiveness of the method in finding epends on A inj . This test was performed for ealisations of the simulated noise. The minirequired to find a peak rising above the noise d also depends on the noise instantiation. We inj ∼ 0.1 can yield a visible peak. This is the ue of A reported for GW150914 in [20]. In one five trials conducted in this first test, however, plitude was necessary to distinguish the sigoise, as shown in Fig. 2, where the noise and t injection have almost identical SNR 2 results. pted us to perform more detailed statistical d injection-recovery analyses, as described be-

PRIOR RANGES AND TEMPLATE SPACING
or each echo parameter are determined from rior range. Each template in the bank is for a specific value of each parameter. The ltering method finds a higher SNR for data the template, but each template can recover h a range of parameter values. The three pa-, t 0 and t echo are determined by maximisation, t 0 kept fixed between the different events. In arameters recovered are defined as the values Recovered γ γ-recovery depending on A inj FIG. 3. Injected and recovered values for γ, correct recovery would be seen as a diagonal. The search method's preference for γ = 1 (dashed line) at lower injection amplitudes is clearly seen.

−20
Amplitude-recovery depending on A inj H1 injection/recovery L1 injection/recovery Figure 30. This plot shows injected and recovered values for γ. The diagonal line is accurate recovery [6]. The preference for γ = 1 (dashed line) at lower injection amplitudes can be clearly seen.

of 84
RIOR RANGES AND TEMPLATE SPACING each echo parameter are determined from ior range. Each template in the bank is r a specific value of each parameter. The ering method finds a higher SNR for data e template, but each template can recover a range of parameter values. The three pat 0 and t echo are determined by maximisation, 0 kept fixed between the different events. In ameters recovered are defined as the values g to the template in the bank which yields SNR. The maximisation is performed over s in the bank and thus over all values in the rid used to create the bank. The boundaries eter grid are determined by a prior range, nges chosen by ADA are displayed in Table   s for γ and t 0 resulting from this maximisand to lie very close to the boundary of their 0.9 and −0.1 respectively [28]. This suggests ay be support for values of these parameters ide of this range. If these values reflect the than the data, then they cannot be reliably s evidence for a detection claim. Furthere greater than unity for γ means that each ho has an amplitude greater than the previch a result would require the echo signal to g energy from the black hole spacetime. whether the preference for these parameter rtifact of the method, again using known siginto simulated noise. We constructed Gausith a PSD estimated from the 4096 seconds ta around GW150914. The injected signals o signals based on the LOSC GW150914various echo parameters. The range of γ to γ ∈ (0.1, 2.0) both in the prior of the he injections. The range of t 0 is widened to )∆t echo,theory in the search. It is not widened

Comment on: "
Low significance of evidence for BH echoes in gravitational wave data" [6] 1. Comment on search strategy: As described in former part in Fig. 30, Westerweck et al. [6] demonstrated the preference of γ = 1 at lower injection amplitudes for ADA model [1]. This is expected result as γ → 1 extends the template to infinity. In other words, it extends the template range to infinite time, which is clearly dominated by noise. However, γ = 1 is still far from the best-fit γ = 0.9 on edge of the prior where at least 90% of energy goes to first 11 echoes. Besides, initial waveform must change significantly in subsequent echoes. Indeed, repeated template that does not damp (γ = 1) is not physical. One solution to this problem might be to use a finite range of data (which is used by Westerweck et al. [6]) making the result strongly dependent on what portion of data has been taken.
It should be also noted that using 1-sigma range for errors in ∆t echo , implies a 32% chance for the signal to be missed that causes reduced significance by diluting SNR 2 from some events.
Westerweck et al. [6] have found that the least significant event LVT151012 which is now called GW151012 has the most contribution to tentative evidence for echoes. This peculiar finding does not disfavor echoes as there is no simply reasonable justification that significance of echoes should be directly related to the significance of main event. Additionally, Wang et al. [148] have shown that by changing only ±20% of frequency of initial condition of echoes the SNR 2 for echoes can change by 3 orders of magnitude. As BBH events have different component spins and mass ratios we might expect a significant diversity in relative echo signal amplitude for each of them. Interestingly, as will be discussed in next part 5.2.3 mass ratio of BBH events appears to show correlation with the echo amplitude [139].

Comment on abstract and conclusion:
The most crucial comments for Westerweck et al. [6] goes to their abstract and conclusion (also provided in [28]). Although ADA [28] strongly acknowledge the analysis by Westerweck et al. [6], which is a careful re-evaluation of ADA analysis, the Abstract/Conclusion of Westerweck et 55   al. [6] misrepresents their finding. The most critical point of this misrepresentation is in Abstract claiming "a reduced statistical significance ... entirely consistent with noise". Contrasted to this claim in their Table I (Table 9 in this paper) they found p-value=0.020 for the noise hypothesis, with the same model and data as in ADA (as opposed to 0.011 in ADA [1]). However, if one follows standard nomenclature (e.g., https://en.wikipedia.org/wiki/P-value#Usage or [149]), p-values < 0.05 disfavour noise hypothesis, providing "moderate to strong" Bayesian evidence for echoes, which is contrary to what they state in their Abstract.
To conclude, considering all the critiques of Westerweck et al. [6], we see NO evidence that their improved analysis with p-value= 0.020 ± 0.009 has reduced the significance of echoes, entirely consistent with p-value = 0.011 of ADA [1]. The fact that completely independent events of O2 also show a low p-value= 0.039 for ADA echoes (see Table 6 [2]) further boosts the statistical evidence for ADA model in LIGO/Virgo data.

Results of Salemi et al. [8]
Another independent group [8] has found similar post-merger GW signals that can be attributed to GW echoes. However, the setup of this methodology, which is based on coherent WaveBurst (cWB) [150] method, was not originally developed to search for echoes. This search, which is independent of the waveform models, has been developed based on coherent excess power in events from the GWTC-1 (catalog of compact binary coalescence). Here, loose bounds on the duration and bandwidth of the signal leads to evaluation of coherent response of independent detectors. This search has focused on detected features as deviations from GR and has presented the method to obtain their significance. It appears that from eleven events reported in the GWTC-1, two of them (GW151012 and GW151226) in Figs. 32 and 33 respectively, show an excess of coherent energy after the merger (∆t 0.2 s and 0.1 s, respectively) with p-values (0.004 and 0.03, respectively). However, [8] have shown that (Fig. 34) the post-merger signal from GW151012 favours different sky location than that of the main event.
In TableI [8]. The upper plot shows the squared coherent network SNR and the plot bellow shows the normalized residual noise energy. The residual plot is given after the reconstructed signal was subtracted from the data. In this plots a secondary cluster occurring 200 ms after the merger (consistent with echo times predicted and seen by ADA, Equation 144). The dashed vertical lines denote coalescence time for GW151012 (the network has used the Livingston detector time as a reference).   [8]. The upper plot shows the squared coherent network SNR and the plot bellow shows the normalized residual noise energy. The residual plot is given after the reconstructed signal was subtracted from the data. In this plots a secondary cluster occurring 100 ms after the merger (consistent with echo times predicted and seen by ADA, Equation 144). The dashed vertical lines denote coalescence time for GW151226 (the network has used the Livingston detector time as a reference).  In this part, we first re-examine the interpretation of Salemi et al. [8] about the signals they found and then give a more conclusive support for echoes hypothesis.
• Comment on: Results of Salemi et al. [8] Salemi et al. [8] disfavoured echoes hypothesis pointing that post-merger signal of GW151012 has arrived from a different sky location than that of the main event. However, the p-value∼ 0.004 of this secondary signal, disfavours two signals being unrelated. We see that all the secondary (post-merger) clusters they claim as signals in Figs. 32 and 33 are nearly monochromatic. That means the waveform of these signals are quasi-periodic, leading to degeneracies in inferred time-delays. So looking again to the secondary signal of GW151012, the null (residual) plot in Fig. 32 shows the peak of the cluster (which is mostly responsible to make a different sky localization) is at ∼ 130 Hz which corresponds to 7.7 m sec time delay. This 7.7 m sec is the same time delay of first peak and second peak in Fig. 34 for post-coalescence signal (green), confirming that the monochromatic degeneracy in Fig. 34 might have caused the different sky localization. Fig. 35 shows better interpretation of cause of this error.
• Dependence of significance of echoes on binary BH mass ratio Consider a system of binary BHs (BBH), with progenitor masses m 1 and m 2 . It is known that these systems consist of two almost equal mass BHs m 1 ∼ m 2 . However, the inevitable diversity in the initial conditions, specifically binary mass ratio, can lead to different echo properties. Here, we review the evidence for correlation between the significance of echoes and the progenitor BBH mass ratio (first presented in [139]): 1. We have used LIGO parameter estimation samples for BBH events provided in [151]. Then obtained mass ratios by weighting all events as equal. We used full m1 vs m2 distribution samples for "Overall_posterior". We fit these posterior points to a straight line shown in Fig. 36.  2. We take p-values reported in [8] for each event post-coalescence signal and plot the best fit line of mass ratio vs (− log(p − value)) 0.5 . We used least square method [152] in fitting a straight line consisting all the posterior points provided by LIGO. Finally, we use the slope of best fit line as our primary parameter in determining significance.
3. Finally, we obtain the significance (Fig. 37) assuming that there is no relation between p-value and mass ratio by taking random p-values for BBH Catalog events within the uniform range 0 < (− log(p − value)) 0.5 < 2.5 (0.0019 < p-values< 1). Indeed, having no relation between p-value and mass ratio of events shall end up with zero slope in large number of random selections. Therefore, in order to find a false detection rate for correlation, number of slopes higher than the actual measured slope is calculated (see Fig. 37). Accounting for the "look elsewhere" effect, we find tentative hint of mass-ratio dependence of echo significance reported in [8] at false detection probability of 1%. Table 12 shows events and p-value of their post-coalescence signal reported in [8] versus expected Planckian echo time delays and average mass ratios. At a glance it is seen that smallest mass ratios go to smallest p-values. This can happen with 1/10 chance out of 10 BBH events. For two most significant events, being also the most extreme BBHs, the random chance becomes p-value= 1 9 × 1 10 = 0.011 which is consistent with the statistics of regression analysis in Fig. 37.

5.3.1.
Template-based gravitational-wave echoes search using Bayesian model selection by Lo et al. [9] Lo et al. [9] found that using a wider range of priors (listed in Table 13) compared to that of ADA model [1], including the main event template in the echo template, and keeping just three echoes, leads to lower significance on echo signals evaluated by the Bayes factor using Bayesian analysis. They have considered two hypotheses as the null hypothesis H 0 and alternative hypothesis H 1 ,  . Histogram of slopes for uniform random choices of 0 < (− log(p − value)) 0.5 < 2.5. We see that only 1.3% of these random realizations, the slope can exceed the observed value [139].
where d and n indicate the GW data, and the instrumental noise respectively and h IMR , and h IMRE are the inspiral-merger-ringdown (IMR) gravitational-wave signal and inspiral-merger-ringdown-echo (IMRE) gravitational-wave signal respectively. The log Bayes factor ln B in their search method gives the detection statistics which determines whether there is an IMRE signal or an IMR signal in data. From the Bayesian perspective if the log Bayes factor, is greater than 0, we can conclude that the data favor the alternative hypothesis.
The relation between p-value and null distribution of detection statistic ln B is given by where ln B detected is the detection statistic which is obtained using a segment of data in analysis, and p(ln B|H 0 ) is called the null distribution of ln B, i.e., the distribution of ln B assuming that H 0 is true. Lo et al. have injected an IMRE injection of template with echo parameters discussed earlier for the event GW150914 into simulated Gaussian noise. The detection statistic for this injection is given as follows, Therefore, this finding indicates that the data slightly favor the null hypothesis from Bayesian analysis point of view. While the p-value and the corresponding statistical significance, for Gaussian noise, show that the data favors echo hypothesis,  Table 15. The detection statistic and its corresponding statistical significance and p-value for O1 events [9]. The ordering of events by their statistical significance is consistent with what reported by Nielsen et al. [7] Table 14 shows the values of the detection statistic ln B vs corresponding statistical significance in Gaussian and O1 backgrounds. Therefore, for a detection of gravitational-wave echoes having statistical significance ≥ 5σ, the detection threshold would be, ln B threshold,Gaussian = 1.9, for Gaussian noise and O1 noise respectively. Table 15 shows the detection statistics and the corresponding statistical significance and p-value for the O1 events. This table also shows that the ordering of the events by their statistical significance is consistent with what has been reported by Nielsen et al. [7]. ln B (cat)

Comments on Lo et al. [9]
In their analysis, Lo et al. [9] have included both the main event, as well as the ADA echo waveform in their template, but they used expanded priors in Table 13. Although expanding the ADA priors covers a larger space of possibilities, it tends to dilute marginal signals and bury them in the noise. For example, there is no good physical interpretation for repeating echoes that do not damp (with γ = 1), as they violate energy conservation.
For these reasons, it is not surprising that Lo et al. [9] find a smaller Bayes factors, due to their expanded priors. However, it is well-known that expanding prior into nonphysical regimes will artificially lower Bayesian evidence for any model, especially since γ = 1 is a (formal) singularity of the likelihood function.  Table 17. P-value for each event and total p-value [2]. Result 1 is the case when the phase shift is fixed to π, and result 2 is the case when the total phase shift is also a parameter.

Results of Uchikata et al. [2]
As discussed in Section 5.1.3 above, Uchikata et al. [2] can approximately reproduce the evidence for ADA echoes in both O1 and O2 events. Here, we present their search using an alternative template that failed to find any evidence for echoes. In order to build the latter, they considered Kerr spacetime, replacing the event horizon with a reflective membrane. They then used the transmissivity of the Kerr angular momentum barrier T BH (ω) to filter the ADA template, which acts as a hi-pass filter (truncating the low frequency part of the ADA phenomenological waveform 147). Moreover, the overall phase shift of the waveform as a free parameter is taken into account contrary to ADA search. Using this template, they have found no significant echo signals in the binary BH merger events. The background estimation, has used the same method provided by Westerweck et al. [6]. Table 16 gives p-values for all events. The combined p-value is well above the critical p-value 0.05. In other words, echo signals using this model do not exist in the data, or their amplitudes are too small to be detected within the current detector sensitivity. 2. Since the phase shift at the membrane (due to the reflection and boundary condition) is model dependent, it is physically reasonable to assume a total phase shift as a parameter (see Fig. 21 that has used extra phase parameter for GW150914). In contrast, former studies [1,6] only considered phase inversion (or Dirichlet boundary conditions) at the reflective membrane. Therefore, the results of two cases, when the phase shift is fixed to π (result 1) and when it is a free parameter (result 2), respectively, in Table 17  Model provided by Uchikata et al. [2] substantially truncates the low frequency part of the GR waveform (which is the basis of ADA template 147). However, one may argue that the GR waveform has already been filtered once by the transmissivity of the angular momentum barrier, T BH (ω) as it is what is seen by observers at infinity. In fact, [29] have shown that both GR signal (main event) and echoes can be constructed from a superposition of the QNMs of the quantum BH, which are essentially the modes trapped between the angular momentum barrier and the quantum membrane (see Section 4.4 above). Therefore, both the GR signal and the echoes pass through same barrier and are thus truncated by the same T BH (ω). This implies that T BH (ω) cancels in the ratio of echo to main signal waveform, and in contrast to Uchikata et al. [2], no truncation is needed.

The results in
Indeed, as evidenced by their own analysis (Section 5.1.3 above), the low-frequency part of the ADA template is necessary to obtain a significant signal. This is physically justified since the high frequencies leak out of the angular momentum barrier quickly, leading to a rapid decay, while echoes can last much longer at lower frequencies [148]. 5.3.5. Results of Tsang et al. [10] Tsang et al. [10,140] proposed a morphology-independent search method which consists of a large number of free parameters for echoes compared to ADA model [1] (49 versus 5). They search for echoes in all the significant events in (GWTC-1), and found that for all the events, the ratios of evidences for signal versus noise and signal versus glitch do not rise above their respective background. Only the smallest p-value=3% goes to the event GW170823. Hence they found no significant evidence for echoes in GWTC-1. The results of search are given in Table 18 and 19. 5.3.6. Comment on: Results of Tsang et al. [10] Tsang et al. [10] have developed a model which consists of a large number of free parameters (49 by our count). Indeed, a larger space of possiblities leads to lower significance. So it is not surprising to get large p-values out of large free-parameter space. Indeed, all the SNRs for echoes reported in ADA would be below the detection threshold reported by [140], given their large number of parameters (see [139] for more detail).

A concordant picture of Echoes
In order to have an optimal search for echoes, one may want to take the following guidelines into consideration: 66 Table 19. Same as Table 18, while for the events that three detectors has been involved [10]. In the case of GW170817, in order to cover 1.0 sec after the merger for echoes found by Abedi and Afshordi [3], additional search as first echo being from this time has been set. For this particular event, latter prior choice has been taken for combined p-values.
indicate that the echo signals found in [1] mostly consist of low frequency modes, which is independently confirmed in [3], and Uchikata et al's own search in both O1 and O2, using the original ADA template.

Non-Gaussianity of backgrounds
In order to search for signals of a given echo model, we need a proper understanding of background behaviour. Only then we might be able to determine the best statistical methodology. Since ADA [1] used different parts of LIGO data to get a combined significance, one may think about what would be the best method of combination of separate sets of data with different background behaviour. Then, we can also ask how the search changes by including three or higher number of detectors.
It was already observed that LIGO noise vary significantly and is very non-gaussian over long time-scales (see Fig's 14-15 in [142]). This either non-stationary or non-Gaussian background makes the interpretation of p-value ambiguous, particularly in finding marginal echo signals which are often near the detection threshold. Therefore, one must examine how much this varying background affects the inferred significance of a detection. This is studied by looking at other different minute-long stretches of data within a minute of the main events [1]. As can be seen in Fig. 38, the variation of p-values at the tail of the distribution is much higher than what is expected from Poisson statistics of the SNR peaks. This becomes more interesting when we see that the smallest p-value is coming from the range which is closest to the main event. This might be because the marginal LVT151012 (now called GW151012) detection is over a minimum of the LIGO (combined) detector noise.  Figure 38. p-value distribution for combined events of different stretches of data within 1 minute of the main events. Surprisingly, the blue line which is closest to the main event, and has used to define p-value in [1] (Fig. 19), happens to give the smallest p-value. The shaded region represents the Poisson error range for blue histogram. This shows that the variation in p-values is clearly much larger. This behaviour is interpreted as non-gaussianity and/or non-stationarity of the LIGO noise. In this plot the y-axis on the left (right) shows p-value (number of higher peaks) within the mentioned range of data. In each histogram the total number of "peaks" is (38 − 9)/0.02 = 1450. 71 of 84

Towards Synergistic Statistical Methodologies
As we summarized in the previous sections, the past fours years have witnessed hundreds of theoretical studies focusing on model-building for echoes, as well as dozens of observational searches and statistical methodologies. However, in spite of remarkable progress on both fronts, the theoretical and observational tracks have largely developed independently. However, it appears that both tracks have become mature enough, so much so that the time is ripe for a synergistic convergence. For example, Bayesian methods developed in [7,9] applied to a superposition of QNMs of quantum BHs (as outlined in [29]) would put coherent methods developed by [4,28] on more sound statistical and physical footings. The analogy will be with helio-or astro-seismology, where modeling a dense spectrum of QNM frequencies can be used to infer the intrenal structure of the compact objects [128].
The real challenge will be in allowing enough freedom in our best physical models, in order to capture all the remaining theoretical uncertainties, but not any more!

Echoes in Numerical Relativity
Most studies of echoes have so far focused on the linear perturbation theory around the final BH for simplicity, but in reality the mergers start with the highly nonlinear binary BH inspiral. Hence, we need a covariant numerical implementation of binary quantum BHs within a highly-nonlinear dynamical spacetime to fully address the entire dynamics, especially the initial conditions. There are several possible approaches borrowed from numerical relativity which can be modified to either include the quantum boundary condition or the full dynamics of binary quantum BHs.
For instance, the effective one body (EOB) formalism [156,157] is a concrete strategy which only needs to solve ordinary differential equations rather than to perform the costly 3d numerical relativity simulations. It uses higher-order post-Newtonian expansion in a resummed form (different from the usual the Taylor-expansion), to include the non-perturbative result using a conservative description of binary BHs dynamics, radiation-reaction and emitted GW waveform. One possible approach, that is currently underway, is to capture the nonlinear effects in echoes by modifying the boundary condition in the EOB codes to implement the quantum BH dynamics.
Another route is to directly modify numerical relativity codes that have successfully produced waveforms for BBH merger events. A concrete strategy could be incorporating the mock fuzzball energy-momentum tensor (Section 3.6) as a source for Einstein equations, directly into the numerical relativity codes. If the fuzzball "fluid" manages to stay just outside the apparent BH horizons in a dynamical setting, then it can potentially generate echoes in a fully nonlinear numerical simulation of quantum BBH merger.
Recently, [158,159] presented the first numerical simulation of BBH mergers in Chern-Simon gravity. They start with the modified action and predict the dynamics order by order. It is possible that a similar iterative approach can be applied to model boundary conditions at apparent horizons, or evolution of mock fuzzballs.

Quantum Gravity, Holography, and Echoes
As we discussed in Section 4 above, any modification of event horizons that could lead to echoes should be a non-perturbative modification of general relativity, and can only be fully captured by a non-perturbative description of quantum gravity. A possible example of this is the fuzzball program in string theory (Section 3 above). But more generally, what can non-perturbative approaches to quantum gravity tell us about BH echoes? 72 of 84 One of our greatest insights into the dynamics of quantum gravity has come from the Holographic Principle, that extending Bekenstein-Hawking area law for entropy of BHs [42], suggests the entire dynamics of a quantum gravitational system should be captured on its boundary. The most concrete realization of this principle was proposed by Juan Maldacena [94], in the form a conjectured duality between quantum gravity in Anti-de Sitter (AdS) spacetime and a Conformal field theory (CFT), commonly known as AdS/CFT correspondence or conjecture. It proposes that CFT in spacetime of d-1 dimension, at the asymptotic boundary of an AdS spacetime is mathematically equivalent to string theory (or quantum gravity) within the bulk AdS in d dimension. This topic has been extremely fruitful over the past two decades, offering many synergies between seemingly disparate notions in geometry and quantum information. For example, the Ryu-Takayanagi conjecture [160] relates the entanglement entropy of boundary CFTs with the areas of extremal surfaces in the bulk AdS, generalizing the notion of Bekenstein-Hawking BH enetropy to arbitrary geometries.
An intriguing connection between AdS/CFT and echoes is the appearance of echo times: as "scrambling time", in the AdS/CFT literature [161]. Here, S BH and T H are the entropy and temperature of the BH respectively. The scrambling time refers to the time it takes to destroy quantum entanglements in a chaotic system, while BHs (and their CFT duals) are conjectured to be fast scramblers, i.e. the most efficient in destroying entanglement (e.g., [162]). Interestingly, Saraswat and Afshordi [163] have recently shown that the scrambling time (computed using Ryu-Takayanagi conjecture in a dynamical setting) is identical to the Planckian echo times, for generic charged AdS BHs. Could this imply that echoes could be a generic property of (possibly a certain class of) quantum chaotic systems? Another possible connection could come in the form of the fluid-gravity correspondence, e.g., in the context of membrane paradigm discussed in Section 2.2. For example, in [23], we have argued that Boltzmann reflectivity of GW echoes, implies that viscosity of the boundary fluid should vanish at small frequencieshω kT. One may also speculate that other holographic manifestations of BH echoes may appear in the Kerr/CFT conjecture [164], Braneworld BHs [165], or as Regge poles of the boundary plasma in AdS/CFT.

Einstein Telescope, Cosmic Explorer
The Einstein Telescope (ET) [166] and Cosmic Explore (CE) [167] are the third-generation ground based GW detectors. The ET consists of three underground detectors with three arms 10 kilometers long and CE will be realized with two arms 40 kilometers long, which are 10 times longer than Advanced LIGO's. These next-generation GW detectors might allow us to observe some Planckian signatures from quantum BHs such as GW echoes from merger events leading to a remnant BH. We plot the spectra of GW echoes and ringdown with the sensitivity curves of Advanced LIGO, ET, and CE in FIG. 39. The detection of GW echoes with the third generation GW observatories are discussed in [117,118], and it may be possible to distinguish ECOs with |R| 0.3 from BHs with at 2σ level when SNR ∼ 100 in ringdown, which would be possible for the third-generation GW detectors. The relative error on the reflectivity of would-be horizon is also investigated in [117,118], and the relative error for measurement of relectivity in ground-based detectors is approximately given by where M = 30M , ρ ringdown is the SNR in the ringdown phase, while the distance between the top of the angular momentum barrier and the would-be horizon is assumed to be longer than 50M in the tortoise coordinate. For comparison, we note that the loudest detected BBH event, GW150914, has ρ ringdown 8. The detectability of GW echoes from failed supernovae, leading to the formation of BHs, with the third-generation GW observatories is also discussed in [128]. Calculating the SNR of GW spectrum consisting of echo and ringdown, ρ ringdown+echo , in the Boltzmann reflectivity model, the horizon distance D h , defined as the distance where ρ ringdown+echo = 8, is estimated. Given the optimistic case in the Boltzmann reflectivity model, T H /T QH = e 15(ā−1) , the horizon distance can be estimated as D h ∼ 10 Mpc for the Advanced LIGO at design sensitybity and D h ∼ 100 Mpc for the third-generation detectors such as ET and CE. Therefore, the authors in [128] argue that the searching for GW echoes, sourced by failed supernovae within our Galaxy and nearby galaxies, may be possible. However, in the case of T QH = T H , the horizon distance is less than or comparable with 10 Mpc and so the echo search with failed supernovae would be restricted to within the Local Group. For the comparison, the strain amplitude of GW echoes in T QH /T H = 1 and T QH /T H = e 15(ā−1) are shown 13 in FIG. 40.

LISA
The Laser Interferometer Space Antena (LISA) is planned to be the first GW observatory in space. It will have three satellites separated by millions of kilometers and their orbits maintain near-equilateral triangular formation. LISA might enable us to reach high-precision detection of ringdown in SNR ∼ O(10 3 ), which puts stronger constraints on the reflectivity of BHs [117,118].
Recently, a novel proposal to discriminate BH horizons based on the tidal heating was proposed in [168]. One of the main targets of the LISA mission is precision measurements of extreme-mass-ratio inspirals (EMRIs), in which the tidal heating could be important. The (partial) absorption of ECOs or BHs plays the role of dissipation at the surface, by which tides back-react on the orbital trajectory. It is argued that this tidal heating is responsible for a large dephasing between the orbits of a BH and ECO. This dephasing accumulates over the timescale of months and the accumulation speed is faster for a higher spin. The authors in [168] also found a proportionality relation between the dephasing δφ and energy reflectivity |R| 2 .
In order to make use of this scheme to put strong constraints on the reflectivity of ECOs, one has to obtain accurate EMRI waveforms by properly taking into account the tidal heating for orbiting objects, which may decrease systematic errors in data analysis.
Not only the tidal heating, but also the tidal deformability contributes to the GW Fourier phase and it can be characterized by the tidal Love number k. The Love number of ECO of mass M may scale as 1/| log δ|, where δ ≡ r 0 − r h , with r h is the BH horizon radius of mass M and r 0 is the radius of the ECO. So the k − δ relation is δ = r h e −1/k , (164) and assuming this relation, one can infer the near-horizon structure characterized by δ from the measurement of the Love number k. For instance, if the Love number of the order of k ∼ 10 −2 is measured by LISA from a supermassive BH binary signal, leading to the formation of a BH of M ∼ 10 6 M , it yields the resolution of δ ∼ l Pl . However, the authors in [169] point out that the statistical and quantum mechanical uncertainties in measurements of near-horizon lead to some difficulty to measure δ precisely. The former one comes from the fact that the statistical uncertainty in δ is proportional to 1/k, and the inferred value of k, where the inferred value of δ is comparable with its statistical uncertainty, is around k ∼ 0.2. Therefore, any inferred value of δ, derived from k that is smaller than ∼ 0.2, would be dominated by the statistical uncertainty. The latter one comes from the uncertainty principle in quantum mechanics. Once precisely measuring δ ∼ l Pl , it may lead to the uncertainty in the mass of the ECO, which then leads to the uncertainty in the binding energy. This results in the uncertainty in the orbital and GW frequencies, which means that one cannot measure δ precisely if it is much shorter than l Pl .

Pulsar Timing Arrays
Following their first discovery in 1968 [170], over 2000 pulsars have now been detected by radio telescopes across the world. The pulsars' intrinsic properties, as well as propagation effects in the interstellar medium, can influence the arrival times of pulsar pulses. Therefore, pulsar timing arrays (PTA) can be used as a detection tool for BH binaries [171], and thus, might be used to detect singatures 75 of 84 of echoes from quantum BHs. In particular, millisecond pulsars stand out for their unparalleled stability (comparable to atomic clocks!) without being subject to starquakes and accretion. To give an explicit example, we show the spectrum of GW echoes predicted by the Boltzmann reflectivity model [23,29] with the sensitivity curve of International Pulsar Timing Array (IPTA) and Square Kilometre Array (SKA) (FIG. 41). The lower curve in Fig. [170] is for a 3 9 M BH merger at D o = 1 Gpc. Given that this mass is comparable to that of M87 supermassive BH, located at 16 Mpc, we expect ∼ 2 × 10 5 of such BHs at < Gpc. Assuming that each BH merges once every Hubble time ∼ 10 10 years, and that echoes last for 20 years (from simple mass scaling), the chances of detecting such a loud event with PTAs at any time is 0.1%. However, fainter events will be more prevalent as their number increases as SNR −3/2 from volume scaling. Furthermore, increase in supermassive BH merger activity observed at high redshifts shall boost this statistics. PTAs are anticipated to detect the low frequency GW signal from supermassive BBH within the next few years [171]. We expect that the first GW detection will be a stochastic background of supermassive BH binaries. With any luck, this shall lead to new insights into the nature of quantum BHs and gravity.

Final Word
In this review article, we provided a comprehensive overview of the theoretical motivations for why quantum black holes in our universe may have different observable properties, in contrast to their classical counterparts in Einstein's theory of general relativity. The most prominent and potentially observable smoking gun for these quantum black holes comes in the form of gravitational wave echoes, which have been the subject of intense theoretical and observational scrutiny over the past few years. We provided a concise account of theoretical predictions, as well as the exciting and confusing state of observational searches for echoes in LIGO/Virgo observations. We closed by article by our vision of the 76 of 84 future of "Quantum Black Holes in the Sky", via a synergy of statistical methodology, quantum gravity, and numerical relativity, and in light of the next generation of gravitational wave observatories.
While this review article focuses on the gravitational wave echoes, as arguably the most concrete and promising signature of quantum black holes, other possible observable signatures can be (and should be) explored. For example, interactions of photons or neutrinos with near-horizon quantum structure could lead to signatures in radio images in Event Horizon Telescope observations [172], or ultra high energy neutrinos in Ice Cube observatory [173], respectively. However, these signals will be suppressed if Boltzmann reflectivity is assumed, as they havehω kT H . Another alternative to echoes may come through non-localities in non-violent unitarization, which would be observable far from the horizon (see Section 3.9). However, it is arguably difficult to pin down concrete predictions in this scenario.
To conclude, the world of Quantum Black Holes remains a wide open and largely uncharted territory, spanning from the dark corners of obscure mathematical structures to the nitty-gritty details of gravitational wave detector noise. It also holds the promise to crack the century-old puzzle of quantum gravity, and yet be imminently testable in the next few years. Therefore, the study of "Quantum Black Holes in the Sky" remains extremely exciting, active, and confusing, and is bound to provide us with new surprises in the new decade, and beyond.