Entanglement at a Two-Dimensional Quantum Critical Point: a T=0 Projector Quantum Monte Carlo Study

Although the leading-order scaling of entanglement entropy is non-universal at a quantum critical point (QCP), sub-leading scaling can contain universal behaviour. Such universal quantities are commonly studied in non-interacting field theories, however it typically requires numerical calculation to access them in interacting theories. In this paper, we use large-scale T=0 quantum Monte Carlo simulations to examine in detail the second R\'enyi entropy of entangled regions at the QCP in the transverse-field Ising model in 2+1 space-time dimensions -- a fixed point for which there is no exact result for the scaling of entanglement entropy. We calculate a universal coefficient of a vertex-induced logarithmic scaling for a polygonal entangled subregion, and compare the result to interacting and non-interacting theories. We also examine the shape-dependence of the R\'enyi entropy for finite-size toroidal lattices divided into two entangled cylinders by smooth boundaries. Remarkably, we find that the dependence on cylinder length follows a shape-dependent function calculated previously by Stephan {\it et al.} [New J. Phys., 15, 015004, (2013)] at the QCP corresponding to the 2+1 dimensional quantum Lifshitz free scalar field theory. The quality of the fit of our data to this scaling function, as well as the apparent cutoff-independent coefficient that results, presents tantalizing evidence that this function may reflect universal behaviour across these and other very disparate QCPs in 2+1 dimensional systems.


Introduction
At a quantum critical point (QCP), the Rényi [1] entanglement entropies contain the muchcelebrated ability to access the central charge of the associated conformal field theory (CFT) in 1 + 1 space-time dimensions [2,3]. This has provided the mainstay in a healthy dialog between field theory and numerical lattice simulations, where unbiased methods such as density matrix renormalization group (DMRG) [4] are able to calculate the Rényi entropies of order α precisely in quantum models in one spatial dimension (1D). The form for systems with open (η = 2) or periodic (η = 1) boundaries in 1D [5,6], can be straightforwardly compared to lattice numerics [7], where measuring both L and x in terms of the lattice spacing a gives this term access to the universal central charge, c.
The simplicity of the entangled boundary in D = 1 means that this and other geometrical factors, such as occur when Rényi indices α > 1 [8], can be compared between field theory and numerics to a high degree of accuracy in 1 + 1 space-time dimensions. A similar success is only beginning to be enjoyed in D > 1. Scalable finite-size lattice simulation methods which can address behaviour at critical points (where the length scale diverges), such as quantum Monte Carlo (QMC), predominantly measure the strongest signal from a looming non-universal leading-order contribution proportional to the boundary length [9,10], L D−1 /a D−1 = (except in the presence of a Fermi surface, where even higher entanglement is expected [11,12]). The subleading terms that exist must be obtained using subtraction of very large statistically-fluctuating contributions from this "area law" (unless they can be calculated separately, as in infinite-lattice linked-cluster expansions [13]). What is left is a universal piece, which may depend on the shape or topological characteristics of the boundary, but is believed to not depend on the entangled volume. This geometry-dependence can be exploited to access different universal numbers in strongly-interacting models, but the types of geometries amenable to comparison between lattice-model simulations and continuum theories has so far been limited.
One promising function to access universal quantities in 2+1 is the logarithmic contribution that comes from vertices in a polygonal-shaped entangled region [14,15,16,17,18]. The universal coefficients, which depend on the vertex angle θ, have been compared in the past to calculations on finite-size lattices where the entangled region is a square, with four independently-contributing corners (θ = π/2). Particularly relevant are recent numerical results on interacting models [13,19,20], however to date, the only field theory calculations that have been performed are on non-interacting theories [15], inhibiting a quantitative check of universality.
However, other shape-dependent contributions, that occur with smooth boundaries but otherwise impart important geometric characteristics of a 2+1 dimensional space-time geometry, are also accessible by a simulation cell cut into two subregions (e.g. a torus bifurcated into two cylinders), and reveal sub-leading contributions to the Rényi entropy. Speculation exists about whether this shape dependence reveals an underlying universality [21,22]. Quantum Monte Carlo simulations on interacting QCPs in 2 + 1 dimensions are uniquely poised to answer this question.
In this paper we study the subleading shape dependence of the Rényi entanglement entropy on L × L toroidal lattices with spatial dimension D = 2, using the critical point of the transverse field Ising model (TFIM), where − → σ i is a Pauli spin operator, accessed by a novel projector QMC algorithm that works strictly at T = 0. Simulating a variety of lattice sizes up to 40 × 40 reveals universal vertex contributions that are close to a non-interacting field theory in 2 + 1. For smooth boundaries bifurcating the torus into two cylinders, subleading contributions show close functional form to Eq. (1) when the two cylinder lengths are associated with x and L − x. However, a cutoff-(or L/a-) independent coefficient is not obtained on finite lattices. Instead, if we look at the well-studied 2+1 dimensional quantum Lifshitz model [14,23,24,25] (the field-theory associated with the Rokhsar-Kivelson (RK) Hamiltonian [26]) there is a proposed functional form [22] for the universal subleading term to the area law that is in good agreement with our data with a size-independent coefficient. This allows us to speculate on the universality of the scaling function, and the interpretation of its coefficient as a measurable witness to the universality class. Such potential demonstrates the importance of obtaining efficient simulation methods for calculation of Rényi entropies in strongly-interacting lattice models, and the continuing need for field theory calculations on interacting fixed points in 2+1 dimensions.

