Matrix product states for gauge eld theories

Boye Buyens,
1
Jutho Haegeman,
1
Karel Van Acoleyen,
1
Henri Verschelde,
1
and Frank Verstraete
1, 2
1
Department of Physics and Astronomy, Ghent University, Krijgslaan 281, S9, 9000 Gent, Belgium
2
Vienna Center for Quantum Science and Technology, Faculty of Physics,
University of Vienna, Boltzmanngasse 5, 1090 Vienna, Austria
The matrix product state formalism is used to simulate Hamiltonian lattice gauge theories. To
this end, we dene matrix product state manifolds which are manifestly gauge invariant. As an
application, we study 1+1 dimensional one avour quantum electrodynamics, also known as the
massive Schwinger model, and are able to determine very accurately the ground state properties
and elementary one-particle excitations in the continuum limit. In particular, a novel particle
excitation in the form of a heavy vector boson is uncovered, compatible with the strong coupling
expansion in the continuum. We also study non-equilibrium dynamics by simulating the real-time
evolution of the system induced by a quench in the form of a uniform background electric eld.
PACS numbers:
Gauge theories hold a most prominent place in physics.
They appear as eective low energy descriptions at dif-
ferent instances in condensed matter physics and nuclear
physics. But far and foremost they lie at the root of
our understanding of the four fundamental interactions
that are each mediated by the gauge elds correspond-
ing to a particular gauge symmetry. At the perturbative
quantum level, this picture translates to the Feynman di-
agrammatic approach that has produced physical predic-
tions with unlevelled precision, most famously in quan-
tum electrodynamics (QED). However the perturbative
approach miserably fails once the interactions become
strong. This problem is most pressing for quantum chro-
modynamics (QCD), where all low energy features like
quark connement, chiral symmetry breaking and mass
generation are essentially non-perturbative.
Lattice QCD, which is based on Monte Carlo sampling
of Wilsons Euclidean lattice version of gauge theories,
has historically been by far the most successfull method
in tackling this strongly coupled regime. Using up a siz-
able fraction of the global supercomputer time, state of
the art calculations have now reached impressive accu-
racy, for instance in the ab initio determination of the
light hadron masses [1]. But in spite of its clear superior-
ity, the lattice Monte Carlo sampling also suers from a
few drawbacks. There is the infamous sign problem that
prevents application to systems with large fermionic den-
sities. In addition, the use of Euclidean time as opposed
to real time, presents a serious barrier for the understand-
ing of dynamical non-equilibrium phenomena. Over the
last few years there has been a growing experimental and
theoretical interest in precisely these elusive regimes, e.g.
in the study of heavy ion collisions or early time cosmol-
ogy.
In this letter we study the application of tensor net-
work states (TNS) as a possible complementary approach
to the numerical simulation of gauge theories focusing
on the Schwinger model. This is highly relevant as this
Hamiltonian method is free from the sign problem and
allows for real-time dynamics. For the Schwinger model
the TNS approach has been studied before by Byrnes et
al [2] and Ba nuls et al [3]. By integrating out the gauge
eld (which one can only do for d=1+1), the model was
reduced to an ordinary spin model, yet with a non-local
Hamiltonian. Our approach is conceptually dierent, as
we keep the gauge eld degrees of freedom, which enables
us to take the thermodynamic limit, with the relevant
global symmetries exact. TNS have been considered also
for the discrete Z2 gauge theory, for d=1+1 by Sugihara
[4], and for d=2+1 by Tagliacozzo and Vidal [5].
Over the last decade the TNS framework has emerged
as a powerful tool for the study of local quantum many
body systems, exploring the fact that physical states (i.e.
ground states and their low energy excitations) only oc-
cupy a tiny corner of the full Hilbert space [6]. This is
exemplied by the relatively small amount of quantum
entanglement that these states possess. TNS are then
trial quantum states that precisely capture this feature,
allowing for relatively low cost numerical variational cal-
culations. In one spatial dimension, they also go by the
name of matrix product states (MPS), underlying the
well known density matrix renormalization group algo-
rithm (DMRG) [7]. At present MPS/DMRG is the state
of the art method in the numerical study of both static
and dynamical properties of d=1+1 strongly correlated
condensed matter systems. And also in higher dimen-
sions the TNS framework [8], although less developed, is
considered to be a promising candidate for the numerical
simulation of strongly interacting quantum many body
sytems.
The essential new ingredient with respect to the usual
MPS applications on quantum many body systems is that
for the Hamiltonian formulation of gauge theories, of the
full Hilbert space only the subspace of gauge invariant
states is actually physical. Although, due to Elitzurs
theorem [9] gauge invariance will not be broken on the
a
r
X
i
v
:
1
3
1
2
.
6
6
5
4
v
1


[
h
e
p
-
l
a
t
]


2
3

D
e
c

2
0
1
3
2
full Hilbert space, the low-energy excitations on the full
space will typically be completely dierent than those on
the constrained physical subspace. It is therefore crucial
to restrict the variational MPS manifold to this physical
subspace. Notice that the very same issue poses itself in
the context of the simulation of gauge theories with ultra-
cold atoms. See [10] for a recent proposal to implement
gauge invariance in that case.
The massive Schwinger model is QED in 1+1 dimen-
sions, with one avor of fermionic particles with mass m,
interacting through a U(1) gauge eld with coupling g
(which has mass dimension one for d=1+1). This model
shares some interesting features with QCD, most notably
the fermions are conned into zero charge bound state.
Furthermore in the continuum it can be studied by a
strong coupling expansion [11, 12], which makes it a per-
fect benchmark model. We will apply our gauge invariant
MPS construction on the Hamiltonian lattice formulation
of the model, focusing on the strongly coupled regime
g/m 1, and extrapolating our results to the contin-
uum. We determine the ground state and stable bound
states. In addition, we use our formalism to simulate
the full real time quantum dynamics induced by a back-
ground electric eld. This set up was recently also con-
sidered by Hebenstreit et al [13], with classical-statistical
simulations for the gauge elds.
The Schwinger Hamiltonian To write down a lattice
Hamiltonian for the Schwinger model, one starts from
the Lagrangian density in the continuum:
/ =

