# Tensor lattice field theory with applications to the renormalization group
and quantum computing

Yannick Meurice1 Ryo Sakai1 Judah Unmuth-Yockey2,3 1 Department of Physics and
Astronomy, The University of Iowa, Iowa City, IA 52242 USA  2 Department of
Physics, Syracuse University, Syracuse, NY 13244 USA  3 Fermi National
Accelerator Laboratory, Batavia, IL 60510 USA

(August 27, 2024)

###### Abstract

We discuss the successes and limitations of statistical sampling for a
sequence of models studied in the context of lattice QCD and emphasize the
need for new methods to deal with finite-density and real-time evolution. We
show that these lattice models can be reformulated using tensorial methods
where the field integrations in the path-integral formalism are replaced by
discrete sums. These formulations involve various types of duality and provide
exact coarse-graining formulas which can be combined with truncations to
obtain practical implementations of the Wilson renormalization group program.
Tensor reformulations are naturally discrete and provide manageable transfer
matrices. Combining truncations with the time continuum limit, we derive
Hamiltonians suitable to perform quantum simulation experiments, for instance
using cold atoms, or to be programmed on existing quantum computers. We review
recent progress concerning the tensor field theory treatment of non-compact
scalar models, supersymmetric models, economical four-dimensional algorithms,
noise-robust enforcement of Gauss’s law, symmetry preserving truncations and
topological considerations. We discuss connections with other tensor network
approaches.

###### Contents

  1. I Introduction
  2. II Lattice field theory
    1. II.1 The “Kogut sequence”: from Ising to QCD
    2. II.2 Classical lattice models and path integral
    3. II.3 Physical applications
    4. II.4 Computational methods beyond perturbation theory
  3. III Quantum computing
    1. III.1 Situations where importance sampling fails
    2. III.2 Qubits and other quantum platforms
    3. III.3 From Euclidean transfer matrices to Hilbert spaces
    4. III.4 Topological and geometrical dualities
    5. III.5 Real-time evolution with qubits
    6. III.6 Lloyd-Suzuki-Trotter product formula
    7. III.7 Dealing with noise in the NISQ era
    8. III.8 Quantum computations and simulations
      1. III.8.1 Ising model
      2. III.8.2 Gauge theories
  4. IV Quantum versus classical
    1. IV.1 Models
    2. IV.2 Phase transitions
    3. IV.3 Tensor networks
  5. V Tensor methods explained with the Ising model
    1. V.1 Tensor formulation
    2. V.2 Remarks on the forms of duality
    3. V.3 Boundary conditions
    4. V.4 Exact blocking
  6. VI Tensor renormalization group
    1. VI.1 Block-spinning through SVD
    2. VI.2 Optimized truncations
    3. VI.3 Higher dimensional algorithms
    4. VI.4 Observables with tensors
    5. VI.5 Niemeijer-Van Leeuwen equation
    6. VI.6 A simple example of TRG fixed point
    7. VI.7 Corner double line structure on tensor network
  7. VII Tensors for spin models with an Abelian symmetry
    1. VII.1 $O(2)$ nonlinear sigma model
    2. VII.2 $q$-state clock models
    3. VII.3 Dual reformulations with unconstrained variables
    4. VII.4 Chemical potential, complex temperature, and importance sampling
  8. VIII Tensors for spin models with non-Abelian symmetries
    1. VIII.1 O(3) nonlinear sigma model
    2. VIII.2 $SU(2)$ principal chiral model
    3. VIII.3 Truncations and asymptotic freedom
  9. IX Tensors for lattice gauge theories
    1. IX.1 Pure gauge $U(1)$
      1. IX.1.1 Discrete Maxwell equations
      2. IX.1.2 Abelian gauge duality
    2. IX.2 The compact Abelian Higgs model
    3. IX.3 $SU(2)$ gauge theory
    4. IX.4 The non-Abelian Higgs model
  10. X Tensors for models with (non-compact) scalars
    1. X.1 Real $\phi^{4}$ theory
    2. X.2  Tensors from Gaussian quadrature
    3. X.3 Additional topics and references
  11. XI Models with fermions
    1. XI.1 Tensor representation for free Wilson fermions
    2. XI.2 Grassmann tensor renormalization group
    3. XI.3 2D Schwinger model with Wilson fermions
    4. XI.4 3D free fermions
    5. XI.5 2D $\mathcal{N}=1$ Wess–Zumino model
    6. XI.6 2D Schwinger model with staggered fermions
    7. XI.7 Additional topics and references
  12. XII Transfer matrix and Hamiltonian
    1. XII.1 Spin models
    2. XII.2 Quantum simulations for the O(2) and O(3) model
    3. XII.3 Gauge models
    4. XII.4 Duality revisited and Gauss’s law
    5. XII.5 Quantum simulation of the Abelian Higgs model with cold atoms
    6. XII.6 Algebraic aspects of the Hamiltonian formulation
    7. XII.7 Additional topics and references
  13. XIII Additional aspects
    1. XIII.1 Symmetries and truncations
    2. XIII.2 Topological considerations
    3. XIII.3 Quantum gravity
  14. XIV Conclusions
  15. A Review of mathematical results
    1. A.1 Character expansions
    2. A.2 Orthogonality & completeness
    3. A.3 Singular value decomposition

##  I Introduction

Quantum field theory models on space or space-time lattices play an important
role in our understanding of strongly interacting particles, nuclei,
superconductivity, condensed matter and phase transitions. In high-energy and
nuclear physics, lattice quantum chromodynamics (QCD) provides an ab-initio
theory of strong interactions. In the QCD context, the lattice is a non-
perturbative ultraviolet regularization which preserves local gauge
invariance. As the lattice is not a physical feature, we need to approach the
continuum limit where the lattice spacing is small compared to the physical
length scales involved in the problem.

In the context of solid-state physics, lattice spacings of the order a few
Angstroms are present, and atomic physicists can create optical lattices with
a lattice spacing on the order of the laser wavelength; however, large
correlation lengths appear near critical points and universal behavior
independent of the microscopic details can be observed. It is also possible to
create actual optical lattices in laboratories by trapping cold atoms in
counter-propagating laser beams, and tune the interaction in order to quantum
simulate lattice models with interactions similar to Hubbard models Bloch  _et
al._ (2008). This can be called an analog computing method or a quantum
simulation experiment. Again, it is possible to tune the parameters to reach
universal behaviors related to quantum phase transitions with large
correlation lengths.

More generally, we are getting better control in the manipulation of small
quantum systems evolving in small Hilbert spaces and the idea of using
physical quantum systems to study theoretical quantum models Feynman (1982)
has generated many exciting developments Georgescu  _et al._ (2014). As
bits—that can be either on or off—can be thought of as the basic building
blocks of classical computers, one can envision qubits that can each be used
as a two-dimensional Hilbert space as the building blocks of a quantum
computer. If we want to use the $2^{N}$ dimensional Hilbert space provided by
$N$ qubits to represent the Hilbert space of a quantum field theory problem we
need to apply discretizations and truncations. Discretization of space can be
achieved by the lattice approximation, while the discretization of continuous
field integration can be done by using character expansion, as will be
discussed extensively in this review, or by other methods which include gauge
magnets/quantum links Orland and Rohrlich (1990); Horn (1981); Brower  _et
al._ (1999a); Wiese (2013), or field digitization, or summation of discrete
subgroups Jordan  _et al._ (2014); Klco and Savage (2019); Lamm  _et al._
(2019); Alexandru  _et al._ (2019a); Hackett  _et al._ (2019).

General arguments Lloyd (1996) show that for local interactions, a quantum
computer will reduce the computational effort, for problems like the real-time
evolution, to a polynomial in the size of the system rather than an
exponential for a classical computer.

From a purely theoretical point of view, studying models with a large number
of strongly correlated degrees of freedom is very challenging. In order to
deal with this situation, L. Kadanoff Kadanoff (1966) suggested to consider
the average field or spin in cells of variable sizes often called “blocks”.
The procedure is often called “blockspinning” and it played a crucial role in
the development of the renormalization group (RG) ideas Wilson and Kogut
(1974). Sometimes great theoretical intuitions can take a long time to be
practically realized. Despite its visual appeal, the blockspinning procedure
is not easy to implement numerically. It typically involves approximations
that are difficult to improve.

Successful applications of the RG idea were made possible without requiring
numerical implementations of the original blocking idea. A well-known example
is the discovery of asymptotic freedom Gross and Wilczek (1973); Politzer
(1973) which initially relied on a one-loop calculation of the Callan-Symanzik
beta function. Typically, the interplay between small and large energy scales
is more easily seen in the momentum representation. However, it is clear that
$if$ it was possible to design practical methods such that each step of the
blocking could be performed within a reasonable amount of time and with a
desired accuracy, the computational cost would scale like the number of
blockings, in other words, the logarithm of the volume, since the size of a
block after each blocking step doubles Kadanoff (1966). Achieving this goal is
nontrivial and not guaranteed in general. Below, we briefly explain the
practical issues that have prevented the use of blocking for quantitative
purposes and how new tensorial methods can be used to make progress in this
direction.

A simple way of blocking consists in introducing 1 in the partition function
in the following generic form:

$\prod_{\mathfrak{B}}\int
d\Phi_{\mathfrak{B}}\delta(\Phi_{\mathfrak{B}}-\sum_{x\in{\mathfrak{B}}}\phi_{x})=1,$
(1)

where the blocks ${\mathfrak{B}}$ form a partition of the original lattice.
For instance, on a three-dimensional cubic lattice, the blocks can be chosen
as cubes with a linear size of two lattice spacings and contain eight sites.
The $\phi_{x}$ are the original lattice fields and $\Phi_{\mathfrak{B}}$ the
block fields which inherit new effective interactions after one performs the
integration over the original fields. Conceptually, this sounds easy, however
in practice it appears to be more complicated than the original problem. As an
example, one can try to write a simple algorithm for the two-dimensional Ising
model on a square lattice by replacing four spins in a $2\times 2$ square
block by a single variable and write an effective energy function (or at least
some effective measure) for the new block variables. The procedure becomes
more intricate as we proceed, and finding the effective energy function is
nontrivial Liu  _et al._ (2013).

For this reason, approximate procedures were developed such as the Migdal-
Kadanoff approximation Migdal (1975); Kadanoff (1976), the approximate
recursion formula Wilson and Kogut (1974) or other hierarchical approximations
Dyson (1969); Baker (1972); Meurice (2007), where no new interactions are
generated by the blocking process. However, in these examples, the lack of
reference to an exact procedure to handle the original model with localized
interactions makes the systematic improvement of these approximations
difficult. Similar issues appear in non-perturbative functional methods based
on the momentum space representation Berges  _et al._ (2002), where the local
potential approximation allows high-accuracy estimates of the critical
exponents Bervillier  _et al._ (2007), but its improvement with methods such
as the derivative expansion remains difficult Bervillier (2013). For Ising
models, it is possible to deal with the proliferation of couplings generated
by the blocking process by starting with the most general set of interactions
Kadanoff (1975); Kadanoff and Houghton (1975); Niemeijer and van Leeuwen
(1976). They introduce the identity in terms of probabilities
$P(\\{\sigma^{\prime}\\},\\{\sigma\\})$ such that

$\sum_{\\{\sigma^{\prime}\\}}P(\\{\sigma^{\prime}\\},\\{\sigma\\})=1,$ (2)

where $\\{\sigma^{\prime}\\}$ are new Ising spins associated with blocks. As
we will discuss later, this special setup allows us to write formal
expressions for the effective couplings as double partition functions and
write RG equations. However, from a computational point of view the locality
of the interactions is lost and additional assumptions are needed to proceed.

In contrast, reformulations of the partition function of classical spin models
as the trace of a product of local tensors provide a new type of blocking
procedure in configuration space called Tensor RG (TRG) Nishino and Okunishi
(1996); Levin and Nave (2007); Gu and Wen (2009); Xie  _et al._ (2009); Gu
_et al._ (2010); Xie  _et al._ (2012).

TRG procedures can be obtained by applying truncations to exact blocking
formulas. The blocking neatly separates the degrees of freedom inside the
block (which are integrated over), from those kept to communicate with the
neighboring blocks Meurice (2013). However, the remaining degrees of freedom
after the blocking are still microscopic and finding a truncation that
captures the low-energy physics and the entanglement is a non-trivial task.

At early stages of the TRG development, Singular Value Decomposition (SVD)
methods were used extensively. This is reviewed in Ref. Efrati  _et al._
(2014). Some of the SVD procedures can be simplified by using character
expansions Liu  _et al._ (2013) when applied to most models studied in the
context of lattice gauge theory Liu  _et al._ (2013); Yu  _et al._ (2014);
Denbleyker  _et al._ (2014); Zou  _et al._ (2014); Shimizu and Kuramashi
(2014a); Takeda and Yoshimura (2015); Bazavov  _et al._ (2015); Shimizu and
Kuramashi (2018); Yoshimura  _et al._ (2018); Nakamura  _et al._ (2019);
Kuramashi and Yoshimura (2019); Kadoh  _et al._ (2018); Unmuth-Yockey (2019);
Kadoh  _et al._ (2019); Bazavov  _et al._ (2019); Kadoh  _et al._ (2020); Butt
_et al._ (2019). Tensorial methods are also used in the context of quantum
gravity Perez (2013); Dittrich  _et al._ (2016); Asaduzzaman  _et al._ (2019).

In general, tensorial methods represent a new approach to lattice field theory
that we will call Tensor Lattice Field Theory (TLFT ). TLFT can be used for
purposes more general than the blocking procedure. In particular TLFT provides
a very convenient, discrete framework to perform quantum computations or
simulations. There are—additionally—continuous tensor network, or tensor-like,
methods such as continuous matrix product states Haegeman  _et al._ (2013);
Campos  _et al._ (2019).

In this review article, we introduce TLFT for lattice models studied in the
context of lattice gauge theory and report progress made for blocking and
quantum computing, two “competing” methods that attempt to reduce the
computing time logarithmically. The models targeted are introduced in Sec. II.
We advocate a road map starting with the Ising model and culminating with QCD,
which we call the “Kogut sequence” Kogut (1979, 1983). This sequence is
sometimes called a “ladder” and has been followed successfully in situations
where importance sampling methods such as the Metropolis algorithm are
effective. Lattice QCD has become a very reliable precision tool to study the
static properties of hadrons. As we are writing this article, we are roughly
in the middle of the sequence, _i.e._ we have roughly half of the models in
the sequence remaining to be studied thoroughly, with improvements and
optimization on previous models still possible. However, we hope that recent
progress on higher dimensional algorithms Adachi  _et al._ (2019); Kadoh and
Nakayama (2019) could be combined with the methods that we describe to deal
with gauge fields, fermions, and non-Abelian (non-commuting) symmetries in
order to attempt calculations directly related to lattice QCD in the coming
years.

In Sec. III, we discuss situations where importance sampling cannot be used
and where quantum computations or simulations could provide alternate ways to
perform computations. This includes unitary, real-time evolution and other
situations where a sign problem is encountered. One important long-term goal
with potential impact on the interpretation of high-energy collider data is
doing ab-initio real-time calculations relevant to fragmentation processes and
parton distribution functions. In other words, starting with lattice QCD, we
would like to be able to perform calculations that would ultimately replace
the use of event generators such as Pythia Sjostrand  _et al._ (2015).

The simplest starting point for the real-time evolution is the evolution
operator $\exp(-i\hat{H}t/\hbar)$ acting on the Hilbert space of the quantum
Hamiltonian $\hat{H}$. We provide a first look at the transfer matrix which
smoothly connects the “classical” Lagrangian approach to the Hilbert space
used in the Hamiltonian formalism. We discuss various types of dualities
(geometrical and topological) that are often used together and taken for each
other.

For the models in the Kogut sequence, the bosonic field variables and the
symmetry groups are compact. General mathematical theorems, namely the
Pontryagin duality Pontryagin (1939) and the Peter-Weyl theorem Peter and Weyl
(1927), guarantee that functions over compact groups can be expanded in terms
of discrete sums of representations. This is called the “character expansion”
and was exploited to calculate strong coupling expansions Balian  _et al._
(1975) or introduce new variables on geometrically dual lattice elements Savit
(1980).