Entanglement at strongly-interacting critical points in 2+1 D
Extensively studied in one spatial dimension, a multi-disciplinary community is beginning to examine the scaling of entanglement entropy in two (spatial) dimensions, thanks to the rapid development of both theory, and numerical methods for calculating entanglement-related quantities in scalable simulations. With the important exception of many-body systems housing a Fermi surface [11,12], the prevailing paradigm for entanglement in ground state wavefunctions is the area (or boundary) law [9,10], where A is a non-universal constant, = L D−1 /a D−1 is cutoff-dependent, and the ellipses indicate a combination of different subleading corrections. Heuristically, the existence of an area law can be related to the finite extent of correlations across the entangling boundary [27,28]. In gapped systems, the subleading term may have contributions from several constants, some universal, some not. The most important universal subleading correction gives a contribution dependent on the topology of the entangled surface in a fractional topological phases, and is called the topological entanglement entropy: where e.g. γ = log(2) for a simple Z 2 spin liquid [29,30]. In 2D critical systems, the behaviour of these subleading corrections to the area law become much more rich. One may naively suspect that, due to the diverging correlation length, a violation of the area law might be possible. However, it can be understood through a course-graining picture that the area law is still obeyed. First, assume that due to scaleinvariance each length scale (in a renormalization group (RG) sense) contributes order O(1) bit of entanglement entropy across the boundary. One may take scale-invariance to mean that, when rescaling the system by some factor b, the number of modes at this new length scale is proportional to the new boundary length /b. Then, using this assumption and summing over all length scales, an entropy proportional to the boundary length is obtained.
‡ At criticality, additional universal subleading terms to this area law are also possible, however they may have a complicated dependence on the geometry of the bipartition. Although typically believed, it is not generally known if particular geometric features, for example the number of vertices or the Euler characteristic [14,31,17,18], give rise to certain universal numbers that can be compared reliably between field theories and quantum lattice models. This would be important in making progress towards developing an analog of c-theorems [32] in 2+1 space-time dimensions, which aim to identify a universal function with monotonic behaviour under RG flows [33]. The fact that most field-theoretic calculations are limited to non-interacting systems hampers progress in this regard. In order to study interacting systems, one must turn to numerical techniques on finite-size lattice. Understandably, the geometries amenable to study on finite-sizes lattice are sometimes different than those that can be studied with continuum field theories, as we now discuss.

Bipartitions with smooth boundaries
In the continuum thermodynamic limit, dividing a systems into two partitions A and B is easily done with a smooth curved boundary. This geometry has become particularly important in entanglement monotonicity studies, where for example, in Lorentz-invariant theories, Casini and Huerta have shown that the entanglement of a smooth circle of ‡ We thank M. Hastings for pointing out this argument, which might possibly be traced back to J. Preskill. circumference scales as S 1 ∼ A − γ , where the universal constant decreases along an RG flow [34]. It is thus a compelling candidate for the monotonic c-function. This may also be related to the entanglement of a three-sphere in odd space-time dimensions, which contains a subleading constant term which changes monotonically along RG flows in holographic theories [35,36]. Such developments may provide a route to a (2+1)-dimensional analog to the c-theorem, especially if the fixed-point value of the monotonic quantity could be determined.
In essentially any interacting theory (and some non-interacting theories) in 2+1, this task would necessarily fall to numerical simulations. Unfortunately, such curved geometries are inaccessible on lattices, where obtaining boundaries with sharp vertices is unavoidable when attempting to draw smooth curves. In the case where the vertices disappear in the continuum, it is unknown how such "pixelization" might affect the approach to the thermodynamic limit. In the case where vertices or corners remain in the geometry in the continuum, they will contribute an additional universal factor, as we discuss in Section 2.3.
Additionally, one may also examine the entanglement entropy across a smooth boundary perturbed away from the critical point, where the correlation length becomes finite. In this case, even if the boundary has no curvature, one expects [16], for spatial dimensions D. Since the definition of flat boundaries is possible in lattice systems, recent numerical calculations on interacting lattice models have been able to make important comparisons of r α calculated perturbatively in interacting field theories. In these numerical works, series expansion techniques are able to capture entanglement contributions across flat boundaries on infinite lattices by systematically including larger cluster sizes [19,13]. An important distinction between these numerical series expansion techniques and quantum Monte Carlo (QMC) is the fact that the latter is typically restricted to periodic finite-size systems (i.e. toroidal lattices of size L × L), meaning that these methods approach the thermodynamic limit in a different way. A smooth spatial boundary in a two-dimensional toroidal lattice is possible only if one bifurcates the torus into two separate cylinders. This geometry can also be studied in certain field theories, which is important for addressing the potential universality of cutoff-independent subleading scaling terms in the Rényi entropies.
In this paper will we study this geometry in the context of the TFIM 2+1 QCP in great detail.

Bifurcated torus: two-cylinder entropies
As discussed above, in two spatial dimensions QMC simulations typically take place on toroidal lattices of size L × L. For the measurement of entanglement between two subregions A and B, a simple geometry is illustrated in Fig 1, where upon bifurcation two cylinders of "length" x and L − x are produced. Past numerical studies of strongly-interacting gapless systems have demonstrated that this geometry is a sensitive probe of the entanglement structure of the wavefunction [37,38,22]. Two possible features are particularly prominent as one varies the length x: a striking even-odd effect which arises due to dimer-like physics in the wavefunction; and, a smooth x-dependent curvature that may contain universal scaling behaviour. This smooth x-dependent curvature is not present in gapped states, making it a sensitive indicator for gapless behaviour in general wavefunctions [38] (see Section 2.2.2).
2.2.1. Even-odd effect: As first observed in Ref. [38], a striking "even-odd" branching effect is observed in S n (x/L) in certain RVB wavefunctions. In Ref. [22], this effect was understood in a free scalar field theory in the 2+1 dimensional quantum Lifshitz model, which is the fieldtheory of a Rohksar-Kivelson (RK) Hamiltonian. It arises due to the underlying dimerization of the wavefunction -and not from any underlying symmetry breaking, as one might naively expect for example in a valence-bond solid phase.
Away from this exactly-soluble free field theory, the even-odd effect serves as a sensitive probe of the degree of dimerization in the low-energy effective description of the wavefunction. For example, the Néel ground state of the 2D spin-1/2 Heisenberg model can be described in an RVB singlet-basis where long singlets, decaying as 1/r 3 , occur [39]. Correspondingly, no even-odd branching effect is observed in the Rényi entropy S 2 [37]. Similarly, for QCPs which are described by theories "sufficiently" far from RK-like Hamiltonians, it's reasonable to expect that no even-odd branching occurs. As we will see in this paper, in the 2+1 dimensional QCP in the transverse-field Ising model (TFIM), no such even-odd branching occurs.

