Qubits on the Horizon: Decoherence and Thermalization near Black Holes

We examine the late-time evolution of a qubit (or Unruh-De Witt detector) that hovers very near to the event horizon of a Schwarzschild black hole, while interacting with a free quantum scalar field. The calculation is carried out perturbatively in the dimensionless qubit/field coupling $g$, but rather than computing the qubit excitation rate due to field interactions (as if often done), we instead use Open EFT techniques to compute the late-time evolution to all orders in $g^2 t/r_s$ (while neglecting order $g^4 t/r_s$ effects) where $r_s = 2GM$ is the Schwarzschild radius. We show that for qubits sufficiently close to the horizon the late-time evolution takes a simple universal form that depends only on the near-horizon geometry, assuming only that the quantum field is prepared in a Hadamard-type state (such as the Hartle-Hawking or Unruh vacua). When the energy difference, $\omega_\infty$, between the two qubit states (as seen by an observer at infinity) satisfies $\omega_\infty r_s \ll 1$ this universal evolution becomes Markovian and describes an exponential approach to equilibrium with the Hawking radiation, with the off-diagonal and diagonal components of the qubit density matrix relaxing to equilibrium with different characteristic times, both of order $r_s/g^2$.


Introduction and Summary
Making reliable predictions can be difficult at the best of times. But reliably predicting behaviour at very late times is notoriously hard. What makes it difficult is the inevitable breakdown of perturbative methods that happens at late times; a huge handicap given that perturbative methods dominate a theorist's intellectual toolbag.
Perturbative methods break down for a simple reason: if a Hamiltonian can be written H = H 0 + gH 1 for some small dimensionless parameter g, then there is always a time beyond which the time-evolution operator U (t, t 0 ) = exp[−iH(t − t 0 )] is not well-described by perturbing in g. The time where this breakdown occurs scales as an inverse power of g, but is eventually exceeded no matter how small g might happen to be. Like death and taxes, perturbative failure is just a matter of time.
This might not be bothersome if such late times were never of interest. However many important physical processes occur on long time-scales like these. For example, even when individual photons interact weakly with individual atoms, phenomena like refraction and reflection (where 100% of photons scatter in one direction or another) occur on time-scales long enough to invalidate perturbing in electromagnetic interactions. Thermalization is another phenomenon whose time-scales are very large and scale inversely with coupling strengths like g.
In this paper we explore similar late-time issues for interacting quantum systems moving in gravitational fields. That similar phenomena must exist -particularly in the presence of horizons -is clear given the thermal nature of quantum fields in these spacetimes [1][2][3][4][5][6]. Test probes should be expected to thermalize in such environments, and any description of this process should share all of the late-time complications that thermalization calculations always have [7][8][9][10][11][12][13][14]. In this paper we show this is true for quantum systems exterior to a Schwarzschild black hole, extending our own earlier work that does so for spacetimes with Rindler [15,16] and de Sitter horizons [17].
The reason for doing so is not because this kind of thermalization is soon likely to be observed. On the contrary, it is worth doing because the tools used are informative in their own right. In particular, they show how standard techniques used to describe late-time behaviour in optics and thermal physics apply equally well in gravitational settings [18][19][20][21][22][23][24][25][26][27][28]. This makes them potentially relevant to late-time puzzles known to occur in gravity, such as the problem of secular growth in cosmological spacetimes  (for reviews see [53,54]) and to problems like information loss [55] or 'firewall' problems [56,57] in black-hole physics (for reviews see [58,59]).
In this paper we compute the late-time evolution of a two-level quantum system (i.e. a qubit or Unruh detector [5,60,61]) that hovers at fixed radius r = r 0 above the event horizon of a Schwarzschild black hole while interacting with a quantum scalar field. We do so perturbatively in the dimensionless coupling strength g with which the qubit interacts with the quantum field. We show that if ω ∞ is the splitting of the two qubit energy levels at distances far from the black hole, and if ω ∞ r s 1 where r s = 2GM is the usual Schwarzschild radius, and if the scalar-field mass also satisfies m r s 1, then such a qubit has universal late-time behaviour (for t > ∼ r s /g 2 ) provided that it sits sufficiently close to the event horizon: 0 < r 0 − r s r s . Not surprisingly, this universal evolution describes the evolution of the qubit towards an asymptotic thermal state whose temperature equals the Hawking temperature T = T H := (4πr s ) −1 . Perhaps more surprisingly we show that this approach to equilibrium is also very robust, occurring exponentially with two different thermalization time-scales proportional to ξ = 4π tanh (2πr s ω ∞ ) g 2 ω ∞ 8π 2 r s g 2 + · · · since ω ∞ r s 1 . (1.1) This evolution is robust in the sense that it depends only on the qubit/field coupling strength, g, and on the background geometry for any quantum state whose Wightman function has the standard 'Hadamard' form [62][63][64] at small field separations: i.e. it satisfies (3.13), reproduced here as where σ(x, x ) is half the square of the geodesic distance between spacetime points x and x . In particular, eq. (1.1) applies equally well if the quantum field is prepared in either the Hartle-Hawking or Unruh vacua, and is independent of the scalar-field mass in the mass range m r s 1. It has been known for some time that Hadamard behaviour suffices for deriving the steady-state Hawking flux around Schwarzschild black holes [65], and our results extend this conclusion to the approach to equilibrium for quantum probes. We remark in passing that our results differ from early -and some recent -calculations of Unruh detectors in Schwarzschild geometries [66][67][68][69][70][71], which often compute qubit excitation rates, finding results that differ when the field is prepared in different states (such as the Hartle-Hawking or Unruh vacua). These calculations usually compute the rate with which a qubit is excited out of its ground state, as opposed to the qubit's latetime approach to its asymptotic thermal state (as is computed here). Although the excitation rate can be accessed perturbatively in g, more effort is required to obtain the approach to equilibrium since the time-scale involved is of order r s /g 2 .
We are able to make reliable predictions using arguments of Open Effective Field Theories (Open EFTs) [24,26,28]. As is explained in more detail in [16], these recast techniques from elsewhere in physics into an effective field theory language that is easily adapted to gravitational systems. In essence these arguments have a renormalizationgroup like structure: one sets up a differential evolution equation for the object of interest (in this case the reduced density matrix for the qubit) whose domain of validity is larger than the integrated evolution from which it is derived. That is, one explicitly evolves the system using perturbation theory starting from an arbitrary initial time, t 0 . Although perturbative evolution can only be used to evolve a limited way forward in time, say from t 0 to t 1 , within this window the result can be differentiated with respect to time to derive a differential evolution equation.
If this evolution equation itself makes no specific reference to t 0 then the same construction could equally well be used to derive the same evolution equation starting at t 1 , with perturbative validity out to t 2 , and again starting at t 2 and so on. Whenever this can be done the solutions to the differential evolution equation can be valid on the union of each of these derivation intervals. If g 1 is the small perturbative expansion parameter then this process ends up resumming all orders in g 2 t, say, but neglecting contributions in the evolution 1 that are of order g 4 t. As a result the solutions found this way can be trusted even when t ∼ O(r s /g 2 ).
One reason to explore the simple qubit systems considered here is to make this construction very explicit, making it easier to understand. The starting point for the argument is the Nakijima-Zwanzig equation [72,73], which is a general evolution equation for the reduced 2×2 density matrix, (t), of the qubit. It is obtained by tracing over the Liouville equation describing the evolution of the full qubit/field system, and then eliminating that part of the density matrix that describes the non-qubit degrees of freedom. The result is an integro-differential evolution equation that is useful because it refers only to the qubit's reduced density matrix and not to the other degrees of freedom, which appear only implicitly through correlation functions of H int . Although the Nakajima-Zwanzig equation does not in itself automatically allow perturbative time-evolution to be extended out to very late times, it provides a useful starting point for identifying situations where this can be done.
As is true for most effective field theories, relative simplicity comes only when there is a hierarchy of scales that can be exploited. The important hierarchy arises in this case if the field correlation function H int (t)H int (t ) falls off to zero for |t − t | > ζ, for some characteristic time-scale ζ. In this case the useful hierarchy arises when exploring time-evolution over much longer time-scales ∆t ζ. Access to late times can happen if the Nakajima-Zwanzig equation remains sufficiently simple once expanded in powers of this ratio ζ/∆t.
The qubit example studied here shows in detail how this can happen: the leading terms in the Nakajima-Zwanzig equation become Markovian, in the sense that ∂ t (t) depends only on (t) and not on the details of its past history prior to time t. Markovian behaviour of this form emerges for qubits near a black hole once ∆t r s (at least this is true when the energy difference ω ∞ between the two qubit energy levels -seen by a static observer far from the black hole -satisfies ω ∞ r s 1), Evolution to all orders in g 2 t is then described by a Lindblad equation [74,75]. (Some implications of Lindblad evolution in Schwarzschild geometries are also explored in [76][77][78][79][80][81][82].) By deriving the Lindblad equation as a limit of the Nakajima-Zwanzig equation for this system, we are able to assess its domain of validity.