The discreteness of the character expansion provides a natural starting point
for building approximate reformulations of lattice models suitable for quantum
computing or quantum simulation experiments. The Ising model is an elementary
example where the Hilbert space of the transfer matrix can be implemented with
a set of qubits, the basic components of actual quantum computers which exist
in a linear superposition of two states, $\ket{0}$ and $\ket{1}$, rather than
being just “on” or “off” like the bits of a classical computer. For models
with continuous fields, character expansions allows us to perform the “hard
integrals” analytically without the need of approximate numerical
discretizations which break the continuous symmetries. Demonstrating the power
of the character expansion is one of the main goals of this review. Examples
of quantum computations and simulations are provided at the end of the
section. In Sec. IV, we clarify the use of “classical” and “quantum” in
various contexts and make connections with other approaches Schollwöck (2011);
Ran  _et al._ (2020); Haegeman and Verstraete (2017); Cirac  _et al._ (2020).

Sec. V introduces the tensor reformulation for the Ising model. SVD,
truncation and the TRG method are discussed in Sec. VI. Spin models with an
O(2) symmetry or with discrete subgroups are discussed in Sec. VII. In Sec.
VIII, we derive expressions for local tensors in the simple case of a non-
Abelian spin model with O(3) symmetry. We also find tensor expressions for
effective theories of gauge theories known as principal chiral models.

Models with local gauge symmetry are introduced in Sec. IX. We first consider
Abelian gauge theories and work up in complexity to tensor expressions for
non-Abelian gauge theories as well.

In Sec. X, tensor network expressions for the real and the complex $\phi^{4}$
theory are derived. For models with non-compact fields such as the scalar
$\phi^{4}$ theory, the Gaussian quadrature rule can be used to extract
discrete degrees of freedom, just as the gauge degrees of freedom are
discretized via the character expansions. The accuracy of the tensor network
approach is shown for the real-field case, and the ability to deal with a
severe sign problem is shown for the complex-field case.

In Sec. XI we present tensor formulations for models with fermionic degrees of
freedom. In general fermions fit in well with the tensor (and discrete)
approach thanks to the nilpotency of the Grassmann variables. In the section,
various models that contain fermions such as pure fermions, gauged fermions,
and fermions combined with scalars are discussed.

In Sec. XII we re-discuss the transfer matrix using the tensor formalism and
broaden the perspective. Recent TLFT developments regarding symmetries,
topological solutions and quantum gravity are discussed in Sec. XIII.

##  II Lattice field theory

###  II.1 The “Kogut sequence”: from Ising to QCD

In the early 70’s, QCD appeared as a strong candidate for a theory of strong
interactions involving quarks and gluons. However the perturbative methods
that provided satisfactory ways to handle the electroweak interactions of
leptons failed to explain confinement, mass gaps and chiral symmetry breaking.
A non-perturbative definition of QCD was needed. In 1974, K. Wilson proposed
Wilson (1974) a lattice formulation of QCD where the $SU(3)$ local symmetry is
exact. As this four-dimensional model is fairly difficult to handle
numerically, a certain number of research groups started considering simpler
lattice models in lower dimensions and then increased symmetry and
dimensionality. This led to a sequence of models, sometimes called the “Kogut
ladder” which appears in J. Kogut’s review articles Kogut (1979, 1983), and
later with small modifications in textbooks by A. Polyakov Polyakov (1987) and
C. Itzykson with J. M. Drouffe Itzykson and Drouffe (1991).

The sequence is approximately the following:

  1. 1.

$D=2$ Ising model

  2. 2.

$D=3$ Ising model and its gauge dual

  3. 3.

$D=2$ O(2) spin and Abelian Higgs models

  4. 4.

$D=2$ fermions and the Schwinger model

  5. 5.

$D=3$ and 4 $U(1)$ gauge theory

  6. 6.

$D=3$ and 4 Non-Abelian gauge theories

  7. 7.

$D=4$ lattice fermions

  8. 8.

$D=4$ QCD

This sequence should not be understood in a rigid way as if each step is
necessary for the next step. For instance steps 3. to 5. could be interchanged
and the problems involving fermions have very specific features that are not
easily compared to those involving only bosonic fields. The message that we
want to convey is that there is an approximate road map that has proven to be
very effective for the “classical” approach of lattice field theory in order
to deal with static problems using importance sampling (Monte Carlo) methods.
We advocate to follow a similar path to develop the quantum versions of these
models and deal with real-time evolution and other problems not accessible
with classical methods. The difference between quantum and classical is
explained more precisely in Sec. IV.1. In addition, a similar path is being
followed to develop numerical coarse-graining.

###  II.2 Classical lattice models and path integral

In this subsection we introduce lattice versions of classical field theory
models. At this point, we would like to point out that while we will provide
definitions of the fields used, notations, and acronyms/initials, more details
on basic quantum field theory and lattice field theory can be found in
textbooks and review articles _e.g._ Peskin and Schroeder (1995); Itzykson and
Drouffe (1991); Kogut (1979); Montvay and Münster (1994). We use a Euclidean
time and treat space and time on the same footing. The metric is simply a
Kronecker delta in $D$ dimensions. We then discretize space and time. We use a
$D$-dimensional (hyper) cubic Euclidean space-time lattice. The sites are
denoted $x=(x_{1},x_{2},\dots x_{D})$, with $x_{D}=\tau$, the Euclidean time
direction. In lattice units, the space-time sites are labelled with integers.
In the following, the lattice units are implicit. The links between two
nearest neighbor lattice sites $x$ and $x+\hat{\mu}$ are labelled by $(x,\mu)$
and the plaquettes, the smallest squares on a square or (hyper) cubic lattice,
delimited by four sites $x$, $x+\hat{\mu}$, $x+\hat{\mu}+\hat{\nu}$ and
$x+\hat{\nu}$ are labelled by $(x,\mu\nu)$. By convention, we start with the
lowest index when introducing a conventional circulation at the boundary of
the plaquette. The total number of sites is denoted by $V$. Unless specified,
periodic boundary conditions are assumed, and they preserve a discrete
translational symmetry. If we take the time continuum limit, we obtain a
quantum Hamiltonian formulation in $D-1$ spatial dimensions.

In the continuum, the Lagrangian density for $N$ real scalar fields with a
$O(N)$ global symmetry reads

$\mathcal{L}^{O(N)}_{Euclidean}=\frac{1}{2}\partial_{\mu}\vec{\phi}\cdot\partial_{\nu}\vec{\phi}\delta^{\mu\nu}+\lambda(\vec{\phi}\cdot\vec{\phi}-v^{2})^{2},$
(3)

with $\vec{\phi}$ a $N$-dimensional vector. Here $\lambda$ is a coupling
constant, whose size determines fluctuations of the $|\vec{\phi}|$ field
around some value, $v$. The potential has degenerate minima on a $N$-1
dimensional hyper-sphere $S_{N-1}$ and a local maximum at $\vec{\phi}=0$. For
$N=2$, the low energy part of the potential has a shape reminiscent of the
bottom of a wine bottle. The degenerate minima form a circle at the very
bottom. We can study the small fluctuations about a given minimum on the
circle. Note that the choice of a minimum breaks the $O(2)$ symmetry. There
are “soft” fluctuations along the circle that restore the symmetry and “hard”
fluctuations in the radial direction.

We can extend this analysis for arbitrary $N$. We have one massive mode
(fluctuations in the symmetry breaking direction) with mass
$2\sqrt{2\lambda}v$ and $N-1$ massless modes (“Nambu-Goldstone” (NG) modes)
which is the number of broken generators Goldstone  _et al._ (1962); Nambu
(1960); Goldstone (1961); Peskin and Schroeder (1995).

We can write the Euclidean action for the NG modes on a $D$-dimensional
lattice with isotropic lattice spacing $a$ as

$S_{\text{NLSM}}=\frac{1}{2}\sum_{x}\sum_{\mu=1}^{D}a^{D-2}(\vec{\phi}_{x+\hat{\mu}}-\vec{\phi}_{x})\cdot(\vec{\phi}_{x+\hat{\mu}}-\vec{\phi}_{x}).$
(4)

This model is called the nonlinear sigma model (NLSM) Peskin and Schroeder
(1995). The constraint $\vec{\phi}_{x}\cdot\vec{\phi}_{x}=v^{2}$ (which
enforces the “nonlinear” part of its name) can be expressed by introducing
unit vectors: $\vec{\phi}_{x}=v\vec{\sigma}_{x}$ such that

$\vec{\sigma}_{x}\cdot\vec{\sigma}_{x}=1.$ (5)

Redefining $a^{D-2}v^{2}\equiv\beta$, we get the simple action

$S_{\text{NLSM}}=\beta\sum_{x,\mu}(1-\vec{\sigma}_{x+\hat{\mu}}\cdot\vec{\sigma}_{x}).$
(6)

These models are often called spin models as well. The first term in the
action $\beta\sum_{x,\mu}1$ is a constant that is often dropped. However, for
large $\beta$, the configurations with almost constant $\vec{\sigma_{x}}$
dominate the partition function and since under these circumstances
$\vec{\sigma}_{x+\hat{\mu}}\cdot\vec{\sigma}_{x}\simeq 1$ it is useful to
subtract the constant in order to just keep the small fluctuations.

The case $N=1$ is the well-known Ising model with $\sigma_{x}=\pm 1$ Ising
(1925). For $N=2$, the terminology “planar model” or “classical XY model” Vaks
and Larkin (1965); Bowers and Joyce (1967) is common and if we use the circle
parametrization

$\sigma_{x}^{(1)}=\cos(\varphi_{x}),{\rm and}\
\sigma_{x}^{(2)}=\sin(\varphi_{x}),$ (7)

then

$\vec{\sigma}_{x+\hat{\mu}}\cdot\vec{\sigma}_{x}=\cos(\varphi_{x+\hat{\mu}}-\varphi_{x}),$
(8)

with $\varphi_{x}\in[0,2\pi)$.

There is another class of models which breaks the $O(2)$ symmetry in Eq. (8)
into a discrete $\mathbb{Z}_{q}$ symmetry, _i.e._ the possible angles are
restricted to those of the $q$th roots of unity in the complex plane. They are
called the $q$-state clock models Potts (1952). These models have the same
action from Eq. (6) with the identification of the angles as being discrete,

$\displaystyle\varphi_{x}=\frac{2\pi n_{x}}{q}$ (9)

with $n_{x}=0,1,2,\ldots,q-1$. With this identification, the $O(2)$ model
emerges as the $q\rightarrow\infty$ limit of the $q$-state clock models, and
the Ising model is simply the $q=2$ model.

For $N=3$, the symmetry becomes non-Abelian and the model is sometimes called
the “classical Heisenberg model”. In the large-$N$ limit, the model becomes
solvable if we take the limit in such a way that $N/\beta(N)=\lambda$ remains
constant Coleman  _et al._ (1974).

It is instructive to rewrite the $O(2)$ model using the complex form

$\Phi_{x}={\rm e}^{i\varphi_{x}}.$ (10)

Dropping the constant terms, the $O(2)$ action reads

$\displaystyle S_{O(2)}$
$\displaystyle=\frac{\beta}{2}\sum_{x,\mu}(\Phi_{x+\hat{\mu}}-\Phi_{x})\cdot(\Phi_{x+\hat{\mu}}-\Phi_{x})^{\star}$
$\displaystyle=-\beta\sum_{x,\mu}\cos(\varphi_{x+\hat{\mu}}-\varphi_{x}).$

The $O(2)$ model has a global symmetry

$\varphi_{x}\rightarrow\varphi_{x}+\alpha.$ (11)

With the complex notation, this transformation becomes

$\Phi_{x}\rightarrow{\rm e}^{i\alpha}\Phi_{x}.$ (12)

We would like to promote this symmetry to a local one

$\Phi_{x}\rightarrow{\rm e}^{i\alpha_{x}}\Phi_{x},$ (13)

_i.e._ one that is site dependent. This can be achieved by inserting a phase
$U_{x,\mu}$ between $\Phi_{x}^{\star}$ and $\Phi_{x+\hat{\mu}}$ which
transforms like

$U_{x,\mu}\rightarrow{\rm e}^{i\alpha_{x}}U_{x,\mu}{\rm
e}^{-i\alpha_{x+\hat{\mu}}}.$ (14)

The procedure can be extended for arbitrary $N$-dimensional complex vectors
${\bf\Phi}_{x}$ with a local transformation involving a $U(N)$ matrix $V_{x}$:

${\bf\Phi}_{x}\rightarrow V_{x}{\bf\Phi}_{x}$ (15)

In addition, we can also introduce $U(N)$ matrices ${\bf U}_{x,\hat{\mu}}$
transforming like

${\bf U}_{x,\hat{\mu}}\rightarrow V_{x}{\bf
U}_{x,\hat{\mu}}V_{x+\hat{\mu}}^{\dagger}.$ (16)

The action

$S_{U(N)}=-\frac{\beta}{2}\sum_{x,\mu}({\bf\Phi}_{x}^{\dagger}{\bf
U}_{x,\hat{\mu}}{\bf\Phi}_{x+\hat{\mu}}+\text{h.c.})$ (17)

has a local $U(N)$ invariance which we call gauge invariance. If we consider
two successive links in positive directions, then the local transformation at
the middle site cancels and

${\bf U}_{x,\mu}{\bf U}_{x+\hat{\mu},\nu}\rightarrow V_{x}{\bf U}_{x,\mu}{\bf
U}_{x+\hat{\mu},\nu}V_{x+\hat{\mu}+\hat{\nu}}^{\dagger}.$ (18)

If the second link goes in the negative direction, we use the Hermitian
conjugate and a similar property holds

${\bf U}_{x,\mu}{\bf U}_{x+\hat{\mu}-\hat{\nu},\nu}^{\dagger}\rightarrow
V_{x}{\bf U}_{x,\mu}{\bf
U}_{x+\hat{\mu}-\hat{\nu},\nu}^{\dagger}V_{x+\hat{\mu}-\hat{\nu}}^{\dagger}.$
(19)

We can pursue this process for an arbitrary path connecting $x$ to some
$x_{\text{final}}$. The transformation on the right will be
$V_{x_{\text{final}}}^{\dagger}$. If we close the path and take the trace, we
obtain a gauge-invariant quantity. We call these traces of products of gauge
matrices over closed loops “Wilson loops” Wilson (1974). In the case where the
loop goes around the imaginary time direction, we often call it a “Polyakov
loop” Polyakov (1978).

On a square, cubic or hypercubic lattice, the smallest path that gives a non-
trivial Wilson loop is a square. We call this square a plaquette. Claude
Itzykson coined this terminology after Ken Wilson’s seminar in Orsay in 1973.
The corresponding matrix is

${\bf U}_{\text{plaquette}}={\bf U}_{x,\mu\nu}={\bf U}_{x,\mu}{\bf
U}_{x+\hat{\mu},\nu}{\bf U}_{x+\hat{\nu},\mu}^{\dagger}{\bf
U}_{x,\nu}^{\dagger}$ (20)

The simplest gauge-invariant lattice model has an action, called Wilson’s
action:

$S_{\text{Wilson}}=\beta_{pl.}\sum_{\langle
x,\mu\nu\rangle}\left[1-\frac{1}{2N}\text{Tr}\left[{\bf
U}_{x,\mu\nu}+\text{h.c.}\right]\right].$ (21)

where $\sum_{\langle x,\mu\nu\rangle}$ indicates a sum over all plaquettes.
Here each $\bf{U}_{x,\mu}$ is related to the vector potential, or gauge field,
in the continuum theory through

$\mathbf{U}_{x,\mu}=e^{i\mathbf{A}_{x,\mu}}.$ (22)

On the lattice, both $\mathbf{U}_{x,\mu}$ and $\mathbf{A}_{x,\mu}$ are located
on a link starting at $x$ and going in the $\hat{\mu}$ positive direction. In
the Abelian case ($N=1$), the matrix reduces to a phase

$U_{x,\mu}={\rm e}^{iA_{x,\mu}},$ (23)

and there is no need to take the trace.

