Entanglement Entropy in Excited States of the Quantum Lifshitz Model

TL;DR

This paper analyzes the entanglement properties of excited states in the quantum Lifshitz model, revealing universal behaviors and introducing entanglement propagator amplitudes that are explicitly calculated and interpreted.

Contribution

It extends the understanding of entanglement entropy from ground states to excited states in the quantum Lifshitz model using geometrical and topological methods.

Findings

Entanglement entropy for excited states expressed via entanglement propagator amplitudes.

EPAs are geometrical probabilities that are universal across lattice and continuum models.

The quantum Lifshitz model allows analytical computation of excited-state entanglement in 2+1 dimensions.

Abstract

We investigate the entanglement properties of an infinite class of excited states in the quantum Lifshitz model (QLM). The presence of a conformal quantum critical point in the QLM makes it unusually tractable for a model above one spatial dimension, enabling the ground state entanglement entropy for an arbitrary domain to be expressed in terms of geometrical and topological quantities. Here we extend this result to excited states and find that the entanglement can be naturally written in terms of quantities which we dub "entanglement propagator amplitudes" (EPAs). EPAs are geometrical probabilities that we explicitly calculate and interpret. A comparison of lattice and continuum results demonstrates that EPAs are universal. This work shows that the QLM is an example of a 2+1d field theory where the universal behavior of excited-state entanglement may be computed analytically.

Equations

$$S_A=F_A+F_B-F_{A\cup B}+\mathcal{O}(1{)}_{\text{topological}}.$$

$$F_A=f_b\left|A\right|+f_{s}L-\frac{c\chi }{6}\ln L+\mathcal{O}(1)$$

$$S_{A}=2f_{s}L-\frac{c}{6}\Big{(}\chi_{A}+\chi_{B}-\chi_{A\cup B}\Big{)}\ln L+{\cal O}(1).$$

$$H=\frac{1}{2}{\int }_M{d}^{2}x\left[{\Pi }^2+{\kappa }^2{\left(\Delta \phi \right)}^2\right].$$

$$Q^{†}(x)=\frac{1}{\sqrt{2}}\left(-\frac{\delta }{\delta \phi (x)}+\kappa \Delta \phi \right),$$

$${\Psi }_0[\phi ]=⟨[\phi ]⟩{\Psi }_0=\frac{1}{\sqrt{Z}}{e}^{-\frac{1}{2}S[\phi ]}\text{ where }Z=\int \mathcal{D}\phi e^{-S[\phi ]},\text{ and }S[\phi ]={\int }_M{d}^{2}x\kappa {\left(\nabla \phi \right)}^2.$$

$$\Delta L_{\lambda }(x)+\lambda L_{\lambda }(x)=0,$$

$$A_{\lambda }^{†}=\frac{1}{\sqrt{\kappa \lambda }}{\int }_M{d}^{2}x{L}_{\lambda }(x)Q^{†}(x).$$

$$H=\sum _{\lambda }\kappa \lambda A_{\lambda }^{†}{A}_{\lambda }.$$

$$⟨[\phi ]⟩(n_{{\lambda }_1},n_{{\lambda }_2},\cdots )=\prod _{\lambda }{\left(Z{2}^{n_{\lambda }}{n}_{\lambda }!\right)}^{-1/2}{H}_{n_{\lambda }}\left(\sqrt{\kappa \lambda }{\phi }^{\lambda }\right)e^{-\frac{1}{2}S[\phi ]}.$$

$${\Psi }_{\lambda }[\phi ]=⟨[\phi ]⟩{\Psi }_{\lambda }=Z^{-1/2}\sqrt{2\kappa \lambda }{\phi }^{\lambda }{e}^{-\frac{1}{2}S[\phi ]}.$$

$$⟨[{\phi }_A^1]⟩{\rho }_A[{\phi }_A^2]=\int \mathcal{D}{\phi }_B⟨[{\phi }_A^1\oplus {\phi }_B]⟩\rho [{\phi }_A^2\oplus {\phi }_B].$$

$$\text{Tr}{\rho }_A^n=\int \prod _{a\in {\mathbb{Z}}_n}\mathcal{D}{\phi }_A^a\mathcal{D}{\phi }_B^a⟨[{\phi }_A^a\oplus {\phi }_B^a]⟩\rho [{\phi }_A^{a+1}\oplus {\phi }_B^a]=\underset{{\mathcal{S}}_n}{\int }\mathcal{D}{\phi }^a⟨[{\phi }_A^a\oplus {\phi }_B^a]⟩\rho [{\phi }_A^{a+1}\oplus {\phi }_B^a],$$

$$\braket{\mathcal{O}(x_{1})\mathcal{O}(x_{2})\cdots}_{\mathcal{S}_{n}}=Z^{-n}\int_{\mathcal{S}_{n}}\mathcal{D}\varphi_{a}\;\Big{(}\mathcal{O}(x_{1})\mathcal{O}(x_{2})\cdots\Big{)}e^{-\sum_{a=1}^{n}S[\varphi^{a}]}.$$

$$\text{Tr}{\rho }_A^n={\left(2\kappa \lambda \right)}^n{⟨{\phi }_{\lambda }^{11}{\phi }_{\lambda }^{21}{\phi }_{\lambda }^{22}\cdots {\phi }_{\lambda }^{nn}{\phi }_{\lambda }^{1n}⟩}_{{\mathcal{S}}_n}$$

