Deformations of the boundary theory of the square lattice AKLT model
John Martyn,
1, 2
Kohtaro Kato,
2
and Angelo Lucia
2, 3
1
Department of Physics, University of Maryland, College Park, MD 20742, USA
2
Institute for Quantum Information and Matter,
California Institute of Technology, Pasadena, CA 91125, USA
3
Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA 91125, USA
The 1D AKLT model is a paradigm of antiferromagnetism, and its ground state exhibits
symmetry-protected topological order. On a 2D lattice, the AKLT model has recently gained at-
tention because it too displays symmetry-protected topological order, and its ground state can act
as a resource state for measurement-based quantum computation. While the 1D model has been
shown to be gapped, it remains an open problem to prove the existence of a spectral gap on the 2D
square lattice, which would guarantee the robustness of the resource state. Recently, it has been
shown that one can deduce this spectral gap by analyzing the model’s boundary theory via a tensor
network representation of the ground state. In this work, we express the boundary state of the
2D AKLT model in terms of a classical loop model, where loops, vertices, and crossings are each
given a weight. We use numerical techniques to sample configurations of loops and subsequently
evaluate the boundary state and boundary Hamiltonian on a square lattice. As a result, we evidence
a spectral gap in the square lattice AKLT model. In addition, by varying the weights of the loops,
vertices, and crossings, we indicate the presence of three distinct phases exhibited by the classical
loop model.
I. INTRODUCTION
The boundary state of a quantum system encodes am-
ple information about the system’s bulk. This concept
has motivated numerous research developments in differ-
ent areas, from the theory of topological insulators
1–3
to
the AdS/CFT correspondence.
A correspondence between a system’s boundary and its
bulk is explicitly manifest in Projected Entangled Pair
States (PEPS)
4,5
. This type of tensor network state has
risen to prominence in the past decade due to its versa-
tility in approximating ground states and expressiveness
in capturing essential physical properties of exotic phases
of quantum matter.
The correspondence we discuss involves the spectral
gap of a quantum many-body Hamiltonian, or the dif-
ference between its two lowest energy levels. Quantum
phase transitions correspond to points in the phase dia-
gram where the spectral gap vanishes. While the prob-
lem of determining the presence of a spectral gap or the
related phase diagram is in general undecidable
6–8
, it is
tractable for some specific cases of interest.
In the case of PEPS models, solid numerical evidence
9
and analytical results
10,11
provide a connection between
the spectral gap of the Hamiltonian and locality prop-
erties of the boundary state. This bulk-boundary cor-
respondence can be used to evidence a spectral gap in
systems that are otherwise difficult to study analytically.
One noteworthy system is the 2D generalization of the
model proposed by Affleck, Kennedy, Lieb and Tasaki
(AKLT)
12,13
, whose ground state is known as a valence
bond state (VBS), or equivalently, a PEPS. Since the 2D
AKLT VBS can be used as a resource for measurement-
based quantum computation
14–16
, the spectral gap has
implications on the complexity of preparing the states
and their stability against noise. However, despite exten-
sive analyses of the 2D AKLT model, its spectral gap re-
mains to be fully understood. Various numerical studies
have suggested that the general AKLT model is gapped,
some of which have explicitly calculated the spectral gap
for small systems and estimated this gap in the thermo-
dynamic limit
17,18
. Very recently, the spectral gap for the
2D hexagonal model has been confirmed by numerically
verifying a rigorous finite size criteria
19
.
In this work, we numerically study the spectral gap of
the 2D AKLT model on a square lattice via the bulk-
boundary correspondence in the PEPS framework. For
the AKLT model, the boundary states are mapped to a
classical loop model, where vertices and loops are each
assigned a specific statistical weight. The corresponding
classical loop model on the hexagonal lattice has been
studied in depth, but not within the context of the AKLT
model
20
. However, unlike loops on a hexagonal lattice,
loops can intersect on a square lattice, and therefore the
analysis is more involved in our scenario.
We use random sampling to numerically construct the
boundary state and the boundary Hamiltonian on lat-
tices of different sizes. We find evidence of locality of the
boundary Hamiltonian and therefore indicate the exis-
tence of a spectral gap, in accordance to previous findings
in the literature
9,21,22
. We obtain further evidence that
the boundary of the AKLT model does indeed show the
locality features of a gapped model by varying the pa-
rameters of the classical loop model. By doing so, we ob-
tain a phase diagram that exhibits multiple phases. The
point corresponding to the AKLT model sits well inside
a non-critical phase, thus confirming that the model is
gapped. While many deformations of the AKLT model
have been proposed in the literature
9,23,24
, to the best
of our knowledge these deformations do not give rise to
boundary states related to these loop models, which ap-
pear instead for some recently proposed AKLT models
arXiv:1912.10327v2 [quant-ph] 7 Jul 2020
2
on decorated lattices
25,26
.
II. BACKGROUND
A. Projected Entangled Pair States (PEPS)
Adopting the conventions of Ref.
10
, PEPS are defined
as follows. First begin with a subset Λ of an infinite
graph on which to model a state. Associate with each
vertex
v
∈
Λ a
d
-dimensional Hilbert space
H
d
, known
as the physical Hilbert space. The total Hilbert space
of the graph is then
H
Λ
=
⊗
v
∈
Λ
H
d
. Next, denote the
set of edges of Λ by
E
Λ
. Associate with each end point
of an edge a Hilbert space
H
D
of dimension
D
, spanned
by an orthonormal basis
{|
j
〉}
D
j
=1
.
D
is known as the
bond dimension of the PEPS. On the Hilbert space of
edge
e
, denote the maximally entangled state by
|
ω
〉
e
=
1
√
D
∑
D
j
=1
|
j,j
〉
e
. Focus now on a vertex
v
of degree
r
, and
associate with this vertex a linear map
T
v
:
H
⊗
r
D
→ H
d
.
T
v
maps from states of the virtual Hilbert space,
H
⊗
r
D
,
to the physical Hilbert space,
H
d
.
With these definitions, a PEPS state on Λ (with no
outgoing edges) is defined, up to normalization, as
|
PEPS
Λ
〉
=
⊗
v
∈
Λ
T
v
⊗
e
∈
E
Λ
|
ω
〉
e
.
(1)
We interpret this as a state where the transformation
T
v
is applied to entangled pairs at the edges of each vertex,
producing a state of the physical Hilbert space. If the
bond dimension
D
grows as a polynomial in the system
size, then this PEPS description is called efficient. If Λ
has outgoing edges, we collect the non-contracted virtual
degrees of freedom into a boundary Hilbert space
H
∂
Λ
.
B. The Bulk-Boundary Correspondence
As explained in Refs.
9,21,27
, PEPS exhibit a bulk-
boundary correspondence, wherein the virtual degrees of
freedom at the boundary correspond to the physical de-
grees of freedom in the bulk. Explicitly, the boundary
theory is governed by the boundary Hamiltonian, defined
as follows. First, focus on a region Λ with all the outgoing
edges. The density matrix on the boundary of Λ is then
obtained by tracing out the bulk physical degrees of free-
dom and keeping only the remaining boundary virtual
degrees of freedom:
ρ
∂
Λ
= tr
Λ
(
|
PEPS
Λ
〉〈
PEPS
Λ
|
)
∈B
(
H
∂
Λ
)
,
(2)
where
B
(
H
∂
Λ
) is the set of bounded operators on
H
∂
Λ
.
The boundary Hamiltonian,
H
∂
Λ
, is defined by equating
ρ
∂
Λ
to a thermal state at inverse temperature
β
= 1:
ρ
∂
Λ
=
e
−
H
∂
Λ
.
(3)
Via the bulk-boundary correspondence, many proper-
ties of the bulk can be learned from
ρ
∂
Λ
. For instance,
the boundary Hamiltonian exhibits the same symmetries
as the bulk Hamiltonian. More importantly, the local-
ity of the boundary Hamiltonian indicates a spectral gap
in the bulk Hamiltonian
9
. This correspondence is also
useful for numerical analyses. Ref.
9
employs the bulk-
boundary correspondence to numerically study the en-
tanglement spectrum of the boundary theory of a modi-
fied 2D AKLT model.
C. The AKLT Model and the VBS
Following Ref.
22
, let us define the AKLT model on a
graph
G
= (
V,E
). With each vertex
k
∈
V
, we asso-
ciate a spin operator
~
S
k
whose spin value is
S
k
=
|
~
S
k
|
=
1
2
deg(
k
), where deg(
k
) is the degree of vertex
k
. The
AKLT model is then defined by the following Hamilto-
nian:
H
AKLT
=
∑
〈
k,l
〉∈
E
A
(
k,l
)
π
S
k
+
S
l
(
k,l
)
,
(4)
where
A
(
k,l
)
∈
R
+
are (arbitrary) coefficients, and
π
S
k
+
S
l
(
k,l
) is the projector onto the space of maximal
total spin of vertices
k
and
l
:
S
k,l
:=
S
k
+
S
l
. Explicit
formulae for
π
S
k
+
S
l
(
k,l
) are provided in Ref.
22
.
The ground state of this model is known as the valence-
bond-state (VBS)
21,22
. It is the unique ground state, up
to degeneracies arising from the states of the boundary
vertices. As shown in Refs.
4,25,28
, the VBS in 2 dimen-
sions has an exact PEPS representation of bond dimen-
sion
D
= 2. This representation has been used to study
the properties of the AKLT model on different lattices,
such as cylinders and hexagonal lattices
9
.
On a graph
G
= (
V,E
), the VBS is obtained by placing
at each edge (
k,k
′
) in
E
a singlet state shared between
k
and
k
′
(so that the virtual vector space at a vertex
k
is (
C
2
)
⊗
deg(
k
)
), and then projecting the virtual space
at each vertex onto its appropriate symmetric subspace,
which has dimension deg(
k
). Formally, the VBS is de-
fined as
|
VBS
〉
=
⊗
k
∈
V
T
2
,
deg(
k
)
sym
⊗
e
∈
E
|
s
〉
e
,
(5)
where
T
2
,
deg(
k
)
sym
: (
C
2
)
⊗
deg(
k
)
→
C
deg(
k
)
is the projector
onto the symmetric subspace, and
|
s
〉
=
1
√
2
(
|
01
〉−|
10
〉
)
is the singlet state. Note that we will ignore the normal-
ization factor in the rest of this paper.
D. The VBS as a PEPS
Eq. (5) is almost in the same form as Eq. (1), ex-
cept that we are contracting against antisymmetric sin-
glet states
|
s
〉
instead of
|
ω
〉
. In order to connect with the
3
existing results in the PEPS literature, one would like to
express the VBS as in Eq. (1). This is easily remedied by
representing each
|
s
〉
in terms of
|
ω
〉
, by applying a mul-
tiple of a Pauli matrix:
|
s
〉
= (
I
⊗
iσ
y
)
|
ω
〉
. This shows
that the VBS is indeed a PEPS.
For bipartite lattices, such as the square or hexagonal
lattice, one can then choose to distribute the
iσ
y
matri-
ces performing the rotation from
|
ω
〉
to
|
s
〉
in such a way
that either all of the virtual indices at a given site are
affected, or none are. With this choice, the tensor
T
v
representing the VBS state as a PEPS depends on which
bipartition the vertex
v
belongs to: in one case it will be
exactly
T
2
,
deg(
v
)
sym
, while in the other it will be given by
T
2
,
deg(
v
)
sym
(
iσ
y
)
⊗
deg(
v
)
. The latter is equal to
U
v
T
2
,
deg(
v
)
sym
for some local unitary
U
v
. When computing the bound-
ary state
ρ
∂
Λ
, these local unitaries
U
v
cancel out under
the partial trace.
On the contrary, at the virtual boundary sites, on
which we are not taking the trace when computing the
boundary state, the unitary matrices
iσ
y
will not cancel
out, so that the PEPS boundary state will be equal to
U
∂
Λ
tr
bulk
(
⊗
v
∈
V
T
2
,
deg(
v
)
sym
⊗
e
∈
E
|
s
〉〈
s
|
e
⊗
v
∈
V
T
2
,
deg(
v
)
†
sym
)
U
†
∂
Λ
,
where
U
∂
Λ
is a tensor product of Pauli
σ
y
acting on the
virtual boundary Hilbert space. Since this transforma-
tion does not affect the locality properties of the bound-
ary state, we will ignore it and compute
ρ
∂
Λ
without it,
with the understanding that what we obtain is the correct
PEPS boundary state for the VBS up to a local change
of basis.
Finally, in order to simplify the notation, we will de-
note by
P
2
,
deg(
k
)
sym
the composition
T
2
,
deg(
k
)
†
sym
T
2
,
deg(
k
)
sym
,
thinking of it as an operator on (
C
2
)
⊗
deg(
k
)
. With this
final simplification, the boundary state we are going to
compute is
ρ
∂
Λ
= tr
bulk
(
⊗
v
∈
V
P
2
,
deg(
v
)
sym
⊗
e
∈
E
|
s
〉〈
s
|
e
)
.
(6)
III. MAPPING THE BOUNDARY VBS TO A
CLASSICAL LOOP MODEL
In this section, we establish the connection between the
boundary state of the AKLT model and a classical loop
model. In order to make the description more accessible,
we apply our method to the VBS on lattices of increasing
complexity: a 1D chain, a 2D hexagonal lattice, and fi-
nally a 2D square lattice. The calculations regarding the
projection operators
P
2
,
deg(
k
)
sym
are presented in Appendix
A.
A. Boundary State of VBS on a 1D Chain
We begin by looking at the boundary state of the VBS
on a 1D chain. We can find the exact boundary state in
this case, and it will serve as a precursor to studying the
VBS on 2D graphs.
On the 1D chain, deg(
k
) = 2
∀
k
∈
V
. Consider an
interval [1
,n
], whose bulk is composed by
n
spin-1 degrees
of freedom, while the boundary consists of the two spin
1/2 degrees of freedom at the outgoing edges. We trace
out the bulk and we denote by
ρ
n
the resulting boundary
state, that is,
ρ
n
= tr
bulk
(
n
⊗
i
=1
P
2
,
2
sym
|
s
〉〈
s
|
⊗
n
+1
)
.
(7)
We can write
P
2
,
2
sym
as (see Appendix A 1 for details):
P
2
,
2
sym
=
3
4
(
I
+
1
3
σ
1
·
σ
2
)
,
(8)
where
σ
1
·
σ
2
:=
σ
x
⊗
σ
x
+
σ
y
⊗
σ
y
+
σ
z
⊗
σ
z
, and the
sub-indices on
σ
denote the virtual sites on which the
operator acts. One can then expand the tensor product
⊗
n
i
=1
P
2
,
2
sym
as a polynomial in the
σ
·
σ
terms:
n
⊗
i
=1
P
2
,
2
sym
=
(
3
4
)
n
∑
J
⊂
[1
,n
]
1
3
|
J
|
⊗
i
∈
J
σ
i,
1
·
σ
i,
2
,
(9)
where the sub-index denotes which of the two virtual sites
of each site the operator acts on. However, it turns out
that when contracted against the
|
s
〉
states in Eq. (7), all
except two of these terms vanish.
To see this relation and straightforwardly calculate the
boundary state, it is convenient to introduce a diagram-
matic representation of the expansion in Eq. (9). For
each term in the expansion in Eq. (9), we say a site is
“turned on” if
σ
·
σ
acts on it, and “turned off” if the
identity
I
acts on it. We then associate a configuration
of
strings
to each term by drawing lines on the edges
adjacent to every turned on site. A string ends when it
reaches a turned off site or one of the boundary sites.
The expansion in Eq. (9) is then an expansion over
all the possible configuration of such strings (where each
string corresponds to a connected component of the set
J
). It is easy to show that if a string begins or ter-
minates in the bulk then it necessarily vanishes, since
〈
s
|
I
⊗
σ
i
|
s
〉
= 0 for
i
=
x,y,z
.
Therefore, there are only two configurations that con-
tribute to
ρ
n
: the configuration with all sites turned
off, which has no strings, and the configuration with all
sites turned on, which has a string connecting the two
boundary qubits. Neither of these configurations contain
strings that begin or terminate in the bulk.
The configuration with no strings yields a contribu-
tion equal to the identity.
By using the fact that
〈
s
|
σ
i
⊗
σ
j
|
s
〉
= (
−
1)
δ
ij
, one can see that the configu-
ration with all sites turned on yields a contribution of
4
(
−
1)
n
+1
3
−
n
σ
0
·
σ
n
+1
(where 0 and
n
+ 1 denote the two
boundary spins). Therefore, after normalizing the state,
the boundary state of the 1D VBS in 1D is
ρ
n
=
I
4
+
(
−
1)
n
+1
3
n
σ
0
·
σ
n
+1
4
.
(10)
B. Boundary State of VBS on a Hexagonal Lattice
We now turn to the VBS on a 2D hexagonal lattice,
which we can study analogously to the 1D case. On this
lattice, deg(
k
) = 3, so we will use (see Appendix A 2 for
details)
P
2
,
3
sym
=
1
2
(
I
+
1
3
(
σ
1
·
σ
2
+
σ
2
·
σ
3
+
σ
1
·
σ
3
)
)
,
(11)
where again the subindex denotes the virtual leg that the
operator acts on. In this case, there are three types of
turned on sites corresponding to different combinations
of outgoing lines determined by the subindices of
σ
’s. We
construct a diagrammatic representation by drawing lines
only on the outgoing edges of each site where the
σ
’s are
acting. By using arguments similar to those in the 1D
case, we find that the only contributions to
ρ
∂
come from
configurations where no strings begin or terminate in the
bulk. Thus,
ρ
∂
is a sum over configurations where strings
either form closed loops in the bulk, or terminate at the
boundary of the hexagonal lattice. Since each vertex has
degree 3 and the number of
σ
’s in the expansion at each
site has to be even, no crossings of the loops or strings
are allowed.
As the only terms contributing to the boundary state
are configurations of loop and strings with endpoints on
the boundary, we can collect all the terms with the same
endpoints, so that
ρ
∂
will be a weighted sum of products
of terms
σ
i
·
σ
j
, where
i
and
j
are virtual sites at the
boundary. Let
i
=
{
(
i
1
,j
1
)
,
(
i
2
,j
2
)
,...,
(
i
k
,j
k
)
}
be a set
of disjoint pairs of boundary qubits and
I
=
{
i
}
be the
set of all possible such disjoint pairs. The boundary state
is then given by
ρ
∂
=
1
N
∑
i
∈I
Z
i
n
∏
k
=1
σ
i
k
·
σ
j
k
,
(12)
where
{
Z
i
}
i
are coefficients and
N
is a normalization fac-
tor. Each coefficient is obtained by summing the contri-
butions from all possible loops and string configurations
with endpoints
i
, which we denote by
C
i
:
Z
i
=
∑
∆
∈C
i
w
(∆)
.
(13)
Here,
w
(
·
) is the weight function determined by the con-
traction of the tensors, which can be computed with the
following rules, starting from the empty diagram with no
strings which is given an arbitrary weight (the choice of
which will only affect the normalization of the boundary
state, so we set it equal to 1 without loss of generality):
1. Each edge contributes a multiplicative factor of
(
−
1);
2. Each loop contributes a multiplicative factor of 3;
3. Each vertex with a string passing through it con-
tributes a multiplicative factor of 1
/
3.
For a given configuration ∆, if we denote by
m
the to-
tal number of edges covered by strings,
`
the number of
loops, and
v
2
the number of bulk vertices with a string
passing through them (we choose the index 2 because this
is the degree of the vertex), then we obtain the following
weight, which we can interpret as a Boltzmann factor up
to a sign:
w
(∆) = (
−
1)
m
3
−
v
2
+
`
= (
−
1)
m
exp[
−
ln(3)(
v
2
−
`
)]
.
(14)
On a bipartite graph like the hexagonal lattice, the
weights of configurations with the same endpoints all
have the same sign, so the sign of these weights only
affects the overall sign of
Z
i
without causing any cancel-
lation. Therefore, up to the sign, each coefficient
Z
i
can
be interpreted as a partition function of a classical loop
model at inverse temperature
β
= 1 with boundary con-
ditions fixed by
i
, and the energy of a configuration given
by ln(3)(
v
2
−
`
). This has been previously analyzed as a
classical loop model on a hexagonal lattice, albeit not in
connection to the AKLT model
20
.
It is also now straightforward to calculate the normal-
ization factor
N
in Eq. (12). Because tr(
σ
) = 0,
N
is
determined entirely by the coefficient of the identity in
ρ
∂
, which corresponds to configurations with no strings
connecting boundary qubits. Denoting this coefficient by
Z
∅
(we use the empty set to signify the identity), and let-
ting
N
∂
be the number of boundary qubits, we see that
N
= 2
N
∂
Z
∅
.
C. Boundary State of VBS on a 2D Square Lattice
Our goal is to study the VBS on a 2D square graph.
We have now built up enough formalism to approach this
problem.
On a 2D square graph, deg(
k
) = 4, so the relevant
projector is now (see Appendix A 3 for details)
P
2
,
4
sym
=
15
48
(
I
+
1
3
∑
(
i,j
)
∈S
4
(
σ
i
·
σ
j
)
+
1
15
∑
(
i,j
)(
k,l
)
∈S
4
(
σ
i
·
σ
j
)(
σ
k
·
σ
l
)
)
.
(15)
Here,
S
4
is the set of permutations of 4 elements,
∑
(
ij
)
∈S
4
denotes the sum over all 2-cycles in
S
4
(there
are 6 of these), and
∑
(
ij
)(
kl
)
∈S
4
denotes the sum over all
disjoint (2
,
2)-cycles in
S
4
(there are 3 of these).
Evidently, this projector contains both 2-body and 4-
body interactions. Similar to the cases on the 1D chain
5
and the hexagonal lattice, we can again find
ρ
∂
and ex-
pand it in terms of
σ
i
k
·
σ
j
k
terms as in Eq. (12). The
coefficients of these terms are again sums over allowed
configurations of loops and strings. However, due to the
presence of the 4-body term in
P
2
,
4
sym
, the loops and strings
are now allowed to cross. In addition, there is now an in-
teraction term between loops and strings passing through
the same vertex, irrespectively of whether they cross or
not: they are weighted by 1
/
15 instead of the weight of
(1
/
3)
2
= 1
/
9 that they would receive as the product of
independent strings.
Analogous to the rules on the hexagonal lattice, the
rules for calculating the weight
w
(∆) of a configuration
on the square lattice are given by:
1. Each edge contributes a multiplicative factor of
(
−
1);
2. Each loop contributes a multiplicative factor of 3;
3. Each vertex with a single string passing through it
contributes a multiplicative factor of 1
/
3;
4. Each vertex with two single string passing through
it contributes a multiplicative factor of 1
/
15.
If we denote by
m
the total number of edges covered by
strings,
`
the number of closed loops,
v
2
the number of
bulk vertices with only one string passing through, and
v
4
the number of vertices with two strings passing through
(again we choose indices 2 and 4 because these are the
degrees of the respective vertices), then the weight of an
arbitrary configuration is given by the Boltzmann weight:
w
(∆) = (
−
1)
m
3
−
v
2
+
`
(15)
−
v
4
=
(
−
1)
m
exp[
−
ln(3)(
v
2
−
`
)
−
ln(15)
v
4
]
.
(16)
As before, because of the bipartite nature of the lattice,
the sign (
−
1)
m
only depend on the boundary sites
i
, so
that we can interpret the coefficients in
ρ
∂
as partition
functions of a classical loop model. As an example, in
Fig. 1 we display a sample configuration and calculate
its weight.
D. Locality at the Boundary
Our goal is to determine the locality of the boundary
Hamiltonian,
H
∂
, on the 2D cylindrical square lattice,
which will evidence the presence of a spectral gap. As
we describe in Appendix B, we study cylindrical lattices
with
N
x
plaquettes in the periodic direction and
N
y
pla-
quettes in the longitudinal direction. We use the nu-
merical procedure outline in Appendix B to estimate
ρ
∂
and subsequently compute the boundary Hamiltonian as
H
∂
=
−
log(
ρ
∂
).
FIG. 1: An example configuration on a square lattice.
The circled vertex with degree 2 (one string passing
through it) contributes a multiplicative factor of
1
3
,
whereas the circled vertex with degree 4 (two strings
passing through it) contributes a factor of
1
15
. We can
compute its weight straightforwardly: we have
m
= 9,
`
= 1,
v
2
= 6, and
v
4
= 1. Hence the weight of this
configuration is
−
(3
5
·
15)
−
1
. Because the strings of this
configuration connect boundary spins 2 and 3, this
weight contributes to
Z
{
(2
,
3)
}
.
1. Probing Locality
To probe the locality of the boundary Hamiltonian, we
look at the amplitude of the Heisenberg interactions:
A
r
=
1
3
·
2
N
x
Tr
(
H
∂
σ
i
·
σ
i
+
r
)
(17)
(this is the same for all
i
due to the rotational symmetry
of the cylinder). This amplitude measures the strength
of interactions between two qubits separated by
r
lattice
spacings. It is the coefficient of the Heisenberg term in
the boundary Hamiltonian:
H
∂
=
∑
r
A
r
∑
i
σ
i
·
σ
i
+
r
+
...
(18)
In addition, we will study the strength of interactions on
subsets of n qubits:
d
n
= Tr
(
h
2
n
)
,
(19)
where
h
n
contains the terms in
H
∂
with interaction range
n
. As defined in Ref.
9
, these are the terms in which
the largest contiguous block of identity operators acts on
6
N
x
−
n
qubits. For instance,
h
0
consists of the term in
H
∂
that is proportional to the identity,
h
2
consists of
all Heisenberg interactions on neighboring qubits, and so
on. Because
H
∂
is comprised of products of Heisenberg
interactions, which can be deduced from Eq. (12),
d
1
=
0. All other
d
n
are nonzero.
Similar quantities are analyzed in Ref.
9
in order to
study the locality of a boundary Hamiltonian. For a sys-
tem of local interactions, we expect both
A
r
and
d
n
to
decay exponentially in their respective argument or van-
ish past a certain argument value. We will use these
features to diagnose locality.
2. Results
In Fig. 2, we plot
A
r
and
d
n
on a lattice of size
N
x
=
N
y
= 10. We see that
A
r
decays exponentially in
r
, and
d
n
decays exponentially in
n
. Similar results were seen
for lattices of other sizes. Clearly, the interactions in the
Hamiltonian decay exponentially in their strength and
size, and thus the boundary Hamiltonian is quasi-local.
To justify this quasi-locality in the thermodynamic
limit, we performed a finite size scaling analysis to ex-
trapolate
A
r
in the limit of an infinitely long cylin-
der (
N
y
→ ∞
). The results are shown below in Fig.
3. Note that we collect data on systems up to size
(
N
x
, N
y
) = (8
,
20).
The finite size scaling analysis indicates that
A
r
(for
r
≥
1) continues to decay exponentially in the limit that
N
y
→ ∞
. Hence, the boundary Hamiltonian is quasi-
local in the thermodynamic limit. In Ref.
10
it was shown
that a strictly local boundary Hamiltonian implies a spec-
tral gap in the bulk, and it was furthermore conjecture
that the same results holds for quasi-local interactions,
which has been partially confirmed by the recent results
of Ref.
11
. The locality estimate we have obtained thus
provides a strong evidence that the 2D AKLT model on
a square lattice is gapped.
IV. PHASE STRUCTURE OF THE LOOP
MODEL
A. Procedure
Here we investigate the phase structure of the loop
model by varying the factors associated with closed loops,
vertices with a single string passing through, and vertices
with two strings. These factors were 3,
1
3
, and
1
15
, respec-
tively, in our previous analysis. By varying these param-
eters we can induce qualitative changes in the boundary
state and reveal a phase diagram for the classical loop
model.
To study the phase structure, we introduce a parame-
FIG. 2: Plots of the Heisenberg interaction amplitudes
(
A
r
), the n-qubit interaction strengths (
d
n
), and the
corresponding lines of best fit. This data was collected
on a lattice of size
N
x
=
N
y
= 10. Also note that these
plots only extend to
r
= 5 and
n
= 6 because the
rotational symmetry of the cylinder inhibits exponential
decay at larger argument values. For instance,
A
1
=
A
9
because of this symmetry.
terized operator
P
(
c
2
,c
4
) =
I
+
c
2
∑
(
i,j
)
∈S
4
(
σ
i
·
σ
j
)
+
c
4
∑
(
i,j
)(
k,l
)
∈S
4
(
σ
i
·
σ
j
)(
σ
k
·
σ
l
)
,
which when contracted against the network of singlet
states
|
s
〉
produces a boundary state described by a loop
model with closed loops weighted by 3, vertices with a
single string passing through by
c
2
, and vertices with
two strings by
c
4
. By replacing the operator
σ
·
σ
by a
7
FIG. 3: Finite size scaling analysis of Heisenberg
interaction amplitudes (
A
r
) as a function of 1
/N
y
on a
cylindrical lattice with
N
x
= 8. Data is displayed for
N
y
= 2
,
4
,
8
,
12
,
20. We note that, although
A
0
,
A
1
,
A
2
, and
A
3
appear constant, they actually grow as a
function of 1
/N
y
, but less noticeably than
A
4
.
sum over only one or two of the Pauli matrices, we can
change the loop weight to 1 or 2, respectively. Higher
integer loops weight can be obtained similarly by allow-
ing a larger bond dimension and choosing a larger set of
orthonormal traceless Hermitian unitary matrices. We
denote this loop weight by
c
`
.
In general,
P
(
c
2
,c
4
) is not a projector. Except for a
few special points of the parameters (such as the one
corresponding to the VBS), the rank of
P
(
c
2
,c
4
) will be
maximal, which implies a growth in the local physical
dimension.
Furthermore,
P
(
c
2
,c
4
) is not necessarily positive semi-
definite either. However, the boundary state computed
with
P
(
c
2
,c
4
) can be obtained from a PEPS contraction
only at points where
P
(
c
2
,c
4
) is positive semi-definite.
In this regime, one can interpret a change in
c
2
and/or
c
4
as a perturbation of the PEPS tensor of the VBS after
renormalizing at the injectivity length. Note that this
perturbation does not break the
SU
(2) symmetry of the
tensor. Outside of the range where
P
(
c
2
,c
4
) is positive
semi-definite, the resulting boundary state is not guar-
anteed to be a valid density matrix. For instance, it may
have negative eigenvalues and therefore not be a physical
quantum state. However, it is still meaningful to study
the classical loop model in this region, as it exhibits an
interesting phase structure.
To probe the behavior of the boundary state when the
weights are changed, we study the spin-spin correlation
functions
T
k
=
〈
σ
i
·
σ
i
+
k
〉−〈
σ
i
〉·〈
σ
i
+
k
〉
, where the average
is computed from the boundary state. (This definition
is independent of
i
by the rotational symmetry of the
cylinder.) At first glance, this quantity may seem ill-
defined when the boundary state is not physical, but
T
k
really only depends on the classical loop model. To see
this, note that Eq. (12) implies
〈
σ
i
〉
= 0 and that
〈
σ
i
·
σ
i
+
k
〉
=
Z
{
(
i,i
+
k
)
}
/Z
∅
, where we have used the value of
N
mentioned in Sec. III B. Therefore,
T
k
=
Z
{
(
i,i
+
k
)
}
/Z
∅
is simply the ratio of two classical partition functions,
which is well-defined even when the boundary state is
not physical.
We also record the average number of loops and the
average perimeter of loops, where these averages are
weighted means computed with the Boltzmann factors
of loop configurations. These loop properties, as well as
the spin-spin correlation functions, are calculated using
the random sampling procedure detailed in Appendix B
to sample configurations on the lattice.
B. Phases
1. Integer
c
`
We first study
T
k
at integer
c
`
and various choices of
c
2
and
c
4
. We find two characteristic phases exhibited
by our model: one in which
|
T
k
|
decays exponentially,
and one in which it decays sub-exponentially, roughly as
a power law. In addition, as indicated by our numerics,
the phases can be indexed by a simple, but non-local,
order parameter given by the ratio
η
:=
Average Number of Loops
Average Total Loop Perimeters
.
(20)
The average number of loops is computed as the weighted
average over all configurations of the number of loops in
a configuration. The average total loop perimeter is an
analogous average of the sum of perimeters of the loops
in a configuration. Qualitatively, this ratio measures the
size and prevalence of loops. Because both the numerator
and denominator are extensive quantities, this ratio can
be finite and nonzero in the thermodynamic limit.
We find that when the ground state is in the phase
of exponentially decaying correlations,
η
≈
1
/
4. On the
other hand, when the ground state is in the phase of
power law decaying correlations, 0
< η <
1
/
4; a typical
such value would be
η
∼
0
.
1. We also discover that
η
≈
0
corresponds to a nonphysical regime, as we discuss in
Sec. IV B 2.
In Fig. 4, we display plots of log
(
η
1
/
4
)
for various
c
2
and
c
4
at
c
`
= 3
,
2
,
1. We also enclose in dashed red
lines the regime in which the boundary state is physical.
We compute this regime by calculating the eigenvalues
of
P
(
c
2
,c
4
) to determine where this operator is positive
semi-definite.
8
FIG. 4: log
(
η
1
/
4
)
for various
c
2
and
c
4
, with integer
c
`
.
We also plot the contour
η
= 0
.
15 as a dashed gray line
to help visualize the phases (other values of
η
∼
0
.
15
produce similar lines). Likewise, the dashed red line en-
closes the region where the boundary state is physical, in
which
P
(
c
2
,c
4
)
≥
0.
(Top)
c
`
= 3. This plot contains
the AKLT model, which we indicate with a red circle.
(Middle)
c
`
= 2.
(Bottom)
c
`
= 1.
The region in Fig. 4 where log
(
η
1
/
4
)
≈
0 (the blue
region) encompasses the phase of exponentially decay-
ing correlations, and the region where log
(
η
1
/
4
)
<
0 (the
white region) encompasses the phase of power law decay-
ing correlations. From these plots, we see that the AKLT
model is well within the phase of exponentially decaying
correlations. Increasing
c
2
and/or
c
4
will cause the model
to enter the phase of power law decaying correlations.
Fig. 4 also indicates that increasing
c
4
sharpens the
transition between the two phases and shifts the transi-
tion to lower values of
c
2
. To study this transition more
closely, we calculated
η
as a function of
c
2
near the tran-
sition. We display our results for different lattice sizes in
Fig. 5. In the top plot, we fixed
c
`
= 3 and
c
4
= 0
.
2;
in the bottom plot, we fixed
c
`
= 3 and
c
4
= 1
.
2. It
is clear that
η
transitions from
≈
1
/
4 to a value around
∼
0
.
1 more sharply when
c
4
= 1
.
2 than when
c
4
= 0
.
2.
In addition, these plots indicate that, in the limit of an
infinite length cylinder,
η
converges to a value near 0
.
1
in the power law phase.
FIG. 5:
η
as a function of
c
2
near the transition. We have
fixed
N
x
= 10 and displayed data for various
N
y
.
(Top)
c
4
= 0
.
2 and
c
`
= 3.
(Bottom)
c
4
= 1
.
2 and
c
`
= 3.