Replica Wormholes and Holographic Entanglement Negativity

Recent work has shown how to understand the Page curve of an evaporating black hole from replica wormholes. However, more detailed information about the structure of its quantum state is needed to fully understand the dynamics of black hole evaporation. Here we study entanglement negativity, an important measure of quantum entanglement in mixed states, in a couple of toy models of evaporating black holes. We find four phases dominated by different types of geometries: the disconnected, cyclically connected, anti-cyclically connected, and pairwise connected geometries. The last of these geometries are new replica wormholes that break the replica symmetry spontaneously. We also analyzed the transitions between these four phases by summing more generic replica geometries using a Schwinger-Dyson equation. In particular, we find enhanced corrections to various negativity measures near the transition between the cyclic and pairwise phase.


Introduction
Replica wormholes have played an important role in recent progress on solving the black hole information problem [1,2]. These wormholes appear as nontrivial saddle points that could dominate gravitational path integrals with replicated boundary conditions. Their appearance leads to nontrivial "island" contributions in the quantum extremal surface (QES) formula for gravitational entropy [3][4][5][6].
So far most of the discussion has been centered on the von Neumann entropy. While obtaining the von Neumann entropy is a good first step, we need more detailed information about the quantum state -such as more general measures of entanglement -to fully solve the black hole information problem.
In this paper, we take a first step towards understanding the structure of entanglement in an evaporating black hole and its Hawking radiation by studying entanglement negativity and its Rényi generalizations in a couple of toy models. Just as the von Neumann entropy is a measure of quantum entanglement in pure states, the negativity is an important measure of entanglement in generally mixed states. Therefore, the negativity provides an interesting probe in diagnosing the structure of multipartite entanglement in systems such as an evaporating black hole.
To understand negativity intuitively, consider a general state on two subsystems that is described by a density matrix, and take its partial transpose on the second subsystem. The partially transposed density matrix could have negative eigenvalues, and the degree to which the eigenvalues are negative is characterized by the negativity and logarithmic negativity. Both of these negativity measures are entanglement monotones, and the logarithmic negativity provides an upper bound on the distillable entanglement [7][8][9]. Negativity has been discussed in a number of interesting prior works .
We now give a short summary of this paper. In Section 2, we review the definition and properties of the negativity and its Rényi generalizations. In Section 3, we start our study of negativity in a toy model of an evaporating black hole in Jackiw-Teitelboim (JT) gravity with an end-of-the-world (EOW) brane. This is a slight generalization of the model studied in [1], with the system describing the Hawking radiation divided into two subsystems so as to study negativity.
As we tune the parameters of the model, we find a rich phase diagram for the negativities consisting of four phases (see Figure 4). Each of the four phases is dominated by a saddle-point geometry of JT gravity (or a set of saddle points). For a black hole before the Page time, we find a phase dominated by a totally disconnected geometry, whereas after the Page time, we find three distinct phases depending on how we divide the radiation system into two subsystems: the first phase is dominated by a cyclically connected geometry (which is the replica wormhole of [1]), the second dominated by an "anti-cyclically" connected geometry, and the third dominated by pairwise connected geometries that are in one-to-one correspondence with non-crossing pairings. These pairwise connected geometries are new replica wormholes that spontaneously break the replica symmetry. Their appearance agrees with the general discussions on holographic negativity in [30].
In Sections 4 and 5, we study the behavior of negativities near the transitions between the four phases. Near these phase transitions, more geometries than the four types described earlier could dominate the gravitational path integral for Rényi negativities, and we need to sum over them. In order to obtain the negativity and logarithmic negativity (as well as related negativity measures such as the partially transposed entropy [20,30]), we need to analytically continue in the replica number. We achieve this by using the resolvent for the partially transposed density matrix to find its eigenvalue distribution (which we call the "negativity spectrum"). To calculate this "negativity resolvent," we organize the sum over geometries into a Schwinger-Dyson equation, which is similar to the method used in [1]. We develop this method for negativity in Section 4 and apply it to both a microcanonical ensemble and canonical ensemble in Section 5.
When the black hole is in a microcanonical ensemble, the Schwinger-Dyson equation simplifies into a cubic equation for the negativity resolvent, leading to concrete results for the negativities near all phase transitions. This is similar to the case of a random mixed state studied in [29].
When the black hole is in a canonical ensemble, the gravitational calculation is technically more difficult. As a result, we study each of the phase transitions separately, for we only need to sum over a subset of geometries near each transition. Near the transition between the "disconnected" phase and "pairwise" phase, the Schwinger-Dyson equation again simplifies, this time into a quadratic equation for the negativity resolvent. Near the transition between the "cyclic" phase and "pairwise" phase, it is difficult to solve the Schwinger-Dyson equation exactly, but we solve it approximately in the semiclassical, or β → 0, limit. From its solution, we find that the negativity spectrum near the phase transition consists of two branches, each of which is approximately a shifted thermal spectrum with a cutoff. One branch consists of positive eigenvalues, and the other has negative eigenvalues. From this we find enhanced corrections to various negativity measures near the phase transition. In particular, a quantity known as the refined Rényi-2 negativity receives an O(1/ √ β) correction, similar to the enhanced corrections to the von Neumann entropy at the Page transition [1,33,34], whereas other negativity measures such as the logarithmic negativity and the partially transposed entropy exhibit O(1/β) corrections, similar to what happens to Rényi entropies S n with n < 1.
Moving beyond the JT gravity model, we study in Section 6 the behavior of negativities in a topological model of 2-dimensional gravity with EOW branes. This is a slight generalization of the model of [35], where we again divide the radiation system into two subsystems to study negativity. We find the situation to be very similar to the case of a microcanonical ensemble in JT gravity described earlier. In particular, the Schwinger-Dyson equation again simplifies into a cubic equation for the negativity resolvent, leading to concrete results for the negativities.
We end with some concluding remarks in Section 7 and several appendices. In Appendix A, we derive the set of dominant geometries in each of the phases and near phase transitions. In Appendix B, we provide a more detailed analysis near the transition between the cyclic phase and pairwise phase in the canonical ensemble. In Appendix C, we study the Rényi entropies near the Page transition in a similar fashion and show that they exhibit corrections analogous to the corrections to the negativities near the phase transition.
Related works appeared recently and have some partial overlap with our results on the study of negativity in JT gravity in the microcanonical ensemble [36] and the canonical ensemble [32,37].

