Black hole information recovery from gravitational waves

We study the classical and quantum black hole information in gravitational waves from a black hole's history. We review the necessary concepts regarding quantum information in many-body systems to motivate information retrieval and content in gravitational waves. We then show the first step in an optimal information retrieval strategy is to search for information in gravitational waves, compared to searching for correlations in Hawking radiation. We argue a large portion of the information of the initial collapsing state may be in the gravitational waves. Using the Zerilli equation for particles falling radially into Schwarzschild black holes, we then describe a method to retrieve full classical information about infalling sources, including masses, infall times and angles.


Introduction
In response to the black hole information paradox, a vast body of literature discussing the problem in the fields of general relativity, string theory, holography and quantum information theory has appeared in the last 50 years (see [1,2,3,4,5,6,7,8,9,10,11,12] and references therein for a non extensive overview). Assuming subtle correlations between evaporation quanta which are not taken into account in Hawking's initial calculation [13], it is often expected that the final state can be made pure and information about the initial state retrieved (e.g. [7]). In doing so black holes are often implicitly portrayed as initially closed (or at least already equilibrated) thermal systems, allowing information to escape only in Hawking radiation. However since black holes are open systems throughout their history, they leak information at the early stages of their formation history, when gravitational waves are emitted. Further, although thermal many-body quantum systems and black hole perturbations have separately been extensively studied (see eg [14,15,16,17,18,19,20] for a review), to our knowledge the two fields are quite separate and there has been no study of the quantum information content of gravitational waves from black holes. We will nevertheless find that a large portion of the information of an initial gravitationally collapsing state is present and easily accessible in gravitational waves. This is because a large share of the information of the state is stored in the classical observables of the state, which form a subset of the total information, usually referred to as quantum information †, and which can be precisely recovered using black hole perturbation theory.
We believe that the gravitational wave community has therefore made extensive progress on accurately determining information from black holes in a way that is relevant to the information paradox, motivating further work to determine its content in GWs. In section 2 we review quantum information conservation in thermal systems and motivate information retrieval in GWs. In section 3, we review black hole perturbation theory in the case of a Schwarzschild black hole. In section 4 we show the algorithm for the inversion to obtain the masses of the infalling particle and black hole from the gravitational wave signal. Finally sections 5, 6 and 7 focus on obtaining the time and angle at which the particles fell in from the signal.

