Density matrix of black hole radiation

Hawking's model of black hole evaporation is not unitary and leads to a mixed density matrix for the emitted radiation, while the Page model describes a unitary evaporation process in which the density matrix evolves from an almost thermal state to a pure state. We compare a recently proposed model of semiclassical black hole evaporation to the two established models. In particular, we study the density matrix of the outgoing radiation and determine how the magnitude of the off-diagonal corrections differs for the three frameworks. For Hawking's model, we find power-law corrections to the two-point functions that induce exponentially suppressed corrections to the off-diagonal elements of the full density matrix. This verifies that the Hawking result is correct to all orders in perturbation theory and also allows one to express the full density matrix in terms of the single-particle density matrix. We then consider the semiclassical theory for which the corrections, being non-perturbative from an effective field-theory perspective, are much less suppressed and grow monotonically in time. In this case, the R\'enyi entropy for the outgoing radiation is shown to grow linearly at early times; but this growth slows down and the entropy eventually starts to decrease at the Page time. In addition to comparing models, we emphasize the distinction between the state of the radiation emitted from a black hole, which is highly quantum, and that of the radiation emitted from a typical classical black body at the same temperature.


Introduction
There has been a recent spike in activity on understanding the implications of black hole (BH) evaporation [1,2]. This can be attributed, in large part, to the controversial proposal that a unitary evaporation process comes with the cost of a "firewall" [3] or "energetic curtain" [4]-these being colorful euphemisms for an apparent tension between general relativity and the unitarity of quantum field theory in a BH setting [5][6][7][8]. While the debate rages on, a consensual mechanism for information release is still lacking.
In Hawking's model, the process of BH evaporation is not unitary [2]. Hawking argued that the correlation functions of the emitted radiation are diagonal in modeoccupancy number, frequency and emission time and, from this, deduced that the density matrix for the radiation is diagonal in the very same quantities. In fact, to good accuracy, the final state of the emitted radiation is thermal, and so it is maximally mixed. Then the evaporation cannot be the result of unitary evolution from a nearly pure state, which lead Hawking to state (a statement which he more recently retracted [9]) that the process of gravitational collapse is not compatible with the standard principles of quantum mechanics. The established benchmark model for describing unitary BH evaporation has been the Page model [10]. In this model, the BH and the emitted radiation are assumed to be in a pure state in some large Hilbert space, which is partitioned into an "in" part representing the BH and an "out" part representing the outgoing radiation. Because the combined state is pure, there must be a special basis in which the full density matrix has only a single eigenvalue. An external observer measures the reduced "out" density matrix in some random basis and treats the transformation matrix U as a random matrix with prescribed statistics [10,11]. (See also [12] for a recent review.) When the "out" system is smaller than the "in" system, it looks to an external observer as if it were thermal. Conversely, when the "out" system becomes the larger one, the deviations from a thermal state grow and indicate that the total state is indeed pure. The critical time in which the midpoint of evaporation is reached is normally called the Page time. Given that the BH is initially in a pure state and that the evolution is unitary, the Page model can plausibly be viewed as setting the minimal rate for purification of the emitted radiation.
Two of the current authors (RB and AJMM) have recently proposed a semiclassical model of BH evaporation [13,14], with the premise of repeating Hawking's seminal calculations [1,2]-which assume a classical background metric-so as to include a fluctuating BH geometry. The incipient BH is endowed with a quantum wavefunction [15,16], leading to expectation values in place of fixed classical parameters. From this perspective, Hawking's model pertains to the limit of an infinitely massive BH with a fixed size, whereas the semiclassical reformulation treats the BH mass as finite with a continually decreasing size due to classical back-reaction effects. In the new model, the correlation functions of the emitted radiation are no longer diagonal and the evaporation process becomes unitary, even though the thermal-like emission spectrum is maintained. This model of unitary BH evaporation obeys the general constraints about the rate of purification. In particular, it has been shown in [17], relying on the results of the current paper, that the Rényi entropy of the radiation decreases monotonically from its peak value at the Page time and does so at a faster rate than the Page model predicts.
From the point of view of an effective field theory in a fixed curved-space background, the correlation functions of the semiclassical model are modified in a nonperturbative way [16]. Contrary to expectations, these modifications could become significant in some situations. It has been shown that, when the leading result in the effective field theory either vanishes or diverges, then the non-perturbative corrections are particularly relevant [18]. In the first case, we expect a correction of the order 1/S BH (the inverse of the BH entropy) to become the leading result and, in the second, the divergence should be replaced by a large but finite S BH . These small non-perturbative quantum effects could be coherent, and so their amplitudes could add up and grow with time; even growing to the point where they become comparable to the classical outcome [13].
The previous studies of the new model have utilized the single-particle density matrix (the two-point function) of the radiation as the primary tool 1 rather than the full multi-particle density matrix. The single-particle density matrix is easier to calculate and, for the radiation emitted by a BH, sufficient to completely determine the full density matrix. One of the primary goals of the current paper is to make explicit the relationship between the single-particle and multi-particle density matrices. The results are presented in a way that can be applied to the Hawking model, the proposed semiclassical model and even to more general modifications of Hawking's framework. As an upshot, we are able to find a closed expression for the Rényi entropy of the full density matrix that is expressed directly in terms of the single-particle density matrix.
It has been argued for BHs in Anti-de Sitter (AdS) space-at first by Maldacena [19] with subsequent elaborations by many others -that non-perturbative contributions from other geometries which do not possess a horizon could be relevant to the state of the Hawking radiation. Maldacena's specific example was the contribution of thermal AdS to the matter correlations for large BHs. In this case, the contribution to the density matrix of the radiation is expected to be exponentially suppressed ∼ e −S BH . However, because the number of off-diagonal elements is exponentially large ∼ e +S BH , one could not be sure as to the relevance of such a contribution. It could then be argued that, without a better control on the non-perturbative contributions from other geometries (which are usually not available), it becomes difficult to trust the Hawking result.
We have two comments on this issue. First, according to our analysis, one requires much larger corrections ∼ e −S BH /2 to challenge the reliability of the leadingorder outcomes (See, for example, [20].) Second, the quantum fluctuations of the background BH geometry induce off-diagonal effects in frequency space that are only power suppressed ∼ 1/ √ S BH and sufficient by themselves to restore unitarity. These statements will be made more precise in due course. The rest of the paper is organized as follows: In the next section, we review the basic framework of the Hawking model of BH evaporation and recall its main results. Then, in Section 3, we formulate a precise relationship between the multiparticle density matrix and the single-particle density matrix. The preceding analysis is applied to the Hawking model in Section 4, with the aim of regularizing formal divergences that arise in the eternal-BH framework. Section 4 also includes a discussion on how the state of the BH radiation differs from that of thermal radiation emitted by a "normal" black body. The formal analysis up to this point sets the stage for Section 5, where we discuss similarities and differences between the three models of BH evaporation: the Hawking model, the Page model and the semiclassical model. There is then a brief conclusion in Section 6, followed by four appendices that fill in some technical gaps in the main text.