This paper
The rest of this paper is organized as follows. The next section, §2, sets up the system whose late-time near-horizon evolution is to be computed. In particular §2 defines our qubit/quantum-field system for static spacetimes, and then briefly explores the properties of qubit trajectories that hover at fixed positions just above a Schwarzschild black hole.
§3 follows this with a brief description of how reduced density matrices are evolved in open systems, describing the Nakajima-Zwanzig equation whose solutions govern the qubit's late-time behaviour. Since at lowest nontrivial order the quantum field enters into the qubit evolution only through its Wightman function, we also summarize in §3 the near-horizon form for this function for field states that satisfy the Hadamard form for small separations.
Finally, §4 shows how the near-horizon limit of the Wightman function allows the Nakajima-Zwanzig equation to be approximated by a Markov process, describing the late-time exponential decay towards a Hawking-temperature thermal state. The time-scale for this approach to equilibrium is computed for qubits asymptotically close to the horizon, and found to be universal in the sense that it is determined only by qubit properties and the black-hole geometry. Provided m r s 1 this rate is largely independent of the details of the quantum field, and assumes only that it is prepared in a Hadamard state. In particular the approach to equilibrium is the same when the field is prepared in either a Hartle-Hawking or Unruh state.

Qubits in Schwarzschild
This section sets up the framework -a qubit/field system and the spacetime through which the qubit moves -that is used to perform the calculations to follow.

