The Edge of Entanglement: Getting the Boundary Right for Non-Minimally Coupled Scalar Fields

In entanglement computations for a free scalar field with coupling to background curvature, there is a boundary term in the modular Hamiltonian which must be correctly specified in order to get sensible results. We focus here on the entanglement in flat space across a planar interface and (in the case of conformal coupling) other geometries related to this one by Weyl rescaling of the metric. For these"half-space entanglement"computations, we give a new derivation of the boundary term and revisit how it clears up a number of puzzles in the literature, including mass corrections and twist operator dimensions. We also discuss how related boundary terms may show up in other field theories.


Introduction
Measures of quantum entanglement, and especially entanglement entropy, have been an inspirational and unifying theme in theoretical physics in recent years. Invented largely in the quantum information community, they have been successfully applied in a much broader context. That entanglement entropy and black hole entropy both have an area law [1,2] has been a source of inspiration in the quantum gravity community. Certain types of entanglement entropy are believed to order quantum field theories under renormalization group flow [3][4][5][6]. These entanglement measures can also serve as order parameters for certain exotic phase transitions [7][8][9][10]. In our field theory context, we are interested in the entanglement between a spatial region A and its complementĀ. We assume, for now, that the Hilbert space may be factorized with respect to these spatial regions: H = H A ⊗ HĀ. Given this nontrivial assumption, one may construct the reduced density matirx ρ A = trĀ ρ by tracing over the degrees of freedom in the complementary regionĀ, where ρ is the initial density matrix. The Rényi entropies are then moments of the reduced density matrix, while the entanglement entropy itself can alternately be defined as the lim α→1 S α of the analytically continued Rényi entropies or as the von Neumann entropy with respect to the reduced density matrix: The modular Hamiltonian is defined to be the logarithm of the reduced density matrix, and plays a key role in computing entanglement and Rényi entropies as well as other measures of quantum entanglement such as the mutual information and the relative entropy.
In general, H is non-local. One simple and important case where H is believed to be known is in computations of the entanglement across a planar interface x 1 = 0 in flat space coordinatized by x µ , µ = 0, 1, . . . , d − 1. In this case, by the Bisognano-Wichmann theorem [11,12], the modular Hamiltonian is given by the generator of boosts in the direction x 1 perpendicular to the interface. In the Euclidean setting, where x 1 and the Euclidean time coordinate ix 0 form a plane, these boosts are simply rotations about the origin. In terms of the stress tensor T µν (defined by varying the action with respect to an external metric), we then naively should be able to claim that the modular Hamiltonian is (1.4) We use the subscript "cov" for covariant. For conformal field theories, this result has additional implications for geometries related to a planar interface in flat space by Weyl rescaling of the metric. For instance, the modular Hamiltonian for entanglement across a spherical interface in flat space would then also be known (as described for example in refs. [13,14]). The problem with the definition (1.4) is that the charge H can be modified by a boundary term [15][16][17]. More generally, given a conserved current J µ , such that ∂ µ J µ = 0, we are free to introduce an antisymmetric tensor field q µν = −q νµ and add a derivative of it to the current J µ → J µ + ∂ ν q νµ without spoiling the conservation condition. The conserved charge is then modified by the boundary integral d d−2 x q 01 . (With regards to H, q µν could for example be a time component of the Belinfante tensor.) In the case of the conformally coupled scalar field φ, refs. [18,19] proposed seemingly distinct methods for fixing the ambiguity. Ref. [18]'s proposal, which we henceforth adopt in this paper, is to add a boundary term such that H = 2π where T µν is the usual "improved" stress tensor for a conformally coupled scalar field, derived by varying the action with respect to the metric, and ξ = d−2 4(d−1) is the conformal coupling. A nice feature of this expression is that each term is separately Weyl invariant; it is straightforward to construct H for geometries related by Weyl rescaling. In contrast, ref. [19] proposed that one should use the "unimproved" stress tensor, dropping the ∇ 2 φ 2 term from the T 00 component of the "improved" stress tensor. In this case, it is less clear how to apply a Weyl rescaling. Ref. [20] later partially clarified the situation by pointing out that at least in flat space, the ∇ 2 φ 2 contribution could be integrated by parts to yield the boundary term in (1.5).
In section 2, we present two partially independent arguments for the choice of boundary term in (1.5). 1 As suggested by refs. [18,19] and as we further clarify here, the φ 2 boundary term resolves a number of puzzles in the literature concerning various perturbative corrections to entanglement entropy of a conformally coupled scalar. The unifying feature of these puzzles is that the perturbative result depends on correlation functions involving H. If the boundary term in (1.5) is not included, these correlation functions take the wrong value, and the perturbative result will disagree with other independent computational methods. Ref. [18] discovered this boundary term affected thermal corrections to entanglement entropy. Ref. [19] noted the importance of this term for computing the entanglement entropy of a massive scalar perturbatively in the mass. Ref. [19] also noted its importance for computing Rényi entropies in the α → 1 limit, or equivalently the dimension of twist operators perturbatively in an α → 1 limit. (The discrepancy in the α → 1 limit of twist operators was noted but not explained in ref. [22].) In sections 3 and 6, we extend to all dimension d ≥ 3 the perturbative mass and twist operator calculations presented in ref. [19] for the case of a dimension d = 3 scalar. (For d = 2, the conformal coupling vanishes, ξ = 0, and there is no discrepancy to be explained.) Indeed we claim for a massive scalar with an arbitrary non-minimal coupling, ξ = d−2 4(d−1) , the improved modular Hamiltonian (1.5) continues to be the correct one to use in "half-space" entanglement calculations. With a mass, we lose the ability to perform Weyl transformations, but we may still consider the mass dependence of the entanglement across the x 1 = 0 interface.
The operator algebra for a non-minimally coupled scalar field in flat space cannot depend on ξ because the curvature vanishes. Therefore, it seems reasonable to assume that the entanglement entropy cannot depend on ξ [20]. For instance in flat space, standard numerical methods which depend on the operator algebra [2] cannot detect a nonzero value of ξ. We show that the boundary term in (1.5) acts precisely to cancel any naive ξ dependence in the area law term of the entanglement entropy for a massive scalar. In this context, the calculations we present are very similar to those in ref. [20], but our interpretation is somewhat different. The possibility and attendant difficulties of ξ-dependence in entanglement entropy is a twenty year old subject [23][24][25][26] which has still not been successfully resolved [27]. We hope our arguments may finally help put the debate to rest.
The paper is organized as follows. In section 2, we give two independent arguments for the necessity of the boundary term in the modular Hamiltonian (1.5). Both arguments rely on specifying appropriate boundary conditions for the entangling surface. The first argument uses the modular Hamiltonian, while the second uses the path integral on a conical space. The next three sections discuss mass contributions to the entanglement entropy of a nonminimally coupled scalar; we discuss in order the small mass limit, the large mass limit, and the case of arbitrary mass. In section 3, we explain how the boundary term resolves a discrepancy in perturbative calculations of the mass dependence of the entanglement entropy for a conformally coupled scalar. As we are in the business of clearing up puzzles in the calculation of entanglement entropy for non-minimally coupled scalars, we decided to resolve a puzzle [28] that turns out not to be associated with boundary terms. This puzzle concerns the large mass limit in 2+1 dimensions. As we describe in section 4, the resolution of this puzzle has a more mundane origin in the dimensional dependence of the conformal coupling. Third, we discuss the mass dependence of the area law contribution for an arbitrary mass, non-minimally coupled scalar in section 5, arguing that the prescription (1.5) should hold also away from the conformal coupling ξ = d−2 4(d−1) . Moving on, in section 6, we use the boundary term to clear up a discrepancy in the calculation of twist operator dimensions (or equivalently Rényi entropy in the limit α → 1). The Discussion in section 7 is an attempt to apply the lessons learned here about boundary terms to gauge fields and more general quantum field theories. Appendix A contains details of the numerical algorithms we use to check our results. Appendix B has two additional versions of the perturbative mass calculations described in section 3.

