Corner contribution to the entanglement entropy of an O(3) quantum critical point in 2+1 dimensions

0.64 MB·

by  Ann B. Kallin

Corner contribution to the entanglement entropy
of an O(3) quantum critical point in 2+1 dimensions
Ann B. Kallin,1 E.M. Stoudenmire,2 Paul Fendley,3 Rajiv R.P. Singh,4 and Roger G. Melko1,2
1Department of Physics and Astronomy, University of Waterloo, Ontario, N2L 3G1, Canada
2Perimeter Institute for Theoretical Physics, Waterloo, Ontario, N2L 2Y5, Canada
3Physics Department, University of Virginia, Charlottesville, VA 22904-4714 USA
4Physics Department, University of California, Davis, CA, 95616 USA
(Dated: February 21, 2014)
Theentanglemententropyforaquantumcriticalsystemacrossaboundarywithacornerexhibits
a subleading logarithmic scaling term with a scale-invariant coefficient. Using a Numerical Linked
ClusterExpansion,wecalculatethisuniversalquantityforasquare-latticebilayerHeisenbergmodel
4 at its quantum critical point. We find, for this 2+1 dimensional O(3) universality class, that it is
1
thricethevaluecalculatedpreviouslyfortheIsinguniversalityclass. Thisrelationgivessubstantial
0
evidencethatthiscoefficientprovidesameasureofthenumberofdegreesoffreedomofthetheory,
2
analogous to the central charge in a 1+1 dimensional conformal field theory.