Concepts in density matrices, entanglement entropy and thermal systems
Naively, black holes at rest are featureless and seem to destroy information as they evaporate. In the original paradox, Hawking showed the density matrix of the Hawking radiation was diagonal [13], which he noted to be information-theoretically inconsistent with an initial "pure" quantum state. Therefore before discussing gravitational waves, we begin by an introduction to quantum information conservation and motivate the study of black hole perturbations.
Compared to wavefunctions, density matrices are generalized representations of states. States with well defined wavefunctions, called "pure states", have N × N density matrices ρ = |ψ ψ|, where |ψ = i a i |n i is the wavefunction in some basis |n i i∈ [1,...,N ] and N is the dimensionality of the Hilbert space. Off-diagonal elements ρ i =j can then appear as the phase coherence of the state. These phases are important as they guarantee the existence of a basis in which we can write ρ = 1 N ×N δ i,1 . Pure states, however, are more often characterised using Tr(ρ 2 ) = 1. More generally, a state with well defined probabilities but partial or absent phase coherence (e.g. classical statistical ensembles) will exhibit a partial or total basis independent absence of offdiagonal matrix elements and is called a partially or totally "mixed" state. Mixedness can be quantified as Tr(ρ 2 ) < 1. Unitary evolution of the system, such as the time evolution of gravitationally collapsing radiation and/or matter, implies mixedness and vice versa purity remain constant since Tr(ρ 2 ) = Tr(U ρU † U ρU † ) .
(1) † Note here that the term quantum information might be somewhat ambiguous as it refers to the total information contained in a system and therefore comprises the classical information as well as the purely quantum information. In this work we are interested in the share of classical information that can be retrieved.
This implies the mixedness/purity is a measure of information since information is assumed preserved in unitarily evolving systems. Further the Von Neuman entropy of a state, defined as S vN = −Tr(ρ ln(ρ)), is also a measure of mixedness of the system and remains constant. It is identically zero for a pure state only. But what of the classical entropy of a system described by a pure state which has thermalized? It is non zero as expected, and in fact can be found as the upper bound in the space of subsystems, which will be mixed, such that S BH (Ω) = max Ω 0 ⊂Ω S vN (Ω 0 ). Indeed, [21,4] showed that random subsystems from a thermal ensemble (such as blackbody radiation escaping a black hole) are mixed because they are entangled subsystems of a globally pure state. This is why Von Neuman entropy is sometimes called "entanglement entropy". This can be illustrated with a simple two particle-two spin state example with states |ψ AB = 1 √ 2 (|+− − |−+ ), taking the partial trace of the density matrix over the A states such that Tr A (AB) = B Tr(A) gives: where we expressed ρ AB in the basis {|−− , |−+ , |+− , |++ }, and summed over values (+ and -) of A for each value of B in the basis {|+ , |− } to obtain ρ B . Note this matches the classical entropy S cl (AB) = ln 2 = S vN (ρ B ) since AB has two possible microstates. Since the particles A and B are entangled, the density matrix of a subsystem can only contain partial information. Let us now assume a black hole is a thermalized body made of entangled quanta, which escape during evaporation. Since the global state is pure, as more particles escape the black hole the mixedness of the subsystem of evaporated quanta will eventually start to decrease after the "Page time" until purity is restored. In the above example, assuming AB evaporates one particle at a time, the detector will initially see S vN (ρ vacuum ) = 0, then S vN (ρ A or B ) = ln 2 and finally S vN (ρ AB ) = 0. This yield a so-called Page curve. However this is in contrast to Hawking's result which states the evaporation particles are maximally mixed, such that mixedness keeps increasing until the end of the evaporation, and emitting an entropy S BH (illustrated in figure 1). Nevertheless recovering entanglement at late times has led the community to believe the information could be recovered if the black hole could be shown to evolve unitarily from a pure state and was referred to as the "central dogma" of black hole information in [11]. However, this requires finding correlations between patches of the black hole metric, which could be transferred to the radiation. These correlations would naturally appear in a microscopic model of gravitational degrees of freedom that have become highly entangled through scatterings and eventually thermalizing. This also introduces Figure 1. We show in blue the evolution of the black hole thermodynamic (Bekenstein-Hawking) entropy derived from classical arguments S BH cl = A/4 [22,23], in orange the classical radiation entropy S HR cl from Hawking's calculation [23], and in green the Von Neuman entropy S vN (Ω 0 ) = Tr(ρ Ω0 ln(ρ Ω0 )) for a thermal pure state, as a function of the size of the subsystem Ω 0 ⊂ Ω.  Figure 2. Evolution of the average entanglement entropy of half of a many-body system as a function of size L (data taken from [39]). S vN is a measure of the local correlations with the rest of the system. As the system thermalizes it converges up to small quantum fluctuations corresponding to the (small) interactions still occurring in the system, and which decrease relatively to the S vN .
along with it the idea of thermal metric fluctuations at the horizon, albeit small, as seen in figure 2. This also revisits the classical assumptions of a featureless background behind Hawking's initial calculation. Nevertheless, since the original formulation of the paradox, many have revisited the classical assumption of a purely classical the background metric (e.g. [ These various works have described black hole horizons as a thermal microscopic ensemble in order to recover correlations. Assuming a pure state, the corresponding off-diagonal elements can be seen as quantum fluctuations [16] and we will assume them to be present in rest of this work, although they have been subject of debate [23,40]. Indeed, similarly to a classical state, a quantum many-body state also exhibits quantum fluctuations in addition to classical fluctuations at late times, when the system has "equilibrated". This is defined as the moment the average local entanglement entropy has reached its maximum, or scrambling time, as illustrated in figure 2. The classical fluctuations can be seen as arising from the distribution of diagonal elements (representing probabilities) whereas quantum fluctuations can be seen as arising from the off-diagonal ones. These fluctuations appear due to interactions of the smaller microscopic degrees of freedom of the system. They can be used to recover the initial state in a closed system and therefore contain information. In the context of metrics, fluctuations are called "hair", and have been studied extensively [41,42,43]. An important classical result in "vanilla" general relativity is that asymptotically in time "black holes have no hair" [44,45,46]. Black hole hair appears therefore purely in models with quantum features.
As we have established quantum and classical fluctuations in an equilibrated system carry classical and quantum information (often referred to as "coarse grained" and "fine grained" respectively), we may study how to recover this information. The more entangled the final state, the more measurements and operations are necessary to recover quantum information (see e.g. [47,48], although the algorithm may depend on the microscopic physics). After the scrambling time any algorithm to recover the information from fluctuations requires a larger number of steps, contrary to the initial state where the information is more readily available. Furthermore the size of fluctuations in the system reach their minimum after this time, and their size can be well modelled using the Eigenstate Thermalization Hypothesis (ETH), which has experimental and numerical success with a variety of generic systems [16]. The ETH describes local operators aŝ where (f (Ē), g(Ē, ω)) are smooth functions of (Ē, ω) = ((E i − E j )/2, E i − E j ) and R ij is a matrix of random variables with zero mean and order one variance. Quantum fluctuations are therefore highly suppressed, and the classical fluctuations will dominate.