Fixing the Boundary Term Ambiguity
The action for a non-minimally coupled, massive scalar field φ in a d-dimensional curved space time M with boundary ∂M is The action is Weyl invariant for the choices m 2 = 0 and .
The quantities R and K are the Ricci scalar curvature and the trace of the extrinsic curvature respectively. We have also introduced the induced metric γ µν on the boundary ∂M. The usual action has been supplemented by a boundary term K φ 2 that preserves the Weyl invariance in the presence of ∂M. The boundary term also improves the behavior of the action under variations with respect to the metric, which yields the usual "improved" stress tensor with a less well-known boundary contribution: Here n µ is an outward pointing unit normal vector to ∂M and K µν is the extrinsic curvature.

Hamiltonian approach
We begin our entanglement entropy computations with the case where M is the Rindler The cut-off > 0 is introduced to control singular behavior in the entanglement entropy in the limit → 0. We would like to compute the entanglement entropy of the half space x 1 > at x 0 = 0. This quantity is controlled by the modular Hamiltonian H, also known as the boost generator in the 01-plane. Starting from the stress tensor and the naive definition (1.4) of the modular Hamiltonian, H should be where the canonical momentum Π = ∂ 0 φ and ∇ 2 = 2 + ∂ 2 0 . In deriving the first line, we have used that K µν − Kγ µν vanishes, because K µν has only one nonzero component on the boundary. We have also assumed that φ is not divergent at x 1 = 0, allowing us to drop a x 1 ∂ x 1 φ 2 term at the boundary. This modular Hamiltonian evolves the theory not in x 0 but in τ = tanh −1 (x 0 /x 1 ). We have multiplied H cov by a factor of 2π associated with the periodicity of the Euclidean time τ E = iτ or correspondingly the inverse temperature. In moving from the first line to the second line, we have integrated by parts twice and again assumed that φ is not divergent at x 1 = 0.
The ambiguity in H can be seen by instead starting in flat space and deriving the modular Hamiltonian using Noether's theorem and boost symmetry. The Noether charge is basically the canonical Hamiltonian multiplied by an extra factor of x 1 : where we need to be careful to evaluate K = 1 before taking the limit → 0. The subscript "can" here means canonical. Clearly the expressions (2.4) and (2.6) differ by a boundary term.
The claim of refs. [18][19][20] can be paraphrased thus -the modular Hamiltonian relevant for an entanglement entropy computation is the one with no boundary term at all, or equivalently eq. (1.5).
One way of justifying this choice is through a consideration of boundary conditions. Hamilton's equations reduce toφ = zΠ , where for ease of notation we have set x 1 = z. The dot means a derivative with respect to the time associated with H, i.e. the modular time τ , not x 0 . In deriving theΠ equation of motion, we have assumed that a boundary term vanishes, where we have allowed for a constant c depending on the strength of the φ 2 boundary term.
Our assumption that φ is finite at the entangling surface sets the z∂ z φ combination to zero. If c = 0, then the value of φ at the entangling surface is fixed. Such a Dirichlet condition is not compatible with the notion of an entanglement entropy computation where the value of the field φ should be free to fluctuate at the entangling surface. The remaining option is to set c = 0 and choose the modular Hamiltonian (2.7). A similar argument was presented in ref. [18] but using a Lagrangian framework and the Weyl rescaled geometries: We would like to elevate these boundary condition considerations to a more general principle: To compute entanglement entropy, we must choose a modular Hamiltonian which has sensible boundary conditions for the quantum fields at the entangling surface. For the case of the scalar, that meant the fields should be free to fluctuate but should also not be divergent at z = 0. If a naively computed modular Hamiltonian, for example H cov or H can above, does not have sensible boundary conditions, one may be able to improve it by adding boundary terms.

Path integral approach
We now present an alternate derivation of the same boundary term, using instead a path integral (or replica trick) approach. The Rényi entropy S α across a planar interface is related to the Euclidean partition function on the conical space C α × R d−2 where C α is an α-sheeted cover of R 2 , branched over the origin: (2.11) The codimension-two conical singularity allows one to consider adding codimension-two boundary terms to the action for a conformally coupled scalar, in particular where c is an as yet undetermined constant and x ⊥ are the transverse coordinates on R d−2 parametrizing the entangling surface. Note that such a term does not spoil the Weyl invariance of the theory. On C α , the Ricci scalar has a delta function distribution at the conical singularity R = 4π(1 − α)δ 2 (x). If we put a metric ds 2 Cα = dr 2 + r 2 dθ 2 on C α where 0 ≤ θ < 2πα, then the distribution can be written as R = −2(1 − 1/α)δ(r)/r. The differential operator in the equation of motion is then (2.13) (There is in principle also a distributional contribution to , proportional to δ(r)r∂ r . However, with our choice of boundary condition, this contribution will vanish.) To carry out the path integral, we look for eigenfunctions of this operator using a separation of variables ansatz, φ(t, x ⊥ , r) ∼ e inθ/α+ik i x i ⊥ ϕ(r). We find (2.14) This ordinary differential equation has been well studied both in the context of cosmic strings and conformal quantum mechanics (see for example ref. [29]). The general solution is Normalizability implies that if c 2 = 0, then n < α. An additional constraint comes from integrating the differential equation near the origin. For the n > 0 case, the simplest solution is to set c 2 = 0 in which case ϕ(0) = 0 and the effect of the delta function vanishes. For the n = 0 case, the delta function produces the constraint lim r→0 rϕ (r) = 2ξ(1 − 1/α)ϕ(0) .