Another generalization of the $N=1$ expression of the complex phase given in
Eq. (10), consists of replacing $\Phi_{x}$ by a $SU(N)$ matrix ${\bf U}_{x}$.
This is called the principal chiral model Gürsey (1960); Green and Samuel
(1981); Samuel and Yee (1985); Campostrini  _et al._ (1995); Peskin and
Schroeder (1995),

$S_{\text{PCM}}=-\frac{\beta}{2N}\sum_{x,\mu}\left[\text{Tr}\left[{\bf
U}^{\dagger}_{x+\hat{\mu}}{\bf U}_{x}\right]+\text{h.c.}\right].$ (24)

This model has a _global_ rotational symmetry under the $U(N)$ group, such
that
$\mathbf{U}_{x,\mu}\rightarrow\mathbf{U}_{x,\mu}^{\prime}=V\mathbf{U}_{x,\mu}V^{\dagger}$,
just as in Eq. (16) in the case of a uniform $V$ in all of space-time.

We can also consider the $N=1$ case for Eq. (3) when the group is no longer
compact, _i.e._ $\phi_{x}$ can take on the values of any real, or complex,
number. The action on the lattice is then,

$\displaystyle
S_{\text{scalar}}=\sum_{x}\left(\frac{1}{2}\sum_{\mu=1}^{D}|\phi_{x+\hat{\mu}}-\phi_{x}|^{2}+\lambda(|\phi_{x}|^{2}-\nu^{2})^{2}\right)$
(25)

or equivalently,

$\displaystyle
S_{\text{scalar}}=\sum_{x}\left\\{\frac{1}{2}\sum_{\mu=1}^{D}\left|\phi_{x+\hat{\mu}}-\phi_{x}\right|^{2}-\frac{\mu_{0}^{2}}{2}|\phi_{x}|^{2}+\frac{\lambda^{\prime}}{4}|\phi_{x}|^{4}\right\\}$
(26)

with the substitution $\mu_{0}^{2}/2=2\lambda\nu^{2}$ and
$\lambda^{\prime}/4=\lambda$, and an overall constant is ignored.

Besides scalar fields, we will also consider fermionic fields (or Grassmann
fields on the lattice). In the case of free fermions, we can write down a
lattice action using a straightforward discretization due to Wilson Wilson
(1974),

$\displaystyle S_{\text{WD}}=\sum_{x}\bar{\psi}_{x}\left(D\psi\right)_{x}$
(27)

where the Wilson-Dirac operator is defined by

$\displaystyle D_{xx^{\prime}}=\left(am+rD\right)\delta_{x,x^{\prime}}$
$\displaystyle+\frac{1}{2}\sum_{\mu=1}^{D}\left\\{\left(r-\gamma_{\mu}\right)\delta_{x^{\prime},x+\hat{\mu}}+\left(r+\gamma_{\mu}\right)\delta_{x,x^{\prime}+\hat{\mu}}\right\\}.$
(28)

where $\gamma_{\mu}$ are the gamma matrices in $D$ dimensions, and $r$ is the
“Wilson parameter” to control species doubling. This so-called “doubling” is
the name given to the existence of extra copies of fermions—16 in the case of
four dimensions—and they must be removed or gapped-out appropriately.
$\psi_{x}$ and $\bar{\psi}_{x}$ are multi-component Grassmann variables that
_anti_ -commute, as opposed to classical bosonic variables which can be
interchanged without sign change.

There is another formulation of lattice fermions, where the different
components of the fermion fields are located at different lattice sites,
called staggered fermions. This comes about from a transformation which mixes
the fermion components and space-time components Kogut and Susskind (1975).
The action for free fermion fields is given by,

$\displaystyle S_{F}=\frac{1}{2}\sum_{x=1}^{N}\sum_{\mu=1}^{D}\eta_{x,\mu}$
$\displaystyle[\bar{\psi}_{x}\psi_{x+\hat{\mu}}-\bar{\psi}_{x+\hat{\mu}}\psi_{x}]$
(29)

where

$\displaystyle\eta_{x,\mu}=(-1)^{\sum_{\nu<\mu}x_{\nu}}$ (30)

with $x_{\nu}$ the coordinate in the $\nu$th direction. One can include gauge
fields in a gauge-invariant manner by inserting $\mathbf{U}_{x,\mu}$ like,

$\displaystyle S_{F}=\frac{1}{2}\sum_{x=1}^{N}\sum_{\mu=1}^{D}\eta_{x,\mu}$
$\displaystyle[\bar{\psi}_{x}\mathbf{U}_{x,\mu}\psi_{x+\hat{\mu}}-\bar{\psi}_{x+\hat{\mu}}\mathbf{U}^{\dagger}_{x,\mu}\psi_{x}].$
(31)

These models will appear in the rest of the review and will be very briefly
reintroduced in each section. As seen above, the matter degrees of freedom are
placed on each lattice site. When one considers a tensor network
representation of a model, integer degrees of freedom will arise at each link
on the lattice (see Fig. 1). As for lattice gauge theories where gauge degrees
of freedom are initially placed on links, the character expansion will be used
for generating integer at each plaquette. These points will be made clear in
the following sections. The simplest example is the Ising model discussed in
Sec. V.

Figure 1: Reinterpretation of physical degrees of freedom ($\phi$) to tensor
indices ($i$). After that, the partition function is expressed as a summation
of tensor indices instead of the path integral.

###  II.3 Physical applications

The sequence of models described in Sec. II.1 is designed to handle lattice
QCD which is so far our best definition of the theory describing strongly
interacting particles observed in a large number of experiments. Some of the
models discussed in Sec. II.2 are also studied in condensed matter. For
instance, the $O(2)$ model with a chemical potential can be seen as an
effective theory for the Bose-Hubbard model Sachdev (2001). $U(1)$ gauge
theories with either scalar of fermion fields are studied in the context of
superconductivity Herbut (2007). Tight-binding approximations for solids lead
to interesting lattice models, for instance a sheet of graphene can be
described with fermions on a hexagonal lattice having a dispersion relation
similar to the one for a massless Dirac fermion at half-filling. Castro Neto
_et al._ (2009). Other interesting applications include spin liquids and
topological phases Anderson (1973); Wen (2004, 2017); Zhou  _et al._ (2017);
Knolle and Moessner (2019).

We would like to point out that in the context of high-energy physics, local
Lagrangian densities have played a crucial role in the development of the
Standard Model. In the continuum limit, models compatible with relativistic
invariance, local gauge invariance and renormalizability have a very small
number of free parameters provided that the interactions among the fields and
their first derivatives are kept local. The standard model of electroweak and
strong interactions has only 18 free parameters if we ignore the QCD vacuum
angle and the additional parameters related to the masses and mixing of the
neutrinos in extensions of the standard model. It is not possible to tweak the
theory each time a new experiment is completed. This makes the standard model
a very predictive theory. In the lattice formulation, the local interactions
involving derivatives in the continuum are replaced by interactions involving
fields located on neighboring sites, links and plaquettes and can be
considered quasi-local. For these reasons, the idea that reformulations of
lattice models and their coarse-grained versions could be kept local seems to
be an important consideration.

###  II.4 Computational methods beyond perturbation theory

Besides the RG methods mentioned in the Introduction, lattice models have been
studied with a variety of analytical and numerical methods. For example,
expansions at small and large coupling for spin and gauge models Kogut (1979,
1983); Polyakov (1987); Itzykson and Drouffe (1991); Parisi (1998) have been
investigated extensively. As power series, they can be used to check numerical
calculations in their respective limits, however, it seems hard to capture the
non-perturbative effects such as the generation of a dynamical mass gap in the
continuum limit. After Wilson’s original proposal Wilson (1974) a very
suggestive Hamiltonian picture was developed by Kogut and Susskind Kogut and
Susskind (1975), however it became clear that the size of the Hilbert space
would make numerical calculations impractical.

In many practical situations, the path integral formulation uses a real
Euclidean action $S_{E}$. This allows importance sampling. A typical way to
proceed is to start with a random field configuration and then apply some
random changes on this configuration. If the new configuration has a lower
action, it is accepted. If the new configuration has an action larger by
$\delta S_{E}$, then it is accepted with a probability $\exp(-\delta S_{E})$.
The fact that the fields are only connected to a few neighboring fields makes
the calculation of $\delta S_{E}$ easy and the exploration of the important
configurations controllable.

Mathematicians are often amazed that it is possible to obtain very reliable
results with this method Villani (2012). So far, this has been the most
reliable way to capture the non-perturbative behavior while taking the
continuum limit. M. Creutz started numerical lattice gauge theory with
simulations of the gauge Ising model on a $3^{4}$ lattice using a HP 9830
calculator, programmed in Basic Creutz (2001). This was rapidly followed by
the study of other Abelian gauge theories Creutz  _et al._ (1979) and other
models. This is reviewed in C. Rebbi’s book Rebbi (1983).

Several decades after its inception, numerical lattice gauge theory has become
an area of research where corrections on the order of a few percent are
considered important. For instance, in heavy flavor physics the determination
of the mixing angle $|V_{ub}|$ is quoted as $(3.72\pm 0.16)10^{-3}$ Bailey
_et al._ (2015). This result was obtained using ensembles of lattices of
various lattice spacings and sizes, some of them as large as $64^{3}\times
192$. For a review of averages of numerical results with estimated errors see,
for instance, the Flavour Lattice Averaging Group (FLAG) summary Aoki  _et
al._ (2020).

Given the success of importance-sampling for static problems in QCD, not much
has been done to develop efficient numerical Hamiltonian methods. On the other
hand, in condensed matter, the idea that it is possible to approximate low-
energy states of the Hamiltonian with a significantly reduced Hilbert space
has been very successful in 1+1 dimensions and has some similarities with
importance sampling. It appeared in the context of the density matrix
renormalization group (DMRG) method White (1992); Schollwock (2005);
Schollwöck (2011); Haegeman and Verstraete (2017) and other methods discussed
in Sec. IV.3. We expect that in the coming years great progress will be made
adapting related methods in the context of lattice gauge theory.

##  III Quantum computing

###  III.1 Situations where importance sampling fails

The numerical successes of lattice gauge theory can be linked to the fact that
when the action of the Euclidean-time path integral $S_{E}$ is real,
importance sampling works remarkably well for the selection of configurations
with a Boltzmann distribution $\exp(-S_{E})$. However, if we return to real
time or introduce a chemical potential that makes the action complex, this
powerful tool becomes ineffective de Forcrand (2010). This is because in these
instances the Boltzmann weights become complex making their interpretation as
probability weights impossible.

Real-time evolution can be set up by using a Hamiltonian acting on a Hilbert
space which can be constructed by noticing that the Euclidean path-integral
can be recast as the trace of a transfer matrix Wilson and Kogut (1974);
Creutz (1977); Fradkin and Susskind (1978); Kogut (1979); Luscher (1990).
Outside of the transfer-matrix method, one can also extract real-time
observables, such as Green’s functions, from generating functional methods in
Minkowski spacetime. There are also other techniques for calculating off-
equilibrium real-time observables requiring the use of the Keldysh contour
Keldysh (1965).

Analysis of the exponential decays of correlations as the Euclidean time is
increased is a standard tool to extract masses and form factors for momenta
that are small compared to the lattice cutoff. On the other hand, real-time
evolution and a large Hilbert space involving states with large momenta are
needed to describe aspects of hadron fragmentation in real space-time and
deep-inelastic scattering Bauer  _et al._ (2019); Mueller  _et al._ (2020);
Lamm  _et al._ (2020). The real-time evolution operator does not provide a
positive measure nor a projection onto a small Hilbert space. Another
situation where sampling methods are challenged is the construction of
interpolating operators for nuclei as the number of Wick contractions grows
rapidly with the number of light quarks involved Detmold and Orginos (2013).

One could consider the possibility of abandoning the Lagrangian formulation
and directly considering the Kogut-Susskind Hamiltonian Kogut and Susskind
(1975) for QCD and/or the standard many-body formalism in condensed matter and
nuclear physics Fetter and Walecka (2003) where the Hilbert space is generated
by creation operators on a Fock-space vacuum. For the sake of argument, let us
consider the simple case of $N$ conjugate pairs of fermionic creation and
annihilation operators. They generate a Hilbert space of dimensions $2^{N}$.
The exponential growth of this number with $N$ rapidly restricts our ability
to store or manipulate the matrix elements of operators using brute force
methods with classical computers. Computations involving spatial lattices with
$64^{3}$ sites, as made possible in the Lagrangian formulations, would be
completely out of the question if we had to set up the whole quantum Hilbert
space with existing classical computers. Note however, that for a broad class
of many-body problems, it has been shown Poulin  _et al._ (2011) that most
states of the Hilbert space can only be reached after an exponentially long
time. This suggests that tensor networks methods should allow expressing
Hamiltonians and ground-state without needing to explore the full Hilbert
space.

###  III.2 Qubits and other quantum platforms

The building blocks of an ordinary (classical) computer memory are bits taking
the values 0 or 1. For instance, in some designs of Dynamic Random Access
Memory (DRAM) devices, this is achieved by using small capacitors that are
either on or off. The typical units of capacitance are femtofarads with
voltages of the order of 1V. These capacitors store a few thousand electrons
and can be charged or discharged in times of order $10^{-15}$ sec. It is
interesting to consider the limit of miniaturization where electrons, atoms or
photons could be manipulated individually and where the peculiarities of the
quantum behavior become prominent Lloyd (1996).

Following the physical examples of the electron spin or the photon
polarization, one could envision some ideal and generic quantum system, where
the on/off concept for bits is replaced by the linear superposition of two
states. This basic unit of quantum computing is often called a qubit
Schumacher (1995) and can be represented as

$\ket{\rm qubit}=\alpha_{0}\ket{0}+\alpha_{1}\ket{1},$ (32)

with two complex numbers $\alpha_{0}$ and $\alpha_{1}$ such that
$|\alpha_{0}|^{2}+|\alpha_{1}|^{2}=1$. A set of $N$ qubits spans a Hilbert
space of dimension $2^{N}$. If we use a classical computer to apply a dense
unitary matrix representing the real time evolution on an arbitrary state, we
need to perform on the order of $2^{2N}$ arithmetic operations. On the other
hand, if the Hamiltonian is such that any qubit is only connected to a
restricted number of other qubits which is fixed by the dimensionality of
space and the internal symmetries, then it is possible to design a method that
performs the evolution with a number of quantum manipulations which only
scales polynomially with $N$ Lloyd (1996). This will be discussed in Sec.
III.6. Nowadays, commercial quantum computers provide sets of qubits and the
possibility of preparing, evolving and measuring quantum states. The
limitations of the current hardware are discussed in Sec. III.7.

More generally, the idea that quantum devices could be used to perform
computations for quantum problems involving many degrees of freedom is very
appealing Feynman (1982). Physical systems involving cold atoms (see Jaksch
_et al._ (1998) and Bloch  _et al._ (2008) for a review of the early
developments) or trapped ions (see Leibfried  _et al._ (2003) for early
developments and Debnath  _et al._ (2016) for a recent example) can be used to
mimic the behavior of simplified many-body models such as various types of
spin chains or Hubbard models. In addition, superconducting circuits Devoret
and Schoelkopf (2013), Rydberg atoms Bernien  _et al._ (2017); Wu  _et al._
(2020) and photonic systems Wang  _et al._ (2019) also provide interesting
opportunities.

###  III.3 From Euclidean transfer matrices to Hilbert spaces

For the lattice models introduced in Sec. II.2, the Euclidean time was treated
on equal footing with the $D-1$ spatial dimensions. In order to discuss real-
time evolution, we first need to single-out the time direction. Evolution then
occurs along this direction according to a transfer matrix. The key ingredient
then is to connect the Lagrangian formulation to the Hamiltonian formalism and
to identify the Hilbert space from the transfer matrix introduced in Eq. (33).
The general idea is to organize the partition function sums or integrals into
operations performed on successive time slices Wilson and Kogut (1974); Creutz
(1977); Fradkin and Susskind (1978); Kogut (1979). A general procedure to
construct the transfer matrix of lattice models in configuration space in the
general context of scattering theory is nicely presented in Luscher (1990). A
“dual” method based on character expansions which are at the heart of TLFT are
discussed in Sec. XII. We will illustrate these two possibilities with
examples. Note that, in this subsection, Euclidean time methods are used to
derive a Hamiltonian that can later be used to do real-time calculations.