The Hawking model of black hole evaporation
This paper uses the notations of Hawking in [2]. In particular, we consider the initial vacuum state |0 − and a final vacuum state |0 + . The initial vacuum contains no particles with respect to the creation and annihilation operators {a + i , a i } as defined at past null infinity, Following Hawking, the Hilbert space at future infinity is a direct product of the Hilbert space of the particles falling into the BH, H in , and of those going out, H out , so that H = H in ⊗ H out . The ingoing particles are created and annihilated by the operators {c + i , c i }, whereas the outgoing ones are created and annihilated by Both are related to {a + i , a i } via Bogolyubov transformations. The final vacuum state can be represented as |0 + = |0 in |0 out and contains neither outgoing nor ingoing particles, The two vacuum states |0 − and |0 + do not coincide with one another due to the particle creation by the BH. In other words, an initially empty state will appear to have a non-zero occupation number in the final state, meaning that b i |0 − = 0. The initial vacuum can be expressed as a linear combination of the Fock states at future infinity as where |a out and |b in are the basis vectors of the Fock spaces which are spanned by the operators {b + i , b i } and {c + i , c i }, respectively. As a consequence, The object of interest is the density matrix of the initial vacuum state, which allows one to calculate the vacuum expectation values of observables according to the standard rule 0 − | O|0 − = tr ρ vac O.
In particular, we will be interested in the expectation values of the observables at future infinity. These will be composed only from the outgoing creation and annihilation operators {b + i , b i } and so can be written as O = O out ⊗ I in , where I in denotes the identity operator in the Hilbert space of the ingoing states. The expectation value of this operator then takes the form Here, ρ out ac are the matrix elements of the reduced density matrix of the Hilbert space H out for the outgoing radiation and are determined by The elements of the density matrix for the outgoing radiation can be found by calculating different moments of the operators {b + i , b i }. In particular, the number operator of the jth outgoing mode is The famous result of Hawking is that these elements are the same as for thermal radiation, where Γ is the grey-body factor. The density matrix (2.10) is diagonal in both j a labeling the mode frequencies and n ja , the occupation number of each mode. The index a (also c) runs through the set of basis vectors of the Fock space that is built by the operators b + i . Hence, the dimension of the density matrix is formally infinite. The result of Hawking indicates that the density matrix for the outgoing radiation is thermal and thus maximally mixed. This is in contradiction to unitary evolution and at the core of the information paradox since the collapsing matter was in a pure state and evolved into a mixed state. Nonetheless, we will eventually show how the presence of off-diagonal elements in the density matrix ρ out can lead to the purification of the outgoing radiation in the course of the BH evaporation. 3 The multi-particle density matrix in terms of the singleparticle density matrix The goal of this section is to expose the connection of the full density matrix of the outgoing radiation ρ out , as defined in (2.7), to the widely used single-particle density matrix Let us begin with the well-known Bogolyubov transformation between the creation and annihilation operators at future infinity {b + i , b i } and at past infinity In terms of the Bogolyubov coefficients, we can express ρ i j as follows: As shown in Appendix A, there is a closed expression for the density matrix ρ out in terms of the b's, where Z is the normalization factor and Ω is a c-number Hermitian matrix. The latter is related to the single-particle density matrix ρ i j as , or e Ω = 1 + 1 ρ . (3.5) From now on, the superscript "out" on ρ will be omitted for brevity, The full density matrix ρ is thus completely defined by ρ i j =β ik β jk , with the normalization factor Z expressed as (3.7) One can define the entropy and related quantities directly in terms of the matrix ρ. In particular, an explicit expression for the Rényi entropy, readily obtained by evaluating From this expression, it would appear that the Rényi entropy does not vanish unless ρ = 0, when no particles have been emitted. However, the entropy for a state with N radiated particles can reach as low as H 2 = ln (2N) ≪ S BH (0), which is below the scale of validity of the approximations that were used in deriving Eq. (3.9).
The expression for the von Neumann entropy is more cumbersome but included for completeness, The state of the emitted particles in the Hawking model Having a formal expression for the density matrix at hand, we can now evaluate the Rényi entropy and the particle number in the state of the outgoing radiation in the Hawking framework. However, as made clear below, these quantities are formally proportional to an infinite sum over all the possible modes in the Fock space of the outgoing radiation. This divergence comes about from implicitly assuming an eternal BH that radiates for an infinite period of time and, hence, emits an infinite amount of particles. The physical situation, however, must be different: The number of particles emitted by the BH during a finite time interval is finite and depends on the initial mass of the BH. The goal of this section is to propose a way of counting the modes emitted by the BH which reach an observer at future infinity during some finite time interval. By implementing this method, we find that the physical quantities indeed become finite. For simplicity, the analysis is initially specialized to the case of the Hawking model.
In the current case, the matrix Ω i j = ω i T δ i j is diagonal, and Eqs. (3.7) and (3.9) for the partition function and the Rényi entropy then read while the total particle number is given by The sums in the above expressions run over the full set of the frequency modes ω i of a quantum field in the BH background. Hence, the dimension of the set of the possible frequencies is infinite, leading to an infinite entropy and particle number.

A reminder: A scalar field in a box
To better understand the correct way of counting the modes so as to render physical quantities finite, let us first recall a simple case: A massless quantum field in a box.
A massless quantum scalar field φ(t, x) in a box of size L with reflecting boundary conditions (i.e., φ(t, 0) = φ(t, L) = 0) can be expanded in terms of the Fourier modes where k x , k y and k z can be any integer. For a large enough box (L ≫ 1/T ), the sum over the modes can be approximated by an integral over a continuous momentum p, kx,ky,kz with the corresponding frequency given by ω( p) = | p|. In this way, the quantities in Eqs. (4.1)-(4.3) are proportional to the volume of the box V = L 3 ; that is, they are extensive. Moreover, since the only other dimensional parameter is the temperature T , any extensive, dimensionless quantity is given by V T 3 times some numerical constant: where x = ω/T . In general, due to the diagonal form of the thermal density matrix, any quantity that is given by a trace over the Hilbert space H out is proportional to a product over the available frequencies (as in Eq. (A.7)). A logarithm of such a quantity is then proportional to a sum over the frequencies and is extensive by virtue of Eq. (4.5). Imposing some physical boundary conditions on the scalar field results in a restriction on the allowed frequency modes-in this case, the condition that the wave numbers have to be integer multiples of 2π/L. The formally infinite quantities (4.1) -(4.3) then become finite, as should be true for any physical system.