(2.16)
We can try to use this constraint to find a relation between c 1 and c 2 . One finds however that the relation can only be solved at a small cut-off scale r and not in the strict limit r → 0: where C = 2ξ(1 − 1/α). Using this procedure, one will find also a bound state (or tachyon) K 0 (irλ) (λ is pure imaginary) with energy proportional to 1/r when α < 1: In the conformal case m = 0 and ξ = ξ c , the introduction of such a scale breaks conformal invariance. The bound state also implies the theory is not stable. One could imagine truncating the n = 0 modes from the spectrum. But that will have the awkward consequence of forcing the field to vanish at the origin, a sort of Dirichlet condition on an entanglement problem where the origin is not a special point and the field should be free to fluctuate there. We propose a different resolution. We use the freedom to add a codimension-two boundary term to the action to remove the delta function distribution from the Laplacian and restore scale invariance to the solutions: The sign is a bit tricky. Eq. (2.19) correctly cancels the contribution from the tip to the Ricci scalar in the Euclidean signature we use in this section. In Lorentzian signature, δI would have the opposite sign. It seems most natural to define theories with arbitrary values of m and ξ as a smooth deformation of the Weyl invariant case. Thus, we should keep the codimension-two boundary term (2.19), with coefficient ξ, that cancels the delta function contribution from the Ricci scalar. On the cone then, our action reduces to 20) and, using boosts and Noether's theorem, the modular Hamiltonian to (2.7). (Note that the space C α × R d−2 has no codimension-one boundary.) The procedure we are advocating here is equivalent to the so-called Friedrichs extension [29] of the Laplacian on the cone. Such an extension is typically used in heat kernel calculations on the cone. See for example ref. [30]. We are arguing that the Friedrichs extension implicitly involves adding a codimension-two boundary term to the action for a non-minimally coupled scalar.

Conformal Perturbation Theory and the Massive Scalar
In this section, we consider how the mass of a scalar field affects the entanglement entropy in the limit where the mass is kept small. We will use conformal perturbation theory.
Let O(x) be a relevant operator of dimension ∆(≤ d) in d-dimensional QFT. The relevant perturbation by the operator induces an RG flow from the UV fixed point described by the CFT with the action I UV to an IR fixed point. Perturbative studies around the UV fixed point (g O = 0) show entanglement entropy has an expansion in the relevant coupling g O of the form [14,31,32] (see also [33][34][35] for related works) where the first order term is given by the integrated two-point function of the modular Hamiltonian H and the relevant operator If we consider half-space entanglement and were to adopt the naive expression (1.4) for the modular Hamiltonian, the first order term should vanish, s 1 = 0. It would vanish because the two-point function of the stress tensor and a scalar operator is zero in CFT, as noted in this context in ref. [32]. However, as we have already discussed, the modular Hamiltonian is not simply an integral over the stress tensor. It also involves a boundary term (1.5). This boundary term can generically give a nonzero result for s 1 .
If the conformal perturbation theory calculation described above is to give the correct answer, it should agree with an equivalent replica trick approach to the computation. In order to reduce the complexity of the problem, we place the CFT on a compact space to remove the IR divergence. Among many choices of the IR regulator, we employ a cylinder type manifold R × S d−1 whose Euclidean metric is and take the entangling surface Σ dividing S d−1 into two subregions at t = 0 and θ = θ 0 as in [28]. Using the replica trick, the entanglement entropy can be extracted from the partition function Z[M α ] on the α-fold cover M α of the manifold (3.4) branched over the codimension-two hypersurface Σ. In particular, we find Instead of working on the intricate geometry M α , it is easier to make use of a Weyl rescaling to a simpler space S 1 × H d−1 with metric whose relation to the original one (3.4) is given by ds 2 Accordingly, the α-fold cover M α is mapped to the space (3.6) with the α times larger period τ ∼ τ + 2πα. The perturbative expansion of the free energy for the relevant perturbation (3.1) is easily seen to take the form of (3.9) Here S 1 α stands for a circle parametrized by τ of the period 2πα. Formal substitution of (3.9) into the formula (3.5) will calculate the entanglement entropy, the first order term (3.3) being 3 In this form, there is no particular reason for s 1 to vanish. For a concrete example, we consider a free scalar theory (2.1) with mass term m 2 φ 2 /2. The theory becomes conformal when ξ = ξ c and m = 0 and we regard the mass term as a relevant perturbation with g O = m 2 and O = φ 2 /2 of ∆ = d − 2. Indeed, we will see below that for a mass perturbation, s 1 does not vanish and that the replica trick agrees with the modular Hamiltonian approach, provided we include the φ 2 boundary term in H.