Turning to the transfer matrix, with generic notations for a lattice model
with $N_{\tau}$ sites in the Euclidean time direction,

$Z=\int{\mathcal{D}}\Phi\,{\rm
e}^{-S[\Phi]_{E}}=\text{Tr}(\mathbb{T}^{N_{\tau}}).$ (33)

If the lattice spacing $a_{\tau}$ in the Euclidean time direction is small
compared to the physical time scales involved, we have

$\mathbb{T}\propto e^{-a_{\tau}\hat{H}}.$ (34)

For $N_{\tau}$ large enough, the use of Euclidean time provides a projection
in the low energy sector of the Hilbert space. This property remains effective
if we insert operators that create and destroy states with non-trivial quantum
numbers.

The simplest possible example is the one-dimensional Ising model,

$\displaystyle S_{\text{Ising}}$
$\displaystyle=\beta\sum_{\tau}(1-\sigma_{\tau+1}\sigma_{\tau})$ (35)
$\displaystyle=\frac{\beta}{2}\sum_{\tau}(\sigma_{\tau+1}-\sigma_{\tau})^{2}$
(36)

with partition function

$\displaystyle Z=\sum_{\\{\sigma\\}}e^{-S}.$ (37)

This is a product of exponentials which each share one spin variable with the
next factor. Then we can write the partition function as

$\displaystyle Z=\text{Tr}[\mathbb{T}^{N_{\tau}}]$ (38)

with

$\mathbb{T}_{\alpha\alpha^{\prime}}=\exp\left\\{-\frac{\beta}{2}(\sigma^{(\alpha)}_{\tau+1}-\sigma^{(\alpha^{\prime})}_{\tau})^{2}\right\\}$
(39)

and $\sigma^{(\alpha)}=1,-1$ for $\alpha=0,1$, respectively. Along the
diagonal of the transfer matrix we see only unity. On the off-diagonal, a spin
flip comes with weight $e^{-2\beta}$. Then to leading order in the temporal
lattice spacing, $\mathbb{T}\simeq 1-a_{\tau}\hat{H}+\ldots$ which allows us
to identify

$\displaystyle\hat{H}=-h_{x}\hat{\sigma}^{x},$ (40)

with $h_{x}\equiv e^{-2\beta}/a_{\tau}$ and $\hat{\sigma}^{x}$ the $x$-Pauli
matrix. For this case, to extract a Hamiltonian from the original Lagrangian
formulation, we required that the coupling, $\beta$, goes to infinity to match
the temporal lattice spacing, $e^{-2\beta}\propto a_{\tau}$.

In the above we found the Hamiltonian and Hilbert space in configuration
space; however, one can use a “dual” method as well by expanding the original
Boltzmann weights. Consider the action for the Ising model in Eq. (35) with
the constant dropped,

$\displaystyle
S_{\text{Ising}}=-\beta\sum_{\tau}\sigma_{\tau+1}\sigma_{\tau}.$ (41)

The Boltzmann weight can be expanded in the form,

$\displaystyle\mathbb{T}_{\alpha\alpha^{\prime}}$
$\displaystyle=e^{\beta\sigma^{(\alpha)}_{\tau+1}\sigma^{(\alpha^{\prime})}_{\tau}}$
(42)
$\displaystyle=\cosh(\beta)\sum_{n=0}^{1}(\sigma^{(\alpha)}_{\tau+1}\sigma^{(\alpha^{\prime})}_{\tau})^{n}\tanh^{n}(\beta).$
(43)

which is simply the Euler identity for imaginary angles. Ignoring the factor
out front as it does not affect the Hamiltonian, in this form we can easily do
the summation over the values of $\alpha$ at all lattice sites. The transfer
matrix becomes a diagonal matrix with matrix elements labeled by the integers
$n$, from each Boltzmann factor,

$\displaystyle\mathbb{T}_{nn^{\prime}}=\begin{pmatrix}1&0\\\
0&\tanh(\beta)\end{pmatrix}.$ (44)

In the literature, the integers $n$ are often called dual variables or
characters indices. The terminology is discussed more systematically in Sec.
III.4. To recover a Hamiltonian, the transfer matrix must have the form
$1-a_{\tau}\hat{H}+\ldots$ for small times. This is found by taking
$\beta\rightarrow\infty$, and recalling $\tanh{\beta}=1-2e^{-2\beta}+\ldots$
for large $\beta$. Then the Hamiltonian in these new variables takes the form,

$\displaystyle\hat{H}=h_{z}(1-\hat{\sigma}^{z})$ (45)

with $h_{z}\equiv e^{-2\beta}/a_{\tau}$, and $\hat{\sigma}^{z}$ the $z$-Pauli
matrix. These are two distinct procedures which give Hamiltonians in the time
continuum limit; one, in the original configuration variables, and the second
in the dual variables.

The transfer matrix of the Ising model in higher dimensions can be constructed
in a similar manner Kaufman (1949). The action for the Ising model in
dimension $D$ is the $N=1$ case of Eq. (6). In higher dimensions there are now
two types of interactions. If we consider a particular time slice, we can
first collect all the time links connected to the next time slice, each with a
representation given in Eq. (35). In addition, we have all the spatial links
in the time slice with nearest-neighbor interactions. We can represent any
spin configuration in a time slice as an element of a tensor product of
eigenstates of the Pauli matrix $\sigma^{z}_{\bf x}$:

$\\{{\rm configurations}\\}=\bigotimes_{\bf{x}}\ket{\pm 1}_{\bf x}.$ (46)

We can introduce operators $\hat{\sigma}^{z}_{\bf x}$ and
$\hat{\sigma}^{x}_{\bf x}$ acting on this Hilbert space that can be identified
with a set of qubits. Following Kaufman (1949), we can collect the two types
of interactions in two matrices. Using the relation

$\displaystyle
e^{\beta\sigma_{\alpha^{\prime}}\sigma_{\alpha}}=(2\sinh(2\beta))^{1/2}(\exp(\tilde{\beta}\sigma^{x}))_{\alpha^{\prime}\alpha},$
(47)

where $\tilde{\beta}$ is the dual inverse temperature introduced by Kramers
and Wannier Kramers and Wannier (1941) which satisfies the relation
$\tanh(\tilde{\beta})=e^{-2\beta}$, the contribution of time links can be
summarized with the matrix connecting the configuration of the two time slices

${\mathbb{V}}_{1}=(2/\sinh(2\tilde{\beta}))^{N_{s}^{D-1}/2}e^{\tilde{\beta}\sum_{\bf
x}\hat{\sigma}^{x}_{\bf x}}.$ (48)

On the other hand, the spatial links can be recast in the diagonal matrix

${\mathbb{V}}_{2}=e^{\beta\sum_{{\bf x},j}\hat{\sigma}^{z}_{\bf
x}\hat{\sigma}^{z}_{{\bf x}+\hat{j}}}$ (49)

with $\hat{j}$ a unit vector pointing along one of the $D-1$ directions on the
lattice. We can now write the transfer matrix as

$\mathbb{T}={\mathbb{V}}_{2}^{1/2}{\mathbb{V}}_{1}{\mathbb{V}}_{2}^{1/2},$
(50)

where the matrix indices label the spatial configurations. Geometrically, the
Hilbert space is located on the time slices. Illustrations of the time-slices
and the location of the Hilbert space can be found in Sec. XII.

Alternatively, we can work in a dual representation where the $\sigma^{x}_{\bf
x}$ are diagonal

$\tilde{\mathbb{T}}={\mathbb{V}}_{1}^{1/2}{\mathbb{V}}_{2}{\mathbb{V}}_{1}^{1/2},$
(51)

where the matrix indices label sets of group characters. Geometrically, the
Hilbert space is located between the time slices. Graphical illustrations of
this situation will be provided in Sec. XII. The construction generalizes
easily for finite Abelian groups and in a non-trivial way for continuous and
compact Abelian groups. The advantage of using the second (dual)
representation is that it remains discrete and as we will see in Sec. XIII.1,
it preserves the symmetry when truncations are applied Meurice (2019, 2020a).

###  III.4 Topological and geometrical dualities

In Sec. III.3 we used the concept of duality in two occasions. The first was
the relation between $\beta$ and $\tilde{\beta}$ which interchanges their low
and large values regimes. The second occasion was the discussion of the two
ways to represent the transfer matrix. In addition, in Secs. VII.3 and IX.1.2,
we will use the geometrical duality reviewed in Savit (1980). Duality is a
general concept used in many branches of mathematics. According to Atiyah
Atiyah (2007) duality gives “two different points of view of looking at the
same object”. In the following we clarify the various usages of the concepts
in the rest of this review.

An important notion of duality is the so-called “Pontryagin duality”
Pontryagin (1939) coming in the study of topological groups. It relates an
Abelian group and its characters (for instance, Fourier modes). It states that
if the former is compact, the later is discrete and vice-versa. The simplest
situation is a finite group which is compact and discrete. In the case of
finite cyclic Abelian groups the characters form a finite group which is
isomorphic to the group itself Serre (1973). A simple example is
${\mathbb{Z}}_{q}\ $the additive group of integers modulo $q$. If $x$ denotes
an element of $\mathbb{Z}_{q}$, the characters (See A.1) have the form

$\chi_{k}(x)=\exp\left(i\frac{2\pi}{q}kx\right),$ (52)

and clearly satisfy the character property

$\chi_{k}(x+x^{\prime})=\chi_{k}(x)\chi_{k}(x^{\prime}).$ (53)

The product of two characters is another character

$\chi_{k}(x)\chi_{k^{\prime}}(x)=\chi_{k+k^{\prime}}(x),$ (54)

and one sees that they also form a $\mathbb{Z}_{q}$ group. They obey the
orthogonality relations

$\frac{1}{q}\sum_{x=0}^{q-1}\chi_{k}(x)\chi^{\star}_{k^{\prime}}(x)=\delta_{k,k^{\prime}},$
(55)

and

$\frac{1}{q}\sum_{k=0}^{q-1}\chi_{k}(x)\chi^{\star}_{k}(x^{\prime})=\delta_{x,x^{\prime}}.$
(56)

A simple example of a continuous (nondiscrete) and compact group is $U(1)$,
the multiplicative group of complex numbers with modulus one. The group
properties can be reformulated in an additive manner by introducing the phases
$z=e^{i\varphi}$ which are added modulo $2\pi$. Topologically it is a circle
which is a compact manifold. Its characters are discrete and labelled by the
integers. They are the usual Fourier modes $e^{in\varphi}$. The orthogonality
relations appear in an asymmetric way

$\int_{-\pi}^{\pi}\frac{d\varphi}{2\pi}e^{in\varphi}(e^{in\varphi})^{\star}=\delta_{n,n^{\prime}}$
(57)

$\sum_{n=-\infty}^{\infty}e^{in\varphi}(e^{in\varphi^{\prime}})^{\star}=2\pi\sum_{m=-\infty}^{\infty}\delta(\varphi-\varphi^{\prime}+2\pi
m)$ (58)

In practice, a deep understanding of the mathematical statements made above is
not necessary and we just need to remember a few character expansions. For the
Ising models we have functions over the multiplicative group $\sigma=\pm 1$,
and we can recall the expansion from Sec. III.3 as a character expansion

$\displaystyle\bra{\sigma_{\alpha^{\prime}}}\mathbb{T}\ket{\sigma_{\alpha}}$
$\displaystyle=\exp(\beta\sigma_{\alpha^{\prime}}\sigma_{\alpha})$
$\displaystyle=\frac{1}{2}\sum_{n=0}^{1}\lambda_{n}(\beta)(\sigma_{\alpha^{\prime}}\sigma_{\alpha})^{n}$
(59)
$\displaystyle=\cosh(\beta)+\sigma_{\alpha^{\prime}}\sigma_{\alpha}\sinh(\beta).$
(60)

A similar representation can be obtained for the ${\mathbb{Z}}_{q}\ $spin
models. If we replace the discrete angle variables $\frac{2\pi}{q}x$ by
continuous ones $\varphi$, we obtain a “matrix” with continuous elements which
can be calculated using the standard Fourier transform

$\displaystyle\bra{\varphi^{\prime}}\mathbb{T}\ket{\varphi}$
$\displaystyle=\exp(\beta\cos(\varphi^{\prime}-\varphi))$ (61)
$\displaystyle=\sum_{n=-\infty}^{+\infty}I_{n}(\beta)\exp(in(\varphi^{\prime}-\varphi)),$
(62)

where $I_{n}(\beta)$ is the modified Bessel function of order $n$. We will
show in Sec. VII.2 that finite versions of this expansion hold for the finite
${\mathbb{Z}}_{q}\ $subgroups.

Generalizations of Pontryagin duality to compact non-Abelian groups appear in
the Peter-Weyl theorem Peter and Weyl (1927). As an example, this translates
into expansions in spherical harmonics for problems involving the O(3)
symmetry. Related methods were used for tensor networks with continuous
symmetries Tagliacozzo  _et al._ (2014); Zohar and Burrello (2015).

From the point of view of quantum computing, we see that using compact fields
guarantees that we can replace the continuous integrals by discrete sums. Our
strategy is to associate the indices of these sums with quantum states. One
important aspect of TLFT is that when the fields are compact, we don’t need to
discretize the integrals using numerical approximations. Instead, we can use
characters expansions as Eq. (62) where some integrals have been done exactly
and result in Bessel functions that we can input with any desired accuracy. In
other words, the difficult part of the classical path-integral approach can be
done efficiently with classical methods. After that, the original integrals
reduce to orthogonality relations and can be performed exactly.

Another notion of duality is of a geometrical nature. It is related to the
Levi-Civita tensor $\epsilon^{\mu_{1}\dots\mu_{D}}$. Its meaning is dimension-
dependent and relates objects of dimension $d$ to objects of dimension $D-d$.
For instance, in $D=4$, a dual field-strength tensor with two indices is
obtained by contracting the original field-strength tensor with the Levi-
Civita tensor. This duality transformation interchanges the electric and
magnetic fields and reduces to the identity when repeated twice. It also
relates sites to four-dimensional hypercubes, and plaquettes to plaquettes. In
$D=3$ it relates the field-strength tensor to a divergenceless pseudovector
$\epsilon^{ijk}\partial_{j}A_{k}$, links to plaquettes and sites to cubes.

The various notions of duality are often used simultaneously. A common example
is the “dual formulation” of the Ising model. As we will explain in Sec. V.2,
the new set of indices from Eq. (60) which is a consequence of Pontryagin
duality, leads to a representation in terms of paths. We can then try to
represent these paths as the boundaries of surfaces, which brings the
geometrical duality and new “dual variables” together Savit (1980).

There are also occasions where the phrase “dual variables” has become
associated with generic integer fields, regardless or their origin or their
relation to the concepts of duality mentioned above Bruckmann  _et al._
(2016); Gattringer  _et al._ (2018b, 2015a); Marchis and Gattringer (2018);
Bruckmann  _et al._ (2015). These integer fields can arise from a Taylor
series of the Boltzmann weight, _e.g._ in the case of the Ising model,

$\displaystyle
e^{\beta\sigma\sigma^{\prime}}=\sum_{q=0}^{\infty}\frac{\beta^{q}}{q!}(\sigma\sigma^{\prime})^{q}.$
(63)

This expansion associates a natural number with the links of the lattice,
similar to the character expansion from before. Also just as before, this
expansion allows one to perform the path-integral sums over the $\sigma$
fields, leaving one with a theory of constrained, positive, integer fields on
the links of the lattice. However, instead of two values, _i.e._ $n=0,1$,
these integers can take on an infinite number of values. So while this
expansion accomplishes similar feats, it can be seen as a less economical
parametrization of the model. Of course, the character expansion from before
can be found within the Taylor expansion by summing the even and odd integers,
respectively, leaving one with two terms (0 and 1).

###  III.5 Real-time evolution with qubits