2.2.2.
Length dependence from conformal field theories: Previous numerical studies of interacting models on two-cylinder entropies show a clear geometry-dependent function, that occurs even in the presence of even-odd branching (described above), and which depends on the cylinder length x. Previous results on the Néel ground state of the Heisenberg model, and the square-lattice RVB wavefunction showed reasonably good numerical fits [38] to the shape-dependence motivated from 1+1 CFT in Eq. (1), specifically, where however b = c in general. The term proportional to log(sin(πy)) is heuristically included in order to account for a strong shape-dependence observed in these gapless wavefunctions in 2+1, in analogy to the 1+1 exact result. Importantly, a non-zero logarithmic term b of order unity was first observed in QMC data on the Heisenberg model in Ref. [37]. This phenomena was subsequently explain by Metlitski and Grover [40] as arising from the presence of Goldstone modes and the restoration of symmetry in a finite-volume simulation cell. The coefficient of this additive logarithm should be universal, b = N G (D − 1)/2, where N G is the number of Goldstone modes. Thus, it is only expected to be present in the case of spontaneously broken continuous symmetry. It has previously been demonstrated not to exist in an exactly-solvable finite-size model of free spinless fermions on a square lattice, with π flux through each plaquette [38]. As we will see in Section 4.1, QMC data for the TFIM quantum critical point is also consistent with b = 0.
Recently, motivated by the study of dimer RVB wavefunctions that in the continuum limit have conformal invariance in 2D space, CFT techniques have been used to derive an alternate scaling form of the shape-dependent piece for Fig. 1 [22,41]: where y = x/L is the fractional width of the strip (referred to as y /L + y in some of the previous literature), η is the Dedekind eta function and θ ν is the Jacobi-Theta function. The above form for J α (y) applies for the Rényi entropies with α ≥ 2, but the definition does not extend to the von Neumann entropy. Through the rest of the paper we use J(y) = J 2 (y) for simplicity. Due to our geometry (always using L×L systems), τ = iL x /L y = i never changes in the above equation. Although this universal function was derived for the quantum Lifshitz field theory, it is interesting to test its universality away from the critical points describing RK and RVB-like wavefunctions. In the present paper, we explore its potential for universality by comparing the functional dependence of the Rényi entanglement at the TFIM QCP to fits of Eq. (6) and Eq. (8).

Polygons on a torus: entanglement due to vertices
It is well-known that, in critical systems in 2+1 dimensions, another cut-off dependent contribution to the entanglement entropy distinct from the boundary ("area"-law) is induced by the presence of vertices or corners. This term can be seen to have a logarithmic dependence on through heuristic mode-counting arguments in a renormalization group framework.
To begin, one assumes that each length scale (in an RG sense) contributes O(1) bit of entanglement entropy for each vertex or corner. Then, summing the contributions from each length scale, one discovers that the corners contribute a constant times the total number of length scales, log( ). § Refinements and generalization of this argument to other geometrical boundaries exist in the literature. Most promising for comparison between analytical field theory and lattice numerics is the simple 90-degree vertex. In Ref. [15], the full angle-dependence of the logarithmic vertex contribution is explored, and found to obey, where n c is the number of vertices, and a α (θ) is a universal quantity independent of lattice cutoff. For the square lattices studied in this paper, it is natural to use a vertex of θ = π/2. The quantity a α (π/2) will be universal between the continuum field theory and the lattice model, provided that the angle in the lattice entangled region approaches a 90 • corner as one increases the number of sites to infinity. Past numerical studies have calculated a α (θ) for a variety of Rényi indices α at the QCP of the TFIM. Of relevance to the present study, the value of a 2 (π/2) has been calculated several times in the past literature. A value for the non-interacting fixed point in a scalar field theory was provided by Casini and Huerta, a 2 (π/2) = −0.0064 [42]. Series expansion studies of the interacting 2+1 dimensional QCP of the TFIM give −0.0055(5) [19] while Numerical Linked-Cluster Expansion (NLCE) gives −0.0053; both results are consistent with each other, and lower than that result of the free field theory. Finally, previous finitetemperature QMC calculations on a single 36 × 36-size lattice report −0.0075 (25) [20]. In the current paper, we aim to improve on these QMC results by employing a new projector method that converges the Rényi entropy S 2 directly at T = 0, as now described. § See footnote 1.