The qubit/scalar system
The system whose evolution we follow consists of a real massive scalar field φ(x) coupled to a single two-level qubit through the action S = S B + S Q + S int , where S B describes a self-interacting quantum field within a background metric 2 given as in (2.20) or (2.22). For this paper we neglect self-interactions (λ = 0), though we briefly comment in the conclusions on how things can change in their presence. The coupling ξ plays no role because Schwarzschild is Ricci flat, and for reasons to be clear below the mass m is assumed to satisfy m r s 1.
The free qubit action is given by 3 where z i (s) are classical Grassman variables (with i = 1, 2, 3) withż i := dz i /ds and z i := δ ij z j . The quantities Z, ω 0 andω and u i are real parameters (with u i u i = 1), with Z eventually absorbed into the z i to obtain a convenient normalization that simplifies later formulae. 4 The integral over d 4 x is trivially done using the delta-function, and reveals that the integration is over a specific timelike world line x µ = y µ (s), along which the qubit moves through the ambient spacetime. 5 Here s is a parameter along this world line, and the quantityẏ is what is required to ensure that S Q is invariant under reparameterizations of s. It is usually convenient to fix this freedom by choosing proper time, τ , along the curve as the parameter, in which caseẏ 2 = −1.
Interactions beween z and φ are assumed to take the form whereĝ and n i are real coupling constants, with n i n i = 1, and our analysis is ultimately performed perturbatively inĝ.