Entanglement negativity and Rényi negativities
The motivation for entanglement negativity comes from the Peres-Horodecki criterion [38,39], also known as the PPT (positive partial transpose) criterion for mixed states, which we review here. Consider a mixed state ρ AB defined on the product Hilbert space for states ρ i A and ρ i B on H A and H B , respectively. Separable states are classical mixtures of product states and thus do not contain quantum entanglement; inseparable states are said to be entangled.
We denote the algebra of operators on H i by A i , and the space of linear maps from maps positive operators to positive operators, and is completely positive if for all nonnegative integer n, is positive, where M n denotes the algebra of n × n complex matrices. For separable states, this condition is clearly satisfied when Λ is a positive map, as ( For inseparable states, this no longer holds in general, so a good diagnostic of entanglement would be a positive but not completely positive map, such that entangled states would have negative eigenvalues under the action of Λ ⊗ I.
The partial transpose is such a positive but not completely positive map. Consider a bipartite system AB with an orthonormal basis {|a } on A and {|b } on B. Given a density matrix ρ AB on AB, we define the partially transposed density matrix as 1 Acting on a reduced density matrix on B, the partial tranpose becomes the usual transpose which preserves the eigenvalues of the original reduced density matrix and is therefore a positive map. Acting on the full density matrix is not guaranteed to preserve positivity. As an example, take an EPR pair of two qubits A and B. The partial transpose of its density matrix has eigenvalues { 1 2 , 1 2 , 1 2 , − 1 2 }. We therefore see that the partial transpose can be a useful tool for differentiating between separable and inseparable states. 2 The entanglement negativity N (ρ) is defined as the sum of the absolute values of the negative eigenvalues of this partially transposed density matrix and can be variably written as Here ||X|| 1 ≡ Tr |X| = Tr √ X † X is the Schatten 1-norm of a matrix X. We see why negativity is such an appealing entanglement measure, as it is computed directly from a trace, as opposed to a variational principle in the case of other entanglement measures. Note that as we are taking a trace, it does not matter which subsystem we take the partial trace over, so choosing ρ T B AB instead of ρ T A AB is merely a convention. The logarithmic negativity is similarly defined by The logarithmic negativity is an upper bound on the distillable entanglement, i.e., the asymptotic number of EPR pairs that can be extracted from a set of identically prepared ρ AB with local operations and classical communication (LOCC).
We can also write a Rényi version of negativity via There is a subtlety in the analytic continuation of the Rényi negativity. As the negativity is defined by the absolute value of the eigenvalues of the partially transposed density matrix and the Rényi negativity is defined without an absolute value, we need to define different analytic continuations for even and odd n such that The logarithmic negativity is then obtained from the even analytic condition: In this paper, we will also be interested in a generalization of the Rényi negativities termed refined Rényi negativities, which are given by (2.10) The refined Rényi negativities are inspired by the refined Rényi entropies defined in [40]. In particular, we will be interested in two measures that descend from the refined Rényi negativities. The first is the partially transposed entropy S T B (ρ AB ) of [20,30], defined as the m → 1 limit of the refined odd Rényi negativity: S T B is so named in analogy with the von Neumann entropy. The other measure is the refined Rényi-2 negativity S T B (2,even) (ρ AB ), which can be written as We will refer to this quantity as S T B (2) for short. It is equivalent to the von Neumann entropy of ρ T B

JT gravity with EOW branes
We start by reviewing the simple model of black hole evaporation studied in [1] (see also [41,42]). This model consists of a black hole in JT gravity, decorated with an end-of-the-world (EOW) brane with tension µ. The action is given by with the JT action being We have set G N = 1, though it can be restored by sending the inverse temperature β → βG N . The parameter S 0 can be thought of as the zero temperature entropy of an eternal two-dimensional black hole. The EOW brane is endowed with k orthonormal states, or "flavors," which are entangled with an auxiliary reference system R. The states on the brane can be thought of as describing the interior partners of the early Hawking radiation R, so by increasing k we can probe later regimes of an "evaporating" black hole. The entangled state of the black hole system B and the "radiation" R can be written as The density matrix of the R subsystem is therefore The inner product ψ i |ψ j B is given by a gravitational path integral with standard Dirichlet boundary conditions on an asymptotic boundary interval and Neumann boundary conditions on the EOW branes: as well as specifying the brane states i and j at the endpoints of the EOW brane.
As was shown in [1], while naively ψ i |ψ j ∝ δ ij , this should be understood as an Figure 1. Boundary conditions for ρ R 1 R 2 and ρ T 2 R 1 R 2 . Blue (dashed) lines denote states in R 1 , and red (dotted) lines denotes states in R 2 . If we take a trace, these two boundary conditions are equivalent. expectation value in an ensemble, and wormhole contributions indicate exponentially small fluctuations of the inner product. We illustrate the boundary conditions for the matrix elements of ρ R as follows: The solid black line denotes an asymptotic boundary interval for the gravitational path integral, while the blue dashed lines are index lines that impose boundary conditions for the brane states. Computing Tr (ρ R ) means contracting the open index lines and summing over all possible geometries respecting the boundary conditions (3.5).
To study negativity, we need to consider a bipartite mixed state. To that end, we split the radiation system into two subsystems H R = H R 1 ⊗ H R 2 consisting of k 1 and k 2 states, respectively, such that k = k 1 k 2 and We will refer to this partitioned density matrix as ρ R from now on. We define our partially transposed density matrix as the partial transpose over R 2 , i.e., We will use the shorthand ρ T 2 R moving forward. This partial transpose affects the boundary conditions for our path integral by swapping the brane flavor index lines corresponding to states in H R 2 . The resulting boundary conditions are illustrated in Figure 1.

Dominant saddles
As in any calculation with a gravitational path integral, our first task is to identify the saddle-point geometries which obey the given boundary conditions and sum over them with the appropriate weight. As our goal is to compute Rényi negativities Tr ρ T 2 R n , our boundary conditions will consist of n copies of the boundary conditions illustrated on the right of Figure 1, with matching brane flavor indices contracted. The set of all classical saddles consists of oriented two-dimensional surfaces which end on the asymptotic boundaries and EOW branes, possibly connecting two or more boundaries.
As our gravitational action (3.1) is independent of brane flavor, we can factorize the flavor contributions so that for some function f of the brane Hilbert space dimensions k 1 and k 2 . The gravitational partition function Z n for a surface connecting n boundaries depends on the Euler characteristic χ of the surface in the schematic form The contribution of a surface with genus g ≥ 1 is therefore suppressed by e −2gS 0 for large S 0 . This means that the only classical geometries we need to consider are disks or disjoint unions of disks, and Z grav is a product of disk partition functions Z n . We will therefore assume S 0 1 throughout the rest of the paper. We illustrate some examples of these disk geometries as well as a higher genus geometry in Figure 2.
More precisely, for a disk connecting n boundaries in JT gravity we have In order, they are the disconnected (g = 1), cyclic (g = X), anti-cyclic (g = X −1 ), and pairwise (g = τ ) geometries. The pairwise geometries spontaneously break replica symmetry. The black lines are oriented asymptotic boundaries, the purple lines are EOW branes, the blue (dashed) lines denote k 1 index loops, and the red (dotted) lines denote k 2 index loops.
where ρ(s) = s 2π 2 sinh (2πs) is the disk density of states 3 in JT gravity. In order to recover (3.10), we take S 0 1 while keeping other parameters fixed. We emphasize this is a schematic approximation that should only be used to motivate the pertinent saddles for our problem; in general, there are parametric corrections from the full expression of Z n which will be discussed in more detail in Sections 4 and 5.
Since we are ignoring higher genus surfaces, the sum over geometries with n replicated asymptotic boundaries is equivalent to a sum over elements of S n , the permutation group on n elements. In particular, the sum takes the form where χ(g) is the number of disjoint cycles of the permutation g, |c i (g)| is the length of the i-th disjoint cycle of g, and X (X −1 ) is the (anti-)cyclic permutation of length n. Unless otherwise specified, we will take k, e S 0 1. Note (3.14) 3 The density of states is more typically written in the E = s 2 energy basis such that ρ(E) = 1 4π 2 sinh 2π √ E . Here s can be thought of as an entropy, and is the more natural variable for our purposes.
The sums (3.13) and (3.12) over elements of the permutation group has a simple geometric interpretation. The permutation g determines how the asymptotic boundaries are connected by EOW branes, while the powers of k 1 and k 2 count the number of index loops. The totally disconnected geometry corresponds to g = 1, while the totally connected geometry corresponds to g = X. What does the g = X −1 geometry look like? We show examples of these three geometries in Figure 3. Two of these geometries, the disconnected and cyclic geometries, belong in the class of non-crossing diagrams discussed by [1]. The third, the anti-cyclic geometry, is in some sense equivalent to the cyclic geometry if one reverses the orientation of the boundary, or equivalently if one exchanges k 1 and k 2 . This statement will be explained in more detail in Section 4.2.
One might naively guess that the anti-cyclic geometry would never dominate the Rényi negativity, for the same reasons as in [1] where crossing partitions in the calculation of the Rényi entropy were suppressed by factors of 1/k 2 . In fact this geometry dominates in a very large parameter regime: as we show in Appendix A, we have the following phases dominated by the corresponding permutation g: Totally disconnected: Cyclically connected: Pairwise connected: To see this intuitively, we calculate the contributions to Tr ρ T 2 R n from these geometries. Consider the totally disconnected phase dominated by g = 1. We have χ(g) = n and χ(g −1 X) = χ(g −1 X −1 ) = 1, so this diagram contributes schematically to the sum in (3.13) as in [1]. It is then unsurprising that this dominates in the parameter regime e S 0 k, since it is the unique diagram which maximizes the power of e S 0 . Note that the contribution only depends on k, and not k 1 and k 2 individually. Now, consider the cyclically connected phase with g = X. Then χ(g) = 1, χ(g −1 X) = χ(1) = n, and 1, n odd, 2, n even. (3.17) Hence, the cyclic diagrams contribute schematically g = X ⇒ e S 0 k n 1 k f (n) 2 (3.18) to the sum in (3.13). This configuration maximizes the power of k 1 , and therefore it is expected to become important in the parameter regimes where k 1 is comparably large.
In fact, as we prove in Appendix A, it is the unique dominant diagram in the regime k 1 k 2 e S 0 . Note that compared to the Rényi entropy calculation in [1] (which can be recovered by setting k 2 = 1 here), the cyclic geometry is suppressed by 1/k . We will also show this diagrammatically in Section 4.
The anti-cyclically connected phase with g = X −1 is similar, so we will not go through the analysis. In the end, the anti-cyclic diagrams contribute schematically to the sum (3.13). It is thus expected to become important in the parameter regimes where k 2 is comparably large, and we can prove that they are the unique dominant diagrams when k 2 k 2 e S 0 . Finally, there is one additional class of dominant geometries we should consider: the pairwise connected phase with g = τ . As we show in Appendix A, these diagrams dominate in a fourth regime satisfying both k 1 k 2 e S 0 and e −S 0 k 1 /k 2 e S 0 , and they are the only diagrams aside from the disconnected, cyclically connected, and anti-cyclically connected diagrams that can dominate in a large regime of the parameter space. These geometries are in one-to-one correspondence with the set of permutations τ known as non-crossing pairings. For even n, a pairwise connected geometry is constructed by choosing an element in τ , for example (12)(34) · · · (n−1, n), and connecting paired asymptotic boundaries by two-boundary wormholes. For odd n, the geometries are given by a similar non-crossing pairings of boundaries, plus a single one-boundary connected component. We show an example of such a geometry in Figure 3. It is evident that such geometries spontaneously break the replica symmetry.
As we show in Appendix A, each pairwise connected geometry contributes schematically to the sum in (3.13), where n 2 and n 2 denote the ceiling and floor function, respectively. A pairwise connected diagram in some sense puts k 1 , k 2 , and e S 0 on the most equal footing by maximizing the sum of the three exponents in (3.13). As with the disconnected geometry, the contribution of a pairwise connected geometry only depends on k, and not k 1 and k 2 individually.

Contributions to negativities
Having worked out the dominant geometries in the four phases and how they contribute to the Rényi negativities schematically (i.e., using Z n ∼ e S 0 ), we now write their

Negativities in Dominant Phases
2m m is the Catalan number which gives the number of non-crossing pairings.
contributions exactly using (3.12). This calculation is straightforward to do, for both even and odd n. We then analytically continue the resulting Rényi negativities and find the values of three special limits that we defined in Section 2: the logarithmic negativity E, partially transposed entropy S T 2 , and refined Rényi-2 negativity S T 2 (2) .
We collect these results for all four phases in Table 1.
From these results, we find a phase diagram for negativity, which we show in Figure 4. Unlike the von Neumann entropy which only has a single phase transition at k ∼ e S 0 , we see that there are two distinct types of phase transitions for negativity, one from the disconnected phase to the pairwise phase and one from the pairwise phase to the cyclic phase. The transition from the pairwise phase to the anti-cyclic phase is similar to the pairwise-to-cyclic transition under the exchange k 1 ↔ k 2 .

Resolvent equation for partial transpose
Having analyzed the negativity measures deep within each of the four phases, we now turn our attention to the behavior of negativities near the phase transitions. Generally speaking, more geometries than the four types studied in the previous section could dominate the Rényi negativities near a phase transition, and we need to sum over them. We would then need to analytically continue the resulting sum to find special limits such as the logarithmic negativity. This is technically difficult to do directly.
Instead, we study the resolvent for the partially transposed density matrix ρ T 2 R , which we refer to as the negativity resolvent or simply the resolvent. From this resolvent, we then extract the negativity spectrum, i.e., the eigenvalue distribution of ρ T 2 R . This allows us to calculate the Rényi negativities and their special limits.
In this section, we derive a self-consistent equation for calculating the resolvent. As we will show, the sum over dominant diagrams in a "planar" regime reduces to a Schwinger-Dyson equation, which can be resummed. This allows us to write down a closed-form equation for the resolvent.
The negativity resolvent is defined in terms of the partially transposed density matrix as 4 (4.1) From the resolvent, the eigenvalue spectrum for ρ T 2 R , which we will denote by D(λ), can be obtained by taking the discontinuity across the real axis as follows From this, we can compute the Rényi negativities via from which all other negativity measures we consider can be derived. It will be useful to consider (4.1) in the following matrix form: where we denote the R 1 subsystem by upper i-type indices and R 2 by lower j-type indices. In the last line, we have expanded the expression in a formal power series in 1/λ. Each term in the series is given by the n-replicated density matrix ρ T 2 R n , which defines a boundary condition with n asymptotic boundaries with the brane indices contracted: where the blue (upper) dashed lines denote i-type index lines and red (lower) dotted lines denote j-type index lines. Each pair of blue/red index lines gives a factor of 1/λ and each asymptotic boundary gives a factor of 1/(kZ 1 ) coming from the normalization of the density matrix. Note that for up to two boundaries, the boundary conditions are the same after taking a final trace, with or without partial transpose. In general, the gravitational path integral can be performed as a sum over bulk geometries satisfying the boundary conditions. In the JT model we introduced in Section 3, higher genus corrections are highly suppressed by factors of e −S 0 1, so we only need to consider disjoint unions of disk geometries connecting any number of asymptotic boundaries. The path integral for the disk can be performed exactly including quantum corrections and is given by (3.11). For n asymptotic boundaries, these disjoint unions of disk geometries are in one-to-one correspondence with elements of the permutation group S n . As we show in Lemma 1 and Corollary 2 of Appendix A, the only geometries that can possibly dominate in the limit e S 0 1 are the planar and anti-planar diagrams, which correspond to certain subsets of permutations in S n . The term "anti-planar" will be explained in detail in Section 4.2 and in Appendix A. In fact, large regions of the phase diagram are dominated by either the planar or anti-planar geometries, which we will call the planar and anti-planar regimes. As we will now show, in each of these two regimes the resulting sum over diagrams can be resummed via a Schwinger-Dyson equation.

Planar regime
We now derive a resolvent equation in the parameter regime k 2 k 1 e S 0 . For reasons that will become clear shortly, we call it the planar regime.
Our strategy is to keep only the subset of geometries which have a possibility of dominating. As we outlined in Section 3.2, the phases in this regime away from phase transitions are defined by the disconnected, cyclic, and pairwise geometries. Closer to phase transitions, we expect more generic geometries which "interpolate" between these three types of geometries to have a chance of dominating. Indeed, as we show in Lemma 1 of Appendix A, the geometries which can dominate are precisely the planar diagrams, which correspond to the non-crossing partitions, i.e., the permutation group elements g ∈ S n lying on a geodesic between 1 and X.
We can write the sum over planar geometries diagrammatically as This sum can be recast as a Schwinger-Dyson equation. Diagrammatically, we have or, as an equation, where we have defined the partial trace of the resolvent matrix over the R 1 subsystem: and repeated indices are summed over 5 . Note that for n = 1 the product of resolvents is simply R i 1 i 2 j 1 j 2 and for n = 2 it isR k 1 k 1 R i 1 i 2 j 1 j 2 . As is evident from the diagrammatics, the i-type indices denoting R 1 form simple self-contractions on all but the last resolvent, while the j-type indices denoting R 2 form a complicated set of contractions. Fortunately, this equation can be solved iteratively: starting with the leading solution R i 1 i 2 As we will see, this allows us to rewrite the complicated product of resolvents as the following simple expressioñ Let us explain this in more detail. In general, the behavior for n odd and even are different so we will need to treat these cases separately. To illustrate the simplification for the odd case, we first consider the contribution at n = 3, which is the first non-trivial diagram under the partial transpose. We can writẽ where no summation on j 1 is implied. Now, recalling the full trace R = jR jj , we findR jj = R/k 2 and we can therefore writẽ We can understand the even case by looking at the first non-trivial contribution at n = 4. By a similar analysis as above, we havẽ Note the additional factor of k 2 compared to n = 3 due to the closed index loop formed from the first and third resolvent factors. More generally, one can show that the even case always has a single index loop and the odd case has no index loops, leading to (4.8). Using this, we can rewrite the Schwinger-Dyson equation (4.7) as Taking the full trace, we find The gravitational partition function of the n-boundary totally connected geometry is given by (3.11), which we repeat here: Since this depends on n only through y(s) n , the sum over n becomes a geometric series and (4.12) can be resummed into where w(s) ≡ y(s)/Z 1 . As a consistency check, when k 2 = 1 this reduces to the resolvent equation for the original (non-partially-transposed) density matrix derived in [1].

Anti-planar regime
We now consider a different parameter regime k 1 k 2 e S 0 . For reasons that will become clear shortly, we call it the anti-planar regime. This anti-planar regime has a large overlap with the planar regime; the overlap region is e −S 0 k 1 /k 2 e S 0 . Together the two regimes cover the entire parameter space.
Once again, we will keep only the subset of geometries which have a possibility of dominating. As we outlined in Section 3.2, the phases in this regime away from phase transitions are defined by the disconnected, anti-cyclic, and pairwise geometries. Closer to phase transitions, we expect geometries which interpolate between these geometries to have a chance of dominating. Indeed, as we show in Corollary 2 of Appendix A, the geometries which can dominate are precisely the set of anti-planar diagrams, which are in one-to-one correspondence with permutations lying on a geodesic between 1 and X −1 . An example of an anti-planar geometry is the anti-cyclic geometry; two other (perhaps less obvious) examples are the disconnected and pairwise geometries. Geometrically, anti-planar diagrams are precisely those that become planar diagrams after reversing the orientation of the asymptotic boundaries: The sum over anti-planar diagrams can similarly be recast as a Schwinger-Dyson equation. Diagrammatically, we have or, as an equation, where we have defined the partial trace over the R 2 subsystem: Note that the n = 1 and n = 2 terms are the same as in the planar case. Compared to the planar case, the j-type indices denoting R 2 now form simple self-contractions on all but the last resolvent, while the i-type indices denoting R 1 form the complicated set of contractions. In other words, the i-type indices now play the role of j-type indices in the planar regime, and vice-versa. As before, we can use the fact that R i 1 i 2 j 1 j 2 ∝ δ i 1 i 2 δ j 1 j 2 to all orders to rewrite the complicated product of resolvents as (4.16) Using this, we can rewrite the Schwinger-Dyson equation (4.15) as Taking the full trace, we find Finally, using (3.11), we resum (4.17) into (4.18) As expected, this is simply the resolvent equation in the planar regime (4.14) with the exchange k 1 ↔ k 2 .

Negativity spectrum near phase transitions
Having derived the resolvent equation for the partial transpose, we now solve it to find the negativity spectrum near phase transitions. For each negativity measure, we will be interested in a neighborhood near one of the phase transitions as we tune the relative sizes of our parameters, and we will analyze the corrections to the negativity measures listed in Table 1. Generically these corrections take the form of fluctuations about a fixed saddle point in the gravitational path integral. Near a phase transition, however, multiple saddles are competing for dominance, so enhanced corrections, i.e. corrections larger than those for any individual saddle, provide additional information about the entanglement structure of the system. We fix e S = 10, k 2 = 2, and move in a horizontal line in the phase diagram by tuning k. The disconnected-to-pairwise transition (left) occurs at k = e S = 10, though we only plot to k = 5 for visual clarity, as the qualitative behavior is the same. The pairwise-to-cyclic transition (right) occurs at k = k 2 2 e S = 40.

Microcanonical ensemble
Before studying the negativity spectrum in more detail, let us take a brief detour and consider the situation where the black hole is in a microcanonical ensemble in the JT model. In this case, we restrict to some small energy window [s, s + ∆s]. We write where ρ(s) ≡ e S 0 ρ(s) is the density of states. In this case, it can be shown that the sum over geometries for the Rényi negativities (3.12) is given by We recognize this as the n th Rényi negativity of a Wishart matrix with e S degrees of freedom. This implies that ρ T 2 R in a microcanonical ensemble can be thought of as the partial transpose of a random matrix drawn from the Wishart distribution. A similar expression for the moments of the partial transpose of a random mixed state was derived in [29]. The similarity between the microcanonical JT model and a random mixed state was also noted in [36].
In the planar regime k 1 k 2 e −S 0 , the resolvent equation (4.14) becomes This matches the resolvent equation derived in [29] for a random mixed state under appropriate rescaling of variables. 6 As was shown there, a closed-form solution to this cubic equation for R can be found, and leads to concrete results for the negativity spectrum and various negativity measures. The resolvent for the anti-planar regime k 2 k 1 e −S 0 can be obtained by k 1 ↔ k 2 . We plot the eigenvalue density in the microcanonical ensemble for various parameter values in Figure 5. The spectrum is approximately a Wigner semicircle distribution in the disconnected phase, continues to be connected in the pairwise phase, develops singularities at the pairwise-to-cyclic transition, and has two branches in the cyclic phase, where it is well approximated by the difference of two disjoint Marchenko-Pastur distributions.

Canonical ensemble: disconnected-pairwise transition
The transitions that involve the cyclic and anti-cyclic phases are complicated, as they involve a sum over diagrams with pieces connecting more than two asymptotic boundaries. Here, we will focus on the transition between the totally disconnected phase and the pairwise connected phase. The disconnected phase involves single-boundary diagrams, while the pairwise phase involves pairwise connected wormholes (plus a single disconnected piece for odd n).
The disconnected-pairwise transition happens within the large overlap e −S 0 k 1 /k 2 e S 0 between the planar regime and the anti-planar regime 7 . Therefore, the dominant geometries are those that are simultaneously planar and anti-planar. As we show in Appendix A, these geometries are disjoint, non-crossing unions of single-boundary disks and pairwise connected wormholes. This result is the content of Lemma 3 in Appendix A.
Intuitively, these dominant geometries interpolate between the disconnected and pairwise geometries. Geometries with pieces connecting more than two asymptotic boundaries are parametrically suppressed. As such, the resolvent equation (4.12) trun- to m = 1/2 [30].
cates at the quadratic order: This quadratic equation can be solved analytically giving the resolvent and eigenvalue density as where A 2 ≡ 4Z 2 /(kZ 2 1 ). Thus the eigenvalue density is a Wigner semicircle distribution supported on λ ∈ −A + 1 k , A + 1 k . For k ≤ 1/A, D(λ) only has support on λ ≥ 0 and the negativity vanishes; for k > 1/A, D(λ) has support on λ < 0 and we find the negativity is non-vanishing. The phase transition thus occurs at k = 1/A ∼ e S 0 , as expected from the schematic analysis in Section 3.2.
We can write an explicit expression for the logarithmic negativity using (2.6) and (5.6): We plot this result in Figure 6, along with the naive answers for logarithmic negativity in the dominant phases. How large is the correction at the transition? It is easy to verify that it is O(1). There are no enhanced corrections here, as we are working in a regime where higher order terms in the Schwinger-Dyson equation are parametrically suppressed. The logarithmic negativity, along with all other negativity measures, never receives contributions from geometries containing pieces with Z n>2 , so corrections are O(1).

Canonical ensemble: cyclic-pairwise transition
In this subsection we will study the richer phase transition between the pairwise phase and the cyclic phase. 8 Our computation is inspired by that of Appendix F of [1].

(5.9)
This means that the integral that defines Z n in (3.11) can be well approximated by the saddle point located at Throughout, we take our parameters S 0 , k 1 , and k 2 to be large before taking the semiclassical limit, such that e.g. log Z n ≈ S 0 as S 0 1/β. We will also need to define s k , the value of s for which We can approximate s k by Note that here we are considering the values of k and k 2 at transition, so the particular values of s k we are interested in will depend on the details of the negativity measure we are computing. In our schematic analysis where Z n ∼ e S 0 , we derived the location of the phase transition between the cyclic and pairwise phase and found that it was independent of n. However, taking into account dependence on β, the Z n 's are distinct for different n, which leads to n-dependent transition points. The Rényi negativities in the cyclic and pairwise phases are given by the contributions of the dominant geometries in each phase (see Table 1), and coincide at transition. Equating their contributions at the transition, we find, up to factors O(1) in β, In terms of the approximation (5.10), we can solve for log (k/k 2 2 ) at transition to obtain (5.14) From this, we can solve for s k at the transition using (5.12) to obtain s (n, even) k As expected, the transition point depends on n. In particular, it is O(1/β) at leading order and bounded below as a function of n.
For this phase transition, we will fix k and tune k 2 . In the phase diagram, this corresponds to moving along a line between the upper left corner and the lower right corner. We need to consider diagrams with pieces made of an arbitrary number of boundaries, and we can restrict ourselves to planar diagrams, as anti-planar diagrams are suppressed by factors of k 2 /k 1 relative to their planar counterparts. The resolvent equation is again (4.18): We are going to split this integral at some transition s t such that We rewrite this simple step to emphasize that no approximations have been used yet.
We are now going to use a set of three assumptions: These assumptions are justified in detail in Appendix B, where we show that the resulting simplifications to the resolvent equation give a self-consistent treatment of the problem. For now we will take these as facts and proceed. The first approximation allows us to simplify the final term in (5.17) such that where we have dropped an R 2 term from the last integral because it can be shown to be much smaller than the leading term k using Assumptions 1 and 3. We define the coefficient of R in the last term to be the constant λ 0 , given by We can now write the resolvent equation as (5.20) Now we turn to Assumption 2, which allows us to treat the second term above as a perturbation to the zeroth order solution Plugging this solution back into (5.20), we obtain the first order iterated solution Now we can find the discontinuity in this expression and extract D(λ). There are three contributions to the spectrum: a simple pole at λ = λ 0 , and a pair of branch cuts given by the poles at λ = λ 0 ± w(s)/k 2 in the integrand. We obtain Let us pause for a second to unpack this equation. The spectrum consists of a delta function located at λ 0 from the simple pole and two regions of nonzero eigenvalue density from the integrated delta functions. We plot a sketch of this eigenvalue density in Figure 7. There are two distinct regions with nonzero eigenvalue density, similar to the spectrum in the microcanonical ensemble. One point we emphasize in Appendix B is the presence of a "controlled" region 0 < s < s c in which our assumptions hold and an "uncontrolled" region s > s c where we claim ignorance about the spectrum. In terms of the eigenvalues, this corresponds to an ignorance in the spectrum for a region We show that our ignorance about the uncontrolled region leads to at most O(1) multiplicative corrections to the Rényi negativities, or O(1) additive corrections to E, S T 2 , and S T 2 (2) , due to the constraint that the total number of eigenvalues must be k. Clearly, λ 0 lies within this uncontrolled region, so we should not take seriously the presence of the delta function at λ 0 . In fact, we will now show that this delta function vanishes if we extend the upper limit of the integral from s t to s k in (5.23) (which only affects the uncontrolled region and therefore causes a small error). Our density matrix has a total of k eigenvalues and is unit normalized, which translates into conditions on the zeroth and first moments of D(λ), namely We see that these conditions are satisfied by (5.23) if we replace s t by s k , compensating for the fact that s t = s k − κ by sending the coefficient of the δ(λ − λ 0 ) piece to zero. We conclude that a good approximation for the spectrum of the partially transposed density matrix at transition is given by 26) In order to simplify our calculations for negativity, we would like that the delta function branch cuts were purely positive or negative, which is equivalent to the condition λ 0 < w(s k )/k 2 . By definition, λ 0 is bounded above by 1/k, and at transition we have k/k 2 2 = e S 0 +2πs k . In the semiclassical limit, we can use the approximation Our condition on the branch cuts becomes kw(s k ) k 2 > 1 ⇒ k 2 e S 0 +2πs k w(s k ) ≈ k 2 e 2πs k −βs 2 k /2−2π 2 /β ∼ k 2 e C/β 1 (5.28) as s k ∼ 1/β for a generic n. This is satisfied even if the unknown order one constant is negative under our previous assumption that we take our counting parameters to be large before taking small β, such that log k 2 1 β . Now that we have the spectrum at transition, we are ready to calculate the corrections to any negativity measure we want! Let us start with the logarithmic negativity, which has a transition located at We find As s k < s (1) for logarithmic negativity, we have approximated λ 0 ≈ 1/k, and the final integral is well approximated by its maximum value ρ(s k )w(s k ). The logarithmic negativity is then In the second line we used our previous approximation (5.12) for s k , and in the third line we used (5.28). As we see, the logarithmic negativity experiences an O(1/β) correction to the naive answer E = log k 2 .
Where do we expect O(1/ √ β) corrections? In the case of even Rényi negativities, this happens at s We can check this explicitly for a negativity measure descending from the even analytic continuation. The simplest such measure, the Rényi-2 negativity N (even) 2 , is related to the second Rényi entropy S 2 by N (even) 2 = e −S 2 (5.33) and therefore comes with O(1) corrections. We instead turn to the refined Rényi-2 negativity S T 2 (2) , defined in (2.12). We can read off the naive answer for S T 2 (2) from Table 1, again using the approximation (5.10). We find To compute S T 2 (2) , we first need to compute i λ 2 i = N (even) 2 . At transition, the naive answer is given by Table 1, where N (even) 2 = Z 2 /Z 2 1 . However, using (5.26), we find Again, s k < s (1) , so λ 0 ≈ 1/k. This integral gives N (even) 2 In the limit k e S 0 , we can safely ignore the first two terms, and the Rényi-2 negativity becomes The factor of 1/2 out front may seem like a problem, as we do not reproduce the naive answer for N (even) 2 . However, as the refined Rényi negativities are functions of log N n , this factor will only contribute an O(1) difference from the true answer for S T 2 (2) , and we are safe in using this approximation. The refined Rényi-2 negativity is therefore given by Again, as λ 0 ≈ 1/k, the dominant contribution to this integral will come from the w(s) 2 term, which is where the enhanced transition should come in. Previously, using (5.9), we showed that an integral of the form ρ(s)w(s) n can be approximated by a sharply peaked Gaussian with mean s (n) and standard deviation 1/ √ nβ, up to normalization.
We use these simplifications to obtain − e S 0 N (even) 2 Our final expression for S T 2 (2) is therefore confirming that there is an O(1/ √ β) correction at transition. The fact that the refined Rényi-2 negativity experiences this particular correction is not surprising due to its close connection to von Neumann entropies. It is known that von Neumann entropies receive O(1/ √ β) or O(1/ √ G N ) corrections at the Page transition [1,33,34], which can be explained using a diagonal approximation with respect to a basis of fixed-area states [33,34,44]. It was shown that the refined Rényi-2 negativity can be written in holography as the sum of the von Neumann entropies of R 1 and R 2 in the state ρ 2 R 1 R 2 (once properly normalized) [30]. 9 Therefore, the O(1/ √ β) correction that we find in the refined Rényi-2 negativity can similarly be explained using the diagonal approximation.
As we show in Appendix C, the Rényi entropy S n with n < 1 experiences O(1/β) corrections in the model of [1], as there too we are computing an entanglement measure with s (n) k < s (n) . In other words, the Rényi index of both the logarithmic negativity and the Rényi entropy with n < 1 is below some "critical" Rényi index at which there exist O(1/ √ G) corrections. For measures descending from odd Rényi negativity, we might not expect O(1/ √ β) corrections, as we never have s (n,odd) k = s (n) . However, we may still expect some enhanced corrections for some negativity measures in this case. The partially transposed entropy S T 2 is one such measure. As s (1,odd) k = 5π/2β, the naive answer for S T 2 is given by Our approximation gives There is however a subtlety here. The dominant contribution to this integral no longer comes solely from the w(s) term. This can be seen by expanding (5.43) using From this we obtain Treating the log w(s) k 2 as negligible compared to the exponential ρ(s), the two terms in the integrand are of the same order, so we should keep them both. If we look at the naive transition point s k = 5π/2β > s (1) , we find and we would conclude that the correction is O(1). However, if we were to find the largest correction to this quantity, we would look not at the naive transition, but at the point where we might find O(1/ √ β) corrections, at s k = s (1) . At this point λ 0 = 1/2k and we capture half of the Gaussian ρ(s)w(s), so we have This looks the same as the naive answer with a O(1/ √ β) correction. However, as we are working at fixed k, we should really be writing everything in terms of k and e S 0 using (5.12), in which case our naive and corrected S T 2 's are Naive: Corrected: and we find an O(1/β) correction.

