Infrared Effects in the Late Stages of Black Hole Evaporation

As a black hole evaporates, each outgoing Hawking quantum carries away some of the black holes asymptotic charges associated with the extended Bondi-Metzner-Sachs group. These include the Poincar\'e charges of energy, linear momentum, intrinsic angular momentum, and orbital angular momentum or center-of-mass charge, as well as extensions of these quantities associated with supertranslations and super-Lorentz transformations, namely supermomentum, superspin and super center-of-mass charges (also known as soft hair). Since each emitted quantum has fluctuations that are of order unity, fluctuations in the black hole's charges grow over the course of the evaporation. We estimate the scale of these fluctuations using a simple model. The results are, in Planck units: (i) The black hole position has a uncertainty of $\sim M_i^2$ at late times, where $M_i$ is the initial mass (previously found by Page). (ii) The black hole mass $M$ has an uncertainty of order the mass $M$ itself at the epoch when $M \sim M_i^{2/3}$, well before the Planck scale is reached. Correspondingly, the time at which the evaporation ends has an uncertainty of order $\sim M_i^2$. (iii) The supermomentum and superspin charges are not independent but are determined from the Poincare charges and the super center-of-mass charges. (iv) The supertranslation that characterizes the super center-of-mass charges has fluctuations at multipole orders $l$ of order unity that that are of order unity in Planck units. At large $l$, there is a power law spectrum of fluctuations that extends up to $l \sim M_i^2/M$, beyond which the fluctuations fall off exponentially, with corresponding total rms shear tensor fluctuations $\sim M_i M^{-3/2}$.


