Perturbative calculations of entanglement entropy

This paper is an attempt to extend the recent understanding of the Page curve for evaporating black holes to more general systems coupled to a heat bath. Although calculating the von Neumann entropy by the replica trick is usually a challenge, we have identified two solvable cases. For the initial section of the Page curve, we sum up the perturbation series in the system-bath coupling κ; the most interesting contribution is of order 2s, where s is the number of replicas. For the saturated regime, we consider the effect of an external impulse on the entropy at a later time and relate it to OTOCs. A significant simplification occurs in the maximal chaos case such that the effect may be interpreted in terms of an intermediate object, analogous to the branching surface of a replica wormhole.


Introduction
While not an observable quantity, entropy is useful as an abstract measure of active degrees of freedom and correlations in the system. Quantum correlations can be elusive, particularly in black holes, where the classical space-time picture is incomplete. There has been a long but ultimately successful chase of correlations in the Hawking radiation. If a black hole forms from an object in a pure quantum state and then evaporates, the resulting radiation must also be in a pure state. Thus, it is strongly (albeit nonlocally) correlated. The general form of such correlations was predicted by Don Page [1], who considered a black hole as a generic quantum system. Still, it long remained unclear how such correlations could emerge in semiclassical gravity. Some important works that contributed to the solution include the Dray-'t Hooft mechanism of gravitational interaction between infalling matter and subsequent radiation [2,3] and the calculation of out-of-time-order correlators (OTOCs) in the black hole setting [4]. However, the OTOC physics is relevant on short time scales and explains correlations that are present not in the radiation itself but relative to a purifying system [5] (that is, under the assumption that the black hole is part of a thermofield double and that we have unrestricted access to the other part). The recent breakthrough in understanding the correlations developing over the Page time [6,7] required a careful formulation of the problem, which we will now summarize.
The problem at hand is a semiclassical one. We do not have a complete theory of quantum gravity, nor should it be required. When working at the semiclassical level, it is not possible to derive long-term evolution from the short-term one. Rather, one should look for a global solution, which may depend on the quantity of interest. We consider the JHEP03(2021)198 entanglement entropy between the black hole and the emitted radiation at a particular time t. So let ρ = ρ(t) be the black hole's density matrix; we want to compute its von Neumann entropy, S(ρ) = − Tr(ρ ln ρ). The latter is expressed as the s → 1 limit of the s-Renyi entropy, For integer s, the expression Tr ρ s may be interpreted as the partition functions of s replicas of the system. Now, it turns out that the transition from the early phase of the black hole evaporation (when the radiation is uncorrelated as the naive theory predicts) to the later phase (when the entanglement entropy equals the black hole's coarse-grained entropy) is first order. The later phase is described by a new type of space-time geometry, the replica wormhole [6,7]. Although choosing the correct solution of the two is a global problem, each of them can be examined locally. We will study some properties of both solutions for general many-body systems, where the geometric description is not applicable.
The von Neumann and Renyi entropies are nonlinear functions of the quantum state, which is why they are not observables. However, the logarithmic nonlinearity is mild, such that in the thermodynamic limit, S(ρ) is determined by typical microstates that contribute to the mixed state ρ. In contrast, Renyi entropies are often dominated by a fraction of microstates of tiny overall weight. This distinction is also evident from the replica wormhole picture. The s-Renyi entropy is related to an s-fold cover of space-time, whose metric is different from the physical one. But when we analytically continue the solution in s and take s to 1, we get the standard metric with an additional piece of data, the branching surface. Thus, the s → 1 limit is essential for compatibility with the usual (non-entropic) physics. Our main technical advance is how to take this limit in some specific cases.

Early phase of entanglement growth
We adopt a simpler variant of the evaporation problem, where instead of radiating energy, the system comes into contact with a heat bath at the same temperature. Turning the system-bath interaction on represents a slight change in the Hamiltonian and results in a brief period of non-equilibrium dynamics. Then a steady state is achieved such that all simple correlation functions are thermal. However, if the system's initial state was pure (though mimicking the thermal state), its von Neumann entropy will grow at a constant rate. We focus on this regime as well as the very beginning of quantum evolution. The entropy growth eventually saturates at the thermal (i.e. coarse-grained) entropy, but that is not captured by our method.
Our calculation is perturbative in the system-bath coupling strength κ. Note that the von Neumann entropy has a logarithmic singularity at the unperturbed state, which is pure. This is reflected by the fact that in addition to terms of order κ 2 (or any constant power of κ), terms of order κ 2s (where s is the number of replicas) play an important role.