Quantization
Working in the interaction picture, quantization of z i and φ is performed as if they did not interact, with interactions included in powers of g (and λ) once time-evolution is evaluated. For z i this quantization goes through as usual, keeping in mind this is a constrained system. To see why, recall that the canonical momenta are given by Terms linear in z i do not appear in this action because we require the classical Grassmann action to be Grassman-even. 4 Hats appear on couplings likeω andĝ to distinguish them from the corresponding quantities once appropriate powers of Z have been absorbed into z i . 5 For real systems y µ (s) is itself a dynamical variable to be quantized, but for simplicity we ignore this complication here (treating the qubit trajectory as being specified), since it does not affect our later discussions.
Because this cannot be solved forż i as a function of the p j it is instead regarded as a constraint: The canonical quantization conditions turn out to imply that the anticommutator of p i and z j is proportional to δ i j , and the parameter Z can be chosen to ensure that the quantum version of the Grassmann condition becomes The space of quantum states furnishes a representation of this algebra, and for a 2-level qubit this representation is two-dimensional. The required operator representations for the z i therefore are the Pauli matrices With this choice, properties of the Pauli matrices ensure that ijk z j z k = 2iz i , and so defining coordinates so that u · σ = σ 3 then allows (2.6) to be written explicitly as where I is the 2×2 unit matrix, and hats are dropped on ω when variables are normalized so that (2.7) holds. Eq. (2.9) reveals the free-qubit energy eigenvalues to be ω 0 ± 1 2 ω and so ω 0 is their mean energy while ω (which we take to be positive) gives their level splitting (as measured by an observer whose time is the qubit's proper time, τ ).
With this representation for z i the qubit/field interaction (2.4) becomes where the second equality specializes to the case where n i is perpendicular to u i (and so can be chosen to lie along the '1' axis).
In the absence of couplings (λ = g = 0) the scalar field is quantized in the usual fashion for a static curved space [6]. The interaction-picture field equation is and so for a static geometry one expands where u n simultaneously satisfies L t u n = −iω n u n and eq. (2.11), where L t is the symmetry generator in the timelike direction along which the metric is static. Canonical commutation relations imply the creation and annihilation operators satisfy [a n , a * m ] = δ mn . As mentioned earlier, for applications to Schwarzschild our interest is in small masses, 6 m r s 1, and Ricci-flatness makes the term ξR drop out of subsequent discussion. The scalar field hamiltonian (including self-interactions) is easily computed in the presence of any spacetime metric of the form where g at = 0 and both g tt = −f and the spatial metric g ab = γ ab are t-independent.
The hamiltonian density (for the generator of evolution in t) is then given by where L is the lagrangian density from (2.1) and the canonical momentum is where (2.1) is again used, together with the metric (2.13), to evaluate the derivative. Therefore the scalar-field hamiltonian density, H, becomes

Total energy
We can now assemble everything to identify the total field/qubit hamiltonian, leading to the following sum: where the 'free' hamiltonian is the sum of the two free hamiltonians constructed above Here I is the unit operator in the scalar-field state-space, while H and h are as given in eqs. (2.16) and (2.9). The factor dτ /dt is required because h generates translations in proper time τ while H 0 (and H) are defined to generate translations in t.
Since the interaction Lagrangian does not involve time derivatives its contribution to the hamiltonian is simple to write down, starting from (2.10) In what follows we compute the implications of this interaction out to second order in g. The nature of this perturbation theory depends on the relative size of ω and the O(g 2 ) corrections to the qubit energy levels, and for simplicity we work in the regime where these corrections are much smaller than the qubit's zeroth-order level splitting, a restriction that eventually leads to the parameter conditions summarized in Table 1.

Near-horizon geometry
Our interest lies in the near-horizon limit of the exterior of a spinless black hole, defined in Schwarzschild coordinates by where r > r s , where r = r s := 2GM defines the event horizon. These coordinates are useful because they fall into the category defined in (2.13), with t being the static time on which the metric does not depend, and so L t = ∂ t . The (Kretchsmann) curvature invariant for this geometry is R ρστ κ R ρστ κ = 12r 2 s /r 6 , and so is nonsingular at r = r s . Schwarzschild coordinates famously break down at the horizon, in the vicinity of which Kruskal coordinates, (T, X, θ, φ), defined by are more useful. In terms of these the line element becomes where now r = r(X, T ) is the implicit function of X and T given by solving and so is given by and W is the Lambert W function defined by W(z) exp[W(z)] = z (which has a unique real solution for z > 0 and two real branches for −e −1 < z < 0). Although these coordinates are well-behaved at the horizon, the spatial geometry at fixed T is T -dependent. This is a reflection of the coordinate-independent statement that the metric is only static outside of the horizon (where it can be rewritten as (2.20)). Inside the horizon the metric's symmetry directions are all spacelike.

Hovering world-lines
The qubits whose late-time evolution we follow are chosen to hover at fixed r = r 0 > r s above the event horizon. This is clearly not a geodesic and so the qubit can be maintained along this trajectory through the action of some non-gravitational force, whose detailed nature need not concern us here.
Consider two events on this world-line that are distinguished by two values, τ 1 and τ 2 , of the qubit's proper time; separated by ∆τ := τ 2 − τ 1 > 0. These two events are separated by a redshifted time ∆t for hovering observers situated at spatial infinity, with (2.20) implying that ∆τ := where y µ 0 (t) denotes the curve along which only t varies, with r, θ and φ all fixed. Notice, for future use, that ∆τ ∆t when 0 < r 0 − r s r s . Both of these intervals differ from the geodesic separation of these two points, ∆s := where the subscript 'g' indicates that integration is evaluated along the geodesic y µ g (t) that satisfiesÿ Because this is a geodesic it must describe the longest time interval as measured along any timelike curve that connects the two events, so ∆s > ∆τ . Such a geodesic is possible if dr/dt(t 1 ) > 0 is chosen appropriately, since then any freely falling body initially moves radially away from the horizon before eventually turning back and falling into the black hole.
In what follows it proves more convenient to work with the Synge world function, σ(x 1 , x 2 ), defined for timelike geodesics by [87][88][89] σ(x 1 , since this has an integral form that is easier to manipulate (see Appendix A). For later purposes we are interested in a limit that simultaneously has a small invariant interval, ∆s r s , but corresponds to late times ∆t r s . We remark in passing that the above formulae show that both of these can be simultaneously true provided we pick r 0 > r s sufficiently close to the horizon so that Being close to the horizon suffices because the curve that hovers at fixed r = r s is a geodesic, although it is a null geodesic -for which ∆s = 0 -rather than a timelike one.

Time evolution in open systems
We return now to the evolution of the qubit that hovers just above the horizon while interacting with the quantum field. Our interest is in how this qubit responds to the fluctuations of the quantum field, and in how this response becomes universal in the late-time limit very near the horizon.
Since it is only the qubit's behaviour that is to be predicted, it is convenient to trace out the scalar field from the system density matrix, and work instead only with the qubit's 2 × 2 reduced density matrix, defined as where ρ(t) is the total -i.e. the combined field/qubit -density matrix, and the trace is over the scalar-field part of the Hilbert space. When needed we assume the field and qubit to be initially uncorrelated, where 0 defines the initial qubit state and Ω is the density matrix for the quantum field. Three commonly made choices for Ω might be the Hartle-Hawking state, Ω H := |H H|, the Unruh state Ω U := |U U| or the Boulware state, Ω B := |B B|. These are all pure states that are candidate vacua for the field, with |H corresponding to the vacuum in the presence of a black hole that is in equilibrium with a bath of radiation prepared at the Hawking temperature, while |U is the late-time vacuum for a black hole that forms in isolation. Time evolution for is in principle determined by the evolution of the full system's density matrix, which in the interaction picture satisfies where V (t) := e iH 0 t H int e −iH 0 t . Eq. (3.3) has a standard perturbative solution given the initial condition ρ(t i ) = ρ 0 .
As is discussed at great length elsewhere -see for example [16,17,28] and references therein -there are two major obstacles to using eqs. (3.3) or (3.4) to predict (τ ) at late times.
1. At first sight one could trace (3.3) over the scalar-field sector to obtain ∂ t , but the result is hard to solve for (τ ), because the dependence of the right-hand side on is only implicitly given through its dependence on the full density matrix ρ. The Nakajima-Zwanzig equation [72,73] provides the solution to problem (1) above, and this is useful because the result shows how to find solutions that are not afflicted by problem (2), in that they allow reliable perturbative predictions even when t is so large that g 2 t cannot be neglected relative to r s .

The Nakajima-Zwanzig Equation
The logic of the Nakajima-Zwanzig equation is to project the full density matrix onto the reduced density matrix and its complement: = P(ρ) and Ξ := Q(ρ) (3.5) for some projection operator P 2 = P and the second definition uses Q := 1 − P = Q 2 .
Since the time-evolution equation (3.3) for ρ is a linear equation it can be turned into a pair of coupled linear evolution equations for the two quantities and Ξ. Eliminating Ξ from this system gives the Nakajima-Zwanzig equation: an evolution equation that involves only , but is nonlocal in time due to the elimination of Ξ. Because this is essentially a linear problem, it can be solved in great generality [72,73].
As applied to the current example, following identical steps as given in [16,17] leads to the following result at second order in the coupling g: where I (τ ) = e +ihτ (τ ) e −ihτ is the reduced density matrix in the interaction-picture representation and so similarly m I (τ ) := e +ihτ σ 1 e −ihτ (and conventions generally follow [16,17]). At this point several features of (3.6) bear explanation.
• First, notice that eq. (3.6) gives the evolution of as a function of proper time along the qubit trajectory, and does so despite its derivation starting from the Liouville equation ( • Second, in expression (3.6) the quantity G Ω (τ, s) represents the scalar-field Wightman functions evaluated at two places along the qubit trajectory, y µ (s), and we use the property G Ω (s, τ ) = G * Ω (τ, s) for the Wightman function of a real scalars. These are fairly complicated functions when evaluated in Schwarzschild spacetime and are usually given implicitly in terms of a sum over mode functions [60,66,76,[90][91][92][93][94][95][96]. In what follows we choose a near-horizon trajectory along which they take a simple form.
• Third, the final term in (3.6) comes from a counter-term interaction, obtained by replacing ω → ω bare = ω + δω in h, with δω regarded as being O(g 2 ). This counter-term is required because the qubit/field interaction shifts the inter-level energy spacing, and so makes the parameter ω appearing in h no longer equal to this spacing. If the parameter in h is therefore instead called ω bare then ω remains the physical spacing of qubit levels if δω is chosen to cancel the order g 2 qubit energy shift. 7 • Finally, notice that (3.6) would agree with the time derivative of (3.4) if in the right-hand side one were to replace I (s) with its initial condition I 0 . Furthermore, such a replacement at face value seems to be compulsory, because the difference between I (s) and I 0 is higher order in g. It is this assumption that I (s) and I 0 are only perturbatively different that breaks down at very late times, and when it does it is (3.6) that is the more reliable equation.
Since is a Hermitian 2 × 2 matrix with unit trace, its elements 21 = * 12 and 22 = 1− 11 can be eliminated from (3.6) to leave the following two decoupled evolution 7 A bonus of this definition is that δω also automatically cancels an ultraviolet divergence that arises in the computed energy shift. We note in passing that it can happen that the appearance in the above equations of the oscillatory factors e ±iωs and e iωτ can complicate the construction of their solutions. Such terms can be removed from an ordinary differential equation by standard changes of dependent variable, which in the present instance amount to returning to the Schrödinger picture. The result in Schrödinger picture is

Near-horizon Wightman function
So far the description of qubit evolution has been quite general, with little said about the specific field state Ω or about the details of the qubit trajectory. Application of this formalism to a qubit near a black hole requires filling in some of this detail, starting with some information about the scalar-field Wightman function in a Schwarzschild geometry.

Hadamard correlation functions
As mentioned earlier, for generic trajectories the scalar field Wightman function can be quite complicated, even for comparatively simple states like the Hartle-Hawking, Unruh or Boulware vacua [60,66,[90][91][92][93][94][95][96]. One of our central points is that the late-time evolution very close to the horizon does not depend on which of these choices for field state is made, with universal predictions relying only on the state being 'Hadamard', in the sense that the Wightman correlation function has -in four spacetime dimensions -the following limit as x → x [64,[97][98][99][100]: 14) and σ(x, x ) the so-called Synge world function [63,87] that is equal to half the square of the geodetic length between x and x (see Appendix A). Here T is any future-increasing function of time, and → 0 + a small-distance regulator with dimensions of length that appears in the above formula so that G Ω (x, x ) satisfies the correct temporal boundary conditions. The quantities ∆(x, x ), V (x, x ) and W Ω (x, x ) are biscalar functions that are symmetric in x ↔ x , and regular in the limit that x → x . The renormalization length scale L > 0 is introduced on dimensional grounds, and different values for L can be absorbed into the precise definition of W Ω (x, x ). The subscript Ω on W Ω is meant to emphasize that its detailed form depends on the state Ω [101]. The same is not true of the functions ∆(x, x ) and V (x, x ), which are universal in the sense that they depend only on the geometry of the spacetime (and -in the case of V (x, x ) -on parameters like the mass of the field).
What this says is that the leading part of the coincident limit of G Ω (x, x ) is universal in curved space, and shares in particular the singularity structure also found in flat space. The Hadamard form expresses the physical condition common to all effective field theories [28] that states that the details of very high-energy field modes are irrelevant provided because for slowly changing backgrounds they are prepared within their adiabatic vacuum. This amounts to a quantum variant of the principle of equivalence: modes with wavelengths much shorter than the local radius of curvature do not 'know' that they are in curved space.
Because they depend only on local properties, there is a general procedure for computing the geometric functions V (x, x ) and ∆(x, x ) in the coincident limit, for which σ(x, x ) → 0 [63,102]. For a real massive scalar field evaluated on a Ricci-flat spacetime (like the Schwarzschild geometry) they have the form where ∂ µ σ = σ µ and σ µ = g µν σ ν obey the relation σ µ σ µ = 2σ. Terms written O(σ 3/2 ) are those containing three or more factors of σ µ . For massive fields it is conventional to choose the form of W Ω (x, x ) so that L 2 = 2/m 2 , so that Of the vacuum states described above, the Hartle-Hawking [103] and Unruh vacua [104] are both Hadamard states, and so share the same values for ∆(x, x ) and V (x, x ) but not for W Ω (x, x )). The Boulware vacuum is not, however, as can be seen from its singular form for the stress-energy tensor at the horizon [96].
In practice the leading behaviour suffices for our purposes, which means we may use ∆(x, x ) 1 and drop V (x, x ) in the applications to follow, leaving the result that applies when |σ(x, x )| is much smaller than both r 2 s and m −2 .

Evaluation for qubits hovering near the horizon
What is special about the small-σ(x, x ) limit is that it applies not just as x → x , but also when x and x are generic points situated sufficiently close to a null geodesic. Small σ(x, x ) should apply in particular for any two points hovering at a fixed position (r, θ, φ) = (r 0 , θ 0 , φ 0 ) just outside the Schwarzschild event horizon, with σ(x, x ) → 0 as r 0 → r s . The function σ(x, x ) is evaluated in this limit in Appendix A for points on such a hovering trajectory as a function of their separation ∆t in Schwarzschild time, with the result [70,81,82]) in the limit σ(x, x ) → 0. It is important that (3.18) remains valid even if ∆t r s , provided that r 0 is chosen close enough to r s to ensure that |σ(x, x )| r 2 s . (The validity of this approximation in the regime ∆t r s is verified numerically in Appendix A.) For separations for which (3.18) applies, eq. (3.17) states that the Wightman function for any Hadamard state has the form

Universal late-time near-horizon evolution
This section ties everything together to obtain a closed-form expression for the two universal thermalization time-scales that arise for qubits hovering asymptotically close to the horizon. The result is surprisingly simple because of a somewhat paradoxical result: the simplicity occurs because in the near-horizon limit one can exploit the Wightman function's small-σ(x, x ) Hadamard form (3.13). This is paradoxical because thermalization occurs in the limit of very long time separations, ∆t r s . The coexistence of these two limits is possible only because of the enormous time-dilation that relates static clocks running very near the horizon and those far from the black hole; two nearhorizon events separated by a small geodesic separation can look to a distant observer like they are separated by very large times.

The near-horizon Nakajima-Zwanzig equation
The starting point is the interaction-picture Nakajima-Zwanzig equations (3.8) and Our later interest is in late times as seen by an observer far from the black hole, so changing coordinates τ = t/γ gives and where for convenience we define the red-shifted qubit gap as seen by observers far the black hole and perform a similar scaling of the O(g 2 ) counter-term δω ∞ := δω(1 − r s /r 0 ). Finally, F denotes the scaled Wightman function In the small-σ(x, x ) limit inspection of (3.19) shows that F has the simple asymptotic form F(∆t) − 1 64π 2 r 2 s (sinh [∆t/(4r s )] − i /(4r s )) 2 , which is identical to the analogous result for the massless Rindler correlation function found in [16] once one replaces r s → 1/(2a). Recall from Appendix A that this asymptotic form for F(t) is valid so long as |σ(x, x )| r 2 s and so applies when and so in particular there is always an r 0 > r s sufficiently close to the horizon for which this is satisfied, no matter how large ∆t/r s happens to be.

The late-time Markovian approximation
From here on the story evolves much as it did in the Rindler example considered in [16], by virtue of the similarity between (4.8) and its counterpart for an accelerated qubit in flat spacetime.
In particular eqs. (4.4) and (4.5) greatly simplify when I is slowly varying (compared with the light-crossing time of the black hole, r s ) and we focus on t r s , because in this case the sharply peaked form for F(t) allows the upper integration limit to be taken to infinity, and implies that a Taylor expansion in the integrand of ij (t − s) in powers of s converges very quickly.
After choosing δω to cancel the field-induced shift in qubit energy -which means picking with D S defined by (4.14) below -these steps lead (at face value) to the following approximate evolution equations (see [16] for details) and 12) in which the quantities C S and D S are defined by where ψ (0) (z) = Γ (z)/Γ(z) is the digamma function [106].

Control over approximations
The words 'at face value' are added above eqs. (4.11) and (4.12) because the term involving D S must actually be dropped in the above if we are consistent. The reasons for this lie in the size of the deviations from the leading approximation, and the assumptions that must be made in order to neglect them. We briefly summarize the issues, following closely the discussion in [16,17]. A side effect of this observation -together with the energy shift (4.10) -is the elimination of all singular dependence 9 in the limit → 0 that enters through eq. (4.14).
8 A flat-space analog of D S is computed in [16] for generic field masses m = 0 and with the replacement a → 2/r s (using the different notation ∆ M there). Equation (4.14) follows as the m → 0 + limit of this function (this same function is evaluated in [105]). 9 For the purposes of estimating the size of different contributions we take here to be much smaller than other scales, but not infinitely small so that logarithms of cannot overwhelm powers of g 2 .
There are two kinds of approximations to consider -one convenient and one essential. The issue of convenience concerns the relative size of the qubit splitting ω and the generic size of field-driven corrections to this splitting. Assuming ω ∞ to be much larger than the corrections to ω induced by the interactions with the field simplifies calculations by allowing use of non-degenerate methods. In terms of the functions C S and D S this condition requires Table 1 displays the asymptotic form for these two quantities in the limit of large and small ω ∞ r s , showing that they require ω ∞ r s not to be taken smaller than g 2 /4π. The large-and small-ω ∞ r s asymptotic forms for the two quantities that must be small to work with nondegenerate perturbation theory (see [16]).
The essential approximation is the one that makes the Markovian evolution dominate the Nakajima-Zwanzig evolution. To see what this involves, recall that the Markovian approximation is derived from the Nakajima-Zwanzig equation by Taylor expanding I ij (t − s) I ij (t) − s˙ I ij (t) + · · · inside the integrands of equations (4.4) and (4.5), where t r s is used to take the upper limit of integration to infinity (given the exponential falloff of f (s) for s r s ). The size of the s˙ I ij (t) term characterizes the size of deviations from the Markovian limit, and we evaluate it to understand what demands are made on the free parameters of the model by the requirement that these be small. Physically this amounts to requiring the evolution time-scale of I ij to be large compared with the domain of support of the rest of the integrand. The quantitative conditions are obtained self-consistently, by evaluating˙ I ij assuming the time dependence is given by (4.11) and (4.12), whose integration implies I 11 (t) = 1 e 4πrsω∞ + 1 and I 12 (t) = e −g 2 C S t 12 (0) + * 12 (0) . (4.18) Differentiating this to find˙ I ij then allows the˙ I term to be computed in equations like (4.16), and requiring the result to be negligible relative to the leading term for all of the integrals appearing in eqs. (4.4) and (4.5) requires the following four quantities all to be negligible: The first two of these are required whenever the derivative in˙ ij is of order g 2 C S while the second two arise when it is order ω ∞ . [Differentiation with respect to ω ∞ arises from use of identities like s cos(ω ∞ s) = (d/dω ∞ ) sin(ω ∞ s) in equations like (4.16).] The large-and small-ω ∞ r s asymptotic forms for the four quantities that must be small to believe the Markovian approximation to the Nakajima-Zwanzig equation (see [16]). Notice that r s ω ∞ 1 is incompatible with Markovian evolution. Table 2 displays the asymptotic behaviour for these four quantities in the limits where r s ω ∞ is very large or very small. This table makes clear in particular that only r s ω ∞ 1 is consistent in the Markovian regime, since otherwise the bounds ω ∞ C S /C S 1 and ω ∞ D S /C S 1 necessarily fail (because the e iω∞t oscillations are too rapid).
The asymptotic forms of Table 2 say more than just this, however. From them we also notice that r s ω ∞ 1 implies 10 D S /ω ∼ D S and so 20) which implies that the g 2 D S /ω ∞ term appearing in the solution (4.18) is negligible relative to the g 2 C S /ω ∞ term. This means that the g 2 D S terms in the Markovian 10 Using ψ (0) (z) 1 z −γ+ π 2 z 6 −ζ(3)z 2 +. . . for |z| 1 (with ζ the Riemann zeta function) [106], D S ω∞ 2π 2 (log( 2rs ) + 4ζ (3) evolution equations (4.12) can be neglected, allowing the Markovian evolution instead to be written as (4.11) and ∂ I 12 (t) ∂t −g 2 C S I 12 (t) + g 2 C S e +2iω∞t I * 12 (t) . (4.21) In particular the divergent quantity D S plays no role in the Markovian limit, apart from shifting the qubit energy levels in the way that is renormalized into the definition of ω. Following the steps discussed at great length in [16,17]) shows that these equations preserve positivity of (t) to O(g 2 ) in the Markovian limit, with no additional approximations necessary. The solutions in the Markovian regime therefore become and I 12 (t) = e −t/ξ 12 (0) + i * 12 (0) where ξ := 1 and the last line follows since the Markovian approximation demands ω ∞ r s 1. These solutions describe the exponential decay towards a thermal distribution (with temperature T = 1/(4πr s ) = T H that equals the Hawking temperature), doing so with the characteristic time-scale ξ 8π 2 r s /g 2 . Notice that the approach to equilibrium takes place twice as fast for the diagonal components of compared to its off-diagonal parts.
We remark in passing that it is also possible to solve the Nakajima-Zwanzig equation at late times using weaker assumptions than those that lead to the above Markovian solutions, using methods similar to those used in [17] (see also [105], where a non-Markovian solution for an accelerated qubit is derived by method of Laplace transforms). The utility of such a solution is less interesting here since the Markovian condition r s ω ∞ = r s ω 1 − r s /r 0 1 is satisfied for any qubit of fixed rest-frame energy splitting that hovers sufficiently close to the black-hole horizon.

Conclusions
In summary, this paper shows how Open EFT methods can lend themselves to late-time resummation in more general gravitational systems than the cosmological examples previously explored. As in the examples of [16,17] simplicity arises near the horizon at late times, even when the underlying geometry tends to makes quantum mechanical calculations difficult. Standard tools for open quantum systems give relatively easy access to times of the order r s /g 2 , at least in the specific instance of an Unruh-DeWitt detector placed very close to a Schwarzschild horizon and interacting with a quantum field. The resulting evolution describes qubit thermalization with the expected Hawking radiation, asymptoting to the Hawking temperature T H = (4πr s ) −1 . The time-scale for thermalizing a hovering qubit can be computed, and in the very-near-horizon limit takes a universal form that relies only on properties of the near-horizon geometry given only the relatively weak assumption that the quantum field is prepared in a vacuum state of Hadamard form (including in particular the Hartle-Hawking and Unruh states).
What makes the late-time evolution easy to resum is its Markovian nature over Schwarzschild times that are long compared with r s . Autocorrelations of the field in a Hadamard state then fall off very robustly for qubits hovering very near the horizon, effectively washing out the past entanglement history. As one might expect from the equivalence principle, the qubit behaviour becomes equivalent to that of a qubit accelerating through flat space in the limit of infinite acceleration. It is the large acceleration (and blueshift) experienced by the qubit which ensures that the quantum field mass eventually becomes negligible in the near-horizon limit, explaining why the mass largely drops out of our result. As a consequence, the Markovian evolution seems likely to be very robust, at least asymptotically close to the horizon (provided r s ω ∞ 1, so that the qubit states are not too split to allow thermal excitation).
The absence of mass dependence (in the m r s 1 limit) also carries information about dependence on the scalar self-coupling, λ. Scalar self-couplings are known to give rise to secular effects for accelerated observers even in flat space [15] (see also [107] for other evidence for secular growth in black-hole geometries), where it is known that they can also be resummed at late times. For the Rindler problem late-time resummation amounts to re-organizing perturbation theory using a small shifted mass, δm 2 ∼ λa 2 , similar to the development of small temperature-dependent masses in thermal environments [15]. Similarity with the Rindler problem makes it is very plausible that a similar resummation can be obtained near the Schwarzschild horizon by shifting the scalar mass by an amount δm 2 ∼ λ/r 2 s , making the m-independence of near-horizon qubit evolution likely also to imply the same for λ-dependence, at least when λ is small and times are late.

Acknowledgements
We thank Viacheslav Emelyanov and Laszlo Zalavari for helpful discussions. This work was partially supported by funds from the Natural Sciences and Engineering Research Council (NSERC) of Canada. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Colleges and Universities.

A The Synge world function
This appendix derives some of the features of the Synge world-function that are used in the main text.

Definitions
To this end consider two points, x and x , that are connected by a timelike geodesic Γ. If λ is an affine parameterization of Γ then it is described by the curve y µ (λ) along whichÿ µ + Γ µ νσẏ νẏσ = 0 (A.1) is obeyed for all λ, withẏ µ := dy µ /dλ. The fact that Γ connects x and x is expressed as the boundary conditions y(λ i ) = x and y(λ f ) = x. For such a geodesic the Synge world function, σ(x, x ), is defined by [87][88][89] σ(x, x ) : where the integral is performed along the geodesic Γ. This integral is actually quite easy to evaluate because the geodesic equation (A.1) implies that the quantity  Performing the same calculation using Schwarzschild coordinates instead gives Although this agrees with (A.12) for ∆t r s , the domain of validity of (A.9) turns out to be larger, applying even when ∆t/r s is not small. This can be seen numerically in σ, as is shown in Figures 1 and 2. Also shown in these figures is how the domain of validity of (A.12) can be extended out to extremely large values of ∆t/r s simply by choosing r 0 to be ever-closer to the horizon itself.