Introduction and summary
Hawking's black hole information loss paradox is one of the most enduring mysteries in theoretical physics: how does information escape from a black hole during its evaporation? [1][2][3]. Great progress has been made on this issue in the past few years, using explicit Euclidean path integral methods. It is now possible to explicitly compute the Page curve that describes the time evolution of the entanglement entropy of the emitted Hawking radiation and the black hole, and to show that it is consistent with unitarity [4][5][6][7][8][9]. In addition, the amount of time taken for the information in a diary thrown into a black hole to return in the Hawking radiation can be reliably computed [5]. Nevertheless, some of the processes and computational prescriptions that arise in the Euclidean domain remain mysterious in the Lorentzian domain, so there is still much to be understood.
A central role in this subject is played by the semiclassical approximation, where the gravitational field is treated classically (aside from linear perturbations that can be treated as free gravitons) and the matter fields are treated quantum mechanically. This approximation excludes macroscopically large quantum fluctuations in the geometry. It is the only approximation in which we can compute the full state of the outgoing Hawking radiation. In addition, the new computational prescriptions for computing the entanglement entropy of the exact state of the Hawking radiation [4][5][6][7][8][9] are expressed in terms of a single semiclassical geometry, as are other similar powerful theoretical tools and results (the Ryu-Taganacki formula [10,11], the quantum focussing conjecture [12] and the covariant entropy bound [13]).
On the other hand, it has been known since the work of Page in the 1980s [14] that the semiclassical approximation actually fails drastically during the course of black hole evaporation. This failure arises as follows: each emitted quantum carries of a momentum ∼ M −1 in a random direction, where M is the mass of the black hole in Planck units, and the corresponding change in the velocity of the black hole is of order ∼ M −2 . This change in velocity causes a net displacement in the center-of-mass of the black hole of order ∼ M after an evaporation time ∼ M 3 . During the evaporation process we have n ∼ M 2 such kicks that accumulate as a random walk, giving a total net uncertainty in the black hole location of order ∼ √ nM ∼ M 2 , much larger than the size of the black hole. Thus we have superpositions of macroscopically distinct geometries.
Several authors have argued for the importance of the center-of-mass fluctuations in understanding the unitarity of black hole evaporation [15][16][17]. They note that unitarity is required only for the evolution in the complete Hilbert space, not in the subspace that corresponds to a single semiclassical geometry, and that the relative phases of different semiclassical geometries in a quantum superposition contain information. However, there are counterarguments [18,19] which suggest that the breakdown of the semiclassical approximation is fairly innocuous. First, there are situations where black holes evaporate in anti-de Sitter space where the center of center-of-mass spreading is suppressed but where there is still an information loss paradox [19]. Second, the dimension of the Hilbert space associated with the center-of-mass motion is negligibly small compared to the relevant scale of the exponential of the Bekenstein-Hawking entropy, since it scales as a power law in the mass M of the black hole 1 .
In Ref. [20] we show that the center-of-mass fluctuations give rise to large corrections to the angular distribution of the Hawking radiation. We also argue there that those corrections remove one of the primary objections to the proposal that soft hair on black holes plays a key role in how unitarity of the evaporation is achieved [21][22][23][24][25], by increasing the number of soft hair modes that can interact with outgoing Hawking quanta.
The purpose of this paper is to study the macroscopic fluctuations in the geometry of evaporating black holes, in more detail than hitherto. A more detailed understanding may be useful for eventually extending some of the theoretical tools discussed above to situations where infrared quantum fluctuations are large. It may also shed light on the role of soft hair. Finally, some of the results derived here were used in the computations of Ref. [20].
The theoretical framework we use to study these fluctuations is as follows. In the classical theory, the geometry of a stationary black hole is determined by the conserved charges on future null infinity, including the soft hair charges associated with extensions of the Bondi-Metzner-Sachs (BMS) group [21,22]. We assume that the this property remains true in the quantum theory, and evolve the charges using a simple model, described in Secs. 2.1 and 3.6 below. We extend similar previous studies [14,15] of the fluctuations in a number of directions: • We extend the computations to late times when the black hole mass M is small compared to its initial mass M i , by making use of the approximation that the fluctuations in charges are small compared to their expected values [Eqs. (2.5) below]. This approximation is valid until M ∼ √ M i . The variance in the center of mass location does not evolve significantly as the black hole shrinks from M ∼ M i to M M i , but remains ∼ M 2 i .
• We compute the evolution of the fluctuations in the mass of the black hole. This is of order unity in Planck units after an evaporation times, but grows at late times according to ∆M ∼ M 2 i /M 2 [Eq. (2.5c) below]. It follows that ∆M ∼ M when M ∼ M 2/3 i . At this epoch, there is an order unity amplitude for the evaporation to be completed (M = 0), but there is also an order unity amplitude for the black hole mass to be macroscopic, M ∼ M 2/3 i .
• We extend previous studies to include the charges associated with an extension of the BMS group [26][27][28][29]. These charges are reviewed in Sec. 3.1 below, and are higher-l analogs of the center-of-mass, momentum and spin that are encoded in the asymptotic metric near future null infinity. We show that for evaporating black holes only some of these charges are independent. The independent charges can be taken to be the so called super center-of-mass charges, or soft hair. These charges can be parameterized in terms of the supertranslation required to set them to zero, a function Φ on the two-sphere with dimensions of length with only l ≥ 2 components. See Secs. 3.1, 3.2 and Appendix B for more details.
The accumulated fluctuations in the supertranslation Φ are relatively small. For multipole orders l of order unity, the fluctuations are of order unity in Planck units [Sec. 3.8]. There is also a contribution to the fluctuations associated with quanta that have only partially arrived at future null infinity by the time the charges are measured, which we compute in Sec. 4. This contribution gives a power law spectrum of fluctuations extending up to l ∼ M 2 i /M which is subdominant at early times, but becomes dominant when M becomes small compared to M The organization of this paper is as follows. The predictions of our model for the evolution of the Poincaré charges are given in Sec. 2. Section 3 extends this model to include soft hair charges. In Sec. 4 we consider some additional contributions to the fluctuations of the soft hair charges that are associated with quanta that have only partially arrived at future null infinity by the time the charges are measured. Preliminary versions of some of the results were presented at the conferences [30,31] [14] and explored in more detail by Nomura, Varela and Weinberg [15]. Page estimated the fluctuations in position of the black hole after an evaporation time. Below we extend his computations to late times and also estimate the mass fluctuations of the black hole. A more sophisticated model which includes all the BMS and extended BMS charges will be given in Sec. 3 below. The results from that more sophisticated model for the Poincaré charges agree qualitatively with those of the simple Page model discussed here. We start by describing the motivation for the model. Consider the evaporation of a black hole of mass M 1, in Planck units with 8πG = = c = 1. Roughly one Hawking quanta per time ∆t ∼ M is emitted; this is explicit in an orthonormal wavepacket mode basis. Each quantum carries an energy ∆E ∼ M −1 and spatial momentum ∆p ∼ M −1 , which change the energy and momentum of the black hole by corresponding amounts. The fractional fluctuations in ∆E and ∆p are of order unity, since they are carried by a single quantum. Also the spatial momentum can be in a random direction.
The classical stochastic process is defined by 1c) x n+1 = x n + p n . (2.1d) Here n = 1, 2, 3 . . . labels the steps, one for each emitted quantum. The variable M n is the mass of the black hole at step n. The quantity n is a random variable which takes on the values 0 and 1 with probability each of 1/2. The mass evolution equation (2.1a) describes the mass of the black hole being reduced, with a probability to occur of order unity, by each emitted quantum. The emission process takes a time of order the current mass of the black hole, as encoded in the time evolution equation (2.1b); t n is the time of step n. The spatial momentum of the black hole after step n is p n , and its evolution is governed by Eq. (2.1c). The magnitude of the emitted spatial momentum is n /M n , the same as the energy, since the Hawking quantum is massless. However the momentum can be in either direction (we are using a one dimensional model of the black hole motion). The directionality is encoded in the random variable δ n , which takes on the values −1 and 1 with equal probability 2 . The change in the position x n of the black hole during step n is the momentum p n divided by the mass M n times the time interval M n , which yields the evolution equation (2.1d). All of the variables n , δ n for n = 1, 2, 3 . . . are uncorrelated. The initial conditions at n = 1 are taken to be M 1 = M i , the initial black hole mass, and p 1 = x 1 = t 1 = 0. The evaporation model (2.1) clearly incorporates a number of simplifications and approximations. However we believe that the key predictions of the model are robust and are insensitive to these simplifications. Some of the approximations are: • It incorporates only one spatial dimension, and treats only the Poincaré conserved charges of the black hole, neglecting the additional charges associated with the BMS algebra and its extensions. These restrictions will be lifted to some extent in Sec. 3 below.
• It treats all the fluctuations classically rather than quantum mechanically. This restriction will be lifted to some extent in Sec. 2.5 below.
• It models a continuous process as a discrete process. However we believe that this idealization does not affect the scale of the late time fluctuations predicted by the model.
• It neglects any initial fluctuations in the black hole's conserved charges. This is acceptable since at late times the fluctuations will be dominated by the cumulative effects of the emitted Hawking quanta, for any reasonable estimate of initial fluctuations 3 .
• Clearly, the the model could be generalized and made more precise by inserting dimensionless parameters of order unity into each of the equations (2.1); this would not change the qualitative predictions.
too literally. In particular, if we quantize a scalar field near future null infinity on a spherical harmonic basis, then the operator (3.12) below that describes linear momentum radiated consists entirely of cross terms between different modes, rather than individual modes carrying linear momentum. However, consider the following slight generalization of the model: at each timestep there are two independent random variables n and n that take on the values 0 and 1 with equal probability, with n = 1 representing the emission of a l = 0 scalar quantum and n = 1 representing an l = 1 quantum. Then Eqs. (2.1) effectively hold with the n in Eq. (2.1a) replaced by n + n and with the n in Eq. (2.1c) replaced by n n . The qualitative predictions of the model are unchanged by this refinement. 3 For example, for a particle of mass M i , the standard quantum limit on the uncertainty in position after a time t is ∆x t/M i [32]. This uncertainty is of order ∼ M i after an evaporation time t ∼ M 3 i , much smaller than the late time uncertainty ∆x ∼ M 2 i due to the Hawking quanta, cf. Eq. (2.5a).
• The motion of the black hole is treated non-relativistically. This is consistent since the motion is still non-relativistic at the time the model breaks down, cf. Eq. (2.5b) below.
The key feature of the model is that independent, uncorrelated fluctuations are introduced into the black holes conserved charges at each timestep or each emission event. Those fluctuations ultimately originate in incoming modes of quantum fields at past null infinity, which are orthonormal and so uncorrelated for the incoming vacuum state.