Topological model with EOW branes
Having studied a toy model of an evaporating black hole in JT gravity, we will now consider entanglement negativity in the context of the topological model of Marolf-Maxfield [35], including dynamical end-of-the-world (EOW) branes. This is a theory of topological two-dimensional gravity in which spacetimes are two-dimensional manifolds endowed only with orientation. In contrast to the JT model of Section 3, there are no metric or dilaton degrees of freedom and EOW brane boundaries can be generated dynamically.
The action for the topological model is given by where S 0 is some arbitrary parameter, χ(M ) is the Euler characteristic of the (possibly disconnected) manifold M , and |∂M | counts the number of boundaries. S ∂ |∂M | is a non-local term that we put in by hand to ensure reflection positivity. As shown in [35], the simplest choice which results in reflection positivity is S ∂ = S 0 . 10 The action then becomes S top = −S 0χ (6.2) whereχ = 2 − 2g for any manifold, with or without boundary. To make contact with black hole evaporation, we can extend the model to include EOW branes, which can take one of k "flavors". Since we are interested in studying negativity, we will allow the branes to be labeled by a set of two of flavor indices {i, j}, where i ∈ {1, . . . , k 1 }, j ∈ {1, . . . , k 2 }, and such that k = k 1 k 2 . This is exactly analogous to our construction in Section 3 and is a slight generalization of the model in [35].
There are three distinct types of boundaries allowed by the theory. The first are circular asymptotically AdS boundaries denoted by Z. These boundaries are associated with an operator Z which acts on the baby universe Hilbert space H BU and creates a Z boundary. Second, there are boundary conditions which we denote by (ψ i 2 j 2 , ψ i 1 j 1 ) composed of an oriented interval of asymptotically AdS boundaries with endpoints labeled by flavor indices {i 1 , j 1 } and {i 2 , j 2 }. The diagram that describes this is the same as in (3.6). These boundaries are associated with an operator (ψ i 2 j 2 , ψ i 1 j 1 ) on H BU . Finally, there are circular EOW brane boundaries labeled by an arbitrary flavor index {i, j}, independent of all boundary conditions. These brane boundaries can be dynamically generated as additional boundaries when performing the gravitational path integral.
Let us consider the simplest quantity one can compute with the gravitational path integral of this theory, namely the partition function associated to a single connected component of spacetime with some number of asymptotic boundaries. The gravitational path integral demands that we sum over all such manifolds with arbitrary genus, weighted by e −S 0χ . Additionally, since the EOW branes are dynamical, there is the possibility of the gravitational path integral generating an arbitrary number of closed brane boundaries, each of which contribute a factor of k and are mutually indistinguishable. The partition function for a single connected component with some fixed number of asymptotic boundaries is therefore More generally, one can consider amplitudes Z m (ψ i 1 j 1 , ψ i 1 j 1 ) · · · (ψ i n j n , ψ injn ) ≡ NB Z m (ψ i 1 j 1 , ψ i 1 j 1 ) · · · (ψ i n j n , ψ injn ) NB (6.4) which are computed using the gravitational path integral by summing over all (possibly disconnected) manifolds with boundary conditions specified by m circular boundaries Z and n oriented intervals (ψ i j , ψ ij ) with endpoints labeled by the corresponding flavor indices and connected to oriented brane boundaries labeled with matching flavors. The brackets in (6.4) can be interpreted the expectation value of the corresponding operators in the no-boundary state NB ∈ H BU . In what follows, we will assume that the no-boundary state is unit normalized, NB NB = 1. 11 Let us now proceed to the calculation of negativity in this model. In analogy to Section 3, we can define the (unnormalized) density matrix which plays the role of the state of the Hawking radiation. We are interested in studying the Rényi negativities N n = Tr ρ T 2 n of this density matrix. The most straightforward method is to use the moment generating function of ρ i 1 i 2 j 1 j 2 ≡ (ψ i 2 j 2 , ψ i 1 j 1 ), which one can show is given by where t can be thought of as the k × k matrix with entries t i 1 i 2 j 1 j 2 by treating {i, j} as a single index of size k and I is the k ×k identity matrix . This is a slight generalization of the result derived in [35]. In principle, one can compute all moments of ρ i 1 i 2 j 1 j 2 , and hence all Rényi negativities, by taking appropriate partial derivatives of (6.6) with respect to with the distribution of a random mixed state and the microcanonical JT model in Section 5.1. Thus, we can immediately write down the Rényi negativities in a fixed Z = d sector which matches the answer for the microcanonical ensemble (5.2) with d playing the role of e S . The results for the negativity spectrum obtained in Sections 3.2 and 5.1 therefore apply with this replacement. However, since d does not correspond to the partition function on some manifold, it is difficult to interpret the result in (6.8) geometrically.
To obtain the Rényi negativities in the full theory, we simply sum over d ∈ N with Poisson weight p d (λ): where B m (λ) = e −x ∞ k=0 λ k k m m! are the Bell polynomials, whose asymptotic behavior is B m (λ) ∼ λ m as λ → ∞. We therefore find This is once again equivalent to the microcanonical ensemble in (5.2), with λ now playing the role of e S , and therefore we can obtain concrete results for the negativity spectrum. Since λ is the gravitational partition function of a single connected component of spacetime, we can in fact find a geometric interpretation for the terms in (6.9).
To understand the geometric origins of the terms in (6.9), let us first look at the case n = 2, which gives the purity Tr ρ 2 = λ 2 k + λk 2 + λk. (6.11) The terms in (6.11) correspond to the following geometries: The first two diagrams are familiar: they are the disk and wormhole geometries, summed over genus and closed brane boundaries. The last diagram represents two disk geometries joined by an arbitrary number of wormholes; we thus call it a joining wormhole.
More generally, the geometries which contribute at leading order in λ in the Rényi negativities (6.9) are in one-to-one correspondence with elements of the permutation group. To be precise, each of these geometries is actually a disjoint union of disks, summed over genus and closed brane boundaries. 12 The subleading contributions in (6.9) can be identified with the same geometries but with arbitrary numbers of joining wormholes between connected components, and thus can not be mapped to elements of the permutation group.
It is clear that log λ plays the same role as S in the microcanonical JT model, namely it is the Bekenstein-Hawking entropy. In analogy to black hole evaporation, we should assume λ 1. The joining wormholes are therefore parametrically suppressed, but disks with handles are not since they instead come with factors of e −2gS 0 and S 0 is not a priori a large parameter (in fact, it may have a small or even negative real part). 13 This is the analogue of the "planar" limit in the topological model. There are thus two distinct classes of higher genus geometries: disks with handles and joining wormholes. The higher genus disk geometries can be systematically included in a Schwinger-Dyson equation as in Section 4, while the joining wormholes can not.
To study the Page curve, we would like to fix the value of λ and tune k. Since λ ∼ e k , this involves scaling the prefactor e 2S 0 1−e −2S 0 with e −k . However, this function has a minimum value for real S 0 , which means the black hole can not evaporate completely. To decrease λ beyond the minimum value, we need to go to complex values of S 0 , namely e 2S 0 ∈ 1 2 + iR. This is a bit strange because it implies a complex action, but is presumably fine because S 0 is not a physical parameter (it is not the Bekenstein-Hawking entropy here). This is simply a quirk of the model, and can be attributed as a consequence of having a non-vanishing S ∂ .