b
e
F I. INTRODUCTION terms involving log((cid:96)/δ). With smooth boundaries, this
logarithmic divergence occurs only for odd d.19,20 For
9
1 Entanglementbetweentwosubregionsofasystempro- even d, it can occur when there is a singularity such as
vides a novel probe of quantum correlations.1 For this a corner or a cone in the surface separating A and B.
] probe to be useful, it is necessary to extract quanti- In d = 2, the case of interest here, the contribution to
l
e ties that are not only universal, but also give intuition γ of a corner of interior angle θ in the one-dimensional
r- into physical properties. One prominent idea is to use boundary is
t
s the entanglement entropy to define a measure of the de- (cid:18) (cid:19)
(cid:96)
t. grees of freedom at and near quantum critical points. γ =a(θ)log + . (2)
a Heuristically, such a measure quantifies the information δ ···
m
lost during RG transformations so that one may con-
Since the cutoff length δ is contained within the loga-
- strainflowsoftheoriesrelevantformanycondensedmat-
d rithm, rescaling it only affects the terms in the ellipses,
ter systems.2,3 While much is known for 1 + 1 dimen-
n not the coefficient a(θ).
o sional systems via the connection of Zamolodchikov’s c-
Such corner contributions in d = 2 have been com-
c theorem4 to entanglement,5–8 studies in higher dimen-
puted in several interesting situations. At a conformal
[ sionality are in their infancy.
quantum critical point21 such as the quantum Lifshitz
Subleadingtermsintheentanglemententropyprovide
2 theory (describing e.g. the square-lattice quantum dimer
v some very intriguing possibilities. The leading contribu-
model22), results from two-dimensional conformal field
4 tiontotheentanglemententropyS(A),betweentwosub-
theory23 have been adapted to give a(θ).24 It has been
0 regions A and B, scales with the ubiquitous area law9–11
5 (with a few important exceptions12–14). In a quantum calculated in free-scalar field theory,25 numerically in
3 interacting lattice models,26–29 and more generally us-
critical system in d+1 space-time dimensions, a shape-
. ing the AdS/CFT correspondence.30 These are consis-
1 dependent subleading term γ is expected, so that15
tent with a conjectured geometrical form valid in any
0
4 (cid:18)(cid:96)(cid:19)d−1 dimension.31
1 S(A)=C + ... + γ + . (1) These computations confirm that a(θ) is indeed uni-
δ ···
: versal. Even more strikingly, they point to a(θ) as a
v
i Thelineardimension(cid:96)characterizesthesubregionAand useful measure of the number of degrees of freedom. For
X δ is the ultraviolet/lattice-scale cutoff, while the ellipses example,fora2+1-dimensionalconformalquantumcrit-
r concealnon-universalconstantsandsubleadingtermsde- ical point, a(θ) is proportional to the central charge c of
a
pending on the length scale to some power. Explicit ex- thetwo-dimensionalconformalfieldtheorydescribingthe
pressions for γ in a cylindrical geometry with smooth ground-state wavefunction.24
boundaries were found for φ4 theory in d = 3 (cid:15).16 The purpose of this paper is to show that this corner
−
At the infrared-unstable free-field fixed point, γ indeed contribution behaves similarly in a less exotic situation,
decreases under the renormalization group flow toward namely an interacting theory familiar to both condensed
the non-trivial Wilson-Fisher fixed point, just like the matter and particle theorists. To this end, we study the
c-theorem requires for the central charge c in d=1.17 entanglement entropy for a critical system of spin-1/2
The precise functional form of γ is expected in gen- particlesonasquare-latticebilayerwithnearest-neighbor
eral to depend on scale invariants such as the aspect ra- antiferromagnetic Heisenberg interactions. This quan-
tios,Eulercharacteristic,andothergeometricfeaturesof tum critical point is in the same universality class as
region A.18 A particularly interesting piece comes from the three-dimensional classical Heisenberg model,32 de-
2
scribed in the continuum limit by the well-studied three- represented pictorially in Fig. 1(c) for a 4 4 2 bilayer
× ×
component O(3)-invariant φ4 field theory.33 system.
Specifically, we numerically compute the universal co- This Hamiltonian supports two different phases in its
efficient a(θ = π/2) in a continuous family of Renyi en- ground-state. At J = 0, the physics is simply that of
tanglement entropies for this critical Heisenberg bilayer. two decoupled squa⊥re-lattice Heisenberg models. These
We compare the results to those for the critical point are N´eel ordered, so that the SU(2) symmetry is spon-
ofthetransverse-fieldIsingmodel(TFIM),28 whosefield taneously broken to U(1). For small J /J, the order
theory description is a single-component φ4 theory. We persists, and the coupling between the la⊥yers relates the
find over a wide range of Renyi index values, a(π/2) for ordering. Thus there are two Goldstone bosons in this
the former is thrice that for the latter, to within numer- N´eel phase (which, incidentally, give rise to a sublead-
ical uncertainties. The factor of three is compelling evi- ing logarithm in the entanglement entropy of a straight
dencethattheuniversalcoefficientofthecorner-induced boundary on a finite-size lattice26,42). At large J /J,
logarithm provides a measure of the number of low-lying there is a “dimer” phase where the pair of spins on⊥each
degreesoffreedom. Thisindeedisbehavioranalogousto bond between the bilayers forms a singlet. In this phase
that of c in 1+1 dimensional CFTs. the SU(2) symmetry is unbroken, and the spectrum is
In an interacting theory in general d, one does not gapped. A single quantum critical point separates the
expect such a measure to literally count the degrees of twophases. Estimatesforthecriticalcouplinghavebeen
freedom;rather,itprovidesawayofunderstandingwhich calculated to high accuracy with series expansion43,44
renormalization-group flows are possible. The factor of and quantum Monte Carlo,32 with the most accurate es-
3 we find is presumably a consequence of the fact that timate (J /J)c =2.5220(2) coming from the latter.
O(N)φ4 theoryis“close”tofree,inthatthe(cid:15)expansion The un⊥iversality class of this phase transition turns
around a free-field theory applies. To be precise, our out to be the simplest possibility consistent with
numerical result implies that the leading contribution to the symmetry. The order parameter describing the
a(θ) in N-component φ4 field theory in 2+1 dimensions N´eel/Goldstone phase is the staggered magnetization.
isproportionaltoN, andthatcorrectionsarewithinour Following the Landau-Ginzburg approach, one defines a
(fairly small) numerical error. three-component vector field φ(cid:126) representing a suitably
The outline of this paper is as follows. In section II averaged order parameter. The simplest action for this
we review the Heisenberg bilayer model and its theoret- fieldconsistentwiththesymmetriesistheO(3)-invariant
ical description. In section III, we describe our numeri- φ4 theory:
calmethod,anovelNumericalLinked-ClusterExpansion
(cid:32) (cid:33)
(NLCE).34–37 Our NLCE uses both Lanczos diagonal- (cid:90) ∂φ(cid:126) ∂φ(cid:126)
S = d2xdt φ(cid:126) φ(cid:126) µ2φ(cid:126) φ(cid:126) g(φ(cid:126) φ(cid:126))2 .
ization and Density Matrix Renormalization Group38–40 ∂t · ∂t −∇ ·∇ − · − ·
(DMRG) simulations to calculate the Renyi entangle-
(4)
ment entropies.41 We present our results in section IV, At µ2 < 0 the O(3) symmetry is spontaneously broken
and discuss some implications of this work in section V.
to SO(2), resulting in two Goldstone bosons. There is a
continuousphasetransitionatµ2 =0toaphasewithno
order in φ(cid:126), as in the bilayer.
II. THE HEISENBERG BILAYER MODEL AND
It has long been known that this Landau-
ITS CRITICAL POINT
Ginzburg/effective field theory describes the phase
transition in the three-dimensional classical Heisenberg
The aim of this paper is to understand the universal magnet between the low-temperature Goldstone and
subleadingtermcomingfromcornercontributionstothe the high-temperature disordered phases. Even though
entanglement entropy. We study a strongly interacting the underlying degrees are fixed-length spins (and so
quantum model whose physics is well understood, and are labeled by two angles), the average value of the
which is amenable to treatment by the powerful NLCE magnetization also fluctuates strongly near the critical
method described in section III. This system is the spin- point, and so this three-component theory describes the
1/2 nearest-neighbor Heisenberg model on a square lat- critical behavior. This has been confirmedby comparing
tice bilayer, or equivalently, two flavors of spins on the detailed numerical simulations33,45 with calculations in
square lattice. Labelling the spin operator on site j and the d=3 (cid:15) expansion.46
layer a by Saj, the Hamiltonian is The sam−e Landau-Ginzburg approach applies to the
Heisenberg bilayer model we study here, with the same
(cid:88) (cid:88)
H =J (S S +S S )+J S S , (3) conclusion. This has also been confirmed convincingly
1i 1j 2i 2j 1i 2i
i,j · · ⊥ i · via numerics,32 so that critical fluctuations correspond-
(cid:104) (cid:105)
ing to the restoration of SU(2) symmetry at the crit-
where i,j denotesapairofnearest-neighborsiteswithin ical point indeed are described by three bosonic fields.
(cid:104) (cid:105)
a layer. We take the couplings J and J as positive, The appearance of the longitudinal mode can also be
hence antiferromagnetic. The interactions⊥between sites understood directly in the quantum theory.47 It is a
onthesamelayer(J)andsitesondifferentlayers(J )are bound state in the Ising limit and therefore neglected by
⊥
3
in Ref. 28, a modified NLCE procedure was defined that
(a)
employs an alternative definition of cluster geometries,
involving only m n rectangles. This restriction
×
significantly simplifies the cluster embedding problem,
passing the computational bottleneck to the calculation
of ground-state properties, via numerical techniques
such as exact diagonalization (Lanczos) or DMRG. As a
(b) (c) J consequence, one is able to achieve significantly higher
6 2
J orders of the expansion than conventional NLCE, and
7 ⊥ therefore significantly improve approximations to the
thermodynamic limit. In this section, we give details of
2 4 3 the NLCE procedure, definition of cluster geometries for
the Heisenberg bilayer model, methods to define order
extrapolations, and procedures to calculate the Renyi
FIG. 1. (a) The lattice constant L(c) denotes the number entanglement entropies using the Lanczos and DMRG
of distinct ways a cluster can be embedded in a lattice. The methods.
boxesillustratethefourpossiblevaluesofL(c)=1,2,4,8for
asquarelattice. (b)Thesubclustersofa2×3cluster,corre-
spondingtorow(6)ofTableI.(c)A4×4×2bilayercluster A. NLCE Overview
(row (10) in Table I) with intralayer coupling J denoted by
the thin grey bonds, and interlayer coupling J denoted by
⊥
The foundation of the NLCE method is based on the
the red vertical bonds.
fact that properties of a lattice model can be expressed
as a sum over all distinct clusters which are embeddable
spin-waveandSchwingerbosontreatments,whichamong in the lattice . Let a cluster c be a set of sites with
L
other things leads to very poor estimates for the critical the connectivity of the underlying lattice, and L(c) be
coupling(J /J) . Thissituationcanberemediedwitha the number of distinct ways it can be embedded. A ro-
c
propertrea⊥tmentofallthreemodes48–50). Theexistence tation or reflection could result in a different embedding
of this additional longitudinal mode has been recently dependingonthesymmetryofc,whereasasimpletrans-
emphasized by a numerical study explicitly demonstrat- lationwillnotleadtoadifferentembeddingintheinfinite
ing that it becomes degenerate with the two Goldstone lattice ; see Figure 1(a)). For a square lattice, L(c) can
L
modes exactly at the critical point, corresponding to a take values 1, 2, 4 and 8. Here, we follow Ref. 28 and
full restoration of SU(2) symmetry.44 consider only m n rectangular clusters, for which L(c)
×
Thus, the quantum critical point of the Heisenberg bi- can only be equal to 1 (if m=n) or 2 (if m=n).
(cid:54)
layer model offers an excellent opportunity to examine A property P per site can then be expressed as
thebehavioroftheentanglemententropyinthepresence (cid:88)
P( )/N = L(c) W(c) , (5)
of multiple bosonic modes in the low-lying spectrum of
L ×
an interacting model.In the next section, we discuss the c
details of the numerical simulation scheme used to ex- where the weight of a cluster for a given P is defined by
tract the corner contribution to the Renyi entanglement
(cid:88)
entropies. W(c)=P(c) W(s), (6)
−
s c
∈
where s is any subcluster of c. In the sum, each subclus-
III. NUMERICAL LINKED-CLUSTER tersisincludedthenumberoftimesthatitcanfitinside
EXPANSION cluster c (see Figure 1(b)). This relation is a generaliza-
tion of the inclusion-exclusion principle, the statement
The Numerical Linked-Cluster Expansion34–37 that the number of elements in the union of two finite
(NLCE) is a method of extending measurements of a sets A B is simply A + B A B . This gives in-
| ∪ | | | | |−| ∩ |
series of finite-sized lattice clusters towards the thermo- tuitionforwhydisconnectedclustersneednotbeconsid-
dynamic limit, cancelling off finite-size and boundary ered in NLCE calculations: their weight W(c) is always
effects using sums and differences of various clusters. In zero, just as the number of elements in the intersection
its general formulation, this method uses measurements of disconnected sets is zero.
of a given property P from all possible clusters of sites TheNLCEprocedureusesEqs.5and6tobuildupthe
that can be embedded in the chosen lattice. Typical value of P/N, starting from the smallest cluster (which
NLCEapproachesinvolvethecomputationallyexpensive has no subclusters, thus W(1) = P(1)) and ending with
task of generating all clusters, and embeddable sub- some maximal cluster size. For rectangular clusters on
clusters. This cluster embedding problem results in an a square lattice, this maximal cluster size is limited in
exponentialbottleneck,restrictingthemaximumnumber practice only by the computational cost of calculating
ofsitesto 16orso,dependingonthelattice. However, the given property P on the cluster.28
∼
4
ThevalueoftheweightW(c)isanindicatorofthecon- Subcluster Multiplicity
vergence of the NLCE. In a system with no broken sym- id# n ×n N L(c) (1) (2) (3) (4) (5) (6) (7) (8) (9)
x y
metriesandafinitecorrelationlength,theweightsshould
(1) 1×1 2 1
decrease exponentially with cluster size, once these sizes
(2) 1×2 4 2 2
exceed the correlation length. This would lead to an ex-
ponential convergence of the NLCE with cluster size (or (3) 1×3 6 2 3 2
“order” as we define below). At (and near) a critical (4) 1×4 8 2 4 3 2
point, where the correlation length becomes large com- (5) 2×2 8 1 4 4
pared to the sizes of clusters we can study, the weights
(6) 2×3 12 2 6 7 2 2
will vary as a power of the cluster size. This will lead
(7) 2×4 16 2 8 10 4 2 3 2
toanalgebraicconvergencewithorderforquantitieslike
(8) 3×3 18 1 9 12 6 4 4
ground-state energy and entanglement entropies, requir-
(9) 3×4 24 2 12 17 10 3 6 7 2 2
ing a careful extrapolation. Note that other quantities,
suchasorder-parametersusceptibilitywilldivergeatthe (10) 4×4 32 1 16 24 16 8 9 12 6 4 4
critical point. But even for divergent properties, the
TABLE I. The properties of all clusters up to a maximum
NLCE can extract useful information if one can reach
n and n of four sites, including: number of sites N, lattice
orderslargeenoughforthepropertytobefittoaknown x y
constantL(c),andthemultiplicityoftheirsubclusters(thisis
scaling relation. The fact that some quantities converge
left blank when the multiplicity is zero). There is an implied
and others diverge is analogous to convergence or diver-
bilayer (n =2) for each cluster.
z
gence of a series expansion at its radius of convergence.
It is instructive to work out the NLCE for one-
dimensionalsystems,whichcanbedoneinfullgenerality.
B. Rectangular Clusters for the Square Lattice
Allsub-clustersofaone-dimensionalsystemareuniquely
Bilayer
labeledbytheirlengthn,andcanbeembeddedonlyone
way, thus L(n) 1 and
≡ In this paper we study a square lattice bilayer system.
Each cluster used in the calculation is an n n 2
(cid:88)∞ x × y ×
P/N = W(n) . (7) array of sites, where n ,n 1. The second layer of the
x y
≥
n=1 bilayer lattice does not change the NLCE calculations
fromthestrictlytwodimensional(2d)NLCE,otherthan
The cluster weights are
doublingthenumberofsitesinagivencluster,increasing
thecomputationalexpenseofmeasuringmostproperties.
W(1)=P(1) (8)
Table I shows the characteristics of all rectangular clus-
W(2)=P(2) 2P(1) (9)
− terswithnx,ny 4,includingthenumberofembeddings
W(3)=P(3) 2P(2)+P(1) (10) ofeachsubcluste≤r,neededforEq.(6). Unfortunatelythe
−
2d cluster weights do not simplify as in the 1d case, and
···
all terms (even those from the smaller clusters) need to
W(n)=P(n) 2P(n 1)+P(n 2) [n 3]. (11)
− − − ≥ be included in Eq. (5).
Defining partial sums up to order n of the NLCE as InanNLCEcalculation,itisnecessarytodefineanor-
der, a way of grouping together clusters of similar sizes.
n
(cid:88) This allows one to assign a length scale related to the
p(n)= W(m) , (12)
largest order included in a calculation. At a quantum
m=1
criticalpointforexample,thisallowsonetostudyNLCE
then for n>1 data as a function of order, giving scaling relationships
that can be extrapolated towards the thermodynamic
p(n)=P(n) P(n 1). (13) limit.
− −
In this paper we consider two methods of grouping
Assuming that P(n) changes ever more slowly with in-
clusters: (1)bytheaverageedge-lengthofacluster =
creasing n, then in the limit of large clusters O1
(n +n )/2 and (2) by the number of sites = N/2
x y 2
O
∂2 in one layer of the cluster. In the two cases we have
W(n) P(n) (14) the average cluster length, (cid:96) = (n +n )/2 and
(cid:39) ∂n2 ∼ O1 x y
∂ (cid:96) √ 2 =√nxny,whichcanbethoughtofrespectively
∼ O
p(n) P(n). (15) asthearithmeticandgeometricmeanoftheedge-lengths
(cid:39) ∂n
for clusters of that order; Fig. 2 attempts to give some
It is interesting to observe from Eq. (13) that in 1d, the intuitionforthetwodifferentdefinitions. Using tends
2
O
NLCE is nothing but the “subtraction trick” used, for toincludelonger1dclusters,whileonlyincludingsmaller
example,inDMRGcalculationstoestimatebulkproper- squareclusters, whereas truncatingusing includesall
1
O
tiesinthethermodynamiclimitfromfinitesystemswith the clusters that will fit inside a diamond of a given size
open boundary conditions.40 defined by that order, thus excluding long 1d clusters.
5
FIG. 3. (a) The entanglement due to a corner is calcu-
latedbyconsideringeachplaquetteofacluster. Inthisfigure
we consider a 4×3 cluster. (b) The corner entanglement is
equal to the difference of the entanglement due to opposing
boundaries with corners (V1, V2) at the given plaquette and
the entanglement across to the corresponding horizontal and
vertical lines (Lx, Ly). (c) This figure shows all four entan-
FIG. 2. Two methods for defining order O. Each tile shows
glement cut geometries for the first plaquette superimposed.
thenumberofsitescontainedinarectanglewithitstopright
The subtracted line cuts are solid, and the corner cuts are
cornerinthatsquareanditsbottomleftcorneratthebottom
dotted.
leftofthefigure. Asanexamplethe2×3rectangleisoutlined
withadashedline,andfromthetoprightcorneronecansee
thatitcontains6sitesperlayer. UsingO thedifferentorders
1 C. The Renyi Entanglement Entropies
aredefinedbythediagonalsofthisdiagram. O =1contains
1
onlythe1×1×2cluster(denotedby“1”inthediagram),O =
1
1.5 contains the clusters denoted by the 2’s, O =2 contains The Renyi entanglement entropies41 are used to mea-
1
3,4,3andsoforth,movingalongthediagonals. O issimpleto sure bipartite entanglement between two spatial regions
2
understandfromthisdiagram. IfO2 =x,thatordercontains of a system, labelled A and B,
clusterswithxsitesperlayer. Theshadedtilesshowclusters
solved by Lanczos and DMRG in this work. 1
S (A)= logTr(ρα), (16)
α 1 α A
−
where ρ = Tr (ρ) is the reduced density matrix of re-
A B
gion A. Here, α is the Renyi index, such that the limit
α=1 gives the familiar von Neumann entropy.
Additionally, each order using includes an increasing
1
O As mentioned above, the measurement of Renyi en-
number of clusters (equal to rounded down to the
1
O tropies in NLCE is non-standard because it requires the
nearest integer), while the number of clusters for each
definition of two spatial regions A and B. In this paper
new order of is generally smaller and determined by
2
O we focus on the entanglement due to a 90 corner in the
the number of factors of . This tends to give the data ◦
2
O boundary between entangled regions. This is extracted
a step-like distribution as a function of while allow-
2
O by defining a boundary with a corner, and subtracting
ing for higher orders to be calculated, resulting in more
offtheentanglementcontributioncomingfromthelinear
data points. Contrastingly, plotting data as a function
portions,leavingonlytheentanglementduetothecorner
of gives a smoother distribution with fewer overall
1
O (Fig. 3).
datapoints. Bothmethodscontainthesameinformation
Entanglement measurements in general are done by
(though excludes long 1d clusters) simply viewed in
1
O consideringeverypossiblewaythatthechosenboundary
different ways.
can intersect a cluster. Thus, to measure the entangle-
ment across a line running in the y direction Ly, for a
The lattice constant, as described in Section IIIA, re-
given cluster we use
flects the number of ways a cluster can be embedded
in the underlying lattice, and in the case of our square
bilayer lattice, is limited to the values L(c) = 1,2. It P(c)=n(cid:88)x−1S (Ly), (17)
α i
is important that measurements considered on a cluster
i=1
obeythesamesymmetry,i.e.ifL(c)=2thenameasure-
ment must give the same results for both orientations where Ly denotes a region A including i columns of the
i
of the cluster, otherwise the two orientations should be cluster. This measurement can be simplified for many
considered separately with L(c) = 1 for each. Standard rectangular clusters. Since S(A) = S(B) for a pure
measurements, such as energy, would never depend on ground-state wavefunction, we have S(A ) = S(A )
the orientation of the cluster. However we will see below which can approximately halve the numbier of meansxu−rei-
thatthisbecomesimportantforsimplifyingthemeasure- ments required. Different types of boundaries will have
ment of entanglement, which depends on spatial bound- different symmetries that can be exploited to reduce the
aries defined in the lattice. total number of measurements required.
6
AsmentionedinSection??,ifaclusterhaslatticecon- (a) (b)
stant L(c) > 1, any measurements on that cluster must
B
share the same symmetry, or else L(c) must be modified
for that cluster. This becomes especially important to
considerforentanglementmeasurements,whereacluster
isdividedintotwospatialregions. Theexampleabove,of A
the entanglement across a vertical line will give different
results for the two orientations of a non-square cluster. FIG.4. Atypicalpath(a)usedwithinDMRGforcomputing
This can be remedied by instead calculating the sum of theground-stateofthe5×4×2cluster. Hereweshowonlythe
theentanglementduetoaverticallineandthatduetoa toplayerofsites;thepathvisitseachbottom-layersiteinturn
horizontal line. A second, but equivalent, method would before going to the next top-layer site. An irregular path (b)
be to always treat n m clusters as distinct from m n is needed to obtain the corner cut dividing the system into
× × regions A and B as shown.
clusters and measure entanglement across only a line in
thexdirection,forexample. Bothtechniquesrequirethe
same computational effort, but use different methods of
where a , a are numerical coefficients such that
bookkeeping for the clusters and measurements. i jk
(cid:80) a 2 = (cid:80) a 2 = 1, ψi are basis vectors for the
The NLCE, unlike many other numerical techniques, | i| | jk| { c}
is able to isolate the entanglement due to a 90◦ corner full cluster of both regions A and B, and {ψcAj} and
by subtracting off the entanglement contributions from ψBk are the basis vectors of region A and B respec-
{ c }
the linear portions of the boundary. This is done by tively. Computationally, this means constructing a ma-
subtractingtheentanglementduetothelinesLx andLy, trix with all region A basis states as the rows and all
fromentanglementacrossopposinglineswith90 vertices region B basis states are the columns. Then, running
◦
(V1 and V2), as shown in Figure 3(b) and (c). To fully through Ψc , each entry is assigned to an element of
| (cid:105)
cancel off any line contributions to the entanglement we Mc where ψci = ψcAj ⊗ψcBk. From there, obtaining the
calculate thecorner entanglement α of acluster c using reduced density matrix simply requires multiplying this
V
matrix by its conjugate transpose,
2Vα(c)=(cid:80)ni=x1−1(cid:80)nj=y−11(Sα(Vi1,j)+Sα(Vi2,j)) (cid:88)(cid:88)
−(ny−1)(cid:80)ni=x1−1Sα(Lyi) (18) ρA =McMc† = aija∗kl|ψcAi(cid:105)(cid:104)ψcBj|ψcBl(cid:105)(cid:104)ψcAk|
whereVi1,j denotesa−re(gnixon−A1)in(cid:80)clnju=yd−1in1gSαa(nLixj×),j-siterect- =(cid:88)i,k (cid:88)j aijai,∗kjj(cid:104)kψ,lcBj|ψcBj(cid:105)|ψcAi(cid:105)(cid:104)ψcAk|. (20)
angle of cluster c. It is more intuitive to say that we
measurethefourtermsinFigure3(b)foreachofthepla- ThemultiplicationMcMc† =ρA orMc†Mc =ρB ischosen
quette in 3(a), which easily extends to clusters of other basedonwhichwillresultinthesmallerreduceddensity
shapesandsizes. Inpractice,duetothesymmetryofthe matrix, since one final diagonalization must be done to
rectangular clusters all V2 measurements will already be extract its eigenvalues. The diagonalization of the re-
done by V1. The number of measurements required can duced density matrix must give the full eigenvalue spec-
be further reduced for a square cluster, where V =V trum; this requires a computationally more expensive
i,j j,i
for both V1 and V2. algorithm than Lanczos which only returns the largest
eigenvalues. The above method is limited by computer
memoryandtime,anditissuitableonlyforsmallerclus-
D. NLCE cluster solvers ters. It is imperative to use a good linear algebra library
asitwillsignificantlyimprovethespeedandperformance
1. Lanczos & Full Diagonalization ofthisalgorithm. ForthisworktheEigenC++template
library51 was used to solve clusters of up to 30 sites with
Lanczos and exact diagonalization.
It has not yet been discussed how the measurements
onthedifferentclustersareobtained. InpreviousNLCE
studies, generally the Lanczos algorithm for diagonaliza-
2. Density Matrix Renormalization Group
tion was used obtain the full ground-state wavefunction
ofthecluster, Ψ . Toextracttheentanglemententropy
c
foreachofthe|diff(cid:105)erententanglementcutgeometries(see For larger clusters we use a complementary method,
Fig. 3 for examples) one must take a partial trace of the the density matrix renormalization group (DMRG).52
fulldensitymatrixρ= Ψ Ψ overthestatesinregion Unlike exact diagonalization and Lanczos, DMRG is not
c c
B. The partial trace is| do(cid:105)(cid:104)ne b|y rewriting the ground- necessarily limited by the total number of sites in the
state vector in matrix form, cluster. Instead, DMRG traverses the system along a
one-dimensionalpath,suchastheoneshowninFig.4(a)
(cid:88) (cid:88)
Ψ = a ψi M = a ψAj ψBk , (19) forthe5 4 2cluster(top-downview). Forapathofthis
| c(cid:105) i| c(cid:105)→ c jk| c (cid:105)(cid:104) c | × ×
type, DMRG scales exponentially in n but scales much
i j,k y
7
more favorably in n (e.g. linearly in n if the system is
x x -1.12
gapped). NLC-
1
A key advantage of using DMRG for computing en- -1.121 NLC-O2
tanglement entropy is that it diagonalizes the reduced SerOies
QMC
density matrix for a different bipartition of the system -1.122
at every step. These bipartitions correspond to cut- e
t
ting the one-dimensional path at each bond. Thus one /si -1.123
y
can obtain the full entanglement spectrum for many of g
r -1.124
the cuts needed for the corner entanglement contribu- ne
E
tion within a single DMRG calculation. For example,
-1.125
the path in Fig. 4(a) provides two inequivalent vertical-
line cuts (on the 4th and 8th bonds of the DMRG path) -1.126
as well as corner cuts separating sites in the first col-
umn from the rest of the system. Other corner cuts, -1.127
such as the one separating the system into the regions A 0 0.1 0.2 0.3 0.4 0.5
and B of Fig. 4(b), require modifying the DMRG path. 1/‘
Suchirregularpathstypicallyincreasethecomputational
cost needed for DMRG to reached a fixed accuracy since
FIG.5. TheenergypersitefortheHeisenbergbilayersystem
short-rangeinteractionsinthetwo-dimensionalHamilto-
at its quantum critical point as a function of 1/(cid:96) along with
nian get mapped to much longer-ranged interactions in fits to b/(cid:96)3 +d for some constants b and d. The horizontal
the one-dimensional model seen by DMRG. lines are predictions for the ground-state energy in the ther-
For the results below, we used DMRG to solve the modynamic limit from series expansion and quantum Monte
4 4 2, 3 6 2, 6 3 2, 4 5 2, and 5 4 2 clus- Carlo.53
× × × × × × × × × ×
ters. Each calculation kept up to 10,000 states or
enough states to obtain a truncation error below 10 12,
−
whichever occurred first. The only exception was one of that the two critical points are “close”, in harmony with
the irregular paths for the 4 5 2 cluster for which we thefactthatcriticalpropertiesoftheWilson-Fisherfixed
× ×
only managed 8000 states (due to memory constraints) point can be computed very accurately in the d = 3 (cid:15)
−
foratruncationerrorof1.4 10 9,whichisstillquiteac- expansion around a free-field theory.46
−
×
curate. One can easily benchmark the accuracy for such AsdiscussedinSectionII,theHeisenbergbilayerquan-
difficult clusters by comparing the energy to the same tumcriticalpointisdescribedbythethree-componentφ4
cluster studied with a path more favorable to DMRG. theory. Thus if a (π/2) provides a measure of the num-
α
For the 4 5 2 cluster with the difficult irregular path, berofdegreesoffreedomanalogoustoc,itshouldbeap-
× ×
for example, we still obtained the energy within a rela- proximately three times the value obtained for the Ising
tive error of 10−9 compared to the 5 4 2 result, which theory. Itneednot(andoughtnot)beexactlythricethe
× ×
is essentially exact to numerical precision. value, sincethisisaninteractingfixedpoint, andcritical
properties depend on the number of components N in a
non-trivial way. Nevertheless, the O(N) Wilson-Fisher
IV. RESULTS fixed point can be reached by perturbing around N free
fields,andsooneexpectsa (π/2)tobeclosetothefree-
α
Inthissectionwepresentourresultsforthesubleading field value, roughly N times the single-component value.
logarithmic term in the Renyi entropy resulting from a As we now describe, our numerical calculations do in-
corner of angle π/2 at the quantum critical point in the deed yield this factor of 3. We perform the NLCE pro-
Heisenberg bilayer. We compare these to the analogous cedure described in the previous section, at the critical
resultsforthequantumcriticalpointinthe2dtransverse- coupling(J /J) =2.5220oftheHeisenbergbilayersys-
c
fieldIsingmodel(TFIM)onthesquarelattice,28 inorder tem,Eq.(3)⊥. Infigure5theNLCEincludesclustersfrom
to study the relation between the two. 1 1 2to4 5 2,whilefigs.6and7excludethe1dclus-
TheTFIMhasonlyaZ spin-flipsymmetry,soapply- te×rs×which d×o n×ot contribute to the corner entanglement
2
ing the Landau-Ginzburg approach to it yields a scalar term(seeEq.(18)). ThelargestNLCEtruncationorders
field theory with a single component. This has the same included are =4.5 and =20.
1 2
O O
O(N) φ4 action given in Eq. (4) for the Heisenberg bi- Before discussing results for the Renyi entropies, we
layer, except for Ising, φ(cid:126) has a single component. The performaninitialcheckofourNLCEprocedure,byusing
critical theory in this N = 1 case is the famed Wilson- it to compute the ground-state energy per site in the
Fisher fixed point, in the same universality class of the Heisenberg bilayer at its critical point. At T = 0, the
classical 3 dimensional Ising model. This theory is not ground-state energy plays the role of free energy, and by
free;theNLCEcalculation28 givesvaluesofa (π/2)very hyperscaling, its singular piece scales as ξ (d+z), where
α −
closeto,butnotexactly,thosecalculatedforRenyiindex forourproblemd+z =3. Thus,weexpectthecorrection
α=1,2,3 in Gaussian free-field theory.25 This indicates at order (cid:96) to scale as 1/(cid:96)3. In Fig. 5, the ground-state
8
0
0.016
0.01
− Gaussian
0.02 2 0.014
V TFIM
−
−0.03 V1.5 0.012 Bilayer O1/3
0.04
Bilayer /3
αV−0.05 V1.25 aα−0.010 O2
−
0.06 0.008
− V1.125
0.07
− 0.006
0.08
− 1
V 0.004
0.09
− 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.0 1.5 2.0 2.5 3.0 3.5
α
log‘
FIG. 6. The entanglement due to a single corner Vα in the FIG.7. Thecoefficient−aα ofthelogarithmduetothepres-
Heisenberg bilayer system for α = 1,1.125,1.25,1.5,2 using enceofa90◦ corner. GreencirclesaretheGaussianfreefield
both definitions of order O1 (circles) and O2 (squares) along calculation.25 The red dashed line is for the transverse-field
withfitstoVα =aαlog(cid:96)+bα. Theresultingcoefficientaα is Ising model.28 The solid lines are results for the Heisenberg
shown in figure 7 for both definitions of order. bilayer divided by 3, with the NLCE fit using the two differ-
ent definitions of the order O and O . The shaded regions
1 2
aroundeachlinecorrespondtoanestimateoftheerrorinthe
energy per site is plotted as a function of 1/(cid:96) for both data, discussed in Section IV.
definitions of the order: =(n +n )/2=(cid:96) and =
1 x y 2
O O
n n = (cid:96)2. Each dataset is fit to the function E ((cid:96)) =
x y 0
1/(cid:96)3+E0(∞),whereE0(∞)isthepredictedground-state the cluster length scale (cid:96), to extract aα and bα. We note
energy per site in the thermodynamic limit. For the two that using results in a systematically higher value of
1
different fits we find, b , thoughOthe difference decreases as α increases.
α
Using this fitting procedure, our results for the log-
:E ( )= 1.12665
O1 0 ∞ − coefficient a (π/2) are plotted in Fig. 7. Here, the coef-
α
:E ( )= 1.12649.
2 0 ficient a is compared between three different theories:
O ∞ − α
−
the single-component φ4 theory via the TFIM, the free
Even though relatively small cluster sizes are included
field theory calculation of Ref. 25, and the present cal-
in this extrapolation, we see that both values of E ( )
0
∞ culation of the Heisenberg bilayer. In Fig. 7, data for
areveryclosetotwoindependentcalculationsfromcom-
the Heisenberg bilayer is divided by 3 to emphasize how
plementary, but very different techniques (see Fig. 5).
remarkably close it is to thrice the Ising data. The im-
First, from series expansions, Pade extrapolations lead
plications of this are discussed in detail in section V.
to a value of E ( ) = 1.1262 at (J /J) = 2.5220.
0 c
Second, from large∞-scale u−nbiased quant⊥um Monte Carlo The error bars shown in Fig. 7 are meant to be a
(QMC)SenandSandvik53 revealahighlyaccuratevalue guide to the reader. They are calculated as the stan-
darddeviationofthedatafromthelinearfitsto ((cid:96))=
ofE ( )= 1.1265201(5)atthequantumcriticalpoint, α
0 V
∞ − a log(cid:96)+b ,examplesofwhichareshowninFig.6. This
which is again consistent with the NLCE results. α α
error is then assigned to a , although strictly speaking
We now turn to our NLCE calculation for the sub- α
it also depends on b . As with any study of this kind
leading scaling term γ induced by a 90◦ corner in the α
which incorporates functional extrapolations using rela-
entanglement boundary. As discussed in the last sec-
tively small cluster sizes, significant uncertainty related
tion, the NLCE can isolate the corner contribution to
totheprecisedataseriesincludedinthefitremains, and
the Renyi entanglement entropy independently from the
is not represented by the error bars in Fig. 7. It is worth
leading “area law” contribution to scaling. Thus for
noting that the NLCE results28 for the second Renyi
each value of the Renyi index α, we extract the value
entropy S (A) in the TFIM were independently bench-
of a (π/2) directly from fits of this corner entropy to 2
α
marked against series-expansion27 and QMC data29,54
((cid:96))=a log(cid:96)+b . (21) obtained through a replica-trick procedure. Both cal-
α α α
V culations yield a coefficient a(π/2) consistent with the
The raw NLCE data for this quantity is shown in Fig. 6, NLCE to within numerical errors. The unbiased QMC
for several values of α. Data is plotted separately for data were obtained from a much different fitting proce-
bothdefinitionsoforder, and ,andseparatefitsare dure involving a square subregion A with four corners
1 2
O O
performedforeachvalueofαtoEq.(21),asafunctionof embedded in a toroidal lattice; thus the match with the
9
NLCE is particularly striking. in higher dimensions also suggest this behavior. Using
It is clear from Fig. 7 that some uncertainty remains the AdS/CFT correspondence, the corner-induced loga-
regarding the relationship between the data for α(cid:46)1.3, rithmfromanangleθanditsanalogsinarbitrarydimen-
obfuscated by the growth in uncertainties in the NLCE sions obey the scaling form31
calculation in this regime. One can see from Fig. 6 that
the data points deviate further above and below the lin- (cid:18) L (cid:19)d (cid:18)(cid:96)(cid:19)
γ = q(θ)log + ..., (22)
earfitforsmallervaluesofα,directlyincreasingtheerror L δ
P
bars. ThepoorernumericalconvergenceoftheNLCEfor
smallerαcouldbeexpected,particularlysincetheα<1 where L is the length scale set by the AdS curvature
Renyi entropies are more sensitive to the tail of the en- and L is the Planck length, while the ellipses include
P
tanglement spectrum. At least for gapped systems, it a log2((cid:96)/δ) term for odd d > 1. The Planck length
hasbeenarguedthatthecloserareduceddensitymatrix arises from Newton’s constant, which appears in Ryu
eigenstate’s eigenvalue is to the tail of the entanglement and Takayangi’s famous formula relating the entangle-
spectrum,thefurtheritprobesthesystemawayfromthe ment entropy to an area of a minimal surface in the
cut.55 Thus, although it seems likely that the factor of AdSspace.19Thefactor(L/L )dscalesappropriatelyfor
P
3 relationship between the TFIM and Heisenberg bilayer counting the degrees of freedom in d-dimensional space
data remains for α(cid:46)1.3, it remains possible that devia- withthePlancklengthasashort-distancecutoff,andthe
tions, corrections, or even a phase transitionin α change cutoff-independent function q(θ) has been computed for
the relationship.16 d<6. In d=1, the coefficient indeed becomes precisely
Despite the uncertainties, Fig. 7 gives substantial sup- c/3,56 and in d > 2 other relations with conformal cen-
port to the hypothesis that a (π/2) for the Heisenberg tral charges can be found.31 Thus it is natural to believe
α
bilayer at its critical point is approximately thrice the that the coefficient of the log term in our case d = 2 is
value for the TFIM at its critical point. also proportional to some universal quantity like a cen-
tral charge giving a measure for the degrees of freedom.
Our results provide support for this idea.
V. SUMMARY AND DISCUSSION Finally, it is important to test this behavior on other
2+1 critical points with different types of effective low-
energy theories. From this work, the conjecture is that
We have studied the Renyi entanglement entropies
the coefficient of the corner-induced logarithm counts a
for a square-lattice Heisenberg bilayer antiferromagnet
very clear signal (i.e. the integer N in the O(N) theory),
withaNumericalLinkedClusterExpansionmethod,em-
and hence may be relatively insensitive to some finite-
ploying both Lanczos and DMRG techniques as cluster
size effects. These finite-size effects will continue to be
solvers. Focussing on the subleading logarithmic scal-
reducedinthefuturewiththefurtheradoptioninNLCE
ing contribution to the Renyi entropies that arises from
of the powerful and general DMRG method, which has
a 90 corner, we calculate the cutoff-independent coeffi-
◦
been shown to provide very accurate results on quasi-2d
cient a (π/2) of the logarithm directly at the quantum
α finite-sizesystems.40 Thus, thecalculationsinthispaper
criticalpointofthemodel. Thiscriticalpointisinthe3D
can immediately and straightforwardly be extended to
O(3) universality class, described by a three-component
φ4 field theory. We find that the universal coefficient is other lattice models.
thrice that computed previously28 for a quantum crit- AnobviousnextcandidateforstudyisanO(2)critical
point in 2+1d, such as occurs in a spin 1/2 XY model
ical point in the 3D Ising universality class, described
by a single-component φ4 field theory. Our calculation with bilayer couplings, or alternatively, in a symmetry-
breaking magnetic field. If these and other conventional
thusprovidessubstantialevidencethatthisuniversalco-
critical points confirm our scenario, then especially im-
efficient provides a measure of the number of degrees of
portant would be the study of critical points with exotic
freedom of the field theory. In this interacting but close-
structures to their low-energy effective theories. An ex-
to-free theory, this is the number of low-lying bosonic
citingprospectwouldbetheSU(2)-invariantHeisenberg
modes present at the quantum critical point.
model with four-site exchange (Sandvik’s J-Q model57),
There have been hints in past literature that the co-
which is believed to contain a critical point in the non-
efficient of the corner-induced logarithm contains valu-
compact CP1 universality class, a field theory that de-
able information about the effective low-energy critical
scribes two flavors of spinons interacting with a non-
theory. Most concretely, for z = 1 conformal quantum
(cid:54) compact U(1) gauge field.58,59
field theories, the coefficient is proportional to the cen-
tralchargeoftheconformalfieldtheoryusedtobuildthe
ground-state wavefunction.24 Since we found analogous
behavior in z = 1 (Lorentz-invariant) field theories, it Acknowledgments
is thus tantalizing to speculate that this coefficient pro-
vides a quantity analogous to a central charge for any We would like to acknowledge crucial discussions with
2+1-dimensional quantum critical point. A.Burkov,A.JohnBerlinsky,J.Carrasquilla,H.Casini,
Holographic calculations of the entanglement entropy A. Ferris, S. Inglis, R. Myers, W. Witczak-Krempa, and
10
A.Sandvik. Thesimulationswereperformedonthecom- Physics. Research at PI is supported by the Govern-
puting facilities of SHARCNET and on the Perimeter ment of Canada through Industry Canada and by the
Institute HPC. Support was provided by NSERC, the Province of Ontario through the Ministry of Economic
Canada Research Chair program, the Ontario Ministry Development&Innovation. PFissupportedbyNational
of Research and Innovation, the John Templeton Foun- ScienceFoundationgrantDMR/MPS1006549andRRPS
dation, and the Perimeter Institute (PI) for Theoretical is supported by NSF grant DMR-1004231.
1 Entanglement entropy in extended quantum systems, J. 29 S. Inglis and R. G. Melko, New Journal of Physics 15,
Phys. A (special issue), Vol. 424 (2009). 073048 (2013).
2 H. Casini and M. Huerta, Phys.Rev. D85, 125016 (2012), 30 T. Hirata and T. Takayanagi, JHEP 0702, 042 (2007),
arXiv:1202.5650 [hep-th]. arXiv:hep-th/0608213 [hep-th].
3 T. Grover, (2012), arXiv:1211.1392. 31 R. C. Myers and A. Singh, (2012), arXiv:1206.5225.
4 A. Zamolodchikov, JETP Lett. 43, 731 (1986). 32 L.Wang,K.S.D.Beach, andA.W.Sandvik,Phys.Rev.
5 C.Holzhey,F.Larsen, andF.Wilczek,Nucl.Phys.B424, B 73, 014431 (2006).
443 (1994), arXiv:hep-th/9403108 [hep-th]. 33 C. Holm and W. Janke, Phys. Rev. B 48, 936 (1993).
6 G. Vidal, J. I. Latorre, E. Rico, and A. Ki- 34 M.Rigol,T.Bryant, andR.R.P.Singh,Phys.Rev.Lett.
taev, Phys.Rev.Lett. 90, 227902 (2003), arXiv:quant- 97, 187202 (2006).
ph/0211074 [quant-ph]. 35 M. Rigol, T. Bryant, and R. R. P. Singh, Phys. Rev. E
7 V. E. Korepin, Phys. Rev. Lett. 92, 096402 (2004). 75, 061118 (2007).
8 P. Calabrese and J. Cardy, J. Stat. Mech.: Theor. Exp. 36 M. Rigol, T. Bryant, and R. R. P. Singh, Phys. Rev. E
P06002 (2004). 75, 061119 (2007).
9 L. Bombelli, R. K. Koul, J. Lee, and R. D. Sorkin, Phys. 37 B. Tang, E. Khatami, and M. Rigol, (2012),
Rev. D 34, 373 (1986). arXiv:1207.3366.
10 M. Srednicki, Phys. Rev. Lett. 71, 666 (1993). 38 S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
11 J.Eisert,M.Cramer, andM.B.Plenio,Rev.Mod.Phys. 39 U. Schollwo¨ck, Rev. Mod. Phys. 77, 259 (2005).
82, 277 (2010). 40 E. M. Stoudenmire and S. R. White, Annual Review of
12 M. M. Wolf, Phys. Rev. Lett. 96, 010404 (2006). Condensed Matter Physics 3, 111 (2012).
13 D.GioevandI.Klich,Phys.Rev.Lett.96,100503(2006). 41 A. Renyi, Proc. of the 4th Berkeley Symposium on Math-
14 H.-H.Lai,K.Yang, andN.E.Bonesteel,Phys.Rev.Lett. ematics, Statistics and Probability 1960, 547 (1961).
111, 210402 (2013). 42 M. A. Metlitski and T. Grover, (2011), arXiv:1112.5166.
15 D.V.Fursaev,Phys.Rev.D73,124025(2006),arXiv:hep- 43 Z. Weihong, Phys. Rev. B 55, 12267 (1997).
th/0602134 [hep-th]. 44 C. J. Hamer, J. Oitmaa, and Z. Weihong, Phys. Rev. B
16 M.A.Metlitski,C.A.Fuertes, andS.Sachdev,Phys.Rev. 85, 014432 (2012).
B 80, 115122 (2009). 45 M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi,
17 Forasimpleexplanationofhowcmeasuresdegreesoffree- and E. Vicari, Phys. Rev. B 65, 144520 (2002).
dom in 1+1 dimensions, see Ref. 59. 46 J. Zinn-Justin, Quantum field theory and critical phenom-
18 E. Fradkin, Field Theories of Condensed Matter Systems, ena (Clarendon Press, Oxford, 2002).
2nd ed. (Cambridge University Press, 2013). 47 A. V. Chubukov and D. K. Morr, Phys. Rev. B 52, 3521
19 S. Ryu and T. Takayanagi, Phys. Rev. Lett. 96, 181602 (1995).
(2006). 48 P. V. Shevchenko and O. P. Sushkov, Phys. Rev. B 59,
20 T. Nishioka, S. Ryu, and T. Takayanagi, J. Phys. A 42, 8383 (1999).
504008 (2009). 49 V. N. Kotov, O. Sushkov, Z. Weihong, and J. Oitmaa,
21 E. Ardonne, P. Fendley, and E. Fradkin, Annals Phys. Phys. Rev. Lett. 80, 5790 (1998).
310, 493 (2004), arXiv:cond-mat/0311466 [cond-mat]. 50 T.Sommer,M.Vojta, andK.W.Becker,Eur.Phys.J.B
22 D.S.RokhsarandS.A.Kivelson,Phys.Rev.Lett.61,2376 23, 329 (2001).
(1988). 51 G. Guennebaud, B. Jacob, et al., “Eigen v3,”
23 J. L. Cardy and I. Peschel, Nuclear Physics B 300, 377 http://eigen.tuxfamily.org (2010).
(1988). 52 U. Schollwo¨ck, Annals of Physics 326, 96 (2011).
24 E. Fradkin and J. E. Moore, Phys. Rev. Lett. 97, 050404 53 A. Sen and A. W. Sandvik, unpublished (2014).
(2006). 54 S. Humeniuk and T. Roscilde, Phys. Rev. B 86, 235116
25 H. Casini and M. Huerta, Nucl. Phys. B 764, 183 (2007). (2012).
26 A. B. Kallin, M. B. Hastings, R. G. Melko, and R. R. P. 55 A. M. Turner, F. Pollmann, and E. Berg, Phys. Rev. B
Singh, Phys. Rev. B 84, 165134 (2011). 83, 075102 (2011).
27 R. R. P. Singh, R. G. Melko, and J. Oitmaa, Phys. Rev. 56 J.D.BrownandM.Henneaux,Commun.Math.Phys.104,
B 86, 075106 (2012). 207 (1986).
28 A. B. Kallin, K. Hyatt, R. R. P. Singh, and R. G. Melko, 57 A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).
Phys. Rev. Lett. 110, 135702 (2013). 58 T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and
M. P. A. Fisher, Science 303, 1490 (2004).
59 R.K.Kaul,R.G.Melko, andA.W.Sandvik,Annu.Rev.
Con. Mat. Phys. 4, 179 (2013).