Black hole quantum information in gravitational waves
In our case, we are interested in fluctuations before the scrambling time, i.e. at times near when a particle fell in, since they are larger and depend more heavily on the disentangled initial state |infalling particle ⊗ |black hole . In the absence of a theory of quantum gravity, these can be calculated as gravitational waves using classical perturbation theory of a black hole metric with a source, which is exempt of the assumptions of the no-hair theorem ‡. Corrections due to fluctuations of the metric from the scrambled past history of the black hole will still be present, but we may reasonably expect them to be smaller in a variety of scenarios. Indeed, assuming the thermodynamic limit of many microscopic gravitational degrees of freedom (quanta) in thermal equilibrium with temperature T BH = T Hawking , the energy fluctuations of the system can simply be estimated as the standard deviation of a Bose-Einstein condensate, which is O(kT BH = M 2 p /(8πM BH )). Meanwhile, the signal energy of early fluctuations of a radially infalling particle has been found to be |u(ω, r)| 2 dωdr ≈ 10 −2 m 2 0 /M BH m 0 [49,50] where u(ω, r) is the GW signal while m 0 and M BH are the particle and black hole masses respectively (similar dependence on m 0 and M BH appear for more general trajectories). The recovered information will therefore have an SNR 2 (∝ χ 2 ) ≈ m 2 0 /M 2 p where M p is the Planck mass. Assuming the black hole and gravitational wave system is thermal such that information is spread evenly throughout the system, the quantum information recovered in this (sub)system can be quantified as the change in the Von Neuman entropy when taking the partial trace of the black hole density matrix over the subsystem such that: · · · |BH n BH n | · · · . . .
where (i, j) and (k, l) run over the Hilbert space of the black hole and the associated gravitational waves, in the columns and rows respectively, and whose states include all possible histories of the black hole given its initial conserved charges of mass, angular momentum, and possible gauge charges. Note this makes the simplest possible assumptions about ρ BH⊗GW . A more accurate estimate would require study of the thermalization and information loss of open quantum systems. Further, although the information recovered by the gravitational wave spectrum will be recoverable up to small quantum fluctuations of the metric, the information in these fluctuations can also be retrieved in principle by measuring them. Assuming Hawking radiation arises from the decay of microscopic gravitational degrees of freedom at the horizon, one expects it contains information that is approximately spread out through the horizon and proportional to the size of the fluctuations which is given by the temperature. These fluctuations represent the late time tail of the black hole history ‡ We will focus on m 0 M BH and use perturbation theory to allow us to show analytically that information is present as we see in the next section, although we believe these conclusions hold also in the strong gravity regime (eg black hole mergers).  . Illustration of the information content of the metric. The black hole metric is assumed to be a thermal bath of microscopic degrees of freedom. As the bath evaporates through (red) blackbody Hawking quanta, the energy stored in the gravitational field is released, causing fluctuations (black) outside the horizon set by the blackbody temperature corresponding to noise in the measurement of the classical gravitational waves from the infalling particle. During the infall the system evolves from the product state |BH I ⊗ |m 0 J to become part of the entangled thermal state i,j |BH i |m 0 j . [26,27,28,30,33,34,35,37,38,51,52]. We may therefore think of information retrieval in two parts: the first is reading large fluctuations in the form of gravitational waves from which information is easily retrievable, and the second is reading the late time tail fluctuations occurring after the scrambling time which encode the black hole history. Again note this tail of quantum and classical fluctuations appears as a consequence of the quantum features attributed to the horizon, since classical metric equations predict exponentially small signal with time in line with the no hair theorem. Further, the slow escape of the information from the horizon during the late time tail is in sharp contrast with early phase gravitational waves, where the infalling particle is still mostly unentangled and coherently scatters with the microscopic degrees of freedom of the black hole, allowing for easy information recovery away from the horizon.
We may also venture to estimate the share of the total information contained in the gravitational wave signal, setting aside for now the difficulty of recovering the information. Assuming the Von Neumann entropy of a small subsystem grows linearly with the energy (which is an extensive quantity), we may estimate the information content of gravitational waves as their share of the total energy of the BH and GW system. This corresponds to taking the entanglement entropy of the black hole after formation but without the gravitational waves. It can be expressed as the average share of energy m 0 emitted as GWs, which we call E GW , during infall: Information in GWs Information in BH Equations Classical picture Table 1. Summary of share of information estimates of the initial collapsing state in GWs in both the fully classical picture (eternal black hole) and quantum case (including Hawking radiation and backreactions).
obtained over the black hole history, assuming most of the signal is recoverable (SNR > 1). However, this fails to account for the potentially large share of black hole history information present in the gravitational waves, which reduces/collapses the Hilbert space further. As we will see in section 4 this is O(1), such that the share of information of an observable can be estimated as O(1) × max(SNR, 1). The total entropy in the system S cl is therefore reduced by a factor of E fluct /(E GW + E fluct ) = 1/(1 + SNR 2 ), and the fraction of the information in GWs becomes 1 − 1/(1 + SNR 2 ) = SNR 2 /(1 + SNR 2 ), which is large for information regarding any constituent of the black hole larger that the Planck mass. Thus gravitational waves contain the black hole information down to a Planck mass granularity, and the share of black hole information may be found as which is quite large for astrophysical black holes. Note we assume an O(1) share of the information is obtained for SNR> 1, although a more conservative estimate would use SNR> 10, and depends on the experimental setup [53]. However the main purpose of this work is to motivate information retrieval using gravitational waves and we leave this for future work.
The scope of black hole information retrieval in gravitational waves in different regimes is summarized in Table 1. In the classical picture ( → 0 such that M p → 0), black holes are eternal and information is lost behind the horizon. However the information of the initial state is also copied and reflected near the horizon in the form of gravitational waves, preserving the initial state information. In the semiclassical picture, where fields are quantized on a smooth metric background, black holes evaporate into maximally mixed states, containing no information besides the conservation of total energy [13]. In this context, gravitational waves provide the only source of information. Finally, in a purely quantum regime, one expects fluctuations of the black hole metric to contain information along with Hawking radiation. However by using only the GW signals (early fluctuations) to extract information we obtain the same share of information as in the semi-classical case ("quantum" picture in Table 1), i.e. a sizeable portion of the initial state information, which converges to the classical result in the limit M p → 0.