The model and general formulas
Let us consider a quantum system (meant to represent a black hole) with some Hilbert space H B and Hamiltonian H B . For an exact analogy with the evaporation problem, we would have to pick a pure state that looks like thermal to all simple measurements. Instead, we double the system and postulate that its initial state is the thermofield double, Only the right part is coupled to the heat bath, but we are interested in the von Neumann entropy of the double system as its density matrix ρ B * B evolves in time. Likewise, the bath is also doubled, so that the initial state of the world is A similar, but not identical, 1 setting was used in [8,9], where the s-Renyi entropy for integer s > 1 was calculated. The full Hamiltonian H = H B + H b + H Bb acts only on the two objects in the middle, i.e. on H B ⊗ H b . We assume that the interaction term has the form In the case of fermionic systems like the SYK model [10][11][12][13], we should multiply the coupling parameter κ by i. For simplicity, we will do the computation for a bosonic system, but the final answer will equally be applicable to fermionic systems.
Thus, the evolution of the world in the interaction picture takes the form where T stands for time ordering. We also assume that O j B = O j b = 0. 2 In the rest of the section, we will compute the s-Renyi entropy of the system's density matrix ρ B * B (t) after tracing out the bath. It is given by the perturbative expansion where the expectation value · · · is with respect to the bath's thermal state and T denotes reverse time ordering. (If u 1 < · · · < u n and u 1 < · · · < u m , then the operators are already ordered.) There is also an implicit sum over repeated indices, with each index going from 1 to N . Note that operators with same indices have the same time argument. 1 In refs. [8,9], the initial state is taken to be the thermofield double of two interacting subsystems rather than the product of two thermofield doubles. 2 In case that O j = 0, one can work with O j − O j .

JHEP03(2021)198
Since the combinatorics might soon get complicated, let us introduce some simplifying graphic notation: (2.5) Then each term in the expansion (2.4) will look like this (where T, T, and the indices are omitted): (2.6) The diagram element in the middle is the Keldysh contour for the heat bath. It consists of a circle at the bottom representing imaginary-time evolution and a stem corresponding to the real-time evolution; the time goes up. For integer s, Tr (ρ B * B (t)) s can be represented by gluing s such diagrams (describing different replicas of the density matrix) in the cyclic order -see figure 1, where the replicas are depicted with different colors. The expectation values should be independently computed for each closed contour, whether it corresponds to the system or the bath.
Let us further assume that the system-bath coupling is sufficiently weak. Then the "radiation quanta" emanating from the system are sparse, which means that dominant diagrams have at most two operators with close times per contour. Therefore, the calculation can be done using Wick contraction. Of course, if the fields O j B , O j b are Gaussian, then no sparseness condition is necessary.
An example of a (subleading) Wick pairing contributing to Tr (ρ B * B (t)) s is as follows: It corresponds to the black loop at the bottom of figure 1 (a). In general, a Wick contraction diagram is a disjoint union of loops that consist of alternating solid and dotted lines. Solid lines represent contractions of fields on the same contour, whereas dotted lines correspond to interaction terms such as O j . Each contraction comes with a Kronecker delta identifying the indices of the contracted fields. The result is nonzero if all the indices on each loop are the same. Therefore, each loop with d solid line segments evaluates to N κ d multiplied by the product of two-point functions.
The diagrams for ln Tr (ρ B * B (t)) s are connected, i.e. contain a single loop. The factor of N that appears here is usual for extensive thermodynamic quantities; thus, the intensive parameter is κ d . If κ is small, we should only keep loops with d = 2 contractions (in one replica) and d = 2s contractions (traversing all the replicas). Both types of loops are shown in purple in figure 1.
There are three types of two-point functions for the system, and similarly for the bath. The expressions for closed paths will be simplified if we think of a two-point function as the matrix element of a bilocal operator G, where f t (u) is a time window function that vanishes outside the interval (0, t). This way, we can extend the time domain to (−∞, ∞) and avoid putting limits on the integrals. We also define the transpose of the operator G, denoted by G , with the matrix element t 1 | G |t 2 = G(t 2 − t 1 ). Note that G T and G T for bosons are symmetric, i.e. equal to their transpose. To illustrate this notation, the expression (2.7) involves a loop with four Wick contractions, As already mentioned, there are two types of loops that contribute to ln Tr (ρ B * B (t)) s to leading order. The loops of length d = 2 can themselves be of two forms, one of which appears in figure 1. They give the following contributions: Here we have used the fact that (±i) 2 from (2.4) cancels i 2 from (2.8). (The extra factors present in the fermionic case also cancel each other.) There are also loops of length d = 2s, which we say to have winding number 1 because they traverse all replicas. The expression for such a loop takes the following form: Since in the end we are interested in the limit s → 1, the quantities (2.12) and (2.11) are of the same order in the coupling.