Discussion
In this paper, we analyzed the behavior of negativity measures in toy models of evaporating black holes in both JT gravity and a topological theory of gravity, with EOW branes. We found four distinct phases dominated by different saddle-point geometries: the disconnected, cyclically connected, anti-cyclically connected, and pairwise connected. The last of these geometries are new replica wormholes that break the replica symmetry spontaneously.
We also studied the negativity resolvent using a Schwinger-Dyson equation that resums the contributions of different geometries, and used it to extract the negativity spectrum and negativity measures. This analysis is valid not only within each of the four phases, but also near phase transitions. For the topological model or a microcanonical ensemble in JT gravity, we found a cubic equation for the resolvent which can be solved exactly. For a canonical ensemble in JT gravity, we found a quadratic resolvent equation near the disconnected-pairwise transition, and we solved a more complicated resolvent equation approximately near the cyclic-pairwise transition. Near this last transition, we found enhanced corrections to various negativity measures: the refined Rényi-2 negativity receives an O(1/ √ β) correction, whereas the logarithmic negativity and the partially transposed entropy receive O(1/β) corrections.
These enhanced corrections to negativities are similar to previously found corrections to the von Neumann entropy at the Page transition [1,33,34]. For the von Neumann entropy, the enhanced corrections can be explained using a diagonal approximation with respect to a basis of fixed-area states [33,34,44]. We argued that the O(1/ √ β) correction to the refined Rényi-2 negativity can be explained in the same way by noting its close connection to von Neumann entropies. It would be interesting to understand further the O(1/β) corrections to the logarithmic negativity and the partially transposed entropy in a similar way. Moreover, it would be useful to study the implications of these O(1/β) corrections for the partially transposed entropy more generally: it was conjectured in [30] that S T 2 (ρ R 1 R 2 ) is given as a sum of von Neumann entropies (S R 1 + S R 2 + S R 1 R 2 )/2 in general non-fixed-area states by assuming a diagonal approximation, but this would imply an O(1/ √ β) correction and seems to be in tension with the O(1/β) correction that we find here.
We focused our study on two specific toy models of evaporating black holes, but it would be interesting to generalize our analysis to other models, including the examples studied in [2].
Finally, it would be very interesting to use these results on negativity to diagnose the structure of multipartite entanglement in a realistic evaporating black hole and learn more about its quantum state. We hope that this will lead to new insights on understanding the interior of black holes and the dynamics of their evaporation.