Free energy method
We use the replica trick to compute the O(m 2 ) correction to the entanglement entropy of a free massive scalar field in d ≥ 3. We assume that the one-point function O(x) S 1 α ×H d−1 does not depend on the position x as the space has translational invariance S 1 α and the homogeneity on H d−1 . Then one can factorize (3.10) to where V d is the integral of the conformal factor The map follows from the coordinate transformation tanh(t/R) = sin θ0 sin τ cosh u + cos θ0 cos τ , tan θ = sin θ0 sinh u cos τ + cos θ0 cosh u . (3.7) We need the expectation value of the one-point function to evaluate s 1 , which can be read off from the partition function of the theory perturbed by the relevant operator defined on Namely, the one-point function is obtained by taking a derivative of the free where we again used the homogeneity assumption. Putting all together, we end up with (3.14) For the mass deformation of the free scalar, we The partition function of a conformally coupled real massive scalar field on S 1 α × H d−1 is given by [36] where µ s (λ) is the Plancherel measure of the real scalar field on H d−1 of unit radius [37,38] This representation of the free energy leads to The λ integral may be evaluated, making the change of variables u 2 = λ and using the relation u sinh(πu)|Γ(iu)| 2 = π and the first Barnes lemma (see (D.1) of [39] Next we need to evaluate V d which is the integral over the conformal factor (3.8) A similar integral was carried out in appendix C of ref. [28]. The τ integral can be done by contour integration, and after a change of variables s = tanh u, the remaining integral can be put in the form Where it converges, this integral reduces to 5 Putting the pieces together, we find that For d > 3, there is a divergence in the s integral from the region near s = 1. Instead of the dimensional regularization we just employed above, we may regulate the divergence by inserting a cut-off = 1 − s max . In the sphere to hyperbolic space mapping, the angular cut-off is related to the cut-off in s via δθ 2 = 2 sin 2 θ 0 . It is useful to make a table of the results for small d to see what is going on (see table 1). In even dimensions, the coefficient of the log is related to the result from dimensional regularization, while in odd dimensions, the δθ independent term is the result from dimensional regularization.
The result in d = 4 dimensions can be written as an area law contribution where Σ = S 2 is the entangling surface and A Σ = 4πR 2 is its area. This result matches the computation in ref. [40] for a minimally coupled scalar. More generally, as we review in section 5, the area law contribution scales as m d−2 log . Thus only in d = 4 is there an overlap region where we can compare the results.
In appendix B, we describe two additional methods for computing s 1 (3.14). Both methods rely on identifying the right hand side of eq. (3.14) as the expectation value :φ(x) 2 : α on the conical space C α × R d−2 . The first method builds on the method of images employed by Cardy in ref. [41] to compute the mutual information of a conformally coupled scalar in the limit where the two regions are far apart. A key quantity in that paper was φ(x)φ(y) α . Using some results from refs. [18,42], it is straightforward to compute φ(x)φ(y) α in the limit x → y and α → 1. The second method computes :φ(x) 2 : α directly from the path integral, assuming the boundary term (2.19) is present in the action. This second method is remarkably efficient.

Modular Hamiltonian method
We now compare this free energy computation with the modular Hamiltonian method. As noted at the beginning of this section, T µν (x):φ(y) 2 : will vanish, and the contribution to s 1 comes entirely from the boundary term in the modular Hamiltonian (1.5): where we used the conformal map from R × S d−1 to S 1 × H d−1 and Wick's theorem in the second equality. We will now derive this two-point function on S 1 × H d−1 , that appears in the second line. We derive the two-point function on S 1 × H d−1 by using a Weyl transformation. We start with the line element (3.6) on S 1 × H d−1 . This metric is related to flat space with line element where the Weyl factor is The explicit coordinate transformations are The two-point function of a free scalar field of dimension ∆ = (d − 2)/2 on R d takes the usual form where we use a standard normalization, . (3.31) This normalization insures that the Laplacian acting on the two-point function gives a Dirac delta function with unit strength.
The Weyl and coordinate transformations map this two-point function to 6 (3.32) With the two-point function on S 1 × H d−1 in hand, we can now evaluate the expression (3.25) for s 1 . We make use of the symmetry on S 1 × H d−1 to fix the point x ∈ ∂H d−1 = S d−2 to the "north pole", Ω i = δ i1 . This choice significantly simplifies the calculation. Plugging (3.32) into (3.25) with (3.8), we get .
The τ integral can be done by contour integration. The s-integral is related to the u -integral by the change of variables s = tanh(u ). This expression is identical to what we obtained using the free energy, (3.23).
We have also computed these mass corrections on the sphere by conformally mapping them to a conical space rather than hyperbolic space. The calculation is far lengthier, and involves Feynman parameters. At the end of the day, the relation between the UV cut-offs in the two calculations is obscured. However, we can still match the universal terms. We will not include the details of this calculation.  demonstrates that the angular dependence of the change in entropy is well approximated by sin θ 0 , as predicted. Additionally, the linearity indicates that the next higher order correction to the change in entropy scales as m 4 . In panel (c), the x-axis is the inverse of the number of grid points on the sphere. The data points were computed from the y-intercept of fits like in panel (b), but for a variety of different grid sizes. The black dot corresponds to the expected coefficient π 3 /32 of the change in entanglement entropy in the continuum limit. The error in the linear extrapolation is about 4 parts in 10 5 .

Lattice calculation
As a third independent check, we calculate the entanglement entropy of a massive scalar numerically, using a variant of Srednicki's method [2]. Further details of our approach can be found in appendix A. In brief, we begin with the canonical Hamiltonian for a massive scalar on a sphere S d−1 . As the sphere has no boundary, issues about codimension one and two boundary terms will not arise. We decompose this Hamiltonian into spherical harmonics on S d−2 , leaving an overall dependence on the polar angle θ of S d−1 . We then compute discretized two-point functions of the field φ and conjugate momentum Π for each angular momentum mode, as a function of the polar angle θ. Given the Gaussian nature of the ground state wave function, from these two-point functions, restricted to a cap on the sphere with opening angle 2θ 0 , we can reconstruct the eigenvalues of the reduced density matrix and thus compute the entanglement entropy.
Given the results of the previous subsections, there is an obvious issue with UV divergences for all d ≥ 4. Only in d = 3 is the check straightforward. In figure 1a, we see that the corrections to the entanglement entropy are indeed linear in m 2 for small m 2 . A more careful analysis of the data, shown in figure 1b, reveals that the linear term has a sin θ 0 dependence, and that the corrections to the linear scaling are O(m 4 ), in agreement with the expectation from conformal perturbation theory. We can do even better and try to compute the π 3 /32 coefficient of the entanglement change numerically. Computing the entanglement entropy for several different grid sizes and extrapolating to get a continuum limit, we are able to obtain the numerical value π 3 /32 = 0.968946 to about four parts in 10 5 . The results are shown in figure 1c.
Already in d = 4, there is a log dependence on the cut-off which makes extracting precise results from the numerics more involved. Figure 2a indicates that the corrections are still linear in m 2 , as predicted. Figure 2b shows that the corrections additionally obey a sin 2 θ 0 dependence, as predicted by table 1. Having varied the lattice spacing, we also get a rough fit for the coefficient of the log. The answer depends slightly on θ 0 . For θ 0 = π/2, we get the best fit of 0.1659 which is very close to the analytic prediction of 1/6. For θ = 41.4 • , we get the worst fit value of 0.158. 6 The invariant distance on H d−1 is known to be Ωi where Ωi (i = 1, · · · , d − 1) are the embedding coordinates of

Large Mass Expansion
Entanglement entropy of an entangling (closed) curve Σ in a gapped (2 + 1)-dimensional system is expected to have a large gap expansion where m is a gap scale and Σ is the length of Σ [43]. An algorithm to fix the coefficients c Σ 2n+1 for free fields was proposed in [44], which relates c Σ 2n+1 to a conformal anomaly in a (2n + 4)-dimensional theory compactified on T 2n+1 . The authors of ref. [44] focused on the flat space case R 2,1 . Attempting to extend the proposal to the cylinder R × S 2 , ref. [28] found a small discrepancy for conformally coupled scalars. The purpose of this section is to explain and remove the discrepancy. The resolution has no direct connection to the boundary terms that form the main topic of this paper; it is instead related to the dimension dependence of the conformal coupling parameter ξ.
Nevertheless, we will be able to perform an interesting cross-check as a result of the computation here. The value was established for a minimally coupled scalar field in refs. [40,45]. One interpretation of the results here is as evidence that β remains equal to −1/12 in the presence of a non-minimal coupling ξ to curvature, consistent with our discussion of area law contributions in the next section.
We begin by reviewing how c Σ 1 can be determined for a free massive scalar on R 2,1 . We start with a free massless scalar field in 3 + 1 dimensions and compactify a spatial direction to a circle of circumference L. Then the theory reduces to an infinite tower of free massive scalars in 2 + 1 dimensions with masses If the entangling surface in 3 + 1 dimensions is topologically a torus Σ 2 = Σ × S 1 , the entanglement entropy S (3+1) Σ 2 becomes a sum over the entropies S (2+1) Σ (m n ) across the entangling curve Σ for the massive free scalars: Taking the L → ∞ limit, the Kaluza-Klein mass becomes continuous (m n → p) and the entropy turns out to be where we introduced a UV cutoff for the KK mass. On the left hand side, there is a logarithmic divergence s 0 log whose coefficient s 0 is fixed by a conformal anomaly, while the corresponding term arises from the order 1/m term in (4.1). Matching these terms, we find Given an entangling surface in CFT 4 , s 0 can be fixed by Solodukhin's formula [5,46] s Solodukhin where a and c are theory dependent central charges, chosen here for the scalar theory, and normalized such that (a, c) = ( 1 180 , 1 120 ) respectively. χ[Σ 2 ] is the Euler characteristic of Σ 2 , R the Ricci scalar, and R aa ≡ a R µν n a µ n a ν , R abab ≡ a,b R µνρσ n a µ n b ν n a ρ n b σ are projected Riemann tensors by the normal vectors n a µ (a = 1, 2) to Σ 2 . k a µν ≡ γ ρ µ γ σ ν ∇ ρ n a σ is the extrinsic curvature with the induced metric γ µν ≡ g µν − a n a µ n a ν . Applying the formula to the entangling surface Σ that is topologically a torus χ[Σ 2 ] = 0 and reducing the theory on S 1 , one can determine c Σ 1 as a function of the extrinsic curvature of the entangling curve Σ on a plane, which is consistent with a numerical calculation [44].
Superficially, nothing about the above argument seems restricted to R 2,1 . We may try to use Solodukhin's formula (4.7) and the relation (4.6) to determine c Σ 1 for an arbitrary entangling surface Σ in an arbitrary two-dimensional manifold M 2 for a spacetime of the form R 1 × M 2 . Ref. [28] made just such an attempt when M 2 = S 2 . Parametrizing a curve Σ by θ = Θ(φ) on S 2 with line element ds 2 S 2 = R 2 (dθ 2 + sin 2 θdφ 2 ), one finds where the extrinsic curvature κ of Σ is (We used the curvatures R = 2/R 2 , R aa = 1/R 2 and R abab = 0.) To make a numerical check, it is useful to restrict to the simple case Θ(φ) = θ 0 of a caplike entangling region with opening angle 2θ 0 . The coefficient (4.8) for a free scalar simplifies to However, a discrepancy was found between the analytic and numerical calculations in the large mass expansion [28]. The numerics suggests an additional contribution to the coefficient (see figure 3)

Resolution by conformal coupling
The issue is that in fact the conformal coupling cannot be ignored on a general spatial manifold M 2 . The discrepancy comes from the dimension dependence of the conformal coupling. In flat space, the mass contribution from the conformal coupling is zero, and there is no discrepancy. On the sphere, however, the dimension dependence leads to mixing between the O(m) and O(1/m) terms in the expansion (4.1). The essential point is that in the expansion (4.1), the mass gap m is defined with reference to the three-dimensional coupling ξ 3 = 1/8 while in the integral (4.5), the mass p should be defined with respect to the coupling ξ 4 = 1/6. We have in general that Given this relation, we can rewrite two terms in the large m expansion in terms of p: In our setup, R (3+1) = R (2+1) = R(S 2 ) = 2/R 2 , Σ = 2πR sin θ 0 , while β = −1/12 for a free scalar field [40,45]. Assembling the pieces, we find that the second term in the square bracket in (4.13) agrees with (4.11) and resolves the discrepancy found in [28]. That a value of β = −1/12 is required for the discrepancy to be resolved can be viewed as evidence that the value of β is independent of the conformal coupling ξ.
To show more clearly how much the shift term (4.11) improves fitting in the large mass region, we plot the numerical and analytic results for a few choices of θ 0 in figure 3. We introduced the renormalized entanglement entropy of Liu-Mezei type [47] F LM (θ 0 , mR) ≡ (R∂ R − 1)S(θ 0 , mR) , (4.14) to remove the UV divergence. 7 The leading large mass expansions with and without the shift term (4.11) are shown in the orange dashed and green dotted lines respectively. The orange lines fit the numerical data beautifully even in the mR ≈ 2 region.

Alternate resolution
An alternative way to fix the large mass discrepancy is to reduce a free massless scalar with non-conformal coupling ξ 3 in four dimensions to an infinite tower of the KK modes with masses m = m n (4.3) in three dimensions. Namely, we uplift the theory of a free massive scalar with conformal coupling in 2 + 1 dimensions to a free scalar theory in 3 + 1 dimensions 7 This renormalized entanglement entropy does not have a monotonicity property under the RG flow in contrast to the flat space case [4]. Another type of renormalized entanglement entropy FC ≡ (tan θ0 − 1)S(θ0) on the cylinder is proposed and shown to be monotonic numerically for a free massive scalar field in [28].
with the action Note that this action is not conformally invariant as the coefficient of the curvature coupling differs from the correct value ξ 4 = 1/6. A manifold such as a cylinder M = R × S 2 has a constant Ricci scalar, and the conformal mass can be treated on the same footing as a usual mass. As we discussed already around eq. (3.24) and as we will discuss in more detail in section 5, one can show that the logarithmic divergence of the entropy for the theory (4.15) has a contribution additional to Solodukhin's formula: where δs 0 is a universal area term. If one naively regards the action (4.15) as a conformally coupled scalar with the mass (ξ 3 − ξ 4 ) R (3+1) , one would get from eq. (3.24) where A Σ 2 is the area of the entangling surface Σ 2 . In our setup, we find which yields a shift for c Σ 1 given by (4.11) using the relation (4.6). Though the discrepancy has been successfully resolved, we emphasize that the formula (4.17) only works for a manifold with a constant curvature. For a general manifold, ref. [48] argues a formula based on several examples where h is the induced metric on the entangling surface Σ 2 , and χ[Σ 2 ] is the Euler number. If this formula were true, we would find δs 0 = 0 as our entangling surface is topologically a torus and come back to the discrepancy. Instead we propose an alternative formula with the Ricci curvature R(g) of the background metric. Our formula differs from theirs by the extrinsic curvature terms due to the Gauss-Codazzi relation, and thus correctly reproduces the result of ref. [48] for the spherical waveguide geometry. Our proposal (4.20) may be testable by putting the theory on an ellipsoid instead of S 2 to see the dependence on the background curvature.

Universal Area Term
In a field theory with a mass gap m, the entanglement entropy is believed to obey a universal area law for any shape of an entangling surface on flat spacetime [40]: where γ d is a dimension-dependent constant and A Σ is the area of the entangling surface.
Employing the heat kernel method on the replica manifold, the constants γ d were determined for free massive scalar and fermion fields [40] with minimal couplings to the background geometry. There are, however, a one-parameter family of non-minimal couplings ξ Rφ 2 for a scalar field when coupled to a curved space and it is not so obvious whether the universal area term depends on the coupling ξ. Ref. [20] argues the minimal choice to be natural on physical grounds as the operator algebra on flat space does not depend on ξ.
The area law dependence can be derived both from the Adler-Zee formula [49][50][51] and from the Rosenhaus-Smolkin formula [52,53]. With some assumptions, ref. [20] explicitly demonstrates these two methods are equivalent. Depending on the precise treatment of boundary terms, contact terms and the equations of motion, both methods can also produce a dependence of γ d on ξ.
In what follows we revisit and address the ξ dependence of γ d from the perspective of the discussion in section 2.

Adler-Zee formula
In general, the effective action induced by matter coupled to gravity has a derivative expansion of the metric, respecting diffeomorphism invariance: The first few expansion coefficients are completely determined by flat space expectation values of the stress tensor as [49][50][51]54] c 0 = − Θ d , where Θ = T µ µ and the variation of the stress tensor in the second term of c 1 is taken for a conformally flat metric. The expression for c 1 is called the Adler-Zee formula, and interpreted by them as a renormalization of Newton's constant.
A familiar process of evaluating the effective action on a conical deficit around the entangling surface Σ yields the leading area term of the entanglement entropy This is a universal area formula valid for any theory defined on any space, but the coefficient c 1 appears to depend on the type of the stress tensor to be used. For the scalar field, in view of the discussion in section 2.2, we should use the stress tensor derived in flat space from the minimally coupled Euclidean action (2.20). This action of course does not depend on ξ and neither as a result will the area contribution to the entropy.
To compare with the literature, we consider a two parameter space of stress tensors. One parameter is associated with the non-minimal coupling ξ in the improvement term. The other is associated with the possibility of using the equations of motion. For a non-minimally coupled scalar with mass, the trace of the stress tensor is We can further modify the trace using the equations of motion where κ can be arbitrary. The value κ = ξ/ξ c is particularly nice because it removes the φ∇ 2 φ dependence from the trace in general and in the conformally coupled case (ξ = ξ c and m = 0) sets Θ (1) = 0. The integral of the correlator of Θ on flat space is performed with Wick contraction where the final result does not depend on ξ as it multiplies a total derivative term in the third line. For a scalar with a non-minimal coupling to the curvature we find We would naively then obtain the (incorrect) universal area term As we have argued, the correct result should correspond to the choice ξ = 0. Only in this case do we match the d = 3 (4.2) and d = 4 (3.24) area law results mentioned earlier in the paper.
We then further compute how the Adler-Zee formula responds to shifting the stress tensor by the equation of motion. The two-point function shifts by , (5.10) while the contact term is modified to In this prescription, the universal area term is shifted from (5.9) by a κ-dependent term: The calculation above duplicates and agrees with one in ref. [20] for the particular choices κ = 0 and κ = ξ/ξ c . Here, however, we see the general κ dependence of the Adler-Zee formula. We have gone through the exercise of computing the κ-dependence of the Adler-Zee formula because of an alternate method of computation that relies on a spectral decomposition of the two-point function of the trace of the stress tensor: where G(x, µ) is the Green's function of a free scalar field on flat space with mass µ and c (0) (µ) is the spin zero spectral function [55]. The area term involving the two-point function Θ Θ is nicely written in ref. [20] as In this spectral decomposition, for the conformally coupled scalar the equations of motion are typically [55,56] used to reduce the stress tensor to a term proportional to m 2 φ 2 . In our notation, these choices correspond to ξ = ξ c and κ = 1. Moreover, in the spectral computation the contact term proportional to δΘ/δR is ignored. One might worry, given the κ-dependence above, that employing the equations of motion and ignoring the contact term will lead to further complications. However, precisely for the choices κ = ξ/ξ c and κ = 1, the κ-dependent shift (5.12) vanishes and the contact term (5.11) can be ignored, respectively. The result then matches the (wrong) calculation (5.9) with ξ = ξ c .