Emission of localized wave packets
We have just seen that, for a field in a box, the frequencies are quantized in units of the inverse size of the box and the total number of modes is proportional to the volume of the box. In order to study how a BH radiates into an infinite space, it is more convenient to use localized wave packets with finite normalization instead of the Fourier modes (as in [1]).
Let us begin here with a scalar field that is propagating in the background of a BH. It can be expanded in terms of the complete set of solutions of the wave equation having only positive frequencies. At past infinity, the expansion contains only ingoing modes and takes the form where {a + i , a i } are the creation and annihilation operators defined at past infinity and the sum runs over a discrete set of modes {f i } with finite normalization. At future infinity, the expansion of the scalar field in terms of positive-frequency solutions contains both ingoing and outgoing modes, where {b + i , b i } are the creation and annihilation operators of the outgoing modes {p i } at future infinity.
In the continuum normalization with frequencies ω (rather than the discrete wave packets i), the ingoing and outgoing Fourier modes for the spherically symmetric solutions are expressible as where l and m are the angular-momentum numbers, and v and u are the advanced and retarded tortoise coordinates (respectively). A complete set of localized wave packets can then be defined as [1] j ε e −2πi n ω/ε p ωlm (u, r, θ, φ) dω . (4.14) In place of the continuous label ω, the wave packets have been labeled by two integer indices j ≥ 0 and n. The wave packet with index j contains waves with frequencies that are localized in the range j ε ≤ ω ≤ (j + 1) ε. Index n labels the ray along which the packet is propagating: The wave packet f jn is peaked around the advanced time v = 2πnε −1 and the wave packet p jn , around the retarded time u = 2πnε −1 , both having a width of 2πε −1 .
Our particular interest is the wave packets which were emitted by the collapsing body and then detected at some fixed distance away from the BH during a given period of time t 0 ≤ t ≤ t 0 + ∆t. Such wave packets will also be localized in the corresponding range of retarded time u 0 ≤ u ≤ u 0 + ∆t. For a long enough time interval, one can always choose the width ε −1 so as to ensure the detection of many wave packets (ε −1 ≪ ∆t), as well as maintain a fine-enough frequency resolution (ε ≪ T ) for each wave packet to be treated as a monochromatic mode of fixed frequency. In such a case, the summation over the wave-packet position n in the interval ∆n = ∆t/(2πε −1 ) can be approximated as whereas the summation over the discrete set of frequencies j = ω/ε can be approximated by the integral As one can now see, in this approximation, the total number of modes which can be detected during the time interval ∆t, does not depend on the choice of the parameter ε and is proportional to ∆t times the integral over mode frequencies, The total number of particles detected during the time interval ∆t in this case is given by where Γ ωlm , the so-called grey-body factor, is determined by the properties of the modes near the horizon. The emission of the modes with high multipoles is highly suppressed, Γ ωlm ≪ Γ ω00 ≡ Γ(ω) for l, m > 0, and the sum (4.18) can be well approximated with the l = m = 0 term only (see, e.g., [21]). In what follows, we will consider only the l = m = 0 modes and omit the angular-momentum labels. Moreover, in the Schwarzschild case, there is only one available dimensional parameter, which can be chosen as the BH temperature T . Therefore, the grey-body factor Γ(ω) depends on the frequency only through the ratio x = ω/T , and the number of emitted particles simplifies as follows: where A ∼ T −2 is the area of the BH horizon. Thus, the rate of BH emission coincides with the thermal emission rate of a body with area A. In reality, due to the non-trivial frequency dependence of the grey-body factors, the spectrum for BH radiation is quite different from the thermal spectrum of a body in an empty box. Some details about the grey-body factors are provided in Appendix B.
And so, in the absence of other dimensional parameters, all the extensive quantities describing the emitted radiation are proportional to the product ∆t T or, equivalently, to the total number of emitted particles (4.19). For example, the Rényi entropy H 2 is given by Hence, by properly accounting for the finite duration of detection, one finds that the number of particles emitted by the BH, as well as the entropy of the radiation, becomes finite as opposed to the idealized case of infinite emission time.

Single-particle density matrix in the wave-packet basis
In order to highlight the connection between the localized wave-packet description and the semiclassical model as presented in [13,14], we shall construct the singleparticle density matrix ρ i j in the finite wave-packet basis. Using the definition of the wave packets (4.14), one can express the matrix elements of ρ in terms of the same Fourier-mode basis, Each wave packet is labeled by the characteristic frequency ω j = j ε and the detection time t n = 2πε −1 n. The choice of ε parametrizes a trade-off between the frequency and the position resolutions. The total number of modes W , within a given range of frequencies ∆ω and an interval of time ∆t, does not depend on ε and is just given by W = ∆j · ∆n = ∆t ∆ω/(2π). In the case of a thermal distribution, the frequency range ∆ω is effectively set by the temperature. Hence, the total number of modes with substantial occupation numbers is of the order W = T ∆t. The matrix ρ jn j ′ n ′ in Eq. (4.21) is, therefore, effectively a W × W matrix. An important feature of BH radiation is that the emission rate Γ is the same as the temperature T . For a general radiating body, this is not true; even if Γ ∝ T , the proportionality constant need not be close to unity. This has important consequences, which we now elaborate on.
From the above discussion on wave packets, we have learned that the total number of emitted modes is given by W = ∆ω∆t, where ∆ω is the range of emitted frequencies and ∆t is the detection time. For the case of a nearly thermal emitter, the frequency range is given by the temperature ∆ω ∼ T , so that W = T ∆t. On the other hand, the total number of emitted particles during this same time is given by N = Γ∆t, where Γ is the emission rate. For a BH, with the emission rate Γ ∼ T , the number of emitted particles N is then of the same order as the total number of occupied frequency modes W ; i.e., N ∼ W . This fact has two important consequences for BH radiation: First, the average occupation number for the modes of radiation is approximately of order unity, as N/W ∼ Γ/T ∼ 1; meaning that, on average, each mode of radiation is occupied by only a few particles. Second, since W ∼ N, the single-particle density matrix ρ is effectively of size N × N. Then, as the average occupation numbers correspond to the diagonal elements of the matrix ρ, these elements go as ρ jn jn ∼ 1. Hence, for BH radiation, the two-point function ρ in the wave packet basis has the form of an N × N identity matrix.
The fact that, on average, each mode is occupied by only a few particles also implies that the radiation field in a BH background cannot be treated as semiclassical. This is not the normal state of affairs. For any other thermal emitter, the radiation rate is proportional to its area, Γ ∼ A · T 3 , and the area for classical emitters is much larger than the inverse temperature squared (i.e., the square of the typical wavelength), so that Γ/T ∼ AT 2 ≫ 1. Consequently, N/W ∼ Γ/T ≫ 1 and the radiation field is highly classical.
In order to make these ideas more precise, let us calculate tr ρ and tr ρ 2 for the Hawking model [1,2]. For a large BH, With this expression, in the wave-packet basis, As one can observe, ρ jn j ′ n ′ depends only on the difference in detection times t n − t n ′ . Because of the finite time and frequency resolution, it has non-zero offdiagonal elements with n = n ′ , which are concentrated in a narrow strip |n−n ′ | ε/T corresponding to |t n −t n ′ | T −1 . The number of non-zero elements of ρ is, therefore, proportional to its dimension W and, consequently, of the order of the number of emitted particles N.
The traces of ρ and ρ 2 can be explicitly calculated: In order to obtain the last line, we have used that (sin(∆nπx)/ sin(πx)) 2 ≈ ∆n δ(x) for large ∆n.
The number of the non-zero eigenvalues, which is just the number of diagonal elements in the Hawking model, can be then estimated by the participation ratio (i.e., the inverse of the purity), The result grows linearly with the number of emitted particles and also provides an estimate of the Rényi entropy (4.20).
For a more general case including off-diagonal elements, one can still expect tr ρ ∼ N, as this result only depends on the dimensionality of the matrix. On the other hand, tr ρ 2 can be much different, but its value can be readily estimated if the number of non-zero eigenvalues is known or vice versa.