$${\phi }^{ab}=\{\begin{array}{ll}{\phi }^a(x)& x\in A\\ {\phi }^b(x)& x\in B\end{array}\text{ and }{\phi }_{\lambda }^{ab}={\int }_M{d}^{2}x{L}_{\lambda }(x){\phi }^{ab}(x).$$

$$\text{Tr}{\rho }_A^n={\left(2\kappa \lambda \right)}^n{\int }_M\prod _{i=1}^{2n}{d}^2{x}_i{L}_{\lambda }(x_i){⟨{\phi }^{11}(x_1){\phi }^{21}(x_2)\cdots ⟩}_{{\mathcal{S}}_n},$$

$$\frac{{⟨{\phi }^{a_1}(x_1){\phi }^{a_2}(x_2)\cdots {\phi }^{a_{2n}}(x_{2n})⟩}_{{\mathcal{S}}_n}}{\text{Tr}{\rho }_A^n(\text{GS})}$$

$$={\left(2\kappa \right)}^{-n}\sum _{\text{“pairs"}}\prod _{i=1}^n\left({\delta }^{a_{c_i}{a}_{d_i}}{G}^D(x_{c_i},x_{d_i})+\frac{1}{n}\left[G^M(x_{c_i},x_{d_i})-G^D(x_{c_i},x_{d_i})\right]\right)+\text{Har}$$

$$\text{“pairs”}=\left\{C=(c_1,c_2,\dots ,c_n),D=(d_1,d_2,\dots ,d_n)\subset {\mathbb{Z}}_{2n}:C\cup D={\mathbb{Z}}_{2n}\right\},$$

$$G^D(x,y)=\{\begin{array}{ll}{G}^A(x,y)& x,y\in A\\ G^B(x,y)& x,y\in B\\ 0& \text{else.}\end{array}$$

$$\text{Tr}{\rho }_A^n(\text{GS})=\mathcal{T}{\left(\frac{Z_A{Z}_B}{Z_M}\right)}^{n-1},$$

$$\frac{\text{Tr}{\rho }_A^n({\Psi }_{\lambda })}{\text{Tr}{\rho }_A^n(\text{GS})}=\sum _{\text{“pairs"}}\prod _{i=1}^n\lambda \int d^2{x}_{c_i}{d}^2{x}_{d_i}{L}_{\lambda }(x_i)L_{\lambda }(y_i)$$

$$\times \left({\delta }_{a_{c_i}{a}_{d_i}}(x_{c_i}{x}_{d_i})G^D(x_{c_i},x_{d_i})+\frac{1}{n}\left[G^M(x_{c_i},x_{d_i})-G^D(x_{c_i},x_{d_i})\right]\right),$$

$$\displaystyle\alpha\$$

$$\displaystyle\beta\$$

$$\displaystyle\gamma\$$

$$\frac{\text{Tr}{\rho }_A^n({\Psi }_{\lambda })}{\text{Tr}{\rho }_A^n(\text{GS})}=\sum _{\begin{array}{c}\text{“pairs"}\\ C,D\end{array}}\prod _{i=1}^n\frac{\gamma }{n}+\{\begin{array}{ll}\alpha & \text{if ci+1≡di(modn) and ci even}\\ \beta & \text{if ci+1≡di(modn) and ci odd}\\ 0& \text{otherwise}\end{array}\}.$$

$${\alpha }_{\lambda \mu }=\sqrt{\lambda \mu }{\int }_M{d}^{2}x{d}^{2}y{L}_{\lambda }(x)G^A(x,y)L_{\mu }(x).$$

$$\frac{\text{Tr}{\rho }_A^1({\Psi }_{\lambda })}{\text{Tr}{\rho }_A^1(\text{GS})}$$

$$\frac{\text{Tr}{\rho }_A^2({\Psi }_{\lambda })}{\text{Tr}{\rho }_A^2(\text{GS})}$$

$$\frac{\text{Tr}{\rho }_A^3({\Psi }_{\lambda })}{\text{Tr}{\rho }_A^3(\text{GS})}$$

$$\frac{\text{Tr}{\rho }_A^4({\Psi }_{\lambda })}{\text{Tr}{\rho }_A^4(\text{GS})}$$

$$+20a{g}^3+b^4+8b^2{g}^2+20b{g}^3+31g^4,$$

Full text

, ,

1 Introduction

The tools of quantum information theory have proven increasingly useful in recent years to understand and characterize many-body quantum systems [1, 2]. A prime example is the use of entanglement entropy to determine if a groundstate wave function is critical or topologically ordered. Remarkably, the entanglement structure of many-body quantum groundstates of systems with local interactions is quite special, and is quite different from that of random states. This result is known as the “area law” [3, 4]: in the ground state of a model, the entanglement between two subsystems is usually proportional to the surface area between them, whereas in highly excited states or for a random state in the Hilbert space it usually becomes proportional to the volume [3, 5, 6]. These are robust results that are not model-dependent, and they also have crucial implications for numerical techniques as they explain the success of variational matrix-product state algorithms such as the density matrix renormalization group [7, 8, 9, 10].

Defining the entanglement entropy requires an ultraviolet cutoff (lattice spacing) and is therefore not completely universal, but the entanglement often contains universal pieces that are characteristic of a model independently of how it is written. Many examples of this are now known. In a 1+1d conformal field theory, for instance, the entanglement entropy scales logarithmically with the size of the entanglement interval times a coefficient proportional to the central charge of the theory [11, 12, 13, 14]. In topological systems, the number of topological degrees of freedom is reflected in a constant term in the entanglement [15, 16]. Even strongly disordered systems can sometimes show universal logarithmic corrections to the area law in one dimension [17, 18]. One can therefore read off many physically interesting quantities and understand a good deal of the structure of the theory by knowing the entanglement.

Although most results are restricted to low-energy equilibrium settings, the scaling of entanglement in highly excited states or after a global quantum quench has attracted a lot of attention in the past few years [19, 20, 21]. In particular, it was realized that entanglement is a key tool to understand how isolated many-body quantum systems self-thermalize (or fail to do so) under their own unitary dynamics [22, 23, 24]. For example, in many-body localized systems, the area law is preserved even in highly excited states, a fact closely related to the non-ergodic behavior of those systems [25, 26]. Results on one dimensional systems also indicate that the entanglement entropy of low-lying excited states may have universal contributions [27, 28, 29, 30, 31].

Unfortunately the promise of using entanglement to understand condensed matter systems is mitigated by the difficulty of computing it. The best “playground” for entanglement is 1+1 dimensional conformal field theories, where a field theory replica technique pioneered by Cardy and Calabrese [13, 14] provides a practical computational framework for the ground states, at finite temperature [13], and for dynamics [21]. These ideas have been extended to free field theories in higher dimensions [32] and to specialized interacting theories using the AdS/CFT correspondence [33, 34]. There are especially few results concerning the entanglement of excited states in dimensions larger than one, both in highly excited states where entanglement should be a prime tool to study questions related to thermalization, and in low-energy excited states that may show interesting universal contributions to the entanglement.

As a step towards this goal, this paper will describe the entanglement properties in excited states of the 2+1d quantum Lifshitz Model (QLM). The QLM arose as a continuum limit of the venerable quantum dimer model [35, 36, 37], which began as model for the resonating-valence-bond theory of high temperature superconductivity with a quantum critical point known as the RK point [35]. (See [38] for a review). The QLM also provides an example of a deconfining critical point [39, 40]. The quantum Lifshitz model is a continuum limit of the dimer model at the RK point via a height function [41, 42]. It is an example of a conformal quantum critical point, whose wavefunctional has 2d conformal invariance, so that techniques of conformal field theory (CFT) are applicable. This special property simplifies the entanglement of the theory and makes analytic progress possible. It was argued in [43] (and subsequently confirmed in [44, 45, 46, 47, 48, 49] with slight corrections) that the entanglement entropy in arbitrary geometries may be expressed in terms of the free energy $F$ of the associated CFT. One can hope that the understanding of entanglement in excited states of the QLM will inform excited-state studies in standard ($z=1$) quantum critical points with conformal invariance in $d+1$ dimensions, similarly to how the corner contributions to entanglement in the QLM can be compared to the corner contributions in the $z=1$ case, where much has been learned even in the absence of a closed form [50].

Let us quickly recall the results of [43]. Suppose the QLM is defined on a 2d spatial manifold $M$, with a partition $A\cup B=M$. Fradkin and one of us showed that the entanglement entropy can be expressed in terms of free energies of the 2d CFT describing the ground state wavefunctional (in most physically relevant cases, a free boson)

[TABLE]

The scaling of the entanglement entropy then follows from a standard results of Cardy and Peschel [51] that states that if region $A$ is finite and simply connected with a smooth boundary, then the free energy scales as

[TABLE]

where $f_{b}$ and $f_{s}$ are non-universal bulk and surface contributions, $\left|A\right|$ is the area, $c$ is the central charge of the CFT ($c=1$ in our case), $L$ is the perimeter of $A$, and $\chi$ is the Euler characteristic of $A$, $\chi=2-2g-b$ where $g$ is the genus and $b$ is the number of boundaries. The $\left|A\right|$ terms cancel in the entropy, yielding a logarithmic correction to the area law determined purely by the geometry and topology:

[TABLE]

Inspired by this concise and general result, this paper explores what may be said about entanglement in excited states of the quantum Lifshitz model. We find that the entanglement of an infinite class of eigenstates is given in terms of geometrical quantities that have a natural probabilistic interpretation, which we dub “entanglement propagator amplitudes” (EPAs) for convenience. Intuitively, the three EPAs (23a)–(23c), can be thought of as the amount of “entanglement scattering” from side $A$ to itself, since $B$ to itself, and between $A$ and $B$.

The remainder of this paper is organized as follows. Section 2 describes the QLM in the wavefunctional picture and computes its excited states. Section 3 introduces the EPAs and shows that arbitrary Rényi entropies are expressed naturally in terms of EPAs through a combinatorial formula, Equation (24). To understand the significance and universality of EPAs, as well as to form an intuitive picture of them, it is necessary to work on the lattice. Section 4 rederives the result (3) on the lattice via a new method employing the renormalization group and provides numerical checks in tetromino geometries. Section 5 extends this to excited states, defining lattice versions of EPAs. In this setting a physical interpretation is natural and a structure of cancelling divergencs is manifest. Section 6 returns to the continuum to compute EPAs in explicit geometries, presented in Equation (62). A comparison of lattice and continuum results shows that EPAs are numerically equal in both settings — strong evidence for universality. The appearance of such universal quantities is evidence that studying entanglement of excited states in higher dimensions reveals novel physical quantities of interest.

2 The Quantum Lifshitz Model

The quantum Lifshitz model (QLM) is a 2+1d field theory which may be defined on a general 2d spatial manifold $M$. Let $\varphi$ be a compact scalar field (i.e. a field whose target is the circle with radius $2\pi R$) with Dirichlet boundary conditions on $M$. Then define [52]

[TABLE]

Here $\Pi(x)$ is the cannonical conjugate momenta, $\Delta$ is the Laplacian on $M$, and $\kappa$ is a real parameter. The value of $\kappa=(8\pi)^{-1}$ is often used to match correlation functions at the RK point of the quantum dimer model. The QLM (with different choices of $R$ and $\kappa$) is also the continuum limit of various interacting generalizations of the dimer model, as well as the quantum six- and eight-vertex models [46].

The QLM is a Gaussian field theory and thus we expect it to be exactly solvable as an infinite collection of harmonic oscillators. It is convenient to use the Schrödinger picture to find the wavefunctionals $\Psi[\varphi]=\braket{[\varphi]}{\Psi}$, which are complex-valued functionals on the space of field configurations, whose eigenfunctionals satisfy $H\Psi[\varphi]=E\Psi[\varphi]$. This treatment expands the ideas of Ardonne et al. [53], and a pedagogical treatment can be found in Fradkin’s textbook [52]. In this picture, the conjugate momentum becomes a functional derivative $\Pi(x)=-i\delta/\delta\varphi(x)$. We can then take the “square root” of the Hamiltonian by defining

[TABLE]

which satisfies $[Q^{\dagger}(x),Q(y)]=-\kappa\Delta\delta(x-y)$. Then the normal-ordered Hamiltonian becomes $H=\int d^{2}x\;Q^{\dagger}(x)Q(x)$.

To find the ground state wavefunctional, note $H$ is positive semi-definite, so any non-trivial eigenfunctional with eigenvalue zero must be the ground state. We thus seek a functional $\Psi_{0}[\varphi]$ in the kernel of $Q$: $Q(x)\Psi_{0}[\varphi]=0$ for all $x$. This is a first-order functional differential equation whose solution is

[TABLE]

This is the action for a compact boson CFT (up to normalization). The QLM is thus an example of a conformal quantum critical point, whose wavefunctional has the form of the free 2d boson (even though the 2+1d theory (4) has dynamical exponent $z=2$ and is not conformally invariant), so one may employ the results of 2d conformal field theory when working with the wavefunctional [53].

The full spectum of (4) now follows from properties of the Laplace operator on an arbitrary conected 2d manifold. The collection of all eigenfunctions (or ‘modes’) $\left\{L_{\lambda}\right\}$ of the Laplace operator on $M$,

[TABLE]

have positive eigenvalues $\lambda>0$ and together form a complete, orthonormal basis for ($L^{2}$) functions on $M$ [54]. For each mode $\lambda$,111Here $\lambda$ or $\mu$ will be used as an index for modes even though finitely many modes may have the same eigenvalue. define a raising operator

[TABLE]

By using the above commutation relation for the $Q$’s, one may show that that, if $\lambda$ and $\mu$ are two eigenvalues, $[A_{\lambda},A^{\dagger}_{\mu}]=\delta_{\lambda\mu}$, $[A_{\lambda}^{\dagger},A_{\mu}^{\dagger}]=[A_{\lambda},A_{\mu}]=0$, and $[H,A_{\lambda}^{\dagger}]=\kappa\lambda A_{\lambda}^{\dagger}$, and thus a general energy eigenstate is indexed by the number of quanta $n_{\lambda}$ in each mode $\lambda$ of the Laplace operator. We may then rewrite the Hamiltonian for the wavefunctionals as

[TABLE]

It is now clear that the wavefunctional of arbitrary excited states is directly analogous to the quantum harmonic oscillator:

[TABLE]

where $H_{n}$ is the $n$th Hermite polynomial and $\varphi^{\lambda}=\int_{M}d^{2}x\;L_{\lambda}(x)\varphi(x)$.

3 Rényi Entropy in Excited States

Now equipped with the excited states of the model, we can look to understand their entanglement. Using a generalization of Wick’s theorem, this section derives a formula, Equation (24), for arbitrary Rényi entropies in excited states.

For definiteness, let us work in the state with one quantum in the $\lambda$th mode, whose wavefunctional is therefore

[TABLE]

This is an artificial restriction to make the calculations easier and the exposition clearer; the method below works for any wavefunctional. However, the combinatorics is prohibitively difficult to make exact statements for arbitrary states.

To set notation, let us say that ${M}=A\cup B$ and let $\partial$ be the border between $A$ and $B$ — see Figure 1. Then we may write a field-configuration as $\varphi=\varphi^{A}\oplus\varphi^{B}$, the sum of fields with support on $A$ and $B$ respectively. The density matrix $\rho=\ket{\Psi_{\lambda}}\bra{\Psi_{\lambda}}$ has a partial trace with matrix elements

[TABLE]

The path integral is over all possible field configurations on the $B$ side of the manifold. Each field (e.g. $\varphi^{2}_{A}$) lives on one sheet on one side (sheet two on side $A$), as shown in Figure 1. To find an arbitrary Rényi entropy $S^{(n)}_{A}=\frac{1}{1-n}\ln\Tr\rho_{A}^{n}$ we must compute

[TABLE]

where $\int_{\mathcal{S}_{n}}$ is the path integral over all field configurations on ${M}$ with the constraint that each of the $\varphi^{a}$’s are equal on $\partial$ [43]. This comes from the fact that a discontinuity of the fields across the boundary gives an infinite contribution to the action, so the fields must agree at the boundary [43]. Henceforth, the product over $a\in\mathbb{Z}_{n}$ will be implicit. One may consider this as a path integral over an $n$-sheeted surface $\mathcal{S}_{n}$ made of $n$ copies of $M$, which are pinned together along $\partial$. This is quite different than the Riemann surfaces that often show up in 1+1d CFT calculations of entanglement [13]; the surface $\mathcal{S}_{n}$ is not a manifold (since it is not homeomorphic to a neighborhood of $\mathbb{R}^{2}$ near $\partial$) and may have the cut anywhere, not just along a spatial axis. Unfortunately, this precludes the possibility of a conformal mapping to a standard geometry.

Define the surface expectation to be

[TABLE]

The problem of finding the Rényi entropies can therefore be rephrased as computing the expectation value

[TABLE]

where

[TABLE]

3.1 The Sides and the Boundary

The task is now to evaluate the correlation functions in (15). By linearity of expectation,

[TABLE]

so all that is necessary is to evaluate spatial $2n$-point functions on $\mathcal{S}_{n}$. A derives Wick’s theorem — suitably modified for this geometry — by applying a series of coordinate transformations first developed in the Appendix of [48] and simplified in [49]. If the $\varphi_{a_{i}}(x_{i})$’s are fields on sheets $a_{i}$, then the result is

[TABLE]

where

[TABLE]

i.e. all possible partitions of $2n$ elements into pairs. The notation $G^{M}$ (resp. $A$, $B$) is the Green’s function for the Laplacian on ${M}$ (resp. $A$, $B$) with Dirichlet boundary conditions, and

[TABLE]

“Har” means this is only true up to an arbitrary harmonic function; this is irrelevant because we are always integrating against eigenfunctions of the Laplacian, so we can remove it with an integration by parts. Lastly, the normalization factor is exactly the ground state Rényi entropy

[TABLE]

a familar result from, e.g. [43, 48, 49], where $Z_{D}=\left(\det\Delta_{D}\right)^{-1/2}$ is the partition function on the domain $D$ [43] and $\mathcal{T}$ represents the ${\cal O}(1)$ topological terms. Our generalized Wick’s Theorem is therefore the normal Wick’s theorem on each side $A$ and $B$ of each sheet, plus a contribution $G^{M}-G^{D}$ which communicates between sides and between sheets, but is suppressed by $1/n$.

3.2 Entanglement Propagator Amplitudes

Combining (17) and (18) then tells us

[TABLE]

where $\delta_{a_{c_{i}}a_{d_{i}}}(x_{c_{i}},x_{d_{i}})$ is unity only if the fields are on the same sheet. For instance, for the pair of fields $\varphi^{11}(x_{1})\varphi^{21}(x_{2})$, it is unity only when $x_{1},x_{2}\in B$. Equation (22) is a polynomial in three quantities, which we dub the “Entanglement Propagator Amplitudes” (EPAs),

[TABLE]

Here $\alpha$ and $\beta$ are intra-sheet terms that have support on $A$ or $B$ respectively. On the other hand, $\gamma$ is an inter-sheet term, connecting $A$ to $B$. These are dimensionless and should be interpreted as probabilities since $0\leq\alpha,\beta,\gamma\leq 1$ and $\alpha+\beta+\gamma=1$.

Applying a dash of combinatorics,

[TABLE]

This allows us to compute arbitrary Rényi entropies for excited states with one quantum in an arbitrary mode. Such a restriction is artificial; the Rényi entropy of an excited state with multiple quanta could be computed with more complex combinatorics and three EPAs for each pair of modes involved. For example if we had two modes $\lambda$ and $\mu$, then the Rényi entropies would involve quantities such as

[TABLE]

Equation (24) is a main result of this work. To leverage (24) to compute entanglement entropy, we must understand both how to compute the EPAs and the combinatorics of how they are combined. This section comments on the combinatorics briefly and we will return to the interpretation and evaluation of the EPAs presently, with a physical interpretation at the end of Section 5.

To find the entanglement entropy, one must analytically continue (24) to arbitrary $n$, in practice by finding a closed form equation for it [13], a non-trivial problem. Explicitly, the first several traces are

[TABLE]

where $a=\alpha+d$, $b=\beta+d$, and $g=\gamma/n$. One can translate the problem of computing the Rényi entropies into a weighted perfect matching problem. Consider the complete graph on $2n$ vertices, with weights $a$ and $b$ alternating around the outside, and weight $g$ for edges through the middle. The $n$th Rényi entropy is the sum of all weighted perfect matchings. (Amusingly, this is equivalent to finding the partition function of a classical dimer problem on a non-planar graph.)

In the special case $\alpha=\beta$, the Rényi entropies are polynomials in a single variable $\gamma$. We conjecture that

[TABLE]

where

[TABLE]

wherein the coefficients $A_{i,j}$ are the diagonals of OEIS sequence A112486 [55] and may be recursively defined by

[TABLE]

We have confirmed Equation (27) explicitly through $n=11$ via computer.

To find the entanglement in general it would be convenient if one of $a$, $b$ or $g$ could be treated as a small expansion parameter. The behavior of (24) to lowest order in $g=\gamma/n$ is particularly simple:

[TABLE]

It is therefore quite tempting to hope that $\gamma$ vanishes, which would make (27) trivial. Superficially, it looks like $\gamma$ might vanish; we may integrate by parts to get $\Delta G^{M}-\Delta G^{D}$, which would give a delta function on $M$ minus delta functions on both $A$ and $B$. It would appear that the inter-sheet EPA therefore vanishes, since the integral kernel is identically zero except on the co-dimension one boundary $\partial$. In this case, $\alpha+\beta=1$ and $\alpha$ would be the integral of the wavefunction over $A$, i.e. the probability the excitation is on the $A$ side, and similarly for $B$. The analytic continuation of (30) would be trivial, and the entanglement entropy would be that of the ground state plus the entropy of an effective qubit for each quanta with density matrix $\rho_{A}^{\text{eff}}=\operatorname{diag}(\alpha,\beta)$. However we will demonstrate below that, contrary to this naive expectation, $\gamma$ is generically non-zero and is not small in any common physical limit. The reasoning of the above paragraph therefore does not apply, though important aspects of it survive, and it is a good first start for intuition.

To see why $\gamma$ is non-zero and what its interpretation is in the continuum involves regularization and subtle arguments. Therefore, we will explore a explicit lattice regularization of the problem, providing a point of view complementary to the field theory picture, but without divergences. Section 4 will reproduce the entanglement entropy in the ground state from the lattice, and Section 5 will examine lattice excited states, showing $\gamma$ is generically non-zero. With this analysis in hand, Section 6 will reveal how this is realized in the continuum. We will conclude that the EPAs $\alpha$, $\beta$, and $\gamma$ are universal and completely geometric quantities which fully determine the entanglement entropy of excited states of the theory. For states with one quantum, it therefore follows that the difference in entanglement with the ground state is $\mathcal{O}(1)$ and universal.

4 Building the Quantum Drum from Scaffolding

To understand and compute the EPAs on the lattice, it is first necessary to reproduce the ground state results in a lattice setting. To that end, this section re-derives (3) with a lattice regularization and confirms it numerically.

Let us make the following assumptions:

The only difference in the entanglement between the compact and non-compact theories is an $\mathcal{O}(1)$ topological term. We saw in Section 4 that this holds for the ground state and A suggests it holds for excited states as well. Therefore we will examine the non-compact lattice theory. 2. 2.

Regularizing the original Hamiltonian (4), is the same as regularizing the emergent free-boson Lagrangian (6).

We will work on the square lattice of length $a$. To make calculations concrete, we adopt the explicit geometry of Figure 2: a rectangle of size $K\times L$ where the $A$ side is $\ell\times L$. Though a specific geometry has been adopted, the Riemann mapping theorem implies these results hold for any geometry where $A$ and $B$ are simply connected. With a variable $\varphi_{pq}\in\mathbb{R}$ at each lattice site, the action becomes (factors of $a$ cancel out)

[TABLE]

Tilde’s will be used to denote the lattice versions of quantities (except for fields) for easy comparison with the continuum equations. Dirichlet boundary conditions are imposed by setting “out-of-bounds” variables such as $\varphi_{0,3}$ to zero in the action.

From the definition of the partial trace, the $n$th Rényi entropy breaks the “bonds” on the boundary between $A$ and $B$ on each sheet and reinstates them with half the weight going to the same sheet and the other half going to the next sheet (see Figure 2). Therefore the lattice Rényi entropies involve the traces111Here the path integral $\int\mathcal{D}\varphi^{c}$ is just shorthand notation for $\prod_{pq}\int_{\mathbb{R}}d\varphi_{pq}$ to show the analogy with the continuum case and not an actual path integral.

[TABLE]

To prevent a plethora of indices, $A$, $B$, $\partial$, $\overline{\partial}$ and $M$ are used as multi-indices to refer to all the lattice sites in those domains, as labelled in Figure (2). So $\varphi_{M}\widetilde{\Delta}^{MM}\varphi_{M}$ should be interpreted as $\sum_{(pq),(rs)\in M}\varphi_{pq}\widetilde{\Delta}^{pqrs}\varphi_{rs}$, etc. Here $\widetilde{\Delta}$ is the discrete Laplacian restricted to the appropriate domain. Due to the coupling between sheets, the Gaussian integrals here are non-trivial to evaluate explicitly. However, by applying a Fourier transform we may decouple the sheets. The action may be rewritten in terms of the circulant matrix

[TABLE]

as

[TABLE]

We can perform a discrete Fourier transform by the unitary change of variables $\psi^{d}=V_{c}^{d}\varphi^{c}$ withthe Vandemonde matrix

[TABLE]

where $\zeta_{n}=e^{-2\pi i/n}$ is an $n$th root of unity. Since $T$ is a circulant matrix, it is diagonalized by this coordinate transformation, giving

[TABLE]

The result of this transformation is that the $n$-sheeted partition function is now decoupled into $n$ single-sheet partition functions where, at the boundary between $A$ and $B$, the strength of the bonds between the two sides is between zero and one. Therefore, we find $\Tr\widetilde{\rho}_{A}^{n}=\prod_{k=0}^{n-1}\frac{\widetilde{Z}_{k}}{\widetilde{Z}_{M}}$ with

[TABLE]

In particular, $\widetilde{Z}_{0}=\widetilde{Z}_{M}$ since the bond between the two sides is “full strength”, and all the rest have weaker-than-normal bonds across the boundary. (See Figure 3.)

We are thus left with the evaluation of the partition function of decoupled free bosons with “weak links” along the entanglement cut $\partial$. This is similar in spirit to the problem considered by Kane and Fisher who argued that a weak potential scatterer (“weak link”) in one dimensional systems of interacting fermions (Luttinger liquids) with repulsive interactions is a relevant perturbation in the renormalization group (RG) sense, so that it effectively cuts the 1D chain at low energy [56]. We emphasize however that a weak link for the fermions does not correspond to a weak link in the (free boson) Luttinger liquid theory, but rather takes the form of a boundary sine-Gordon perturbation after bosonization. Nevertheless, we expect a similar picture to hold in our case: the weak links with strength $0<\lambda=\cos\left(k\frac{2\pi}{n}\right)<1$ should flow under renormalization to conformally invariant boundary conditions, either $\lambda=0$ (decoupled systems) or $\lambda=1$ (‘healed’ uniform system). To determine which way the RG flow goes, we can study the perturbative stability of these two fixed points. From the lattice analysis above, the perturbation to the homogeneous free boson action is

[TABLE]

with $q$ the coordinate along the cut, which in the continuum limit becomes

[TABLE]

with $y$ the coordinate along the entanglement cut. (Recall we are working with a non-compact boson here.) Starting from an homogenous system with $\lambda=1$, we see that weakly weakening links along the cut is a strongly relevant perturbation that acts as a local mass term pinning down the field to $\varphi(x=0,y)=0$ along the cut $\partial$ at low energy. This indicates that the homogeneous fixed point with $\lambda=1$ is unstable, and flows under renormalization to the free boson theory that vanishes (satisfying Dirichlet boundary conditions) along the cut (corresponding to $\lambda=0$).

One can then self-consistently check that this fixed point with $\phi(x=0,y)=0$ is stable, as the perturbation $V$ around this fixed point $V_{3}=\int_{\partial}dy\left(\lambda_{3}(\partial_{x}\varphi(x=0,y))^{2}+\dots\right)$ has scaling dimension $\Delta=2>1$ and is indeed irrelevant as boundary perturbation. For large systems, we may thus treat the $n-1$ sheets with weak links as if they had no link across the boundary, i.e. Dirichlet boundary conditions on $\partial$, consistent with the result (cf. Equation (21)) of Fradkin and Moore

[TABLE]

Applying the replica trick, the entanglement entropy is

[TABLE]

This is, of course, precisely the same as the continuum result (1) up to the topological term. One can again apply the result (2), whose logarithmic term is universal, to conjecture that entropy takes the same form as in the continuum:

[TABLE]

The appearence of the Euler characteristic comes from an application of the Gauss-Bonnet Theorem in the context of a manifold with smooth boundary. In the case of a shape with straight edges and sharp corners, like the rectangle of Figure 2, there is an extra contribution that depends on the angles of the corners [51]

[TABLE]

Conical singularities on the boundary give a slightly different contribution. When writing the entanglement entropy, if the boundary $\partial$ between $A$ and $B$ has a corner of angle $\gamma$, we must include the contribution from both the angle $\gamma$ and $2\pi-\gamma$. Corners on the boundary of $M$ are cancelled out and do not contribute.

This can be explicitly confirmed in special cases with asymptotics on the lattice. David and Duplantier showed [57] that on a square lattice of size $(K-1)\times(L-1)$,

[TABLE]

where $G$ is Catalan’s number, $\eta$ is the Dedekind $\eta$-function, and $\zeta=K/L$. For the rectangular geometry of Figure (2), this implies that the logarithmic term in the entanglement has coefficient $-1/4$, which is consistent with (43) for 4 angles of $\pi/2$. The result (44) was extended by Kenyon to any simply connected planar shape by mapping the problem to counting states of the classical dimer model [58]. The connection is suggestive that perhaps this lattice regularization is closely linked to the quantum dimer model but beyond the scope of this paper.

4.1 Numerical Confirmation

The discrete setting here gives us a chance to confirm (43). To numerically extract the sub-leading logarithmic behavior, it is necessary to peel back the linear term of (42). One way to do this is by considering the difference in entanglement between two regions $A_{1},A_{2}\subset M$ with the same area and the same perimeter, but a different number of corners. A set of shapes with exactly this property are the tetrominos, shown in Figure 4. (This is conceptually similar to methods to detect topological order in a ground state wavefunction [15, 16].) To do this, the entanglement entropy for the tetrominos shown in Figure 4 was calculated by explicitly computing $\widetilde{F}=\frac{1}{2}\ln\det\widetilde{\Delta}$ on the lattice.

This was done for four tetrominos with area area of $4L^{2}$ and a perimeter of $10L$. Each entropy was evaluated via numerical determinants with the shape centered in a $12L\times 12L$ grid for $1\leq L\leq 20$. Differences in the entropies of pairs are shown in Figure 4. These are fit to

[TABLE]

using linear least squares. There are finite size effects at small $L$, so fits only used $L\geq 5$. All fits had $R^{2}>0.99$. The results confirm (42) to precision of $\sim 10^{-4}$. From this demonstration and the above RG argument, we may conclude that the entanglement entropy in the ground state of the QLM model has the same form on the lattice as in the continuum.

[TABLE]

5 Lattice Excited State Entanglement

The previous section showed that the entanglement entropy in the ground state is the same on the lattice and in the continuum, up to non-universal contributions. It is reasonable to expect that this exact correspondence carries over to excited states as well. Indeed this is the case. This section parallels Section 3 to find arbitrary Rényi entropies, but on the lattice. The correspondence is so close that many of the details are repetitive and are given in B. Here, however, all definitions can be made precise and all quantities are divergence-free and finite. This allows explicit numerical computations of the EPAs. We conclude with a physical interpretation of the EPAs.

Consider the excited lattice wavefunctional

[TABLE]

the direct analogue of (10) on the lattice. Here $\widetilde{L}^{\lambda}$ is the eigenvector of the lattice Laplacian on ${M}$ with eigenvalue $\widetilde{\lambda}$.

Unlike the ground state, where the $A$ and $B$ sides contributed independently to the entanglement, here the two sides are explicitly linked. This can be seen already by examining the reduced density matrix. By definition, for the state $\widetilde{\Psi_{\lambda}}$ it is

[TABLE]

where $\varphi^{21}$ takes the values of $\varphi^{2}$ on $A$ and $\varphi^{1}$ on $B$, and similarly for $\varphi^{11}$. As in the previous section, the path integral is a notational shortcut for the integrals of $\varphi_{pq}$ over all lattice sites and not a literal path integral. After computing the Gaussian integrals one finds

[TABLE]

where $\widetilde{G}^{BB}$ is the inverse of the discrete Laplacian on $B$, and as above $\widetilde{\Delta}^{\partial\overline{\partial}}$ is the restriction of the discrete Laplacian to the boundary, which acts as a “tunneling” matrix between the right-most column of $A$ and the left-most column of $B$. The first and last terms of (49) depend purely on the $A$ side or the $B$ side and will become the analogues of $\alpha$ and $\beta$ in the Rényi entropy. The second term involves tunneling from one side to the other and will manifest as the inter-sheet term $\gamma$. It is therefore not sufficient to consider separately the density matrix restricted to the $A$ and $B$ sides, but rather one must consider the behavior on and across the boundary as well. Indeed, the boundary behavior is crucial in what follows. This is in contrast to the behavior in the ground state, where the two sides contributed to the entanglement independently.

The Rényi entropies for excited states on the lattice involve

[TABLE]

It is again sufficient to compute 2-point functions on the $n$-sheeted geometry, then apply Wick’s theorem and integrate against the eigenvectors of the Laplacian. This is carried out explicitly in B. The upshot is one can write all the Rényi entropies of $\widetilde{\Psi}_{\lambda}$ in terms of lattice EPAs

[TABLE]

These are the direct analogues of (23a) – (23c). Since $\widetilde{L}^{\lambda}$ has dimensions of inverse length, the factors of $a^{2}$ is necessary to preserve dimensions.

These may be simplified. By definition, $\widetilde{\Delta}\widetilde{L}^{\lambda}=-\widetilde{\lambda}\widetilde{L}^{\lambda}$ and $\widetilde{\Delta}^{M}\widetilde{G}^{M}=\text{Id}^{M}$. What one would like to do is use the factors of $\widetilde{L}^{\lambda}$ to produce extra Laplacians, to get

[TABLE]

then reduce the Green’s functions to Kronecker deltas. Unfortunately, this is not possible because the action of the Laplacian does not commute with restricting the domain, so $\widetilde{\Delta}^{M}\widetilde{G}^{A}\neq\text{Id}^{A}$ (resp. $B$ and $D$). Nevertheless, one can do this away from the boundary. As previously, let $\widetilde{\Delta}^{\partial\overline{\partial}}$ be the Laplace operator restricted to the links between $A$ and $B$, i.e. the tunneling operator between $\partial$, the sites on $A$ closest to $B$ and $\overline{\partial}$, the sites on $B$ closest to $A$ (see Figure 2). Then

[TABLE]

Combining these and (53),

[TABLE]

As a concrete example, consider the $K\times L$ rectangular geometry where $A$ is an $\ell\times L$ rectangle, as shown in Figure 2. Numerically, $\widetilde{\Delta}^{\partial\overline{\partial}}$ is minus the identity matrix on $L$ elements, which gives

[TABLE]

The lattice form for the EPAs suggests a physical interpretation in terms of scattering processes. We stress that this is an analogy suggested by the form of the equations to promote intuition and not a genuine physical scattering process. For $\widetilde{\gamma}$ we may consider each entry of the matrix in (57) as a scattering process. The upper right entry, $-\delta_{qs}$ represents “tunneling” from the $\ell$th row to the $\ell+1$th row (Figure 5 (a)). Meanwhile, the upper left entry, $\widetilde{G}^{B}_{\ell+1,q,\ell+1,s}$, represents starting on the $\ell$th row, tunneling across, scattering around anywhere in $B$, and tunneling back across the boundary (Figure 5 (b)). Hence, $\widetilde{\gamma}$ is the sum over all scattering from the $A$ to $B$, minus the scattering back — a “current” flowing across the boundary. Alterntively, $\widetilde{\gamma}$ can be interpreted as the scattering on $M$ sourced by $\widetilde{L}^{\lambda}$ minus the scatter strictly within $A$ and within $B$, which leaves the scattering across $\partial$. Though similar logic, one may interpret $\widetilde{\alpha}$ as the scattering within the $A$ side, minus what scatters out through $\partial$, and likewise for $B$.

It is instructive to examine the scaling of these separate processes. The “tunneling” process scales as the number of lattice sites along the boundary, i.e. it diverges as ${\cal O}(L)$. However, for each lattice site on the edge, there is a “tadpole” process that tunnels across at a site and then tunnels back at the same site, shown in (Figure 5 (c)). The weight of this also diverges (due to the logarthmic divergence of the Green’s function) as ${\cal O}(L)$ with equal but opposite coefficient to the tunneling process. In the thermodynamic limit, therefore, the divergences precisely cancel and $\widetilde{\gamma}$ is a constant. Visually, processes (a) and (c) in Figure 5 cancel each other out and leave only the (b) processes which “tunnel across and then tunnel back”. This is shown numerically for the rectangular geometry in Figure 5 (d). Thus we may conclude that the inter-sheet EPA, $\widetilde{\gamma}$, is generically non-zero and comparable in magnitude to $\widetilde{\alpha}$ and $\widetilde{\beta}$. By applying conformal mappings, this holds whenever $A$ and $B$ are simply connected and it is reasonable that the same conclusions extend to more general domains.

6 Continuum Redux

As with the ground state, the lattice and continuum results are closely analogous. The same structure of cancelling divergences also appears in the continuum and thus necessitates careful regularization of any integrals involved, which may be quite subtle. One can think of the lattice results as a guide for interpreting the continuum quantities.

A key tool to do this is Green’s second identity, which we recall here for convenience: suppose $f:U\subset\mathbb{R}^{2}\to\mathbb{R}$ and $G^{U}$ is the Green’s function on $U$ with Dirichlet boundary conditions (so $\Delta^{U}G^{U}=\text{Id}^{U}$), then for $w\in U$,

[TABLE]

where $\frac{\partial}{\partial n_{z}}=\hat{n}\cdot\nabla_{z}$ is the outwards-pointing normal derivative with respect to $z$.

For our case, $G^{U}(z,w)=0$ identically when $w\in\partial U$, so this reduces to

[TABLE]

where $f^{U}(w)=f(w)$ if $w\in U\setminus\partial U$ and [math] if $w\in\partial U$. The crucial point here is that if $f(w)$ does not vanish on the boundary, then the values of $f$ on the boundary are not preserved by $\text{Id}^{U}$, which only acts as the identity for functions which vanish on $\partial U$; the information on the boundary is lost.

Let us apply this to $\alpha$. Since $\Delta^{M}L_{\lambda}+\lambda L_{\lambda}=0$,

[TABLE]

Applying (59) for both $z$ and $w$ yields, after a calculation,

[TABLE]

Let $d_{B}$ be the same as $d_{A}$ with $A\leftrightarrow B$, so

[TABLE]

Then $\alpha+\beta+\gamma=\int_{A}d^{2}z\,L_{\lambda}^{2}(z)+\int_{B}d^{2}z\,L_{\lambda}^{2}(z)=1$ due to the fact that $L_{\lambda}$ is normalized.

The interpretation of $\int_{A}d^{2}z\,L_{\lambda}^{2}(z)$ and (and $A\leftrightarrow B$) is clear, but $d_{A}$ and $d_{B}$ are strange expressions: they involve the normal derivatives of $L_{\lambda}^{A}$ at the boundary of $A$ where $L_{\lambda}^{A}$ is discontinuous. The double normal derivative of the Green’s function on the boundary is subtle because the Green’s function vanishes at the boundary except when $z=w$ where it is singular. This hypersingular integral kernel arises in several other places, such as in the engineering literature of fracture analysis and always requires careful regularization [59].

From the lattice, we are motivated to interpret these quantities as

[TABLE]

Here $\hat{n}$ is the outwards pointing unit normal, $z,w\in\partial$, and $z-\varepsilon$ is a point $\varepsilon$ inside the boundary. So our regulated version is

[TABLE]

Both terms of this are, of course, superficially divergent quantities. But just as on the lattice, the divergences cancel to leave a finite quantity. Roughly this is because the normal derivative of the Green’s function restricted to the boundary is a delta function, so

[TABLE]

A more precise version of this is shown in C. Because of the minus sign, the $\mathcal{O}(\varepsilon^{-1})$ terms cancel

[TABLE]

This regulation scheme means $\alpha$, $\beta$, and $\gamma$ are well-defined and finite for any choice of geometry. In principle they can be calculated whenever the Green’s function is known. In practice, this requires either an explicit equation for the Green’s function on both $A$ and $B$. When the boundary is smooth, there is a general form for the $\mathcal{O}(1)$ in (66) in terms of Hadamard regularized integrals. Just as in Cardy and Peschel’s expression for the free energy [51], there appear to be special cases for corners, which generalize (43) in some sense.

As an example, consider the rectangle of side lengths $L_{x}$ and $L_{y}$ with an entanglement cut at $x=\ell_{x}$. The modes of the Laplacian are precisely Fourier modes in this geometry,

[TABLE]

with $\lambda=\left(\frac{\pi k_{x}}{L_{x}}\right)^{2}+\left(\frac{\pi k_{y}}{L_{y}}\right)^{2}$. In C, it is shown that

[TABLE]

where, for $L_{x}=L_{y}=L$,

[TABLE]

Here $\mathcal{S}(x)=\int_{0}^{x}\frac{\sin t}{t}\,dt$ is the sin integral. The approximations made, described in C, are accurate for large $k_{y}$. Figure 6 shows the comparison between numerical calculations of $\gamma$ on the lattice and (72), where good agreement can be seen for large $k_{y}$. One can also consider an artificial half-plane geometry with no corners, in which case the agreement is essentially perfect, shown in Figure 7. This is described in detail in C. This is strong evidence that EPAs are numerically equal in the continuum and on the lattice, and are hence universal.

7 Conclusions

In this paper we studied the entanglement of excited states in the quantum Lifshitz model. We found that there is an infinite class of excited states $\Psi_{\lambda}$ (reaching arbitrarily high energies) whose entanglement is completely specified by three geometrical, albeit state-dependent, “entanglement propagator amplitudes” $\alpha$, $\beta$, and $\gamma$ as defined in Equation (23a) – (23c). In terms of these, the second Rényi entropy is particularly simple:

[TABLE]

where $a=\alpha+\gamma/2$, $b=\beta+\gamma/2$ and $g=\gamma/2$. A combinatorial formula for arbitrary Rényi entropies is given in Eqution (24). This same procedure works for any state, but the combinatorics is generally complex. A particularly simple case are the states $\Psi_{K_{\lambda}}$ with $K$ quanta in mode $\lambda$, whose second Rényi entropy has the multinomial-esque expression

[TABLE]

The EPAs $\alpha$, $\beta$, and $\gamma$ have a natural interpretation as probabilities since they are all between zero and one and sum to unity. It is quite possible that these constants, which can be defined purely from Green’s functions of the Laplacian, may have already been studied in the mathematics literature under another name. The lattice results of Section 5 suggest that one can think of $\gamma$ as the total “flow” across the boundary and as $\alpha$ as the “energy” on the $A$ side, minus what flows out, and similarly for $\beta$. The continuum results of Section 6 show that these quantities are geometric but contain subtle cancelling divergences. Naively, since $\gamma$ can be defined as a quantity within a lattice site of the boundary, it appears to be a microscopic quantity. However, the comparison between lattice and continuum results gives strong evidence that $\gamma$ gives a leading order universal contribution to the entanglement.

As a by-product of working on the lattice, Section 4 reconstructed the ground state entanglement entropy on the lattice (41). Previous derivations of this fact [48, 49] made use of careful changes of variables in path integrals, whereas the lattice derivation relies renormalization group arguments. One can therefore consider this as independent confirmation of the non-topological part of the result (1).

In the special case where $\alpha=\beta$, we may use the replica trick [13] to determine the entanglement entropy. Starting from (27), and differentiating yields

[TABLE]

where $\Gamma$ is the incomplete gamma function. The infinite sum converges quickly, even when $\gamma=1$, and is $\mathcal{O}(1)$. With a sufficient application of combinatorics, it is undoubtedly possible to extend this to $\alpha\neq\beta$ and perhaps even to more general states. In all excited states with one quantum, the difference with the ground state entanglement is fully determined by the EPAs and is therefore $\mathcal{O}(1)$ and universal.

The notion of EPAs is entirely general and, in principle, works for any state whose wavefunctional is known. For instance, an eigenstate which populates $M$ modes can be written down in terms of $3\binom{M+1}{2}$ generalizations of $\alpha$, $\beta$ and $\gamma$ such as (25). Of course, more complicated wavefunctions have more complex combinatorics, making it harder to extract information. One possible direction would be to find a state with a particularly simple wavefunctional and examine the entanglement dynamics after a global or local quench. This was studied in the special case of a local quench by a vertex operator in [60]. The “propagating” nature of the EPAs is reminiscent of the light-cone picture of entanglement spreading, but since our results are at finite time, we see no connection between the two at present.

Another future direction is to examine entanglement at finite temperature. Since entanglement entropy is not a good measure of entanglement in mixed states, it is standard in field theory to instead calculate the logarithmic negativity which is an entanglement monotone and bounds several other common measures of entanglement [61]. It may be computed with a modified replica trick [62] which amounts to sightly modifying the combinatorics. At low temperature where only a few modes are populated, it may be possible to find an analytic continuation and determine the negativity. At high temperatures, one may be able to leverage the strong mathematical results on the asymptotic distribution of eigenvalues of the Laplacian into an asymptotic expression for the entanglement.

Acknowledgments.

We thank J. Cardy for helpful discussions and Adam Scherlis for suggesting tetrominos. We acknowledge support from the Department of Energy through the LDRD program of LBNL (RV), from NSF DMR-1507141 and a Simons Investigatorship (JEM), and center support from CaIQuE and the Moore Foundation’s EPiQS initiative.

Appendix A A Modified Wick’s Theorem for the Many-Sheeted Geometry

This Appendix will derive Equation (18), a version of Wick’s Theorem suitably modified for fields defined on the many-sheeted surface $\mathcal{S}_{n}$. The first uses a series of coordinate transformations introduced in [48] to decouple the sheets and the second evaluates the path integrals. The derivation has been written to be clear at the expense of concision.

Recall we are dealing with a compactified field $\varphi$, so $\varphi\equiv\varphi+2\pi R$ for some $R>0$. We can therefore rewrite the “surface expectation” (14) concretely as a constrained path integral

[TABLE]

where $\int_{\text{con}}$ means that the fields obey the constraint

[TABLE]

First, let us make the change of variables $\varphi_{a}=\varphi_{\alpha_{a}}+\widetilde{\varphi}_{a}$ where the first term $\varphi_{\alpha_{a}}$ is topological and encodes the zero modes and vortices, obeying $\Delta\varphi_{\alpha_{a}}=0$ with Dirichlet boundary conditions on $M$. The new field $\widetilde{\varphi}_{a}$ is non-compact. By applying Green’s first identity, the action decouples: $S[\varphi_{a}]=S[\varphi_{\alpha_{a}}]+S\left[\widetilde{\varphi}_{a}\right]$, and

[TABLE]

where the sum is over all possible zero modes and con’ is the constraint

[TABLE]

Adopting a matrix notation, one may express this as $M\bm{\widetilde{\varphi}}=\mathbf{c}_{w}(x)$ where one may choose (other choices are possible)

[TABLE]

where $w_{k}\in\mathbb{Z}$ and $\sum_{k}w_{k}=0$. Note that $\Delta_{a_{k}}\varphi(x_{k})=\Delta_{a_{k}}\left(\varphi_{\alpha_{a}}+\widetilde{\varphi}_{a}\right)=\Delta_{a_{k}}\tilde{\varphi}_{a}$, so a product of $\varphi$’s is equal to the same product of $\widetilde{\varphi}$’s, plus terms that are harmonic in at least one of the variables. These additional terms will not concern us because they can be removed with an integration by parts, and will be denoted Har. So after this first change of coordinates, the surface expectation is

[TABLE]

We now perform a second change of coordinates to decouple the sheets. Choose some $T\in O(n)$ such that $T_{1}^{a}=n^{-1/2},\forall a\in\mathbb{Z}_{n}$, and put $\bm{\psi}=T\widetilde{\varphi}$. There are many possible such matrices $T$. One explicit family [48] is

[TABLE]

Since the action is quadratic in $\widetilde{\varphi}$, $S[\widetilde{\varphi}_{a}]=S[\psi_{a}]$, and the constraint becomes

[TABLE]

This has the special property that $\psi_{1}$ is unconstrained, and the constraint on the other $(n-1)$ fields is an invertible matrix. In particular, they are constrained by

[TABLE]

Therefore

[TABLE]

The factor of $\sqrt{n}^{-L/a}$ comes from the Jacobian determinant and is non-universal.

To make it possible to actually evaluate the path integrals, we make one last coordinate transformation. For $a=2$ to $n$, let $C_{a}$ be the solution to the Dirichlet problem $C_{a}=c(x)$ on $\partial$, $C_{a}=0$ on $\partial M$, and $\Delta C_{a}=0$ on $M\setminus\partial$. This specifies $C_{a}(x)$ uniquely on $M$. Then define $\psi^{\prime}_{a}=\psi_{a}+C_{a}$ for $a>1$ and $\psi^{\prime}_{1}=\psi_{1}$. Again note that $\Delta\psi^{\prime}_{a}=\Delta\psi_{a}$, so the integrands of the path integrals are the same up to harmonic terms, and the actions again decouple, so

[TABLE]

where $\int_{M}$ means there is no constraint on $\partial$ and $\int_{D}$ means that the fields obey the Dirichlet boundary condition $\psi_{a}^{\prime}=0$ on $\partial$. At this point, the topological part, in parentheses, is completely separate from the path integrals. This is the end of the first part: the sheets have been decoupled and the path integrals may now be evaluated.

The path integrals now describe expectations of a free boson with two sets of boundary conditions: either nothing happens at $\partial$, or the fields are set to zero there. As in the above paragraph, these will be referred to as $M$ or $D$ respectively. Let us recall the standard (normalized) Wick’s theorems for these geometries

[TABLE]

where $Z_{D}=Z_{A}Z_{B}$, $G^{M}$ is the Green’s function for the Laplacian on $M$ with Dirichlet boundary conditions and if $G^{A}$ and $G^{B}$ are the corresponding Green’s functions for $A$ and $B$, then

[TABLE]

In the $D$ case, the fields are constrained to zero along $\partial$, so correlations may not be “transmitted” through $\partial$, but they can for the $M$ case. The sum over “pairs” may be written explicitly as the sum over the set of all possible partitions of $2k$ integers:

[TABLE]

This is necessary to keep track of the explicit indices below.

The first field $\psi^{\prime}_{1}$ is evaluated with the “$M$” Wick’s theorem, and the rest use the “$D$” version. Since the fields are now decoupled, the expectation of the product of $\psi^{\prime}_{a}$’s is the product of the expectations. Hence

[TABLE]

This serves as a definition for $S^{a}$. The overall prefactor is the ground state Rényi entropy

[TABLE]

where the first parenthetical is topological. All that remains is to rearrange terms to simplify the combinatorics:

[TABLE]

This implies Equation (18) for $N=n$.

Appendix B Arbitrary Rényi Entropies on the Lattice

In this Appendix we evaluate excited state entropies on the lattice. By definition, the $n$th Rényi entropy is $\widetilde{S}_{A}^{(n)}=\frac{1}{1-n}\ln\Tr\widetilde{\rho}_{A}^{n}$. We must compute

[TABLE]

The calculation proceeds in close analogy to the continuum case in Section 3. It is again sufficient to compute 2-point functions on the $n$-sheeted geometry, then apply Wick’s theorem and integrate against the eigenvectors of the Laplacian.

First, apply the the discrete Fourier Transform $\psi^{c}=V_{d}^{c}\varphi^{d}$ from (35), the $n$-sheets are decoupled as described in the previous section. The RG arguments from Section 4 may now be applied, whose result is that the fields $\psi_{2},\psi_{3},\dots,\psi_{n}$ have Dirichlet boundary conditions at $\partial$, so two-point functions on those sheets can only scatter from $A$ to $A$ or from $B$ to $B$, each according to the Dirichlet Green’s function for $A$ or $B$ respectively. Meanwhile, $\psi_{1}$, the center-of-mass sheet, has nothing to distinguish $\partial$, so scattering proceeds by the Green’s function on all of $M$. Therefore

[TABLE]

Here $\widetilde{G}^{M}$ is the Dirichlet Green’s function for the lattice Laplacian on $M$ and

[TABLE]

We can thus compute the two-point function back on the original sheets

[TABLE]

From the definition of $V$ and the fact it is unitary, it follows that

[TABLE]

Physically, the inter-sheet term term occurs because the “center of mass” Fourier mode connects the $A$ side to the $B$ side. If you start with a field on sheet $c$, then $1/\sqrt{n}$ is part of the center of mass sheet, which may scatter to the other side, and then $1/\sqrt{n}$ of that is on sheet $d$. Unlike in the ground state, we cannot cancel out the inter-sheet term with the normalization, because here it is additive instead of multiplicative.

Because this is a Gaussian theory (even on the $n$-sheeted surface), one may apply Wick’s theorem and find that $2N$-point functions are

[TABLE]

One can check that this is precisely the form of Equation (18) under the replacement $\widetilde{G}\leftrightarrow G$ (without the topological prefactor for the ground state entropy, since we have not compactified on the lattice). By the same reasoning as Section 3,

[TABLE]

where

[TABLE]

Appendix C Boundary Term in Explicit Geometry

C.1 Behavior of the Green’s function near the boundary

A key quantity to evaluate the probabilities $\alpha,\beta,\gamma$ is the Green’s function of the Laplacian with Dirichlet boundary conditions very close to the boundary (entanglement cut). Consider the Green’s function on the half-plane,

[TABLE]

Denote the Green’s function a distance $\varepsilon=x=x_{0}\ll 1$ away from the boundary by $G_{\epsilon}(u)$, with $u=\left|y-y_{0}\right|$ the coordinate along the boundary,

[TABLE]

It is easy to show that the normal derivative at the boundary is simply a delta function

[TABLE]

but we are actually interested in the subleading contribution to this delta function

[TABLE]

where the divergent contribution $-\frac{\delta(u)}{\epsilon}$ will cancel out another delta function contribution in (65). We would like to understand what $g(u)$ is (as a distribution). Since the logarithm is hard to deal with, it is convenient to consider the Taylor expansion of the derivative instead $\frac{1}{\epsilon}\frac{dG_{\epsilon}(u)}{d\epsilon}=-\frac{\delta(u)}{\epsilon}+2g(u)+{\mathcal{O}}(\epsilon)$, which differs from (104) by a factor of 2. Let us compute this quantity explicitly. It reads:

[TABLE]

that we conveniently rewrite as

[TABLE]

The first piece goes to zero as $\epsilon^{2}$ if $u\neq 0$, and its integral is $-1/\epsilon$, so it simply gives a delta function

[TABLE]

while we recognize the well-known expression of the Cauchy principal value in the second term

[TABLE]

Therefore, we find

[TABLE]

where ${\cal H}\frac{1}{u^{2}}=-\frac{d}{du}{\cal P}\frac{1}{u}$ is the derivative of the principal value, also known as the Hadamard regularization. In the following, we will denote the Hadamard regularization of an integral by a double-dashed integral

[TABLE]

with $x_{0}\in[a,b]$. We conclude that $g(u)=-{\cal H}\frac{1}{\pi u^{2}}$, so that

[TABLE]

This was derived explicitly for the half-plane, but conformal invariance implies this result should be true to leading order for simply connected geometries with smooth boundary.

C.2 Explicit calculation of the boundary term in the rectangle geometry

Let us now compute explicitly the value of $\gamma$ for a rectangle $L_{x}\times L_{y}$ with the entanglement cut in the $x$ direction between regions of size $\ell_{x}\times L_{y}$ and $(L_{x}-\ell_{x})\times L_{y}$. Our starting point will be the formula

[TABLE]

where the factor of $2$ in the second integral comes from the fact that the region $A$ and $B$ add up. The modes of the Laplacian read

[TABLE]

and $\lambda=(\frac{\pi k_{x}}{L_{x}})^{2}+(\frac{\pi k_{y}}{L_{y}})^{2}$. Both terms in $\gamma$ look badly divergent as $\epsilon\to 0$, but as explained in the main text, these divergent contributions cancel out. In principle, one should consider the full expression of the Green’s function on the rectangle, but for simplicity, we will use the simplified formula (102) which should hold far away from the boundaries of the rectangle ${M}$. This is equivalent to only considering the red poles in Figure 8. We do not expect this to change the answer significantly since the modes of the Laplacian all vanish near the corners for our geometry (the leading corrections due to these corners will be considered below). Let us expand in $\epsilon$. The first term of 112 is easily dealt with

[TABLE]

Note that there is no ${\cal O}(1)$ term in the integral. For the second term, we will use (111) (which once again should hold away from the corners where the entanglement cut $\partial$ meets the boundary of ${M}$) to find

[TABLE]

We see that in the limit $\epsilon\to 0$, the divergent contributions indeed cancel out, and we get

[TABLE]

Plugging in the expression for $L_{\lambda}$ and changing variables to $u=y-y_{0}$ and $v=y_{0}+y$, this yields

[TABLE]

The integrals can be computed explicitly using (110), and we find

[TABLE]

where we restricted ourselves to $L_{x}=L_{y}=L$ for simplicity.

This expression for the rectangle is approximate since we ignored the contribution of the corners to simplify the calculation. Let us now evaluate how the corners modify this result. We will focus on the leading correction coming from the corner at $x=0$, $y=0$ (the blue poles in Figure 8) so that the Green’s function near the vertical boundary ($x=x_{0}=\epsilon$) reads

[TABLE]

which vanishes for $y=0$ or $y_{0}=0$. The correction from the corner therefore amounts to changing $u$ to $v$ and comes with a minus sign:

[TABLE]

Evaluating the integrals for $L_{x}=L_{y}=L$, this yields

[TABLE]

The other corner at $x=0$, $y=L_{y}$ gives the same contribution, so that we obtain a modified expression taking into account the leading contributions of the corners as $\gamma_{\rm with\ corners}\approx\gamma+2\Delta\gamma_{\rm corner}$, so

[TABLE]

C.3 Numerical Verification

Equation (122) may be verified numerically by comparing it to results on the lattice. In particular, (122) was compared to (57) for a rectangle with height $L=9,18,\dots,273$ for $k_{x}=1$ and $k_{y}=1,2,3,4$. To extract the asymptotic behavior, the numerical results were fit to

[TABLE]

which is a good heuristic model for the results with $R^{2}>0.999$ in all cases. The coefficient $a_{1}$ should match the continuum result.

[TABLE]

It can be seen in Figure 6 that the agreement is excellent for large $k_{y}$.

The main source of the error comes from the approximations made to deal with the corners. To determine this, let us consider a geometry without corners. Consider an artificial half plane: use the half-plane Green’s functions but with a finite height $L$ for the boundary $\partial$. In the continuum, the answer is then given by (118). On the lattice, we need the half-plane lattice Green’s function, which is [63]

[TABLE]

where $L$ is a discrete version of the logarithm

[TABLE]

with $\gamma_{E}$ Euler’s constant, $\cos\lambda+\cosh\mu=2$ and $\mu/\lambda\to 1$ as $\lambda\to 0$. Using the same fitting procedure as before, the numerical results match the analytic ones closely, as shown in Figure 7 and the table below.

[TABLE]