Projector Quantum Monte Carlo
Due to a diverging correlation length, the numerical study of quantum phase transitions in strongly-interacting models requires extraordinary care [43]. Every effort must be taken to converge data on lattices of as large a size as possible, such that reliable finite-size extrapolations may take place. Quantum Monte Carlo (QMC) on sign-problem free models provide an unbiased numerical procedure to systematically approach the thermodynamic limit. For lattice spin models, such as the transverse-field Ising model (TFIM) studied here, the standard procedure used to access quantum critical behaviour involves using highlyefficient finite-temperature algorithms, such as continuous world-line [44] or Stochastic Series Expansion (SSE) [45,46], operating at sufficiently low temperatures. Then, the temperature T is either decreased until it is converged to its T → 0 behaviour for each lattice size (and parameter set) of interest; or, the inverse temperature β = 1/T is set proportional to the linear lattice size, β ∝ zL (where z is the dynamical scaling exponent), for each lattice size studied.
Recently, a new flavour of QMC algorithm has emerged that combines the simplicity and efficiency of SSE QMC, with the ability to study ground state properties of model Hamiltonian directly at T = 0 [39]. Such "projector" methods share many of the features of their T > 0 counterparts, such as: the direct numerical coding of a D-dimensional quantum system to a D + 1 dimensional classical configuration, which is represented and sampled on a computer; and the existence of the prohibitive "sign-problem" for fermonic and frustrated systems. Unlike T > 0 simulations, these methods do not operate in a D + 1 simulation cell that is periodic in the temporal direction; rather, they operate with a Hamiltonian on a trial state, repeatedly, projecting out the ground state, as explained below. These projector methods have been widely adopted for the study of T = 0 properties of SU(2) (and SU(N ) [47]) invariant Hamiltonians with singlet ground states, providing a simple and efficient method to converge ground state properties [39,48,49]. In the next section, we describe a procedure, first developed by Sandvik [50,51] (and recently generalized [52]), for the efficient simulation of the U(1) symmetric TFIM Hamiltonian, using an adapted projector QMC method which employs non-local cluster updates. In Section 3.2, we describe how this algorithm may be adapted to calculate the Rényi entanglement entropies using a straight-forward adaptation of the "replica" trick [53]. It is this method which allows us to accurately study the finite-size scaling of S 2 at the QCP of the TFIM, and to access the universal subleading quantities of interest.

Algorithm for the transverse-field Ising model
In 2005, Sandvik introduced a ground state projector QMC method for SU(2) quantum spins, using the so-called valence-bond basis [39]. At T = 0, QMC methods are tasked with calculating the operator expectation value, Here, one aims to use some procedure to sample ψ 0 , the ground state wavefunction of a Hamiltonian, where the normalization is Z = ψ 0 |ψ 0 . The transition probabilities of the Metropolis algorithm are based on this overlap; the non-orthogonality of the valence-bond basis ensures that this is trivially non-zero. However, as we will see below, this trial state is not strictly required to be an overcomplete basis. In the case of the TFIM, an orthogonal σ z basis can be used, provided that wavefunctions are sampled such that this overlap is non-zero, using a non-local cluster algorithm. In a projector QMC representation, the ground state wavefunction is estimated by a procedure where a large power of the Hamiltonian is applied to a trial state, call it |α . To see the projection of the ground state wavefunction, one can write the trial state in terms of energy eigenstates |ψ n , n = 0, 1, 2 . . ., |α = n c n |ψ n , so that a large power of the Hamiltonian will project out the ground state, Here, we have assumed that the magnitude of the lowest eigenvalue |E 0 | is largest of all the eigenvalues. To achieve this, one may be forced to add a sufficiently large negative constant to the overall Hamiltonian (that we have not explicitly included). Then, from this expression, one can write the normalization of the ground state wavefunction, Z = ψ 0 |ψ 0 with two projected states (bra and ket) as, for large m. In a procedure that will be familiar to any SSE aficionado [54], the Hamiltonian is written as a (negative) sum of elementary lattice interactions the indices t and a referring to the operator "types" and lattice "units" over which the terms will be sampled. In order to represent the normalization as a sum of positive-definite weights we can insert a complete resolution of the identity between each H t i ,a i , Where this equation has been cast in a form similar to that for finite-T SSE. Note that, the sum over the set {α} and the operator list S m = [t 1 , a 1 ], [t 2 , a 2 ], · · · [t 2m , a 2m ] must be done with importance sampling. As we will see below, an update procedure can be constructed that efficiently samples both the list of operators H t,a , and (separately) the left and right basis states. Thus, for sufficiently large m, any trial state |α can be chosen, which is equivalent to using the equal superposition of all spin states σ z i [51]. Other choices, such as a variationally optimized state, may also be used [51].
Turning to the Hamiltonian Eq. (2), a convenient definition of operator types is, Note that the index a has two different meanings: a site or a bond, depending on the operator type. It is evident that some simple constants have been added to the Hamiltonian (Eq. (2)): the diagonal operator H 0,a , and also the +1 in Eq. (17). The first results in matrix elements with equal weight for both one-site operators. Denoting each matrix element in the standard basis: the σ z i = +1 eigenstate is | • i and σ z j = −1 is | • j , the non-zero matrix elements are, The D + 1 dimensional projected simulation cell is built such that 2m operators of the type 15 to 17 are sampled between the "end points" (i.e. the trial states). Then, sampling occurs via two separate procedures, as follows. First, the diagonal update where one traverses the list of all 2m operators in sequence, e.g. from |α l to |α r . If an off-diagonal operator H −1,a is encountered, the σ z spin associated with that site is flipped but no operator change is made. If a diagonal operator is encountered, the Metropolis procedure is: (i) The present diagonal operator, H 0,a or H 1,a , is removed.
(ii) A new operator type is chosen, t = 0 or t = 1, corresponding to the insertion of either a diagonal h or a diagonal J operator, i.e. Eq. (19) or (20). The transition probability to add H 0,a is, where N and N b are the number of sites and bonds in the lattice, respectively. Note, P (H 1,a ) = 1 − P (H 0,a ).
(iii) If H 0,a is chosen, a site a is chosen at random, and the operator is placed there.
(iv) If H 1,a is chosen, a random bond a is chosen. The configurations of the two spins on this bond must be parallel for the matrix element to be nonzero. If they are not, then the insertion is rejected. Steps (ii) to (iv) are repeated until a successful insertion is made. One can see that this diagonal update is necessary in order to change the topology of the operator sequence in the simulation cell. However, in order to get fully ergodic sampling of the TFIM Hamiltonian operators, one must employ non-local updates in addition to these simple diagonal updates. For T = 0 projector methods, non-local updates have been discussed previously in the literature in the context of the valence-bond basis [49]. In the present case, the TFIM Hamiltonian does not conserve σ z , and a different type of branching non-local update, called a "cluster" update, must be used. We use cluster updates adapted from the finite-T SSE procedure described in Ref. [50]. A crucial observation is the judicious choice of H −1,a and H 0,a to both have the weight h, which allows for unrestricted sampling between the two operator types. One can easily see that a functional definition of non-local clusters will be an Ising spin forming a space-time group, bounded by either single-site operators, or by spin states of the end point trial states |α l and |α r (see Fig. 3). If all spins within a cluster are flipped, the total weight of the configuration remains unchanged. One is then free to build all clusters deterministically, and flip each with a Swendsen-Wang algorithm, i.e. a probability of 1/2. In this way, one sees how a fully ergodic sampling of both the operator types and basis state σ z is sampled in the projector QMC.
In Fig. 3 special care has been taken to note the mid-point of the simulation projection, since operator expectation values are measured there: It is particularly important to reiterate, especially in anticipation of the next section, that the wavefunction overlap ψ 0 l |ψ 0 r must be non-zero (indeed, the Metropolis sampling is set up to force this). This is ensured by requiring that spin states are the same on the left and  Figure 4. An illustration of the "SWAP A " operator acting on a replicated simulation cell, for the calculation of S 2 . Basis states and TFIM operators follow the same legend as Fig. 3. Here, each simulation has three physical spins, and m = 3 operators in each non-interacting replica, represented one on top of the other. At left, clusters within each replica that cross the mid-point of the 1+1 simulation cell are highlighted in different colours. At right, a SWAP A operator has been applied to permute the basis states corresponding to region A, between each replica. The reconfiguration of the clusters crossing the mid-point can be seen from their colours. the right of the "mid-point" of the D + 1 dimensional cell. Since spins are only constrained to be the same within a given cluster, a proper normalization for expectation values could be defined as, where N 0 is the number of independent clusters that cross the boundary. This follows from the fact that each connected cluster in Fig. 3 has two possible spin orientations, σ z = 1 and σ z = −1, independently. Various relevant measurements, such as the ground state energy, can be constructed for this simulation [51]. However, for the purposes of this paper, we are interested specifically in the Rényi entanglement entropies, which only require knowledge of the cluster structure of the simulation at the middle of the propagated simulation cell; albeit in a non-trivial replicated (or multi-sheet) geometry. Thus, in the next section, we describe the specific implementation of the Rényi entropy estimator for the projector QMC for the TFIM.