Off-diagonal elements of the density matrix
With the previous framework in hand, we will now consider the structure of the off-diagonal elements in both the single-and multi-particle density matrices. This will be carried out for the three distinct models of an evaporating BH: (1) Hawking's model [1], (2) the Page model [10] and (3) the semiclassical model [13,14]. 3 In this section, all dimensional quantities are expressed in units of the BH temperature, N = N(t) is the total number of emitted particles at time t and S BH ≡ S BH (N) = S BH (0) − N(t) is the BH entropy at "time" N. This relation is modified, semiclassically, by fluctuations in the emission time of the Hawking modes. However, the fluctuations were shown to be small [13], being proportional to inverse powers of S BH .
One of the main issues to be clarified is to what extent the off-diagonal elements in terms of mode-occupation number of the full density matrix are suppressed. If they are indeed highly suppressed, then the relationship between the density matrix ρ and the single-particle density matrix ρ is reliable.

The Hawking model
In Hawking's model, the outgoing radiation is described to a very good approximation by the diagonal density matrix given in Eq. (2.10). What we would like to address here is how accurate an approximation this is by quantifying the off-diagonal contribution. 4 The Hawking density matrix ρ H is only approximately diagonal in two obvious ways: Elements between states differing in total occupation number are only approximately vanishing, as are those between states with some modes differing in frequency or in the time of emission. We would like to estimate the magnitude of the off-diagonal elements of the density matrix. As far as we know, this issue has not been addressed directly in the literature, although Hawking did expect such corrections to be exponentially small [1].
What we find below is an overall suppression for off-diagonal elements in modeoccupation number that is stronger than e −S BH . This is significant because it implies that the results of Hawking are not modified in a meaningful way by off-diagonal corrections. 5 We also find that the off-diagonal elements of the single-particle density matrix ρ i j are smaller than 1/S BH , leading to corrections to physical quantities, such as the entropy, that are suppressed by at least a factor of 1/S 2 BH . These are negligible at all times, as originally anticipated by Hawking. Corrections to the diagonal elements of ρ i j are similarly suppressed by at least 1/S BH and readily absorbed into the normalization.
The fact that the Hawking density matrix is diagonal can be understood by looking at Hawking's derivation of the single-particle density matrix ρ H , as defined in Eq. (4.22). The relevant factor goes as ρ H (ω, ω ′ ) ∼ I − (ω, ω ′ ) such that [1] where y = lnω has been used. This equation makes it clear that the matrix elements of ρ H are diagonal in frequency space. The off-diagonal elements of ρ H for different occupation numbers are also vanishing, as follows from an estimation of the products The products αβ andᾱβ arise in the expectation values of operators like and, thus, define the elements of the density matrix between states that differ in total occupation number by ∆N = 2k. The products αβ andᾱβ also appear in the computation of the generating function (A.2). Consequently, the dependence of the density matrix on ρ is only truly valid in the case when the products αβ andᾱβ are negligible. In Hawking's Fourier-space analysis, such terms contain an integral that is similar to the one in Eq. (5.1) but now the two frequencies appear with the same sign, Since both frequencies are positive, the argument of the delta function can only vanish if both of the frequencies vanish. This leads to one or more factors of δ(ω + ω ′ ) when contracting pairs of operators in the corresponding density-matrix element, as there would necessarily be unequal numbers of creation and annihilation operators within the expectation value.
To help reveal where corrections to off-diagonal elements can appear, let us consider the integration limits in Eqs. (5.1) and (5.5). These limits are an idealization, as a frequency can never be truly infinite nor exactly zero in physically realistic situations. For instance, the upper limit onω should be an exponentially large number, representing ultra-high frequencies on past null infinity that need to be red-shifted in order to produce the thermal radiation at future null infinity. An estimate of this upper limit in units of the BH mass M is ∼ e τ /M ∼ e (M/M P l ) 2 ∼ e S BH (0) , where τ ∼ M 3 denotes the lifetime of the BH (see, e.g., [22]). As for the lower limit oñ ω, a natural choice is to fix it as the inverse of the BH lifetime (τ /M) −1 ∼ S −1 BH (0). However, in Hawking's model, the BH is regarded as eternal, and so this limit should rather be the exponential of a large negative number.
In the following, we replace the idealized limits with suitable ultraviolet and infrared cutoffs, e ymax and e −y min , which are assumed to be of order y max ∼ y min ∼ S BH (0). We will further set y * ≡ y max = y min , as this symmetry will simplify the calculations without affecting the conclusions. This estimate determines the strength of the off-diagonal corrections of the density matrix in mode-occupation number. If, for some reason, the limits are such that y max ∼ y min ≪ S BH (0), then the original Hawking calculation is significantly modified and the original conclusion about the nature of the density matrix needs to be revised.
The tails of this "regulated" delta function provide the off-diagonal correction of interest. Notice that this comes about for both I ± (ω, ω ′ ), and, hence, the off-diagonal elements of the density matrix can now be non-vanishing, even when connecting states that differ in total occupation number.
To better understand the implications of these off-diagonal corrections, the singleparticle density matrix (and its related products) can be convolved with the same wave packets as in Eq. (4.21). This analysis is carried out in detail in Appendix C, from which there are two main conclusions: First, we find that the off-diagonal elements of the single-particle density matrix ρ jn j ′ n ′ ≡ (ββ) jn j ′ n ′ are suppressed relative to the diagonal elements by a factor of µ −2 * J −1 , where µ * ≡ y * ε/(4π 2 T ) and J = j ′ −j. Second, it is found that all the elements of the products αβ,ᾱβ are subleading and suppressed as µ −2 * J −1 αβ . Here, J αβ = j ′ +j +1 and can be treated as some typical average value of the frequency label j referring to the frequency range ω j ∈ [jε, (j + 1)ε]. These products set the order of magnitude of the off-diagonal elements of the full density matrix. In particular, this result tells us that the leading-order contribution to a typical off-diagonal element of the density matrix 6 is relatively suppressed by a factor of µ Here ∆N is the difference in total occupation number of the corresponding states and j * can be viewed as a typical value for the integer j.
Let us elaborate on this last claim by, first, understanding the structure of the density matrixρ. A natural way of organizing the elements of the density matrix which connect different Fock states is to classify them by the total number of particles in each state. A density-matrix element ρ N N ′ then denotes a block of elements connecting states with total occupation numbers N and N ′ . Thus, each element where d N is given by the size of the Fock subspace spanned by the states with total occupation number equal to N . If we denote the number of available frequency modes by W , then d N coincides with the number of distinct elements in a symmetric tensor of rank N in a W -dimensional vector space: Suppose that our interest is the strength of the density-matrix element ρ N N ′ such that ∆N = N − N ′ is an even, non-zero integer. (There are no contributions when this difference is odd; cf, Eqs. (5.3) and (5.4).) The strength of such an element is determined by a higher-order moment than the two-point function; rather, by an (N + N ′ )-point function with N annihilation operators and N ′ creation operators. But, because the theory of interest is non-interacting, the relative strength of such an element can be assessed by looking at products of two-point functionsββ, αβ,ᾱβ. Here, we will be able to correctly match a creation operator with an annihilation operator for all but |∆N | 2 of the total number of operator pairs N +N ′ 2 . Each mismatched pair-two creation operators or two annihilation operators-will then contribute one suppression factor of µ −2 * j −1 * . Formally, the density matrix can be represented in terms of the products ββ, αβ, andᾱβ aŝ Hence, a typical density matrix element between states with a difference in the total occupation number ∆N is suppressed by the product of ∆N 2 factors of µ −2 * j −1 * . As already stated, the off-diagonal contribution to ρ H -the single-particle density matrix-is suppressed by a factor µ −2 Since the final result, when traced over all the frequency modes, should not depend on the parameter ε, we can chose it to be of the order of the BH temperature, ε ∼ T . This amounts to saying that all the modes have the same frequency and, hence, j ∼ O(1). It is then clear that the off-diagonal contribution to ρ H is suppressed by at least y −2 * ∼ S −2 BH (0) < N −2 with respect to the diagonal elements and is, therefore, inconsequential to the moments of ρ H .
Nevertheless, the question of whether the off-diagonal elements can make a substantial correction to the moments of the full Hawking density matrix ρ H is well founded, as the number of off-diagonal elements is certainly very large. To address this concern, let us consider the regime in which N ∼ S BH because this is when significant changes to the nature of the BH radiation could be expected. To estimate the importance of the corrections for the Rényi entropy tr ρ 2 H (tr ρ H ) 2 we shall replace the hierarchical structure of the density matrix (5.8) by a uniform estimate of the typical value of the off-diagonal corrections. We do so by setting |∆N | ∼ N ∼ S BH for all the off-diagonal elements, so that the relative suppression factor goes as (µ 2 * j * ) − N 2 . The overall off-diagonal contribution to the Rényi entropy can now be calculated as the product of the dimensionality of the matrix d and the square of the relative suppression factor, d × (µ 2 * j * ) −N . The size d of the full density matrix for N emitted particles is given by where d n is the size of the n-particle subspace defined in Eq. (5.7), and W is the number of available frequencies. Since W ∼ N for BH radiation (see Subsection 4.3), this gives d ∼ e 2N . Thus, the overall off-diagonal contribution, relative to that of the diagonal, is e 2N (y 2 * ) −N ∼ e 2N N −2N ≪ e −N ∼ e −S BH , where we have again used that ε ∼ T and j * ∼ O(1), leading to the estimate µ 2 * j * ∼ y 2 * ∼ S 2 BH ∼ N 2 . It can be concluded that the off-diagonal contribution to the Hawking density matrix is exponentially suppressed. The physical reason behind this result is that the thermal radiation at future null infinity originates from ultra-high frequency modes at past null infinity which are highly red-shifted due to the presence of the horizon. We can therefore deduce that the existence of a region of high redshift-the horizon-induces a characteristic exponential suppression of the off-diagonal elements in mode-occupation number of the full density matrix.
Finally, let us use Eqs. (3.9) and (4.25), to determine the Rényi entropy for the Hawking model. As off-diagonal corrections are suppressed at least by order e −S BH at all times and the diagonal is approximately uniform, the single-particle density matrix always has N non-zero eigenvalues, each of which is unity up to negligible corrections. Hence, the estimate (H 2 ) Hawking ≃ N ln 3 (5.10) is valid at all times.