JHEP03(2021)198
As an exercise, let us sum up the leading diagrams in Tr (ρ B * B (t)) s -we should get the exponential of a sum of single loops with certain coefficients. It is sufficient to only keep track of the system's replica, leaving the bath implicit as in figure 1 (b). Let m 1 , . . . , m s and n 1 , . . . , n s be the numbers of fields in the time ordered and anti-time ordered branches of the Keldysh contours. The number of ways to break them into k loops of length 2s with winding number 1 and some loops of length 2 (with winding number 0) is given by Defining m r − k = 2p r and n r − k = 2q r , after manipulation we will get Using the fact that , the final answer is as follows: (2.15) While the above equation was derived for bosonic systems, it is the same for fermionic systems like the SYK model.

Short initial period vs. linear growth
Equation (2.15) allows one to compute the von Neumann entropy S(ρ B * B (t)) for t < t Page .
Although the exact answer is model-dependent, there are two universal regimes: very early times, just after the system-bath coupling is turned on, and intermediate times, when the entropy grows linearly.
Very early times. Let t t UV such that the effect of the Hamiltonians H B and H b is negligible. For example, t UV = J −1 for the SYK model. More exactly, we assume that the Green functions G B (t) and G b (t) may be approximated by some constants. Then the expression (2.15) and the von Neumann entropy take this form: Intermediate times. For systems with continuous excitation spectrum, connected correlators decay in time. Exponential decay is typical; for example, in a conformal system at finite temperature, the correlator of fields with scaling dimension ∆ decays as exp − 2π∆ β t if t β ∆ . Let us assume that both G B (t) and G b (t) decay exponentially at t t * .

JHEP03(2021)198
If t t * , then Tr G B • G b can be approximated as follows. This expression is as integral over u, u ∈ (0, t), but the integrand is negligible unless |u − u| t * . Therefore, we may remove the limits on u , and then use the Fourier transform: The second term in (2.15) can also be approximated in such a way. Thus, where (2.20) After analytically continuing to s = 1, the von Neumann entropy will be given by The integrals in A(s) and A (1) converge because a possible peak at ω = 0 is broadened to have a width ω min = t −1 * . There is also a natural UV cutoff at ω max = t −1 UV . An interesting case is where both the system and the bath are conformal at ω ω max . Let us first assume that the temperature is zero; . (2.23) A good example is two SYK models at large βJ. (A bath with ∆ b = 1 2 can also be realized by a critical Majorana chain.) The Renyi entropies in this case were studied in [9] using the effective action method, which is generally more powerful than perturbation theory. However, the analytic continuation to s = 1 was not obtained. The computed growth rate of the s-Renyi entropy for integer s > 1 is consistent with our estimate,

Perturbations to the saturated phase
We now consider the system at later times, such that its von Neumann entropy has reached the coarse-grained (thermodynamic) entropy. The entanglement entropy in this phase shows interesting behavior under perturbations. For example, a short impulse increasing the system's energy (similar to throwing a rock into a black hole) will cause a resurgence JHEP03(2021)198 of entropy growth. Indeed, such an action can be described by some unitary operator V . It increases the coarse-grained entropy and effective temperature, though the true microscopic entropy does not change. Letting the system interact with the bath, we should see a behavior similar to the cusp in the Page curve. Specifically, we expect the von Neumann entropy to grow until it becomes equal to the coarse-grained entropy. Since the growth rate is constant while the perturbation can be arbitrarily weak, the resurgence can be short -just slightly longer than the scrambling time. It will be followed by a thermal equilibration period, when both the coarse-grained and microscopic entropies decrease, see figure 2. Thus, the Page curve cusp is accessible in this setting. However, to actually produce a cusp, the perturbation should be sufficiently strong, likely beyond the Taylor expansion. We will study a simpler problem, calculating the effect in the lowest order. While our formal goal is to compute the von Neumann entropy in a rather general setting, the key result pertains to systems that saturate the chaos bound [14]. The expression obtained in this case admits a holographic interpretation, which will be discussed section 4.