Measuring Rényi entropies through the SWAP A operator
In a development crucial to the study of entanglement at interacting quantum critical points, QMC methods have recently been introduced that are capable of measuring the degree of entanglement in a ground state wavefunction, specifically using the Rényi entropies [1], for integer α ≥ 2. Here, ρ A is the reduced density matrix, ρ A = Tr B {ρ}, and A is a subregion of a lattice system (with B being its complement) -see Figs. 1 and 2. The von Neumann entanglement entropy corresponds to the limit α → 1. Unlike a linked-cluster expansion or other methods based on Lanczos diagonalization [13], the direct measure of Rényi entropies can not be done in QMC using conventional estimators. However, recent work has demonstrated that it is possible to measure S α for integer α ≥ 2 using a replica trick [6,14,55,56,57]. For T = 0, the calculation of S α is done via the following procedure: (i) The system is copied into α independent "replicas".
(ii) Each replica is independently projected to sample the ground state. The total wavefunction for the system of all replicas is denoted by |ψ 0 .
(iii) Following Eq. (22), the operator O is replaced by the "SWAP A " operator for α = 2 [53] or permutation operator for α ≥ 3 [37]. This operator literally swaps (or permutes) basis states in the region A between the α copies. See Fig. 4 (iv) The expectation value where N A is the number of independent clusters crossing the middle of the re-connected ("swapped") partition function, and N 0 is the number of clusters before the swap took place. Note: here N 0 is the number of clusters in all α replicas or copies.
(v) Finally, Eq. (24) is applied: e.g. S 2 = − log SWAP A Note that, Step (iv) comes from the fact that Eq. (22) is used, with O = SWAP A as the operator being measured. This results in the simple procedure of measuring the overlap in the numerator and denominator -both of which are simply the number of spin states per cluster (2) raised to the power of the number of independent clusters crossing the middle of the D + 1 dimensional simulation cell.
As first pointed out in Ref. [53], as the lattice size (and particularly the size of region A) grows, sampling statistics become exponentially poor when using the naive SWAP A operator as described above. Hence, a slight adaptation called the ratio trick must be used in order to improve statistics. The ratio trick involves calculating the Rényi entropy in several steps, each step being a separate simulation that involves sampling a ratio: where X denotes a subregion that is spatially smaller than A. In other words, the weight of a simulation becomes (un-physically) related to a partially-swapped simulation cell. In a practical simulation, many spatial sub-regions X i , each employed as the weight of a separate simulation, are used to build up towards the physical bi-partition A of interest.
The Rényi entropy is then built up by multiplying the contribution of Eq. (26) from each spatial subregion. For the TFIM simulation discussed in this paper, the procedure for efficiently calculating the Rényi entropy closely follows that used in another T = 0 projector QMC -the valencebond basis QMC for the spin-1/2 Heisenberg model [39,53], which has a detailed description (including the ratio trick) in Ref. [37]. Remarkably, the only algorithmic difference between the two QMC algorithms is the structure of the space-time clusters employed in the projector method. In Eqs. 25 and 26, the N -numbers (N A , N 0 , N X ) count the number of branching clusters, spanning both the spatial and propagation directions, that cross the middle of the simulation cell. In the spin-1/2 Heisenberg model [37], this number counts the non-branching loop structures that cross this middle point, due to the different nature of Hamiltonian operators in that model. As in the case of the TFIM clusters, each Heisenberg loop in the valence-bond representation has two spin states associated with it (for SU (2); this is modified to N for SU(N )). It is remarkable that Hamiltonians with different symmetries, and completely different basis-state representations in the projector QMC, end up with equivalent measurement procedures for the Rényi entanglement entropies.