The Page model
In the Page model of BH evaporation [10], the BH and the emitted radiation are assumed to be in a pure state in some large Hilbert space H. This space is partitioned into an "in" part, H BH , representing the BH (including the ingoing radiation) and an "out" part, H out , representing the outgoing Hawking radiation, so that H = H BH ⊗ H out . The two Hilbert spaces are characterized only by their dimensionality. Let us label the states in H BH by i = 1, . . . , m, where m = e S BH , and the states in H out by A = 1, . . . , n, where n = e S rad . A state in H can then be written as |A, i . The working assumption is that both Hilbert spaces are large, n, m ≫ 1. Therefore, the Bekenstein-Hawking entropy S BH of the BH and the radiation entropy S rad ∼ N (which is calculated as if the radiation was in a thermal state) characterize the dimensionalities of the corresponding Hilbert spaces rather than the real entropies of the subsystems. Since the total state is pure, the actual entropy of the radiation and that of the BH are entirely due to entanglement and, thus, equal to one another at all times.
The partition between "in" and "out" is meant to mimic the horizon separating the interior of the BH from its outside. The variation in time of the dimensions of the Hilbert spaces H BH and H out is meant to model the evaporation of the BH, so that H BH shrinks and H out grows while the number of emitted particles increases. However, additional physical effects resulting from the existence of a horizon are not taken into account.
The density matrix of the Page model ρ P is given by The reduced density matrix for the radiation is then obtained by tracing over the BH Hilbert space, just as in Eq. (2.7), ρ out = tr BH ρ P . (5.12) Because the combined state is pure, there must be a special basis in which ρ P has only a single eigenvalue. This can be chosen, without loss of generality, to be given by ρ 1,1;1,1 . It follows that the matrix elements ρ A,i;B,j can be expressed in terms of a basis transformation matrix U acting on the original state vector |1, 1 . Denoting the resulting vector as V A,i = U A,i;1,1 , we have and 14) The vectors V A,i are meant to be treated statistically and assumed to be random vectors of unit size with a uniform distribution on an mn-dimensional sphere. The initial investigation of such a random system was conducted by Lubkin [11] and then improved by many subsequent investigations [23][24][25]. The final result is that the distribution of eigenvalues of ρ out was found to obey the Marchenko-Pastur (MP) law [26], which is the eigenvalue distribution for a certain ensemble of semipositive-definite, random matrices.
The MP distribution is generically divided into two parts, one set of large eigenvalues and another set of vanishingly small eigenvalues. It depends on only two parameters, the total number of the eigenvalues and the participation ratio corresponding to the number of large eigenvalues. The distribution of the eigenvalues (up to a normalization convention) is given by where The early times of BH evaporation, when n ≪ m, correspond to c ≫ 1 and P R ∼ n. Essentially, all the eigenvalues are about equal and their value is determined by the normalization convention. This means that the Rényi entropy of the radiation is the same as that of the Hawking model, (H 2 ) Page ∼ N ln 3. But well after the Page time, when n ≫ m, one finds that c ≪ 1 and P R ∼ m ≪ n. In this case, about n − m ∼ n of the eigenvalues become zero and both reduced density matrices for the radiation and BH are becoming pure. The Rényi entropy of the radiation for the Page model coincides with that of BH and is given by This is depicted in Fig. 1 as a function of the number of emitted particles. Since the entropy of the subsystems is entirely due to entanglement, it is symmetric if one exchanges the Hilbert spaces of the BH and radiation. As a consequence the evolution of the entropy is symmetric under reflections around the Page time, when the sizes of H BH and H out are equal; that is, when half of the particles have been emitted, N P age ≈ S BH (0)/2.