Early time predictions: large position fluctuations
At times small compared to an evaporation time, t M 3 i , simple random walk arguments can be used to estimate the fluctuations in the black hole charges, as first done by Page [14] and explored in more detail by Nomura, Varela and Weinberg [15]. In this section we review these early time predictions of the model.
We define τ to be the time since black hole formation in units of the evaporation time, τ = t/M 3 i . The fluctuations in the various quantities are (1 + δ nm ), gives p 2 n = (n − 1)/(2M 2 i ). Combining this with t n = (n − 1)M i from Eq. (2.1b) yields the momentum fluctuation estimate (2.2b). Similarly the mass evolution equation (2.1a) in this approximation yields M n = M i − n j=1 j /M i , and an analogous argument yields the mass fluctuation estimate (2.2c). Finally solving the position evolution equation (2.1d) yields x n+1 = n j=1 j k=1 k δ k /M i . Squaring and taking the expectation value gives x 2 n = n 3 /(6M 2 i ), up to fractional corrections ∼ 1/n, which yields the estimate (2.2a). 4 discussed in the introduction, this implies that there is an order unity amplitude for the evaporation to be completed at the same time as there is an order unity amplitude for the black hole mass to be macroscopic, M ∼ M 2/3 i . There is an elementary argument for the large, late time mass fluctuations, which is as follows [33]. Consider a black hole with initial mass M i . After an evaporation time ∼ M 3 i , when the mass is M i /2, the spread in mass is of order unity, from the random walk estimates of Sec. 2.2. Consider now two different histories, one with mass M i /2 at this time and one with mass M i /2 + 1. For the subsequent evolution of these two histories, we consider just the evolution of the mean mass for simplicity, neglecting fluctuations produced by the subsequent Hawking emission. The corresponding mean masses a time t later are [(M i /2) 3 − 3t/2] 1/3 and [(M i /2 + 1) 3 − 3t/2] 1/3 , and when the mass of the first history is zero, the mass of the second history is ∼ M

Numerical simulation of model
It is also straightforward to numerically simulate the stochastic model

Beyond the classical stochastic model
So far, the motion of the black hole has been treated classically. However, generalizing to a more detailed quantum mechanical treatment does not qualitatively change the results, as we now outline. The black hole center-of-mass motion can be described by its Wigner function W(p, x), with variance-covariance matrix for any function f . As before, we idealize the evolution as a series of n steps, each of which has two parts. First, the black hole emits a quantum, under which by momentum conservation Σ transforms as where After n ∼ M 2 steps, the effect of the initial value of Σ is negligible, and the predicted scalings of position uncertainty and momentum uncertainty agree with Eqs. (2.2) above.

Evolution and fluctuations of black hole extended Bondi-Metzner-Sachs charges
In this section we generalize the Newtonian model of the previous section to include all of the Bondi-Metzner-Sachs (BMS) charges of the black hole and in addition the charges associated with the extended BMS algebra [34][35][36][37][38].
As discussed in the introduction, in this paper we focus on charges of the black hole measured at future null infinity. One could instead consider the symmetries and charges defined on the black hole horizon, which has a different symmetry algebra [21,22,[39][40][41][42][43][44][45][46][47]. These charges are related to charges at future null infinity by global conservation laws, as detailed in Ref. [21], assuming that one can find the appropriate identification between horizon symmetry generators and asymptotic symmetry generators (currently known in some special cases). However our approach here follows the perspective of a distant asymptotic observer.

Review of BMS and extended BMS charges
We start by reviewing the nature of the BMS and extended BMS charges of asymptotically flat spacetimes. For more details on this topic see the review by Strominger [26] and the expositions [29,48]. The BMS group is the group of asymptotic Killing vectors that act on spacetimes which are asymptotically flat at future null infinity [49][50][51]. Associated with each asymptotic Killing vector ξ or generator of the group, and with each cut of future null infinity, there is a conserved charge Q( ξ) [52,53].

2)
A, B = 1, 2, and all the functions that appear in the metric are are functions of u and θ A . The metric h AB is the unit round metric on the two-sphere and is used to raise and lower capital Roman indices, and D A is the associated covariant derivative. There are three important, leading-order functions in the metric's expansion coefficients [34,37,49,[54][55][56][57][58]: the Bondi mass aspect m(u, θ A ), the angular-momentum aspect N A (u, θ A ), and the symmetric tensor C AB (u, θ A ) whose derivative is the Bondi news tensor. Evolution equations for these metric functions in terms of retarded time are given by (FN,2.4), (FN,2.11a) and (FN,2.11b). The leading order components of the stress energy tensor are given by (FN,2.6) and involve functionŝ The BMS algebra is the algebra of infinitesimal diffeomorphisms on null infinity that map from one Bondi frame (u, θ A ) to another. A general BMS generator can be written as [Eq. (FN,2.13)] where Y A is a globally smooth conformal Killing vector on the 2-sphere (the set of which is isomorphic to the Lorentz algebra), and α is an arbitrary smooth function that parameterizes the supertranslation transformations. Two different extensions of the BMS algebra have been proposed. Barnich and Troessaert [34,54] suggested an extension that includes all local conformal Killing fields Y A on a 2-sphere, allowing isolated singular points. This replaces the Lorentz algebra of vector fields with an infinite dimensional Virasoro algebra; see Ref. [26] for more details. Campiglia and Laddha [27,28] suggested extending the Lorentz transformations to include all smooth infinitesimal diffeomorphisms Y A on a 2-sphere.
We focus here on the second extension, whose status can be summarized as follows. It arises as the symmetry group of an extended phase space of general relativity at future null infinity, in which fewer 6 of the diffeomorphism degrees of freedom are fixed than is usual [27,28,59]. While the presymplectic current of general relativity diverges at null infinity in the extended phase space, by exploiting a redefinition freedom [53] in the presymplectic current it can be made finite. Compère, Fiorucci and Ruzziconi use this method to construct a finite presymplectic current and derive charges Q( ξ) associated with each symmetry generator ξ in the extended algebra [29]. However, their construction uses a specific coordinate system, and so is not obviously local and covariant (which would be necessary for uniqueness). Indeed it can be shown there is no redefinition of the presymplectic current that is local and covariant throughout the spacetime and which makes the presymplectic current finite at null infinity [59]. Nevertheless, the charges defined by Ref. [29] can be shown indirectly to be covariant and unique [60]. They have also been shown to be consistent with the leading and subleading soft graviton theorems [29].
The charges associated with a symmetry of the form (3.4), with Y A an arbitrary smooth vector field on the two sphere, on a cut u = constant of a stationary region of future null infinity can be written as [29] In a given Bondi frame 8 , the charges we consider are (Sec. III of FN): • The Bondi four momentum P α which is encoded in l = 0, 1 pieces of the the Bondi mass aspect m(u, θ A ) and conjugate to normal translations.
• The supermomentum charges which are encoded in the l ≥ 2 pieces of m(u, θ A ) and are conjugate to supertranslations. They encode a separate energy conservation law at each angle [55].
• The angular momentum J αβ which is encoded in the l = 1 piece ofN A (u, θ A ) and is conjugate to the Lorentz generators Y A (conformal Killing vectors). As usual this can be split into intrinsic angular momentum, and orbital angular momentum or center-of-mass charge (center-of-mass minus velocity times time).
• The superspin charges which are encoded in the magnetic parity piece of the l ≥ 2 piece ofN A (u, θ A ), and are conjugate to l ≥ 2 magnetic parity symmetry generators Y A in the extended algebra. They encode a separate conservation law for intrinsic angular momentum at each angle [58].
• The super center-of-mass charges which are encoded in the electric parity piece of the l ≥ 2 piece ofN A , and are conjugate to l ≥ 2 electric parity symmetry generators Y A in the extended algebra. They encode a separate conservation law for orbital angular momentum or center-of-mass charge at each angle [38,61]. In the context of black holes they are also called soft hair [21,22,37].
We can parameterize these charges in terms of a set of symmetric, tracefree tensors J ij , J ijk , . . . as follows. For any symmetric, tracefree Cartesian tensor Y i 1 ...i l , we consider the symmetry generator (3.4) with α = 0 and with where L is the multi-index i 1 . . . i l , n L = n i 1 . . . n i l , and n i is the unit vector (sin θ cos ϕ, sin θ sin ϕ, cos θ). We define the symmetric tracefree tensor J L by demanding that the corresponding charge (3.5) is Y L J L . The tensor J L is related to the lth multipole electric parity piece ofN A bŷ with g l = 2(2l + 1)!!/(l(l + 1)l!), from Eq. (3.5) and using the identity (C5) of Ref. [48].