Information recovery in the gravitational wave literature
With this estimate we are left to quantify the other advantage of gravitational waves: the complexity, or relative simplicity, of recovering information from gravitational wave signals, which we do for black hole-particle mergers in the following sections. First however we comment on their role within the current body of literature on gravitational waves.
Gravitational wave signals from black hole perturbations have been studied extensively (e.g. [17,18,49]), and obtaining them to good numerical/computational accuracy is non trivial. However as we will see, once a few key quantities which can be computed in advance are known the dependence on the various quantum/classical numbers of the infalling particles and initial black hole is trivial, such that with full knowledge of the signal (i.e. in the absence of noise) one can determine them exactly. This is in fact a common implicit assumption behind many gravitational wave detection works [53,54,55,56]. However remarkably it has not been explicitly shown before, although mentioned in passing in [57], partly due to its irrelevance in the face of traditional experimental constraints such as intrinsically noisy detectors, irreducible background signals and numerical accuracy. The information is therefore usually determined with a Bayesian approach using catalogues of waveforms (e.g. [58,59,60]) such that the widths of the significance contours go linearly with these uncertainties: where (θ S , θ N ) represent vectors of parameters describing the signal and noise respectively and d is the detector data. Nevertheless in an experimentally finite SNR regime, we may recast our previous statement regarding the accuracy of an observable by instead claiming that Bayesian contours do not give disjoint regions of high statistical significance when in the limit of highly sensitive/low noise experiments.

Zerilli perturbation equation
We will consider here a Schwarzschild black hole formed from particles radially falling in, perturbing the metric and emitting gravitational radiation. Other types of perturbation have been considered, and many are solvable to first order [61]. Schwarzschild black hole perturbations are separable into ten scalar degrees of freedom using a tensorial harmonic basis, which can be reduced to four polar (P-even) and two axial (P-odd) degrees of freedom by an appropriate choice of gauge [62]. We may then solve the Einstein field equations of each set separately, resulting in the Zerilli and Regge-Wheeler equations [62]. Furthermore Chandrasekhar proved these equations were equivalent (isospectral) [63] through a supersymmetry transformation. Chandrasekhar also proved [64] that the limit of the Teukolsky equations at zero spin, where Kerr spacetimes become Schwarzschild, is equivalent to the Regge-Wheeler and Zerilli equations, as expected. Thus we expect our method to be generalizable to Kerr perturbations. We also expect the method to generalize to Reissner-Nortstrom black holes as it reduces to a Schwarzschild metric in the limit of zero charge, and since perturbations on this background follow at Zerilli-like equation. The Zerilli equation is obtained by perturbing a background Schwarzschild metric g µν = diag(−f (r), 1/f (r), r 2 , r 2 sin 2 θ), where f (r) = (1 − 2M/r) and use M = M BH for notational brevity, with a small parity-even perturbation h even µν such that g µν → g µν + h even µν and solving the field equations. Writing where Y lm are spherical harmonics and r * = r + 2M ln(r/2M − 1) is the tortoise coordinate §, one can express the resulting equations in terms of just two out of the initial four scalar degrees of freedom, which are chosen to be H 1 and K such that we discard H 0 & H 2 . With the following choice of variablê one can write the radial Zerilli equation [65]: where with η = (l − 1)(l + 2)/2 and as usual (M, ω, l) are the black hole mass, frequency (eigenvalue of the Zerilli equation) and spherical harmonic l-number. Therefore as expected from the e.g. de Donder gauge, we see perturbations of a Schwarzschild metric (gravitational waves) are solutions of a wave equation with admit two real degrees of freedom: polarizations K and H 1 . In what follows we will assume that there is a oneto-one mapping between experimental GW data, H 1 and K, andû. We will further assume a gedanken experiment which allows us to obtain GW signals andû from any point x µ and with arbitrary accuracy.

Source term
The right hand side of 12 describes the source of the perturbations: where m 0 is the infalling particle mass m 0 M . The first sum term of I is called the "source term" denoted S whereas the second term represents perturbations present at t 0 when we start observing GWs, and appears in the frequency domain as an artifact of the Laplace transformation of the time domain Zerilli equation: where the Laplace transform is: For now we may ignore the second term in I(ω, r) if the signal is detected such that u ≈ 0 for t 0 → −∞. T (r) describes the coordinate time at which the particle can be observed at radius r, which for radial infall from rest is [49] T (r) (17) While the expression for T (r) and hence I(ω, r) appear cumbersome, remarkably m 0 only appears as an overall factor in I(ω, r), which allows us to extract information about the falling mass relatively easily, as we will see in the next section.