(

(i

+gA

) m)
1
4
F

. (1)
One can then perform a Hamiltonian quantization in the
time-like axial gauge (A
0
= 0), which can be turned
in a lattice system by the Kogut-Susskind spatial dis-
cretization [14] with the two-component fermions sited
on a staggered lattice. These fermionic degrees of free-
dom can then nally be converted to spin 1/2 degrees of
freedom by a Jordan-Wigner transformation, leading to
the gauged spin Hamiltonian (see [2] for more details):
H =
g
2

x
_

nZ
L(n)
2
+

2

nZ
(1)
n
(
z
(n) + (1)
n
)
+x

nZ
(
+
(n)e
i(n)

(n + 1) +h.c.)
_
. (2)
Here we have introduced the parameters x 1/(g
2
a
2
)
and 2

xm/g, with a the lattice spacing.

The spins live on the sites of the lattice, with

z
(n) [s
n
= s
n
[s
n
(s
n
= 1), and

= 1/2(
x
i
y
)
the spin ladder operators. Notice the dierent second
(mass) term in the Hamiltonian for even and odd sites.
This can be traced back to the staggered formulation,
with the even sites being reserved for the positrons
and the odd sites for the electrons. For the even
positron sites s
2n
= +1 can be viewed as an occu-
pied state, while s
2n
= 1 corresponds to an empty
state, and vice versa for the odd electron sites. The
gauge elds (n) = agA
1
(na/2), live on the links be-
tween the sites. Their conjugate momenta L(n), with
[(n), L(n

)] = i
n,n
, correspond to the electric eld,
gL(n) = E(na/2). Since (n) is an angular variable,
L(n) will have integer charge eigenvalues p
n
Z. The
local Hilbert space, spanned by the corresponding eigen-
kets [p
n
is therefore innite, but in practice we will do
a truncation and consider [p
n
[ p
max
in a numerical
scheme. For our calculations we take p
max
= 3.
The Hamiltonian (2) is invariant under T
2
, a trans-
lation over two sites, and the corresponding eigenvalues
read T
2
= e
2ika
, where k is the physical momentum of the
state. Another symmetry that will be useful is CT, ob-
tained by a translation over one site, followed by a charge
conjugation, C [s
n
, p
n
= [s
n
, p
n
. Since C
2
= 1, we
will have CT = e
ika
. The states with positive sign then
correspond to the scalar sector, while the negative sign
corresponds to the vector sector.
In addition, the Hamiltonian is invariant under the
residual time-independent local gauge transformations,
generated by
G
n
= L(n) L(n 1)
1
2
(
z
(n) + (1)
n
) . (3)
It is this gauge invariance that sets the Hamiltonian
quantization of gauge theories apart from the Hamilto-
nian quantization of ordinary systems. For gauge theo-
ries only the subspace of gauge-invariant states will be
physical: G
n
[
phys
= 0 for every n. This is called the
Gauss law constraint, as G
n
= 0 is indeed the discretized
version of
x
E = . And we will now show how one can
tailor the MPS formalism towards a constrained varia-
tional method on this physical gauge-invariant subspace.
Gauge invariant MPS. A general, not necessarily gauge-
invariant MPS for the lattice spin-gauge system (2) has
the form:

sn,pn
Tr[B
s1
1
C
p1
1
B
s2
2
C
p2
2
. . . C
p
2N
2N
] [s
1
, p
1
, s
2
, p
2
. . . , p
2N
,
(4)
where for now we consider a nite lattice of 2N sites.
Here, each B
sn
n
(and C
pn
n
) is a complex D D matrix
with components [B
sn
n
]

, that constitute the variational

parameters of the trial state. The indices , = 1, . . . D
are referred to as virtual indices, and D is called the bond
dimension.
Gauss law (see (3)) prescribes how to update the elec-
tric eld L(n) at the right link of a site n: either staying
with the value L(n 1) at the left in case there is no
charge at the site, or adding/subtracting one unit in case
there is a positive/negative charge at the site. This can
be conveyed by the matrix multiplications in an MPS
by giving the virtual indices a multiple index structure
3
(q,
q
), where q labels the charge, and taking the
matrices of the form:
[B
sn
n
]
(q,q),(r,r)
= [b
sn
n,q
]
q,r

q+(sn+(1)
n
)/2,r
[C
pn
n
]
(q,q),(r,r)
= [c
pn
n
]
q,r

q,pn

r,pn
. (5)
One can readily verify that an MPS (4) with matrices
of this form, indeed obeys the Gauss law constraint at
every site. Conversely, we show in the appendix that
every gauge-invariant state [, obeying G
n
[ = 0 for
every n, has an MPS representation of the form (5).
We now show how to apply the time-dependent varia-
tional principle (TDVP) of [20] to obtain a ground-state
approximation in the thermodynamic limit (N ),
taking into account the gauge invariance and anticipating
CT = 0 (see the appendix for the details). It will be use-
ful to block a site and link into one site with local Hilbert
spaced spanned by the states [q
2n1
= [s
2n1
, p
2n1

and [q
2n
= [s
2n
, p
2n
. The ground-state ansatz then
takes a form similar to a uniform MPS (uMPS) [20]:
[(A) =

qn
v

L
_

nZ
A
qn
_
v
R
[q
c
, (6)
where [q
c
= [(1)
n
q
n

nZ
, v
L
, v
R
C
D
, and A
q

C
DD
as follows from (5) of the general form (see ap-
pendix):
[A
{s,p}
]
(q,q);(r,r)
= [a
s,p
]
q,r

p,q+(s+1)/2

r,p
. (7)
The variational freedom of the gauge invariant state
[(A) then lies within the matrices a
s,p
C
DqDr
and the total bond dimension of the uMPS equals D =

qZ
D
q
.
The TDVP method evolves the Schrodinger equation
(SE), i
t
[
_
A
_
= H[
_
A
_
, within the manifold of
uMPS. To this end the right-hand side of the SE is re-
placed by [
_
B
H
(A), A
_
, where [(B, A) is given by

mZ

qn
v

L
_

n<m
A
qn
_
B
qm
_

n>m
A
qn
_
v
R
[q
c
, (8)
with B
q
also of the block structure form (7), B

H
(A) =
arg min
B

[(B, A) H[(A)

and [(B
H
(A), A)
[(A). The SE then boils down to an ordinary dieren-
tial equation for the variational parameters, i a = b
H
(a),
where b
H
(a) can be calculated in O
_
(2p
max
+1) max
p
D
3
p
_
time. Starting from a random state [(A) we will then
evolve towards the ground state by an imaginary time
evolution = it, that we stop once the state has con-
verged.
Once we have a good approximation for the ground
state, we can use the method of [21] to obtain the one-
particle excited states. The excitations are labelled by
their (physical) momentum k [/2a, /2a[ and their
0.4 0.2 0 0.2 0.4 0.6
3.9
3.95
4
4.05
4.1
4.15
k
3 2 1 0 1 2 3
18
16
14
12
10
8
6
4
2
0
p
FIG. 1: Results for m/g = 0.75. Left (a): distribution of the
logarithm of the Schmidt-coecients in every charge sector.
Right (b): Fit of the Einstein-dispersion relation E
2
= k
2
+
M
2
v,1
(x) for x = 100, 300, 800 (dashed lines) to the data (small
circles). The stars represent the estimated continuum values,
the full line (lowest lying curve) is the curve E
2
= k
2
+M
2
v,1
.
CT quantum number = 1. For a given ground-state
approximation we then take the following ansatz state
[
k,
(B, A) for the one-particle excitations:

mZ
e
ikma

qn
v

L
_

n<m
A
qn
_
B
qm
_

n>m
A
qn
_
v
R
[q
c
,
(9)
with B
q
again of the gauge-invariant form (7) with gen-
eral matrices b
s,p
. These are determined variationally
by minimizing their energy in the ansatz subspace which
leads to a generalized eigenvalue problem. For a given
momentum and CT quantum number we typically nd
dierent local minima of which only one or two are stable
under variation of the bond dimension D. It are these
stable states that we can interpret as approximations to
actual physical one-particle excitations.
Results. The continuum limit a 0 of the Schwinger
model corresponds to the limit x in (2). To obtain
the energies of the ground state and of the one-particle
excitations in this limit, we have calculated these quanti-
ties for values of x = 100, 200, 300, 400, 600, 800. At every
x we considered dierent values of D till convergence was
reached at some D
max
. We estimated the truncation er-
ror on D from comparison of the result for D = D
max
with the result for the next to largest value of D. Larger
values of x typically required larger values of D for the
same order of the error. For instance for m/g = 0.5
our maximal D varied from 185 for x = 100 to 358 for
x = 800. This scaling of D is not surprising, as it is
well known that MPS representations require larger D
for systems with larger correlation lengths (in units of
the lattice spacing) [6]. For the Schwinger model indeed
diverges in the x limit.
It was also important to choose the distribution of D
q
wisely, according to the relative weight of the dierent
charge sectors. As illustrated in g.1a, this is done by
looking at the Schmidt coecients for an arbitrary cut,
and demanding that the smallest coecients of each sec-
tor coincide more or less. The resulting distribution of
4
m/g 0 Mv,1 Ms,1 Mv,2
0 -0.318320(4) 0.56418(2)
0.125 -0.318319(4) 0.789491(8) 1.472(4) 2.10 (2)
0.25 -0.318316(3) 1.01917 (2) 1.7282(4) 2.339(3)
0.5 -0.318305(2) 1.487473(7) 2.2004 (1) 2.778 (2)
0.75 -0.318285(9) 1.96347(3) 2.658943(6) 3.2043(2)
1 -0.31826(2) 2.44441(1) 3.1182 (1) 3.640(4)
TABLE I: Energy density and masses of the one-particle ex-
citations (in units g = 1) for dierent m/g. The last column
displays the result for the heavy vector boson, compatible
with the prediction of Coleman [11, 12]
D
q
is peaked around q = 0, and justies our p
max
= 3
truncation that corresponds to D
q
= 0 for [q[ > 3.
To extrapolate towards x we used a third order
polynomial t in 1/

x through the largest ve x-values.

Similar to [2] our extrapolation error is then estimated by
considering a third and fourth order polynomial through
all six points, taking the error to be the maximal dier-
ence with the original inferred value.
In table 1 we display our resulting values for the ground
state energy density and the mass of the dierent one-
particle excitations. For m/g = 0 this can be com-
pared with the exact result that follows from bosoniza-
tion [15]. In this limit the model reduces to a free the-
ory, of one bosonic vector ( = 1) particle with mass
M
v,1
= 1/

= 0.56419 and with a ground-state energy

density
0
= 1/ = 0.318310 (both in units g = 1).
Furthermore, in the strong coupling expansion g/m
1 on this exact result, it is found that the vector boson be-
comes an interacting particle, leading to two more stable
bound states. There appears one scalar boson that is a
stable bound state of two vectors and one more vector bo-
son, that is best interpreted as a bound state of the scalar
and the original lowest mass vector [11, 12]. For g/m ,= 0
we also nd three excited states, one scalar and two vec-
tors, with the hierarchy of masses M
v,1
< M
s,1
< M
v,2
matching that of the strong coupling result. But notice
that for our values of g/m, the strong coupling expan-
sion result is not reliable anymore, making a quantita-
tive comparison useless. One can also show that in the
continuum limit the ground-state energy is independent
of g/m which is compatible with our ndings.
This is the rst time that the second vector excitation
has been found numerically. For the energy density and
the two lowest mass excitations our results are consistent
with the previous most precise simulations [2, 3], with a
similar or sometimes better accuracy.
A nice cross-check of our method also follows from cal-
culating the excitation energies for non-zero momenta k.
The Schwinger model is Lorentz invariant in the contin-
uum limit, so we should have an approximate Einstein
dispersion relation at nite lattice spacing a, for small
momenta ka 1. As shown in g.1b, this is precisely
what we nd.
0 5 10 15 20
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
El ectri c el d ( = 0. 75)
t
! " #! #" $!
!%&
!%'
!%$
!%#
!
!%#
Current densi ty ( = 0. 75)
(
FIG. 2: The electric eld and current density as function of
time for = 0.75 (in units g = 1). The results are plotted
for two dierent values of the bond dimension D = 239 (blue)
and D = 207 (green).
Finally, we have investigated the non-equilibrium dy-
namics induced by applying a uniform electric eld E
on the ground-state [
0
at time t = 0. Here we show
some rst results, a more detailed analysis will be pre-
sented elsewhere [16]. Physically, the situation corre-
sponds to the so called Schwinger particle creation mech-
anism [17], but now for a conning theory. This sub-
ject was recently also explored in the AdS/CFT set-up
[18]. In our set-up this can be simulated by applying
a uniform quench, replacing L(n) L(n) + in the
Hamiltonian (2), where E = g. Again we used TDVP,
but now for real-time evolution with the quench Hamil-
tonian. As the background eld breaks CT invariance
our ansatz is now dened simply by blocking two sites
and two links into one site and taking a gauge-invariant
form for A
q
that follows from (5). In g. 2 one can
see our results for the evolution of the electric eld,

0
(t)[L(2n 1) +L(2n)[
0
(t) and for the current den-
sity,

x
0
(t)[
+
(2n 1)e
i(2n1)

(2n) +h.c.[
0
(t);
and this for = 0.75. We observe the typical plasma os-
cillations that are damped over time, which is a typical
feature of thermalization. This is corroborated by the lin-
ear growth of the half-system entanglement entropy that
we nd in this case. It is precisely this growth of the
entropy that will at some point invalidate the MPS ap-
proximation for a given D. We can determine this point
self-consistently by looking at variations of the result for
dierent D. From g.2 we can infer then that at t 15,
our D = 239 MPS-result starts to become less reliable.
Conclusions. In this letter we have demonstrated the
potential of MPS as numerical method for gauge the-
ories. It is clear that we have only scratched the sur-
face of this approach and that even within the Schwinger
model there are many other types of calculations one
could do, like for instance the construction of two-particle
scattering states. Looking further aeld, one can eas-
ily generalize our gauge invariant MPS ansatz to other
5
gauge groups like SU(N) and also to higher dimensions.
Explicitly for d=2+1, the gauge-invariant 2d PEPS [8]
construction now involves ve-leg tensors with four vir-
tual indices and one physical index (c = charge) on
the sites, of the form [B
c
]
(q
l
,q
l
),(qr,qr
),(q
d
,q
d
),(qu,qu
)
=
[b
c
q
l
,qr,qu
]
q
l
,qr
,q
d
,qu

q
l
+q
d
+c,qr+qu
, while on the links
we get a three-leg tensor with two virtual indices
and one physical index (p = electric eld unit)
[C
p
]
(q
l
,q
l
),(q,qr
)
= [c
p
]
q
l
,qr

q
l
,p

qr,p
.
While preparing our manuscript the paper [19] ap-
peared with an approach that is conceptually close to
ours. There the authors use a quantum link model
to write down gauge invariant MPS for the Schwinger
model.
Acknowledgements. We thank Mari-Carmen Ba nuls for
suggesting us to look at real-time quench dynamics.
Furthermore we acknowledge very interesting discus-
sions with Mari-Carmen Ba nuls, David Dudal and Lucca
Tagliacozzo. This work is supported by an Odysseus
grant from the FWO, a PhD-grant from the FWO (B.B),
the FWF grants FoQuS and Vicom, and the ERC grant
QUERG.
[1] D urr, S., Fodor, Z., Frison, J., et al. 2008, Science, 322,
1224; Fodor, Z., & Hoelbling, C. 2012, Reviews of Mod-
ern Physics, 84, 449
[2] Byrnes, T. M. R., Sriganesh, P., Bursill, R. J., &
Hamer, C. J. 2002, Nuclear Physics B Proceedings
Supplements, 109, 202; Byrnes, T. M., Sriganesh, P.,
Bursill, R. J., & Hamer, C. J. 2002, Phys. Rev.
D , 66, 013002; Byrnes, T. M. R, google-books,
http://books.google.co.jp/books?id=FwaiXgJWM2MC
[3] Ba nuls, M. C., Cichy, K., Jansen, K., & Cirac, J. I. 2013,
JHEP 11 (2013) 158; Ba nuls, M. C., Cichy, K., Cirac,
J. I., Jansen, K., & Saito, H. 2013, arXiv:1310.4118
[4] Sugihara, T. 2005, Journal of High Energy Physics, 7, 22
[5] Tagliacozzo, L., & Vidal, G. 2011, Phys. Rev. B , 83,
115127
[6] Cirac, J. I., & Verstraete, F. 2009, Journal of Physics A
Mathematical General, 42, 4004
[7] White,S.R., Phys. Rev. Lett. 69, 2863-2866, 1992
[8] Verstraete, F., & Cirac, J. I. 2004, arXiv:cond-
mat/0407066
[9] Elitzur, S. 1975, Phys. Rev. D , 12, 3978
[10] Zohar, E., Cirac, J. I., & Reznik, B. 2013, Phys. Rev. A
, 88, 023617
[11] Coleman, S. 1976, Annals of Physics, 101, 239
[12] Adam, C. 1996, Physics Letters B, 382, 383
[13] Hebenstreit, F., Berges, J., & Gelfand, D. 2013, Phys.
Rev. D , 87, 105006
[14] Kogut, J., & Susskind, L. 1975, Phys. Rev. D , 11, 395
[15] Schwinger, J. 1962, Physical Review, 128, 2425
[16] Buyens, B., Haegeman, J., Verstraete F., Van Acoleyen
K., in preparation.
[17] Schwinger, J. 1951, Physical Review, 82, 664
[18] Kawai, D., Sato, Y., & Yoshida, K. 2013, arXiv:1312.4341
[19] Rico, E., Pichler, T., Dalmonte, M., Zoller, P., & Mon-
tangero, S. 2013, arXiv:1312.3127
[20] Haegeman, J., Cirac, I., Osborne, T.J., Pizorn, I., Ver-
schelde, H.,, Verstraete, F.,Phys. Rev. Lett. 107, 070601
(2011)
[21] Haegeman, J., Osborne, T.J., Verstraete, F., Phys. Rev.
B, vol. 88, Issue 7, id. 075133 (2013)
[22] Perez-Garcia, D., Verstraete F., Wolf, M.M., Cirac, J.I.,
Quantum Inf. Comput. 7 (2007), 401
[23] Haegeman, J, Pirvu, B., Weir, D.J., Cirac, I., Osborne,
T.J., , Verschelde, H.,, Verstraete, F., Phys. Rev. B 85,
100408(R) (2012)
6
A1: Gauge invariant states
Consider a lattice with 2N sites. The basis of the total Hilbert space 1 is [q [s
n
, p
n

1n2N
: s
n
= 1, p
n
Z.
A general state [ 1 can be written as a MPS in the canonical form (see [22], theorem 1):
[ =

{s}=1

{p}Z
N

n=1
B
s2n1
2n1
C
p2n1
2n1
B
s2n
2n
C
p2n
2n
[q , (10)
where B
s
k
k
C
D
k
D
k
, C
p
k
k
C
D
k
D
k+1
and D
1
= D
2N+1
= 1. By being in its canonical form we mean that

s=1
B
s
n
B
s
n

= 11,

pZ
C
p
n
C
p
n

= 11 (11)
and that there exist positive denite diagonal matrices l
C
n
and l
B
n
such that

s=1
B
s
n

l
B
n
B
s
n
= l
C
n
,

pZ
C
p
n

l
C
n
C
p
n
= l
B
n+1
. (12)
Because QED is a gauge theory we have to restrict to 1
phys
, the set of all gauge invariant states. This means that
every [ 1
phys
has to satisfy
[ = exp(i
l
G
l
) [ , l = 1, . . . , 2N. (13)
The right-hand side can also be written as a MPS with the same bond dimension:
exp(i
l
G
l
) [ =

{s}=1

{p}Z
tr
__
N

n=1

B
s2n1
2n1

C
p2n1
2n1

B
s2n
2n

C
p2n
2n
__
[q (14)
where if l > 1:

B
s
k
k
= B
s
k
k
for k ,= l,

C
s
k
k
= C
s
k
k
for k ,= l 1, l and

C
p
l1
= e
i
l
p
C
p
l1
,

B
s
l
= e
i
l
[
s+(1)
l
2
]
B
s
l
,

C
p
l
=
e
i
l
p
C
p
l
and if l = 1:

B
s
k
k
= B
s
k
k
,

C
s
k
k
= C
s
k
k
for k ,= 1 and

B
s
1
= e
i1[
s1
2
]
B
s
1
,

C
p
1
= e
i1p
C
p
1
. Because the two MPS
with identical bond dimension have to represent the same state [ and (10) is assumed to be in the canonical form,
it follows by [22], theorem 2, that there exists invertible square matrices U
n
and V
n
such that

B
sn
n
= U
1
n
B
sn
n
V
n
,

C
sn
n
= V
1
n
C
sn
n
U
n+1
where U
1
= 1 and U
2N+1
= 1. Note that U
n
and V
n
depend on
l
and l.
The matrices U
n
and V
n
are unitary matrices. Indeed, its not hard to check that

B
s
n
and

C
p
n
also obey (11) and
(12). For n = 2N we have that

C
s
2N
= V
2N
C
s
2N
which implies that V
2N
V

2N
= 11. Using this,

B
p
2N
= U
1
2N
B
p
2N
V
p
2N
and the fact that B
p
2N
and

B
p
2N
obey (11) it follows that U
2N
U

2N
= 11. Proceeding in the same way from n = 2N till
n = 1 one sees that all the matrices U
n
and V
n
are unitary.
Now we will prove that U
n
, V
n
= 11 for n ,= l. If n < l we may assume that l > 1. Note that U
1
= 1 and that
B
s
1
=

B
s
1
= B
s
1
V
1
. Using (12) it follows that V
1
= (l
C
1
)
1

s
B
s
1

l
B
1
B
s
1
= 11. Assume now V
n1
= 11 (n < l 1) then
C
p
n1
=

C
p
n1
= C
p
n1
U
n
, which implies U
n
= (l
B
n
)
1

p
C
p
n1

l
C
n1
C
p
n1
= 11, i.e. V
n1
= 11 implies that U
n
= 11
for n < l 1. Using the same ideas one proofs that U
n
= 11 implies V
n
= 11 (n < l). This concludes the case
n < l. For n > l one starts from

C
p
2N
= C
p
2N
= V
2N
C
p
2N
. From (11), we obtain that V
2N
= 11. As a consequence

B
s
2N1
= B
s
2N1
= U
2N
B
s
2N
holds. By (11) it follows that U
2N
= 11. One can now repeat this reasoning and see that
U
n
, V
n
= 11 for all n > l.
So the MPS (10) is gauge invariant i for every l = 1, . . . , 2N there exist unitary matrices U
l
and V
l
(depending on

l
) such that
U

l
B
s
l
V
l
= e
i
l
[
s+(1)
l
2
]
B
s
l
, C
p
l1
U
l
= e
i
l
p
C
p
l1
(l > 1), V

l
C
p
l
= e
i
l
p
C
p
l
. (15)
Consider now the case
l
= 1, then the matrices U
l
and V
l
do not depend on
l
anymore. The unitary matrices can
be diagonalized (as exponential of a Hermitian matrix): U
l
= W

l

U
l
W
l
, V
l
= X

l

V
l
X
l
, where W
l
, X
l
are unitary
matrices and
U
l
and
V
l
are diagonal matrices where all the diagonal-elements have modulus one. If we perform
the following MPS-gauge transformation:
B
s
l


B
s
l
W
l
B
s
l
X

l
, C
p
l


C
p
l
X
l
C
p
l
W

l+1
, W
1
= W
2N+1
= 1 (16)
7
the MPS (10) is unaected and the conditions (15) now read

U
l

B
s
l

V
l
= e
i(s+(1)
l
)/2

B
s
l
,

C
p
l1

U
l
= e
ip

C
p
l1
(l > 1),

V
l

C
p
l
= e
ip

C
p
l
. (17)
The property (11) will also hold for

B and

C, however the property (12) is modied in the sense that l
B
n
and l
C
n
are
not diagonal anymore (but they remain positive denite). We will denote this matrices with l

B
n
and l

C
n
. As already
mentioned, the entries of the diagonal matrices
U
l
and
V
l
are complex phase factors. Let e
i
l,j
, j = 1, . . . , n
U
l
,
respectively e
i
l,j
, j = 1, . . . , n
v
l
be the eigenvalues of
U
l
with multiplicity m(
l,j
) respectively of
V
l
with
multiplicity m(
l,j
),

U
l
=
nu
l

j=1
m(
l,j
)

j=1
e
i
l,j
[
l,j
,
j

l,j
,
j
[,
V
l
=
nv
l

j=1
m(
l,j
)

j=1
e
i
l,j
[
l,j
,
j

l,j
,
j
[, (18)
then we can write

B and

C as

B
s
l
=
nu
l

j=1
nv
l

k=1
m(
l,j
)

j=1
m(
l,k
)

k
=1
[

B
s
l
]
(
l,j
,j);(
l,k
,
k
)
[
l,j
,
j

l,k
,
k
[, l > 1 (19a)

C
p
l
=
nv
l

j=1
nu
l+1

k=1
m(
l,j
)

j=1
m(
l+1,k
)

k
=1
[

C
p
l
]
(
l,j
,j);(
l+1,k
,
k
)
[
l,j
,
j

l+1,k
,
k
[, l < 2N (19b)

B
s
1
=
nv
1

k=1
m(
1,k
)

k
=1
[

B
s
1
]
1;(
1,k
,
k
)

1,k
,
k
[,

C
p
2N
=
nv
2N

j=1
m(
2N,j
)

j=1
[

C
p
2N
]
(
2N,j
,j);1
[
2N,j
,
j
. (19c)
Using (17) it follows that
(e
i(p
l+1,k
)
1)[

C
p
l
]
(
l,j
,j);(
l+1,k
,
k
)
= 0, (e
i(p
l,j
)
1)[

C
p
l
]
(
l,j
,j);(
l+1,k
,
k
)
= 0. (20)
(e
i(p
2N,j
)
1)[

C
p
2N
]
(
2N,j
,j);1
= 0, (21)
so
[

C
p
l
]
(
l,j
,j);(
l+1,k
,
k
)
=
p,
l,j

p,
l+1,k
[c
p
l
]
j,
k
, [

C
p
2N
]
(
2N,j
,j);1
=
p,
2N,j
[c
p
2N
]
j,1
, (22)
Note that
l,j
and
l,j
are only unique up to a multiple of 2. By writing
p,
l,j
we mean that we must take for
l,j
up to a multiple of 2 the value p. Of course this will not inuence the eigenvalue e
i
l,j
.
Assume now that there would exist a
l+1,k0
(l < 2N) with
l+1,k0
,= p, p Z. Then it follows by (22) that
[

C
p
l
]
(
l,j
,j);(
l+1,k
0
,
k
0
)
= 0, (23)
p Z, j = 1, . . . , n
v
l
,
j
= 1, . . . , m(
l,j
). If we now consider the non-singular matrix l

C
l
, see (12), then
_

pZ
(

C
p
l
)

C
l

C
p
l
_
(
l+1,k
0
,
k
0
),(
l+1,k
,
k
)
= 0, (24)

k0
= 1, . . . , m(
l,k0
),k = 1, . . . , n
u
l+1
,
k
= 1, . . . , m(

l+1,k
). By (12) this would mean that l

B
l+1
has a zero-row
and would be singular which is a contradiction because l

B
l+1
is positive denite. As a consequence all the
l,k
are
integers. In the same way, but now by using the condition (11) one proves that all the
l,j
are integers.
We can write (18) as

u
l
=

qZ
D
l
q

q=1
e
iq
[q,
q
q,
q
[,
v
l
=

qZ
D

l
q

q=1
e
iq
[q,
q
q,
q
[, (25)
8
and expand

B,

C:

B
s
l
=

q,rZ
D
l
q

q=1
D

l
r

r=1
[

B
s
l
]
(q,q);(r,r)
[q,
q
r,
r
[,

C
p
l
=

q,rZ
D

l
q

q=1
D
l+1
r

r=1
[

C
p
l
]
(q,q);(r,r)
[q,
q
r,
r
[. (26)

B
s
1
=

rZ
D

l
r

r=1
[

B
s
1
]
1;(r,r)
r,
r
[,

C
p
2N
=

qZ
D

2N
q

q=1
[

C
p
l
]
(q,q);1
[q,
q
(27)
where D
l
q
respectively D

l
q
denotes the multiplicity of the eigenvalue q in the matrix u
l
respectively v
l
. Note that
D
l
=

q
D
l
q
and D

l
=

q
D

l
q
. We have already proven, see (22), that
[

C
p
l
]
(q,q);(r,r)
=
q,p

q,r
[c
p
l
]
q,r
, [

C
p
l
]
(q,q);1
=
q,p
[c
p
l
]
q,1
, (28)
where c
p
l
C
D

l
p
D
l+1
p
. Finally, if we substitute (25) in (19a), we obtain
(e
i[(s+(1)
l
)/2+qr]
1)[

B
s
l
]
(q,q);(r,r)
= 0, (l > 1), (e
i[(s1)/2r]
1)[

B
s
l
]
1;(r,r)
= 0 (29)
meaning that
[

B
s
l
]
(q,q);(r,r)
=
r,q+(s+(1)
l
)/2
[b
s
l,q
]
q,r
, [

B
s
1
]
1;(r,r)
=
r,(s1)/2
[b
s
1,0
]
1,r
(30)
where b
s
l,q
C
D
l
q
D

l
q+(s+(1)
l
)/2
is random.
We have now proven that every MPS that is invariant under local gauge transformations with
l
= 1 can be brought
in the form (28) and (30) by a MPS gauge transformation. A state in this form is also invariant under any gauge
transformation. Indeed, according to (17), we need to nd unitary matrices U
l
and V
l
such that
e
i
l
p

C
p
l1
=

C
p
l1
U
l
, e
i
l
p

C
p
l
= V

C
p
l
, e
i
l
(s+(1)
l
)/2

B
s
l
= U

B
s
l
V
l
, (31)
where

B equals (30) and

C equals (28). Taking
[U
l
]
(q,q);(r,r)
=
q,r

q,r
e
i
l
q
, [V
l
]
(q,q);(r,r)
=
q,r

q,r
e
i
l
q
, (32)
solves this problem. This proves that every gauge-invariant state can be brought in the form (28) and (30) by a
MPS-gauge transformation and, conversely, that every MPS in the form (28) and (30) is gauge invariant.
A2: A MPS ansatz for CT-invariant systems in the thermodynamic limit
Consider a one-dimensional lattice of size 4N where every site n, n 2N + 1, . . . , 2N, contains a ddimensional
Hilbert space 1
n
spanned by the basis [q
n

n
: q
n
= 1, . . . , d. The total Hilbert space is spanned by [q
[q
2N+1
. . . q
2N
: q
n
= 1, . . . , d. We will take the thermodynamic limit (N +). Let H be a Hamiltonian which
can be written as
H =

nZ
(CT)
n
m

k=1
_
h
(k)
1
(Ch
(k)
2
C)
_
(CT)
n
, (33)
where h
(k)
i
has only support on one site (i = 1, 2; k = 1, . . . , m), C =
nZ
C, C is an idempotent Hermitian operator
which induces a permutation c on the basis vectors (C [q
n

n
= [c(q
n
)
n
, c
2
= 11) and T is the translation operator. One
can think of C being the charge conjugation. It is clear that the Hamiltonian is invariant under the transformation
CT CT : H = (CT)H(CT)

. Further one notes that the Hamiltonian is invariant under translations over an even
number of sites. As a consequence it is possible to label the eigenstates of the Hamiltonian by the quantum numbers
k [, ) and 1, +1: H[k, = E
k,
[k, where CT [k, = e
ik/2
[k, , T
2
[k, = e
ik
[k, . The
number k corresponds to the momentum of the excitation (for translations over two sites). The excitations with
9
quantum number = +1 will be referred to as scalar particles and the excitations with quantum number = 1
will be referred to as vector particles.
When the ground state of the Hamiltonian H does not suer from spontaneous symmetry breaking of the
CT-symmetry, we can write down an ansatz which resembles a uniform MPS but is CT-invariant instead of
translation invariant. Thereto we dene [q
c
[c(q
2N+1
), q
2N+2
, . . . , c(q
2n1
), q
2n
, . . . , c(q
2N1
), q
2N
(N +)
which can be obtained by letting C act on the odd components of [q. The MPS [(A) is now dened as
[(A)

{qn}
nZ
v

L
_

nZ
A
qn
_
v
R
[q
c
, A C
DdD
, v
R
, v
L
C
D1
. (34)
One easily veries that this state is CT-invariant and that it can be obtained from the uniform MPS
[
u
(A)

{qn}
nZ
v

L
_

nZ
A
qn
_
v
R
[q , (35)
by letting C act on the odd sites.
Once we have obtained (to sucient accuracy) the ground state [(A) corresponding to the ground state energy
density E
0
one can look for the excited states. For one-particle excited states an ansatz with quantum numbers k and
(k [, ), 1, +1) is
[
k,
(B, A) =

nZ

n
e
i(k/2)n
d

{q}=1
v

L
_

m<n
A
qm
_
B
qn
_

m>n
A
qm
_
v
R
[q
c
. (36)
It is not hard to see that CT [
k,
(B, A) = e
ik/2
[
k,
(B, A) and T
2
[
k,
(B, A) = e
ik
[
k,
(B, A). These
states can be obtained by letting C act on the odd sites of the states [
u
([k+(1)]/2)
(B, A) where
[
u
l
(B, A) =

nZ
e
iln
d

{q}=1
v

L
_

m<n
A
qm
_
B
qn
_

m>n
A
qm
_
v
R
[q , l [, [. (37)
The states [
u
l
(B, A) were introduced in [23] as ansatz for momentum-l particles for translational invariant systems.
To obtain the ground state of H with ansatz [(A), we will use the TDVP method [20]. It can be ver-
ied that this is equivalent with applying the TDVP method to the translational invariant Hamiltonian
H
u

nZ
T
n

m
k=1
_
h
(k)
1
h
(k)
2
_
T
n
where we take as ansatz [
u
(A). The procedure is explained in [21],
[23]. The A that we will obtain as the tensor corresponding to the ground state [
u
(A) of H
u
will also correspond
to the ground state [(A) of H.
To nd the excited states we will apply the Rayleigh-Ritz method and nd B in such a way that it minimizes

k,
(B

, A

)[H[
k,
(B, A) /
k,
(B

, A

) [
k,
(B, A). By noting that

k,
(B

, A

)[H[
k,
(B, A)

k,
(B

, A

) [
k,
(B, A)
=

u
([k+(1)]/2)
(B

, A

)[H
u
[
u
([k+(1)]/2)
(B, A)

u
([k+(1)]/2)
(B

, A

) [
u
([k+(1)]/2)
(B, A)
(38)
this problem is mapped to an analogue problem for uniform MPS. It is extensively discussed in [20], [21] how to deal
with this and how to obtain an approximation for the excited states and their energy by minimizing the right-hand side.
For the Schwinger model, the Hamiltonian can be written as
H =
g
2

nZ
(CT)
n
_
L(1)
2
+

2
_

z
(1) + 1) +x(

(1)e
i(1)
[C

(2)C] +h.c.
_
_
(CT)
n
, (39)
where C is the charge conjugation: C [s, q = [s, q. This implies that we have to apply the TDVP method for
uniform MPS to the translational invariant Hamiltonian
H
u
=
g
2

nZ
T
n
_
L(1)
2
+

2
_

z
(1) + 1) +x(

(1)e
i(1)

(2) +h.c.
_
_
T
n
. (40)
10
We will now construct an ansatz of the form (34) which is gauge invariant. We start from a MPS invariant under
translations over an even number of sites and perform a charge conjugation on the odd sites:

{sn}=1

{pn}Z
v

L
_
N

n=N+1
B
s2n1
1
C
p2n1
1
B
s2n
2
C
p2n
2
_
v
R
[s
2n1
, p
2n1
, s
2n
, p
2n
, (41)
where N +. To make the state gauge invariant, we will require that they have the form (28) and (30):
[B
(1)
l
s
l
]
(q,q);(r,r)
=
r,q+(s+(1)
l
)/2
[b
s
l,q
]
q,r
, [C
(1)
l
p
l
]
(q,q);(r,r)
=
q,p

q,r
[c
p
l
]
q,r
, (42)
where l = 1, 2. We will now perform the following MPS-gauge transformation: B
s
1


B
s
1
= UB
s
1
, C
p
1


C
p
1
= C
p
1
,
B
s
2


B
s
2
= B
s
2
and C
p
2


C
p
2
= C
p
2
U

where [U]
(p,);(q,)
=
p,q

,
. It follows that
[

B
s
1

C
p
1
]
(q,);(r,)
= b
s
1,q
c
q(s+1)/2
1

p,q+(s+1)/2

r,q(s+1)/2
, (43)
[

B
s
2

C
p
2
]
(q,);(r,)
= b
s
2,q
c
q+(s+1)/2
2

p,q+(s+1)/2

r,q(s+1)/2
. (44)
By taking b
s
1,q
= b
s
2,q
b
s
q
, c
p
1
= c
p
2
c
p
and dening [A
s,p
]
(q,q);(r,r)
= [a
{p,s}
]
q,r

p,q+(s+1)/2

r,q(s+1)/2
where
a
{p,s}
= b
s
p(s+1)/2
c
p
, the state (41) can be written as
[(A) =

{sn}=1

{pn}Z
v

L
_
2N

n=2N+1
A
sn,pn
_
v
R
[s
2k1
, p
2k1
, s
2k
, p
2k
, (45)
which is a gauge- and CT-invariant ansatz. Once we have obtained the ground state, we can use the ansatze (36)
to approximate the excited states. If we put [B
s,p
]
(q,q);(r,r)
= [b
{p,s}
]
q,r

p,q+(s+1)/2

r,q(s+1)/2
they will also be
automatically gauge invariant. The variational freedom then lies within the matrices b
{p,s}
. Note that in (36) the
energy corresponding to the physical momentum k is obtained by substituting k 2ka.