Charges in stationary regions of future null infinity
We will idealize the evaporation of a black hole a sequence of transitions between stationary states: after each Hawking quantum is emitted, the black hole settles down to a stationary state, then the next quantum is emitted, and so on.
In regions of future null infinity that are stationary, the various charges discussed above are not all independent. This follows from the fact that there exists a canonical Bondi frame associated with the stationary region in which the metric functions take a simple form that encode the mass and intrinsic spin (see, for example, Sec. II.D of FN): Now a general Bondi frame will be related to the canonical frame (3.8) by a nonlinear BMS transformation of the form (FN,2.12), parameterized by a Lorentz transformation and a supertranslation, and therefore the metric functions m, C AB and N A in the general frame encode just one infinite family of charges, and not three (see Appendix B). We will focus here on the independent charges, which we take to be the momentum P α , angular momentum J αβ , and super center-of-mass charges; the other charges can be determined from these. We also show in Appendix B that the super center-of-mass charges are determined to a good approximation by the shear tensor C AB , so we focus on this quantity in subsequent sections rather than on N A . Specifically, we decompose C AB in terms of an electric parity potential Φ and a magnetic parity potential Ψ via where we take the l = 0, 1 pieces of Φ and Ψ to vanish. We also expand the electric parity potential Φ in terms of a set of symmetric, tracefree tensors The relation between the tensors Q L and the super center-of-mass charges J L is given in Eqs. (B.9) and (B.10) of Appendix B.

Stationary to stationary transitions and changes in the charges
In our model of black hole evaporation, each emission of a Hawking quantum is idealized as a stationary to stationary transition as viewed at future null infinity. Specifically, this means that the spacetime at some early retarded time u 1 is vacuum near future null infinity, and is also approximately stationary there. There is subsequently a burst of gravitational waves and/or matter energy flux to infinity, and the spacetime is again vacuum and approximately stationary near future null infinity at some later retarded time u 2 . In this section we will give formulae for the changes in the BMS and extended BMS charges in such transitions, in terms of fluxes to null infinity of mass-energy or gravitational-wave energy. These formulae will be one foundation of our model of black hole evaporation of Sec. 3.4 below, and are derived in Appendix D.
The changes in the linear momentum P α and angular momentum J αβ have the same form as in special relativity, but with the stress-energy fluxes supplemented by gravitational wave terms. The total energy radiated per unit solid angle in either matter or gravitational waves is Similarly we define whereT uA = lim r→∞ r 2 T uA and T uA is a kind of gravitational wave angular momentum flux given in terms of C AB and N AB by Eq. (FN,3.23). The quantity ∆E A can be interpreted as angular momentum radiated per unit solid angle in either matter or gravitational waves. We also define the quantity 14) The components of angular momentum are given in terms ofN A by Finally we turn to the super center-of-mass charges. As discussed in Sec. 3.2 above, the center-of-mass charges are encoded in the electric parity potential Φ for the shear tensor C AB defined in Eq. (3.9). The change ∆Φ in this potential (which encodes the gravitational wave memory) is given by 9 Here m 0 is the rest mass of the initial Bondi 4-momentum, ∆P is the momentum change (3.12), P is the projection operator that sets to zero the l = 0, 1 pieces of functions on the sphere, and D is the angular differential operator The formula (3.17) is valid in initial rest frames, i.e., Bondi frames in which the spatial components of the initial Bondi 4-momentum vanish.