Rosenhaus-Smolkin formula
Recalling the definition of the modular Hamiltonian H ≡ − log ρ, the entanglement entropy is its expectation value S = H . The variation of the entropy with respect to the coupling g O of a relevant operator O(x) is then given by Rosenhaus-Smolkin in ref. [57]: Note that the expectation value on the right hand side is taken for a state perturbed by the operator. (Comparing this result to our earlier (3.3), one may draw an analogy to the relation between time-independent first order perturbation theory in quantum mechanics and the Feynman-Hellman Theorem.) As in the preceding sections we add the boundary term in the modular Hamiltonian that contributes to the right hand side of (5.15) For a flat entangling surface, the modular Hamiltonian of a non-minimally coupled scalar is improved to be the minimal one by adding the boundary term: Thus the universal area term should not depend on ξ.
We can be more explicit and see how the boundary term δH acts to cancel the ξ dependent contribution of the entropy. For the mass perturbation for a conformally coupled scalar with g O = m 2 and O = φ 2 /2, we find Integration of both sides yields δS = ξ Adding it to (5.9), we obtain the same result as the minimally coupled scalar, i.e. (5.9) with ξ = 0. Interestingly, the boundary term contributes oppositely to the contact term in the Adler-Zee formula without imposing the eom: Thus our prescription to add the boundary term may be restated as removing the contact term.
Assuming for the moment the naive definition (1.4) of H that does not include an extra boundary term, we can relate the result (5.15) directly to the Adler-Zee formula (5.9). Assuming O is the only relevant operator in the theory, the first step is to make use of a Ward identity to replace O(x) H cov with Θ(x) H cov : where we have kept a contact term dropped in [57] but continued to ignore anomaly terms as we work in flat space. The correlator Θ(x) H cov can then be directly related to the trace-trace correlator in the Adler-Zee formula [20]. The contact term can be rewritten as a functional derivative of Θ with respect to R on the conical space. The Rosenhaus-Smolkin formula with the naive H cov (1.4) actually is then equivalent to the Adler-Zee formula with the contact term, yielding the same universal area term (5.9) as before [20]. Leaving a more detailed discussion of the equivalence between Θ Θ and Θ H cov to ref. [20], we discuss instead how to rewrite the contact term.
The variational derivative with respect to g 00 will produce terms proportional to δ(x − y) and ∂ 1 derivatives acting on it. Terms proportional to ∂ 1 δ(x − y) do not appear. The term we care about looks like ∂ 2 1 δ(x − y) and will lead, after a double integration by parts, to a contact term evaluated at the entangling surface y 1 = 0. Such a term can only come from the variation of the Ricci curvature. We know δR = −R µν δg µν + ∇ µ (∇ ν δg µν − g λρ ∇ µ δg λρ ) . (5.22) In flat space, we may rewrite (ignoring the δ(x − y) term) and find where the second inequality follows from the assumption that the stress tensor has a curvature dependence through the coupling R O to the relevant operator. By dimensional analysis, δΘ/δR must scale as g . Integrating this result with respect to the coupling g O , we produce a factor of d−2 d−∆ and match the contact term contribution to (5.4) on the nose. Finally, let us comment on the δ(x − y) term ignored in (5.24). It is independent of the position from the translational invariance, e.g., δΘ(x)/δg 00 (y) ∼ c δ(x − y) with constant c, and contributes to (5.24) an infrared divergent constant. Thus such a divergence cannot be universal as it needs to be regularized and renormalized in some way.