Statement of the problem
For the study of the saturated phase, it is sufficient to consider one copy of the system and the bath rather than the thermofield double. The bath can be integrated out, giving rise to the interaction function σ(τ 1 , τ 2 ) = κ 2 G b (τ 1 , τ 2 ) on the Keldysh contour. Let us simplify the model a bit and replace the interaction σ, which is constantly on, with a superoperator R acting at a specific time. (The results we will obtain in this setting can be easily generalized to the original model.) So, the exact problem involves a quantum system at thermal equilibrium subjected to a sequence of two instantaneous perturbations, V and R.
Let V = e −ixX , where X is a Hermitian operator and x is a small parameter. We consider the action of V on the thermal state ρ 0 and expand the resulting density matrix ρ 1 (x) to the second order in x: In fact, a nontrivial effect will be seen in the second order, and only when combined with a subsequent interaction with the environment. The latter is described by a physically realizable (i.e. completely positive, trace-preserving) superoperator R. Suppose that R is

JHEP03(2021)198
close to the identity such that it can be expanded to the first order in some parameter : where A · B stands for the superoperator that takes ρ to AρB. The first term in L (which involves a Hermitian operator C and act as ρ → −i[C, ρ]) may be neglected because it represents an infinitesimal unitary transformation, and thus, does not change the entropy. The sum over j (known as Lindbladian) corresponds to tracing out the environment. As will be justified later, we may replace L with j A j ·A † j so that the final density matrix becomes Our goal is to compute ∂ 2 ∂x 2 ∂ ∂ S(ρ(x, )), where S(ρ) = − Tr(ρ ln ρ). We assume that V acts at time 0, whereas R acts at a later time t. Thus, A j is understood as A j (t) = e iH 0 t A j (0)e −iH 0 t , where A j (0) is some simple (e.g. one-or twobody) operator. The calculation will be done by the replica method for a general large N system in the early time regime, i.e. before the scrambling time. However t is taken to be sufficiently large such that OTOCs are parametrically greater than correlators with non-alternating times. Note that ρ(x, ) involves only non-alternating operators such as A j Xρ 0 XA † j . However, OTOCs appear due to the use of replicas. The "unimportant terms" in (3.3) are exactly those that do not generate any OTOCs.
In the next section, we study partial derivatives of S(ρ), assuming that ρ depends on parameters in some particular way. This setting does not directly include the function ρ(x, ) given by equation (3.3). To cover this case, we will use a trick called "locking two operators in the same replica", see section 3.5.

Thermodynamic response theory for the replicated system
Let us recall the standard definition of connected correlators. We begin with the partition function Z = Tr W , where W is the imaginary-time evolution operator: Without perturbation, we have H(τ ) = H 0 . The insertion of operators X 1 , . . . , X n at times τ 1 , . . . , τ n is described by perturbing the Hamiltonian: where x j are infinitesimal numbers. We generally assume that the operators X j are bosonic.
(If any of them is fermionic, the corresponding variable x j should be anti-commuting.) Thus,