Evolution model
We now describe the evolution model for the black hole evaporation process, which is based on the same philosophy as the simple Newtonian model of Sec. 2 above. As before the evaporation is treated as a series of discrete steps, and each step is a classical stochastic event. The only generalization is that all of the BMS charges are included instead of just the Poincaré charges.
After the black hole is first formed, with initial mass M i , it rapidly settles down to a stationary state, on a timescale ∼ M i . We will call the canonical Bondi frame associated with this initial stationary state the initial Bondi frame. After each Hawking quantum is emitted, we assume that the black hole settles down again to a new stationary state, with a new associated Bondi frame which we call the instantaneous Bondi frame, before the next quantum is emitted. The changes in the BMS charges, that is, the charges carried off by the Hawking quantum, will have a simple universal form in the instantaneous Bondi frame, but their form in the initial Bondi frame will be more complicated and will depend on the values of all the BMS charges.
We now describe the model in more detail. We denote by x α = (u, r, θ A ) the initial Bondi frame. At the nth step, the black hole BMS charges in this frame are 4-momentum P α n , angular momentum J αβ n , and the tensor C n AB (θ A ) which encodes the super center-of-mass charges. Our goal is to derive a formula for the charges at step n + 1 in the initial Bondi frame, in terms of the corresponding values at the nth step, and also the changes in the charges in the instantaneous Bondi frame.
We compute the BMS transformation from the initial frame to the nth instantaneous frame in two stages. First, we make a supertranslation parameterized by a function β on the 2-sphere (see Appendix B of FN), to a Bondi frame xα = (û,r, θÂ). We can write this function as β = t 0 −t i n i +β 2 , where n i = (sin θ cos ϕ, sin θ sin ϕ, cos θ), t µ is a 4-vector associated with the normal translation piece of the transformation, and β 2 is purely l ≥ 2. The charges in the hatted Bondi frame are [Eqs. (FN,B7) and (FN,B8)] Here δJ αβ [β 2 ] is given by Eqs. (FN,B7) with β replaced by β 2 . Next, we choose the supertranslation β 2 to make C n AB = 0, (3.20) which determines β 2 uniquely as shown in Sec. II.D of FN. We also choose the translation to make P nα Jαβ n = 0, (3.21) which makes the hatted frame be a center-of-mass frame. A translation which achieves this is where M 2 n = − P 2 n . Next, we perform a boost Λᾱα from the hatted Bondi frame to the instantaneous Bondi frame xᾱ = (ū,r, θĀ). From Eqs. (FN,B3) and (FN,B6) the charges transform as

23a)
Jᾱβ n = ΛᾱαΛββJαβ n , (3.23b) Here ϕ : S 2 → S 2 is the conformal isometry of the two sphere associated with the boost, as described after Eq. (FN,B3), and ω ϕ is given by Eq. (B.2). The boost is determined in the usual way by the requirement that the new frame be a rest frame, ie P¯i n = 0. Jᾱβ n+1 = Jᾱβ n + ∆Jᾱβ n , (3.25b) The prescription we use for the changes ∆Pᾱ n , ∆Jᾱβ n and ∆C n AB is discussed in Sec. 3.6 below.
We now transform the new charges back to initial Bondi frame, and express the results in terms of the initial Bondi frame components of the old charges. The final result is 10 P α n+1 = P α n + Λ ᾱ α ∆Pᾱ n , (3.26a) Here the right hand sides are functions of the changes in the charges in the instantaneous Bondi frame, and of the charges at step n: the Lorentz transformation Λ ᾱ α and associated map ϕ are determined as a function of P α n by Eq. (3.24), while the quantity δJ αβ is given as a function of C n AB by Eqs. (3.19c), (3.20) and Eq. (FN,B7) with β replaced by β 2 .
From the structure of these evolution equations we see that the evolution (3.26a) of the 4-momentum is uncoupled from that of the angular momentum and super centerof-mass. In particular, this implies that the results of the simple model of Sec. 2 above for the 4-momentum evolution should still be valid. We also see that the super center-of-mass evolution (3.26c) is uncoupled from the angular momentum, and can be computed once the 4-momentum evolution is known. Finally, the angular momentum evolution (3.26b) depends on both the 4-momentum evolution and the super center-ofmass evolution.

Slow motion approximation
The velocity of the black hole is of order ∼ 1/M at late times, up to a logarithmic factor, from Eq. (2.5b) above. This is small compared to unity until the Planck scale M ∼ 1; in particular it is small compared to unity when M ∼ M 2/3 i , the epoch when ∆M ∼ M . Therefore in the evolution equations (3.26) it is a good approximation to treat the velocity as small. Expanding to linear order in the velocity v n and splitting into space and time components and writing P n = (E n , P n ) = M n (1, v n ) yields Here we have parameterized the electric quadrupole piece of C n AB in terms of a tensor Q n ij using the definitions (3.9) and (3.10). Equations (3.27) show explicitly the leading coupling of the super center-of-mass to the angular momentum evolution, which occurs through the quadrupole Q n ij . The velocity terms in these equations are suppressed relative to the other terms by a factor ∼ M , so we can take the zero velocity limit. This yields M n+1 = M n + ∆E n , (3.28a) The angular momentum evolution is more transparent if we switch to a different set of variables, namely the center of mass X i n = (P i n u n − J 0i n )/M n and the intrinsic angular momentum S ij n . Here u n is the value of the retarded time coordinate u along future null infinity I + at step n. We also note from Eq. (3.16b) that for u large compared to ∆u n = u n+1 − u n ∼ M , the change in the space-time components of the angular momentum can be written as where ∆J 0i n is given by the first term on the right hand side of Eq. (3.16b). Using these new variables to rewrite Eqs. (3.28c) and (3.28d) and making use of Eqs. (3.27a) and (3.27b), we obtain our final result for the evolution prescription in the slow motion approximation:

