Shape Dependence of Entanglement Entropy in Conformal Field Theories

We study universal features in the shape dependence of entanglement entropy in the vacuum state of a conformal field theory (CFT) on $\mathbb{R}^{1,d-1}$. We consider the entanglement entropy across a deformed planar or spherical entangling surface in terms of a perturbative expansion in the infinitesimal shape deformation. In particular, we focus on the second order term in this expansion, known as the entanglement density. This quantity is known to be non-positive by the strong-subadditivity property. We show from a purely field theory calculation that the non-local part of the entanglement density in any CFT is universal, and proportional to the coefficient $C_T$ appearing in the two-point function of stress tensors in that CFT. As applications of our result, we prove the conjectured universality of the corner term coefficient $\frac{\sigma}{C_T}=\frac{\pi^2}{24}$ in $d=3$ CFTs, and the holographic Mezei formula for entanglement entropy across deformed spheres.


Introduction
Entanglement entropy has a central position in the study of quantum field theories. It is a powerful tool to probe the structure of quantum states, primarily because: (i) it is sufficiently non-local to capture certain global properties, and (ii) it is geometric by definition and hence universal in its applicability. As a result, entanglement entropy has provided great insights in a wide class of systems such as relativistic field theories [1][2][3], conformal field theories (CFTs) [4][5][6], topologically ordered phases of matter [7][8][9][10], strongly-coupled theories with holographic duals [11][12][13], etc. It has also become clear that entanglement will play a crucial role in understanding the emergence of geometry in the AdS/CFT correspondence [14,15]. Despite this, computing entanglement entropy for arbitrary shaped regions in general dimension still remains a non-trivial task, especially outside the arena of quantum field theories with classical gravitational duals. While much progress can be made in special symmetric cases such as the entanglement entropy across planar surfaces in relativistic quantum field theories, spherical surfaces in CFTs, etc., it is desirable to develop a larger theoretical toolkit.
x i Figure 1: The set-up for the planar case: the original subregion A is the half-space x 1 > 0, with the entangling surface x 1 = 0 (dashed line). We deform the entangling surface to x 1 = − χ(x i ) (bold line) by glueing on the area elements δA a,b at points x a,b along the entangling surface.
In this paper, we study entanglement entropy for deformed half-spaces and ball-shaped regions in the vacuum state of a conformal field theory on R 1,d−1 . To be concrete, we first explain our construction for deformed half-spaces (see Fig. 1). Let us pick global coordinates x µ = (x 0 , x 1 , · · · , x d−1 ) on R 1,d−1 , where x 0 is the time coordinate. Pick the Cauchy surface x 0 = 0, and consider the reduced density matrix ρ 0 on the half space A given by The entanglement entropy between A and its complementĀ is defined as the Von Neumann entropy of ρ 0 . Next, deform the region A slightly to where χ(x 2 , · · · , x d−1 ) is a smooth function of the (d − 2) transverse spatial coordinates (parametrizing the entangling surface) which we denote collectively as x i = (x 2 , · · · , x d−1 ), and is a positive infinitesimal parameter. This corresponds to deforming the entangling surface to x 1 = − χ(x i ) within the original Cauchy surface. We can also generalize this and deform the entangling surface by the infinitesimal vector field ζ i (x i ), which we take to lie in the plane perpendicular to the original surface x 1 = 0 (i.e., the overlined indices run over i = 0, 1), and which also includes, for instance, time-like deformations. The entanglement entropy across the deformed surface can be written as a perturbative expansion in ζ i ij (x a , x b )+· · · (3) The quantity S (2) ij (x a , x b ) is known as the entanglement density [16][17][18] 1 , and will be the primary focus of the present paper. In the case where ζ i is spacelike and ζ 1 > 0, there is another nice way of thinking about the entanglement density: start with the half space A and glue on to it small area elements δ a A i and δ b A i at the points x a and x b on the entangling surface such that δ a A and δ b A are non overlapping. Then to lowest order in δ a,b A, the entanglement density is proportional to the conditional mutual information between δA a and δA b given the state on A δ a A i δ b A j S where I(X; Y ) is the mutual information between the regions X and Y . The strong subadditivity property then implies ij ≤ 0.
Consequently, the entanglement density provides a natural notion of a metric on the space of geometries of the entangling surface. In theories with holographic duals, the Ryu-Takayanagi proposal maps this space into the space of mimimal-area surfaces in the bulk, and so the 1 Sometimes entanglement density is defined with an extra minus sign to make it a naturally positive quantity, see (5). It is also not clear why one should think of it as a densityentanglement susceptibility would probably be a more appropriate name; however we will follow [16][17][18] in using the term entanglement density. A similar quantity was studied in [19]. entanglement density provides a natural metric on the latter space as well (see [20] for more details in the AdS 3 /CFT 2 case). It has also been argued in [18] that in holographic theories, equation (5) applied to a special class of deformations maps to the integrated null-energy condition on the bulk minimal-area surface.
In general, entanglement density in conformal field theories can contain two types of terms: (1) contact terms which arise in the coincident limit x a → x b , and (2) a non-local term which is finite and well-defined when x a and x b are separated. For the most part, we will be interested in the latter. This non-local term is isolated via the definition (4) in terms of the conditional mutual information which makes it is clear this term should be independent of the UV cutoff. The main result of the present paper is as follows: for any conformal field theory, the non-local term in the entanglement density for a planar surface is universal and given by S (2) ij,non−local (x a , x b ) = − where C T is the numerical coefficient appearing in the two-point function of stress tensors in the CFT. Equation (6) was obtained in [17] for a class of holographic theories using the Ryu-Takayanagi formula. However, we emphasize that in this paper we are working with completely general CFTs. 2 We will employ purely field theoretic techniques (developed in [6,[21][22][23]) to prove equation (6), thus extending the validity of this formula to arbitrary conformal field theories with or without holographic duals, and further providing a nontrivial check on the Ryu-Takayanagi and the Hubeny-Rangamani-Takayanagi proposals for entanglement entropy in holographic theories.
An analogous formula also holds in the case where we take A to be a ball-shaped region of radius R (see Fig. 2). Let x = (x 1 , · · · , x d−1 ) denote the spatial coordinates on the Cauchy slice x 0 = 0 and take A to be the region x 2 ≤ R 2 . For x a = R Ω a and x b = R Ω b two well separated points on the entangling surface, the non-local term in the entanglement density is given by In fact, since a ball-shaped region can be mapped into a half-space by a conformal tranformation, we will argue that equation (7) follows as a direct consequence of equation (6) in a CFT.
A number of results follow as corollaries: (i) in [24,25], it was conjectured based on holo-2 Actually it is enough to invoke the SO(1, d−1)×SO(1, 1) conformal symmetries that leave the entangling surface ∂A fixed, including the boosts in the transverse plane, to argue that the cut-off independent part of the entanglement density should take the form as in (6). Here we will be tracking down the overall coefficient.
A δA a δA b Figure 2: The set-up for the spherical case. We deform the original entangling surface x 2 = R 2 (dashed line) by glueing on the area elements δA a,b at points Ω a,b along the entangling surface.
graphic and numerical evidence that the coefficient a(θ) of the corner term contribution to the entanglement entropy in d = 3 CFTs has the universal behaviour where θ is the opening angle of the corner. We will prove this conjecture as a special case of our results. (ii) We also prove the Mezei formula for the universal part of the entanglement entropy across deformed spheres which was conjectured in [26] based on holographic calculations in a large class of theories. In (9), a ,m 1 ,··· ,m d−3 are the coefficients of the expansion of the shape deformation in terms of real hyperspherical harmonics on the entangling surface. The Mezei formula is meant to apply to the universal term in CFT entanglement entropy for a deformed sphere, and the positivity of the overall coefficient demonstrates that the sphere locally minimizes this universal term in the space of shapes, suggesting that the sphere is somehow the optimal measure of degrees of freedom in a CFT for use as an RG monotone. Further, the above formula was used in [27] to compute universal corner contributions to entanglement entropy in higher dimensions. Therefore, our proof of the Mezei formula also completes the proof of universality of corner contributions in higher dimensions. In this way, our CFT calculation fits nicely into the triangle of recent studies and conjectures [17,[23][24][25][26][27][28][29][30] (see also [31][32][33][34][35]) on entanglement density, corner contributions and the Mezei formula.
The rest of the paper is organized as follows: in section 2, we review some elementary facts about entanglement across planar and spherical surfaces, which will be relevant for our subsequent calculations. In section 3, we present the CFT calculation of the universal non-local term in the entanglement density for planar and spherical surfaces. In section 4, we will then use our result for the entanglement density to prove the universality of corner contributions in d = 3 CFTs and the Mezei formula. Finally, we will end with some discussion about prospects for future work.

Preliminaries
Entanglement entropy is defined as follows -consider the density matrix |Ψ Ψ| corresponding to a pure state defined on the Cauchy surface Σ. In this paper, we will take |Ψ to be the ground state of a conformal field theory. Let us partition Σ into two subregions A and A. For local quantum field theories, we expect the Hilbert space h Σ to factorize into the tensor product h Σ = h A ⊗ hĀ. If this is the case, we can trace over hĀ to obtain the reduced density matrix ρ 0 = tr hĀ (|Ψ Ψ|) (10) which contains all the relevant information pertaining to the subregion A. Then the entanglement entropy between A andĀ is defined as the von Neumann entropy of ρ 0 In this context, the boundary ∂A of A is referred to as the entangling surface. It is also useful to define the modular Hamiltonian (also known as the entanglement Hamiltonian) In general, the modular Hamiltonian is not a local operator, in the sense that the modular evolution U (s) = ρ is 0 does not map local operators to local operators. However, there are a few special cases where symmetry forces the modular Hamiltonian to be local. The simplest such case is when we take A to be the half-space x 1 > 0. In this case, the modular Hamiltonian takes a simple form in terms of the CFT defined on Euclidean space R d : (13) where K is the generator of rotations around the entangling surface in the (x 0 E , x 1 ) plane (x 0 E is Euclidean time) and the constant in (13) is chosen such that tr h A ρ 0 = 1. This is known as the Bisognano-Wichmann theorem [36]. The fact that the modular Hamiltonian for planar entangling surfaces in the vacuum state of a conformal field theory on R 1,d−1 is local, and can be written as an integral over the stress tensor will play a crucial role in our calculation of the entanglement density. In fact, the statement of the Bisognano-Wichmann theorem is true for the vacuum state of any relativistic quantum field theory, irrespective of conformal symmetry, and so it should be possible to extend our calculation to the more general class of relativistic quantum field theories. However, in this paper we will restrict ourselves to CFTs, because the calculation simplifies greatly in this case.
In conformal field theories, the modular Hamiltonian for a ball-shaped region (of radius R) is also local [4]. This happens because the conformal transformation with C µ = (0, 1 2R , 0, · · · , 0), maps the half-space x 1 > 0 to the ball-shaped region x 2 ≤ R 2 . Since conformal transformations are symmetries in a CFT, such a map leaves the ground state invariant, and the reduced density matrix on the ball-shaped region can be related to the reduced density matrix on the half-space by a unitary transformation. Additionally, one can transplant the modular Hamiltonian from the half-space to the ball-shaped region by pushing forward the modular flow by ψ, which gives H E,sphere = 2π For this reason, the calculation of the entanglement density for ball-shaped regions is no more difficult than the calculation for half-spaces in CFTs.

The CFT computation
Let us now delve into the calculation of the entanglement density in conformal field theories. For simplicity, we will describe the computation for half spaces in some detail, and then derive the corresponding result for ball-shaped regions by using the conformal transformation mentioned previously. So take A to be the half-space x 1 > 0. Consider now the entanglement entropy of the deformed region A + δA given by x 1 > − χ(x i ). In order to compute this entropy, we can use the coordinate transformation to map the deformed entangling region A + δA back to the half-space A. However, we must bear in mind that such a coordinate transformation has a non-trivial action on the metric. In terms of the new coordinates, the metric is given by Therefore, to compute the entanglement entropy for A+δA in flat space, we may equivalently compute the entanglement entropy for the half-space A but with the deformed metric g µν [37] 3 For our purpose it suffices to keep only the term linear in ζ µ in equation (18) because we are interested in computing the non-local contribution to the entanglement entropy at second order in the perturbation series, while the O( 2 ) terms in (18) can at most generate a local contribution at this order. The shape deformation in (17) is somewhat special in that it preserves the Cauchy surface x 0 = 0. In our calculation we will relax this and consider the more general deformation which also includes time-like deformations of the entangling surface, and of which equation (17) is a special case. 4 The advantage of trading the original problem of computing S EE [A + δA, η µν ] with that of computing S EE [A, η µν + 2∂ (µ ζ ν) ], is that it is possible to use conformal perturbation theory to write an expansion for the latter in terms of the deformation δg µν = 2∂ (µ ζ ν) [21][22][23]. To see how this works, consider the reduced density matrix ρ on A in the presence of the metric deformation δg µν . A straightforward calculation shows (see Appendix A for details) are now coordinates on Euclidean space R d , and 3 Note that (19) is true (even for the UV divergent terms) if we use a "covariant" regulator to define EE [38][39][40]. However since we are ultimately interested in a UV finite quantity the regulator used at intermediate stages in the calculation should not matter. 4 We need not include components along the transverse directions ∂ i because these simply amount to reparametrizations of the entangling surface, which do not change the entanglement entropy.
is the original reduced density matrix on A in the absence of the metric perturbation. 5 Further, T is the angular-ordering operator in the (x 0 E , x 1 ) plane, i.e., if θ ∈ [0, 2π) is the angular coordinate in the (x 0 E , x 1 ) plane, then where H(θ a − θ b ) is the Heaviside step function. The next step is to perturbatively expand the entanglement entropy S EE = −tr( ρ ln ρ) using eqaution (21). However, care must be taken in expanding the logarithm, because ρ 0 and δ ρ = ρ − ρ 0 do not commute in general.
In order to deal with this, we use the following integral representation for the entanglement entropy Expanding this out to second order in δ ρ, we obtain Substituting equation (21) in (25), we find that δS EE can be written as a sum of two terms coming respectively from the first and the second term in equation (25). The first term (after performing the β integration) is given by where H E is the modular Hamiltonian corresponding to ρ 0 and tr conn. is the connected trace. From the above equation, we see that δS (1) EE can be interpreted as the change in the expectation value of the (original) modular Hamiltonian; we will henceforth refer to this term as the modular Hamiltonian term. Given that all the operators inside the trace are naturally T -ordered, we can rewrite the above traces in terms of connected Euclidean correlation functions 5 From now on, by tr we mean tr h A unless otherwise specified. Now, the second term in (26) is given by This term is in fact the negative of the relative entropy of ρ with respect to ρ 0 at second order in δg µν δS and will henceforth be referred to as the relative entropy term. The non-negativity of relative entropy then implies δS (2) Unfortunately, the operators appearing in equation (29) are not manifestly T -ordered, and so the trace in this form cannot be written as a Euclidean correlation function. However, it is possible to perform the β integral and manipulate this expression further to bring it to a more convenient form (see Appendix B) where (R(θ)) µ ν is a rotation in the (x 0 E , x 1 ) plane by the angle θ. This manipulation essentially involves steps similar to passing from old-fashioned perturbation theory to time dependent perturbation theory in quantum mechanics. This is usually achieved using Schwinger parameters and the s appearing above can be thought of as such.
The trade-off however is the additional s integral with the attendant measure. Interestingly, note that the way s appears in the above correlation function corresponds to a relative boost between the two stress tensor insertions, with s being the boost angle/rapidity. Equivalently, from the point of view of the modular Hamiltonian, we are forced into "real time" evolution. Indeed, we can rewrite the above equation in the following way to make this point manifest Having written this term in the above form we can now use the T -ordering to rewrite the trace in terms of the Euclidean two-point correlation function to obtain .
From equations (28) and (34) we see that the entanglement density can be computed in terms of the two-point and three-point Euclidean correlation functions of the stress tensor. Indeed, in any conformal field theory, these correlators are universal and fixed by conformal invariance modulo finitely many parameters [41]. The two-point function in particular takes the form and is determined entirely by specifying the single parameter C T . The three-point function is more complicated, and in general dimension depends on three independent parameters. Nevertheless, it is clear from the above discussion that the (non-local part of the) entanglement density in a CFT is uniquely determined in terms of the parameters appearing in the two-and three-point correlators.
All that remains now is to explicitly evaluate the integrals in equations (28) and (34). Doing so, one encounters the following surprising result -the modular Hamiltonian term (28) does not contribute to the non-local part of the entanglement density. Since the explicit computation is somewhat tedious, we will defer the details to Section 3.2. We also give a quicker more sketchy proof of the vanishing of the modular Hamiltonian term, using a slightly different setup, in Appendix E. The non-trivial contribution to the entanglement density then comes entirely from the relative entropy term. Indeed, this is why the result (6) for the non-local part of the entanglement density depends only on the single parameter C T . We now proceed to compute the relative entropy term.

Relative entropy term
In order to compute the integrals in (34), it is much more efficient to use the conformal transformation from H = S 1 × H d−1 to R d to pull-back and evaluate the integrals on H. To see how this works, let us coordinatize S 1 × H d−1 by y α = (τ, z, x i ), where τ is periodic with period 2π, and (z, x i ) are Poincaré coordinates on the hyperbolic space H d−1 . The metric on H in these coordinates is given by The map ϕ : H → R d given by is a conformal transformation, i.e.
with Ω(y) = z being the Weyl factor (and ϕ * being the pullback). This implies that the stress tensors on the two spaces are related by where S µν denotes additional Schwarzian derivative-type terms, which vanish in odd dimensions, but are present in even dimensions. The integral (34) then pulls back to where we have defined where the vector field ξ α on H is the push-forward of the vector field ζ µ on R d by ϕ −1 Note that the Schwarzian terms S µν have dropped out of equation (41) because of the connectedness of the correlation function. Additionally, we note that the second term in (43) can also be dropped by the tracelessness of the stress tensor (more precisely, the trace Ward identity) in conformal field theories. 6 So we obtain Since the above integrals include integration over hyperbolic space, there are potential divergences coming from the conformal boundary of hyperbolic space at z = 0. These divergences in the entanglement entropy of course correspond to the short-range entanglement coming from the region close to the entangling surface. One way to regulate such potential divergences is to put a cut-off at z = 1 Λ (which corresponds to cutting out a tubular neighbourhood around the entangling surface in the original description on Euclidean space). We denote the resulting regulated space as H Λ , and rewrite the above integral as Next, integrating by parts and using the diffeomorphism Ward identity, 7 we arrive at where ∂H Λ is the boundary of the regulated space H Λ at z = 1 Λ , n = 1 Λ ∂ z is the outward pointing unit-normal on the boundary, and dμ is the measure induced on the boundary The next thing to compute is the two-point function of stress tensors on H. This can be done efficiently using the embedding space formalism developed in [42]. In this formalism (see Section 2 of [6] for a review relevant for this calculation), one considers the larger embedding space (or ambient space) R 1,d+1 on which the (Euclidean) conformal group acts linearly. Let us pick global coordinates P A = (P I , P II , P 0 , · · · , P d−1 ) on this space, with the coordinate P I being time-like. One then embeds H (more generally, any space which is conformally equivalent to R d ) as a section of the upper light-cone P 2 = 0, P I > 0. Here, we pick the embedding Now the two-point function of stress tensors in equation (42) can be computed using the embedding space formalism following [42], where Z a,b are auxiliary variables, and the right hand side above is evaluated on the section (49). The index-free function G (2) is defined as, 7 Which says that ∇ a α Π αβγδ (y a , y b ) = 0 and ∇ b γ Π αβγδ (y a , y b ) = 0 for separated points.
Further, the projector P AB αβ is defined as (see equations (82) and (83) for explicit expressions) Using this formalism, we compute the required two-point correlation functions terms in the above expressions do not contribute in the limit Λ → ∞, so we will drop them henceforth. Substituting equations (53)-(55) in (47) and using (42), we obtain where the overlined indices run over i, j = (τ, z), and We have also defined The factor of 1 Λ 2 out front in equation (56) apparently suppresses δS (2) EE in the limit Λ → ∞. However, we must be careful in taking this limit because there is a possible enhancement from the s integration inside M ij . Indeed, naively sending Λ → ∞ inside the integral in equation (58), we see that the s integral diverges as ds e s in the limits s → ±∞. We can extract these divergences by zooming in on the integral in these limits; this gives two contributions The contribution from s → ±∞ can be extracted by changing variables to β = Λ −2 e ±s (61) Finally substituting the above into equation (56) and integrating over τ a , τ b , we get Reverting back to Lorentzian signature, we obtain which is our primary result. It still remains to be shown however, that the modular Hamiltonian term δS EE does not give additional contributions -we show this in the next section.
To end this section we would like to give some insight into the above calculation and in particular where the main contribution to the non-local part of the entanglement density is coming from. In words, the two stress tensor insertions start their lives close to the boundary of the entangling region (a distance 1/Λ from ∂A in the flat space metric.) However when we boost one of these operators by a rapidity of order 2 ln Λ, then the stress tensor gets liberated from ∂A and moves into one of the null generators of ∂D(A), the boundary of the domain of dependence of A -otherwise known as the Rindler horizon. Here the relevant integrated correlation function receives an enhancement of order Λ 2 in such a way that Λ drops out of the final expression. Thus the main contribution to the entanglement density is coming from the correlation function of stress tensors inserted along null generators of ∂D(A). We find this result intriguing and intend to study this further in future works.

The modular Hamiltonian term
Now we return to the modular Hamiltonian term -in particular, the second term in equation (28) (the first term in (28) was studied in ref. [23] where it was shown to have no universal contributions). An alternative, less constructive, proof that this term vanishes is given in Appendix E. Using the map ϕ : H → R d , this term pulls back to δS (1) Further, integrating by parts and using the diffeomorphism and trace Ward identities allows us to rewrite this as (66) In Appendix C, we will argue that the contact terms vanish in the limit Λ → ∞, and so we focus here on the three-point function term. The modular Hamiltonian H E written on where the integral above is on the constant τ = τ c slice. The constant term above drops out of all connected correlators. So the relevant correlation function in the present calculation is the three-point function of stress tensors on H

δS
(1) Once again, it is efficient to use the embedding space formalism to obtain this correlation function (70) The {A m }'s are conformally invariant structures -polynomials made from six basic building blocks [42] H IJ = −2 (Z I · Z J )(P I · P J ) − (Z I · P J )(Z J · P I ) (71) where the triplet (I, J, K) is a cyclic permutation of (a, b, c). The allowed structures are 9 The coefficients α m in (70) are not all independent -imposing the conservation condition on the stress tensors for non-coincident points gives the constraints This fixes two of the coefficients in terms of the rest, leaving three independent coefficients. 10 In general, computing the integrals in (68) over the hyperbolic slice at τ = τ c is a difficult task. However, the following observation makes this computation tractable -the modular Hamiltonian is a conserved charge and so we are free to move it in τ . One therefore expects δS (1) EE to be independent of τ c . One might worry about potential crossing contributions to δS (1) EE when we move the modular Hamiltonian across one of the other stress-tensor insertions, but it can be checked explicitly that these vanish in the limit Λ → ∞ (see Appendix C). Therefore, in the complex w = e iτc plane, δS (1) EE can be extended to an analytic function which is constant along the unit circle |w| = 1, and hence a constant on the entire w plane. We can therefore use this to our advantage by computing δS (1) EE at a special point such as w = 0 or w = ∞. Physically, these two-points correspond to light-cone limits: w → 0 corresponds to writing the modular Hamiltonian as an integral over the past null boundary of the domain of dependence D(A) (or the past Rindler horizon for brevity), while w → ∞ corresponds to writing it as an integral over the future Rindler horizon (see Fig. 3). We take w → ∞ in what follows. 9 In three dimensions there is also potentially a parity odd structure [43,44] that we did not write down.
However it is easy to argue that such a term cannot contribute to the non-local part of the entanglement density based on symmetries and unitarity -there is no parity odd term that we could add to (6) or (7) which preserves the appropriate conformal symmetries and the strong subaddativity constraint. 10 Imposing the conservation equation in the coincident limit fixes a linear combination of these three coefficients in terms of C T [41]. The computation simplifies dramatically in this limit. To see this is more detail, define the points E ± = (0, 1, ±i, 0, · · · , 0) , E 0 = (1, 0, 0, −1, 0 · · · , 0) in the embedding space. These points have an interesting physical interpretation -(the rays corresponding to) E ± form the future and past tips of the light cone comprising the boundary of the domain of dependence of A, while (the ray corresponding to) E 0 constitutes the point at spacelike infinity, or equivalently the point at infinity in the Poincaré coordinates of hyperbolic space H d−1 . We note the relations The relevant projectors in (69) (with z a,b = 1 Λ ) can be written in terms of these points as 11 The function G (3) in (70), by construction, satisfies the transversality conditions P a,b,c · ∂ ∂Z a,b,c G (3) = 0, and we have used these to simplify the projectors.
where P ± a,b = P a,b · E ± . The only other ingredient required to compute δS (1) EE is the following generic integral where dY c is the appropriate integration measure over H d−1 . Precisely in the limit w → ∞, this integral simplifies greatly and can be written in terms of a single integral (see Appendix D) (86) where Y a,b are the embedding space coordinates for H d−1 , and the constant k is given by Putting everything together, one finds where the matrix N ij can be explicitly computed as a series expansion in Λ 2 . For instance, (89) Now comes the surprising part: the terms which could potentially survive in the Λ → ∞ limit come multiplied by the function C 1 (α m ) defined in (78), which vanishes by the conservation constraints. The same is true for all the components of the matrix N -the potentially non-trivial terms in the Λ → ∞ limit are all proportional to linear combinations of C 1 (α m ) and C 2 (α m ). Therefore, lim This completes our proof that the modular Hamiltonian term does not give additional contributions to the non-local part of the entanglement density.

Spherical case
So far, we have presented the calculation for the non-local part of the entanglement density in the case of planar entangling surfaces. It is possible to repeat the above calculation for spherical entangling surfaces, but here we will obtain the corresponding result more directly by making use of the conformal transformation ψ : R 1,d−1 → R 1,d−1 given by where we have used X µ = (X 0 , X 1 , X i ) as coordinates on the domain of ψ, and C = 0, 1 2R , 0, · · · , 0 . This is a conformal transformation because ψ * η = ω 2 η, with the conformal factor ω given by If we use global coordinates x µ = (x 0 , x) to cover the image of this map, then it is a simple matter to check the following statements: (i) ψ maps the Cauchy surface X 0 = 0 to the Cauchy surface x 0 = 0, (ii) if A is the half-space X 1 ≥ 0 on the Cauchy surface X 0 = 0, then B = ψ(A) is the ball-shaped region x 2 ≤ R 2 on the Cauchy surface x 0 = 0, (iii) ψ maps the domain of dependence of A to the domain of dependence of B. Consequently, we can compute the entanglement density for the ball-shaped region B by pushing forward the deformation vector field ζ ∂B (Ω i ) (where Ω i are coordinates on the sphere ∂B) by ψ −1 and then computing the corresponding entanglement density for the half-space A (94) Since the map ψ is a conformal transformation, the Jacobian factor in (93) can be written as where R is a rotation. Since ζ ∂B lies in the plane perpendicular to the entangling surface x 2 = R 2 , it follows that ζ ∂A also lies in the plane perpendicular to the surface X 1 = 0. Further, by rotation symmetry along ∂B (or equivalently, translation symmetry along its inverse image ∂A), we deduce Finally, using the relations we obtain which is the result for the entanglement density for spheres.

Applications
In this section we will present some applications of our formula for the entanglement density. In section 4.1, we will prove the conjectured universality of the corner term contribution to entanglement entropy in d = 3 [24,25]. In section 4.2, we will prove the Mezei formula for the shape dependence of entanglement entropy across deformed spheres, which was conjectured based on holographic calculations in [26]. In [27], the Mezei formula was used to deduce further universality results for corner terms in higher dimensions. Our proof of the Mezei formula thus also establishes the universality of corner terms in higher dimensions.

Corner terms in d = 3
The entanglement entropy of a general subregion in the vacuum state of a d = 3 CFT takes the general form where δ is a short-distance cutoff and is a length scale associated with the size of the subregion. The first term above is the area-law term, while the second term, which is universal, only appears in cases when the subregion has a sharp corner with opening angle θ, henceforth referred to as the corner term. It has been conjectured based on holographic, free-field and numerical calculations, that in the smooth limit θ → π, the corner term in any d = 3 CFT behaves as where C T , once again, is the coefficient appearing in the two-point function of stress tensors in that CFT. Here, we will show that our formula for the entanglement density directly reproduces equation (101).
We start with our formula for the planar case in d = 3: and consider the special shape deformation: which has two sharp corners at x = ±L with opening angle θ = (π −α). The two corners will both lead to independent logarithmic divergences which we can then isolate. We choose this form because it does not have any IR issues and because it is easy to work with analytically. 12 Note that in order to do the integrals in (102) we are forced to confront UV divergences when the two-points x a and x b come together. This is then related to local contact terms in the entanglement density that we have so far avoided discussing. These terms are also related to the usual UV divergence of EE (that is the area law piece for d = 3 shown in (100).) An efficient way to deal with these contact terms is to use dimensional regularization where the absence of a scale in the regulator means that we will only ever see logarithmic divergences (which would then show up as 1/(d − 3) poles.) Since we do not expect logarithmic divergences in d = 3 in the absence of sharp corners, this is then a good way of isolating the term of interest. To this end, we consider an entangling surface in a d-dimensional CFT with the shape determined by (103) independent of the other d − 3 transverse directions. At second order the change in entanglement entropy is then: where V d−3 = µ d−3 is the volume of the transverse space. This last integral can easily be done (and converges for d < 0) giving: Taking the limit d → 3 we find the desired pole and logarithmic behavior: Since we had two corners with equal opening angles we can infer that: which proves the conjecture of [24,25]. We have thus shown that in d = 3 the logarithmic corner term in the entanglement entropy is entirely captured by the non-local part of the entanglement density. Presumably similar remarks/proofs hold for higher dimensional cones as conjectured recently in [27] using the Mezei formula.

Mezei formula
Next we turn to proving the Mezei formula [26] for the entanglement entropy across deformed spheres at second order in the shape deformation. Once again in this section, we will only be interested in spatial deformations. Let us expand the shape deformation χ(Ω) in terms of real hyperspherical harmonics on the entangling surface S d−2 Based on holographic calculations in a large class of models, it was conjectured in [26] that the universal contribution to the entanglement entropy at second order in the deformation is given by (109) The positivity of the overall coefficient implies that the sphere is a local minimum for (the universal part of) entanglement entropy across all shapes with the same topology. Having derived the entanglement density for spheres from purely CFT considerations, we are now in a position to prove this conjecture. We start with our expression for the sphere entanglement density (setting R = 1 for convenience) Substituting equation (108) , we obtain From rotation invariance, it is evident that the integral above takes the form The constant on the right hand side can in turn be written as where H is the space of all harmonics of order . In order to explicitly compute c( ), we can use the higher-dimensional analog of the addition theorem for spherical harmonics 13 where C (n) (x) is the Gegenbauer polynomial. This allows us to perform all but one of the integrals in (113) to obtain Note that we have also introduced the regulator z above to control the divergences which arise in the θ → 0 (coincident) limit. Fortunately, there exists a closed form expression for the above θ integral [45,46] π 0 dθ (sin θ) D−1 C (117) Using these expressions, we obtain 14 13 We have normalized the hyperspherical harmonics as S d−2 d d−2 Ω Y ,m1,··· ,m d−3 (Ω) 2 = 1.
14 The careful reader might observe that the functional form (in z) appearing above is very closely related to the (deformed) Ryu-Takayanagi surface, with sin Θ = 1 z playing the role of the bulk coordinate defined in [26]. This motivates the identification = z 2 −1 All that remains to be done is to take the limit z → 1 + . Let us first consider d odd; in this case the hypergeometric function behaves as Going back to (118), we see that c( ) is divergent in the limit z → 1 + . However, these divergences, as before, are associated with the coincident limit Ω a → Ω b . A proper treatment of these divergences would require knowledge of contact terms in the entanglement density, which we are not in control of. However, we can extract the universal (cutoff independent) term in (118), which comes from the d/2 term in the expansion of the hypergeometric function close to z = 1, where the corresponding coeffcient b 0 is given by For d even, the expansion of the hypergeometric function contains a logarithmic term d/2 ln which then gives the universal contribution to the entanglement entropy, and whose coefficient is given byb Finally, putting everything together and simplifying gives S (2) 2 ln R δ d even (122) which is precisely the formula conjectured in [26].
As mentioned previously, in [27] the Mezei conjecture was used to compute the universal corner term contributions to entanglement entropy in higher dimensions. Since we have now explicitly proved the Mezei conjecture, this also completes the derivation of the higher dimensional corner terms in [27].

Discussion
We have presented a proof of the universality of the non-local part of entanglement density in any CFT in any dimension. The form of the entanglement density is fixed by conformal invariance and the overall coefficient is determined by C T . We have also shown that this universality fits into a triangle of results that have been the focus of recent studies/conjectures on the shape dependence of CFT entanglement entropy and that we summarize in Fig. 4. We have a good understanding of the mechanism behind this universality. The calculation Entanglement density Mezei formula Corner contributions presented here and previous works [6,21] studied entanglement essentially in conformal perturbation theory, writing the answer in terms of n point correlation functions on flat space up to some order n. Conformal invariance fixes the low point correlation functions that go into the calculation. As we push these calculations to higher order in the expansion parameter, four-point functions and higher will appear and we expect universality to break down. At this point one might impose more restrictive conditions on the CFT that are expected of a theory with a gravitational dual -large-N factorization and an appropriate sparseness condition on the low-lying spectrum of operator dimensions [47] -after which we would expect universality to re-emerge. Ryu-Takayangi taught us the surprising result that all CFTs with classical gravity duals have the same vacuum entanglement structure. We are far from being able to prove such a statement 15 , but these calculations represent a first step.
It is interesting to note that even before we reach the stage where universality breaks down, entanglement in this perturbative framework already displays rich features that are expected of a CFT with a holographic dual [6]. The relative entropy contribution studied in Section 3.1 arises from an integral over an operator located inside the domain of dependence of A, which reminds us of smearing functions that are used to construct local bulk fields [49]. It is this feature of the CFT calculation that we think is responsible for probing deep into the bulk of an emergent AdS dual to measure the change in area of the minimal surface.
We now discuss some possible generalizations that one may pursue. Firstly it would be of interest to study the Renyi generalization of entanglement density (see [50] for related discussion). For example using the mutual information definition of entanglement density (4), we can imagine simply generalizing this by taking I → I n where I n is called the Renyi mutual information. This captures the non-local part of the Renyi entanglement density which is then a UV finite quantity. It is not hard to argue, based on conformal invariance alone, that Renyi entanglement density should take the same form as regular entanglement density: S for some unfixed coefficient e n depending on the Renyi index n. As a function of n we expect that e n is non-universal and theory dependent. We do know that as n → 1 it should equal e 1 = 2π 2 C T /(d + 1). In d = 3 we can use this result to make a prediction for the logarithmic term in the Renyi entropies in the presence of a corner. Following the same steps as in Section 4.1 we find the coefficient of the log is a n (θ) = e (d=3) n 12 (π − θ) 2 + . . .
for opening angles close to π. Quite a bit is known about such contributions to the corner terms in Renyi entropies which we may then use to make predictions about e n . For example the conjecture in [28] would lead to the relation: where h n is the (higher dimensional) twist operator dimension which can be related to a one-point function of the stress tensor for the CFT living on the space H d−1 × S 1 (n) where the radius of the circle S 1 has been enlarged by a factor of n relative to the conformally flat version. Thus if we could establish (123) for the entanglement density then it would prove the conjecture of [28]. Our field theory approach applied to this problem would naively suggest e n is related to an integrated (connected) stress tensor two-point function for the CFT living on H d−1 × S 1 (n) and it is not at all clear (without putting too much thought into it) how this could be related to a one-point function on the same space.
Further generalizations, that we hope to pursue in the future, include the computation of higher order terms in the perturbative expansion of the shape deformation as well as studying entanglement density in relativistic non-conformal theories.
Finally we comment on previous studies of second order shape deformations of EE in [22] using similar CFT techniques (see also [51] for the Renyi entropy case). These authors considered more general metric deformations that contain the shape deformation as a special case and attempt to access the universal logarithmic divergences in EE for d = 4 -first written down using different arguments in [52]. This should be contrasted with our approach of examining the non-local (finite) shape dependent part of EE. Certain issues with this CFT perturbative approach were identified in [22], which we suspect would be resolved by a more careful analysis of the relative entropy term along the lines in this paper. However it is not clear that these universal ln terms can be extracted using a non-local finite term in EE, say deformed by a more general metric, which was the origin of the many simplifications that occurred in our calculation.

A Perturbative change in the reduced density matrix
In this appendix, we prove the formula (21) for the perturbative expansion of the reduced density matrix on the subregion A in terms of the metric perturbation δg µν . Consider a generic relativistic quantum field theory that admits a path-integral description in terms of the field φ, which collectively denotes all the fields over which we integrate. The density matrix corresponding to the ground state wavefunction on the Cauchy surface Σ is given by where H is the Hamiltonian, and h Σ is the entire Hilbert space. In the path integral language, we can describe a matrix element of this density matrix as a product of path integrals over the regions x 0 E > 0 and x 0 E < 0 of Euclidean space R d with the appropriate boundary conditions where we have denoted the spatial coordinates collectively as x = (x 1 , x i ), and Z is the partition function of the theory on R d .
Let us now denote the reduced density matrix for the half space A = {x A = (x 1 , x i )|x 1 > 0} by ρ 0 . The matrix element φ A − | ρ 0 |φ A + is then given by gluing the above path integrals along the complementary spaceĀ at x 0 where x A are spatial coordinates on A, and φ A ± denotes a field configuration restricted to A. By slicing this path integral along the angular direction θ in the (x 0 E , x 1 ) plane, it becomes immediately clear that the reduced density matrix can be written in operator form as [36] where h A is the Hilbert space on A, and K is the generator of θ-rotations From equation (A.4), we see that up to an overall shift coming from the normalization, the entanglement Hamiltonian in this case is given by Next, we turn on a small (background) metric deformation δg µν . The new reduced density matrix ρ is given by where we have introduced x µ = (x 0 E , x) to collectively denote the coordinates on R d . One might worry about an extra term in the exponential coming from a change in the stress tensor upon introducing δg. However, at second order in the δg-expansion, such a term can at best give a local contribution, and so we drop it. Therefore to quadratic order in δg, we have From equation (A.4), we can then infer the following operator expression for the change in the reduced density matrix Switching to polar coordinates (r, θ) in the (x E 0 , x 1 ) plane, the operator T (θ, r, x i ) above is to be interpreted (from the point of view of the reduced density matrix) as where R(θ) is the appropriate rotation matrix in the vector representation. Further, T is the angular-ordering operator in the (x 0 E , x 1 ) plane where H(θ a − θ b ) is the Heaviside step function.

B Angular ordering in the relative entropy term
Recall that the relative entropy term is given by (B.1) where we have used the short notation Unfortunately, the above expression is not T -ordered, and cannot be written in terms of a Euclidean correlation function. To resolve this problem, we will perform the β-integral, and then manipulate the expression further to bring it into a T -ordered form. Let's begin by rewriting it as δS (2) Let us denote the first line of (B.3) as − dμ µνλσ for convenience. Using the eigenstates of K defined by K|ω = ω|ω to carry out the above trace, and writing ρ 0 = ce −2π K , we get The β integral can be performed to obtain where ν = 2π(ω a − ω b ). Next, using the formulae allows us to revert back from the spectral representation to the operator-trace form. To proceed, let's split the integral in eq (B.5) into two parts: θ a > θ b and θ b > θ a . For the first integral, we use equation (B.8) and for the second integral we use (B.7)

δS
(2) Now making the replacements x a ↔ x b and s → −s in (B.10), and adding (B.9) and (B.10), we obtain Finally, once again making the replacements x a ↔ x b and s → −s in the above integral, and adding to itself, we obtain δS (2) where the contour is given by C = R + iε sign(θ a − θ b ). This then gives the result (32) used in the main text.

C Contact & Crossing terms
In this appendix, we analyse (i) the contact terms in equation (66), (ii) crossing terms in the three-point function term in (66).

C.1 Contact terms
The contact terms in (66) are given by where we have defined Ξ = ξ α ∂ α ln Ω 2 . The modular Hamiltonian can be written as the following integral over the constant τ = τ c slice We will need to use the diffeomorphism and trace Ward identities for three-point functions of stress tensors, which can be found in [41] (for Euclidean space). On a general manifold with the metric g µν , these Ward identities take the following form Using these identities, we see that most of the contact terms in (C.1) drop out trivially because y a and y b are well separated (i.e., because we are only keeping terms which contribute to the non-local part of the entanglement density). The only potentially non-trivial terms are Before proceeding to consider these terms, we first pause to observe that both the above terms are independent of the time coordinate τ c at which the modular Hamiltonian is placed. This is because H E is a conserved charge, and we can freely move it in τ c as long as we don't cross other operators. When we do in fact cross another operator, say T αβ (y a ), then using the commutator H E , T αβ (y a ) ∼ ∂ τa T αβ (y a ), we generate extra terms involving the two-point functions of stress tensors. However, using the Ward identities for the 2-point functions, it is straightforward to check that such crossing terms in (C.5) vanish for y a and y b well-separated. Thus, we conclude that both the terms in (C.5) are independent of τ c .
In light of the above discussion, we can simplify the integrals in (C.5) by integrating over τ c (and dividing by 2π). For instance, the first term in (C.5) upon using the diffeomorphism Ward identity and integrating over τ c gives where in the second line we have performed the y c integration, and once again we have taken y a and y b to be well-separated. The two-point function appearing above can be computed efficiently using the embedding space formalism. Having done so, one finds that the above term is suppressed by a factor of 1 Λ . The only thing to check is whether the z a integral inside dµ a is divergent, because such divergences could give potential enhancements. Happily, one finds that the z a integral is finite, and thus the above term vanishes in the limit Λ → ∞. Similarly, one can check that the second term in (C.5) also vanishes as Λ → ∞.

C.2 Crossing terms
Next, we argue that the three-point function term in (66) is independent of the time τ c at which we place the modular Hamiltonian. Since the modular Hamiltonian is a conserved charge, we are indeed free to move it around in τ , as long as we don't cross another operator insertion. However, when we do cross another operator, we pick up an extra contact (or commutator) term, which we will refer to as a crossing term. For instance, let us take τ c from τ a − to τ a + ; in this case we pick up the crossing term = 1 8 ∂H Λ dμ a n α (y a )ξ β (y a ) (C.10) Similar to our previous discussion, the two-point function appearing above can be computed using the embedding space formalism. Having done so, one finds that the above term is suppressed by a factor of 1 Λ 2 . There are no other enhancements to cancel this factor, and the above term simply vanishes in the limit Λ → ∞. Therefore in this limit, the crossing terms can be ignored.

D Integral
In this section, we wish to evaluate the generic integral I(n + , n − , n 0 |m a , m b ) = which appears in the calculation of the modular Hamiltonian term, in the limit where we send the modular Hamiltonian to the Rindler horizon. Using Schwinger parameters, we can rewrite this integral as We will use embedding space coordinates Analytically continuing the above integral in the w = e iτc plane and sending τ c → ∓i∞, we find where β c = e |τc| . Now we partition n 0 into two integers α + β = n 0 , and rewrite the above integral as where in the last line we have rotated Y c to align W = t a Y a + t b Y b with |W |(1, 0, · · · , 0). We have also defined P ± a,b = P a,b · E ± . (D.4) Then, with the change of variables z = z /|W |, x = x /|W |, the integral becomes (dropping the primes) Let us focus on I for the moment. We now change the order of the (t a , t b ) and (z, x i ) integrals, and rescale t a,b = √ zt a,b (D.6) Finally, performing the x i integrals, and redefining z = t 2 c , we get where and the matrix A is given by Rescaling t c → β −1 c t c , we obtain (D.9) We can perform the t c integral in the limit β c → ∞. The integral converges for τ a , τ b in a neighborhood of π, and we can then continue the expression outside this region: (D.10) Finally switching to new integration variables t a = σλ, t b = σ λ (D.11) and performing the σ integration, we obtain (D.12) So the full integral becomes .
Using the fact that E 0 · E 0 = 0, we can further simplify this to obtain the final expression used in the main text.

E Quicker argument for the vanishing of the modular Hamiltonian term
We would like to give a quick argument that the second line of (28) vanishes. We will do this without passing to the hyperbolic coordinates as in the main text. We will also not make any attempt to explicitly calculate the modular Hamiltonian integral. Rather our argument here will be based on scaling symmetry and the operator product expansion in the CFT. We cut off the integrals over the stress tensor close to the entangling surface by cutting out a tubular region around ∂A of radius δ and only integrate over the remaining region U δ . Our goal will be to show that this term vanishes as we remove the cutoff δ → 0. This cutoff is related to the cutoff in hyperbolic space used in Section 3.2 via δ = 1/Λ. So recall that the term of interest is Integrating by parts on x a and x b and using the diffeomorphism and trace ward identities we can rewrite this as: (E.2) In Appendix C, we have already shown that the contact terms vanish in the limit δ = 1/Λ → 0, and so we focus here on the remaining term coming from the boundary of the tubular region.
Of course had we not cut off the integral around the tubular region then we would be done -the diffeomorphism Ward identity would leave us with just the contact terms. However we choose to worry about potential divergences around the entangling surface for several reasons. Firstly such terms are generic in entanglement entropy calculations -although, since we are calculating a finite quantity (the non-local part of the entanglement density), one might expect this to not be an issue. Secondly, when we computed the relative entropy term, an enhancement occured in this region which ruins the naive argument that this term should vanish at least as δ 2 ; here we are checking that such an enhancement does not occur for the modular Hamiltonian term.
At this stage it is convenient to rewrite the boundary term in (E.2) as a bulk integral inside the tubular region U δ δS (1) where again there is a contact term we have dropped based on the analysis in Appendix C. This form is now convenient because we can argue that the leading contribution as δ → 0 is: where we have integrated over the tubular region assuming the relevant integrated threepoint function is a constant over this region. If the term multiplying (πδ 2 ) 2 above can be shown to be finite as δ → 0 then this assumption is true and further we can argue there are no enhancements from this contribution to the modular Hamiltonian term. Unfortunately this is not quite correct; instead, we will use the OPE of two stress tensors to show that there is at most a logarithmic divergence coming from the x 1 integral which we should then cut off at small x 1 ≈ δ close to ∂A. This should only lead to a mild enhancement such that the overall scaling of this term is δ 4 ln δ.
x a x b x c Figure 5: The mild logarithmic enhancement comes from the stress tensor in the modular Hamiltonian coming close to one of the other two stress tensor insertions inside U δ .
The issue comes about for x i → x i a,b and x 1 → 0 where the stress tensor in the modular Hamiltonian comes close to either one of the two other stress tensor insertions (recall that x a is well separated from x b so we need never consider these two operators colliding) -see where σ = (0, x 1 , x i − x i a ) andσ = σ/|σ|. The function A depends on several conformally covariant structures with three unfixed theory-dependent parameters. These are the same parameters that appear in the stress tensor three-point function.
Plugging this into the three-point function we find the leading behavior: This integral is log divergent close to σ = 0 which we cut off by hand at x 1 ∼ δ. We justify this since for x 1 ∼ δ we cannot replace the integral over the tubular region by an integral over ∂A. Then we can argue that C 00;µν;αβ ∼ α ln(δ/|x i a − x j b |) + β, where the log is cut off at long distances when the OPE expansion breaks down. This argument does not fix β. It is also possible to show using the same OPE argument that the log is indeed the only enhancement possible for x 1 ∼ δ and this argument also leads to an explicit and finite expression for β. We leave this as an exercise to the reader.