Twist Operator of Spherical Entangling Surface
The partition function on the replica manifold can be regarded as a correlation function of a codimension-two twist operator σ α located on an entangling surface. Their characterization is far from our reach in general though. Some aspects have been investigated [22,58] for a planar twist operator in CFT. The tracelessness and conservation of the stress tensor determine the correlation function T µν σ α to be 8 up to a factor h α which will be called the conformal dimension of the twist operator spanning over the coordinates x a (a = 3, · · · , d) and sitting on the origin y 1 = y 2 = 0. n i (i = 1, 2) is the unit normal vector to the entangling surface. By the conformal map to S 1 × H d−1 the conformal dimension takes the form with the vacuum energy on H d−1 at temperature T = 1/(2πR α) The minus sign is due to the convention for the Euclidean stress tensor and correspondingly the bulk modular Hamiltonian has minus sign, H = −2π H d−1 d d−1 x T τ τ . As we have already discussed, the modular Hamiltonian H is not simply an integral over the energy density T τ τ (1.4) but involves also a boundary term (1.5).
Just as in the case of massive perturbation theory considered in the previous section, we have two methods at our disposal for calculating h α . The first involves the replica trick and the Plancherel measure on H d−1 . The second involves calculating correlation functions of H with other fields, in this case T τ τ .
In fact we will not compute h α but the first few coefficients in a Taylor series expansion of h α near α = 1: The coefficients are closely related to the derivative of the Rényi entropy around α = 1 [59] S α,a = 2πR a + 1 E α,a , (6.5) 8 Our convention for the stress tensor in Euclidean signature is where we used the same notation as (6.4) and E α is the integrated energy The homogeneity of H d−1 allows us to factor out the volume in the right hand side of (6.6). Then we find a simple relation between h α,a and S α,a For completeness we give the inverse relation This is the same relation as the one in [22] up to a = 2 and corrects the expressions for a ≥ 3. Corrected relationships can be found in appendix E of [60] along with further discussion of the h α . With these relations, what we are going to show is equivalent to the extension to general dimensions of ref. [19] where a related discrepancy was resolved examining S α,a for d = 3 and 4 dimensions.