The semiclassical model
The semiclassical model [13,14] improves upon Hawking's framework by taking into account the BH's quantum fluctuations, as well as its time dependence due to the emission of the radiation. In this setup, the single-particle density matrix ρ SC is no longer diagonal and, as a result, the evaporation process becomes unitary even though the thermal-like emission spectrum is kept. The basic prescription is to assign a quantum wavefunction to the collapsing shell of matter in Hawking's model and then recalculate all relevant quantities as expectation values. The main outcome is that ρ SC picks up off-diagonal contributions that are uniform in terms of frequency but suppressed relative to the diagonal elements by C 1/2 BH (N), where C BH (N) = S −1 BH (N) ≪ 1 is a classicality parameter-a "time-dependent " that keeps track of how close the system is to the limit of a classical spacetime.
The elements of ρ SC do have a non-uniform suppression in terms of emission time; modes emitted at different times tend to decohere. Nonetheless, if the radiation is being regularly monitored at intervals of ∆N ∼ √ S BH or less, then this suppression can be compensated. We will specifically be considering this case, which has been called the "tracking case" in [17]. Hence, the off-diagonal elements of the singleparticle density matrix ρ SC can be regarded as uniform in magnitude with respect to both frequency and emission time.
As explained in Subsection 4.3, ρ can be viewed as an N × N matrix, with the indices running over the wave-packet modes with non-vanishing occupation number and with the diagonal elements given by the average occupation number for each mode. The elements of the semiclassical single-particle density matrix are found to be [13,14] (ρ SC ) ii = 1 , where the phases θ ij can be treated as random for most purposes (but see below). Let us next consider the full density matrix ρ SC for this model. In [13] (see, in particular, Sect. 2.4), it was argued that the off-diagonal elements in modeoccupation number of the density matrix will be suppressed much in the same way as for the Hawking model, as discussed in Subsection 5.1. The physical reason for this is the same as before-the presence of the horizon implies that the frequencies observed at future null infinity are highly red-shifted and, thus, determine the magnitude of the suppression factor of the off-diagonal elements. The technical explanation is as follows: In the semiclassical model, the terms that lead to off-diagonal elements in mode-occupation number are, again, terms of the form αβ andᾱβ as in Eq. (5.2). These contain a classical term, which is equal to the Hawking-model contribution, along with a semiclassical correction. Irrespective of the details, it is clear that the semiclassical corrections to αβ andᾱβ vanish as some power of C BH in the limit C BH → 0. This is enough to guaranty the exponential suppression of the semiclassical contribution to the off-diagonal elements in mode-occupation number.
Hence, there will be a hierarchical structure in the suppression factors of the off-diagonal elements in mode-occupation number, as shown in (5.8). Like before, the off-diagonal elements can be expressed in terms of (N + N ′ )-point functions, which will factorize into a product of N +N ′ 2 two-point functions. The strength of an (N +N ′ )-point function will then be suppressed by a factor similar to the suppression factor of the Hawking model.
In some sense, this semiclassical model can be viewed as the "middle ground" between the other two models; on one hand, remaining almost thermal like Hawking's with the associated hierarchical structure of the density matrix but, on the other, evolving over time like Page's. At the early stages of BH evaporation, the semiclassical model is essentially Hawking's plus small corrections. However, at later times (after the Page time), the dominant contributions will come from those elements that are off-diagonal in frequency. This can be attributed to the effective perturbative parameter being NC BH for this framework [13], as this parameter grows monotonically throughout the evaporation process and finally becomes large (NC BH > 1) at times later than the Page time.

Semiclassical Rényi entropy
The objective here is to estimate the Rényi entropy of the BH radiation H 2 (N) using the single-particle density matrix ρ SC as specified in Eq. (5.21). However, we are currently lacking a precise knowledge of the eigenvalue distribution for ρ SC , which is necessary for the use of Eq. (3.9), This deficiency can be attributed to our lack of specification as to the nature of the phases e iθ ij . These phases are not entirely random but rather dictated by the effect of the quantum horizon on correlations in the emitted radiation [14]. We do expect that the eigenvalue distribution is similar in some aspects to the previously discussed MP distribution. The moments of the distribution function of the eigenvalues of ρ SC are calculated and compared to the moments of the MP distribution in Appendix D. Although differing in details, the two sets of moments follow a similar pattern. In this subsection, we will instead call upon the known properties of ρ SC in order to draw conclusions about the entropy of the radiation.
We expect that ρ SC has a number of large eigenvalues and a number of vanishing eigenvalues, with the large eigenvalues distributed about the mean in a way that depends on their total number. We can gain some understanding about the distribution of the eigenvalues by using the participation ratio P R along with some symmetry arguments.
The participation ratio provides, for a given matrix, a measure of the number of non-vanishing eigenvalues. For ρ SC , it is given by [17] Then, since tr ρ SC = N (as shown in Subsection 4.3), the average value of a nonvanishing eigenvalue can be estimated as This finding reveals that, up to the Page time (NC BH < 1), there are about N eigenvalues of order one. However, after the Page time (NC BH > 1), the average λ is larger than one; signaling that the matrix ρ SC has a certain number of large eigenvalues while the majority vanish.
Let us now establish what must be true about the Rényi entropy for this semiclassical model. What is known with certainty is that H 2 (ρ SC ) (1) must tend to the Hawking value of N ln 3, as in Eq. (5.10), at early times and (2) must be symmetric under an exchange between the radiation and the BH, since the two reduced density matrices necessarily share the same set of eigenvalues. In particular, all physical quantities are required to be symmetric under the exchange of N ↔ S BH (N), as long as these quantities characterize the sizes of their respective Hilbert spaces H out and H BH .
The participation ratio (5.23) is manifestly symmetric in just this way, and of the same form as the participation ratio of the MP distribution (5.17), given that one identifies the parameters n and m of the Page model to be N and S BH (N), respectively. 7 This may seem peculiar because, as discussed in Appendix D, the actual eigenvalue distribution of ρ SC differs from the MP law. But, as far as we can tell, N S BH (N ) S BH (N )+N is the simplest non-trivial, symmetric function which tends to N for N ≪ S BH (0), suggesting that such an identification is valid.
Motivated by these observations, our expectation is that the semiclassical Rényi entropy can be expressed in terms of the participation ratio as follows: where the coefficients are dimensionless numbers that are determined by the higher moments of the eigenvalue distribution. To see that c 0 = ln 3 is indeed the correct leading-order coefficient, consider that, at early times, N ≃ P R to leading order in a 1/P R expansion. Thus, the average eigenvalue λ = N P R ≃ 1 and the logarithm in H 2 then expands as ln (1 + 2λ) ≃ ln 3.
The dependence (5.26) of the Rényi entropy on the number of emitted particles for a semiclassical BH is presented on Fig. 1. Because the Rényi entropy depends only on the participation ratio N 1+N C BH , one can determine the evolution of the former from that of the latter. The participation ratio steadily grows until reaching a maximum at the Page time and then steadily declines for the remainder of the evaporation process. This is qualitatively similar behavior to the Rényi entropy for the Page model although, in the semiclassical model, the rate of purification-which is the rate of deviation from the linearly growing result of the Hawking's model -is actually faster than the Page-model rate [17].