Changes in the charges in the instantaneous Bondi frame
To complete the model we need to give a prescription for the changes ∆E n , ∆P n , ∆J ij n , ∆J 0i n and ∆C n AB to the charges in the nth instantaneous Bondi frame. For the first four of these charges, we can use simple order of magnitude estimates as we did in the Newtonian model of Sec. 2 above. However, for the super center-of-mass charges, the relevant physics is less familiar, which is why we derived general formulae for the changes in the charges in terms of fluxes to infinity in Sec. 3.3 above. We now use those formulae as a guide to develop a prescription.
For a single Hawking quantum, Eqs.
lines of the Newtonian model of Sec. 2, is given by Here as before n = 0 or 1 is a random variable, and Ξ n , χ n and Θ n are independently and randomly distributed on the unit sphere.
For the super-center-of-mass charges, we parameterize ∆C n AB in terms of a potential ∆Φ n as in Eq. (3.9), which in turn is given by Eq. (3.17). We can neglect the second term on the right hand side, since it is smaller than the first term by a factor ∼ M 2 . The remaining term is the l ≥ 2 piece of the time integral of the energy flux, which by Eqs. (3.11) and (3.31) scales as ∼ M −1 . We therefore adopt the simple model 13 where ϕ n is the random process on the twosphere given by with ϕ lmn = 0 and ϕ lmn ϕ * l m n = c 2 l δ ll δ mm δ nn . Here c l are dimensionless constants which are of order unity for l of order unity. We will take the coefficients c l to fall off exponentially with l as l → ∞, since this is the behavior of the transmission coefficients that enter into the amplitudes for emitted Hawking quanta 14 . Here the Kronecker delta δ nn enforces the fact that successive emitted quanta are uncorrelated, while the other two Kronecker delta factors ensure isotropy.

Results for Poincaré charges
The evolution prescription for the Poincaré charges given by Eqs. (3.30a) -(3.30d) is similar to the simple Newtonian model of Sec. 2 above; here, however, it has been derived from the full BMS kinematics. There are also some differences from the model of Sec. 2, aside from the trivial generalization to three dimensions. First, there is the evolution equation (3.30d) for the intrinsic angular momentum, which was not tracked in Sec. 2. We can make order of magnitude estimates of the terms in Eq. (3.30d), using X i n ∼ M 2 , P i n ∼ 1, ∆P i n ∼ ∆E n ∼ M −1 , and ∆J ij n ∼ ∆J 0i n ∼ 1. The first two terms on the right hand side are of order ∼ 1 while the third term is ∼ M −1 and can be neglected. The random walk estimate then gives that the fluctuations δS ij in S ij at late times are of order ∼ √ n ∼ M i , where n ∼ M 2 i is the number of steps. The dimensionless angular momentum parameter is then of order 15 (3.37) Hence the spin fluctuations are unimportant until M ∼ √ M i , which occurs after the regime M ∼ M 2/3 i of interest in this paper. A similar conclusion was reached in Appendix C of Ref. [15].
A second difference from the model of Sec. 2 is the angular momentum term (third term) on the right hand side of Eq. (3.30c) for the position evolution, which does not appear in the corresponding Eq. (2.1d). This term is of order ∆J 0i n /M n ∼ M −1 , and is statistically independent of the second term which is ∼ 1 at late times, so it gives a subdominant contribution to the displacement fluctuations.
To summarize, we have argued that our model given by Eqs.

Results for super center-of-mass charges
Turn now to the super center-of-mass charges. Combining the definition (3.9) of the potential Φ n , the result (3.30e) for how the charges are updated at each step, and the model (3.34) for the charges carried away by the nth quantum, we obtain for the potential Φ n at late times Φ n =  The calculations so far in this paper have neglected some transient effects that are important for the super center-of-mass fluctuations at high multipole orders. Specifically, consider the black hole charges evaluated on the cut S of I + given by u = u 0 in the initial Bondi frame. We have assumed that each emission event impacts I + either completely before S, or completely after S, since we have counted only complete emission events and not partial events. In fact, there will be approximately ∼ M 2 i /M emission events or outgoing quanta for which the associated stress energy impacts I + partially before S, and partially after S, whose contributions to the charges have not been correctly accounted for. In this section we will estimate the contribution from these events.
We start by summarizing the results. We denote by Φ tr the "transient" contribution from the partially counted quanta to the potential Φ for the shear tensor at u = u 0 . We parameterize the spectrum of angular fluctuations of Φ tr by a quantity dΦ 2 tr /d ln l defined so that 17 where the angular brackets denote expected value. Here l denotes multipole order, which we treat as a continuous variable at large l. Our result for the fluctuations of Φ tr at large l is dΦ 2 where the symbol ∼ means that we have dropped constant factors of order unity. Here M is the expected mass of the black hole and σ 2 ∆ the variance in the center-of-mass location at retarded time u 0 . Using the late-time result (2.5a) for this variance 18 we can rewrite the power spectrum as Thus the transient contribution to the spectrum is a power law up to a critical angular scale given by l crit ∼ M 2 i /M , and at higher l it is exponentially suppressed. The transient contribution (4.3) is small compared to the previously computed contribution (3.39) for l of order unity, but dominates for 1 l l crit where the previous contribution is exponentially suppressed.
The spectrum (4.3) characterizes the fluctuations in the potential Φ, or equivalently of the tensors Q L , at large l. The super center-of-mass charges J L are related to these tensors by a factor of the black hole mass, as discussed in Sec. 3.8 above.
Finally, we note that the shear tensor (3.9) is related to Φ by two angular derivatives. Hence the spectrum of fluctuations of C AB is given by the right hand side of Eq. (4.3) multiplied by l 4 , which scales ∝ l −5 . The total rms fluctuation in C AB is of order Here the first term is the contribution (3.39) previously computed, and the second term is the transient contribution (4.3). The first term dominates at early times, while the second term begins to dominate when M becomes small compared to M 2/3 i .