We are now ready to discuss the possibility of using quantum devices to
represent states of the Hilbert space emerging from the transfer matrix
construction and to design methods to apply unitary transformations
corresponding to the real time evolution. For the sake of definiteness, we
assume that we have at our disposal a set of qubits. For the spin Ising model,
the construction of the transfer matrix leads to a Hilbert space with an
obvious qubit structure given by Eq. (46). The method to take the time
continuum limit and identify a Hamiltonian using Eq. (33) is well-known
Fradkin and Susskind (1978); Kogut (1979): we deform the original transfer
matrix by increasing $\beta$ in the time direction (which makes the dual value
$\tilde{\beta}$ small) and decreasing $\beta$ in the spatial directions. The
arguments of the exponentials in ${\mathbb{V}}_{1}$ and ${\mathbb{V}}_{2}$
become infinitesimal and provide the two non-commuting pieces of the “quantum
Hamiltonian”. The role played by ${\mathbb{V}}_{1}$ is quite special in the
intuitive picture that we are trying to convey: it only acts on “single
qubits” without connecting them together. Consequently, working in the
representation where $\hat{\sigma}^{x}_{\bf x}$ is diagonal is a good starting
point. Next, we can “turn on” the terms in ${\mathbb{V}}_{2}$. At lowest order
in the time lattice spacing they only connect qubits which are nearest
neighbors. We will show in Sec. III.6 that this type of situation permits
quantum computation in a time scaling polynomially with the size of the system
Lloyd (1996).

Character expansions and TLFT provide the natural tools to perform similar
constructions for the models presented in Sec. II.2. A first step consists in
isolating “building blocks” which are localized in space and have a very
simple real-time evolution. The models have interactions associated with links
and plaquettes. As a first approximation we set the interactions on spatial
links and space-space plaquettes to zero.

For spin models, this results in a collection of $N_{s}^{D-1}$ isolated one-
dimensional spin models. These isolated models are the building blocks. They
are solvable and it is easy to calculate their evolution at real time. For
gauge models we have, in addition, electric degrees of freedom which can be
associated with the spatial links of a given time slice of the lattice. They
are required to satisfy a constraint called “Gauss’s law”, however, when this
condition is satisfied, the real-time evolution in the isolation limit
described above is straightforward as we will explain in Sec. XII.3. The
Hilbert space of the isolated building blocks depends on the model considered.
For the Ising spin models, a single qubit is all we need. For models with
continuous symmetries, the exact treatment requires an infinite Hilbert space,
however, small size truncations provide good approximations and preserve the
symmetries of the models Meurice (2019, 2020a). This is a very attractive
feature of TLFT which is discussed in Sec. XIII.1. Having setup a finite
Hilbert space with isolated building blocks, our next step is to restore the
interactions associated with the spatial links and the space-space plaquettes.
This will be done for a variety of models in Sec. XII. Independently of the
model-specific aspect of this procedure, it is clear that each building block
is only connected to a limited numbers of other building blocks. For instance,
for a spin model, there are two connections for each spatial direction. This
quasilocality is crucial to implement real-time evolution with a quantum
computer.

###  III.6 Lloyd-Suzuki-Trotter product formula

An important motivation for using a quantum computer is to calculate the real-
time evolution for systems with many degrees of freedom having quasilocal
interactions in the sense discussed in Sec. III.5. Ideally the time to perform
computations should scale polynomially with the size of the system rather than
exponentially. A general argument leading to these conclusions has been put
forward by S. Lloyd in Lloyd (1996) where he states that: “Feynman’s 1982
conjecture, that quantum computers can be programmed to simulate any local
quantum system, is shown to be correct”. The proof is based on the basic idea
behind the Suzuki-Trotter product formula Trotter (1959); Suzuki (1976); Reed
and Simon (1980), namely that for two non-commuting operators $\hat{A}$ and
$\hat{B}$ and sufficiently small $\epsilon$

$e^{i\epsilon(\hat{A}+\hat{B})}\simeq
e^{i\epsilon\hat{A}}e^{i\epsilon\hat{B}}+\mathcal{O}(\epsilon^{2}).$ (64)

In the standard construction of the path integral in quantum mechanics, it is
applied to the kinetic and potential energy, but it can also be applied to all
the quasilocal parts of the Hamiltonian.

The argument goes as follows Lloyd (1996). Consider a system composed of $N$
variables with Hamiltonian

$\hat{H}=\sum_{j=1}^{\ell}\hat{H}_{j}$ (65)

where each $\hat{H}_{j}$ acts on a space that involves at most $k_{max}$ of
the variables. It is assumed that $\ell$ increases linearly with $N$ but that
$k_{max}$ is fixed by the dimension and the symmetries and independent of $N$.
The individual $H_{j}$ can be represented as finite matrices in their local
subspace. Under these assumptions, it is shown Lloyd (1996) that the error
associated with the approximation

$e^{i\hat{H}t}\simeq(e^{i\hat{H}_{1}t/n}\dots\
e^{i\hat{H}_{\ell}t/n})^{n}+\dots$ (66)

can be controlled by taking $n$ large enough. In addition, once the accuracy
goal is determined, the computing time scales linearly in $N$ and $t$.

In order to fix the ideas, for most of the available quantum computers, the
$\hat{H}_{j}$ act on one or two qubits and can be represented by $2\times 2$
or $4\times 4$ matrices in this restricted space. A simple quantum circuit for
the quantum Ising model used in Refs. Gustafson  _et al._ (2019, 2021a) is
displayed in Fig. 2. The basic elements are rotations generated by
$\hat{\sigma}_{j}^{x}$ or $\hat{\sigma}_{j}^{z}$ and acting on the $j$-th
qubit only and CNOT gates acting on a pair of qubits and flipping the target
qubit when the control qubit is in the $\ket{1}$ state. The circuit can be
repeated in the space (vertical) and time (horizontal) directions and conveys
the linear scalings stated above. More evolved circuits can be designed with
the purpose of creating energy eigenstates Verstraete  _et al._ (2009);
Cervera-Lierta (2018).

It is instructive to compare the computational resources to perform the
unitary rotation $\exp(i\theta\hat{\sigma}^{x}_{j}\hat{\sigma}^{x}_{k})$ for
which the operator in the exponential flips the $j$ and $k$-th qubit and does
not act on the rest of the system. With a quantum computer, one expects that
this operation keeps the coherence of the qubits with the rest of the system,
present in the original state, but the cost is ideally independent on the size
of the system. On the other hand, the same operation on an arbitrary state
performed with a classical computer involves a matrix with $2^{N}$ non-zero
and non-diagonal matrix elements and the resources necessary increase
exponentially with the size of the system.

Figure 2: Circuit for 4 qubits with open boundary conditions used in Gustafson
_et al._ (2019). Here $RZ$ and $RX$ are rotations about the $z$ and $x$-axes,
respectively.

###  III.7 Dealing with noise in the NISQ era

The discussion in Sec. III.6 considers algorithmic aspects in an idealized
situation where computer errors or noise can be neglected. Quantum computing
technology is still in its early development. The current Noisy Intermediate
Scale Quantum (NISQ) computer hardware platforms can only accommodate small-
depth circuits and sources of errors need to be understood in detail Preskill
(2018). Understanding these errors and their noise, adjusting, and contriving
algorithms to lessen these errors is the topic of error mitigation in quantum
computing.

Various types of noise affect the single-qubit gates in a way that can be
parameterized Nielsen and Chuang (2000) in terms of the density matrix
$\hat{\rho}$:

$\displaystyle\mathcal{E}(\hat{\rho};p_{x},p_{y},p_{z})=(1-p)\hat{\rho}$
$\displaystyle+p_{x}\hat{\sigma}^{x}\hat{\rho}\hat{\sigma}^{x}$
$\displaystyle+p_{y}\hat{\sigma}^{y}\hat{\rho}\hat{\sigma}^{y}+p_{z}\hat{\sigma}^{z}\hat{\rho}\hat{\sigma}^{z}.$
(67)

Here $\mathcal{E}$ is a quantum operation, and the values $p_{x}$, $p_{y}$,
and $p_{z}$ correspond to the probabilities of a $\sigma^{x}$, $\sigma^{y}$,
or $\sigma^{z}$ error, respectively, occurring and $p=p_{x}+p_{y}+p_{z}$. The
error channel for two qubit gates is given by
$\mathcal{E}^{(2)}=\mathcal{E}\bigotimes\mathcal{E}$. Modeling the error in
this way neglects spatial and temporal correlations, 2-qubit correlations, and
assumes that the errors are identical and evenly distributed. In practice, if
a classical simulation of the qubit evolution is performed, each unitary
evolution operation needs to be followed by applying one of the four
possibilities $(\hat{1},\hat{\sigma}^{x},\hat{\sigma}^{y},\hat{\sigma}^{z})$
with respective probabilities $1-p$, $p_{x}$, $p_{y}$, and $p_{z}$ on the
qubits involved. The probability distribution of errors is only manifest after
averaging over an ensemble of runs.

In addition, readout errors (misidentifying a $|1\rangle$ for a $|0\rangle$ or
vice-versa) should be taken into account for all current quantum computing
platforms. Given estimates of the probabilities for these errors, it is
possible to correct the actual measurements by a multiplicative factor Kandala
_et al._ (2018).

A common NISQ strategy for error mitigation is to increase the source of error
in a controllable way and then extrapolate to the limit where the error is not
present. Examples with superconducting qubits are given in Refs. Temme  _et
al._ (2017); Kandala  _et al._ (2019); Klco  _et al._ (2018); Gustafson  _et
al._ (2019, 2021a). A simple way to increase the error is to insert two
successive CNOT gates. Their exact multiplication is the identity, however for
a NISQ devices it increases the chances of errors. Note that mitigation
methods can also be used for quantum variational methods Li and Benjamin
(2017).

In this context, the choice of Trotter step $\delta t$ is crucial because the
number of steps is limited by the loss of coherence and the noise. A concrete
discussion can be found in Gustafson  _et al._ (2021a) for the quantum Ising
model in one spacial dimension with four sites. If we pick a small $\delta t$
with a good control on the $(\delta t)^{2}$ error, we may not be able to reach
a time scale relevant for what we want to learn. However, it appears that by
picking a $\delta t$ significantly larger, the rigorous bound is not sharp and
the empirical bound is much tighter as shown in Fig. 3 where we plot the
operator norm of the Trotter error

$\displaystyle\Delta_{2}U$ $\displaystyle\equiv$ $\displaystyle
e^{-i(h_{T}\hat{H}_{T}+J\hat{H}_{NN})\delta t}-e^{-ih_{T}\hat{H}_{T}\delta
t}e^{-iJ\hat{H}_{NN}\delta t}$ (68) $\displaystyle\simeq$
$\displaystyle\frac{h_{T}J}{2}[\hat{H}_{T},\hat{H}_{NN}](\delta t)^{2},$ (69)

with

$\hat{H}_{T}=-\sum_{j=1}^{4}\hat{\sigma}^{z}_{j},\ {\rm and}\
\hat{H}_{NN}=-\sum_{j=1}^{3}\hat{\sigma}^{x}_{j}\hat{\sigma}^{x}_{j+1}$ (70)

The nonlinear aspects of the error are quite interesting and may involve
resonances Gustafson  _et al._ (2021a). For other recent developments see
Sewell and Jordan (2021); A Rahman  _et al._ (2021); Mishra  _et al._ (2019);
Gustafson (2020); Gustafson  _et al._ (2021c, b); Gustafson and Lamm (2021);
Funcke  _et al._ (2021); Honda  _et al._ (2021); Gustafson  _et al._ (2021a).

Figure 3:  $||\Delta_{2}U||$ versus $\delta t$ for the quantum Ising model
with $h_{T}=1$ and $J=0.02$. The vertical lines are at $\pi/2$, $3\pi/4$ ….
For details see Ref . Meurice (2020b)

###  III.8 Quantum computations and simulations

The idea of using quantum devices to mimic or study theoretical quantum models
has been a subject of intense activity in recent years. Here we mention
theoretical and experimental studies of the Ising model, and lattice gauge
theories, in the context of quantum computation and quantum simulation.

####  III.8.1 Ising model

As stressed before, the 1+1 dimensional Ising model is the quintessential
model to use to check computational tools and proposed algorithms. This is
because the model possesses non-trivial features, but has been solved exactly.
It is not surprising then that this is one of the first models to be tried and
tested in different forms of quantum computation. There have been various
investigations into how to simulate the model on a quantum computer, and what
interesting observables to measure (_e.g._ phase shift, and thermodynamics)
Lidar and Biham (1997); Mostame and Schützhold (2008); Gustafson  _et al._
(2019); You  _et al._ (2013). The following approaches have been considered:
super-conducting qubit machines, trapped-ion machines, cold atoms trapped in
optical lattices, and Rydberg atom simulators.

On the super-conducting qubit front, Refs. Hebenstreit  _et al._ (2017); Zhang
_et al._ (2017) made initial simulations, and in some cases at relatively
large system sizes. Their calculations of various spin observables matched the
corresponding quantities in exact diagonalization well. In Cervera-Lierta
(2018) a simulation of the Ising model using a few spins was also carried out,
and a comparison between theory and computation was made for the average
magnetization. In Refs. Gustafson  _et al._ (2021a); Smith  _et al._ (2019),
the Ising model was simulated using a few qubits on IMB’s machines. It was
shown that the Richardson extrapolation could be used to mitigate the noise in
the regime where the nonlinear effects are not too large Gustafson  _et al._
(2021a). In addition, in Chen  _et al._ (2019) a quantum-classical approach
was taken which uses a variational algorithm to compute the ground-state wave
function; the so-called Variational Quantum Computing, where the quantum
computer prepares a state with a circuit ansatz depending on some number of
parameters. The expectation value of the energy is computed from this state,
and circuit parameters are then tuned classically based on the quantum
expectation value in order to minimize the energy.

Using trapped ions, the first work was Friedenauer  _et al._ (2008), followed
by extensive investigations from student’s work Friedenauer (2010); Korenblit
(2013). These are pioneering investigations into the trapped-ion platform.
Refs. Kim  _et al._ (2010); Edwards  _et al._ (2010); Islam  _et al._ (2011);
Kim  _et al._ (2011) look at the phase structure of the model by simulating a
few spins. They calculate the phase diagram of the model using the probability
of the state to be in a ferromagnetic state, or moments of the magnetization.
There has also been an investigation into measuring the Rényi entropy using
digital quantum gates. In Linke  _et al._ (2018) the authors consider a two-
site antiferromagnetic Ising model, and using the SWAP gate measure the parity
of two copies of the system Johri  _et al._ (2017). Ref. Islam  _et al._
(2015) used bosonic many-body states and interfered the copies to extract sub-
system parities. They then calculate the Rényi entropy from sub-system
parities.

Finally, another promising approach is to use highly excited (Rydberg) states
of atoms, which allows for strong atom-atom interactions across relatively
large distances. In Kim  _et al._ (2017) the group uses a chain of 19 Rubidium
atoms whose interactions they control through tuning lattice parameters to
simulate the model. In Guardado-Sanchez  _et al._ (2018) out-of-equilibrium
dynamics are explored through a quench using an array of Lithium atoms, again
placed in a Rydberg state. Reference Simon  _et al._ (2011) uses trapped
Rubidium atoms to simulate the anti-ferromagnetic Ising model. They are able
to identify a phase transition between para- and anti-ferromagnetic phases,
and observe magnetic domains using a site-resolved atomic microscope and noise
correlations measurements. In Bernien  _et al._ (2017) the authors demonstrate
the use of configurable and programmable arrays of atoms, and simulate an
Ising-like model on 51 qubits. They observe a phase transition between
symmetric and ordered phases and discuss out-of-equilibrium properties of spin
models Keesling  _et al._ (2019). For reviews of this topic see Schauss
(2018); Browaeys and Lahaye (2020).

####  III.8.2 Gauge theories