JHEP03(2021)198
and Z(β, x n , . . . x 1 ) = Tr W (β, x n , . . . x 1 ). The full correlator is simply The corresponding connected correlator is defined as follows: 3 Now, let us introduce s replicas of the system, such that the partition function becomes We may think of the parameter s as being associated with a branching operator B, which commutes with everything. It is not defined by itself but only through its connected correlators: The branched correlator (3.11) is related to the entropy S = S(ρ) of the density matrix Thus, the entropy derivative with respect to x 1 , . . . , x n is given by − B, X n (τ n ), · · · , X 1 (τ 1 ) + X n (τ n ), · · · , X 1 (τ 1 ) . It is usually the easiest to compute the derivative of the relative entropy, S(ρ||ρ 0 ) = Tr(ρ(ln ρ − ln ρ 0 )): ∂ n S(ρ||ρ 0 ) ∂x 1 · · · ∂x n x 1 =···=xn=0 = (B +βH 0 ), X n (τ n ), · · · , X 1 (τ 1 ) − X n (τ n ), · · · , X 1 (τ 1 ) (3.13) For integer s, the expression Tr (W (β/s, . . .)) s may be interpreted in terms of gluing s intervals of length β/s to make a circle of length β. The operators X n (τ n ), . . . , X 1 (τ 1 ) are distributed along that circle. Thus, the number in question is, essentially, a correlation function at the given β. An important caveat is that there is no natural definition of the full correlator B A as a function of A such that one could compute B Y Z by substituting Y Z for A. If such a function (with the usual relation to the connected correlator) existed, we would have this corollary of equation (

Branched two-point correlator
Suppose the ordinary correlation function Y (τ ), X(0) is known on the imaginary axis, τ = it, and let us use its Fourier transform in t. In these terms, (3.14) The corresponding branched correlator is expected to have a similar form, The goal of this section is to find the function h Y,X . Let us consider the Fourier modes of the operators Y and X, for example, Y ω = Y (it)e iωt dt. Their connected correlator is and we also have We now calculate the branched correlator of Y ω and X ω , equal to h Y,X (ω) · 2πδ(ω + ω ).
When the number of replicas s is a positive integer, each of the operators in question can be inserted in any replica, so the calculation involves a double sum. Since each replica's length is β/s, putting Y ω in the k-th replica is described by Y ω (kβ/s) = Y ω e −kβω/s . With this in mind, we get:

Branched correlator related to early-time OTOCs
Let us recall the original problem of computing S(ρ(x, )) with ρ(x, ) given by equation (3.3).
In this section, we calculate an analogous branched correlator B, A j (β + it), X(β), X(0), A † k (it) and, more generally,  One can eliminate β from the time arguments by cyclically permuting X 4 , . . . , X 1 . As already mentioned, the replica calculation involves OTOCs, which are dominant for sufficiently large t. Neglecting all terms with non-alternating times, we get: (Equation (3.25) is illustrated by figure 3.) In order to make further progress, we will use the single-mode ansatz for early-time OTOCs [13], combined with the Fourier representation The result has this general form:  (3.33) and hence, (3.34) The function f − is obtained from f + by replacing w with w −1 . Adding both terms together, we get: ln u ln v where u = e −βω 23 , v = e −βω 14 , and w = e −iβκ/2 . A great simplification occurs in the maximal chaos case: f (ω 14 , ω 23 ) = 2π (1 + e −βω 23 )(1 + e βω 14 ) if κ = 2π β . (3.36) Importantly, the function f (ω 14 , ω 23 ) splits into two factors. They may be interpreted in terms of interaction of the fluctuating horizon (which corresponds to B) with incoming and outgoing radiation, see section 4.

Locking two operators in the same replica
We now adapt the obtained result to express the entropy of the density matrix ρ(x, ). The latter is a normalized version of the operator Note that we have made the time explicit and will follow the convention that A j = A j (0). The corresponding entropy could be calculated as in the previous section while restricting JHEP03(2021)198 X 1 = A † j and X 4 = A j to be in the same replica. However, instead of imposing this restriction by hand, we will modify the problem so that it fits the branched correlator setting. First, we replace the set of operators A j with a single operator Y . This is achieved by extending the physical system with an auxiliary one, comprising a ground state |0 with zero energy and a set of excited stated state |j with energy Ω. We denote the Hamiltonian of the extended system by H(Ω) and set Although the transformation just described alters the operator W (β, x, ) in a nontrivial way, we will find an agreement in the Ω → ∞ limit. For the time being, let us construct some operators acting on the extended system that correspond to the two terms in (3.37) as closely as possible:

JHEP03(2021)198
(B is not an independently defined object but rather, part of the notation.) While the calculation of branched correlators is rather involved, a simplification occurs for maximally chaotic systems. We now argue that this special case is consistent with a holographic picture in a very broad sense. For comparison, consider a Euclidean black hole in a hyperbolic space (say, in two dimensions). The replica geometry involves an s-fold cover of both the circle and the disk it bounds, with a branching point at the center. Inserting boundary fields slightly deforms the space. In the s → 1 limit, the geometry is given by a smooth metric on the disk and the position of the branching point. The branching point is a special case of a quantum extremal surface [15] (where "surface" means a codimension 2 submanifold). Its position is determined by an extremum of entropy. Instead of the entropy S, we may consider ln Z −S. Indeed, the partition function Z depends only on the space-time metric, which should be fixed before finding the extremal surface.
In the Lorentzian case, the branching point is described by null coordinates (u + , u − ). (We set aside the ambiguity in the choice of origin due to the deformation of space-time relative to AdS 2 .) The entropy can be expanded to the second order in u + , u − : where p + and p − depend on the inserted field. 5 Solving the extremum problem, we get Now, let us forget about geometry. The only property we need is that if there is large time separation, then p + and p − depend only on the fields inserted in the past and the future, respectively. Thus, the change in the entropy should factor into two quantities dependent on the corresponding fields. This is exactly what we observed in the maximal chaos case, see equation (3.36). We leave the interpretation of these quantities to future research.