Solution to the wave equation
One may solve the Zerilli equation 12 using a Green's function such that: whereû 1 (respectivelyû 2 ) is a solution to the homogeneous Zerilli equation, i.e. I = 0, that respect the desired boundary condition at the horizon (respectively infinity) and W = ∂û 1 ∂r * û 2 − ∂û 2 ∂r * û 1 . The solution to 12 is then: The Green's function is a solution to the Zerilli equation with a point-like source function I located at ζ. Since solutions to the Zerilli equation are linear in the perturbations, one may sum over a distribution of point-like potentials. The solution to 19 is then the weighting of the Green's function against a distribution I, which is a plane wave with time dependent phase for radially infalling point particles.
More specifically, the r * → ±∞ boundary conditions on theû i (ω, r * ) solutions allow us to write in a convenient plane wave basis: (20) and lim r * →+∞û 2± (ω, r * ) = e ±iωr * , lim We define quasinormal modes (QNMs) ω nl as having A in (ω nl ) = 0. Note that we define twoû 2 functions for later use but will only useû 2 =û 2+ in the Green's function since we assume no incoming radiation from infinity . From 20 & 21 we find that W = 2iωA in (ω). We can now write the solution in ω-space: If the source has compact support near the horizon and the experiment, located at r * , is far from the horizon and the source, this simplifies to: This may further simplify if the source is sufficiently close to the horizon to replacê u 1 ∼ e −iωζ . Although this may not be our case, in what follows we will focus on this integral for simplicity, since the result and techniques can be straightforwardly generalized to the case of 22. At first glance the integral in 23 is typically divergent at the horizon. Nevertheless we may regularize it by imposing appropriate boundary conditions (see Appendix).
We now want to understand the relation with the time series which is measured at our gedanken experiment. We must therefore return to the time domain via an inverse Laplace transform. This integral is non-trivial and requires defining a contour in the complex ω-plane (see figure 4). The contour contains poles at the quasinormal frequencies ω nl where A in (ω nl ) = 0, which following the residue theorem contribute to the result: Also note normalizing the solutions and taking the r * → −∞ limit gives B out (ω) = A in (ω) and Often a branch cut (hashed region) is introduced for values with ω R = 0 which don't oscillate and thus aren't QNMs (non-radiative), however we include them in our work without loss of generality (24). Their physical significance still requires careful study, but may represent higher order contributions to the signal or late time tail effects, since both contribute late-time tails (see e.g. [61] and references therein).
Note the inverse derivative of A in is an intrinsic property of the black hole metric, and combined with the source integral is often called "quasinormal excitation coefficient" (QNEC). The observed radiation time series is thus a sum over QNMs ω nl and a continuous spectrum of half circle modes and the branch cut. We discuss in the following sections the effect of the contour on information recovery.

Finding the QNMs and A in 's
Although the A out 's cancel, the residue coefficient requires knowledge of A in in the neighborhood of ω nl . We follow [61,66] in determining these using the following series ansatz in ω-space and Schwarzschild coordinates with 2M = 1 ¶: where G j+ν and F j+ν are the regular and irregular Coulomb wave functions with ν chosen to yield a minimal solution for the series b j (i.e. with a fixed initial value a 0 , which 20 fixes to b 0 = 1) and ¶ A simpler analytical approach is possible from [61] in the limit of large n.
These expansions respect the correct boundary conditions e.g.
Although the above expansion holds for all r * , the series a j diverges as r → ∞, and is therefore impractical for direct evaluation of A out (ω nl ). In 20 & 21, we have defined three solution to the homogeneous Zerilli equation, namely u 1 , u 2− and u 2+ , each respecting a different boundary condition. Any pair of these three solutions will form a basis of the space of solutions to the homogeneous Zerilli equation, which is a second order linear ODE. Hence, any one of these solutions can be expressed as a linear superposition of the other two (the combination is oftentimes referred to as a Bogoliubov transformation). The coefficient of the Bogoliubov transformation can be found from the boundary conditions that defines the various solutions, in particular we can write Thus we can solve for A out and A in for all ω and some chosen finite r * + : from which we can calculate dA in /dω nl once we obtain (û 1 ,û 2± ). To do this we plug the series into the Zerilli equation. The recursion relation obeyed by the coefficient of the series is given by (re-establishing mass units where 2M = 1): α 0 a 1 + β 0 a 0 = 0 α nl a n+1 + β nl a n + γ nl a n−1 = 0 (31) where n corresponds to the series index. Building on the three term recurrence relation, Leaver showed that the solution to the Zerilli equation will converge [68,66,69] if the (a n ) n∈N is a minimal solution to the above recurrence relation. The existence of a minimal solution implies that the following continued fraction holds: where we have used the standard notation for continued fraction in the above equation. 33 will be valid for a discrete set of complex frequencies which are the QNM frequencies.
The index labeling the element of this set is called the "overtone" number. In practice one numerically solves 33 or its n th inversion to find the n th QNM frequencies.
This method of finding QNM frequencies, commonly called Leaver's method has been improved and generalised to a wide range of space-time beyond the Schwarzschild metric considered here. QNM solutions (A in (ω) = 0) occur at fixed discrete solutions of u i (ω, r * ) and {α nl , β nl , γ nl } n∈N . By inspection of 32, the coefficients remain constant if ω nl is simply related to the mass: The proportionality coefficients are determined solely by the structure of the black hole perturbation equations and (n, l) and can be determined up to arbitrary accuracy (i.e. given sufficient computational resources) by inverting the continued fraction. For the "fundamental" frequency (n, l) = (0, 2) they are often quoted as (12.07, 18.06) kHz M . Having explained how to calculate each of the quantities of 24 our gedanken GW experiment may find the mass of the black hole and the rest of the parameters with arbitrary accuracy, as we will see in the following section.