A Derivation of dominant saddles for negativity
In this appendix, we derive the set of saddle-point geometries that give dominant contributions to the Rényi negativity in various regimes of the parameter space. This includes each of the distinct phases and near phase transitions.
Our derivation uses facts about geodesics on the permutation group, which we review first. Let S n be the symmetric group of order n, which is the set of permutations on n elements. For any permutation g ∈ S n , we define (g) as the minimum number of swaps from the identity 1 = (1)(2) · · · (n) to g and χ(g) as the number of disjoint cycles in g, including 1-cycles. These quantities satisfy the relations (g) + χ(g) = n, (A.1) As an example, the permutation 14 g = (12)(345) ∈ S 5 has (g) = 3 and χ(g) = 2. We can define the distance between two permutations g and h by which satisfies the usual properties of a distance measure. In particular, given any sequence of permutations (g 1 , · · · , g m ), the distance satisfies the triangle inequality A sequence of permutations that saturates (A.4) is said to be on a geodesic. We denote a geodesic between two permutations g and h by G(g, h). We say that a permutation g is on G(g, h), or equivalently g ∈ G(g, h), if the sequence (g, g , h) saturates the triangle inequality (A.4).
Our goal is to identify the permutations g that dominate the sum in the Rényi negativity (3.13), in different regimes of the parameter space labeled by e S 0 , k 1 , and k 2 . We repeat the sum in (3.13) here: (A.5) For reasons that will become clear shortly, it is useful to first identify the permutations on one or more of the following geodesics: G(1, X), G(1, X −1 ), and G(X, X −1 ). Here X = (12 · · · n) is the cyclic permutation of n elements, and X −1 = (n · · · 21) is the anti-cyclic permutation.
For a given permutation g, let us use m, p, q to denote the three exponents in the sum (A.5): They satisfy three triangle inequalities, which can be obtained from (A.1), (A.3), and (A.4): where f (n) is a useful function defined as f (n) ≡ 1, n odd, 2, n even, (A. 10) and we have used χ(X) = χ(X −1 ) = 1, χ(X 2 ) = f (n). We now identify the permutations g on one or more of the three geodesics.
Permutations on G(1, X): These are known to be in one-to-one correspondence with non-crossing partitions, so we say that the corresponding geometries are planar.
We can write such a element as a product of m non-crossing cycles (including 1-cycles): It is clear that such an element exists for every m ∈ [1, n]. Since it saturates (A.7), we immediately find p = n − m + 1.
Moreover, it is straightforward to derive where |c i | is the length of the i-th cycle c i .
Permutations on G(1, X −1 ): These can be obtained by simply taking the inverse of the permutations on G(1, X), sending c i in (A.11) to c −1 i . We say that these correspond to "anti-planar" geometries.
Permutations on G(1, X) and G(1, X −1 ): Their cycles c i must be their own inverses, so the length of each cycle is at most 2. Therefore, these permutations are precisely those non-crossing partitions that consist of only 1-cycles and 2-cycles. Such permutations exist for every m ≥ n 2 , with the lower bound saturated by non-crossing pairings consisting of n 2 pairs and at most one 1-cycle. Permutations on G(1, X) and G(X, X −1 ): As they saturate (A.9), it is straightforward to use (A.12) and (A.13) to show that these permutations are precisely those non-crossing partitions with at most one odd cycle. Here we define an odd cycle as one of odd length and an even cycle as one of even length. For even n, these permutations consist of only even cycles, whereas for odd n, they have exactly one odd cycle. Such permutations exist for every m ≤ n 2 , with the upper bound saturated by non-crossing pairings.
Permutations on G(1, X −1 ) and G(X, X −1 ): These are obtained by taking the inverse of the permutations on G(1, X) and G(X, X −1 ).
Permutations on G(1, X), G(1, X −1 ), and G(X, X −1 ): It is clear by combining the previous cases that these permutations are those non-crossing partitions that consist of only 2-cycles and at most one 1-cycle. Therefore, they are in one-to-one correspondence with non-crossing pairings [30]. These all have m = n 2 and p = q = n 2 + 1. We denote these non-crossing pairings by τ , and say that they correspond to pairwise geometries. A simple example is τ = (12)(34) · · · (n − 1, n) for even n and τ = (12)(34) · · · (n − 2, n − 1)(n) for odd n.
These results are illustrated schematically in Figure 8. We now state and prove the main points of this appendix. Lemma 1. In the regime k 1 /k 2 e −S 0 , planar geometries dominate the sum (A.5). In other words, for any g / ∈ G(1, X), there exists g ∈ G(1, X) such that g dominates over g. Proof. As g / ∈ G(1, X), the difference is positive. However, this difference must be even, regardless of whether g is an even or odd permutation. Therefore, we denote this difference by 2r with a positive integer r, and use it to rewrite the triangle inequality (A.7) as Our goal is to choose a more dominant g . For g , we define m , p , q similarly as Let us discuss m+r ≥ n 2 and m+r < n 2 separately. For m+r ≥ n 2 , we choose g to be on G(1, X) and G(1, X −1 ), with m = m + r. As we discussed earlier, such a permutation exists -in particular, (A.15) guarantees m + r < n. As g saturates (A.7) and (A.8) (after primes are added), we find where the second equality for p comes from (A.15) and the inequality comes from (A.8) for g. We thus find that g gives a more dominant contribution to the sum (A.5) than g, as e m S 0 k p 1 k q 2 e mS 0 k p 1 k q 2 ≥ e S 0 k 1 k 2 r 1.
(A. 18) In the other case with m + r < n 2 , we choose g to be on G(1, X) and G(X, X −1 ), with m = m + r. Again, such a permutation exists. As g saturates (A.7) and (A.9) (after primes are added), we find where the inequality comes from (A.9) for g. We again find that (A.18) holds and therefore g dominates over g.
It is clear that Lemma 1 is tight in the sense that every planar geometry could give a dominant contribution to the sum (A.5) at some point in the regime k 1 /k 2 e −S 0 . In particular, at the point where k 1 = e S 0 and k 2 = 1, they give an equal contribution e (n+1)S 0 .
From Lemma 1, we immediately obtain the following corollary by taking the inverse of all permutations and switching k 1 ↔ k 2 .
Combining Lemma 1 and Corollary 2, and recalling that the permutations on G(1, X) and G(1, X −1 ) are precisely those that consist of only 1-cycles and 2-cycles, we immediately obtain the following corollary (which is useful for studying the disconnectedpairwise transition in Section 5.2). It is again clear that Corollary 3 is tight in the sense that every permutation on G(1, X) and G(1, X −1 ) could give a dominant contribution to the sum (A.5) at some point in the regime e −S 0 k 1 /k 2 e S 0 . In particular, at the point where k 1 = k 2 = e S 0 /2 , they all give an equal contribution e (n+1)S 0 . Lemma 4. In the regime k 1 k 2 e S 0 , the disconnected geometry dominates the sum (A.5).
Proof. The disconnected geometry is represented by the identity 1 and contributes e nS 0 k 1 k 2 . For any other permutation g, we have m ≤ n − 1. From (A.7) and (A.8), we obtain p, q ≤ n − m + 1. We thus find that 1 gives a more dominant contribution to the sum (A.5) than g, as Lemma 5. In the regime k 1 /k 2 e S 0 , the cyclic geometry dominates the sum (A.5).
Proof. The cyclic geometry is represented by X and contributes e S 0 k n 1 k . For any other permutation g, we have p ≤ n−1. From (A.7) and (A.9), we obtain m ≤ n−p+1 and q ≤ n + f (n) − p. We thus find that X gives a more dominant contribution to the sum (A.5) than g, as From Lemma 5, we immediately obtain the following corollary by taking the inverse of all permutations and switching k 1 ↔ k 2 . Corollary 6. In the regime k 2 /k 1 e S 0 , the anti-cyclic geometry dominates the sum (A.5).
We now show that pairwise geometries dominate a fourth phase. We first derive the following lemma as a useful intermediate step.
Lemma 7. In the regime k 1 k 2 e S 0 , the permutations on G(X, X −1 ) dominate the sum (A.5). In other words, for any g / ∈ G(X, X −1 ), there exists g ∈ G(X, X −1 ) such that g dominates over g.
Proof. As g / ∈ G(X, X −1 ), the triangle inequality (A.9) must fail to saturate by a positive but even integer, which is at least 2: p + q ≤ n + f (n) − 2 = 2 n 2 .
(A. 22) Therefore, one of p, q must be no greater than n 2 . Without loss of generality, we consider the case of p ≤ n 2 . We then choose g to be on G(1, X) and G(X, X −1 ), with p = p + 1. As g saturates (A.7) and (A.9) (after primes are added), we find where the two inequalities comes from (A.7) and (A.22), respectively. As we discussed earlier, such a permutation g exists, as m = n − p ≥ n 2 . From this, we find that g gives a more dominant contribution to the sum (A.5) than g, as Combining Lemmas 1, 7 and Corollary 2, and recalling that the permutations on all three geodesics G(1, X), G(1, X −1 ), and G(X, X −1 ) are precisely non-crossing pairings that lead to pairwise geometries, we immediately obtain the following corollary.
Corollary 8. In the regime satisfying both k 1 k 2 e S 0 and e −S 0 k 1 /k 2 e S 0 , the pairwise geometries dominate the sum (A.5).
can therefore rewrite the integral with an i prescription as e S 0 dsρ(s)f (s) = P V e S 0 dsρ(s)f (s) ±iπk 2 2 e S 0 dsρ(s) w(s)R(k + w(s)R) ∂ s (k 2 k 2 2 − w(s) 2 R 2 ) δ(s−s * ), (B.4) where PV denotes the Cauchy principal value. We choose the sign of i arbitrarily, as we are only looking to bound the absolute value of this integral.
Let us now look at the second term in (B.4). We have where we have used |w(s * )R| = kk 2 k to simplify the numerator. We can rewrite this in terms of s k using (5.11) and the asymptotic form of ρ(s): As stated previously, s k ∼ O(1/β), so for this to be much smaller than k we require s k − s * being at least O(1) but large, and by proxy κ ≡ s k − s t being at least O(1) but large, as stated in assumption (3). What is stopping κ from being much larger, say O(1/β)? Now we check the validity of assumption (1), that is |w(s t )R| kk 2 . Under our approximation (B.2), this assumption translates into the condition However, at the boundary of our spectrum located at s = s t , we should also have |λ − λ 0 | = w(s t )/k 2 . The way to reconcile these two assumptions is to state that our approximation only holds for some range of s between 0 and some control parameter Under our previous assumption that κ is O(1) in β such that s t ∼ 1/β, κ too must be O(1) but large. This answers our previous question about the size of κ. There is something of an inverse relationship between s t and κ : our goal is to design an approximation in which s c is as close as possible to s t , as well as one in which s t is as close as possible to s k . Under the assumption that s t is as close as possible to s k , that is κ is O(1) but large, we also have κ O(1) but large. If s t was far from s k so that, say, s t ∼ 1/ √ β, we would also require κ ∼ O(1)/ √ β with a large O(1) constant, thus missing a large part of the spectrum in our approximation.
Our conclusion is that there exists a region of size O(1) in s where assumption (1) does not hold. As the spectrum is over a region of size w(0)−w(s k ) k 2 , the size of a region O(1) in s is exponentially suppressed in 1/β, and we conclude that very few eigenvalues are in the uncontrolled region.

C Rényi entropies near the Page transition
We draw an analogy between the O(1/β) corrections to the logarithmic negativity and the partially transposed entropy and the O(1/β) corrections to the Rényi entropy S n with n < 1. Here we show this result explicitly in the model of [1]. We recall many of their results, which can equivalently be obtained from ours by sending k 1 to k and k 2 to 1. We consider the model of Section 3, but without partitioning the radiation system. The approximation for the density of states (5.26) is now D(λ) = e S 0 s k 0 dsρ(s)δ(λ − λ 0 − w(s)), (C.1) where w(s) and λ 0 are defined as in the main text. Here s k is defined as The transition at the Page time can be thought of as the transition from the fully disconnected phase to the cyclic phase along the x axis of our phase diagram ( Figure  4). Using our results from Table 1  As s (n) = 2π nβ , for n < 1 we have s (n) k < s (n) . Part of our derivation relied on s k scaling like 1/β, which remains true for n of O(1).
The Rényi entropy is given by As s k > s (1) , λ 0 will be exponentially suppressed in 1/β, so this integral is dominated by the w(s) term, and as s We conclude that there are enhanced corrections of the form O(1/β) to the Rényi entropy for n < 1.