Conclusion
We have exposed the relationship between the density matrix of the outgoing radiation from a BH and the corresponding single-particle density matrix in the case that the density matrix is approximately diagonal in mode-occupation number. It was then shown that the presence of the horizon leads to a high suppression of the offdiagonal elements in mode-occupation number of the density matrix for the emitted radiation. We have therefore concluded that the density matrix is approximately diagonal in mode-occupation number at all times. We have also shown how to regularize the infinities which arise as a consequence of Hawking's idealized picture of an eternal BH. It was explained how the state of the emitted radiation from a BH is different from that of a standard black body of the same temperature. This analysis was then applied to three models of BH evaporation as a means of clarifying their differences and contrasting their main features. Let us briefly summarize the main observations: The Hawking model: The off-diagonal elements of the single-particle density matrix are highly suppressed for different frequencies and emission times, and those of the full density matrix are exponentially suppressed for differences in occupation numbers. The Rényi entropy (4.20) scales with the number of emitted particles N at all times, and so there is no possibility for the radiation to be purified.

The Page Model:
For a randomly chosen basis whose relationship to Hawking's basis is unspecified, the off-diagonal elements of the density matrix are uniform in magnitude with random phases. There is no hierarchy between different elements of the matrix. The Rényi entropy of the radiation (5.20) increases linearly as N until the Page time and decreases linearly as S BH (0) − N afterwards. Hence, purification is inevitable.

The semiclassical model:
The off-diagonal elements of the single-particle density matrix are suppressed by a power of S −1/2 BH in a uniform way, whereas those of the density matrix are exponentially suppressed for different mode-occupation numbers in a similar way to the Hawking model. At the Page time, the contribution from the off-diagonal elements with regard to frequency and emission time becomes significant because of a perturbative parameter that grows monotonically with time. The Rényi entropy (5.26) scales with the participation ratio of the single-particle density matrix, which increases until the Page time and decreases thereafter. This suggests that the radiation purifies and does so at a rate which is faster than that of the Page model.
A Expression for the full density matrix in terms of ρ i j In order to find an expression for the density matrix, it is useful to calculate a vacuum expectation value; namely, The resulting function of µ i and λ j can be used as a generating function for the expectation values of any normal-ordered powers of the creation and annihilation operators b + i and b j . The expectation value (A.1) can be calculated by using (3.2) and applying the Baker-Campbell-Hausdorf (BCH) formula for the exponents: The exponential factors on the outside of the average arise from the commutators. By using Hawking's expressions in [2] for the Bogolyubov coefficients α ij and β ij , one can show that and so the external exponentials are equal to unity as their exponents are vanishing. Relations (A.3) also ensure that only the matrix elements of ρ out between states with the same total occupation number are non-vanishing. But, as discussed in Section 5, these relations hold only approximately for realistic BHs.
Returning to Eq. (A.2), one can see that the two outer factors within the expectation value disappear after acting on the vacuum |0 − . The remaining expression can be further simplified by applying the BCH relation once again: Finally, the answer for the generating function (A.1) reads where ρ i j is the single-particle density matrix as defined in Eq. (3.3). For Hawking, ρ i j is a diagonal matrix but, in more general setups, it can also have non-zero offdiagonal elements. In any case, ρ i j is Hermitian and can thus be diagonalized by a unitary transformation.
We next look for a closed expression for the density matrix ρ out in terms of the b's. In analogy with the harmonic oscillator, let us try the ansatz where Ω is some c-number Hermitian matrix and Z is the normalization factor. The matrix Ω can be diagonalized by a unitary transformation, Ω = U + DU, with D being a real diagonal matrix. This allows us to introduce new annihilation operators are also a good set of creation and annihilation operators.
Moreover, since this transformation does not mix the creation and annihilation operators between each other, the subspaces with a fixed total occupancy of b-and d-particles coincide, and the traces over the Fock spaces of b-and d-particles are equal. In particular, one can find the normalization Z by summing over the states with a fixed occupation number of d-particles |n i : The last equality is due to the invariance of the determinant under unitary transformations of the matrix.
According to ansatz (3.4), the corresponding generating function (A.1) is If it is possible to choose Ω so that the result matches (A.5), then ansatz (A.6) gives the correct expression for the density matrix.
As before, the generating function (A.8) can be calculated by tracing in the d-particle basis: where each of the factors in the bottom line corresponds to a single oscillator. In order to evaluate the above expression, we will use the identity 1 − e −D tr e a + µ e λ a e −D a + a = exp λµ e D − 1 , (A. 10) which is derived in [27]. (For a web resource, see, e.g., [28].) The expectation value (A.8) then reads Comparing the last expression to (A.5), we find that the right-hand sides coincide if the matrix Ω is related to ρ as , or e Ω = 1 + 1 ρ . (A.12)

B Grey-body factors
From the properties of the grey-body factors, one can see that the state of the field outside a large BH is quite different from the usual thermal state in an empty box, ρ = 1 Z e −β H . The actual state of the radiation is that of a field in thermal equilibrium with the BH, which is very selective in what it absorbs and emits. For example, the occupation number of the modes with high angular momentum will be much less than what one would naively expect from just the Boltzmann suppression.
One way to present this modification away from thermality is to notice that, by definition, the grey-body factor Γ ωlm is equal to the modulus square of the scattering amplitude of the corresponding mode on the BH. Therefore, it is physically intuitive to introduce a corresponding effective scattering cross-section, The rate of emission can then be written as This is a three-dimensional thermal emission rate, where the geometric area of the emitter is different for each mode and is given by the absorption cross-section. In this language, the suppression of the high-l modes comes from the fact that the effective area of the BH as seen by them is much smaller than its geometric area. Another way to see the difference from the usual thermal state is to notice that the density matrix of the equilibrium state has the form ρ = Z −1 e −b + Ωb with Ω = ω/T only if Γ ωlm = 1. Indeed, the relation (3.5) between the single-particle density matrix ρ = Γ ωlm e ω/T −1 and the matrix Ω implies that Therefore, the density matrix for the radiation takes the form [cf. Eq. (2.10)] For each given mode, the ratio of the probability of having (N + 1) particles to the probability of having N particles is constant and given by Meaning that, apart from the overall Boltzmann suppression factor, the emitted particle has a probability of (1 − Γ) to be scattered back into the BH,