Convergence of S 2 at a quantum critical point
In this paper, we examine the second Rényi entropy, setting α = 2 in Eq. (24). Like any other observable measured with the T = 0 projector method, the value of S 2 will have its own unique convergence properties, as a function of the projector length m, for each value of h/J studied. When focussing on the critical point, h/J = 3.044, where the correlation length diverges, one must be particularly careful to ensure that the simulation is converged in m for each lattice of size N = L × L.
In Fig 5, we examine the value of S 2 , using the two-cylinder geometry of Fig. 1, at the point when the two entangled cylinders are of equal size, x = L/2. Of the geometries studied in this paper, this point is expected to be the most difficult to get good statistics on, and we use the ratio trick, Eq. (26) to converge it, building up each subregion X i with at most 1 × L sites (less for larger sizes). See also the discussion in the next paragraph for a comparison to data taken without the ratio trick. For each of the four system sizes that we have studied in detail, we see that the value of S 2 at the centre-point x = L/2 converges for sufficiently large m -requiring a slightly larger value of m/N as N increases. For the largest system size that we collect detailed convergence data on, L = 22, the value of S 2 saturates between 50000 < m < 60000 operators, i.e. m/N slightly larger than 100. We continue to collect data (see Fig 5) for larger m, but find that the value of S 2 is essentially converged within error bars. For the data collection in the rest of the paper, we settled on a fixed m/N = 160, which is more than sufficient to converge L = 22 and larger. A comprehensive study of the m-dependence of the Rényi entropy for different values of α and N would be an interesting topic of future study.
In the inset of Fig. 5, we see raw data for the x-dependence of S 2 using the geometry in Fig. 1, on an N = 24 × 24 lattice with m/N = 160. There, we have compared data obtained from a single simulation using a bare measure of the SWAP A operator, Eq. (25), with that obtained from a procedure where the ratio trick, Eq. (26), is used to build up each entangled region A. Clearly, naive measurement of the bare SWAP A operator is insufficient to get controlled statistical sampling for large x values. We find that, in the case of the TFIM at h/J = 3.044, the ratio trick is absolutely necessary for simulation of size L ≥ 16.

Simulation results on finite-size lattices
In this section, we report results for the second Rényi entropy, S 2 , for the Hamiltonian Eq. (2) precisely at the QCP (h/J = 3.044) using the QMC simulation method discussed in the last section. We discuss two entangling geometries: a bipartition between A and B that smoothly cuts each torus into two cylinders (Fig. 1), and a square bipartition with four 90-degree corners (Fig. 2).