The use of optical lattices Bloch  _et al._ (2008) to quantum simulate lattice
gauge theories has been developed extensively. Proposals for the quantum
simulation of lattice gauge theories beginning with early work on Abelian
models includes Zohar and Reznik (2011); Zohar  _et al._ (2012); Banerjee  _et
al._ (2012); Tagliacozzo  _et al._ (2013b); Zohar  _et al._ (2013b), and for
digital quantum devices in Zohar  _et al._ (2017b, a). In the case of non-
Abelian models: Banerjee  _et al._ (2013); Zohar  _et al._ (2013a);
Tagliacozzo  _et al._ (2013a). For reviews see Wiese (2013); Zohar  _et al._
(2015b); Bañuls and Cichy (2020); Bañuls  _et al._ (2020). A useful early
reference on quantum computing for non-Abelian gauge theories is Byrnes and
Yamamoto (2006). For recent developments combining condensed matter and gauge
theory ideas see Li  _et al._ (2016); Kasper  _et al._ (2017); Clark  _et al._
(2018); Schweizer  _et al._ (2019). Trapped ions Leibfried  _et al._ (2003)
have provided new opportunities to approach lattice gauge theory models
Davoudi  _et al._ (2020a). Rydberg atoms offer a versatile platform for gauge
theories Celi  _et al._ (2019); Surace  _et al._ (2020); Meurice (2021).

The Schwinger model is often the first target to develop new approaches
Martinez  _et al._ (2016); Kharzeev and Kikuchi (2020); Davoudi  _et al._
(2020a); Kasper  _et al._ (2017); Klco  _et al._ (2018); Surace  _et al._
(2020); Magnifico  _et al._ (2019). For recent work on non-Abelian models see
Silvi  _et al._ (2019); Raychowdhury and Stryker (2020b, a); Davoudi  _et al._
(2020b); Dasgupta and Raychowdhury (2020).

##  IV Quantum versus classical

In this section, we discuss the meaning of “classical” and “quantum” for
models, phase transitions and tensor networks.

###  IV.1 Models

In textbooks on quantum mechanics, it is a common procedure to start with the
Hamiltonian formulation and derive a path-integral representation which can be
extended to field theory. The path-integral formalism allows for formulations
that are manifestly gauge-invariant and treats space and time on equal
footing. For these reasons, it can be argued that the fundamental definition
of relativistic models should be done in terms of the action and the measure
of integration over all the possible configurations without reference to a
Hamiltonian in the first place. Examples of such actions are given in Sec.
II.2.

It is very common to call models formulated with the path-integral “classical”
while the corresponding formulation using a Hamiltonian acting on a Hilbert
space is called “quantum”. However, except for possible discretization
artifacts, the two formulations describe the same quantum behavior. In the
path-integral formulation for bosonic fields, the action is calculated in
terms of c-numbers as in the classical formulation but the sum over all the
configurations provides a quantum description which amounts to more than the
classical equations of motion. In other words, the path-integral is an
alternate method of quantization which is very convenient in Euclidean time.

Starting with a classical action, a Hamiltonian can be constructed from the
action by using the transfer matrix formalism. This was briefly demonstrated
in Sec. III.3, and will be discussed in detail in Sec. XII for the models
introduced in Sec. II.2. A discussion with detailed references on the
connection between statistical mechanics in $D$ dimensions and quantum
Hamiltonians in $D-1$ dimensions, can be found in the classic work of Wilson
and Kogut Wilson and Kogut (1974).

###  IV.2 Phase transitions

The actions for spin and gauge models introduced in Sec. II.2 contain
parameters $\beta$, or $\beta_{pl.}$ for pure gauge theories, that are often
called the “inverse classical temperature”, or “coupling constant”, and can be
associated with “classical phase transitions”. For instance the $D=2$
classical Ising model has a spontaneous magnetization when
$\beta>\beta_{c}=(1/2)\ln(1+\sqrt{2})$.

In contrast, given a Hamiltonian $\hat{H}$, we can define a thermal quantum
partition function with temperature $T_{qu.}$ in the usual way

$Z(T_{qu.})=\text{Tr}[e^{-\hat{H}/T_{qu.}}],$ (71)

where $T_{qu.}$ has in general a different meaning than $1/\beta$ in the
classical formulation. With the lattice formulation at Euclidean time as a
starting point, we have the identification

$1/T_{qu.}=N_{\tau}a_{\tau},$ (72)

and the non-zero temperature is associated with the finite extent of the
temporal dimension. A typical situation of interest is to start in an ordered
phase at $T_{qu.}=0$, corresponding to the infinite Euclidean time limit, and
induce a finite-temperature phase transition into a disordered phase by taking
a sufficiently small temporal extent. In this way, the temporal extent of the
lattice is responsible for a finite-temperature phase transition, as opposed
to the coupling $\beta$, which is unrelated to the temperature in the quantum
partition function. Transitions related to the couplings at $T_{qu}=0$ are
quantum phase transitions. Sometimes, the transition can be understood in
terms of the classical phase diagram in $D-1$ dimensions. A more detailed
discussion can be found in Cardy’s monograph Cardy (1996)

###  IV.3 Tensor networks

In Secs. V to XI we will introduce “classical tensors” in order to reformulate
classical models as defined in Sec. IV.1. The partition function of these
models can be visualized as an assembly obtained by “wiring” (tracing)
together objects carrying multiple “legs” (tensor indices) and attached to the
sites, links or plaquettes of a Euclidean space-time lattice. This type of
classical constructions has been inspired Levin and Nave (2007); Gu and Wen
(2009) by various quantum tensor networks Fannes  _et al._ (1992); Vidal
(2003, 2004); Verstraete and Cirac (2004a); Shi  _et al._ (2006); Perez-Garcia
_et al._ (2006); Vidal (2007); Verstraete  _et al._ (2008), to mention a few
references, developed in various contexts often related to the density matrix
renormalization group (DMRG) method White (1992). There is abundant literature
on the subject reviewed for instance in Schollwock (2005); Verstraete  _et
al._ (2008); Cirac and Verstraete (2009); Schollwock (2011); Orus (2014);
Silvi  _et al._ (2017); Montangero (2018); Ran  _et al._ (2020).

One important idea is the representation of quantum states by matrix product
states (MPS) which appear in several references mentioned above. As an
example, for a one-dimensional quantum chain problem with $N_{s}$ sites, an
arbitrary element of the Hilbert space can be written as

$\ket{\psi}=\sum_{i_{1},\dots,i_{N_{s}}}c_{i_{1},\dots,i_{N_{s}}}\ket{i_{1},\dots,i_{N_{s}}},$
(73)

where

$\ket{i_{1},\dots,i_{N_{s}}}=\ket{i_{1}}\otimes\dots\otimes\ket{i_{N_{s}}},$
(74)

and each of the indices runs over a local Hilbert space of dimension $d_{H}$
attached to a site. The dimension of the Hilbert space is $d_{H}^{N_{s}}$. It
represents an exponential growth with the size of the system that rapidly
becomes computationally unmanageable. For a MPS state, one assumes the form

$c_{i_{1},\dots,i_{N_{s}}}=\text{Tr}[A_{i_{1}}\dots A_{i_{N_{s}}}],$ (75)

where the $A_{i_{j}}$ are $d_{B}\times d_{B}$ matrices for each value of
$i_{j}$. $d_{B}$ is called the bond dimension. A graphical representation of
such a state is shown in Fig. 4 for open boundary conditions. The filled
circles represent the matrices, the vertical lines represent open indices with
$d_{H}$ values and the horizontal lines traced indices with $d_{B}$ values.

Figure 4: Illustration of states in the MPS approach.

The size of the MPS “subset” 111The sum of two MPS is not a MPS, so we will
not call it a “subspace”. is only growing like $N_{s}\times d_{B}^{2}\times
d_{H}$, so linearly with the size of the system. Similarly, one can represent
operators using the trace of $N_{s}$ $d_{B}\times d_{B}$ matrices
$A_{i_{j}}^{i^{\prime}_{j}}$ with two indices in the one-site Hilbert space at
a computational cost scaling like $N_{s}\times d_{B}^{2}\times d_{H}^{2}$.
This is illustrated in Fig. 5.

Figure 5: Illustration of operators in the MPS approach.

Objects with similar shapes will appear as “times slices” of the classical
construction. In the case where the indices take an infinite number of values
and the sums are truncated in a way that the number of indices kept in the
time direction (similar to $d_{H}$), is the same as the number in the space
directions (similar to $d_{B}$), we will use the isotropic notation
$D_{\text{cut}}$ (one can of course leave the number of states in the temporal
and spatial directions independent of each other). On the other hand, we will
discuss the anisotropic situation where the time direction is singled out to
define the transfer matrix and the Hamiltonian in Sec. XII.

MPS and the DMRG have been used in studies of 1+1 dimensional models such as
the Bose-Hubbard model Bonnes  _et al._ (2014), the Schwinger model Byrnes
_et al._ (2002); Bañuls  _et al._ (2013); Buyens  _et al._ (2014, 2016), see
Sec. XI.7 for more information, SU(2) gauge theory Bañuls  _et al._ (2017b);
Kühn  _et al._ (2015), and the O(3) nonlinear sigma model Bruckmann  _et al._
(2019). Fermionic tensor network studies, and the Hubbard model were discussed
in Barthel  _et al._ (2009); Corboz  _et al._ (2010). Tensor network
techniques for lattice gauge theories are also discussed in Tagliacozzo  _et
al._ (2014); Silvi  _et al._ (2016, 2019); Rico  _et al._ (2014); Pichler  _et
al._ (2016); Zohar and Burrello (2015) and reviewed in Bañuls  _et al._
(2018); Bañuls and Cichy (2020). For a review of matrix product operators and
their relations with the transfer matrix see Haegeman and Verstraete (2017).
The reader is directed to Sec. XII.7 for references on tensor network studies
in 2+1, and 3+1 dimensions using generalizations of MPS.

The MPS framework can also be used to perform real-time calculations based on
the Suzuki-Trotter approximations. An example of method is called the time-
evolving block decimation reviewed in Paeckel  _et al._ (2019).

The success and limitations of the MPS approach can be analyzed in terms of
“area laws” Plenio  _et al._ (2005); Verstraete  _et al._ (2006). Following a
short pedagogical discussion Schollwock (2011), if a bipartition of a system
is introduced, one expects that the entanglement entropy between the two parts
scales like the size of their boundary with possible logarithmic corrections.
In one spatial dimension, we can separate a MPS in two parts by cutting a
single bound carrying a maximal entropy $\ln_{2}d_{B}$. Consequently, $d_{B}$
only needs to increase like the size of the system to capture a possibly
logarithmic entanglement. On the other hand, if we fill a two-dimensional
surface with a one-dimensional MPS, the entanglement grows at least like the
linear size of the system which forces $d_{b}$ to grow exponentially with this
size. Projected Entangled Pair states (PEPS) were proposed Verstraete and
Cirac (2004b) to try to overcome this difficulty.

It should also be mentioned that tensor methods have been applied Nishino and
Okunishi (1996) to the transfer matrix treatment of classical statistical
model such as the problem of monomers and dimers on a rectangular lattice
Baxter (1968) where variational methods can be applied and compared to the
DMRG approach. This approach is called the corner transfer matrix
renormalization group method. Pioneering work connecting the DMRG and the
transfer martix of the Ising model can be found in Nishino and Okunishi
(1996). Related results and their connections with MPS are reviewed in
Haegeman and Verstraete (2017).

For a very recent and comprehensive review on MPS and related topics, we
recommend Cirac  _et al._ (2020).

##  V Tensor methods explained with the Ising model

###  V.1 Tensor formulation

In this section we will construct a tensor formulation for the Ising model in
$D$ dimensions. The partition function for the Ising model is

$Z_{\text{Ising}}=\prod_{x}\sum_{\sigma_{x}=\pm
1}e^{\beta\sum_{x,\mu}\sigma_{x+\hat{\mu}}\sigma_{x}}.$ (76)

For each link $(x,\mu)$ we use the expansion

$\displaystyle e^{\beta\sigma_{x+\hat{\mu}}\sigma_{x}}=$
$\displaystyle\cosh(\beta)\sum\limits_{n_{x,\mu}=0,1}[\sigma_{x+\hat{\mu}}\sqrt{\tanh(\beta)}\sigma_{x}\sqrt{\tanh(\beta)}\
]^{n_{x,\mu}},$ (77)

This attaches an index $n_{x,\mu}$ at each link $(x,\mu)$. It is then possible
to pull together the various factors of
$(\sqrt{\tanh(\beta)}\sigma_{x})^{n_{x,\mu}}$ from links coming from a single
site $x$, and sum over $\sigma_{x}$,

$\displaystyle\sum_{\sigma_{x}=\pm
1}\prod_{\mu=1}^{D}(\sqrt{\tanh(\beta)}\sigma_{x})^{n_{x-\hat{\mu},\mu}+n_{x,\mu}}$
$\displaystyle=(\sqrt{\tanh(\beta)})^{\sum_{\mu}n_{x-\hat{\mu},\mu}+n_{x,\mu}}$
$\displaystyle\times 2\delta({\rm
mod}[\sum_{\mu}n_{x-\hat{\mu},\mu}+n_{x,\mu},2])$ (78)

using the relation

$\sum_{\sigma_{x}=\pm 1}\sigma_{x}^{n}=2\delta({\rm mod}[n,2]).$ (79)

at every site on the lattice. The expression $\delta({\rm mod}[n,2])$ is 1
when $n$ is even (0 modulo 2) and 0 otherwise. We can rewrite the partition
function as the trace of a tensor product,

$Z=2^{V}(\cosh(\beta))^{VD}\;{\text{Tr}}\left[\prod_{x}T^{(x)}_{n_{x-\hat{1},1},n_{x,1},\dots,n_{x,D}}\right].$
(80)

Tr means contractions (sums over 0 and 1) over the link indices (the
$n_{x,\mu}$s). The local tensor $T^{(x)}$ has $2D$ indices. The explicit form
is

$\displaystyle
T^{(x)}_{n_{x-\hat{1},1},n_{x,1},\dots,n_{x-\hat{D},D},n_{x,D}}=$
$\displaystyle(\sqrt{\tanh(\beta)})^{n_{x,{\rm in}}+n_{x,{\rm out}}}$
$\displaystyle\times\delta({\rm mod}[n_{x,{\rm in}}+n_{x,{\rm out}},2])$ (81)

with the definitions

$\displaystyle n_{x,{\rm in}}$
$\displaystyle\equiv\sum_{\mu=1}^{D}n_{x-\hat{\mu},\mu}$ (82) $\displaystyle
n_{x,{\rm out}}$ $\displaystyle\equiv\sum_{\mu=1}^{D}n_{x,\mu}.$

The notions of “in” and “out” refer to position of the link with respect to
the positive directions. The basic tensors and their assembly in two and three
dimensions are illustrated in Figs. 6-9.

Any kind of boundary condition can be accommodated by adapting the method of
integration to the link configuration at the boundary. This will be discussed
in subsection V.3.

Figure 6:  Basic tensor for $D=2$. The indices of the tensor are shown as red
“legs” of the purple “body”. Diagrams of tensors this way are ubiquitous in
the literature. The lattice is shown in black, with the centers of the
plaquettes, and the half-way marks on the links, indicated by black dots.
Figure 7:  Tensor assembly for $D=2$. The crosses mean contraction. Figure 8:
Basic tensor for $D=3$. The six indices of the local tensor are shown as the
red legs, and live along the links of the lattice. Figure 9:  Tensor assembly
for $D=3$. The crosses mean contraction.

The Kronecker delta in Eq. (V.1) implies the discrete conservation law

$\sum_{\mu}(n_{x,\mu}-n_{x-\hat{\mu},\mu})=0\ {\rm mod}\ 2,$ (83)

which we also call the “tensor selection rule”. It implies that only an even
number of $n$s are allowed to take on the value one. For instance, for $D=2$,
there are in principle 16 tensor values, however, only eight are nonzero, one
with all four indices as zero (zeroth order in $\tanh(\beta)$), six with two
zeros and two ones (linear in $\tanh(\beta)$), and one with four ones
(quadratic in $\tanh(\beta)$). Note that if a symmetry breaking term like a
magnetic field coupling to the total spin is introduced, then all the 16
tensors elements will be nonzero generally.

There are also tensor formulations which use the singular value decomposition
(See A.3) on each nearest-neighbor Boltzmann factor to factorize the spins.
The nearest neighbor interactions can be represented as a matrix, whose
indices are the spin variables themselves,