Free energy method
We can evaluate h α defined by (6.2) using the energy density on S 1 × H d−1 which can be obtained by taking a derivative of the free energy (3.15) with respect to 1/T = 2πR α: Here µ s (λ) is the Plancherel measure (3.16). It follows immediately that (6.10) These integrals are very similar to one (3.17) that we did before in computing mass corrections to entanglement entropy. Again, it is useful to make a change of variables λ = x 2 and employ the integral formula (3.18). Before applying this integral formula, we have to replace the hyperbolic trigonometric functions with Γ-functions. The following two relations are useful, the first to evaluate h α,1 and the second to evaluate h α,2 : We obtain then the first two nonzero coefficients of the Taylor series expansion of h α near α = 1: . (6.14) Using a related heat kernel approach, ref. [22] computed h α,1 in general and h α,2 for even d, 4 ≤ d ≤ 14. Our results (6.13) and (6.14) agree with theirs. However, ref. [22] noted a discrepancy when they attempted to compute h α,2 using the modular Hamiltonian because they did not include the boundary term. We will now see in detail how the boundary term fixes the discrepancy.

Modular Hamiltonian method
The boundary term in question amounts to shifting the naive modular Hamiltonian (1.4) by We will repeat the analysis done by [22] and see if the shift resolves the puzzle. We are going to calculate the first and second derivatives of the conformal dimensions at α = 1 where H 0 = H cov of the introduction with the opposite sign and the correlators are in Euclidean signature. The conformal symmetry imposes the constraint for an (local) operator O, leading to δH T τ τ = 0. Thus the first derivative ∂ α h α | α=1 is not affected by the boundary term. Indeed, ref. [22] compute h α,1 using only H 0 and the result agrees with our replica calculation (6.13). However, the second derivative deviates from the naive value To evaluate the three-point functions involving the modular Hamiltonian and the counter term, we find it convenient to put the stress tensor at (τ, u) = (πR, u 1 1) and the others at τ = 0 following [22] where they obtain For a free conformally coupled real scalar, the parameters are given by [61] and we obtain the corresponding conformal dimension (correcting a typo in (3.26) of ref. [22]) which evidently does not agree with the result (6.14) from the previous subsection. We now proceed to evaluate the three-point functions δH H 0 T τ τ and δH δH T τ τ on S 1 × H d−1 by using a Weyl transformation from flat space. The first step is then to calculate the necessary three-point functions in flat space.

The three-point functions for a conformally coupled scalar
It is useful to put the stress tensor in the form ((5.1) of ref. [61]) eliminating the (∂φ) 2 term. In general, (∂ i φ)(∂ j φ) terms are more difficult to deal with because they require point splitting in the analysis. In this form, in considering T τ τ , the (∂ τ φ) 2 term will not contribute to the three point functions because all of the time coordinates of the insertion points have been set to zero. As we have done before, we use the normalization of the two-point function (6.24) We can then write the three-point functions in the form The factor of 8 comes from contracting the φ fields in various equivalent ways. Note we have used the equations of motion to eliminate a φ∂ 2 φ term from the stress tensor (6.23). In general, this substitution could lead to troublesome contact terms in the three-point functions. In the case here, we will always move x 1 far away from x 2 and x 3 , which eliminates most of the contact term ambiguity. There remains a possible issue with the limit x 2 → x 3 in the T τ τ (x 1 )T τ τ (x 2 )φ 2 (x 3 ) correlation function. However, this contact term is proportional to T τ τ (x 1 )φ 2 (x 2 ) , which vanishes in CFT.

Evaluation of δH δH T τ τ
Starting with the second term, we use a conformal map from S 1 × H d−1 to R d with a spherical entangling ball region B of radius R and reduce it to an integral of the three-point function on a flat space where Ω is the Weyl factor (3.28) at t = 0 (6.28) We take a limit of u 1 1 that corresponds to an r 1 R limit, where the three-point function (6.25) simplifies where The spherical symmetry allows us to fix one of the positions r 2 to the north pole on a sphere and factor out the volume, leading to δH δH T τ τ (τ = π, u 1 ) = Evaluation of δH H 0 T τ τ We proceed to evaluate the first term in a similar way In the r 1 R limit, the three-point function (6.26) takes the form where Plugging into (6.32), we find where we utilized the formulae in appendix A of ref. [22].

(6.36)
Adding the correction to the naive result (6.22), we find perfect agreement with the result (6.14) from the free energy computation. Our numerics confirm that the twist operator dimensions (6.13) and (6.14) are correct. See