Bifurcated torus: two-cylinder entropies
The first geometry we consider is that of an L × L torus, where the two entangled regions result from smoothly cutting the torus into two cylinders, as shown in Fig. 1, where the length of each cylinder is x × L and (L − x) × L.
Before examining the full x/L dependence of the two-cylinder geometry, we discuss the possibility of a non-zero b in Eq. (6) by examining regions A of a fixed x/L embedded in different finite-size lattices L. In Eq. (6), one may eliminate the x-dependence of the term proportional to c by fixing x/L to be a constant; e.g. if x is L/2 or L/4, this contribution will be absorbed into the additive constant d. Fig. 6 shows the residuals, that is S 2 (L) − aL − d, for two different choices of x/L as a function of system size for the half and quarter cylinder partitions, with b explicitly set to zero. From this, we conclude that the fit is acceptable within statistical errors. If, instead, we allowed a b = 0 as a fit parameter, a very small negative coefficient is found, b ≈ −0.01, two orders of magnitude smaller than that found in systems with continuous symmetry breaking [37]. As such, we argue that this value is inconsistent with both physical expectations (see Sec. 2.2) and the results below (see Fig. 8). We conclude that our data is consistent with the absence of a subleading log( ) dependence in the case of smooth boundaries.
In the following, we assume that no additive logarithms are present in the entanglement scaling at the TFIM quantum critical point; b = 0 in Eq. (6). Fig. 7 then shows the entanglement entropy as a function of the cylinder size, x, for our three largest system sizes, L = 28, 32 and 36. Recall that in some other gapless wavefunctions, e.g. the RVB  (Fig. 1), along with fits to (orange) Eq. (27) and (teal) Eq. (28). Notice the lack of any even-odd effect in the entanglement as a function of cut length.
wavefunction studied in Refs. [38] and [22], a very prominent branching effect is apparent which produces separate entanglement curves for even-x and odd-x (see Sec. 2.2.1). As discussed by Stephan et al. [22], this even-odd effect is a measure of how RVB-like a wavefunction is. In the present case, it is clear that if any even-odd effect occurs, it is beyond our ability to detect in these QMC simulations.
Next, we are interested in examining the functional dependence of the shape of the curves in Fig. 7. To do this, we look at a variety of system sizes, L, and a variety of entangled-region lengths, x, and examine the fit of the entropy to two equations. The first is the form log(sin(πx/L)), Eq. (1), which is relevant for 1+1 dimensional systems, and which has been used in an ad hoc way in the past analyze some 2+1 dimensional entanglement entropy data [21] including systems with a 2+1 dimensional fermi surface [38] where Eq. (1) gives a reasonable (but not perfect) approximation to the fit. The second is the RVB shape function derived by Stephan et al. [22], Eq. (8). For the purpose of a numerical comparison, we fit to functions of the form, S log = a + c L log(sin(πy)) + d, with y = x/L and where the constants a, c L , and d may be different for the above equations. Note that we allow the freedom that c L can vary with system size. However, a and d are fit by considering all systems sizes together. In this way, we can use the variation of c L with L to examine the quality of fit for each of the two functions. Fig. 8 shows the residual, S(L, x)−S log and S(L, x)−S RVB , of the fit of the entanglement entropy to Eqs. (27) and (28). In the right panel a comparison of the minimizing c L for both functional fits is shown. As is evident from the left-hand plot, the residuals are comparable for the two functional forms, with Eq. (28) giving a slightly better fit, using this metric. However, an important point to notice on the right-hand plot is that, when we attempt to fit the data for all system sizes L to the form Eq. (27), there is a clear system-size dependence in c L . In Ref. [21], using geometries amenable to DMRG simulation, it was argued that this increasing trend of c L may level off at a moderate system size value; this appears to support the notion of a central charge c = 1 in the limit of large system size. As is evident in Fig. 8, no such conclusion can be drawn from the toroidal QMC geometries for the second Rényi entropy -although, the QMC is not able to probe the behaviour for the von Neumann entropy, making a direct comparison difficult. In QMC there is no evidence that c L is converging to a constant value for Eq. (27), leading one to instead conclude that this cannot be a system-size independent (and hence universal) quantity in 2+1 dimensions The data is subtracted in such a way that only the log term should remain, with the line showing the best fit to this assumption. Note that the best fit of the entire data set gives a positive a 2 (π/2) and there is systematic curvature away from a log, that is to say the deviation of the data seems correlated in x rather than random and gaussian. here is defined as the (average of inner and outer) boundary length dividing the two regions.
for the Rényi entropy.
In contrast, for the case of Eq. (28) there is a much less-pronounced size dependence for the coefficient of J(y). As evident in Fig. 8, a much better case can be made that c L levels off to a system-size (ie. cutoff) independent value in the limit of large lattice sizes. This fuels speculation that J(y) may be a universal scaling function relevant for all fixed points in 2+1 dimensions, as discussed more in Section 5.

Polygons on a torus: entanglement due to vertices
The second geometry we examine is that of a rectangle embedded in a torus, as shown in Fig. 2. As discussed in Section 2.3, this geometry is particularly promising for accessing universal subleading corrections to the area law that may be computed in both continuum field theories and lattice models as in this work. Since these universal corrections arise due to vertices (or corners) in the entangled region, any geometry (or aspect ratio) of rectangle should have four separate vertex contributions, 4a 2 (π/2), from Eq. (9), assuming that each rectangle is large enough that interactions between corners may be neglected. In our QMC simulation geometries, the use of differently-shaped rectangles gives us more than one avenue to extract the coefficient of the subleading log, a 2 (π/2), thus providing an independent test for universality. Square, a 2 (π/2) = −0.006791 ± 0.000817 Small Square, a 2 (π/2) = −0.005306 ± 0.000645 Figure 10. The comparison of scaling an L/2 × L/2 square for two definitions of the boundary length. "Square" assumes the boundary is the average of the number of sites on the inside and the outside of the boundary, or simply the length of each edge times four. "Small Square" only counts the number of unique sites bordering in the inside edge of the boundary, which is smaller than the previous definition by four sites. By excluding smaller system sizes the extracted coefficient using these two methods converges, but the error in the extraction also increases.
There are two scalings that we test: Fig. 9 shows a square scaled in a fixed size simulation, while Fig. 10 shows an L/2 × L/2 rectangle scaled in an L × L system. Scaling of an L/2 × L/4 rectangle was also examined, and results are consistent with the L/2 × L/2 rectangle within error. We generate the data using the different geometries, then fit it to Eq. (9) with an additional allowance for a constant term in the fit: Repeating a previous analysis [20], we first perform a preliminary fit using a fixed lattice size, L = 36, and examine the Rényi entropy of all possible squares up to size L/2 × L/2. This data is fit to the form Eq. (9), and the scaling of the subleading logarithmic term as a function of boundary length is extracted. The result of this analysis is plotted in Fig. 9; we find that the constant is actually of the opposite sign from that seen in previous work [20] (reflecting a significant decrease in the error bars on S 2 with the present data set). If we exclude smaller rectangles from the fit, the sign of a 2 (π/2) remains opposite to previous work. In addition, if we look at the quality of the fit in Fig. 9, there is a large systematic deviation from a log( ) form, suggesting this geometry may be influenced by other factors not accounted for by this analysis, in particular, terms which depend on the aspect ratio of region A. In Fig. 10, we attempt another type of analysis aimed at eliminating the systematic shape-dependence observed in Fig. 9. There, the size of region A is fixed to be a square or rectangle with a perimeter proportional to the system's size L × L. Then, many different L×L toroidal lattice sizes are individually simulated (L = 10,12,14,16,20,24,28,32,36,40). Fig. 10, which shows the subleading logarithm dependence isolated from the area law (and additive constant), demonstrates that this analysis produces a very good fit to the expected function log( ). Analyses of two geometries (L/2 and L/4) of region A produce a consistent value for the universal coefficient of this logarithm, with the value extracted for squares of size L/2 × L/2 being more accurate (due to the availability of more system sizes). There are two additional pieces of analysis required for completeness: size exclusion and the definition of . For all results presented thus far on the polygonal entangling geometry, we have taken the definition of to be the average of the number of sites on the inside and outside of the boundary, except where specifically indicated otherwise. An alternate definition of the boundary length uses the minimum of the number of sites on the inside and outside of the boundary. These definitions of the boundary length only differ when considering our rectangular regions, differing by four lattice spacings (assuming a square of at least 2×2 in size). With very large sets of data, size exclusion can be used to eliminate bias when fitting data for which subleading corrections are not known (here, those proportional to 1/L). Since the number of available sizes in our study is limited, we can only do a limited exclusion study. The result of this analysis is that the corner term tends to become larger as smaller systems are excluded, and with smaller systems included the smaller definition of suggests a smaller (in magnitude) value of a 2 (π/2). It should also be noted that with a less comprehensive analysis, it is possible to get values with an artificially lower error bound, and part of the reason for this work is to illuminate the possible pitfalls in the extraction of these universal subleading terms.
All of this analysis suggests a universal corner contribution of a 2 (π/2) = −0.006(2), with caveats mentioned in the previous paragraph, for the TFIM at criticality, i.e. the universality class of the Wilson-Fisher fixed point. Remarkably, this value is very close to the value calculated in a continuum field theory at the non-interacting Gaussian fixed point by Casini and Huerta, a 2 (π/2) = −0.0064 [15,42]. Previous numerical series expansion studies of the interacting 2+1 dimensional QCP of the TFIM give −0.0055(5) [19] while Numerical Linked-Cluster Expansion (NLCE) gives −0.0053 [13]; both results are consistent, and lower than the result from the free field theory. Our QMC result gives independent confirmation of the validity of the series and NLCE results; however, due to statistical error bars, it is impossible for the QMC to distinguish between the non-interacting fixed-point value and the previous estimates for the interacting theory.