C Off-diagonal corrections to the Hawking model
In order to estimate the off-diagonal corrections to the Hawking density matrix, we first need to consider the corrections to the products where the second index is summed over. The first product defines the standard single-particle density matrix ρ i j ≡β ik β jk and, in the case of Hawking's calculation, is diagonal and completely determines the density matrix of the outgoing radiation (3.4). The products αβ andᾱβ arise in the expectation values of operators like and, thus, define the elements of the density matrix between states that differ in total occupation number, ∆N = 2k.
The products αβ andᾱβ also appear during the computation of the generating function (A.2). Consequently, the sole dependence of the density matrix on the single-particle density matrix ρ i j is only truly valid in the case when the products αβ andᾱβ vanish. This assumption is true for Hawking's idealized calculation but, as shown below, modified for physically realistic BHs.
In order to find the corrections to the products (C.1), we rewrite the definition of the coefficients α ωω ′ and β ωω ′ in the basis of Fourier modes, and v 0 denotes the position of the BH horizon in advanced time, κ = 2πT is the BH surface gravity and t ω is the transmission coefficient for which |t ω | 2 = Γ ω .
The products (C.1) can be rewritten in the wave-packet basis of Section 4. For instance, There are similar expressions for (αβ) jn j ′ n ′ and ᾱβ jn j ′ n ′ . The product (C.8) and its analogues can be expressed in terms of the Fourier-mode basis as where the integration variable has been changed to y = ln(ω).
Substituting the previous set of relations into the corresponding expressions for the products in the wave-packet basis and then integrating over the frequencies ω and ω ′ , we arrive at whereω ≡ ε (j + 1/2) andω ′ ≡ ε (j ′ + 1/2) are the mean values of the frequencies in the given range of integration,ñ ≡ n − n 0 ,ñ ′ ≡ n ′ − n 0 with n 0 defined by v 0 = 2πn 0 ε −1 , and the integration variable has been changed to µ = yε/(2πκ). We see that the three integrals differ only by the power in the exponent, and so it is convenient to define Let us now introduce infrared and ultraviolet cutoffs for the integration range of the frequencies, as discussed in the Subsection 5.1, (C. 16) Due to the similarity of the three integrals above, we need only to consider one, and then interpret the results for the different choices of J in Eq. (C.15). The integral does need, however, to be treated differently for the cases with ∆n ≡ n ′ − n = 0 and ∆n = 0, as well as for the cases with J = 0 and J = 0.

C.1 Off-diagonal elements in frequency
When ∆n = 0, the integral (C.17) becomes For the case J = 0, the previous integral can be evaluated to give where Si(x) = x 0 dt sin t/t is the sine integral. It has the large-x expansion Hence, in the limit when µ * ≫ñ, the integral becomes When J = 0, the integral can be split as (C. 23) We find that the last integral is expressible as In the limit when µ * ≫ñ, the function h(J) then becomes The integral I 0 (J) at a given value of J can then be evaluated by combining Eqs. (C.21), (C. 22) and (C.26). For example, I 0 (1) works out to be

C.2 Off-diagonal elements in mode-occupation number
When ∆n = 0, the integral (C.17) can be split as If J = 0, the integral can be computed exactly. For J = 0, it is possible to redefine the integration variables in the two summands so as to move the dependence on the expansion parameter µ * from the limits to the integrand: One can then expand the integrand in the limit when µ * ≫ñ,ñ ′ and evaluate the integral to obtain On the other hand, at J = 0, the integral I ∆n differs from I 0 (0) in Eq. (C.21) only by a zeroth-order, diagonal term: Hence, the corrections in both cases, ∆n = 0 and ∆n = 0, coincide and are uniform in ∆n.

Interpretation
The matrix elements of the product β β jn j ′ n ′ are proportional to β β jn j ′ n ′ ∼ I ∆n (Jβ β ) , Jβ β = j ′ − j . (C.36) Since the corrections are found to be uniform in ∆n = n ′ − n, we can represent the matrix β β jn j ′ n ′ as a tensor product of j's and n's. Then each matrix element that is labeled by j and j ′ is itself a uniform matrix in terms of n and n ′ . In this way, the product can be represented in the following form: where each depicted entry represents a uniform block and only the order of magnitude of the leading term is shown. For the products αβ andᾱβ, the corresponding index is J αβ = ±(j ′ + j + 1), respectively. Since both j and j ′ are positive, these products have no zeroth-order contributions and only the matrix elements with j = j ′ = 0 have corrections of the order O(1/µ * ). All the other entries of the matrices αβ andᾱβ receive only second-order corrections O (1/(µ 2 * J)). How these findings impact upon the off-diagonal elements of the full density matrix is discussed in Subsection 5.1 of the main text.

D Higher moments of the Marchenko-Pastur and ρ SC distributions
Here, we will quantify more precisely the differences between the Marchenko-Pastur (MP) distribution and the eigenvalue distribution of ρ SC . This entails assigning the distributions the same participation ratio and then determining how their higher moments are different. Comparing their respective participation ratios (with n = N) in Eqs. (5.17) and (5.23), one can readily identify the parameter c for the semiclassical two-point function, . (D.1) Let us consider the case NC BH 1. This is really the regime of interest, since smaller values of NC BH correspond to the case in which the Hawking model is valid (up to small corrections). We will further assume that the phases of the off-diagonal terms in ρ SC can be treated as random. One can then estimate the higher moments of ρ SC up to combinatorial factors and sub-leading terms in small 1/NC BH , which leads to 3) The basic idea behind these estimates is that the off-diagonal parts of the matrices are dominant when NC BH 1 and the randomness of the phases requires these parts to sum coherently (i.e., restricted to sums of the form ij M ij M ji ). A simple example should suffice to illustrate the point. Let γ and η be the diagonal and offdiagonal parts respectively of a matrix ρ. Applying the rule of coherent summation and, otherwise, insisting on the maximum power of η, we have, for the p = 2 case, tr ρ 4 = ijkl ρ ij ρ jk ρ kl ρ li ∼ This follows from the observation that the MP distribution (5.15) has, for small values of c, about Nc ∼ 1/C BH large and (roughly) equal-valued eigenvalues λ ∼ 1/c. This discrepancy between the MP distribution and the semiclassical distribution is a consequence of the square root of C BH appearing in the off-diagonal elements of ρ SC . Hence, ρ SC does not precisely conform to an MP distribution nor should it necessarily be expected to. However, when NC BH ∼ 1, both expressions for the moments of the distributions scale in the same way, tr(ρ n ) (tr ρ) n ∼ (C BH ) n−1 . Therefore, we do expect that the two distributions share the same general features; in particular, once c = 1 (NC BH = 1) is reached, both eigenvalue distributions begin to develop support for zero eigenvalues.
Yet, when NC BH ≫ 1 is true, the higher moments of the eigenvalue distribution of ρ SC are much more suppressed than those of the MP -distribution. For instance, the n th semiclassical moment is smaller by a relative factor 1/(NC BH ) n/2 than its MP counterpart. Nevertheless, we do expect that, in this case, both distributions have only a few large eigenvalues but apparently differ in the detail. It would be interesting to find out what is the actual eigenvalue distribution of the semiclassical matrix.