Discussion
There are a number of entanglement computations in conformal field theory where one naively might hope to take advantage of the universality of two-and three-point functions involving the stress tensor. We considered two important examples in this paper. One was perturbation by a relevant operator O (in our case a mass term). The second was a computation of the Rényi entropies in the limit α → 1. In both cases, we found that the boundary term for H, in other words the fact that H was not simply a spatial integral over T 00 , spoiled this expectation.
In the case of perturbation by a relevant operator, the first order corrections are controlled by H O . Because T µν O = 0 in flat space, one might hope to claim in general that first order corrections due to relevant operators vanish, that entanglement entropy is in some sense stationary at a UV conformal fixed point [32]. To make this more precise in 3d, we can introduce the quantity F = (L∂ L − 1)S(L) for entanglement across a circle of radius L [47]. At a conformal fixed point, F corresponds to the subleading constant term in a large L expansion of the entanglement entropy. Ref. [4], assuming sub-additivity of the entanglement entropy, demonstrated that F decreases monotonically as L increases, like the Zamolodchikov c-function for 2d-CFTs [62]. However, refs. [63,64] later pointed out that, unlike the c-function, F was not stationary at the UV fixed point for a conformally coupled scalar perturbed by a mass term. Here, we see why. The boundary term in the modular Hamiltonian means that even though T µν O is zero, H O is not.
In the case of Rényi entropies, we saw that the first couple of derivatives of S α in a Taylor series expansion near α = 1 were controlled by the two-and three-point functions H T 00 and H H T 00 . Again, since the two-and three-point functions of the stress tensor are unique up to a couple of undetermined constants, one could hope that ∂ α S α | α=1 and ∂ 2 α S α | α=2 have a universal structure. Since the boundary term does not contribute to H T 00 for a conformally coupled scalar, this expectation was born out for ∂ α S α | α=1 . On the other hand, for the next derivative, ∂ 2 α S α | α=2 , it was crucial to include the boundary term. A third example where the boundary term spoils universality is in thermal corrections to entanglement entropy. We did not discuss that case in detail here because the issue was already successfully resolved by one of the authors [18]. The correction is governed by H O O where O is a primary operator that creates the first excited state, which in turn is naively given by a spatial integral over T 00 O O , which has a universal form. However, in this third case as well, the boundary term changes the result.
A burning question is to what extent these boundary terms could be an issue beyond the case of a conformally coupled scalar. We saw already in section 5 that the same boundary term continues to plays a role for a massive scalar with an arbitrary non-minimal coupling. Ref. [20] takes the point of view that these boundary terms are very special, that they only can occur when there is an operator in the theory with dimension exactly ∆ = d − 2 that can be added as a contact term on the entangling surface without spoiling Weyl invariance. Any additional interaction in the scalar case, the argument goes, would shift the dimension of φ 2 away from d − 2 and remove the contact term.
While we are sympathetic to this argument, we are troubled that two of the simplest conformal field theories have issues with boundary terms. The second example is a U(1) Maxwell field in 3 + 1 dimensions. If we perform BRST quantization and choose Feynman gauge, the BRST action can be written (on a manifold M without boundary) as where b and c are the usual anti-commuting ghost fields. If we consider a conical space M = C α × R 2 , then the Ricci tensor has a distributional support at the conical singularity, R ij = 2π(1 − α)δ ij δ 2 (x) where i, j = 1, 2 index the cone C α . Analogous to the scalar case we considered in section 2, one can find eigenfunctions of the vector Laplacian on the cone that involve Bessel functions. To keep the J 0 (λr) eigenfunctions, one has to eliminate the delta function at the origin. The simplest way to do so is to add a codimension-two boundary term to the action Such a boundary term has been discussed before, for example in ref. [30,65]. Like the boundary term for the scalar, this boundary term for the gauge field is also Weyl invariant. It is not gauge invariant, but we already broke gauge invariance by choosing Feynman gauge. Whether the boundary term can be made BRST invariant in a nontrivial way is a subtle question. In fact boundary conditions in general are a subtle question in computing entanglement entropy. The usual boundary conditions chosen in quantizing a gauge field on a manifold with boundary do not appear to be compatible with an entanglement computation [66][67][68][69]. More generally, it seems a boundary term |A n−1 | 2 could appear for free field theories in even dimension d = 2n involving an anti-symmetric n-form F n = dA n−1 with Lagrangian density F ∧ F . In the CFT context, modular Hamiltonians can be modified by a boundary counter term when there exists an operator O d−2 of dimension d−2 that can have a coupling R O d−2 . Moving away from CFT, one no longer needs to ensure that the boundary term is Weyl invariant. Conceivably any relevant perturbation of the form R O might necessitate a boundary counter term. In the case of half-space entanglement, the consequences of such a boundary term can be dealt with simply -the boundary term simply removes the added R O bulk term, leaving one with the original theory.
We have kept the polar angle on the S d−1 explicit, θ = cos −1 u, in order to be able to do entanglement computations. Because S d−1 has no boundary, there is no issue with possible codimension one or two boundary terms in the definition of the full Hamiltonian on this space.
The Rényi entropies can then be constructed from two-point functions of φ l and π l restricted to a cap on the sphere. The essential observation is that these two-point functions should be the same, whether they are computed with the full or the reduced density matrix. In this relatively simple case of a scalar field, it is then possible to reconstruct the eigenvalues of the reduced density matrix from the restricted two-point functions. In technical detail, we define the matrix C (u 1 , u 2 ) 2 ≡ 1 u du φ l (u 1 )φ l (u) π l (u)π l (u 2 ) , (A.5) where the range of C is restricted such that u ≤ u i ≤ 1, i = 1, 2. The Rényi entropy contribution from H l to S α is then In practice, we do not use this formula itself, but instead use the first few coefficients in a Taylor series near α = 1: The zeroth order term is the entanglement entropy contribution, which we use in the numerical check of the mass computations in sections 3 and 4. The first and second order terms are useful in computing ∂ α S α and ∂ 2 α S α near α = 1 in the numerical check of the twist operator dimensions in section 6.
The total Rényi entropy is given by the infinite sum where dim( ) is the number of spherical harmonics l with l d−2 = , e.g. 2 in d = 3, 2 + 1 in d = 4, etc. This infinite sum has to be treated with care. In the numerics, we perform the sum explicitly up to some max where dim( )R α ( ) is of order 10 −7 . Then we fit a dozen values of S with > max to a power law a b , and compute a correction to the finite sum by integrating the power law out to = ∞.
To discretize u, we choose a grid with lattice points at u j = −1 + (j − 1 2 ) , j = 1, . . . , N , and = 2/N . The operator D is discretized using valid at second order in . With this choice, the contributions from the ghost points u 0 and u N +1 vanish.

B.1 Method of images (Cardy's method)
The quantity can alternately be interpreted as a two-point function on the hyperbolic space near α = 1: We have a closely related quantity from [42], namely the two-point function G α,d (2θ) = φ(y)φ(y ) on the conical space C α × R d−2 for particular insertion points. The insertion points y and y are at the origin of R d−2 , at a distance r/2 from the origin of C α , and at an angular separation of 2θ. We will add a factor of Vol(S d−1 )(d−2) to change the normalization conventions of the two-point function to the one considered here. Near α = 1, the Green's function has the form G (α,d) (θ) = G (1,d)  It follows from this recursion relation that (B.8) The corresponding two-point function on S 1 ×H d−1 should be independent of position. Therefore, the conformal transformation should act to cancel the 1/r dependence. Indeed it does. We find from this Green's function that, , (B.9) This result agrees with the earlier free energy computation of this two-point function.

B.2 Contact term from the path integral
There is yet one more way of deriving lim α→1 ∂ α :φ 2 : S 1 ×H d−1 . Without the additional boundary term (2.19), the quantity :φ 2 : S 1 ×H d−1 should vanish at O(1 − α) for the same reasons that T µν :φ 2 : vanishes on the plane. The O(α − 1) correction to :φ 2 : S 1 ×H d−1 will thus come purely from this boundary term.