Derivation
We now turn to the derivation of the spectrum (4.2). Our derivation is based on the same kind of heuristic model as used in earlier sections of the paper. A more rigorous derivation based on the two point function of the flux operator yields qualitatively the same result and will be given elsewhere. As previously discussed, each outgoing quantum is characterized by an outgoing flux ∼ M −2 over a timescale ∼ M . For simplicity we will assume that the dependence on time is identical for each outgoing quantum. Thus for the nth quantum we assume, consistently with Eqs. (3.33) and (3.34), where u n+1 − u n = M n [cf. Eq. (2.1b)] and F is a fixed smooth nonnegative function with F(x) = 0 for |x| > 1 and dxF(x) = 1. Also ϕ n is the random process given by Eq. (3.35) but now with the l = 0, 1 terms included. Next, we sum over all the quanta and transform from the instantaneous Bondi frame to the initial Bondi frame, neglecting the relative boost in accordance the slow motion approximation of Sec. 3.5. This givesT where ∆ = ∆(u) is the center of mass location of the black hole. Finally we evaluate the net change in the potential Φ for the shear tensor from u = −∞ to u = u 0 by applying Eq. (3.17) in the initial Bondi frame, neglecting the subdominant second term, and using Eq. (3.11). This gives Here we have defined the functionF(x) = x −∞ dx F(x ), which satisfieŝ There are two types of terms that arise in the sum (4.7). For sufficiently early quanta for which the argument ofF is larger than 1, we can drop theF factor by Eqs. (4.8), and the computation reduces to that of Sec. 3.8. For later quanta for which the argument ofF is less than one in absolute value, the computation is modified. These are the terms with |u 0 − u n | ∆ ∼ M 2 , that is, the final ∼ M quanta in the sum. These terms will enhance the fluctuations at high multipole orders l.
We now make an order of magnitude estimate the spectrum of fluctuations as a function of angular scale of the expression (4.7). We specialize to high multipole orders l, for which the sum over quanta will be dominated by the late quanta just discussed. For these terms we can drop the factor ϕ n (θ A ), since its dependence on θ A is exponentially small at large l, and we have ϕ n ∼ 1 at low l. We also approximate ∆(u) by its final value ∆(u 0 ), and M n by its final value at u = u 0 which we denote simply by M . We assume for simplicity that there is a valuen of n for which un = u 0 , and define j = n −n, so that to a good approximation we have u n = u 0 + jM . We assume initially that ∆ is fixed with ∆ = |∆| M ; later we will consider the effect of fluctuations in ∆. We will also approximate ∆/M by the integer l * that is closest to it, and define µ to be the cosine of the angle between n and ∆. With these definitions and approximations we have that the relevant terms in Eq. (4.7) are where the subscript "tr" denotes the transient contribution to Φ from the late incomplete quanta.
(A.2a); we will restore this fractional error estimate at the end of the computation, when we have computed δM 2 n . Next from the expected value of the time evolution equation (2.1b) and the definitions (A.1) we findt n+1 = n k=1M k . Converting this to an integral 20 and using the expression (A.3) gives We now turn to computing the fluctuations. From Eq. (A.2b) we obtain The product inside the square brackets can be evaluated by taking the logarithm, converting the sum to an integral, and using Eq. (A.3), which yields where q n+1 = n k=1 δ j . Squaring and taking the expected value gives δM 2 n+1 = n/(4M 2 n ). Using this expression to evaluate the error estimate in Eq. (A.6) and eliminating n in favor of M =M n using (A.3) finally yields Note that this is the fluctuation in mass at fixed n, to be distinguished from the more physically relevant fluctuations in mass at fixed time t [cf. Eq. (2.4c) above], which we compute below. We next compute the fluctuations δt n . From the time evolution equation (2.1b) we obtain δt n+1 = n k=1 δM k , and squaring and taking the expected value using Eq. (A.6) gives Converting the sums to integrals as before yields Now in this simple discrete model of the black hole evolution, the black hole mass M (t) at a given time t is obtained by evaluating M n at the value of n for which t n is closest to t. Hence the fluctuations in mass at fixed time t are given by Squaring and taking the expected value gives and evaluating the integral gives the expression (2.4a).

B Independent charges in stationary regions of future null infinity
In this appendix we derive the relationships between the various charges of the extended BMS algebra that apply in stationary regions of future null infinity, and show that the independent charges can be taken to be the 4-momentum P α , the angular momentum J αβ , and the super center-of-mass charges. Equivalently, we show that the supermomentum and superspin charges are determined in terms of the other charges, and so can be neglected for our purposes.
In the canonical Bondi frame associated with the stationary region, the metric functions take the simple form (3.8). We now make a nonlinear BMS transformation to a general Bondi frame, of the form (FN,2.12), following Appendix B of FN. Quantities in the new frame will be denoted with overbars. The transformation is parameterized in terms of a conformal isometry ϕ : S 2 → S 2 of the 2-sphere into itself, and a function β on the two sphere [denoted by α in Eq. (FN,2.12)]. The Bondi mass aspect in this general frame is given from Eqs. (FN,B5) and (3.8) as where ω ϕ (θ A ) is defined by ϕ * h AB = ω −2 ϕ h AB and ϕ * is the pullback. The quantity ω ϕ is determined by the boost part of the Lorentz transformation ϕ and is given explicitly by ω ϕ = cosh ψ − n · m sinh ψ, where n = (sin θ cos φ, sin θ sin φ, cos θ), ψ is the rapidity parameter of the boost and m is a unit vector giving the direction of the velocity of the general frame with respect to the canonical frame. The 4-momentum in the general frame is now given from Eqs.
Here the overbar denotes the value of this function in the general Bondi frame. The transformation law (B.4) can be simplified by defining the new quantity where the potential Φ for the electric parity piece of the shear tensor C AB is defined in Eq. (3.9). Now in the canonical BMS frame, S A will coincide withN A and will be purely l = 1 and of magnetic parity, encoding the intrinsic spin of the spacetime. Combining Eqs. (B.4), (B.5), (FN,B5), (FN,B6) and (FN,B8) yields that in the general BMS frame this function will bē where β 0 consists of the l = 0, 1 components of β. Combining this with a barred version of Eq. (B.5) giveŝ The last three terms in this equation are determined from the intrinsic spin, 4-momentum and from the Poincaré transformation (ϕ, β 0 ) relating the two frames, so they are determined by the linear and angular momentumP α andJ αβ . It follows thatN A and the superspin and super center-of-mass charges are determined fromJ αβ andP α and the electric parity potentialΦ, in a general Bondi frame in a stationary region 21 . We next derive the explicit form of the super center-of-mass charges, expanding to second order in the velocity of the boost. Since the metric is stationary, shifting the supertranslation β by a constant times ω ϕ does not affect the metric in the general Bondi frame, from Eqs. (FN,2.12). Therefore without loss of generality we take the l = 0 component of β to vanish, and we parameterize the translation β 0 as β 0 = β i n i . From Eq. (B.7) we can now read off the orbital angular momentum (3.15b) and the l = 2 super center-of-mass charge (C.4) in the general frame, using Eqs. (B.1) and (B.2) to expand to expand in powers of the velocity v = m tanh ψ of the boost: Here we have used the definition (3.10), J ij is the intrinsic angular momentum in the canonical Bondi frame, and the angular brackets < . . . > denote the symmetric tracefree projection. Combining Eqs. (B.8) to eliminate β i , dropping the intrinsic angular momentum terms and using Eq. (B.3) gives A similar calculation for l ≥ 3 using Eq. (3.7) gives