Discussion
Using a novel "projector" quantum Monte Carlo (QMC) method that operates at zero temperature [51], we have performed a detailed numerical study of the Rényi entanglement entropy, S 2 , at the quantum critical point of the transverse-field Ising model in 2+1 spacetime dimensions. We have focussed on two different entangling geometries amenable to large-scale simulation on toroidal lattices of size L × L.
First, when the entangling region is a square or rectangular polygon with four 90-degree corners, we confirm the expected scaling form of where A is the non-universal coefficient to the area (or boundary) law, and, a 2 (π/2) = −0.006 (2), is our best estimate of the universal coefficient of the subleading logarithmic term, which arises from each corner. It is a very non-trivial success that this value is within error bars of the value calculated at the same quantum critical point by numerical series and linked-cluster expansions [19,13], both of which take very different approaches to the thermodynamic limit. This value of the universal coefficient is also very close to the value calculated in a continuum field theory at the non-interacting Gaussian fixed point, a 2 (π/2) = −0.0064 [42]. The reason for this coincidence may be the fact that the Wilson-Fisher fixed point describing criticality in the TFIM may be reached perturbatively from the Gaussian fixed point. We hope that this fact motivates future analytical calculations of a α (θ) in continuum field theory at the interacting Wilson-Fisher fixed point. Second, we consider the lattice divided into two cylinders of size x × L and (L − x) × L, separated by two smooth (vertex-less) boundaries of length L. With this geometry, our data is consistent with the absence of any "additive logarithm", i.e. no logarithmic divergence depending on L/a. This is the same conclusion found in a previous calculation of free spinless fermions on finite-size square lattices, with π-flux through each plaquette [38]. In addition, no even-odd branching effect, like that observed for RVB-like phases [38,22], is seen in this geometry at the TFIM quantum critical point. Instead, we find excellent functional fit to where y = x/L is the aspect ratio of a cylinder, and J(y) (Eq. 8) is a function first derived for the quantum Lifshitz fixed point which describes certain Rokhsar-Kivelson (RK) Hamiltonians. We argue that this function gives qualitatively better fits to our QMC data, and additionally is consistent with a system-size independent coefficient c, which is not the conclusion one can draw if the data is fit to the familiar form derived for 1+1 CFTs, log(sin(πy)). It is interesting to speculate that this new function, J(y), may be a universal scaling function relevant for all fixed points in 2+1 dimensions, a conjecture that could be addressed with other numerical and field-theoretic studies. In particular, QMC simulations should have the ability to calculate the numerical value of the coefficient c at different interacting quantum critical points, as demonstrated by the current study. Finally, we stress the fact that the RVB-shape function J(y) is only applicable for Rényi entropies of order α ≥ 2, i.e. not the von Neumann entropy, thus emphasizing the practical utility of S 2 in the study of universality at quantum critical points. We hope that this work, which was obtained with a modest amount of CPU time (≈ 300 CPU-years), will motivate the calculation of similar quantities related to Rényi entanglement entropy at a variety of critical points in 2+1 and higher dimensions in the near future. In addition, we hope that this demonstration of the most convenient geometries for the study of entanglement in finite-size lattices via quantum Monte Carlo simulations will lead to fieldtheoretical studies of similar entanglement quantities at a variety of fixed points. If such a synergy could be accomplished between analytical theory and numerical simulation, the past successes in studying universality through entanglement in 1+1 may soon be translated to 2+1 and higher-dimensional quantum critical points.