Inversion
We are now ready to invert the GW signal to retrieve physical information (semi-) analytically. As mentioned previously, the GW signal will contain information regarding the initial black hole and perturbing infalling mass. To our knowledge the method for extracting the information in the case of a gedanken experiment with arbitrarily high SNR and spatial access to the signal outside the horizon has not been shown explicitly before. In practical cases the detector will be noisy and cannot be easily placed at several points around the black hole or take data for arbitrarily long times with arbitrarily high spatial and temporal resolution. The signal noise will then propagate to the amplitudes in the Laplace transform, and subsequent parameters deduced from them, and justifies instead the Bayesian approach used by gravitational wave astronomers. However it is of classical information-theoretical importance to know whether there exist unresolvable degeneracies in the signal, which next generation GW detectors might encounter in the near future.
To invert the signal we assume in this subsection no a priori knowledge of the mass of the black hole or the perturbing particle, other than the assumption of the black hole being Schwarzschild, although generalization to Kerr metrics should be possible and is reserved for future work. We wish to reconstruct both m 0 and M BH to arbitrary precision, assuming arbitrary precision of our GW data can first be obtained (i.e. SNR→ ∞).
Before t = t 0 we must obtain r * , the location of the detector with respect to the singularity, as it will soon become useful. During this time one can match the readings of two or more detectors around the black hole whose measurements are found to coincide when rotated, assuming each detector is of finite non-zero size and orientable. By taking lines aligned with these detectors and which pass through the detectors' centers one finds the singularity as the point the lines intersect. Alternatively, for idealized point-like nonorientable detectors we may find the distance to the singularity by finding the equation of sphere of four points with the same signal using e.g. [70].
When measuring the signal at t ≥ t 0 , assuming the signal comes from a mass infalling at angle (θ, φ) and that the GW detector acts on a sphere of arbitrary radius which we assume near the horizon, this angle must first be found (as shown in Section 6) to divide out the spherical harmonic from which the radial component u(r * , t) can be deduced. The radial signal u(r * , t) comes in the form of 24, where we have measured the half circle modes as well. We then may Laplace transform the signal to obtain the frequency space representationû(ω, r * ) 23.
Next, we observe the ω's at which the divergence of 23 occurs: these are our QNMs frequencies ω nl . The frequency with the lowest imaginary part ω I , or "fundamental" frequency, immediately yields the black hole mass, as discussed around 34. To get the mass m 0 , we first remove the radial phase e iω nl r * in the complex mode amplitudes with knowledge of the detector's distance from the horizon r * and the QNM values ω nl . With the black hole mass in hand, we may also calculate the coefficients (α nl , β nl , γ nl ), which depend on M BH once physical units are reintroduced, * and in turn (û 1 (ω, r * ),û 2 (ω, r * )) according to 25 & 31. Using (û 1 (ω, r * ),û 2 (ω, r * )) we may then determine A in (ω) and dA in /dω nl using 30.
We remind ourselves of the expression for the source, i.e. 14 and ignoring the boundary term. One may now perform the integral in 23: using the regularization technique described in the Appendix, such that it converges for all ω. The source term being simply proportional to m 0 , the mass of the infalling particle we wish to recover, we may write: where I 0 = I/m 0 . Note that this equation is valid for QNMs as a limit ω → ω nl . We require QNMs as the series u i (ω, r * ) do not converge or respect the boundary conditions at other frequencies [68].
We have now peeled away each layer of 24 semi-analytically, revealing the quantities (M BH , m 0 ) determined to arbitrary accuracy assuming the GW signal from the particle's infall could be detected with arbitrary accuracy. Quantities which don't depend on m 0 such as A in (ω) or 35 may be evaluated numerically using M BH found from the signal u(ω, r * ). As we will see in the following section we may also extract the time and angle of infall. * Note here we are talking about (α nl , β nl , γ nl ) for all ω, not just minimal solutions.

Time recovery
In this section we briefly discuss determining the time of the infall which we will call t 0 . Ignoring any angular momentum in the history of the black hole means the classical information should consist entirely of the masses, times, and angles from which the particles radially fell in. Having seen in the previous section how to obtain the masses m 0 and M BH , we turn our attention to the time.
Infall time t 0 may be defined with reference to any point during the infall, since the observer never sees the particle fall in and we are not provided a natural reference time t 0 . However we may choose an arbitrary radial coordinate r * 0 at which we wish to know the time t 0 when the particle crossed by looking only at the GW signal. We start by assuming a Schwarzschild reference frame with an arbitrary moment t = 0, where t is the coordinate time. t 0 is then the coordinate time at which the particle is considered to fall, or be dropped, into to the black hole. Following 19 we assume a radial infall from infinity with zero energy. We assume without loss of generality the observation window includes a period where the particle has been en route to the black hole for enough time to start creating a GW signal larger than the arbitrarily small (but larger than zero) noise of the experiment, which we may choose to be at the level of the metric fluctuations. Although we do not need to observe the entire signal to obtain information, we should always try to include the period of largest GW signal in the observation window such that the quantity/quality of information is optimized.
The GW spectrum in ω-space is characterised by 23. We recall that the only difference in the particle infalling earlier or later is an overall time translation, which we call t 0 , such that "earlier" → e iω nl (t−r * ) "later" → e iω nl (t+t 0 −r * ) .
This e iω nl t 0 phase is carried over toû(ω, r * ), notably in 36 where we had implicitly assumed t 0 = 0. Generally this introduces a phase, which may be recovered by finding t 0 is some frequency up to a periodic degeneracy where A nl represents the computed quantity 36 and n ∈ Z, and which is bounded above by the inverse of the experiment's observation period 2π/ω min . We may however hope improve the degeneracy by obtaining A nl for two frequencies and observing that is solved when (m ω 1 )/ω 0 is an integer, which may require (n, m) 1 or may have only one solution m = n = 0 depending on (ω 0 , ω 1 ) and would either reduce the degeneracy by increasing it to larger periods (n, m) × 2π/ω, or remove it entirely by fixing (n, m) to zero. In a gedanken experiment we may assume the observation time is infinite, thus avoiding the degeneracy in 39. However since we are not guaranteed QNMs with arbitrary small ω R 1, one could search among QNMs for QNM ratios which require large multiplicative factors to become integers and reduce or eliminate the degeneracy. We leave to future work the study of the ratios of QNMs.
Concerning our resolution when determining t 0 , typically in an experiment this is given by the period of the highest accessible frequency, since any error due to a finite SNR in estimating t 0 becomes more visible at higher modes: where we set t 0 = 0 for convenience and expanded to O( ) 1. However since we expect QNEC(ω) decreases faster than ω −1 from the finiteness of the total GW energy, increasing ω does not yield the desired effect of enhancing the error to make it more visible. As discussed in Section 1, since the maximum SNR is fixed due to background fluctuations, we may in principle calculate a time resolution and share of temporal black hole information accessible in the signal, similarly to our estimate for the mass m 0 . Finally, this argument also applies to the angular resolution of the black hole history, as we will see in the next section.