$\displaystyle
e^{\beta\sigma_{x+\hat{\mu}}\sigma_{x}}=\begin{pmatrix}e^{\beta}&e^{-\beta}\\\
e^{-\beta}&e^{\beta}\end{pmatrix}_{(x,\mu)}.$ (84)

One can then of course perform the singular value decomposition to get,

$\displaystyle
e^{\beta\sigma_{x+\hat{\mu}}\sigma_{x}}=\sum_{\alpha,\beta}U_{\sigma_{x+\hat{\mu}}\alpha}\lambda_{\alpha\beta}U^{T}_{\beta\sigma_{x}},$
(85)

and in the case of the Ising model the singular value decomposition has the
same left and right unitary matrices. This factorizes the spins and allows for
the definition of the matrix,

$\displaystyle W_{\sigma_{x}\alpha}\equiv
U_{\sigma_{x}\alpha}\sqrt{\lambda_{\alpha}}.$ (86)

The main local tensor is then the contraction of all $W$ matrices which have a
common spin,

$\displaystyle T_{ij\cdots
N}=\sum_{\sigma_{x}}W_{\sigma_{x}i}W_{\sigma_{x}j}\cdots W_{\sigma_{x}N}.$
(87)

In fact, the singular value decomposition in this case can be completely
related to the expansion in Eq. (77). Consider the character expansion from
Sec. V.1 for the nearest neighbor interaction,

$\displaystyle
e^{\beta\sigma_{x+\hat{\mu}}\sigma_{x}}=\sum_{n_{x,\mu}=0}^{1}H_{x+\hat{\mu}}C_{n_{x,\mu}}(\beta)H_{x}$
(88)

with matrices,

$\displaystyle H_{x}\equiv\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\\
1&-1\end{pmatrix}$ (89)

whose columns are the normalized eigenvectors of the Pauli-$x$ matrix, and
$C_{0}=2\cosh\beta$ and $C_{1}=2\sinh\beta$. The elements of this matrix are
from the four values $n_{x,\mu}$ and $\sigma_{x}$ can take, normalized to be
unitary. Since the variables are compact, operators of them have a discrete
spectrum, and the character expansion is exactly the spectrum decomposition
for this matrix. For any matrix, $M$, by definition the singular value
decomposition is found from the eigenvalues and eigenvectors of $MM^{\dagger}$
and $M^{\dagger}M$. In this case (we drop any specific spacetime lattice
indices for these steps since they are completely general),

$\displaystyle\sum_{\sigma_{j}}e^{\beta\sigma_{i}\sigma_{j}}e^{\beta\sigma_{j}\sigma_{k}}$
$\displaystyle=\sum_{\sigma_{j}}\sum_{n,m}\sigma_{i}^{n}C_{n}\sigma_{j}^{n}\sigma_{j}^{m}C_{m}\sigma_{k}^{m}$
$\displaystyle=\sum_{n,m}\sigma_{i}^{n}C_{n}C_{m}\sigma_{k}^{m}\delta_{n,m}$
$\displaystyle=\sum_{n}\sigma_{i}^{n}C_{n}^{2}\sigma_{k}^{n}.$ (90)

Then the singular values are given by $\lambda_{1}=2|\cosh\beta|$ and
$\lambda_{2}=2|\sinh\beta|$, as one would expect from Eq. (85) and Eq. (88).

###  V.2 Remarks on the forms of duality

The tensor representation can be used to reproduce the set of closed paths
appearing in the expansion in powers of $\tanh(\beta)$ Itzykson and Drouffe
(1991); Parisi (1998) for the Ising model. The links associated with the set
of indices $n_{x,\mu}=1$ form a graph (a set of sites connected by links). The
selection rule means that each site is attached to an even number of nonzero
links. These graphs are closed paths with specific connectivity which can in
principle be enumerated order by order in their length using geometric
constructions and combinatoric techniques.

One can try to construct these closed paths by assembling the most elementary
contributions, namely closed loops on a single plaquette. The way they can be
assembled depends on the dimension. For instance for $D=2$, we can decide that
when two loops around two plaquettes share a link, this link is “erased” from
the path. Alternatively, one can introduce “dual variables” originally
concieved by Kramers and Wannier Kramers and Wannier (1941) which are spins on
the dual lattice located at the centers of the plaquettes. Each dual spin is
then associated with a closed plaquette loop around it on the original
lattice.

Furthermore, the Pontryagin dual variables $n_{x,\mu}$ can be expressed as
$(1-\tilde{\sigma}\tilde{\sigma}^{\prime})/2$ with $\tilde{\sigma}$ and
$\tilde{\sigma}^{\prime}$ the two dual spin variables connected by a dual link
crossing the original link. A picture showing the dual variables in their
different locations can be seen in Fig. 10. This illustrates that the notions
of duality are often combined in a way that may be seen as confusing at first
sight. Note that the dual domains with a given dual spin have boundaries that
can be interpreted as the closed paths of the of the original model.

Figure 10: The distinct dual variables shown on the original lattice (black
lines), and the dual lattice (gray lines). Two dual spins are denoted with a
tilde on the dual sites, an original spin, $\sigma$, on the original lattice,
and the Pontryagin dual shown in red crossing the dual link. A closed loop of
Pontryagin duals is shown all in red generally.

The questions of completeness and multiplicity need to be addressed with
specific boundary conditions.

Similarly in higher dimensions it is possible to introduce dual spins with
interactions involving $2(D-1)$ spins. For $D=3$, this leads to a gauge theory
with plaquette interactions Wegner (1971). It is also possible to introduce
dual variables for Ising models with arbitrary all-to-all spin interactions
Meurice (1994). Duality questions related to Gauss’s law will be discussed in
Sec. XII.4.

###  V.3 Boundary conditions

In Eq. (80), the trace is a sum over all the link indices. We need to specify
the boundary conditions. Periodic boundary conditions (PBC) allow us to keep a
discrete translational invariance: the tensors themselves are translation
invariant and assembled in the same way at every site. Open boundary
conditions (OBC) can also be implemented by introducing new tensors that can
be placed at the boundary. Their construction is similar to the tensors in the
bulk. The only difference is that the “outside links” which would be attached
at sites on the boundary have an index set to zero. With the normalization
introduced in Eq. (77) the indices carrying a zero index carry a unit weight.
This construction can be understood as decoupling the system from a larger
environment by setting $\beta$ on the links connecting to this environment to
zero because $\tanh(0)^{0}=1$ and $\tanh(0)=0$. This is illustrated in Fig.
11.

Figure 11: Assembling the translation invariant tensor with PBC (top), or
using new tensors at the boundary for OBC (bottom). Tensors are assumed to be
put on each lattice site.

###  V.4 Exact blocking

An important feature of the tensor representation presented in Sec. V.1 is
that it allows an exact local blocking procedure in contrast to what can be
done in configuration space 222We want to make clear that this section is not
describing a RG transformation but rather an exact reorganization of the
computation of the partition function in a way that performs the integration
of some degrees of freedom corresponding to increasing distance scales. As
explained in Sec. VI, truncations need to be applied to define a RG
transformation.. We divide the original lattice into cells having a linear
size of two lattice spacings (“blocks”) in such a way that the boundaries are
half-way between lattice sites. In other words, the boundaries are normal to
links and cut them in the middle. As an example for $D=2$, the blocks are
square enclosing four sites and their four sides cross eight links. This cell
partition separates the link degrees of freedom into two disjoint categories:
those completely inside the block, which can be integrated over, and those
shared by neighboring blocks and kept to communicate between the blocks. This
is a generic feature of the method Meurice (2013) which motivated a systematic
study of lattice models Liu  _et al._ (2013). Note that when translation
invariance is present, all the blocks are identical and we only need to do one
calculation.

Figure 12: Graphical representation of the block (dotted square) and
$T^{\prime}_{XX^{\prime}YY^{\prime}}$ as in Meurice (2013).

We now describe explicitly this exact blocking for $D=2$ using generic
notations inspired by those of Refs. Xie  _et al._ (2012); Meurice (2013). The
extension to higher dimensions on hypercubic lattices is straightforward. We
contract the four tensors located at the four sites inside the block along the
four indices located strictly inside the block as well. The remaining eight
indices associated with the eight links piercing the block are left as new
degrees of freedom. By taking the tensor product between the two indices in
each of the four directions coming out of the block, we obtain a new rank-4
tensor $T^{\prime}_{XX^{\prime}YY^{\prime}}$. In the case of the Ising model
(See Eq. (V.1)), each index now takes four values. This provides a simple
isotropic formula:

$\displaystyle
T^{\prime}_{X({x_{1}}{x_{2}})X^{\prime}(x_{1}^{\prime}x_{2}^{\prime})Y(y_{1}y_{2})Y^{\prime}(y_{1}^{\prime}y_{2}^{\prime})}=$
(91)
$\displaystyle\sum_{x_{U},x_{D},x_{R},x_{L}}T_{x_{1}x_{U}y_{1}y_{L}}T_{x_{U}x_{1}^{\prime}y_{2}y_{R}}T_{x_{D}x_{2}^{\prime}y_{R}y_{2}^{\prime}}T_{x_{2}x_{D}y_{L}y_{1}^{\prime}},$

where $X(x_{1}x_{2})$ is a notation for the product states _i.e._ , if we
regroup the indices with the same parity together, $X(00)=1,\ X(11)=2,\
X(10)=3,\ X(01)=4$. This contraction and redefinition, relative to the block,
is illustrated in Fig. 12.

After this blocking, the partition function can again be written as

$Z=\text{Tr}\prod_{2x}T^{\prime(2x)}_{XX^{\prime}YY^{\prime}}\ ,$

where $2x$ denotes the sites of the coarser lattice with twice the lattice
spacing of the original lattice. This coarse-graining provides an exact
representation of the original partition function. However, the number of
states associated with each index is the square of the number of states in the
original tensor. If this exact procedure had to be carried numerically, this
rapid growth would quickly run into practical limitations. Truncations are
thus necessary. It is important to appreciate that truncations are the only
approximations that will be needed. We now proceed to discuss truncations.

##  VI Tensor renormalization group

###  VI.1 Block-spinning through SVD

Once partition functions and physical quantities are expressed as tensor
networks, one needs to contract them to obtain numerical values. However,
contracting exactly requires an extraordinary amount of computational
resources. Then, in this subsection we work through the original idea of how
one can perform tensor contractions approximately, and how truncations appear
in this approximation. To do that a coarse-graining algorithm for tensor
networks is introduced. The original idea of such an approach was proposed by
Levin and Nave in Levin and Nave (2007), where tensor networks are simply
coarse-grained by using the SVD (see Sec. A.3). This method is similar in
spirit to the real-space renormalization group approach, and, in this sense,
it is called the tensor renormalization group 333 There are several numerical
renormalization group methods that can be regarded as ancestors of the TRG;
e.g. Wilson (1975); White (1992, 1993); Nishino (1995); Nishino and Okunishi
(1996); Wang and Xiang (1997). See also references the end of Sec. IV.3 and
Ueda  _et al._ (2014).

In the standard renormalization group procedure Wilson and Kogut (1974), the
blocking process is supplemented by a sorting of the resulting information
according to their degree of relevance. As far as universal properties
characterizing the continuum limit, it is acceptable to discard the
information which only reflects the microscopic details of a specific lattice
formulation. In the context of the tensor formulation discussed above, it
means that, possibly after a certain number of exact contractions, we need to
restrict the number of states associated with the tensor indices to a fixed
number $D_{\text{cut}}$. We are then mapping a tensor with
$D_{\text{cut}}^{2D}$ entries into another tensor of the same shape and the
question of fixed points becomes meaningful. An important goal of the
renormalization procedure is to identify fixed points. Note that the updated
tensor remains a local object which supersedes the notion of action or
Hamiltonian.

Here we assume that a partition function is expressed as a two dimensional
tensor network with bond dimension $D_{\text{cut}}$:

$\displaystyle Z=\text{Tr}\prod_{x}T^{(x)}_{xx^{\prime}yy^{\prime}},$ (92)

and that periodic boundary conditions are imposed in all directions. Hereafter
in this subsection, we explain the algorithm of the original TRG proposed by
Levin and Nave.444Coarse-graining approaches for tensor networks are
generically referred to as TRG. Then “TRG” does not necessarily identify a
specific algorithm.

First, the tensor $T$ can be regarded as a $D_{\text{cut}}^{2}\times
D_{\text{cut}}^{2}$ matrix, and can be can be approximately decomposed using
the SVD in two ways:

$\displaystyle
T_{x^{\prime}yxy^{\prime}}=M^{[13]}_{(x^{\prime}y)(xy^{\prime})}\approx\sum_{m=1}^{D_{\text{cut}}}S^{[1]}_{(x^{\prime}y)m}\lambda^{[13]}_{m}S^{[3]}_{m(xy^{\prime})},$
(93) $\displaystyle
T_{y^{\prime}x^{\prime}yx}=M^{[24]}_{(y^{\prime}x^{\prime})(yx)}\approx\sum_{m=1}^{D_{\text{cut}}}S^{[2]}_{(y^{\prime}x^{\prime})m}\lambda^{[24]}_{m}S^{[4]}_{m(yx)},$
(94)

where $\lambda^{[13]}$ and $\lambda^{[24]}$ are the singular values (assumed
to be the descending order:
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{D_{\text{cut}}^{2}}\geq 0$),
and $S^{[1]}$, $S^{[2]}$, $S^{[3]}$, $S^{[4]}$ are unitary matrices. Here the
degree of the approximation is set to be $D_{\text{cut}}$ for simplicity. Of
course one can freely choose this parameter, and it becomes the bond dimension
of the coarse-grained tensors.

Using the decomposed tensors $S^{[i]}$ ($i=1,2,3,4$), a coarse-grained tensor
is defined by

$\displaystyle T^{\mathrm{new}}_{xx^{\prime}yy^{\prime}}$
$\displaystyle=\sqrt{\lambda^{[13]}_{x}\lambda^{[24]}_{y}\lambda^{[13]}_{x^{\prime}}\lambda^{[24]}_{y^{\prime}}}$
$\displaystyle\times\sum_{a,b,c,d=1}^{D_{\text{cut}}}S^{[3]}_{x(cd)}S^{[4]}_{y(da)}S^{[1]}_{(ab)x^{\prime}}S^{[2]}_{(bc)y^{\prime}}.$
(95)

The number of tensors on the network is now reduced by $1/2$, and the bond
dimension of $T^{\mathrm{new}}$ is a free parameter and here set to be the
same as that of $T$, $D_{\text{cut}}$. Repeating this procedure, one can make
an effective tensor network that consists of a few tensors, and then one can
take the contraction of the tensor indices. A graphical explanation of the TRG
is given in Fig. 13.

Figure 13:  A coarse-graining step for the tensor network. Circles represent
tensors, and closed indices should be contracted. The definitions of the unit
vectors for the original and the coarse-grained network are also shown. The
tensor indices are shown in the same manner as Eq. (VI.1).

After a coarse-graining step, the network is rotated by $45$ degrees, and then
how one defines the new unit vectors is one’s choice. One possible way is to
define them by $\hat{1}^{\star}=\hat{1}+\hat{2}$ and
$\hat{2}^{\star}=\hat{1}-\hat{2}$, where $\hat{1}$ $\left(\hat{2}\right)$ and
$\hat{1}^{\star}$ $\left(\hat{2}^{\star}\right)$ are the unit vector along the
$\hat{1}$\- ($\hat{2}$-)direction of the original lattice and that on the
coarse-grained lattice, respectively (see Fig. 13). Using this definition, the
orientation of the network is recovered after every two coarse-graining steps;
i.e. $\hat{1}^{\star\star}=\hat{1}^{\star}+\hat{2}^{\star}=2\cdot\hat{1}$ and
$\hat{2}^{\star\star}=\hat{1}^{\star}-\hat{2}^{\star}=2\cdot\hat{2}$.

In this procedure the relevant degrees of freedom are decided by the SVD
during the decomposition of the $T$ tensor into an intermediate sum-over-
states. In the next section we will see there are improved methods to pick out
the relevant states during truncation.