C Choice of basis of algebra of charges
In this appendix we discuss two different versions of the super center-of-mass charges.
To explain these versions, it is useful to distinguish between two different kinds of time evolution of the charges. The first is just the kind discussed in Sec. 3.3, associated with evaluating the charges as surface integrals on cuts of I + of the form u = constant, and varying u.
A second kind of time evolution is associated with the choice of basis in the algebra of asymptotic symmetries. Consider for example the orbital angular momentum J 0i that is associated via Eqs. (3.4) and (3.5) with the boost symmetry generator ξ = D A n i ∂ A − un i ∂ u . This choice of boost symmetry is associated with a particular choice of origin of the retarded time coordinate u (which in the body of the paper we took to be the time of formation of the black hole, see Sec. 3.4). However, by conjugating the symmetry generator with a time translation u → u − u 0 where u 0 is a constant, we can obtain the new boost symmetry generator (C.1) We denote the corresponding charge by J 0i (u, u 0 ), where the first argument reflects the dependence on the cut of I + and the second argument the choice of generator. The dependence on u 0 is given by changing u 0 amounts to a change of basis in the algebra of symmetry generators or equivalently in the algebra of charges. The center of mass at retarded time u, given by encodes both types of time dependence. It is this quantity and not the charge J 0i (u, 0) that most directly enters into the metric at retarded time u, from Eqs. There is an exactly analogous story for the super center-of-mass charges [48]. We specialize for simplicity to the quadrupole l = 2 case. Consider the superboost symmetry generator given by (3.4) for α = 0, Y A = D A (n i n j ), ξ = D A (n i n j − δ ij /3)∂ A − 3u(n i n j − δ ij /3)∂ u . (C.4) We denote the corresponding charge (3.5) by J ij (u), a symmetric traceless tensor [cf. the discussion around Eq. (3.7) above]. As before by conjugating with a time translation we can obtain a new superboost symmetry generator, given by Eq. (C.4) with u replaced by u − u 0 , and we denote the corresponding charge by J ij (u, u 0 ). The dependence on u 0 is given by, from Eqs. is the l = 2 supermomentum charge. As before we can define a super center-of-mass quantity that incorporates both types of time evolution, and which is the quantity that appears most directly in the metric, via We will call the charge (C.5) the super center-of-mass charge, and the quantity (C.7) the comoving super center-of-mass, following Compère [48].
In the special case of stationary regions of I + the charges J ij (u, u 0 ) and J 0i (u, u 0 ) are related by the formula (B.9) derived in Appendix B: J ij (u, u 0 ) = − 3 5 m 0 Q ij (u) 1 + O(v 2 ) + 12 5m 0 J 0<i (u, u 0 )P j> (u) 1 + O(v 2 ) . (C.8) Here m 0 is the rest mass associated with the Bondi 4-momentum, v is the velocity of the Bondi frame with respect to the canonical Bondi frame, Q ij is defined in Eq. (3.10), and the angular brackets denote symmetric tracefree projection. This formula is consistent with the transformation laws (C.2) and (C.5) because we have in stationary regions

D Derivation of changes in charges in stationary-to-stationary transitions
In this appendix we derive the formulae (3.12), (3.16) and (3.17) for the changes in BMS and extended BMS charges in stationary to stationary transitions in terms of fluxes to null infinity of mass-energy or gravitational-wave energy. A similar analysis in a different notation can be found in Ref. [48]. The change (3.12) in Bondi 4-momentum is obtained by multiplying Eq. (FN,4.3) by (1, n), integrating over solid angles, and noting that the l = 0, 1 components of Φ vanish by definition.
For general superspin and super center-of-mass charges, we denote by Q Y the charge (3.5) specialized to α = 0. To derive the change in this charge we first define the subleading memory observables ∆Φ = We now multiply by Y A and integrate over solid angles and over u. The left hand side then becomes the change ∆Q Y in the charge, from Eq. (3.5). On the right hand side, the last term gives a vanishing contribution since we assume the stress energy tensor vanishes at u = u 1 and u = u 2 . The second and third terms can be simplified using the definition (3.14) of u-weighted energy flux ∆E, while the first, fourth, fifth, sixth and seventh terms can be similarly simplified using the definition (3.13) of angular momentum flux ∆E A , making use of Eq. (FN,3.23). The eighth, ninth and tenth terms can be written in terms of the subleading memory observables (D.1) using Eqs. (3.3) and (3.9), and using R ABCD = h AC h BD − h AD h BC . The final result is 23 where the differential operator D is given by Eq. Finally we turn to deriving the change (3.17) in the electric parity potential Φ for the shear tensor C AB , which encodes the super center-of-mass charges. We multiply Eq.