Angle recovery
We now wish to recover the angle at which a particle fell radially into the black hole. We remind ourselves that the gravitational perturbations of a spherically symmetric background metric can be separated as such: where u nlm (r) is expressed in 22 and we assume the integration half-circle radius r HC → ∞. This is because the full differential operator L = L Zerilli + L Legendre is separable into a radial and angular piece. We assume data from the GW wave signal of a radially infalling particle at (θ 1 , φ 1 ). We do not assume a reference frame with θ 1 = 0 since we do not know θ 1 . Instead we assume an arbitrary reference frame (θ 0 , φ 0 ) at some fixed radius . A radially infalling particle will create only m = 0 modes in a frame with θ 1 = 0, since the differential equation Lu = S is symmetric around the axis passing through the singularity and the particle and |m| > 0 modes are not. We may then write the solution 24, assuming without loss of generality r HC → ∞ and writing ω nlm → ω nl for simplicity, in the form: We may think of the more complicated case without varying radius without loss of generality. We may also think of partial angular data of the GW. In this case we expect bandwidth considerations analogous to the time resolution (see previous section) to appear.
where (θ , φ ) are the coordinates in the (θ 1 , φ 1 ) frame, and unprimed coordinates mean we are in the (θ 0 , φ 0 ) frame. Note that while we express angles of (θ 0 , φ 0 ) and (θ 1 , φ 1 ) reference frames inside a single variable, these angles admit different representations in different frames which must be accounted for in computations, as we will see below.
Nevertheless it may prove convenient to express the angles in various frames, and thus we use the prime/unprimed convention.
We now wish to recover (θ 1 , φ 1 ). To do so, we may solve (minimize) for (α, β) such that: is maximized when integrated over the sphere, as we will see below. Note we dropped the sums over n which factors out, and we wish to integrate (θ , φ ). This allows us to maximize Equation 44 for every l ∈ N because spherical harmonics are normalized ( Y Y * dΩ = 1) when Y and Y * are identically oriented and centered. Here the second harmonics (unprimed) angles expressed in the (θ 0 , φ 0 ) reference frame are (θ + θ 2 , φ + φ 2 ) = (θ + θ 1 + α, φ + φ 1 + β), where (θ 2 , φ 2 ) is the guess for the angle at which the particle is falling. Using Wigner's D-matrix we may write the expression as: where we used m D lm m (θ i , φ i )D lmn (θ j , φ j ) = D lm n (θ i + θ j , φ i + φ j ) (closure of the D-matrix), m D * mm (θ, φ)D mn (θ, φ) = δ m ,n (orthogonality of the D-matrix) and D l00 (θ, φ, 0) = P l (cos θ) where P l is the l th Legendre polynomial. The expression then gives l P l (cos θ (α, β, θ 1 , φ 1 )) = l P l (cos θ 1 cos(θ 1 + α) + sin θ 1 sin(θ 1 + α) cos β) where we used the spherical law of cosines in the first line then expanded to O(α 2 , β 2 ), and where (l min , l max ) are the lowest and highest angular resolution (i.e. angular bandwidth) of the experiment. We see that, for any finite signal-to-noise ratio and (α, β) > 0, l max needs to be large enough to reject incorrect choice of (θ 2 , φ 2 ), and for an arbitrarily large SNR and (β, α) 1, we need arbitrarily large l max . This is analogous to the determination of the time of infall t 0 for a single mode with ω max → l max .
Regarding angle periodic degeneracies in determining (θ 1 , φ 1 ), we note although it is 2π/l min for |m| > 0, remarkably for m = 0 no degeneracy is present. Finally, for multi-particle sources, preliminary separation of their signals is necessary before finding the angle of infall of each source and can be performed using highly-damped modes, as we will see in the next section.