###  VI.2 Optimized truncations

The previous section discussed the first work at a renormalization group
procedure using a tensor formulation. Here we will discuss refinements which
occurred later. These refinements incorporate an environment tensor from which
relevant states are determined and kept.

We begin with assuming that one has completed contraction of a single block.
Each tensor index of the blocked tensor now posses $D_{\text{cut}}^{2}$
states. The next step is to find a way to decide which of the
$D_{\text{cut}}^{2}$ states (or possibly a linear combination of them)
corresponding to the tensor product associated with each index should be kept.
Ideally, we would like to address this question in terms of the environment of
a single tensor. We can write

$\text{Tr}\prod_{x}T^{(x)}=\sum_{XX^{\prime}YY^{\prime}}E_{XX^{\prime}YY^{\prime}}T^{\prime(0)}_{XX^{\prime}YY^{\prime}},$
(96)

where $E_{XX^{\prime}YY^{\prime}}$ is obtained by tracing all the indices
except for four pairs of indices associated to a single block that we have
located at the origin, and $T^{\prime(0)}$ is the first block of four tensors.
This is illustrated in Fig. 14.

Figure 14: A blocked tensor and its environment.

Constructing an environment from which important states are chosen is
reminiscent of the finite-size density matrix renormalization group.
Obviously, the environment tensor $E_{XX^{\prime}YY^{\prime}}$ is very close
to what are ultimately trying to calculate and so it is not immediately
available. The only purpose of using $E_{XX^{\prime}YY^{\prime}}$ is to rank
order the states of the tensor product and its exact form does not appear in
later calculations. Consequently, the exact form may not be too important and
as a first step, we can try approximations.

A very simple approximation is to ignore the details of the environment and
use Meurice (2013)

$E^{app.}_{XX^{\prime}YY}=C\delta_{XX^{\prime}}\delta_{YY^{\prime}},$ (97)

for some positive constant $C$. We can then optimize the truncation by
maximizing the approximate partition function expressed in terms of the trace
of a matrix $G$ such that

$\text{Tr}\
G=(1/C)\sum_{XX^{\prime}YY^{\prime}}E^{app.}_{XX^{\prime}YY^{\prime}}T^{\prime(0)}_{XX^{\prime}YY^{\prime}},$
(98)

which can be achieved with

$G_{XX^{\prime}}=\sum_{Y}T^{\prime(0)}_{XX^{\prime}YY}.$ (99)

By looking at the expression in terms of the original tensors one realizes
that $G_{XX^{\prime}}$ is in fact the square of another matrix Meurice (2013)
and if the eigenvalues of this matrix are real then all the eigenvalues of $G$
are positive and we can optimize the truncations by selecting the states
corresponding to the largest eigenvalues of $G$.

A more refined approximation is to assume that the environment is a “mirror
image” of the tensor itself Xie  _et al._ (2012):

$E^{app.}_{XX^{\prime}YY}=C^{\prime}T^{\prime(0)\star}_{XX^{\prime}YY^{\prime}}.$
(100)

The trace of $G$ can then be identified with the tensor norm:

$\text{Tr}\
G=\sum_{XX^{\prime}YY^{\prime}}T^{\prime(0)}_{XX^{\prime}YY^{\prime}}T^{\prime(0)\star}_{XX^{\prime}YY^{\prime}}=||T^{\prime(0)}||^{2},$
(101)

which is clearly a sum of positive terms. This can be accomplished with the
Hermitian matrix

$G_{XX^{\prime}}=\sum_{X^{\prime\prime}YY^{\prime}}T^{\prime(0)}_{XX^{\prime\prime}YY^{\prime}}T^{\prime(0)\star}_{X^{\prime}X^{\prime\prime}YY^{\prime}}.$
(102)

The problem is then reduced to selecting the states that provide the best
approximation of $\text{Tr}G$ which is obvious when all the eigenvalues are
positive.

The procedure that we just described is isotropic. It is however possible to
coarse grain one direction at a time Xie  _et al._ (2012) in order to reduce
the size of the summed expressions. For instance, the summation over the
tensor product indices $Y$ and $Y^{\prime}$ in Eq. (102) could be first
replaced by summations over single indices,

$\displaystyle
M_{X(x_{1},x_{2})X^{\prime}(x^{\prime}_{1},x^{\prime}_{2})yy^{\prime}}=\sum_{a}T_{x_{1}x^{\prime}_{1}ya}T_{x_{2}x^{\prime}_{2}ay^{\prime}}.$
(103)

The tensor $M$ can then be used in Eq. (102) in place of $T^{(0)}$ as the
“blocked” tensor to find the most relevant states for the $X$ indices. This
provides a coarse-graining in the first direction. It is then necessary to
coarse-grain in the second direction using sums over single indices in the
first direction. As our discussion is focused on $D=2$, these two steps
constitute a coarse-graining that doubles the lattice spacing in all
directions. In Xie  _et al._ (2012), these calculations were conducted using
higher-order generalizations of the SVD method called the higher-order tensor
renormalization group (HOTRG).

A better description of the environment can be reached by following a local
truncation procedure as described above for a sufficiently large but finite
number of times at which point it is assumed that there is no environment and
Eq. (97) can be used. In other words, by working with a finite lattice, the
procedure is terminated by approximating the partition function as a trace of
the last coarse-grained expression for the tensor. It is then possible to move
“backward” Xie  _et al._ (2009); Zhao  _et al._ (2010) and reconstitute the
approximate environment of a single tensor coarse-grained one less time.
Explicit expressions based on the HOTRG construction can be found in Xie  _et
al._ (2012). This can be pursued recursively until we reach the first coarse
graining level illustrated in Fig. 14.

The analogy with the backward propagation used in machine learning has been
exploited to design new algorithms recently Chen  _et al._ (2020). More
generally, algorithmic improvement in the TRG context is a subject of active
investigation Bal  _et al._ (2017); Fishman  _et al._ (2018); Morita and
Kawashima (2020). The idea that the bound dimension can be regarded as as a
relevant direction which can be used to obtain data collapse has been
discussed in Vanhecke  _et al._ (2019). For earlier developments see
Tagliacozzo  _et al._ (2008); Pollmann  _et al._ (2009); Pirvu  _et al._
(2012). We expect that these considerations will provide a more systematic
understanding of the truncation errors.

###  VI.3 Higher dimensional algorithms

While the Levin–Nave type TRG can be applied to two dimensional systems,
higher dimensional systems are dealt with using other algorithms. One such
algorithm is the higher-order TRG (HOTRG) Xie  _et al._ (2012) mentioned in
the previous section. Using the HOTRG, in principal, any dimensional tensor
network can be coarse-grained. When a $D$ dimensional tensor network is build
by tensors with the bond dimension $D_{\text{cut}}$, the computational
complexity of the HOTRG is $\mathcal{O}\left(D_{\text{cut}}^{4D-1}\right)$,
and the memory complexity is $\mathcal{O}\left(D_{\text{cut}}^{2D}\right)$.

Recently, cheaper algorithms have been invented. The anisotropic tensor
renormalization group (ATRG) Adachi  _et al._ (2019), whose graphical
description is given in Fig. 15, achieved the time and the memory complexities
$\mathcal{O}\left(D_{\text{cut}}^{2D+1}\right)$ and
$\mathcal{O}\left(D_{\text{cut}}^{D+1}\right)$ that are significant reductions
from the HOTRG. The ATRG introduces an approximation of an approximation, and
indeed, when the bond dimensions are the same, the ATRG is less accurate
compared to the HOTRG. However, thanks to the cheaper complexity, the ATRG
leads to better accuracy with fixed CPU time. Using the ATRG, four dimensional
systems, where the HOTRG is much more expensive, have begun investigation
Akiyama  _et al._ (2019, 2020a, 2020b).

Figure 15:  Graphical description of the ATRG. (a) Decompose $T$ into ($A$ and
$B$) and ($C$ and $D$) with the SVD. (b) Swap the bonds of $B$ and $C$. The
swapping can be done also by the SVD:
$\sum_{i}B_{xyi}C_{ix^{\prime}y^{\prime}}=(BC)_{xyx^{\prime}y^{\prime}}=\sum_{i}B^{\prime}_{x^{\prime}yi}S_{i}C^{\prime}_{ixy_{\prime}}$.
(c) Contract $A$, $B^{\prime}$, $C^{\prime}$, $D$ while truncating the
dimensions of horizontal bonds.

Another approach is coarse-graining on a triad tensor network representation
Kadoh and Nakayama (2019). This paper compares the triad tensor
renormalization group approach to the original HOTRG and the ATRG algorithms.

###  VI.4 Observables with tensors

With the tensor formulation from Sec. V.1, along with the coarse-graining
algorithms from Sec. VI it is possible to calculate derivatives of $\ln Z$, as
well as compute $n$-point correlation functions.

Because the tensor renormalization group process requires renormalizing the
tensor during iterations, in order to compute the logarithm of the partition
function these normalizations must be stored. For a coarse-graining that is
isotropic with $N_{1},N_{2},\ldots N_{D}$ iterations in each direction and
$N=\sum_{i=1}^{D}N_{i}$, let us say the normalizations at each iteration are
$\mathcal{N}^{(0)}$, $\mathcal{N}^{(1)}$, $\mathcal{N}^{(2)}$ etc. starting by
normalizing the initial tensor by $\mathcal{N}^{(0)}$. Then during the coarse-
graining process, $\mathcal{N}^{(0)}$ appears $V=2^{N}$ times, giving an
overall factor of ${\mathcal{N}^{(0)}}^{V}$. Likewise each subsequent
normalization appears $2^{N-n}$ where $n$ denotes the iteration number, _i.e._
$1,\ldots N$. After the final step the total normalization on the effective
tensor is given by

$\displaystyle{\mathcal{N}^{(0)}}^{V}\cdot{\mathcal{N}^{(1)}}^{V/2}\cdot{\mathcal{N}^{(2)}}^{V/4}\cdots\mathcal{N}^{(N)}.$
(104)

The logarithm of the partition function is then given as the logarithm of this
normalization, added to the logarithm of the trace of the final tensor,

$\displaystyle\ln
Z=\sum_{n=0}^{N}2^{N-n}\ln\mathcal{N}^{(n)}+\ln\text{Tr}[T],$ (105)

where $T$ is the final normalized effective tensor, and the trace is the
tensor trace.

Expectation values of $N$-point correlation functions Gu  _et al._ (2008);
Nakamoto and Takeda (2016) are given by ratios of partition functions which
are calculated separately using the tensor renormalization group,

$\displaystyle\langle\sigma^{(1)}\sigma^{(2)}\cdots\sigma^{(N)}\rangle=\frac{Z^{(N)}}{Z}$
(106)

with

$\displaystyle
Z^{(N)}=\sum_{\\{\sigma\\}}\sigma^{(1)}\sigma^{(2)}\cdots\sigma^{(N)}e^{-S}$
(107)

and the positions of the $\sigma$ fields have been suppressed. In terms of
tensors this amounts to $N$ “impure” tensors whose namesake comes from their
altered local constraint. Since $Z^{(N)}$ contains additional spin fields
located at specific sites, those spin fields alter the sum at that site over
the field states, and give

$\displaystyle\sum_{\sigma_{x^{*}}}\sigma_{x^{*}}^{1+\sum_{\mu=1}^{D}n_{x^{*}-\hat{\mu},\mu}-n_{x^{*},\mu}}$
(108) $\displaystyle=\begin{cases}2&\text{ if
}1+\sum_{\mu=1}^{D}n_{x^{*}-\hat{\mu},\mu}-n_{x^{*},\mu}\text{ is even}\\\
0&\text{ otherwise.}\end{cases}$

where $x^{*}$ is the location of the additional spin. These additional spins
act as sources and sinks for the vector fields $n_{x,\mu}$ on the surrounding
links as can be seen from the new constraint where the divergence of the
surrounding $n$ fields must now be one modulo two. To compute the ratio, if
both $Z^{(N)}$ and $Z$ are normalized identically throughout the calculation,
only the ratio of the trace of the final effective tensors is necessary for
the expectation value, since the normalizations would cancel.

As illustrated in Fig. 16 the method can be used for the average energy and
the specific heat, and are compared to the exact solutions by Kaufman’s
formula. The energy is calculated using the TRG with two neighboring impurity
tensors, and the specific heat is given as the numerical $\beta$ derivative of
the energy. To show a graphical example of the contractions performed using
impure tensors, four impurities on a plaquette are decomposed and contracted
to form four impurities again, is illustrated in Fig. 17. Note that the energy
can also be obtained by the numerical derivative of the logarithm of the
partition function, and then the specific heat is the second numerical
derivative. In general, the numerical derivatives cause a loss of
significance, so obtaining the energy as a primary output of the TRG using the
impurity tensor method helps to improve the numerical accuracy. The impurity
tensors destroy the translational invariance, but the effect is local so that
the computational complexity does not drastically increase.

Figure 16:  Energy and specific heat of the two dimensional Ising model on a
$32\times 32$ lattice. $D_{\mathrm{cut}}=32$.  Figure 17:  TRG process for a
tensor network with four impurities on a plaquette.

###  VI.5 Niemeijer-Van Leeuwen equation

In Sec. V.4 we have constructed a coarse-grained tensor that can be used to
give an exact expression for the partition function. Despite the fact that we
integrated over microscopic degrees of freedom, the number of tensor indices
needed for this exact representation grows exponentially with the size of the
blocks. In Secs. VI.1 and VI.2, we introduced truncations where a fixed number
of indices $D_{\text{cut}}$ was kept at each step. This procedure discards
some information but allows us to compare the tensors before and after the
coarse-graining. Typically, the tensors tend to grow exponentially with the
number of coarse graining steps and it is important to renormalize their
absolute size or to only consider their ratios. After such a renormalization
takes place, we obtain a RG transformation. The fixed points of this
transformations are the central objects of the RG approach and it is useful to
compare the TRG equations with standard RG equations due to Niemeijer and van
Leeuwen (NvL)Niemeijer and van Leeuwen (1976).

NvL were aware of the difficulty of controlling the new couplings generated by
the blocking procedure, so they started immediately with the most general
Ising interactions in a finite volume $V$. If no conditions are imposed, there
are as many couplings as Ising configurations so this is not suitable for
numerical purposes. They then introduced 1 in the partition function as in Eq.
(1), in order to define new Ising spins $\\{\sigma^{\prime}\\}$ in a volume
$V^{\prime}=V/b^{D}$ with a new lattice spacing rescaled by the linear size of
the blocks $b$. They were able to give a formal expression for the new
couplings in terms of the original ones as ${\bf K}^{\prime}({\bf K})$.
Strictly speaking there are less couplings after the coarse graining because
they considered the most general case involving the products of spins in
arbitrary domains, but they assumed that only a certain number of quasi-local
couplings were important and had the same form before and after the coarse
graining. In addition they assumed that the dependence of the free energy
density $f({\bf K})$ on these couplings is the same after the coarse graining.
This led to the NvL equation

$f({\bf K})=g({\bf K})+b^{-D}f({\bf K}^{\prime}).$ (109)

The function $g=G/V$, comes from

$G=\sum_{\\{\sigma^{\prime}\\}}\ln\Big{(}\sum_{\\{\sigma\\}}P(\\{\sigma^{\prime}\\},\\{\sigma\\})\exp\big{(}\mathcal{H}(\\{\sigma\\})\big{)}\Big{)},$
(110)

and is defined by the condition

$\sum_{\\{\sigma^{\prime}\\}}\mathcal{H}^{\prime}(\\{\sigma^{\prime}\\})=0,$
(111)

where $\mathcal{H}$ and $\mathcal{H}^{\prime}$ are the Hamiltonians before and
after coarse graining.

Even though computing the new couplings and the functions may be very
difficult in practice (the new couplings are double partition functions), NvL
succeeded to obtain a formal relation which can be iterated and linearized
near a fixed point. This allows us to identify the relevant directions and it
is often taken as the starting point for the introduction of the RG method in
textbooks Cardy (1996).

For the TRG, we can factor out the increasing size of the tensors, for
instance, by imposing the normalization condition

$T_{0000}=1,$ (112)