We extend previous efforts to quantify the entanglement generated in Hawkings evaporation process by including rotation and thermal environments e.g. the cosmic microwave background. Both extensions are needed to describe real black holes in our Universe. Leveraging techniques from Gaussian quantum information we find that the black holes ergoregion is an active source of quantum entanglement and that thermal environments drastically degrade entanglement generation. Our predictions are suitable to be tested in the lab using analog platforms and also provide tools to assess the fate of quantum information for black holes in more generic settings.
['General Relativity and Quantum Cosmology', 'High Energy Physics - Theory', 'Quantum Physics']
Entanglement from rotating black holes in thermal baths
{'Entanglement from rotating black holes in thermal baths': "Ivan Agullo, 1, 2, ∗ Anthony J. Brady, 3, † Adri'a Delhom, 1, 4, ‡ and Dimitrios Kranas 1, §  \n1  \nDepartment of Physics and Astronomy, Louisiana State University, Baton Rouge, LA 70803, U.S.A. 2 Perimeter Institute for Theoretical Physics, 31 Caroline Street North, Waterloo, Ontario, Canada N2L 2Y5  \n3 Department of Electrical and Computer Engineering,  \nUniversity of Arizona, Tucson, Arizona 85721, USA  \n4 Laboratory of Theoretical Physics, Institute of Physics,  \nUniversity of Tartu, W. Ostwaldi 1, 50411 Tartu, Estonia  \nWe extend previous efforts to quantify the entanglement generated in Hawking's evaporation process by including rotation and thermal environments (e.g. the cosmic microwave background). Both extensions are needed to describe real black holes in our universe. Leveraging techniques from Gaussian quantum information, we find that the black hole's ergoregion is an active source of quantum entanglement and that thermal environments drastically degrade entanglement generation. Our predictions are suitable to be tested in the lab using analogue platforms and also provide tools to assess the fate of quantum information for black holes in more generic settings.", 'I. INTRODUCTION': "In 1974, Stephen Hawking showed that black holes emit radiation as hot bodies [1, 2], providing an unforeseen relation between relativity, thermodynamics, and quantum mechanics. The quantum nature of this process becomes manifest by noticing that the emitted Hawking radiation is entangled with modes that fall into the interior of the black hole.  \nA milestone in the characterization of entanglement generated by non-rotating black holes was achieved by D. Page in Refs. [3, 4] (see also Ref. [5]). Page considered quantum fields initially prepared in the vacuum state. Since this is a pure quantum state, the von Neumann entropy of the radiation emitted to infinity quantifies the entanglement with the radiation falling inside the black hole. These calculations generated great interest; specifically because they provide the tools to discuss the information loss paradox [6-10] in quantitative terms [3, 4].  \nThe goal of this article is to extend Page's calculations in two important directions: (i) incorporating rotation and (ii) accounting for thermal environments. Both generalizations are essential for a comprehensive description of realistic black holes in our universe. For instance, rotation-prevalent for black holes of astrophysical origin-gives rise to an ergoregion encompassing the black hole horizon. It is well understood that the ergoregion contributes to Hawking radiation by amplifying its intensity via the phenomenon of superradiance [2, 11-13].  \nIn this article, we elucidate that entanglement generated by a rotating black hole contains two contributions of different physical origin-a thermal contribution from the horizon and a non-thermal part from the ergoregion.  \nBoth contributions are described mathematically in similar terms, by two-mode squeezers. A notable distinction is that the ergoregion is seeded by the thermal Hawking radiation from the horizon. The total entanglement results from the interplay of these two phenomena. To the best of our knowledge, there has been no previous investigation into characterizing the entanglement generated by the ergoregion.  \nThat the evolution of field modes across the ergoregion is formally described by a two-mode squeezer was first noticed in [14]. This implies that the vacuum state undergoes a transformation into a two-mode squeezed vacuum made of entangled pairs of particles-this is a vacuum instability, first pointed out by Starobinski [12] and confirmed by Unruh [13] (see also [15] for a discussion on a related instability of rotating compact stars that feature an ergoregion.).  \nThe scenario considered in this article differs from the Starobinski-Unruh framework in that, if the black hole originates from gravitational collapse, quantum fields in the vicinity of the ergoregion are not in the vacuum state, but instead contain a thermal flux of particles (Hawking radiation) emanating from the outer horizon. The ergoregion amplifies these thermal Hawking quanta by the process of superradiance, imprinting additional entanglement in the radiation. This contribution to entanglement differs from the one that the ergoregion would create if Hawking radiation were not present, as in the idealized Starobinski-Unruh scenario. Although the total radiation emitted by rotating black holes has been extensively studied, including the combined contribution of the horizon and the ergoregion [2, 11, 16], no analysis of entanglement-further decoupling each contributionhas been explored so far.  \nWe furthermore demonstrate that, for black holes immersed in thermal baths, the exterior thermal radiation alters the entanglement generated by the black hole, appreciably degrading entanglement if the temperature of the environment is larger than the Hawking temperature. This fact has been recently pointed out by us in the  \ncontext of analogue horizons created in the lab [17, 18] (and qualitatively discussed for non-rotating black holes in the essay [19]), and is of obvious physical interest for astrophysical black holes, considering that all of them are immersed in the cosmic background radiation. Furthermore, thermal environments are ubiquitous in analogue black holes created in laboratory settings.  \nThe quantitative analysis of these effects requires other tools than the ones used by Page. In particular, thermal environments are described by mixed quantum states and, consequently, the von Neumann entropy of the outgoing radiation is no longer a quantifier of entanglement. Weuse a novel way of deriving the Hawking process based on general techniques used to describe Gaussian bosonic quantum systems. These tools make it possible to derive the Hawking effect in a remarkably simple yet powerful manner, allowing us to isolate and quantify the contribution of the ergoregion.  \nOn the other hand, our tools only apply to states of the fields that are Gaussian (either pure or mixed). Although this is an intrinsic limitation, the formalism includes many of the states that are physically interesting: vacuum, coherent, squeezed and thermal states. Hence, the limitation is not significant in practice.  \nIn what follows, we neglect contributions from particle species with rest-mass energies greater than the thermal energy of Hawking radiation, which, for solar-mass black holes or heavier, practically restricts our calculations to the known massless particles, namely photons and gravitons [4]. We use units in which ℏ = c = k B = 1.", 'II. MODES AND THEIR EVOLUTION': "Figure 1 shows the Penrose diagram of the formation of a rotating, axisymmetric black hole (only the region exterior to the outer horizon is shown). Following Hawking's original derivation, the evolution of quantum fields can be obtained by solving the classical wave scattering. Let us first identify the wave modes, describing either electromagnetic or gravitational waves, involved in the Hawking effect (see Fig. 1).  \nConsider an 'out' mode at future null infinity ( I + ). This is described by a wave packet sharply centered around a frequency ω , angular momentum numbers ℓ and m , and helicity s , and supported mostly at late values of the retarded time u along I + . We will denote this mode by φ (out) Q , where Q = ω, ℓ, m, s , and we have suppressed vectorial and tensorial indices in the modes. The evolution of this mode backward in time to past null infinity ( I -) produces a combination of three modes, which, following [20], we denote as φ (in) Q , φ (p) Q and φ (d) Q . The mode φ (in) Q is a wave packet with support centered at late values of the advanced time v in I -, in the region where the spacetime is described by Kerr's metric. Due to the stationary and axisymmetric character of the underlying geometry in this region of the spacetime, φ (in) Q has  \nFIG. 1. Spacetime representation of the Hawking process (only the region exterior to the outer horizon is shown). For each choice of ω , ℓ , m , s , the Hawking effect results from the evolution of three modes at the Cauchy hypersurface Σ i to three modes at the Cauchy hypersurface Σ f . The intermediary mode 'up' is introduced to describe the evolution through the Cauchy hypersurface Σ int and to separate the total evolution in two steps.  \n<!-- image -->  \nthe same frequency and angular momentum numbers as φ (out) Q .  \nThe modes φ (p) Q and φ (d) Q have support only around the time v H at which the horizon forms. Although they have the same ℓ and m numbers as φ (out) Q , they do not have a well-defined frequency at I -. Due to the blueshift induced by the horizon, these modes have support on exponentially large frequencies. Hence, the label ω of these modes should not be confused with their physical frequency at I -. The explicit expression of these modes can be found in Appendix A. These modes in I -, which describe the progenitors of the 'out' mode, were introduced in [11], and are known as Wald basis modes (see [20] for a detailed account). The use of Wald basis introduces a drastic simplification-factorizing the evolution of infinitely many modes into decoupled Q -sectors and reducing the analysis to evolving three modes at a time.  \nIf we start with an arbitrary combination of φ (p) Q , φ (d) Q and φ (in) Q at Σ i = I -, the evolution produces, at the final Cauchy hypersurface Σ f = I + ∪ H ( H denotes the event horizon), a combination of three modes, φ (out) Q and two extra modes which we denote by φ (dn) Q and φ (down) Q (see Fig. 1). The evolution is made of two contributions of distinct physical origin, which we discuss separately. The  \nfirst contribution corresponds to the evolution from Σ i to Σ int (the form of the Cauchy hypersurface Σ int is depicted in Fig. 1; 'int' stands for 'intermediate'). This part of the evolution occurs in the region of spacetime that is time-dependent, in which the horizon is formed. The evolution from Σ int to Σ f takes place in the stationary Kerr geometry.  \nFor the free field theories we consider in this paper, the evolution is linear and can be written as matrix multiplication as follows. Let  \n⃗ A (i) Q ≡ ( a (in) Q , a (p) Q , a (d) Q , a (in) † Q , a (p) † Q , a (d) † Q ) ⊤ (1)  \nbe the column vector of annihilation and creation operator associated with the modes defined in Σ i , and let  \n⃗ A (int) Q ≡ ( a (in) Q , a (up) Q , a (dn) Q , a (in) † Q , a (up) † Q , a (dn) † Q ) ⊤ (2)  \nbe a similarly defined vector associated with the modes at Σ int . The modes 'dn' and 'up' are depicted in Fig. 1, and their explicit form is written in the Appendix A3. Let us consider a black hole with mass M and angular momentum J . The surface gravity is κ = ( r + -r -) / [2( r 2 + + a 2 )], where r ± = M ± √ M 2 -a 2 , a ≡ J/M and Ω H = a ( r 2 + + a 2 ) -1 is the 'angular velocity' of the horizon.  \nThe Heisenberg evolution from Σ i to Σ int is encoded in the relation ⃗ A (int) Q = S (H) Q · ⃗ A (i) Q , where S (H) Q = ( A Q B Q B Q A Q ) and  \nA Q = ( 1 0 0 cosh z H I 2 ) , (3)  \nB Q = ( 0 0 0 sinh z H ¯ I 2 ) , (4)  \nThe label H in S (H) stands for horizon. Here, z H ( ω, m ) = tanh -1 e -2 π | ˜ ω | / (2 κ ) , and ˜ ω ≡ ω -m Ω H . These expressions are valid for both photons and gravitons. We have denoted by I n and ¯ I n the diagonal and anti-diagonal unit matrices in n dimensions, respectively. We note that cosh z H →∞ as ˜ ω → 0; however, this mathematical divergence is tamed by the evolution across the ergoregion because the transmission probability across the ergoregion vanishes in the limit ˜ ω → 0. The elements of S (H) Q , as well as those of other evolution matrices written below, are called Bogoluibov coefficients. These coefficients are derived by solving the classical equations of motion, as detailed in the Appendix A.  \nInterestingly, the form of S (H) Q reveals that this operation precisely corresponds to a two-mode squeezing process between modes 'p' and 'd', to produce 'dn' and 'up' (and the identity operation for 'in'). A two-mode squeezer converts the vacuum state into a 'two-mode squeezed vacuum', which is an entangled Gaussian state such that the reduced state of each individual mode is a mixed thermal state. The occupation number of either mode is equal to ¯ n = sinh 2 z H = ( e 2 π | ˜ ω | /κ -1) -1 ,  \ncorresponding to a black body spectrum at temperature κ/ 2 π .  \nThe evolution from Σ int to Σ f describes modes propagating across the ergoregion encompassing the horizon. This evolution is starkly different for modes with ˜ ω > 0 (non-superradiant modes, NSRM) versus modes with ˜ ω < 0 (superradiant modes, SRM). For NSRM, we have ⃗ A (f) Q = S (NSR) Q · ⃗ A (int) Q , where  \n⃗ A (f) Q ≡ ( a (out) Q , a (down) Q , a (dn) Q , a (out) Q † , a (down) Q † , a (dn) Q † ) ⊤ (5)  \nand S (NSR) Q = ( C Q 0 0 C Q ) , with  \nC Q =   sin θ Q cos θ Q 0 cos θ Q -sin θ Q 0 0 0 1   . (6)  \nThe term cos 2 θ Q represents the transmission coefficient of the potential barrier (i.e., the 'greybody' factor) and depends on ω , ℓ , m and the spin s . Hence, these coefficients differ for electromagnetic and gravitational waves.  \nThe transformation (6) corresponds to a beam splitter . Beam splitters are passive transformations , meaning they do not amplify waves (no particle creation). Their role is simply to divide the incoming waves in transmitted and reflected portions.  \nFor SRM, the evolution across the ergoregion has a completely different character , and is described by  \nS (SR) Q = ( A SR Q B SR Q B SR Q A SR Q ) with  \nA SR Q = ( cosh z erg I 2 0 0 1 ) , (7)  \nB SR Q = ( sinh z erg ¯ I 2 0 0 0 ) . (8)  \nThis is a two-mode squeezing process between modes 'in' and 'up' to produce 'out' and 'down' [14]. The squeezing intensity z erg ('erg' stands for ergoregion) depends on ω , ℓ , m and s . This implies that, for SRM, the ergoregion acts like a quantum amplifier-able to create entangled pairs of particles-very much like the horizon does. Note that the evolution matrices are real; this is possible by exploiting the freedom to introduce a timeindependent phase to the aforementioned creation and annihilation operators.  \nThe transition of the ergoregion-from acting like a beam splitter for NSRM to a two-mode squeezer for SRM-is of vital importance for most results presented in this article. From a physical perspective, this transition arises from the fact that the Killing field describing time translations at infinity becomes spacelike within the ergoregion. Consequently, the modes labeled 'down' having ˜ ω < 0 (although ω > 0) carry negative energy across the horizon from the standpoint of inertial observers at infinity. The appearance of these negative energy modes within the ergoregion provides a mechanism  \nfor energy extraction, effectively transforming the ergoregion into an amplifier. This line of reasoning builds upon the ideas of Penrose and Zeldovich [21-24]. Our goal is to extend their implications into the quantum realm by demonstrating that this mechanism for energy extraction also acts as a source of entanglement.  \nThe complete evolution, from Σ i to Σ f , is described by the matrix S (tot) Q , obtained by simply multiplying the evolution matrices of the two processes. For NSRM, S (tot) Q = S (NSR) Q · S (H) Q ; while for SRM, S (tot) Q = S (SR) Q · S (H) Q . Given a quantum state for electromagnetic and gravitational perturbations at Σ i , S (tot) is all we need to compute observable quantities at Σ f . We have calculated both θ Q and z erg (thereby characterizing S (tot) Q ) by numerically solving Teukolsky's equation in the Kerr geometry for electromagnetic and gravitational waves, following Teukolsky-Press [25] and Page [16, 26]. Some details of our numerical calculations can be found in the Appendix A 6. This is the only place where numerics enters in our analysis (which includes the computation of the angular eigenvalues). The rest, including calculations of entanglement, are obtained from analytic formulas.", 'III. GAUSSIAN STATE EVOLUTION.': "The advantage of restricting to Gaussian states is that they can be described using vectors and matrices in the classical phase space, which are significantly more amenable than density matrices in the (infinite dimensional) Hilbert space. The reason is that Gaussian states are uniquely determined by their first and second moments, since higher order moments can all be derived from them. In field theory, the phase space is infinite dimensional. Though, for our problem, the evolution factorizes in decoupled Q -sectors, reducing the analysis to an effective six-dimensional phase space.  \nLet us fix a value of Q . Given a Gaussian state (pure or mixed) at Σ i , the information in its density matrix ρ Q is encoded in the vector of first moments ⃗ µ (i) Q ≡ Tr[ ρ Q ⃗ A (i) Q ] and the covariance matrix σ (i) Q ≡ Tr[ ρ Q { ⃗ A (i) Q -⃗ µ (i) Q , ⃗ A (i) Q -⃗ µ (i) Q } ], where the curly brackets denote anti-commutators. The covariance matrix contains the symmetric part of the second moments, after subtracting the first moments; the anti-symmetric part is state independent and fully determined by the canonical commutation relations. For Gaussian states, the pair ( ⃗ µ (i) Q , σ (i) Q ) completely specifies the density matrix ρ Q . Linear evolution preserves the Gaussian character of the state and, at Σ f , the state is a Gaussian state fully described by ⃗ µ (f) Q = S (tot) Q · ⃗ µ (i) Q and σ (f) Q = S (tot) Q · σ (i) Q · S (tot) ⊤ Q . From the pair ( ⃗ µ (f) Q , σ (f) Q ) we can obtain particle number, entropies, entanglement, etc. in a simple manner. These techniques are well known (see, e.g., [27, 28]) but a summary is included in the Appendix B.  \nFor the initial state, we consider populating the 'in'  \nFIG. 2. Entanglement rate for the partition separating the interior and exterior modes of the black hole. T H ≡ κ/ (2 π ) is the Hawking temperature.  \n<!-- image -->  \nmode with thermal quanta. The 'p' and 'd' are set in vacuum because they are supported on ultra-high frequencies and, for temperatures well below the Planck scale, it is an excellent approximation to consider them unpopulated. The initial state is mixed and described by  \n⃗ µ (i) Q = ⃗ 0 and σ (i) Q = ( 0 D Q D Q 0 ) , where  \nD Q = ( (1 + 2 n env ω ) 0 0 I 2 ) , (9)  \nand n env ω = ( e ω/T env -1) -1 is the mean number of quanta in the mode 'in' for environment temperature T env .  \nParticle emission. It is straightforward to read from σ (f) Q the mean particle number on each of the modes. We discuss NSRM and SRM separately. For NSRM modes, the number of 'out' quanta is encoded in the component (1 , 4) of the matrix σ (f) Q and given as  \n¯ n (out) Q = sin 2 θ Q n env ω + cos 2 θ Q sinh 2 z H . (10)  \nThe second term accounts for the Hawking quanta, where sinh 2 z H is the thermal part of the spectrum and cos 2 θ Q is the greybody factor. The first term corresponds to the portion of the thermal radiation approaching the black hole from outside that is reflected back to I + .  \nFor SRM, expression (10) is replaced by  \n¯ n (out) Q = n env ω +sinh z 2 erg (1 + n env ω +sinh z 2 H ) . (11)  \nThe first term corresponds to the quanta already present in the initial state, while the second terms is proportional to the quanta generated by the ergoregion-squeezer, sinh z 2 erg . The last contribution describes the interplay of the two squeezing processes. We find that NSRM dominate the emission for a ≲ 0 . 9 M .  \n<!-- image -->  \nFIG. 3. (a) Entanglement rate (quantified by LN per unit time) for the partition between the interior and exterior of the black hole, separated into SRM and NSRM contributions for T env = 0. (b) Entanglement rate for out | down modes (a 'superradiant pair') for various T env .  \n<!-- image -->", 'IV. ENTROPY AND ENTANGLEMENT GENERATION.': "The von Neumann entropy of the 'out' modes can easily be obtained from the covariance matrix σ (f) Q (see Appendix B)  \nS (out) Q = ( ¯ n (out) Q +1 ) ln ( ¯ n (out) Q +1 ) -¯ n (out) Q ln ¯ n (out) Q . (12)  \n̸  \nIf we restrict to a = 0 and to vacuum input ( T env = 0), Eqn. (12) reproduces Page's results [4]. For SRM, the entropy of the Hawking radiation is again given by expression (12), using now Eqn. (11) for ¯ n (out) Q . These expressions extend Page's entropy calculations [4] to the rotating case and for T env = 0.  \nFor mixed quantum states, the entropy of the outgoing radiation is no longer adequate for quantifying entanglement. Instead, we employ the Logarithmic Negativity [29-31], since it is a measure of entanglement that is applicable to both pure and mixed states. Let ρ be the density matrix of a quantum system made of two subsystems, A and B. The Logarithmic Negativity with respect to the partition A | B is defined as  \nLN( ρ ) = log 2 || ρ ⊤ B || 1 , (13)  \nwhere ρ ⊤ B represents the partial transpose of ρ with respect to subsystem B, and || · || 1 denotes the trace norm. A non-zero LN value indicates a violation of the Positivity of the Partial Transpose criterion for quantum states [32]. For Gaussian states, and when one subsystem contains a single mode, LN is non-zero if and only if the state is entangled. Furthermore, LN is as a faithful quantifier of entanglement, signifying that a higher LN value corresponds to a greater degree of entanglement.  \nFrom the covariance matrix σ (f) Q , it is a simple task to compute the Logarithmic Negativity between any bipartition A | B of the three modes (see Appendix B). An  \ninteresting partition is the interior-exterior one, where subsystem A is made of the 'out' mode and subsystem B of the two modes falling in the horizon ('dn' and 'down'). Figure 2 shows the entanglement rate (measured via the Logarithmic Negativity generated per unit of retarded time u ). This plot communicates two important messages: (i) rotation enhances the entanglement of the black hole with the exterior; (ii) ambient thermal radiation inhibits the generation of entanglement, reducing entanglement by ∼ 90%. Intuitively, this is because thermal fluctuations can overwhelm quantum coherence within the system, leading to reduced entanglement. The degradation of entanglement, in turn, illustrates that the entropy of the Hawking mode is very different from entanglement: entropy grows with the environment temperature, while entanglement decreases. It is also interesting to note that, when T env ≫ T H , even when the black hole does not evaporate, it keeps creating entanglement, although at a much slower rate than that in isolation.  \nNext, we want to understand how much of the entanglement is contained in SRM and NSRM. The answer is revealed in Fig. 3. Figure 3(a) shows the entanglement rate for the interior-exterior partition, for T env = 0, together with the contribution from SRM and NSRM. We observe that, as for energy and angular momentum, NSRM dominates for black holes with a ≲ 0 . 9 M . In other words, for a ≲ 0 . 9 M , most of the entanglement generated can be attributed to the horizon.  \nFinally, we want to know what portion of the entanglement in SRM can be attributed to the ergoregion squeezer. This is not obvious a priori because SRM experience both squeezers. For this, we compute the entanglement between the 'out' mode at I + and the 'down' mode at H . This entanglement cannot be generated by the horizon-squeezer; hence, it is a measure of the contribution of the ergoregion to the final entanglement.  \nFigure 3(b) shows that entanglement between 'out' and 'down' modes is very small for a ≲ 0 . 9 M . On the contrary, for a > 0 . 9 M , the entanglement between the 'out' and 'down' pair is of the same order as the total entanglement between exterior and interior of the black hole, confirming that the ergoregion plays a significant role in generating entanglement. Notice that the previous plots report entanglement carried out by both photons and gravitons together. We find that gravitons dominate the generation of entanglement for a > 0 . 7 M (see Fig. 4). This aligns with the findings on energy emission in [16].  \nFIG. 4. Entanglement rate (quantified by LN per unit time) for the partition between the interior and exterior of the black hole at T env = 0. Similar behavior holds for T env > 0, albeit with parametrically suppressed values of LN. We show separately the contribution from photons and gravitons.  \n<!-- image -->", 'V. DISCUSSION': "When the initial state of the quantum fields is assumed to be the vacuum, it suffices to look at the outgoing Hawking radiation to understand the entanglement generated by an evaporating black hole. For the goals of this article, however, a detailed characterization of the final state of modes that fall inside the horizon-and their individual correlation with Hawking modes-is desired. The Gaussian formalism employed here provides a remarkably good balance between simplicity and power, from which we have been able to quantify the impact of ambient thermal radiation and the contribution of the ergoregion to the process of Hawking radiation and entanglement generation.  \nThe input to the horizon squeezer is vacuum fluctuations, whereas the ergoregion is seeded by the thermal Hawking quanta emanating from the horizon. This in-  \nput stimulates the creation of quanta at the ergoregion, in a similar vein to 'classical' superradiant phenomena. However, the thermal character of the input from the horizon partially inhibits the generation of entanglement by the ergoregion alone. It is only for rapidly spinning black holes ( a > 0 . 9 M ) that the ergoregion is able to contribute significantly to the entanglement between the black hole and its exterior, completely dominating the emission in the extremal limit ( a → 1 M ).  \nThe results of this article have significant ramifications. In the realm of black hole physics, our results emphasize that radiation enhancement by the ergoregion of a Kerr black hole is not a purely classical effect, even when the ergoregion is stimulated with classical thermal radiation. We have proven this by showing that the ergoregion is an active source of entanglement.  \nOur findings are also important in addressing information-related questions, such as those pertaining to generalized Page curves for rotating black holes in thermal baths. For instance, by employing the techniques outlined here, one can monitor the entanglement between the escaping radiation and the interior of the rotating hole during the course of the black hole's lifetime.  \nFurthermore, our results can be extended to other types of ergoregions, broadening the implications beyond black hole physics. This presents an exciting opportunity for validating our predictions in analogue models, where radiation from ergoregions can actually be observed [33]. These controllable systems are appealing because they grant access to all modes, in contrast to real black holes. Quantum fluid analogues [34, 35] are particularly promising platforms that provide a low-noise environment amenable to correlation measurements.", 'ACKNOWLEDGMENTS': "This article has benefited from discussion with A. Ashtekar, M. Jacquet and L. Giacomelli. I.A. and D.K. are supported by the NSF grant PHY-2110273, by the RCS program of Louisiana Boards of Regents through the grant LEQSF(2023-25)-RD-A-04, and by the Hearne Institute for Theoretical Physics. A.D. is supported by the NSF grant PHY-1903799, PHY-2206557, the Hearne Institute for Theoretical Physics, the European Regional Development Fund through the Center of Excellence TK133 'The Dark Side of the Universe', and by the Estonian Research Council through the grant PRG356. I.A.'s research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science, and Economic Development, and by the Province of Ontario through the Ministry of Colleges and Universities. A.J.B. is supported by the DARPA Young Faculty Award (YFA) Grant No. N660012014029.  \n- [1] Hawking, Stephen W, Black hole explosions?, Nature 248 , 30 (1974).\n- [2] S. W. Hawking, Particle Creation by Black Holes, Commun. Math. Phys. 43 , 199 (1975), [Erratum: Commun. Math. Phys. 46, 206 (1976)].\n- [3] D. N. Page, Information in black hole radiation, Phys. Rev. Lett. 71 , 3743 (1993).\n- [4] D. N. Page, Time dependence of Hawking radiation entropy, JCAP 2013 (09), 028.\n- [5] W. H. Zurek, Entropy evaporated by a black hole, Phys. Rev. Lett. 49 , 1683 (1982).\n- [6] S. W. Hawking, Breakdown of Predictability in Gravitational Collapse, Phys. Rev. D 14 , 2460 (1976).\n- [7] Unruh, William G and Wald, Robert M, Information loss, Rep. Prog. Phys. 80 , 092002 (2017).\n- [8] D. Marolf, The black hole information problem: past, present, and future, Rep. Prog. Phys. 80 , 092001 (2017).\n- [9] A. Ashtekar, Black Hole Evaporation: A Perspective from Loop Quantum Gravity, Universe 6 , 10.3390/universe6020021 (2020).\n- [10] A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian, and A. Tajdini, The entropy of Hawking radiation, Rev. Mod. Phys. 93 , 035002 (2021).\n- [11] R. M. Wald, On Particle Creation by Black Holes, Commun. Math. Phys. 45 , 9 (1975).\n- [12] A. A. Starobinsky, Amplification of waves reflected from a rotating 'black hole', Sov. Phys. JETP 37 , 28 (1973).\n- [13] W. G. Unruh, Second quantization in the Kerr metric, Phys. Rev. D 10 , 3194 (1974).\n- [14] A. L. Matacz, P. C. W. Davies, and A. C. Ottewill, Quantum vacuum instability near rotating stars, Phys. Rev. D 47 , 1557 (1993).\n- [15] A. Ashtekar and A. Magnon, C. R. Acad. Sci. Ser. A 281 , 875 (1975).\n- [16] D. N. Page, Particle Emission Rates from a Black Hole. 2. Massless Particles from a Rotating Hole, Phys. Rev. D 14 , 3260 (1976).\n- [17] A. J. Brady, I. Agullo, and D. Kranas, Symplectic circuits, entanglement, and stimulated Hawking radiation in analogue gravity, Phys. Rev. D 106 , 105021 (2022).\n- [18] I. Agullo, A. J. Brady, and D. Kranas, Quantum Aspects of Stimulated Hawking Radiation in an Optical Analog White-Black Hole Pair, Phys. Rev. Lett. 128 , 091301 (2022).\n- [19] I. Agullo, A. J. Brady, and D. Kranas, Event horizons are tunable factories of quantum entanglement, J. Mod. Phys. D 31 , 2242008 (2022).\n- [20] V. Frolov and I. Novikov, Black hole physics: Basic concepts and new developments , Vol. 96 (Springer Science & Business Media, 2012).\n- [21] R. Penrose, Gravitational collapse: The role of general relativity, Riv. Nuovo Cim. 1 , 252 (1969).\n- [22] R. Penrose and R. M. Floyd, Extraction of rotational energy from a black hole, Nature 229 , 177 (1971).\n- [23] Y. B. Zel'dovich, Zh.Eksp.Teor.Fiz. 14 , 270 (1971).\n- [24] Y. B. Zel'dovich, Zh.Eksp.Teor.Fiz. 62 , 2076 (1972).\n- [25] S. A. Teukolsky and W. H. Press, Perturbations of a rotating black hole. III - Interaction of the hole with gravitational and electromagnet ic radiation, Astrophys. J. 193 , 443 (1974).\n- [26] D. N. Page, Particle Emission Rates from a Black Hole:  \n- Massless Particles from an Uncharged, Nonrotating Hole, Phys. Rev. D 13 , 198 (1976).\n- [27] C. Weedbrook, S. Pirandola, R. Garc'ıa-Patr'on, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Rev. Mod. Phys. 84 , 621 (2012).\n- [28] A. Serafini, Quantum continuous variables: a primer of theoretical methods (CRC press, 2017).\n- [29] G. Vidal and R. F. Werner, Computable measure of entanglement, Physical Review A 65 , 032314 (2002).\n- [30] M. B. Plenio, Logarithmic negativity: A full entanglement monotone that is not convex, Phys. Rev. Lett. 95 , 090503 (2005).\n- [31] R. F. Werner and M. M. Wolf, Bound entangled gaussian states, Phys. Rev. Lett. 86 , 3658 (2001).\n- [32] A. Peres, Separability Criterion for Density Matrices, Phys. Rev. Lett. 77 , 1413 (1996).\n- [33] T. Torres, S. Patrick, A. Coutant, M. Richartz, E. W. Tedford, and S. Weinfurtner, Rotational superradiant scattering in a vortex flow, Nature Physics 13 , 833 (2017).\n- [34] F. Claude, M. J. Jacquet, R. Usciati, I. Carusotto, E. Giacobino, A. Bramati, and Q. Glorieux, HighResolution Coherent Probe Spectroscopy of a Polariton Quantum Fluid, Phys. Rev. Lett. 129 , 103601 (2022), arXiv:2112.09903 [cond-mat.quant-gas].\n- [35] M. J. Jacquet, M. Joly, L. Giacomelli, F. Claude, Q. Glorieux, A. Bramati, I. Carusotto, and E. Giacobino, Analogue quantum simulation of the Hawking effect in a polariton superfluid, Eur. Phys. J. D 76 , 152 (2022), [Erratum: Eur.Phys.J.D 76, 247 (2022)], arXiv:2201.02038 [quant-ph].\n- [36] S. W. Hawking and G. F. R. Ellis, The Large Scale Structure of Space-Time , Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2023).\n- [37] Fabbri, Alessandro and Navarro-Salas, Jos'e, Modeling black hole evaporation (World Scientific, 2005).\n- [38] S. A. Teukolsky, Perturbations of a rotating black hole. 1. Fundamental equations for gravitational electromagnetic and neutrino field perturbations, Astrophys. J. 185 , 635 (1973).\n- [39] E. Berti, V. Cardoso, and M. Casals, Eigenvalues and eigenfunctions of spin-weighted spheroidal harmonics in four and higher dimensions, Phys. Rev. D 73 , 024013 (2006), [Erratum: Phys.Rev.D 73, 109902 (2006)], arXiv:gr-qc/0511111.\n- [40] P. L. Chrzanowski, Vector potential and metric perturbations of a rotating black hole, Phys. Rev. D 11 , 2042 (1975).\n- [41] R. M. Wald, Construction of Solutions of Gravitational, Electromagnetic, Or Other Perturbation Equations from Solutions of Decoupled Equations, Phys. Rev. Lett. 41 , 203 (1978).\n- [42] P. Candelas, P. Chrzanowski, and K. W. Howard, Quantization of Electromagnetic and Gravitational Perturbations of a Kerr Black Hole, Phys. Rev. D 24 , 297 (1981).\n- [43] D. V. Gal'tsov, Radiation reaction in the Kerr gravitational field, J. Phys. A 15 , 3737 (1982).\n- [44] D. V. Gal'tsov, Particles and fields in the vicinities of the black holes (Moscow University Press, 1986).\n- [45] J. A. H. Futterman, F. A. Handler, and R. A. Matzner, SCATTERING FROM BLACK HOLES , Cambridge  \nMonographs on Mathematical Physics (Cambridge University Press, 2012).  \n- [46] I. Agullo, J. Navarro-Salas, G. J. Olmo, and L. Parker, Hawking radiation by Kerr black holes and conformal symmetry, Phys. Rev. Lett. 105 , 211305 (2010), arXiv:1006.4404 [hep-th].\n- [47] We thank Adrian del Rio for useful guidance on this calculation in the early stages of this project.  \n- [48] http://blackholes.ist.utl.pt/?page=Files https: //web.uniroma1.it/gmunu https://www.darkgra.org/ .\n- [49] R. Simon, Peres-Horodecki Separability Criterion for Continuous Variable Systems, Phys. Rev. Lett. 84 , 2726 (2000).\n- [50] X. Wang and M. M. Wilde, Cost of Quantum Entanglement Simplified, Phys. Rev. Lett. 125 , 040502 (2020).\n- [51] X. Wang and M. M. Wilde, α -logarithmic negativity, Phys. Rev. A 102 , 032416 (2020).", '1. Kerr geometry': "We collect here a few expressions which have been used in the manuscript. A detailed account of Kerr's geometry can be found, for instance, in [20, 36, 37].  \nThe Kerr line element in Boyer-Lindquist coordinates takes the form  \nds 2 = -∆ ρ 2 ( dt -a sin 2 θ dϕ ) 2 + sin 2 θ ρ 2 [( r 2 + a 2 ) dϕ -adt ] 2 + ρ 2 ∆ dr 2 + ρdθ 2 , (A1)  \nwhere a is the angular momentum of the black hole (BH) per unit mass and ρ 2 ( r ) = r 2 + a 2 cos 2 θ , ∆ = ( r + r + )( r -r -), where  \nr ± = M ± √ M 2 -a 2 , (A2)  \nare the locations of the outer and inner horizons, respectively. The surface gravity and angular velocity of the horizon of a Kerr BH are  \nκ = r + -r -2( r 2 + + a 2 ) and Ω H = a r 2 + + a 2 . (A3)  \nFor a = M we have an extremal BH, whose surface gravity vanishes, and for a/M > 1 the Kerr line element contains a naked singularity. For this reason, we restrict to a/M ≤ 1. The line element (A1) is degenerate at r = r ± , but these are unphysical coordinate singularities that can be resolved by a change of coordinates.  \nThe Kerr line element is stationary and axi-symmetric, meaning it has two independent Killing vector fields (KVF), ∂ t and ∂ ϕ . The vector field ∂ t is time-like for r > r 0 ( θ ) = M + √ M 2 -a 2 cos 2 θ , null at r = r 0 ( θ ) and space-like in the interval r 0 ( θ ) > r > r + . This latter region is called the ergoregion, and it plays an important role in this article since it gives rise to the phenomenon of superradiance. From the point of view of a static observer at infinity, this phenomenon is a result of the fact that ∂ t is not time-like inside the ergoregion. Indeed, the KVF that becomes null at the horizon (at r = r + ) is the combination ∂ t +Ω H ∂ ϕ . This vector field is timelike inside the ergoregion, and it is also the generator of the event horizon. Hence, this KVF can be use to define a natural notion of frequency, ˜ ω , for modes propagating near the outer horizon.  \nFinally, some coordinates which are useful at some points in our calculations are the tortoise coordinate r ⋆ , and the angular coordinate ˜ ϕ defined as  \ndr ⋆ = r 2 + a 2 ∆ dr and ˜ ϕ = ϕ -Ω H t. (A4)  \nThe former has the effect of pushing the horizon to r ⋆ →-∞ . The latter can be thought of as an angular coordinate adapted to the rotation of the spacetime at the event horizon.", '2. Perturbations in a Kerr spacetime: Some facts about Teukolsky equation': "As shown by Teukolsky [38], perturbations of minimally coupled massless fields of spin | s | in a Kerr geometry can be described by a general scalar equation, known as Teukolsky's equation, which, quite remarkably, admits separation  \nof variables as  \ns ϕ ( t, r, ϕ, θ ) = ∑ ℓ,m ˆ dω s R ωℓm ( r ) e -i ω t s Z ω ℓm ( θ, ϕ ) , (A5)  \nwhere s = 0 , ± 1 / 2 , ± 1 , ± 3 / 2 , ± 2 , ... . The angular dependence is described by spin-weighted spheroidal harmonics s Z ω ℓm ( θ, ϕ ) (see e.g. [20]). They satisfy the orthogonality relations  \nˆ d Ω s ¯ Z ω ℓm s Z ω ℓ ' m ' = δ ℓ ' ℓ δ mm ' , (A6)  \nwhere d Ω is the standard integration measure on the sphere and we use a bar to indicate complex conjugation. The spheroidal harmonics can be defined to satisfy  \ns ¯ Z ω ℓm ( θ, ϕ ) = ( -1) m s Z -ω ℓ -m ( θ, ϕ ) . (A7)  \nThe ϕ dependence of s Z ω ℓ,m ( θ, ϕ ) is given by a phase e imϕ , and in the limit of zero rotation, s Z ω ℓ,m ( θ, ϕ ) reduces to the standard spin-weighted spherical harmonics (which reduce to the spherical harmonics for s=0). The spin-weighted spheroidal harmonics are solutions of a Sturm-Liouville problem. The corresponding eigenvalues s E ωℓm are known only numerically, although analytical approximations exist in some regimes (see e.g. [39]).  \nThe radial part of Teukolsky equation can be written as follows. Define a new radial variable  \ns χ ωℓm ≡ ( r 2 + a 2 ) 1 / 2 ∆ s/ 2 s R ωℓm ( r ) . (A8)  \nIn terms of this variable, the radial equation takes the form  \nd 2 dr 2 ⋆ s χ ωℓm + s V ωℓm ( r ) s χ ωℓm = 0 , (A9)  \nwhere r ⋆ is the tortoise coordinate defined in (A4), and the effective potential s V ωℓm has the form  \ns V ωℓm ( r ) = K 2 -2 is ( r -M ) K +∆(4 isωr -s λ ωℓm ) ( r 2 + a 2 ) 2 -G 2 -dG dr ⋆ , (A10)  \nwhere  \ns  \nK = ( r 2 + a 2 ) ω -am, (A11) G = r ∆ ( r 2 + a 2 ) 2 + s ( r -M ) r 2 + a 2 ,  \nλ  \n=  \nE  \n-  \n2  \namω  \n+  \na  \n2  \nω  \n2  \n-  \ns  \n(  \ns  \n+1)  \n.  \nωℓm  \ns  \nωℓm  \n(A12)  \nAlthough the form of this potential is complicated, it simplifies in the asymptotic regimes  \ns V ωℓm ( r ) ∼ { ω 2 + 2 isω r r ⋆ →∞ , ˜ ω 2 -2 is ˜ ωκ -4 s 2 κ 2 r ⋆ →-∞ , (A13)  \nwhere ˜ ω ≡ ω -m Ω H and κ is the surface gravity given in (A3). This implies that linearly independent solutions behave asymptotically as r ± s e ± i ω r ⋆ as r ⋆ →∞ , and ∆ ± s/ 2 e ± i ˜ ωr ⋆ in the near horizon region r ⋆ →-∞ .", '3. Modes for scalar waves': "The modes that are relevant to understand the Hawking effect are similar for all spins. Hence, the discussion of these modes and their properties can be done for scalar fields, which simplifies some technicalities which are inessential for our main purposes. Hence, for pedagogical purposes, we discuss first scalar modes and their properties, and discuss later how to generalise to higher spin. Most of the material that is summarized in this section can be found in greater detail in [20].  \nIn the following, we write the explicit form of the modes used in the main body of the paper in the vicinity of the Cauchy hypersurface where they are defined. An important aspect of these solutions is the sign of their norm, defined  \nfrom the so-called Klein-Gordon, or symplectic product. This product can be defined on any Cauchy hypersurface, and it is preserved in time for scalar fields satisfying the Klein-Gordon equation □ g ϕ = 0.  \nIt is given by  \n⟨ φ 1 , φ 2 ⟩ = i ˆ Σ d Σ µ ( ¯ φ 1 ∂ µ φ 2 -φ 2 ∂ µ ¯ φ 1 ) (A14)  \nwhere d Σ µ is the oriented volume element of the Cauchy hypersurface Σ.  \na. The 'out' modes  \nThese are the modes carrying Hawking radiation to future null infinity. At the Cauchy hypersurface Σ f , these modes have support only at I + , and oscillate as e -iωu . More concretely, their form in the vicinity of I + is  \nφ (out) Q ∣ ∣ ∣ I + ∼ 1 √ 4 πω e -iωu √ r 2 + a 2 0 Z ω ℓm ( θ, ϕ ) , (A15)  \nwhere Q = ( ω, ℓ, m ). These modes have positive norm for ω > 0. The 'out' modes used in the main text are wave packets constructed from them, centered at late times uκ ≫ 1, when Hawking radiation reaches I + . We can define normalized wave packets following Hawking's original derivation [2]. For integers j ≥ 0 and n  \nφ (out) jnℓm = 1 √ ϵ ˆ ( j +1) ϵ jϵ e 2 πinω/ϵ φ (out) ωℓm dω. (A16)  \nThese wave packets are composed of frequencies in the interval [ jϵ, ( j +1) ϵ ], and are peaked at retarded time u = 2 πn/ϵ with width 2 π/ϵ . Wave packets for other modes will be defined similarly. To avoid introducing new notation, we will keep using the symbol φ (out) Q to denote these wave packets, where the label Q is now a shorthand for Q = ( j, n, ℓ, m ).  \nb. The 'down' modes  \nThese are modes that enter the horizon at late times. In the vicinity of the Cauchy hypersurface Σ f , these modes are defined to have support only at the horizon (no support at I + ), and to oscillate as e -i ˜ ωv . In more detail, they are defined as follows. First, consider the (normalized) plane waves  \nf (down) Q ∣ ∣ ∣ H = 1 √ 4 π ˜ ω e -i ˜ ωv √ r 2 + + a 2 0 Z ω ℓm ( θ, ˜ ϕ ) , (A17)  \nwhich have a norm with a sign given by sign(˜ ω ). From these functions, we define the down modes as  \nφ (down) Q ∣ ∣ ∣ H ≡ { f (down) Q ω > 0 , ˜ ω > 0 ¯ f (down) Q ω > 0 , ˜ ω < 0 . (A18)  \nThese have a positive norm for ω > 0. (This is due to having defined them from the complex conjugate modes ¯ f (down) when ˜ ω < 0.) The 'down' modes used in the main text, and denoted also as φ (down) Q , are normalized wave packets constructed from (A18), peaked at late times.  \nc. The 'in' modes  \nThe 'in' modes are modes that propagate on the stationary part of the spacetime, once the BH geometry is described by the Kerr metric. They represent radiation approaching the BH from infinity. At Σ i , they take the form  \nφ (in) Q ∣ ∣ ∣ I -∼ 1 √ 4 πω e -iωv √ r 2 + a 2 0 Z ω ℓm ( θ, ϕ ) . (A19)  \nThey have a positive norm for ω > 0. The 'in' modes used in the main text are normalized wave packets constructed from (A19), peaked at late times vκ ≫ 1 along I -.  \nThe 'up' modes represent outgoing radiation trying to escape the BH. At Σ int , these modes have support only at the v =constant part of Σ int , i.e., they have no support at the horizon or at the I -portions of Σ int . They oscillate as e -i ˜ ωu , and can be defined as follows. First, we define plane waves modes as  \nf (up) Q ∣ ∣ ∣ Σ int = 1 √ 4 π ˜ ω e -i ˜ ωu √ r 2 + a 2 0 Z ω ℓm ( θ, ˜ ϕ ) . (A20)  \nThese functions have a norm with a sign given by sign(˜ ω ). In analogy to down modes, using these plane waves we define  \nφ (up) Q ∣ ∣ ∣ Σ int ≡ { f (up) Q ω > 0 , ˜ ω > 0 ¯ f (up) Q ω > 0 , ˜ ω < 0 , (A21)  \nwhich have a positive norm for ω > 0. The 'up' modes used in the main text, denoted also by φ (up) Q , are normalized wave packets constructed from these modes, and supported at late retarded times uκ ≫ 1 (close to the horizon).", "e. The 'dn' modes": "These are the partner modes of the 'up' modes. The term 'partner' refers to the following. In the spontaneous Hawking radiation, 'dn' modes are the modes that purify the thermal Hawking radiation carried in the 'up' modes, in the sense that they are entangled with each other and with nobody else. One can find the 'dn' modes as follows. Propagate the 'up' modes until I -, and decompose the result in its positive and negative frequency contributions (frequency defined by the coordinate v along I -)  \nφ (up) Q ∣ ∣ ∣ I -= φ (+) Q + φ ( -) Q , (A22)  \nwhere we have denoted the positive and negative-frequency parts as φ ( ± ) Q , respectively. Let us denote by α 2 Q and -β 2 Q the norms of φ ( ± ) Q , and let us define from them positive-norm modes  \nφ (p) Q ≡ φ (+) Q α Q and φ (d) Q ≡ ¯ φ ( -) Q β Q , (A23)  \nwhere α Q and β Q can be defined real. Bringing these definitions to (A22), we obtain φ (up) Q = α Q φ (p) Q + β Q ¯ φ (d) Q . From this, the partner modes of 'up' is  \nφ (dn) Q = α Q φ (d) Q + β Q ¯ φ (p) Q . (A24)  \nThe form of this mode was first obtained in [11], by reflecting the mode 'up' about v = v H at I -, followed by complex conjugation (the name 'dn' is chosen because it is the reflection of 'up'). The propagation of f (up) Q to I -produces  \nf (up) Q ∣ ∣ ∣ I -∼ -1 √ 4 π ˜ ω e -i ˜ ω [ v H -1 κ ln κ ( v H -v )] √ r 2 + a 2 0 Z ω ℓm ( θ, ˜ ϕ ) Θ( v H -v ) , (A25)  \nwhere v H is the advanced time at which the horizon forms and Θ is the Heaviside function. From this, one finds [11, 20]  \nf (dn) Q ∣ ∣ ∣ I -∼ -1 √ 4 π ˜ ω e i ˜ ω [ v H -1 κ ln κ ( v -v H )] √ r 2 + a 2 0 ¯ Z ω ℓ,m ( θ, ˜ ϕ ) Θ( v -v H ) . (A26)  \nFrom f (dn) Q , we define the 'dn' plane wave modes as  \nφ (dn) Q ∣ ∣ ∣ Σ int = { f (dn) Q ω > 0 , ˜ ω > 0 ¯ f (dn) Q ω > 0 , ˜ ω < 0 . , (A27)  \nAt the Cauchy hypersurface Σ int , these modes have support only at the horizon, and oscillate as e -i ˜ ωv . They have a positive norm for ω > 0. From these plane wave modes, we define the 'dn' normalized wave packets used in the main text, which are supported near v H .  \nThese modes are the progenitors of the modes 'up' and 'dn', and can be defined from them inverting the expressions written above:  \nφ (p) Q = α Q φ (up) Q + β Q ¯ φ (dn) Q , (A28) φ (d) Q = α Q φ (dn) Q + β Q ¯ φ (up) Q ,  \nwhere α Q = (1 -e 2 π | ˜ ω | /κ ) -1 / 2 and β Q = e -π | ˜ ω | /κ (1 -e 2 π | ˜ ω | /κ ) -1 / 2 . In the main text, we have written these coefficients as α Q = cosh z H and β Q = sinh z H , with z H ( ω, m ) = tanh -1 e -| ˜ ω | / (2 κ ) . Both modes 'p' and 'd' have positive unit norm for ω > 0. As mentioned above, both are made of a combination of positive frequency modes with respect to the coordinate v in I -, with no contribution from negative frequency modes. This can be checked directly by computing their Fourier transform in the coordinate v [11, 20]. This automatically implies that these modes are unentangled when the field is in the vacuum state at I -. In that case, the reduced quantum state describing these two modes is a pure state. This in turn proves that the state describing the pair 'up' and 'dn' is also pure -when the initial state is the vacuum- since evolution from 'p' and 'd' to produce 'up' and 'dn' is unitary. Consequently, the 'dn' modes are the modes that purify 'up'-i.e., these modes are the partners.  \nIn summary, we have constructed an orthonormal wave packet mode basis of the sector of the vector space of solutions to the Klein-Gordon equation that is involved in the Hawking process, at each of the relevant Cauchy hypersurfaces: Σ i , Σ int , and Σ f . These basis elements are labeled by Q , with ω > 0, and can be collected together in (row) vectors as  \n⃗ Ψ (i) Q = ( φ (in) Q , φ (p) Q , φ (d) Q , ¯ φ (in) Q , ¯ φ (p) Q , ¯ φ (d) Q ) ,  \n⃗ Ψ (int) Q = ( φ (in) Q , φ (up) Q , φ (dn) Q , ¯ φ (in) Q , ¯ φ (up) Q , ¯ φ (dn) Q ) ,  \n⃗ Ψ (f) Q = ( φ (out) Q , φ (down) Q , φ (dn) Q , ¯ φ (out) Q , ¯ φ (down) Q , ¯ φ (dn) Q ) .  \nThese are called Wald bases [11, 20]. For each Q , the elements of these bases are the modes that are involved in the Hawking effect. Evolution does not mix modes labeled with different Q . It is straightforward to extend this set of modes to have a complete basis in the space of solutions. However, since such an extension is irrelevant to describe the Hawking effect, we do not write it down.", '4. Evolution matrices': "In this section, we will implicitly exploit the freedom to introduce a time-independent phase to the aforementioned modes, to make the evolution matrices real.  \nThe relation between wave packets ⃗ Ψ (i) Q and ⃗ Ψ (int) Q encodes the dynamics from Σ i to Σ int , and is given by  \n⃗ Ψ (i) Q = ⃗ Ψ (int) Q · S (H) Q , where S (H) Q = ( A Q B Q B Q A Q ) ,  \nwith  \nA Q =   1 0 0 0 cosh z H 0 0 0 cosh z H   and B Q =   0 0 0 0 0 sinh z H 0 sinh z H 0   ,  \nwhere z H ( ω, m ) = tanh -1 e -π | ˜ ω | /κ . These matrices can be read off from the definitions of 'up' and 'dn' modes in terms of 'p' and 'd' modes, Eqn. (A28), adding the identity transformation for the 'in' modes. Subsequently, the evolution from Σ int to Σ f describes modes propagating across the ergoregion surrounding the horizon. This evolution is qualitatively different for modes with ˜ ω > 0 (non-superradiant modes, NSRM) and modes with ˜ ω < 0 (superradiant modes, SRM). For NSRM, the evolution is given by,  \n⃗ Ψ (int) Q = ⃗ Ψ (f) Q · S (NSR) Q , where S (NSR) Q = ( C Q 0 0 C Q ) , (A30)  \n(A29)  \nwith  \nC Q =   sin θ Q cos θ Q 0 cos θ Q -sin θ Q 0 0 0 1   , (A31)  \nwhere cos 2 θ Q is the transmission coefficient of the potential barrier (the greybody factor), and is obtained by solving numerically the radial Teukolsky equation (A9). θ Q depends on ω , ℓ , m , and the spin of the perturbation | s | . Therefore, it is different for photons and gravitons. For SRM, the evolution across the ergoregion has a completely different character and is described by  \n⃗ Ψ (int) Q = ⃗ Ψ (f) Q · S (SR) Q , where S (SR) Q = ( A SR Q B SR Q B SR Q A SR Q ) , (A32)  \nwith  \nA SR J =   cosh z erg 0 0 0 cosh z erg 0 0 0 1   and B SR J =   0 sinh z erg 0 sinh z erg 0 0 0 0 0   . (A33)  \nThe matrix S (SR) Q represents a two-mode squeezer between modes 'in' and 'up' to produce 'out' and 'down', and the identity evolution for the mode 'dn'. The squeezing intensity z erg depends on ω , ℓ , m , and the spin of the perturbation | s | , and is also obtained by solving numerically the radial Teukolsky equation (A9). This implies that, for SRM, the ergoregion acts like a quantum amplifier, creating entangled pairs of particles, very much like the horizon does.  \nMathematically, the distinction between the evolution matrix for SRM and NSRM arises from the fact that the Killing field describing time translations at infinity becomes spacelike within the ergoregion. For superradiant modes with ˜ ω < 0 but ω > 0, the change in the character of ∂ t reverses the sign of the norm of the modes labeled as 'up' and 'down' that interact with modes 'in' and 'out'. The latter two modes always have positive norm for ω > 0. This, in turn, implies that, for superradiant frequencies, the evolution from Σ int to Σ f mixes modes with different signs of their norm. At the quantum level, this translates into Bogoliubov transformations that mix creation and annihilation operators, corresponding to a quantum amplifier. This phenomenon does not occur for NSRM, as all modes with ω > 0 have positive norm.  \nThis fact is also reflected in that, while both evolution matrices S (SR) Q and S (NSR) Q belong to the symplectic group, only S (NSR) Q belongs to the orthogonal subgroup of the symplectic group; i.e. S (NSR) Q · S (NSR) ⊤ Q = I 6 . This automatically implies that the transformation S (NSR) Q leaves the vacuum and the particle number invariant. This is not the case for the scattering matrix of superradiant modes S (SR) Q .  \nWe can now write the total evolution matrix describing the evolution from Σ i to Σ f . It is given by  \n⃗ Ψ (i) Q = ⃗ Ψ (f) Q · S (tot) Q , (A34)  \nwhere, for non-superradiant modes  \nwhile for SRM  \nS (tot) Q = S (SR) Q · S (H) Q . (A36)  \nSince evolution factorizes in decoupled Q sectors, this is enough to compute the evolution for the field. As a consistency check, one can verify that the total evolution matrices are both linear canonical transformations, meaning that they leave invariant the symplectic structure of the phase space:  \nS (tot) Q · Ω 6 · S (tot) ⊤ Q = Ω 6 , where Ω 6 = ( 0 3 I 3 -I 3 0 3 . ) . (A37)  \nThis is equivalent to saying that S (tot) Q belongs to the symplectic group Sp( C , 6). Symplectic matrices have the general form  \nS (tot) Q = ( A B ¯ B ¯ A ) , (A38)  \nS (tot) Q = S (NSR) Q · S (H) Q , (A35)  \nwhere A and B are 3 × 3 matrices and the bar denotes complex conjugation of their elements. The condition (A37) imposes the following constraints on the components of these matrices  \nA · B ⊤ -B · A ⊤ = 0 3 , (A39)  \nA · A † -B · B † = I 3 . (A40)  \nThe components of A and B are usually called Bogoliubov coefficients, and these two equations are the familiar constraints satisfied by them.  \nIn the quantum theory, we can write the evolution in terms of creation and annihilation operators. First, let us define the column vectors of creation and annihilation operators defined from modes ⃗ Ψ (i) Q , ⃗ Ψ (int) Q and ⃗ Ψ (f) Q , respectively  \n⃗ A (i) Q = ( a (in) Q , a (p) Q , a (d) Q , a (in) † Q , a (p) † Q , a (d) † J ) ⊤ ,  \n⃗ A (int) Q = ( a (in) Q , a (up) Q , a (dn) Q , a (in) † Q , a (up) † Q , a (dn) † Q ) ⊤ ,  \n⃗ A (f) Q = ( a (out) Q , a (down) Q , a (dn) Q , a (out) † Q , a (down) † Q , a (dn) † Q ) ⊤ .  \nThe relation between 'i' and 'f' creation and annihilation operators, can be derived from the relation between the basis modes ⃗ Ψ (i) Q and ⃗ Ψ (f) Q using the following argument. The field operator ϕ ( x ) can be expanded in any of these two sets of modes  \nϕ ( x ) = ∑ Q ⃗ Ψ (i) Q ( x ) · ⃗ A (i) Q = ∑ Q ⃗ Ψ (f) Q ( x ) · ⃗ A (f) Q . (A41)  \nIf we start from the 'i' representation, ϕ ( x ) = ∑ Q ⃗ Ψ (i) Q ( x ) · ⃗ A (i) Q and replace ⃗ Ψ (i) Q = ⃗ Ψ (f) Q ( x ) · S (tot) Q , we find  \nϕ ( x ) = ∑ Q ⃗ Ψ (f) Q ( x ) · S (tot) Q · ⃗ A (i) Q .  \nComparing this result with the last part of (A41), we conclude  \n⃗ A (f) Q = S (tot) Q · ⃗ A (i) Q .  \nThe commutation relation satisfied by ⃗ A (i) Q can be compactly written in vector notation as  \n[ ⃗ A (i) Q , ⃗ A (i) Q ] = Ω 6 . (A42)  \nThe fact that S (tot) Q is a symplectic transformation, guarantees that the operators ⃗ A (f) Q satisfy the same commutation relations:  \n[ ⃗ A (f) Q , ⃗ A (f) Q ] = [ S (tot) J · ⃗ A (i) Q , S (tot) Q · ⃗ A (i) Q ] = S (tot) Q · [ ⃗ A (i) Q , ⃗ A (i) Q ] ︸ ︷︷ ︸ = Ω 6 · S (tot) ⊤ Q = Ω 6 (A43)", '5. Electromagnetic and gravitational waves': "The strategy to construct the Wald basis at each of the Cauchy hypersurfaces for massless higher spin fields is conceptually similar to the scalar case. The technical details are, however, more involved, in part due to the presence of polarization tensors and gauge redundancies. Fortunately, this problem was solved already in the seventies and early eighties in [40-45]. A complete summary, containing many useful details, can be found in Appendix G of [20].  \nElectromagnetic and gravitational perturbations are also described by the Teukolsky equation for | s | = 1 and | s | = 2, respectively. A basis of solutions to this equation can again be labeled by ω, ℓ, m, s , where the sign of s denotes the two possible polarization states (helicity). These modes are described by covariant vectors a µ ( ω, ℓ, m, s, x ) for electromagnetic fields, and symmetric covariant tensors h µν ( ω, ℓ, m, s, x ) for gravitational linear perturbations. In  \nan analogous way as done in the previous section, one defines modes 'out', 'in', 'down', 'up', 'dn', 'p' and 'd'. In Kerr geometry, explicit expressions for these modes are known. These expressions are, however, complicated and not very illuminating. Furthermore, they are not necessary for our goals. Rather, the important aspect of these solutions is that their u , v and ϕ dependence on the three Cauchy hypersurfaces Σ f , Σ int and Σ i is the same as for the scalar case. In particular, the evolution matrix describing the propagation from Σ i to Σ int , is precisely the same for electromagnetic and gravitational perturbations as for scalar fields. Heuristically, this can be understood as a consequence of the exponentially large blueshift near the event horizon, which makes conformal invariance to emerge in the near horizon region [46]. The evolution from Σ int to Σ f , in contrast, depends on the spin of the field, and is dictated by the Teukolsky equation.  \nFollowing Teukolsky and Press [25], as well as Page's calculations [16], we have obtained the evolution matrix by solving numerically the radial Teukolsky equation for spin | s | = 1 and | s | = 2 for different values of the black hole spin a . As explained by Page, the dependence of the matrix elements on the BH parameters is through the dimensionless quantities Mω and a ⋆ = a/M . Also, for massless particles, their dependence on the particle species, spin, and helicity is only through | s | and the number of polarizations, which is 2 for both photons and gravitons. Aside from that, the matrix elements also depend on the angular momentum labels ℓ and m of the modes. We have solved Teukolsky equation and obtained the evolution matrix for a range of frequencies ω ∈ [0 , ω max ] and within each frequency for ℓ ∈ [ | s | , ℓ max ] ⊂ N , and for m ∈ [ -ℓ, ℓ ] ⊂ Z . We have chosen ω max and ℓ max in such a way that the contribution of larger frequencies or ℓ to the energy, angular momentum and entanglement radiated away by the BH is negligible within the precision of our numerical calculation. Such maximum values depend on the spin parameter a ⋆ . In particular, we have make sure that the upper value of the frequency is high enough as to explore the SR to NSR transition for each of the modes that are solved.", '6. Our numerical methods': "The evolution matrix that relates modes 'up' and 'in' with modes 'out' and 'down' (corresponding to evolution from the Cauchy hypersurfaces Σ int to Σ out ) is not analytically solvable. We briefly summarize in this subsection the numerical methods we use to obtain the evolution matrix and provide references to earlier work upon which we have built our calculations.  \nThe goal is to solve the radial Teukolsky equation, which is solved numerically with purely ingoing boundary conditions at the horizon. This equation depends on eigenvalues of the spin-weigthed spheroidal harmonics, which are obtained through the small-c expansion method in section II.C of [39], taking up to the 5th order term in an expansion in ωa . In our calculations, we have checked that the first term being neglected in the expansion - the 6th order term - would lead to negligible corrections to our results for total entanglement radiated even in the worst case scenario when the black hole spin approaches a = 1. The sources of error in our calculations and their estimated contribution are provided below.  \nTo solve the radial equation, we use standard integration algorithms included in the software Mathematica through the in-built function NDSolve, with specifications AccuracyGoal= 10, PrecisionGoal= 10 and MaxStepSize= 0 . 1. This function uses an adaptative method, which might make use of different numerical integrator algorithms in each interval, based on the problem to solve.  \nIn order to avoid that numerical errors in the initial conditions excite unphysical outgoing modes at the horizon, we follow Teukolsky and Press (section VII [25]). Thus, we use a modified radial coordinate which vanishes at the horizon, and is defined by x ≡ ( r -r + ) /r + , and a modified radial function that differs from the original variable s R ω,ℓ,m in a polynomial in the variable x of order | s | + 1, where | s | is the spin of the field. This spin-dependent polynomial is computed by requiring that the 'up' (outgoing) modes approach zero asymptotically faster than the 'down' (ingoing) modes when approaching the horizon. This guarantees the numerical stability when integrating outwards for ingoing boundary conditions, which are provided near the horizon, at x = 10 -15 , following [25]. In this way we obtain the sought modes of the system, from which the greybody factors can be computed by requiring conservation of the flux(particle numbers and entanglement remain independent of the phases of the transmission and reflection coefficients, requiring us only to compute their amplitudes). This calculation was carried out in [25] as outlined in section VII, leading to formulas A.17 of that work, which relate the greybody factors and the amplitudes of the modes far from the black hole. We extract these numbers at x = 700 (i.e. r = 701 r + ). Although we have created our own code [47], publicly available resources exist to do this calculation [48].  \nWe have solved Teukolsky equation and obtained the evolution matrix for a range of frequencies ω ∈ [0 . 01 M -1 , ω ℓ max ], where ω ℓ max is an ℓ -dependent cutoff frequency, and ℓ 's within the interval ℓ ∈ [ | s | , ℓ max ] ⊂ Z .  \nWe have chosen ℓ max = 5 for photons and ℓ max = 6 for gravitons. The cutoff frequencies ω ℓ max used in our calculations are shown in Table I. The contribution above these cutoffs as well as the error in the angular eigenvalues are the main sources of error in our calculations. For these choices of ℓ max and ω ℓ max , the error introduced in the  \ndetermination of the angular eigenvalues is at the percent level at most. Concretely, for a = 0 . 99, we have checked that the correction of the first term ignored in the small-c expansion, namely the 6th order term, corrects the eigenvalues below the 0 . 25%. Furthermore, this error decreases for lower values of a and ω . Since errors are smaller for the low frequencies -which carry most of the energy radiated- the actual error in the determination of the angular eigenvalues weakly impacts integrated quantities. In particular, we have verified that the relative errors in the total power emitted in each angular multipole by frequencies above ω ℓ max are below 10 -5 % with respect to the total energy radiated by modes with that value of ℓ .  \nTo quantify the total error in the integrated spectrum we would have to take into account the contribution of the missing modes. In particular, the missing higher multipoles ℓ > ℓ max will be the main source of errors, since ω > ω ℓ max contribute well below the 0 . 17% contribution of the ℓ = 5 photon modes to total emission (the lowest of the considered contributions). This is difficult to compute exactly, but taking into account the exponential decay of the emitted power with the value of the angular multipole, both for gravitons and photons, we estimate that the corrections to total entanglement radiated are below the percent level.  \nTABLE I. These tables contain information about the emission in each multipole for a = 0 . 99. The second row in each table shows the ℓ -dependent cutoff frequency ω ℓ max , in units of M -1 , for photons and gravitons. The third row provides the power radiated in modes lying in the interval ω ∈ [0 . 01 , ω ℓ max ]. The fourth row shows the contribution to the total power radiated in each multipole ℓ -summed for all allowed m 's and for photons and gravitons- normalized to a 100% (the large relative contribution of gravitons is due to the large value of a used in this table). Note that the contribution to the total power decays exponentially with ℓ , so the contribution of higher ℓ 's is negligible. For smaller black hole spin a the main contribution to the power radiated is further shifted toward small ℓ .", 'Appendix B: Gaussian states and their evolution': 'We summarize some general properties of Gaussian bosonic systems.', '1. Gaussian states': 'Let us consider a dynamical system containing N classical degrees of freedom. Quantum mechanically, the system is described by N pairs of canonically conjugate operators ( x I , p I ) with I = 1 , ..., N or, equivalently, by N pairs of creation and annihilation operators ( a I , a † I ). The creation and annihilation operators are related to the canonical variables x and p by  \n( I  \na I a † ) = 1 √ 2 ( 1 i 1 -i )( x I p I ) and ( x I p I ) = 1 √ 2 ( 1 1 -i i )( a I a † I ) . (B1)  \nThe discussion below can be made using either sets of variables. We will use creation and annihilation operators. Let us define a vector ⃗ A by putting together all creation and annihilation operators in the following order  \n⃗ A = ( a 1 , · · · , a N , a † 1 , · · · a † N ) ⊤ , (B2)  \nwhere the transposition symbol indicates we define this as a column vector (in the literature, other orderings can be found, such as ( a 1 , a † 1 , ..., a N , a † N ) ⊤ . It is straightforward to translate expressions written below from one choice of order to another.)  \nGaussian states ρ , either pure or mixed, are states for which the quantum moments ⟨ A i 1 · · · A i n ⟩ satisfy the same relations as the statistical moments of a Gaussian multi-variable probability distribution. In particular, all moments can be determined from the first and second moments. In other words, the first and the second moments contain complete information about the state, and they can be used as to describe a quantum Gaussian state. The first moments can be encoded in a 2 N dimensional vector ⃗ µ = ⟨ ⃗ A ⟩ , while the second moments are conveniently encoded in the so-called covariance matrix, defined as σ ij = ⟨{ A i -µ i , A j -µ j ) ⟩ , where i, j = 1 , ..., 2 N and the curly bracket denotes the anti-commutator. We use the symbol ⟨· , ·⟩ to indicate expectation value, for both pure and mixed states. In the definition of the covariance matrix, the subtraction of ⃗ µ guarantees that the result is independent of the first moments. The anti-commutator selects the symmetric part of the second moments, since the anti-symmetric part, i.e., the commutator, is state-independent and completely determined from the canonical commutation relations.  \nMany properties of a Gaussian state can be easily extracted only from its covariance matrix σ . As an example, one such property is the purity P of the state, which is completely independent of the first moments. For a Gaussian state, it is obtained from its covariance matrix as P ( σ ) = 1 / √ det σ (it is one for pure states and smaller than one for mixed states). Another quantity that is completely determined only by the covariance matrix is the Logarithmic Negativity, which serves as a faithful entanglement quantifier in some cases (see below).  \nFor classical systems, if the Hamiltonian dictating the dynamics is a quadratic function of the canonical variables, the evolution is a linear canonical transformation S . Linear evolution preserves the Gaussianity of states. If ( ⃗ µ (in) , σ (in) ) describes a quantum Gaussian state, evolution under a quadratic Hamiltonian maps it to another Gaussian state, described by  \n⃗ µ (out) = S · ⃗ µ (in) , (B3)  \nσ (out) = S · σ (in) · S ⊤ , (B4)  \nwhere S agrees with the evolution matrix describing the classical dynamics.', 'a. Entropy of Gaussian states': "The von Neumann entropy of a Gaussian state ( ⃗ µ, σ ) describing a system with N degrees of freedom can be easily computed from the N symplectic eigenvalues of σ . These are N real numbers, denoted by ν I , with I = 1 , ..., N , defined as the absolute value of the eigenvalues of the matrix σ · Ω -1 where Ω -1 is the inverse of the symplectic form. [Note there are technically 2 N eigenvalues, since σ · Ω -1 is a 2 N × 2 N matrix, but the eigenvalues come in pairs ± ν I .] In terms of ν I , the von Neuman entropy of a Gaussian state reads  \nS [ σ ] = N ∑ I [ ( ν I +1 2 ) ln ( ν I +1 2 ) -( ν I -1 2 ) ln ( ν I -1 2 ) ] . (B5)  \nAs an example, for a Gaussian state with ⃗ µ = ⃗ 0 and σ = ⊕ I (1 + 2¯ n I ) ¯ I 2 -which holds for outgoing radiation of an evaporating black hole (even with thermal inputs)-we have ν I = (1 + 2¯ n I ), and we can write the entropy as S [ σ ] = ∑ N I =1 (¯ n I +1)ln(¯ n I +1) -¯ n I ln(¯ n I ), which is just the standard 'thermal entropy'.", 'b. Entanglement in Gaussian states': "Consider a partition of the system of N modes into two subsystems A and B , defined by two disjoint subsets of pairs of creation and annihilation operators ( a I , a † I ). If the total state of the system is Gaussian ( ⃗ µ AB , σ AB ), the reduced states describing each subsystem are also Gaussian. The vector of first moments and the covariance matrices describing each of the Gaussian subsystems can be obtained from ( ⃗ µ AB , σ AB ) by simply restricting to the components associated with each subsystem.  \nIf the total state ( ⃗ µ AB , σ AB ) is pure, then the von Neumann entropies of the reduce systems satisfy S [ σ A ] = S [ σ B ]. Moreover, the entropy in this case quantifies the entanglement between A and B -i.e., the von Neumann entropy is the entanglement entropy. On the other hand, if the total state is mixed, the von Neuman entropy of the subsystems can be different, and the entropy does not quantify entanglement.  \nA convenient entanglement measure, useful for pure and mixed states alike, is the Logarithmic Negativity, LN. It is defined as  \nLN( ρ ) = log 2 || ρ ⊤ B || 1 , (B6)  \nwhere ρ is the density matrix of the system, ρ ⊤ B is its partial transpose with respect to subsystem B (the same result is obtained if we partial transpose in subsystem A instead), and || · || 1 is the trace norm. The Logarithmic Negativity quantifies the violation of the PPT criterion (Positivity of Partial Transpose) for quantum states [30, 32, 49]. This criterion is satisfied by all separable states and, as a consequence, its violation signals the existence of entanglement. For a Gaussian state made of two Gaussian subsystems A and B, Eqn. (B6) can be rewritten in terms of the covariance matrix σ AB as [28]  \nLN[ σ AB ] = ∑ I Max[0 , -log 2 ˜ ν I ] , (B7)  \nwhere ˜ ν I are the symplectic eigenvalues of the partially transposed state. The eigenvalues ˜ ν I are computed as follows. Let N A and N B denote the number of modes within subsystems A and B, respectively. Hence, σ AB is a 2 × ( N A + N B ) matrix. It is not difficult to notice that the covariance matrix of ρ ⊤ B can be obtained from the covariance matrix of ρ by interchanging creation and annihilation operators in each mode within subsystem B. Furthermore, this interchange can be operationally implemented via  \n˜ σ AB = T σ AB T (B8)  \nwhere  \nT =    I N A 0 N A N B 0 N A 0 N A N B 0 N B NA 0 N B 0 N B N A I N B 0 N A 0 N A N B I N A 0 N A N B 0 N B N A I N B 0 N B N A 0 N B    . (B9)  \nIn this expression, I N is the N -dimensional identity matrix, 0 N is the N -dimensional zero matrix, and 0 N A N B is a N A × N B matrix with all components equal zero.  \nWith this, the N A + N B positive real numbers ˜ ν I featuring in the expression for LN[ σ AB ] are the absolute values of the eigenvalues of the matrix ˜ σ AB · Ω -1 .  \nWith this information, the calculation of Logarithmic Negativity generated in Hawking's process is performed as follows: Starting from the covariance matrix of the initial state σ in for each wave packet mode Q = ( j, n, ℓ, m, s ), we obtain the 'out' covariance matrix as σ out = S (tot) Q · σ in · S (tot) ⊤ Q , where S (tot) Q is the scattering matrix containing the scattering coefficients for the mode Q . From σ out , we compute the reduced covariance matrix of the set of modes we are interested in (e.g., mode out and dn), then compute the partial transpose with respect to one of the subsystems (e.g., mode out), evaluate its symplectic eigenvalues, and substitute them into (B7).  \nFor a Gaussian state, non-vanishing LN is a sufficient condition for entanglement, but it is not necessary in general. Therefore, there exist Gaussian states which are entangled but have vanishing LN. However, if either of the Gaussian subsystems is made of a single degree of freedom, (e.g. N A = 1), LN > 0 is also a necessary condition for the presence of entanglement, regardless of the size of the other subsystem. LN is also an entanglement monotone and, therefore, can be used to quantify entanglement (see [28] for further details). In the main body of this article we have used LN only for Gaussian states in which one of the subsystems has a single mode, so that it is a legitimate entanglement quantifier for our purposes. As a last remark, Logarithmic Negativity has an operational meaning for Gaussian states: it is the exact cost (measured in Bell pairs or entangled bits, ebits) that is required to prepare or simulate the quantum state under consideration [50, 51]. It is also an upper bound for distillable entanglement."}