Resolving multiple particles and/or matter distributions
We have now determined the mass, time and angle of an infalling massive particle. However there remains the question whether two particles, or more generally a distribution of matter, can create degeneracies in the signal or can be distinguished.
Naively the sum of two two signals can be disentangled if one can solve for the masses and infall times of each particle (m 0 , m 1 , t 0 , t 1 ): where A nl represents 36 (obtained from observation). In principle one may solve for the four variables using only two QNMs, thus forming four real equations (more generally, one needs N QNMs to solve for the time and mass of N sources), however the algebra is non-trivial as can be seen in the sum of two modes. Instead we may look at the set of QNMs to separate these equations using mode damping. Assuming an arbitrary desired accuracy in determining m 0 , one may choose ω nl with large enough ω I (i.e. overtone n 1) such that, assuming t 1 > t 0 without loss of generality, the term m 1 e iω nl (t 1 −t 0 ) in the first line is suppressed by a factor of e −ω I nl (t 1 −t 0 ) , where ω I nl is taken as positive, such that one can extract m 0 with an error of O(m 1 e −ω I nl (t 1 −t 0 ) ). Indeed, for any l, ω R nl will converge to a constant while ω I nl diverges as n → ∞, such that there is no limit to the arbitrary accuracy we demand. The closer the infall times t 1 and t 0 are, the larger n will need to be to meet this criteria, and increasingly high n-modes will contain the information about increasingly smaller distances between particles. Assuming N > 2 we may split the sources/particles into two groups of nearest neighbor particles such that all the particles of the first group fall in earlier than those in the second. We may then distinguish between the signal from each groups similarly to the case of two particles, and repeat the process of dividing each group into its components. Thus our reasoning generalizes to more complex structures, including classical objects formed of individual quanta, at which point a prescription for quantizing the gravitational effects would be necessary.
Looking at more general stress-energy tensors, this resolution should to apply to lightlike (massless) particles as well, although the corresponding signals are solved using different geodesics which we do not consider here, as long as the difference in mode amplitudes the potential provokes can be compensated by a choice of QNM mode with sufficiently large exponential suppression. Also note signals sourced by combinations of massive and massless particles (the most realistic sources) may be disentangled using this approach, with appropriate choices of QNMs, but is left for future work. Finally, considering particles with same infall times, but different infall angles, we may solve the angular analogue of 47, calculated as sums of spherical harmonic contributions to the signal, to obtain the N angles of infall, and associated amplitude of N sources. Although we may not use the decay of modes as in the time domain to help us, numerical investigation has confirmed that the angular analogue of 47 admits a single solution at the correct values.

Conclusion
In this work we determined a complete classical information retrieval process using gravitational waves from a Schwarzschild black hole's history, and found that it is in principle possible to obtain to arbitrary accuracy in the absence of noise and with arbitrary computational resources. Indeed, from the GW frequencies one can recover the mass of the black hole, and the mode amplitudes contain the masses of the perturbers. Meanwhile, the infall times are contained in the phases of the amplitudes and can be separated using multiple modes, while the angle at which the particles fell in are contained in the angular phases of the GW signal modes. Therefore not only has all the information been recovered, but we have used every degree of freedom the signal contained contained.
Further, we may account for noise which is naturally present in the semiclassical regime as the black hole has a temperature. Nevertheless we show that the history of the mass of the black hole can be retrieved to O(1) when SNR> 1 for granularity above the Planck scale (similarly the share of angular and temporal information can be determined). The share of quantum information in GWs can therefore be quite large and can be estimated as the share of entanglement entropy between the black hole and GWs. Accounting for the measuring of the gravitational waves, the information recoverable in the black hole is then reduced from the initially large set of initial state possibilities with identical macroscopic number M BH to those consistent with the gravitational waves. Although the background fluctuations limiting our measurement have information as well, it is harder to retrieve/untangle and requires waiting until the end of the evaporation to obtain fully. This information is more akin to the "closed system" information of the black hole studied in e.g. [7,8,33]. Gravitational wave information is therefore a natural first place to look for information as it is more immediately available and relatively unentangled with other signals and the black hole.
We may put the importance of black hole information in gravitational waves in context by looking at three different regimes. In the purely classical regime ( → 0), black holes don't evaporate and information is considered lost behind the horizon and destroyed at the singularity. However initial state information is preserved in gravitational waves, where the information is copied and reflected at/near the horizon. In the semiclassical picture of quantized fields on a smooth metric background, black holes evaporate into maximally mixed states as shown by Hawking, containing no information besides the conservation of total energy of the initial state. In this context gravitational waves are the only source of information. Finally, in a purely quantum regime, one expects fluctuations of the black hole metric to contain information along with Hawking radiation. However a sizeable portion of the initial state information may still be in the gravitational waves emitted during the black hole history. Our finding of a Planck scale resolution is consistent with the expected quantum noise of gravitational waves below this scale, and implies Hawking radiation is induced by Planck scale fluctuations of the horizon.
Finally, looking ahead future work is needed in several directions, including finding the share of recoverable temporal and angular information, or generalizing our results to information recovery in rotating and electrically charged spacetimes. We may further wish to clarify the role of the branch cut modes, which cannot be described by QNMs, e.g. whether they can describe a late time perturbative increase in mass of the black hole. Also our work is performed to first order in perturbation theory, in principle a self consistent description of the geodesic and perturbations is possible to arbitrary order. We hope that this work can stimulate further studies of gravitational waves in these directions.