Journal of Physics:
Condensed Matter



TOPICAL REVIEW • OPEN ACCESS                                                                     You may also like
                                                                                                     - Grassmann corner transfer-matrix
Tensor renormalization group for fermions                                                              renormalization group approach to one-
                                                                                                       dimensional fermionic models
                                                                                                       Jian-Gang Kong, , Zhi-Yuan Xie et al.
To cite this article: Shinichiro Akiyama et al 2024 J. Phys.: Condens. Matter 36 343002
                                                                                                     - Variational Corner Transfer Matrix
                                                                                                       Renormalization Group Method for
                                                                                                       Classical Statistical Models
                                                                                                       X. F. Liu, , Y. F. Fu et al.

View the article online for updates and enhancements.                                                - Building projected entangled pair states
                                                                                                       with a local gauge symmetry
                                                                                                       Erez Zohar and Michele Burrello




                             This content was downloaded from IP address 103.151.173.99 on 05/07/2026 at 10:21
                                                                                                                         Journal of Physics: Condensed Matter

J. Phys.: Condens. Matter 36 (2024) 343002 (31pp)                                                                     https://doi.org/10.1088/1361-648X/ad4760



Topical Review


Tensor renormalization group for
fermions
Shinichiro Akiyama1,2, Yannick Meurice3,∗ and Ryo Sakai4
1
  Center for Computational Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan
2
  Graduate School of Science, The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan
3
  Department of Physics and Astronomy, The University of Iowa, Iowa City, IA 52242, United States of
America
4
  Jij Inc., Bunkyo-ku, Tokyo 113-0031, Japan

E-mail: yannick-meurice@uiowa.edu and akiyama@ccs.tsukuba.ac.jp

Received 6 February 2024, revised 2 April 2024
Accepted for publication 3 May 2024
Published 28 May 2024


Abstract
We review the basic ideas of the tensor renormalization group method and show how they can
be applied for lattice field theory models involving relativistic fermions and Grassmann
variables in arbitrary dimensions. We discuss recent progress for entanglement filtering, loop
optimization, bond-weighting techniques and matrix product decompositions for Grassmann
tensor networks. The new methods are tested with two-dimensional Wilson–Majorana fermions
and multi-flavor Gross–Neveu models. We show that the methods can also be applied to the
fermionic Hubbard model in 1+1 and 2+1 dimensions.
Keywords: tensor networks, lattice gauge theory, relativistic lattice fermions,
Fermi Hubbard model, Grassmann path integrals, sign problems


Contents                                                                                  3.2. The Schwinger model                                        18
                                                                                       4. Improved TRG methods for fermions                               19
1. Introduction                                                               2           4.1. CDL structure on tensor network                            19
2. Formalism for the Grassmann TRG                                            3           4.2. Removal of CDL from network                                19
   2.1. Grassmann tensor network representation                               3           4.3. Bond-weighting technique                                   20
   2.2. Exact contraction                                                     6           4.4. Multilayered tensor network formulations for
   2.3. Model-independent notation                                            8                N f -flavor fermions                                       22
   2.4. Extension to lattice gauge theories                                   9        5. Relativistic models with fermion interactions                   23
   2.5. Approximate contraction by the Levin–Nave                                         5.1. Gross–Neveu model                                          23
        TRG                                                                 11            5.2. QCD in the infinite-coupling limit                         24
   2.6. Approximate contraction by the HOTRG                                14            5.3. NJL model                                                  25
   2.7. Practical remarks                                                   16            5.4. N = 1 Wess–Zumino model                                    26
3. Examples of numerical calculations                                       17            5.5. Non-abelian lattice gauge theories coupled to
   3.1. Wilson–Majorana fermions                                            17                 fermions                                                   26
                                                                                       6. The Hubbard model                                               26
∗
                                                                                          6.1. (1+1)-dimensional model                                    28
    Author to whom any correspondence should be addressed.
                                                                                          6.2. (2+1)-dimensional model                                    28
                                                                                       7. Conclusions                                                     29
                   Original Content from this work may be used under the
                                                                                       Data availability statement                                        29
                   terms of the Creative Commons Attribution 4.0 licence. Any
further distribution of this work must maintain attribution to the author(s) and       Acknowledgment                                                     29
the title of the work, journal citation and DOI.                                       References                                                         29

                                                                                   1                     © 2024 The Author(s). Published by IOP Publishing Ltd
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                   Topical Review



1. Introduction                                                             Fermionic tensor network studies with applications for
                                                                        the Fermi–Hubbard model were discussed in [75, 76]. The
The renormalization group (RG) approach of lattice models               application of the TRG method to fermionic systems rely-
has been crucial to identify their universal critical behavior,         ing on Grassmann numbers was proposed by Gu, Verstraete,
construct their phase diagram and understanding their con-              and Wen in [30, 31] in the context of the tensor product
tinuum limits. The general idea [1–3] consists in integrating           states (TPS) [77, 78] and projected entangled pair states
over some of the microscopic degrees of freedom in order                (PEPS) [79]. TPS and PEPS are known as variational wave-
to obtain an effective theory with a larger lattice spacing and         functions that can efficiently represent the ground state of a
repeat the procedure until one reaches a macroscopic size. This         gapped local Hamiltonian in higher dimensions. They exten-
leads to an RG map connecting effective theories at increas-            ded TPS to fermionic systems and proposed to construct the
ing lattice spacing. The fixed points and relevant directions of        TPS using the Grassmann variables based on the idea of the
the RG maps have universal properties such as similar critical          path-integral formalism of fermionic systems. This generic
exponents observed in very different microscopic setups (e.g.           variational wavefunction is referred to as the Grassmann TPS
magnets and solids).                                                    (GTPS). The TRG, which explicitly includes the Grassmann
   The generic properties of RG maps are very well-                     variables, was devised as a method for approximate contrac-
understood [2–4]. However, the numerical implementation of              tion of the Grassmann tensor network derived from the GTPS.
the partial integration over some microscopic degrees of free-          Note that [30] showed that the fermionic PEPS (fPEPS) pro-
dom (coarse-graining) can be challenging. This requires to              posed in [80] can be classified as a special subclass of GTPS.
parameterize some ‘space of theories’ in terms of effective             fPEPS and GTPS can be considered as an origin of the cur-
couplings and find numerical methods to calculate the ‘new’             rent Grassmann TRG approach [30, 31, 80]. Note that sev-
couplings in terms of the ‘old’ couplings. This provides the            eral kinds of higher-dimensional tensor network ansatzes with
RG map. Simplified RG maps such as various majority rules,              variational methods, some of which are combined with the
bond moving [5, 6] approximate recursions [2], hierarchical             Monte Carlo method, have recently been applied to various
approximations [7–9], local potential approximations [10–12],           lattice gauge theories [81–92]. For a set of lectures notes on
can be easily implemented numerically and demonstrate the               the applications of PEPS in the context of lattice gauge theory
generic validity of the RG approach. However, these approx-             see [93].
imations provide critical exponents which are different from                A practically important aspect of tensor network includ-
the exponents of the approximated models and it is difficult to         ing the TRG is that in principle, the sign problem does not
systematically improve such approximations.                             exist. Although Quantum Monte Carlo simulations are very
   The basic RG ideas were incorporated in variational                  efficient for bosons, they face the negative sign problem when
algorithms designed to construct the ground state wavefunc-             applied to fermions [94]. Introducing auxiliary bosons, the
tion of Hamiltonians for one-dimensional spatial lattices often         Monte Carlo methods can simulate lattice fermion models
called the density matrix RG method [13–15]. This led to the            but it is usually computationally challenging because the res-
development of tensor network ansatzes (e.g. matrix product             ulting bosonic theory becomes completely non-local. For the
states (MPS)) for these wavefunctions. It became clear that the         numerical approaches for lattice fermions based on the Monte
success of the method was linked to its ability to handle effi-         Carlo method, see references [95–98] and references therein.
ciently the entanglement entropy in one spatial dimension [14,          In contrast, the tensor network approach allows us to deal with
16–25]. The tensor technology was also used to reformulate              local fermionic actions because we can directly manipulate the
and coarse-grain classical lattice models [26–34]. This new             Grassmann numbers. Moreover, if the system has translational
approach is called the tensor RG (TRG). The TRG approach                invariance on a lattice, we can access its infinite-volume limit
allows us to interpret and quantitatively realize the traditional       (zero-temperature limit) by using the TRG. As we will see, this
RG ideas in the language of the tensor network [20, 35–41].             advantage of the TRG has been demonstrated in actual numer-
   Stimulated by a certain number of interdisciplinary work-            ical calculations for lattice fermions. Therefore, the tensor net-
shops (see for instance, [42]), tensorial methods became of             work approach has been a powerful candidate for investigating
interest to the lattice gauge theory community. MPS were                systems where the sign problem makes Quantum Monte Carlo
used for the Schwinger model [43–51], non-Abelian gauge                 simulation difficult to use in the first place.
theories [52–55], and the O(3) nonlinear sigma model [56].                  In this topical review, we discuss recent TRG applications
Tensor network techniques for lattice gauge theories are also           for fermionic models. The starting point is a path integral
discussed in [57–62] and reviewed in [63, 64]. At the same              involving Grassmann variables. We will consider models in
time, it was shown that character expansions used in the con-           Euclidean spacetime and the statistical weights in the path
text of the strong coupling expansion [65, 66] could be used            integral will have the generic form e−S for some classical
for TRG approaches of most models studied in lattice gauge              action S. The first step of the tensorial approach consists in
theory [67–71]. This led to a new unified way to understand             expanding the exponential of the non-local parts of S (attached
global, local, continuous, and discrete symmetries [72, 73]             to links or plaquettes of the lattice) in terms of discrete tensor
more generally to the development of tensor lattice field theory        indices and then perform the integration over the original vari-
recently reviewed in [74].                                              ables. With this procedure, the path integral becomes a sum



                                                                    2
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                     Topical Review



of product of local tensors with their indices contracted. For         functions without the Grassmann variables, we will see that
bosonic variables, the tensors are just ordinary functions of          any TRG algorithm can be applied to fermionic path integrals.
the coupling that can be collected in a straightforward man-
ner. However, for Grassmann variables, one needs to keep
                                                                       2.1. Grassmann tensor network representation
track of the signs resulting from various orderings. As we will
show in section 2, this can be handled by introducing auxiliary        For the time being in this section, we use the following simple
Grassmann variables which will be incorporated in the tensors.         quadratic model to explain how to derive the Grassmann tensor
We then proceed to coarse-grain the reformulation of the path          network representation,
integral by performing approximately partial contractions fol-
lowing the early TRG method of Levin–Nave (illustrated in                                   ∑∑
                                                                                             2
                                                                                               [                                                     ]
figure 5) and their higher-order TRG (HOTRG) version (illus-                     S = −t                 ψ̄ (n + ν̂) ψ (n) + ψ̄ (n) ψ (n + ν̂)
                                                                                            n∈Λ ν=1
trated in figure 7). In section 3, we review simple examples                                 ∑
of numerical calculations for free Wilson–Majorana fermions                          +m             ψ̄ (n) ψ (n) ,                                        (1)
(equivalent to the Ising model) and 1+1 QED (the Schwinger                                      n
model).
    An important aspect of the TRG coarse-graining is that it          where we assume that ψ(n) and ψ̄(n) are single-component
can in principle be performed exactly [34, 67]. However, in            Grassmann fields for simplicity. We consider the model on a
practice, the computational cost still scales exponentially with       two-dimensional square lattice Λ with periodic boundary con-
the size of the system because despite the partial integration         ditions. The path integral on the lattice is defined via
over some of the microscopic degrees of freedom, effective                                  ˆ ∏
tensors with more indices are generated. The scaling is less                           Z=          dψ (n) dψ̄ (n) e−S .         (2)
severe, but nevertheless exponential. For this reason and also                                             n
in order to get RG maps relating same size tensors, trunca-
tions are necessary. It has been argued that the inverse of the        Our goal is to rewrite equation (2) introducing local tensors.
truncation size can be considered as a relevant direction which        There are multiple ways to define these local tensors so that
can interfere with the study of fixed points [99]. It has also         the tensor network representation is in general not unique. This
been known [27, 28] that short-range entanglement can remain           situation is the same as in constructing tensor network repres-
present during TRG iterations and generate unphysical fixed            entations for theories that do not include fermions.
points. Improvements of this situation for fermionic theories              One methodology was given by Shimizu and
are discussed in section 4. We introduce the corner double line        Kuramashi [117] and subsequently refined by Takeda and
(CDL) fixed point and then discuss methods to remove them.             Yoshimura by explicitly introducing auxiliary Grassmann
This includes the tensor network renormalization (TNR) [20],           variables [102]. Here, we review the formalism in [102].
the loop-TNR [100], and the gilt-TNR [101] algorithm.                  Firstly, the hopping terms are decomposed introducing auxil-
    Applications for relativistic models are presented in              iary Grassmann variables as
section 5 which covers the Gross–Neveu model [102, 103],
QCD at infinite coupling [104], the Nambu–Jona-Lasinio                 e tψ̄(n+ν̂)ψ(n)
(NJL) model [105], the N = 1 Wess–Zumino model [106]                              ∑ (ˆ                                               )iν (n)
                                                                                  1
                                                                                       √
and non-abelian gauge theories with fermion [107]. Finally,                  =                       tψ̄ (n + ν̂) dΦ̄ν (n + ν̂)
the Fermi-Hubbard model which has been extensively studied                       iν (n)=0
                                                                                     (ˆ                          )iν (n)
in condensed matter community [108–114] is considered in                                    √                              (                         )iν (n)
section 6 following the approach of [115] in 1+1 dimensions                      ×              tψ (n) dΦν (n)                 Φ̄ν (n + ν̂) Φν (n)             ,
and [116] in 2+1 dimensions.                                                                                                                              (3)
                                                                           tψ̄(n)ψ(n+ν̂)
                                                                       e
2. Formalism for the Grassmann TRG                                                ∑ (ˆ                                  )jν (n)
                                                                                  1
                                                                                       √
                                                                             =                       tψ̄ (n) dΨ̄ν (n)
We begin by explaining how to express fermionic path integ-                      jν (n)=0
rals in the language of tensor networks. Because of the nil-                         (ˆ                                         )jν (n)
                                                                                            √
potency of the Grassmann variables, they can be straightfor-                     ×              tψ (n + ν̂) dΨν (n + ν̂)
wardly rewritten with finite-dimensional tensors. We show
                                                                                  (                    )jν (n)
how to construct the exact contraction between these funda-                      × Ψ̄ν (n) Ψν (n + ν̂)         ,                                          (4)
mental tensors, which is necessary to reproduce original path
integrals. However, the exact contraction is usually prohibited        where Ψν , Ψ̄ν , Φν , and Φ̄ν are single-component Grassmann
in practice because it requires exponentially large computa-           variables along the ν-directional link. The bit iν (n) (jν (n))
tional resources. One of the promising approaches is the TRG           labels the Taylor expansion of the forward (backward) hop-
which carries out these contractions approximately. Although           ping term on the link (n, n + ν̂). These bits are nothing but
the TRG algorithms are often used to compute the partition             the occupation numbers. We use capital Greek letters for the


                                                                   3
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                                                 Topical Review



auxiliary Grassmann variables (Ψν , Ψ̄ν , etc) as opposed to the                                                 equation (6) is provided in figure 1(notice the left-right asym-
original microscopic Grassman variables which are lowercase                                                      metry). Z in equation (2) is reproduced by summing over all
(ψ, ψ̄). In the rest of this section,´ in order to have more com-                                                bits and integrating over all auxiliary Grassmann variables.
pact notations, we will drop the ’s in equations (3) and (4).                                                    This situation is symbolically expressed as
                                                                                                                                            ∑ ˆ ∏
In other words, when iν (n) = 1 or jν (n) = 1, the differentials
of the auxiliary´ Grassmann variables such as dΨν , should be                                                                        Z=                   Tn ,                (10)
understood as dΨν . In the right-hand sides of equations (3)                                                                                       {i1 ,j1 ,i2 ,j2 }     n
and (4), the original fields living on different sites are decom-
posed into different Grassmann-even pairs. Therefore, ψ(n)                                                       where
and ψ̄(n) in equation (2) can be easily integrated at each lat-
tice site n independently. At each lattice site n, we consider the                                                                        ∑                  ∏ ∑
                                                                                                                                                               1              ∑
                                                                                                                                                                              1
                                                                                                                                                         =                          ,               (11)
following integral
                                                                                                                                     {i1 ,j1 ,i2 ,j2 }       n,ν iν (n)=0 jν (n)=0

 Tn;(i1 j1 )(i2 j2 )(i ′ j ′ )(i ′ j ′ )                                                                              ´
                                                                                                                 and ’s are over the auxiliary Grassmann variables, which is
         ˆ            1 1       2 2
                                         ∏ (√           )iν (√             )j ν                                  implicit in equation (8).
    = dψdψ̄ e−mψ̄ψ                            tψdΦν (n)        tψ̄dΨ̄ν (n)                                           This kind of Grassmann tensor network formulation has
                                                ν
                                                                                                                 been widely applied: the Schwinger model [117–119]5 ,
                (√                             )j ′
                                           )iν′ (√
           ×          tψ̄dΦ̄ν (n) tψdΨν (n) ν                                                                    the Gross–Neveu model [102, 120], free Wilson fermi-
            (                    )iν (                    )jν                                                    ons [121, 122], the N = 1 Wess–Zumino model [106], the
           × Φ̄ν (n + ν̂) Φν (n)       Ψ̄ν (n) Ψν (n + ν̂) ,                                           (5)       NJL model [105], Wilson–Majorana fermions [120], infinite-
                                                                                                                 coupling QCD [104], and SU(2) lattice gauge theory with
where we have introduced several shorthand notations such as
                                                                                                                 reduced staggered fermions [107].
iν := iν (n), jν := jν (n), iν′ := iν (n − ν̂), and jν′ := jν (n − ν̂).
                                                                                                                     There are alternative ways to formulate the tensor network
We regard equation (5) as a fundamental tensor describing the
                                                                                                                 with Grassmann variables. Meurice introduced a multilinear
path integral in equation (2). Equation (5) can be written as
                                                                                                                 combination of Grassman variables to define a fundamental
                                                                                                                 tensor [123] and Bao summarized the contraction, decom-
           Tn;(i1 j1 )(i2 j2 )(i ′ j ′ )(i ′ j ′ )
                                     1 1    2 2                                                                  position, and conjugation for the tensors associated with the
                 = Tn;(i1 j1 )(i2 j2 )(i ′ j ′ )(i ′ j ′ ) Gn;(i1 j1 )(i2 j2 )(i ′ j ′ )(i ′ j ′ ) ,   (6)       Grassmann variables [124]. Akiyama and Kadoh then gave a
                                                1 1       2 2                      1 1      2 2
                                                                                                                 general methodology to derive the tensor network representa-
with                                                                                                             tion for fermionic path integrals based on the concrete defini-
                                                                                                                 tion of the Grassmann tensor [125]. Here, we follow the form-
 Tn;(i1 j1 )(i2 j2 )(i ′ j ′ )(i ′ j ′ )                                                                         alism provided in [125].
                        1 1         2 2
                                                                                                                     They also decompose the hopping terms but use
       √ ∑ (iν +jν +iν′ +jν′ )               j i +j ′ +j ′ +j j ′ +j ′ +i ′ j ′
      = t ν                         (−1) 1 ( 2 1 2 ) 2 ( 1 2 ) 1 2                                                                     ˆ
        [                                                                                                             e tψ̄(n+ν̂)ψ(n) = dΦ̄ν (n) dΦν (n) e−Φ̄ν (n)Φν (n)
       × −mδi1 +i2 +j1′ +j2′ ,0 δi1′ +i2′ +j1 +j2 ,0
                                                  ]                                                                                                √                          √
       + δi1 +i2 +j1′ +j2′ ,1 δi1′ +i2′ +j1 +j2 ,1 ,                                                   (7)                              × e tψ̄(n+ν̂)Φ̄ν (n) e tψ(n)Φν (n) ,                        (12)
                                                                                                                                        ˆ
 Gn;(i1 j1 )(i2 j2 )(i ′ j ′ )(i ′ j ′ )
                        1 1         2 2                                                                                e tψ̄(n)ψ(n+ν̂) = dΨ̄ν (n) dΨν (n) e−Ψ̄ν (n)Ψν (n)
                         i1                 j                   i            j                j′
      = dΦ1 (n) dΨ̄1 (n) 1 dΦ2 (n) 2 dΨ̄2 (n) 2 dΨ1 (n) 1                                                                                          √                         √

                              i1′                   j2′             i2′
                                                                                                                                            ×e           tψ̄(n)Ψν (n)
                                                                                                                                                                        e−    tψ(n+ν̂)Ψ̄ν (n)
                                                                                                                                                                                                ,   (13)
           × dΦ̄1 (n) dΨ2 (n) dΦ̄2 (n)
             ∏(                       )iν (                    )jν                                               instead of equations (3) and (4). Again, Ψν , Ψ̄ν , Φν , and
           ×      Φ̄ν (n + ν̂) Φν (n)       Ψ̄ν (n) Ψν (n + ν̂) . (8)
                                                                                                                 Φ̄ν are the single-component Grassmann variables. We are
                 ν
                                                                                                                 now ready to carry out the integral over ψ(n) and ψ̄(n) in
The Kronecker deltas in equation (7) imply that                                                                  equation (2) at each site n independently as before. The fun-
                                                                                                                 damental tensor is defined via
             i1 + j1 + i2 + j2 + i1′ + j1′ + i2′ + j2′ mod 2 = 0.                                      (9)                      ˆ                ∏ √
                                                                                                                           Tn = dψdψ̄ e−mψ̄ψ        e tψ̄Φ̄ν (n−ν̂)
In other words, the fundamental tensor Tn in equation (6) is                                                                             √          √
                                                                                                                                                        ν
                                                                                                                                                                √
                                                                                                                                                      tψΦν (n) − tψ Ψ̄ν (n−ν̂)
always Grassmann-even. Note that the Kronecker deltas in                                                                           ×e     tψ̄Ψν (n)
                                                                                                                                                           e             e                  .       (14)
equation (7) are not mod 2 because one may not have more
than one power of the Grassman numbers. Following [102],
we call T n in equation (7) as the bosonic part of Tn and Gn as
the Grassmann part. Note the order of the Grassmann meas-                                                        5 Although the previous study in [117] does not appear to introduce auxili-
ures in equation (8). Although Φν (n)iν and Ψ̄ν (n)jν (ν = 1, 2)                                                 ary Grassmann variables explicitly, it gives a formulation that is equivalent
can be integrated within equation (8), the integration shall not                                                 to introducing auxiliary variables. Consequently, their fundamental tensor in
be performed at this stage. The graphical representation of                                                      equation (31) in [117] has the same structure with equation (6).

                                                                                                             4
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                    Topical Review




Figure 1. Graphical representation of the fundamental tensor Tn . (Top) The fundamental tensor in equation (6). Each line denotes the
single-component auxiliary Grassmann measure and the occupation number. Grassmann-even factors are associated with forward lines.
(Bottom) The fundamental tensor in equation (15). Each line denotes the single-component auxiliary Grassmann variable.




                                                                     5
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                                                                   Topical Review



Expanding the integrand, one obtains                                                                     2.2. Exact contraction

              Tn;Φ1 Ψ1 Φ2 Ψ2 Ψ̄1 Φ̄1 Ψ̄2 Φ̄2                                                             Let us figure out how to carry out the contraction between the
                                                                                                       fundamental tensors derived above. As an example, we con-
                       ∏ ∑                                                                               sider the contraction between Tn and Tn+1̂ , which reproduces
                 =                       T
                                             n;(i1 j1 )(i2 j2 )(i ′ j ′ )(i ′ j ′ )
                                                                               1 1    2 2                the hopping terms on the link (n, n + 1̂). When we use the
                            ν iν ,jν ,iν′ ,jν′
                                                                                                         fundamental tensor given in equation (6), the contraction is
                                                     j′     i′       j′   i′
                      × Φi11 Ψ1j1 Φi22 Ψj22 Ψ̄11 Φ̄11 Ψ̄22 Φ̄22 ,                             (15)       defined by
                                                                                                           ∑ˆ
with                                                                                                              Tn+1̂ Tn
                                                   ∑    ′   √ ∑ν (iν +jν +iν′ +jν′ )                             ∑
                                                     ν jν
  Tn;(i1 j1 )(i2 j2 )(i ′ j ′ )(i ′ j ′ ) = (−1)             t                                                =         Tn+1̂;(i1 j1 )(i2 j2 )(α1 β1 )(i ′ j ′ ) Tn;(α1 β1 )(k2 l2 )(k ′ l ′ )(k ′ l ′ )
                     1 1     2 2                                                                                                                        2 2                           1 1       2 2
                                                                                                                 α1 ,β1
                                                 j i +j ′ +j ′ +j j ′ +j ′ +i ′ j ′
                                        × (−1) 1 ( 2 1 2 ) 2 ( 1 2 ) 1 2                                            ˆ
                                          [                                                                      × Gn+1̂;(i1 j1 )(i2 j2 )(α1 β1 )(i ′ j ′ ) Gn;(α1 β1 )(k2 l2 )(k ′ l ′ )(k ′ l ′ ) .
                                        × −mδi1 +i2 +j1′ +j2′ ,0 δi1′ +i2′ +j1 +j2 ,0                                                                                    2 2                                  1 1   2 2
                                                                                   ]                                                                                                                                   (20)
                                        + δi1 +i2 +j1′ +j2′ ,1 δi1′ +i2′ +j1 +j2 ,1 . (16)
                                                                                                                     ´
As for equations (7) and (16) shows explicitly that the funda-                                           Note that          on the right-hand side means the integration
mental tensor in equation (15) is Grassmann-even. Compared                                               over the auxiliary Grassmann variables labeled by repeated
with the previous formulation based on [102], all the bits                                               Greek indices. Firstly, we integrate out (Φ̄1 (n + 1̂)Φ1 (n))α1
introduced by the Taylor expansion are summed within the                                                 and (Ψ̄1 (n)Ψ1 (n + 1̂))β1 originating from Gn , which results in
fundamental tensor Tn . Instead, we can identify the auxili-                                               ˆ
ary Grassmann variables as the indices of Tn . This is why                                                    Gn+1̂;(i1 j1 )(i2 j2 )(α1 β1 )(i ′ j ′ ) Gn;(α1 β1 )(k2 l2 )(k ′ l ′ )(k ′ l ′ )
                                                                                                                                                           2 2                              1 1        2 2
we have introduced the notation as in the left-hand side of
                                                                                                                                 (α1 +β1 )(i2′ +j2′ )+α1
equation (15). Following [125], we refer Tn in equation (15)                                                      = (−1)                                           Q(i1 j1 )(i2 j2 )(i ′ j ′ )(k2 l2 )(k ′ l ′ )(k ′ l ′ ) ,
                                                                                                                                                                                        2 2                   1 1   2 2
as the Grassmann tensor and T n in equation (16) as the coef-                                                                                                                                                          (21)
ficient tensor of Tn . The difference between the bosonic part
in equation (7) and the coefficient
                          ∑ ′        tensor in equation (16) is                                          with
just the sign factor (−1) ν jν in equation (16), which origin-
ates from the negative sign in the last exponential factor in                                            Q(i1 j1 )(i2 j2 )(i ′ j ′ )(k2 l2 )(k ′ l ′ )(k ′ l ′ )
                                                                                                                              2 2              1 1       2 2
equation (13). The graphical expression of equation (15) is                                                                                   j′  i′                 l′   k′    l′   k′
shown in figure 1. The path integral is reproduced from the                                                   = dΦi11 dΦ̄j11 dΨi22 dΨ̄j22 dΨ22 dΦ̄22 dΨk22 dΨ̄l22 dΨ11 dΦ̄11 dΨ22 dΦ̄22
fundamental tensor Tn in equation (15). Introducing the fol-                                                      (          )i (          )j (       )i (         )j (       )k (      )l
                                                                                                                × Φ̄1 Φ1 1 Ψ̄1 Ψ1 1 Φ̄2 Φ2 2 Ψ̄2 Ψ2 2 Φ̄2 Φ2 2 Ψ̄2 Ψ2 2                                                        .
lowing abbreviation,                                                                                                                                                                                                   (22)
                    ˆ       ˆ
                          = dΦ̄dΦ e−Φ̄Φ ,                 (17)                                           We have omitted the site dependence in the auxiliary
                               Φ̄,Φ                                                                      Grassmann variables in equation (22) because it can be read
                                                                                                         from the bit indices immediately. The sign in equation (21) is
and                                                                                                      then taken into account by contracting the bosonic parts. As a
                       ( ˆ                             ˆ                         )                       result of equation (20), we obtain
                        ∏
       gTr [ · ] =                                                                   [ · ],   (18)           ∑ˆ
                           n,ν     Φ̄ν (n),Φν (n)         Ψ̄ν (n),Ψν (n)
                                                                                                                     Tn+1̂ Tn = M(i1 j1 )(i2 j2 )(i ′ j ′ )(k2 l2 )(k ′ l ′ )(k ′ l ′ )
                                                                                                                                                                         2 2              1 1         2 2

Z in equation (2) is given by                                                                                                                     × Q(i1 j1 )(i2 j2 )(i ′ j ′ )(k2 l2 )(k ′ l ′ )(k ′ l ′ ) ,          (23)
                                                                                                                                                                               2 2              1 1         2 2
                                                 [               ]
                                                     ∏                                                   where
                                   Z = gTr                 Tn .                               (19)
                                                     n
                                                                                                             M(i1 j1 )(i2 j2 )(i ′ j ′ )(k2 l2 )(k ′ l ′ )(k ′ l ′ )
                                                                                                                                2 2               1 1       2 2
The symbol ‘gTr’ in equation (18) means the Grassmann                                                                 ∑                  (α1 +β1 )(i2′ +j2′ )+α1
tensor trace that is analogous to the tensor trace symbol ‘tTr’                                                =              (−1)                                   Tn+1̂;(i1 j1 )(i2 j2 )(α1 β1 )(i ′ j ′ )
                                                                                                                                                                                                                    2 2
                                                                                                                         α1 ,β1
common in the tensor network formulation in the spin systems.
   This formulation has been applied recently for the                                                                   × Tn;(α1 β1 )(k2 l2 )(k ′ l ′ )(k ′ l ′ ) .                                                    (24)
                                                                                                                                                          1 1      2 2
free Wilson and staggered fermions [125], Hubbard mod-
els [115, 116], Nf = 1, 2, 3 Gross–Neveu model [103, 126],                                               It is important to realize that the auxiliary Grassmann variables
Zn and U(1) gauge theories with Nf = 1, 2, 4 Wilson fermi-                                               should be integrated before the bosonic parts are contracted
ons [127], and several public codes for the Grassmann TRG                                                when we use the fundamental tensor given in equation (6).
methods [128, 129].                                                                                      Figure 2 illustrates the exact contraction in equation (20),

                                                                                                     6
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                          Topical Review




Figure 2. Graphical representation of the exact contraction between Tn+1̂ and Tn . (Top) Illustration of equation (20). Internal lines denote
the integration over the auxiliary
                               ´ Grassmann variables
                                                ´          and summation over α1 and β 1 . (Bottom) Illustration of equation (30). Internal lines
denote the weighted integrals Φ̄1 (n),Φ1 (n) and Ψ̄1 (n),Ψ1 (n) defined by equation (17). We have used the common Greek indices for the
contracted auxiliary Grassmann variables because of equation (32).




                                                                        7
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                                                                                 Topical Review



explicitly showing the site dependence of auxiliary Grassmann                                                                      should be understood as a result of the integrals on the auxili-
variables.                                                                                                                         ary Grassmann variables. From equation (17), we find that
   Note that the sign factor in equation (21) can be fur-                                                                                               ˆ
ther simplified by rearranging the Grassmann measures in                                                                                                       Θi Θ̄j = δij ,                 (32)
equation (22). Instead of equation (23), one can find                                                                                                                         Θ̄,Θ

    ∑ˆ                                                                                                                             where Θ and Θ̄ are the Grassmann variables. The integra-
            Tn+1̂ Tn = M(′i j )(i j )(k l )(k ′ l ′ )(i ′ j ′ )(k ′ l ′ )                                                          tion of the auxiliary variables naturally introduces contractions
                                                 1 1       2 2     2 2         1 1      2 2        2 2

                                          × Q(′i                                                                                   between the corresponding coefficient tensors in this formal-
                                                       1 j1 )(i2 j2 )(k2 l2 )   (k1′ l1′ )(k2′ l2′ )(i2′ j2′ ) ,        (25)
                                                                                                                                   ism. Figure 2 illustrates the exact contraction in the left-hand
                                                                                                                                   side of equation (30), without omitting the site dependence of
where                                                                                                                              auxiliary Grassmann variables.
                                                                                                                                      Rearranging the auxiliary Grassmann variables, the right-
M(′i1 j1 )(i2 j2 )(k2 l2 )(k ′ l ′ )(i ′ j ′ )(k ′ l ′ )                                                                           hand side of equation (30) can be
                            1 1       2 2       2 2
          ∑
   =                       α1
                  (−1) Tn+1̂;(i1 j1 )(i2 j2 )(α1 β1 )(i ′ j ′ )                                                                         ∑
          α1 ,β1
                                                                         2 2
                                                                                                                                             M(′i j )(i j )(k l )(k ′ l ′ )(k ′ l ′ )(i ′ j ′ ) Φi11 Ψj11 Φi22 Ψj22 Φk22 Ψl22
                                                                                                                                                 1 1   2 2   2 2   1 1       2 2       2 2

         × Tn;(α1 β1 )(k2 l2 )(k ′ l ′ )(k ′ l ′ ) ,                                                                    (26)                       l′    k′      l′     k′       j′     i′
                                         1 1         2 2                                                                                     × Ψ̄11 Φ̄11 Ψ̄22 Φ̄22 Ψ̄22 Φ̄22 ,                                                      (33)
Q(′i1 j1 )(i2 j2 )(k2 l2 )(k ′ l ′ )(k ′ l ′ )(i ′ j ′ )
                            1 1       2 2       2 2
                                                                                                                                   where
                                                                    l′         k′        l′      k′        j′      i′
     = dΦi11 dΦ̄j11 dΨi22 dΨ̄j22 dΨk22 dΨ̄l22 dΨ11 dΦ̄11 dΨ22 dΦ̄22 dΨ22 dΦ̄22
         (          )i (          )j (         )i (       )j (       )k (      )l                                                             M(′i                          (k1′ l1′ )(k2′ l2′ )(i2′ j2′ )
       × Φ̄1 Φ1 1 Ψ̄1 Ψ1 1 Φ̄2 Φ2 2 Ψ̄2 Ψ2 2 Φ̄2 Φ2 2 Ψ̄2 Ψ2 2 .                                                                                   1 j1 )(i2 j2 )(k2 l2 )
                                                                                                                                                        ∑
                                                                                                                        (27)                       =             (−1)
                                                                                                                                                                              α1 +β1
                                                                                                                                                                                         Tn+1̂;(i1 j1 )(i2 j2 )(α1 β1 )(i ′ j ′ )
                                                                                                                                                                                                                           2 2
                                                                                                                                                        α1 ,β1
The identity used here is
                                                                                                                                                        × Tn;(α1 β1 )(k2 l2 )(k ′ l ′ )(k ′ l ′ ) .                                 (34)
                                                                                                                                                                                             1 1     2 2
         i ′ +j ′ k +l +k ′ +l ′ +k ′ +l ′ i ′ +j ′ (α +β )
    (−1)( 2 2 )( 2 2 1 1 2 2 ) = (−1)( 2 2 ) 1 1 , (28)
                                                                                                                                   The same identity as in equations (28) and (29) has been
                                                                                                                                   utilized.
because the fundamental tensor Tn is Grassmann-even;
                                                                                                                                       Comparing equations (24) and (31), or equations (26)
Tn;(α1 β1 )(k2 l2 )(k1′ l1′ )(k2′ l2′ ) takes a non-zero value only if
                                                                                                                                   and (34), we see that the difference between the signs obtained
                                                                                                                                   in the two formulations is (−1)β1 . This can be explained by the
         α1 + β1 + k2 + l2 + k1′ + l1′ + k2′ + l2′ mod 2 = 0                                                            (29)       different ways to decompose the hopping terms employed in
                                                                                                                                   the two formulations. In equations (3) and (4), the last factors
is satisfied. This kind of identity is useful to simplify the sign                                                                 built just by auxiliary Grassmann variables inherit the struc-
factor arising in the exact contraction between the fundamental                                                                    ture of original hopping terms. On the other hand, there is no
tensors.                                                                                                                           such difference between equations (12) and (13). The corres-
    Let us consider the same contraction but with equation (15).                                                                                               ′
                                                                                                                                   ponding sign factor (−1)j1 has already been included in the
Since Tn+1̂ and Tn are connected via Φ1 (n), Φ̄1 (n), Ψ1 (n), and                                                                  coefficient tensor as in equation (16).
Ψ̄1 (n), they should be integrated. The exact contraction gives                                                                        If we performed the exact contractions using the methods
   ˆ                       ˆ                                                                                                       described so far, the path integral shown in equations (10)
                                                     Tn+1̂;Φ1 Ψ1 Φ2 Ψ2 Ψ̄1 (n)Φ̄1 (n)Ψ̄2 Φ̄2                                       and (19) could be obtained exactly. However, in practice, it
      Φ̄1 (n),Φ1 (n)           Ψ̄1 (n),Ψ1 (n)                                                                                      is impossible to keep carrying out the exact contractions when
     × Tn;Φ1 (n)Ψ1 (n)Φ2 Ψ2 Ψ̄1 Φ̄1 Ψ̄2 Φ̄2                                                                                        we consider the model in an arbitrarily large volume. This is
         ∑                                                                                  j′ i′                                  because the number of fundamental tensors in equations (10)
      =       M(i1 j1 )(i2 j2 )(i ′ j ′ )(k2 l2 )(k ′ l ′ )(k ′ l ′ ) Φi11 Ψj11 Φi22 Ψj22 Ψ̄22 Φ̄22                                and (19) is the same as the number of lattice sites in Λ.
                                               2 2               1 1       2 2

                           l′ k′ l′ k′                                                                                             Therefore, we need to consider performing the contraction not
             × Φk22 Ψl22 Ψ̄11 Φ̄11 Ψ̄22 Φ̄22 ,                                                                          (30)
                                                                                                                                   exactly but approximately. Several TRG schemes are reviewed
                                                                                                                                   in sections 2.5 and 2.6.
with

 M(i1 j1 )(i2 j2 )(i ′ j ′ )(k2 l2 )(k ′ l ′ )(k ′ l ′ )                                                                           2.3. Model-independent notation
                    2 2               1 1       2 2
          ∑                  (α1 +β1 )(i2′ +j2′ )+α1 +β1
   =              (−1)                                   Tn+1̂;(i1 j1 )(i2 j2 )(α1 β1 )(i ′ j ′ )                                  In this subsection, we introduce a new notation for the funda-
            α1 ,β1
                                                                                                                        2 2
                                                                                                                                   mental tensor
            × Tn;(α1 β1 )(k2 l2 )(k ′ l ′ )(k ′ l ′ ) .                                                                 (31)                                          Tn = Tn;xtx ′ t ′ Gn;xtx ′ t ′ ,                              (35)
                                               1 1         2 2



The right-hand side of equation (30) defines a new Grassmann                                                                       instead of equation (6). The variable x is defined by two bits
tensor. Here, the contraction between the coefficient tensors                                                                      via x = (i1 j1 ) and t, x′ , and t′ are defined similarly. Gn;xtx ′ t ′ is

                                                                                                                               8
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                                                                     Topical Review



defined by the right-hand side of equation (8). The expression                                                         Similarly, we introduce the following expression for
in equation (35) may seem fine because it does not depend on                                                        equation (15),
the details of our model, the number of components in the ori-
                                                                                                                                                                ∑                                  ′    ′
ginal fermionic field or the structure of hopping terms, and it                                                                        Tn;XTX̄T̄ =                           Tn;xtx ′ t ′ Xx T t X̄x T̄ t .                (42)
only depends on the lattice geometry. However, this notation                                                                                                  x,t,x ′ ,t ′
can be problematic when considering their exact contractions.
As we have observed in equation (24), or equation (26), the                                                         As in equation (35), x is defined by two bits via x = (i1 j1 )
sign factor from the auxiliary Grassmann integrals does dis-                                                        and t, x′ , and t′ are defined similarly. X can be regarded
tinguish the forward and backward hopping terms in the for-                                                         as a two-component auxiliary Grassmann variables via X =
mulation in [102]. One of the ways to resolve this issue is to                                                      (Φ1 , Ψ1 ) and X x as Xx = Φi11 Ψj11 . X̄ can also be regarded like
modify equation (4) as                                                                                                                                                                                                j′    i′
                                                                                                                    X̄ = (Φ̄1 , Ψ̄1 ), but X̄x should be understood as X̄x = Ψ̄11 Φ̄11 .
                                                                                                                                        ′
                                                                                                                    T (T̄) and T t (T̄ t ) are defined in the same way. Equation (30)
e tψ̄(n)ψ(n+ν̂)
                                                                                                                    can be equivalently expressed as
              ∑
              1
                (√                )jν (n) (√                         )jν (n)
      =           tψ̄ (n) dΨ̄ν (n)           tψ (n + ν̂) dΨν (n + ν̂)                                                      ˆ
          jν (n)=0                                                                                                                  Tn+1̂;XT1 Θ̄T̄1 Tn;ΘT2 X̄T̄2
               (                                )jν (n)                                                                      Θ̄,Θ
                                                                                                                                           ∑
          × −Ψν (n + ν̂) Ψ̄ν (n)                            ,                                           (36)                                                                            t′   t′    ′   t′     t′
                                                                                                                               =                              Mxt1 t2 x ′ t1′ t2′ Xx T11 T22 X̄x T̄22 T̄11 ,               (43)
                                                                                                                                    x,t1 ,t2 ,x ′ ,t1′ ,t2′
and absorbing the extra factor (−1)jν (n) into the bosonic tensor
in equation (7). Note that this modification makes the bosonic
                                                                                                                     where the coefficient tensor Mxt1 t2 x ′ t1′ t2′ is defined in the exactly
tensor in equation (7) exactly the same as the coefficient tensor
                                                                                                                    same way with equation (40) using the Grassmann parity
in equation (16). From now on, we identify the expression
                                                                                                                     function defined in equation (38). Note that Θ and Θ̄ in
in equation (35) with this modification. The exact contraction
demonstrated in section 2.2 is then denoted by                                                                      ´equation (43) are two-component Grassmann variables and
                                                                                                                      Θ̄,Θ
                                                                                                                           in the left-hand side is defined by
∑ˆ                             ∑                                          ˆ
                                                                                                                                             ˆ
              Tn+1̂ Tn =              Tn+1̂;xt1 αt ′ Tn;αt2 x ′ t2′           Gn+1̂;xt1 αt ′ Gn;αt2 x ′ t2′ .                                                   ∏ˆ
                                 α
                                                    1                                       1
                                                                                                                                                         =                    dΘ̄i dΘi e−Θ̄i Θi ,                          (44)
                                                                                                        (37)                                     Θ̄,Θ             i


Now, we introduce the Grassmann parity function f x for x =                                                         which is a natural extension of equation (17). We can also see
(i1 j1 ) such that                                                                                                  that f x in equation (38) counts the Grassmann parity of X x .
                                                                                                                    The right-hand side of equation (43) is the Grassmann tensor
                                      fx = i1 + j1 mod 2,                                               (38)        that can be written as MXT1 T2 X̄T̄2 T̄1 in our notation. Extend the
                                                                                                                    definition of the Grassmann parity function in equation (38),
and f t , fx ′ , and ft ′ are done in the same way. Using these parity                                              and the notation in equation (42) immediately allows us to deal
functions, equation (37) is evaluated as                                                                            with the case where the original fermionic model is described
                                                                                                                    by the multi-component Grassmann variables.
                        ∑ˆ                                                                                              So far, we have confirmed that the bosonic tensor derived
                                   Tn+1̂ Tn = Mxt1 t2 x ′ t1′ t2′ Qxt1 t2 x ′ t2′ t1′ ,                 (39)        by [102] and the coefficient tensor by [125] can be equivalent
                                                                                                                    if we slightly modify the formulation of [102]. It has also been
where                                                                                                               confirmed that the exact contraction can be described similarly
                                                                                                                    using the same Grassmann parity function with either formu-
                                           ∑
                   Mxt1 t2 x ′ t1′ t2′ =
                                                        f
                                               (−1) α Tn+1̂;xt1 αt ′ Tn;αt2 x ′ t2′ ,                   (40)        lation. Therefore, these two formulations are fully equivalent
                                           α
                                                                              1
                                                                                                                    in terms of memory usage and computational time for con-
                                                                                                                    tractions. Figure 3 graphically shows the Grassmann tensor
and                                                                                                                 network representation of equation (2). Figure 3(A) can be
                                                                                                                    interpreted in two ways following the notation rule in figures 1
Qxt1 t2 x ′ t2′ t1′                                                                                                 and 2. Figure 3(B) assumes the model-independent notation.
                                                                l′   k′       l′     k′     j′     i′               The following demonstration always assumes the formalism
     = dΦi11 dΦ̄j11 dΨi22 dΨ̄j22 dΨk22 dΨ̄l22 dΨ11 dΦ̄11 dΨ22 dΦ̄22 dΨ22 dΦ̄22                                      in [125] with the notation introduced in this section 2.3.
         (          )i (          )j (         )i (       )j (       )k (      )l
       × Φ̄1 Φ1 1 Ψ1 Ψ̄1 1 Φ̄2 Φ2 2 Ψ2 Ψ̄2 2 Φ̄2 Φ2 2 Ψ2 Ψ̄2 2 .
                                                                                                        (41)
                                                                                                                    2.4. Extension to lattice gauge theories

It should be emphasized that the exact contraction between the                                                      Current Grassmann tensor network formulations are easily
fundamental tensors expressed as in equation (35) is straight-                                                      combined with the usual tensor network formulations for pure
forwardly generalized to other models just modifying the                                                            gauge or bosonic theories. See [67] and several reviews [74,
definition of the Grassmann parity function in equation (38).                                                       130–132] for more details on how to construct fundamental

                                                                                                                9
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                       Topical Review




Figure 3. Diagrammatic representation of Grassmann tensor network on a two-dimensional square lattice. (A) Grassmann tensor network
representation of equation (2). Interpretation of each line is given in figure 1. Background dotted lines denote the real-space lattice. (B)
Model-independent description of two-dimensional Grassmann tensor network. The shape of each fundamental tensor is determined only by
the lattice geometry.


tensors in these theories. Here, let us consider the two-                           We parameterize Uν (n) by an integer qν (n) mod N via
dimensional ZN gauge theory on a periodic square lattice                            Uν (n) = exp [2πiqν (n)/N]. Following section 2.1, one can
defined by                                                                          immediately obtain a Grassmann tensor with some integer
                                                                                    indices. We begin with decomposing the hopping terms via
                                    S = Sg + Sf ,                       (45)
                                                                                                                  ∗

as an example. We assume the standard Wilson action,                                    eην (n)χ̄(n+ν̂)Uν (n)χ(n)/2
                                                                                               ˆ
              ∑ [         (      ) (        )       ]                                      = dΦ̄ν (n) dΦν (n) e−Φ̄ν (n)Φν (n)
 Sg = −β       ℜ U1 (n) U2 n + 1̂ U∗1 n + 2̂ U∗2 (n) , (46)                                                                    √              ∗
                                                                                                                                                             √
              n∈Λ                                                                               × eχ̄(n+ν̂)Φ̄ν (n)/             2
                                                                                                                                    eην (n)Uν (n)χ(n)Φν (n)/  2
                                                                                                                                                                   ,     (49)
                                                                                        e−ην (n)χ̄(n)Uν (n)χ(n+ν̂)/2
and the staggered fermion action,                                                             ˆ
                                                                                          = dΨ̄ν (n) dΨν (n) e−Ψ̄ν (n)Ψν (n)
        ∑ ην (n) [                                                        ]                                                             √                      √
 Sf =                  χ̄ (n) Uν (n) χ (n + ν̂) − χ̄ (n + ν̂) U∗ν (n) χ (n)
               2                                                                                × eην (n)Uν (n)χ̄(n)Ψν (n)/              2
                                                                                                                                              eχ(n+ν̂)Ψ̄ν (n)/  2
                                                                                                                                                                    .    (50)
        n,ν
              ∑
        +m         χ̄ (n) χ (n) .                                       (47)
              n
                                                                                    Integrating out the original staggered fields at each site, one
                                                                                    obtains
β and m denote the inverse gauge coupling and mass, respect-
ively. The staggered sign function ην (n) is defined by η1 (n) =
                                                                                        ˆ                             ∏                   √                          √
1 and η2 (n) = (−1)n1 at each site n = (n1 , n2 ). Our goal is to
derive the Grassmann tensor network representation for the                                  dχdχ̄ e−mχ̄χ                    eχ̄Φ̄ν (n−ν̂)/    2 ην (n)Uν (n)χ̄Ψν (n)/ 2
                                                                                                                                               e
                                                                                                                      ν
following path integral,                                                                                  ∗
                                                                                                                              √                    √
                                                                                         × eην (n)Uν (n)χΦν (n)/ 2 eχΨ̄ν (n−ν̂)/ 2
                    ˆ ∏                  ∏                                                      ∑ (f )                                           ′     ′

              Z=               dUν (n)        dχ (n) dχ̄ (n) e−S .      (48)              =            Txf tf x ′ t ′ ,q2 (n)q1 (n) Xxf T tf X̄xf T̄ tf ,                (51)
                                                                                                                      f f
                         n,ν              n                                                    xf ,tf ,xf′ ,tf′




                                                                               10
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                                                                                                Topical Review



where we have employed the model-independent notation in                                                                               MXT1 T2 X̄T̄2 T̄1 ,xg tg 1 tg 2 xg′ tg′ tg′
section 2.3 and the coefficient tensor is given by                                                                                                    ∑                       1 2

                                                                                                                                        =                                 Mxf tf1 tf2 xf′ tf′ 1 ,tf′ 2 ,xg tg 1 tg 2 xg′ tg′                   t′
                                                                                                                                                                                                                                              1 g2
  (f)
Txf tf x ′ t ′ ,q2 (n)q1 (n)                                                                                                                       xf ,tf1 ,tf2 ,xf′ ,tf′ 1 ,tf′ 2
        f f
                                                                                                                                                                                     ′     t′     t′
     = T(i
              (f)
                                                                                                                                                   × Xxf T1tf1 T2tf2 X̄xf T̄2f 2 T̄1f 1 ,                                                                 (58)
              (i1′ j1′ )(i2′ j2′ ),q2 (n)q1 (n)
               1 j1 )(i2 j2 )

      √ − ∑ν (iν +jν +iν′ +jν′ )                    ′   ′         ′   ′    ′ ′
     = 2                               (−1)j1 (i2 +j1 +j2 )+j2 ( j1 +j2 )+i1 j2 +n1 (i2 +j2 )                           whose coefficient tensor is defined by
                    2π i
                           ∑
         ×e          N         ν qν (n)( jν −iν )                                                                                                                                        ∑
          [                                                                                                             Mxf tf1 tf2 xf′ tf′    ,t ′ ,x t t x ′ t ′ t ′          =                 (−1) fαf Tn+1̂;xf tf1 αf t ′                                  ′
                                                                                                                                              1 f 2 g g1 g2 g g g                                                                              f 1 ,xg tg 1 αg tg 1
         × −mδi1 +i2 +j1′ +j2′ ,0 δi1′ +i2′ +j1 +j2 ,0                                                                                                                 1    2
                                                                                                                                                                                         αf ,αg
                                                    ]
         + δi1 +i2 +j1′ +j2′ ,1 δi1′ +i2′ +j1 +j2 ,1 .                                                      (52)                                                                         × Tn;αf tf2 xf′ tf′       ,αg tg 2 xg′ tg′       .               (59)
                                                                                                                                                                                                               2                      2


                                                 ′
Note that Xxf , T tf , X̄xf , and T̄ tf in the right-hand side of
                                                                    ′                                                   The Grassmann parity function fαf is given in the same way
equation (51) have been defined in the same way with                                                                    with equation (38). Introducing a super index p = (pf , pg ) for
equation (42). We also introduce a four-leg tensor to describe                                                          p = x, t, x ′ , t ′ , equation (59) reads
the plaquette interaction term in e−Sg via                                                                                                                           ∑                      Fα
                                                                                                                                        Mxt1 t2 x ′ t1′ t2′ =                   (−1)              Tn+1̂;xt1 αt ′ Tn;αt2 x ′ t2′ ,                         (60)
                                                                                                                                                                                                                      1
                     (g)                                                                                                                                                α
                    Tq n+1̂ q n+2̂ q (n)q (n)
                      2(   ) 1( ) 2 1
                                 [     { (                (     )
                                         2π                                                                             where we have defined a new parity function for the super
                           = exp β cos       q1 (n) + q2 n + 1̂                                                         index α = (αf , αg ) by
                                          N
                                   (     )         ) }]
                              − q1 n + 2̂ − q2 (n)      ,                                                   (53)                                                                F α = fα f .                                                              (61)

following [133]. Regarding qν ’s as tensor subscripts, we now                                                           The function F α tells us the Grassmann parity of the super
define the fundamental tensor associated at the lattice site n as                                                       index α: when Fα = 0 (1), the super index α is describing the
                             ∑                                                                                          Grassmann-even (odd) contribution.
                                                                            ′     ′
   Tn;XTX̄T̄,xg tg xg′ tg′ =   Tn;xf tf xf′ tf′ ,xg tg xg′ tg′ Xxf T tf X̄xf T̄ tf , (54)                                   Therefore, we reach the important conclusion that the
                                         xf ,tf ,xf′ ,tf′                                                               structure of the Grassmann tensor network and the contrac-
                                                                                                                        tion rule are not affected by the gauge fields. The resulting
where                                                                                                                   coefficient tensor in equation (60) has the same expression
                                                                                                                        as equation (40), which was the coefficient tensor for the
                                                                (f )                   (g)
                           Tn;xf tf xf′ tf′ ,xg tg xg′ tg′ = Txf tf x ′ t ′ ,xg tg · Txg tg xg′ tg′ .       (55)        pure fermionic model in equation (1). The exact contractions
                                                                        f f
                                                                                                                        among the Grassmann tensors result in those among the coef-
The path integral in equation (48) is now represented by                                                                ficient tensors, with or without the lattice gauge fields. One
                               [       ]                                                                                can use the tensor network diagram in figure 3(B) instead
                                 ∏                                                                                      of figure 4(B) to represent equation (56). Although we have
                      Z = gTr       Tn .                 (56)                                                           assumed ZN as a gauge group for simplicity, the above con-
                                                                   n
                                                                                                                        struction can be combined with other tensor network formu-
Here, ‘gTr’ stands not only for the integrations over all aux-                                                          lations for various lattice gauge theories including the non-
iliary Grassmann variables but also for the summations over                                                             Abelian fields [67, 134–139].
all integers corresponding to the link variables. Analogously,
we refer Tn and T n in equation (54) to the Grassmann tensor                                                            2.5. Approximate contraction by the Levin–Nave TRG
and its coefficient tensor, respectively. See figure 4 for the dia-
grammatic explanation.                                                                                                  The usual TRG algorithms are intended to be applied to the
    Let us briefly see how the exact contraction is carried out                                                         classical systems without any Grassmann variables. Still, any
between these fundamental tensors. As in section 2.2, we                                                                TRG algorithm can be utilized to evaluate the Grassmann
consider the contraction between the fundamental tensors at                                                             path integrals. As we have seen, the contractions between
n + 1̂ and n for demonstration. Since a link variable is shared                                                         the Grassmann tensors always result in those between the
between these two fundamental tensors, the exact contraction                                                            corresponding coefficient tensors with some sign factors
should be                                                                                                               arising from the integrals over the auxiliary Grassmann vari-
      ∑ˆ
                                                                                                                        ables. Therefore, we can use the TRG algorithms to carry
                Tn+1̂;XT1 Θ̄T̄1 ,xg tg αg t ′ Tn;ΘT2 X̄T̄2 ,αg tg 2 xg′ tg′ . (57)                                      out the contractions among the coefficient tensors approx-
                                                            1      g1                                   2
          αg           Θ̄,Θ                                                                                             imately. In this subsection, we demonstrate how to extend
                                                                                                                        the Levin–Nave TRG [27] for the Grassmann tensor net-
One will immediately find that the above contraction results in                                                         works. The case of the HOTRG [32] is discussed in the next
a new Grassmann tensor such that,                                                                                       subsection.

                                                                                                                   11
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                           Topical Review




Figure 4. Diagrammatic representation of Grassmann tensor network formulation for two-dimensional lattice gauge theories. Background
dotted lines show a real-space square lattice. (A) Structure of the fundamental tensor defined in equation (54). The blue symbol is located on
the site and denotes the right-hand side of equation (51). Note that the coefficient tensor depends on n1 as shown in equation (52). The red
                                                                      (g)
symbol is located on the plaquette and denotes the four-leg tensor Txg tg x ′ t ′ in equation (53). Each diamond shows the link variable. (B)
                                                                           g g
Grassmann tensor network in equation (56). Each internal line represents the contraction between auxiliary Grassmann variables and
integration over the shared link variable.


   Let us begin with reviewing the original Levin–Nave TRG.                                which approximates the original Z as
The algorithm aims to evaluate the partition function or path
                                                                                                                            [                 ]
integral represented by the tensor network                                                                                        ∏
                              [       ]                                                                           Z ≃ tTr                  Tn ′ .             (66)
                                ∏                                                                                               n ′ ∈Λ ′
                      Z = tTr       Tn ,                 (62)
                                                   n∈Λ
                                                                                           Therefore, Tn ′ is a new fundamental tensor defined on a
where we assume the two-dimensional periodic square lat-                                   site n′ in the coarse-grained lattice Λ ′ . Since equation (66)
tice Λ and T n is a four-leg tensor defined on the site n. The                             describes the tensor network whose geometry is the same as
algorithm employs the singular value decomposition (SVD)                                   that in equation (62), we can easily repeat the above decima-
to decimate the four-leg tensor T n into three-leg tensors:                                tion procedure. Repeating this procedure N times, 2N original
                                                                                           fundamental tensors are approximately contracted. Thanks to
                                            D∑
                                             LNTRG                                         this property, the algorithm allows us to evaluate the parti-
                      T   n;xtx ′ t ′   ≃            An;xta Bn;ax ′ t ′ ,    (63)          tion function in the( thermodynamic
                                                                                                                        )         limit. The Levin–Nave
                                                                                                                                                (        )
                                             a=1                                           TRG requires the O D6LNTRG complexity and the O D4LNTRG
                                            D∑                                             memory    cost.
                                             LNTRG
                                                                                              (        ) Note that the cost can(further be) reduced to the
                      Tn;xtx ′ t ′ ≃                 Cn;xt ′ a Dn;ax ′ t .   (64)          O D5LNTRG complexity and the O D3LNTRG memory using
                                             a=1                                           the randomized SVD [140]. The algorithm is graphically sum-
                                                                                           marized in figure 5.
Each three-leg tensor is defined as a unitary matrix multi-
                                                                                               We now need to define the SVD for the Grassmann tensor
plied by the square root of its singular value. In the above
                                                                                           to extend the algorithm for the Grassmann tensor network.
expressions, we have assumed that the singular values are
                                                                                           Suppose OΦΨ is a Grassmann tensor defined via
in descending order. The truncation parameter DLNTRG is
called the bond dimension. Since we have used the SVD,                                                                      ∑
equations (63) and (64) give the best approximation in terms                                                      OΦΨ =               Oij Φi Ψj ,             (67)
of the Frobenious norm of T n under the fixed bond dimension.                                                                   i,j

Then, we define a new four-leg tensor via
                    ∑                                                                      where we have assumed the notation in section 2.3, say i =
  Tn ′ ;xtx ′ t ′ =   Cn+2̂;x2 t1 x An;x1 t1 t Dn+1̂;x ′ x1 t2 Bn+1̂+2̂;t ′ x2 t2 ,        (i1 · · · im ) and j = (j1 · · · jn ) with Φi = Φi11 · · · Φimm and Ψj =
                 x1 ,x2 ,t1 ,t2                                                            Ψj11 · · · Ψjnn . We also assume that O is Grassmann-even, in other
                                                                             (65)          words, the matrix elements are zero unless the sum of all the

                                                                                      12
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                     Topical Review




Figure 5. Schematic illustration of the Levin–Nave TRG algorithm. Background dotted lines show a real-space square lattice. (A) Initial
tensor network on the lattice. (B) Two kinds of SVD as shown in equations (63) and (64). (C) New tensor network by contracting four
three-leg tensors.


indices (i and j) is even. Since the coefficient tensor is normal,                      introducing the low-rank approximation based on the SVD.
we can apply the SVD for Oij ,                                                          In the following, we use the approximation

                     min(2m ,2n )
                                                                                                                    ˆ    DLNTRG
                        ∑                             ∑
             Oij =                  Uia σa Vaj†   =         Aia δab Bbj .   (68)                         OΦΨ ≃                    AΦΞ BΞ̄Ψ ,           (70)
                                                                                                                        Ξ̄,Ξ
                        a=1                           a,b

                                        √            √                                  which means that the SVD of the coefficient tensor is truncated
In the second equality, we set A = U σ and B = σV † .                                   up to the bond dimension DLNTRG as
Recalling the identity in equation (32), we can express δ ab
by introducing new auxiliary Grassmann variables. Plugging                                                          D∑
                                                                                                                     LNTRG
equation (68) into equation (67), we now have the SVD for the                                               Oij ≃              Aia δab Bbj .           (71)
Grassmann tensor O as                                                                                               a,b=1
                            ˆ
                    OΦΨ =        AΦΞ BΞ̄Ψ ,             (69)                            We will make some practical remarks on the SVD of the
                                        Ξ̄,Ξ                                            Grassmann tensor in section 2.7
                                                                                          We now extend the algorithm for evaluating
where Ξ and Ξ̄ are the min(m, n)-component auxiliary
Grassmann variables introduced via equation (32). The                                                                      [          ]
                                                                                                                               ∏
Grassmann tensors A and B in the right-hand side have A and                                                   Z = gTr                Tn ,              (72)
B in equation (68) as their coefficient tensors. Since we have                                                                 n∈Λ
assumed that O is Grassmann-even, the resulting A and B are
also Grassmann-even, as we will see in section 2.7. It is no                            with the Grassmann tensor Tn;XTX̄T̄ . We can formally
exaggeration to say that the SVD for a Grassmann tensor is                              write the SVD for the Grassmann tensor corresponding to
the SVD for its coefficient tensor. There is no difficulty in                           equations (63) and (64) as



                                                                                   13
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                                                                       Topical Review



                                                                                                                2.6. Approximate contraction by the HOTRG

                                                                                                                Let us next focus on the HOTRG algorithm that applies to any
                                                                                                                       ( −1 system.
                                                                                                                d-dimensional   )                                  ( 2d scales
                                                                                                                                        The computational complexity        )
                                                                                                                with O D4dHOTRG and the memory cost does with O DHOTRG ,
                                                                                                                where DHOTRG is the bond dimension in the HOTRG. For
                                                                                                                simplicity, we consider the two-dimensional tensor network
                                                                                                                on a periodic square lattice again. The original algorithm
                                                                                                                aims to evaluate equation (62) not decimating each local
                                                                                                                fundamental tensor but inserting the projectors that describe
                                                                                                                the coarse-graining transformation in the tensor-network lan-
                                                                                                                guage. Unlike the Levin–Nave TRG, the HOTRG performs the
                                                                                                                contraction between two adjacent fundamental tensors along
                                                                                                                each direction sequentially.
                                                                                                                   We begin with reviewing the normal HOTRG. As an
                                                                                                                example, we consider the contraction between Tn+1̂ and T n .
                                                                                                                Introducing
                                                                                                                                                           ∑
                                                                                                                               Mxt1 t2 x ′ t1′ t2′ =                Tn+1̂;xt1 αt ′ Tn;αt2 x ′ t2′ ,                      (77)
                                                                                                                                                                                         1
                                                                                                                                                             α
      Figure 6. Diagrammatic representation of equation (75).

                                                                                                                the HOTRG provides the coarse-graining transformation such
                                        ˆ   DLNTRG                                                              that
                  Tn;XTX̄T̄ ≃                             An;XT Ξ Bn;Ξ̄X̄T̄ ,                       (73)
                                          Ξ̄,Ξ                                                                                                         ∑
                                         ˆ DLNTRG                                                                           Tn ′ ;xtx ′ t ′ =                       Pt1 t2 t Mxt1 t2 x ′ t1′ t2′ Qt ′ t1′ t2′ ,          (78)
                  Tn;XTX̄T̄ ≃                             Cn;XT̄ Ξ Dn;Ξ̄X̄T .                       (74)                                         t1 ,t2 ,t1′ ,t2′
                                            Ξ̄,Ξ
                                                                                                                where the three-leg tensors P and Q are projectors to decimate
When Tn is Grassmann-even, then A, B, C, and D are                                                              the degrees of freedom. They play roles to map the original
also Grassmann-even. Therefore, we can easily define a new                                                                              ( ′) ( ′)                           ′
                                                                                                                degrees of freedom (t1 t2 ) to the coarse-grained one t( ) ,
Grassmann tensor from the contractions among them,
                                                                                                                whose size is restricted by DHOTRG , the bond dimension in the
                      ˆ           ˆ             ˆ         ˆ                                                     algorithm. We can also regard that the HOTRG is inserting the
     Tn ′ ;XTX̄T̄ =                                                     Cn+2̂;X2 T̄1 X An;X1 T1 T               following four-leg tensor,
                        X̄1 ,X1       X̄2 ,X2   T̄1 ,T1       T̄2 ,T2
                      × Dn+1̂;X̄X̄1 T2 Bn+1̂+2̂;T̄X̄2 T̄2 .                                         (75)                                                        D∑
                                                                                                                                                                 HOTRG

                                                                                                                                        Wt1 t2 t1′ t2′ =                      Pt1 t2 t Qtt1′ t2′ ,                       (79)
The above equation is analogous to equation (65) in the usual                                                                                                       t=1
Levin–Nave TRG. Equation (75) is diagrammatically shown
in figure 6. The new Grassmann tensor approximates the path                                                     into the tensor network. With sufficiently large DHOTRG ,
integral by                                                                                                     Wt1 t2 t1′ t2′ should be equivalent to δt1 ,t1′ δt2 ,t2′ and the algorithm
                                                                                                                gives the exact contraction. Tn ′ generates the approximated
                                                [                       ]
                                                     ∏                                                          tensor network representation of the partition function. The
                             Z ≃ gTr                           Tn ′ .                               (76)        coarse-graining transformation in the HOTRG can be easily
                                                    n ′ ∈Λ ′                                                    repeated. As in the case of the Levin–Nave TRG, the HOTRG
                                                                                                                allows us to contract 2N original fundamental tensors approx-
The resulting Grassmann tensor network in equation (76) is                                                      imately, just in N times of iteration. Figure 7 graphically
identical to the previous one in equation (72), including the                                                   demonstrates the procedure of the algorithm.
ordering of the Grassmann measures in gTr. Therefore, one                                                          There are several ways to determine P and Q. Here, we just
can easily repeat the decimation procedure toward the ther-                                                     follow the original proposal in [32] assuming the translational
modynamic limit as in the usual Levin–Nave TRG. We can                                                          symmetry for the tensor network in equation (62). We consider
understand that the SVD of the Grassmann tensor introduces                                                      the following two reduced density matrices,
new auxiliary Grassmann variables on the coarse-grained lat-
tice Λ ′ . Therefore, the algorithm for the Grassmann tensor net-                                                                                         ∑
                                                                                                                             ρ(t1 t2 )(t̃1 t̃2 ) =                       Mxt1 t2 x ′ t1′ t2′ M∗xt̃1 t̃2 x ′ t ′ t ′ ,    (80)
work perfectly corresponds to that for the normal tensor net-                                                                                                                                                 1 2
                                                                                                                                                       x,x ′ ,t1′ ,t2′
work as shown in figure 5. In terms of coefficient tensors, sev-                                                                                           ∑
eral sign factors should be included in equations (74) and (75)                                                            ρ(t ′ t ′ )(t̃ ′ t̃ ′ ) =                     Mxt1 t2 x ′ t1′ t2′ M∗xt1 t2 x ′ t̃ ′ t̃ ′ .    (81)
                                                                                                                                1 2     1 2                                                                   1 2
as we will see in section 2.7.                                                                                                                         x,t1 ,t2   ,x ′




                                                                                                           14
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                                                      Topical Review




Figure 7. Schematic illustration of the HOTRG algorithm. (A) Initial tensor network on a square lattice. (B) Inset projectors into the
network. Red and blue symbols show P and Q in equation (78), respectively. (C) New tensor network by contracting adjacent two
fundamental tensors with two projectors.


These matrices can be decomposed as                                                                    where the coefficient tensor of M has been already derived
                                                                                                       in equation (40) or equation (60). The HOTRG for the
                                               ∑
                      ρ(t1 t2 )(t̃1 t̃2 ) =               U(t1 t2 )i λi U†i(t̃1 t̃2 ) ,    (82)        Grassmann tensor network should provide us with the coarse-
                                                    i                                                  graining transformation such as
                                               ∑                                                                        ˆ         ˆ         ˆ           ˆ
                 ρ(′ t ′ t ′ )(t̃ ′ t̃ ′ ) =            V(t ′ t ′ )i λi′ Vi† t̃ ′ t̃ ′ ,   (83)
                      1 2       1 2                        1 2             ( 1 2)                      Tn ′ ;XTX̄T̄ =                                                   PT̄2 T̄1 T MXT1 T2 X̄T̄ ′ T̄ ′ QT̄T ′ T ′ ,
                                                i                                                                                                                                                2   1       1   2
                                                                                                                        T̄1 ,T1   T̄2 ,T2   T̄1′ ,T1′       T̄2′ ,T2′
                                                                                                                                                                                                             (86)
where U and V are unitary matrices and λ and λ ′ denote the
singular-value matrices with descending order. Defining the                                            where P and Q are the Grassmann projectors defined by
following quantities,                                                                                                            ∑
                                                                                                                    PT̄2 T̄1 T =   Pt1 t2 t T̄2t2 T̄1t1 T t , (87)
                                 ′                  ∑            ( ′)                                                                       t1 ,t2 ,t
                               ϵ( ) =                           λi ,                       (84)                                                 ∑                             ′   ′ ′   ′ ′
                                                                                                                                                                                   t     t
                                               i>DHOTRG                                                                     QT̄T1′ T2′ =                    Qt ′ t1′ t2′ T̄ t T1 1 T2 2 .                    (88)
                                                                                                                                            t1′ ,t2′ ,t ′
we choose Pt1 t2 t = U∗(t1 t2 )t and Qt ′ t1′ t2′ = U(t1′ t2′ )t ′ if ϵ < ϵ ′ and
choose Pt1 t2 t = V(t1 t2 )t and Qt ′ t1′ t2′ = V∗(t ′ t ′ )t ′ if ϵ > ϵ ′ . This                      In analogy with equation (79), we can identify that the
                                                    1 2                                                algorithm inserts
procedure is equivalent to the higher-order SVD (HOSVD),
which is a tensorial extension of the SVD as explained in [32].                                                                        ˆ DHOTRG
See references [141–143] for other ways to derive the optimal                                                       WT̄2 T̄1 T1′ T2′ =          PT̄2 T̄1 T QT̄T1′ T2′ , (89)
                                                                                                                                                    T̄,T
P and Q without assuming the translational invariance on the
tensor network.                                                                                        into the Grassmann tensor network. P and Q, or P and Q
   The HOTRG algorithm is straightforwardly extended to                                                in other words, are determined via the same procedure in
evaluate the Grassmann tensor network in equation (72).                                                the normal HOTRG. Firstly, we define the conjugation of the
Equation (77) is correspondingly denoted by                                                            Grassmann tensor by
                                        ˆ                                                                                  †
                                                                                                                                        ∑
                                                                                                                                 ∗
            MXT1 T2 X̄T̄2 T̄1 =                     Tn+1̂;XT1 Θ̄T̄1 Tn;ΘT2 X̄T̄2 ,         (85)                      (OΦΨ ) = OΨ̄  Φ̄ =   O∗ij Ψ̄j Φ̄i ,     (90)
                                           Θ̄,Θ                                                                                                                         i,j

                                                                                                  15
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                                                                          Topical Review



where we have assumed that OΦΨ is defined by equation (67).                                                                      Grassmann-even. Therefore, we can always convert Oij into a
Note that for Φi = Φi11 · · · Φimm , we define Φ̄i in equation (90)                                                              block-diagonalized form by the matrix elementary operations;
as Φ̄i = Φ̄imm · · · Φ̄i11 . Ψ̄j is defined in the same way. Then, we                                                            regarding i as a row and j as a column, Oij can be
define the Grassmann reduced density matrices as
                       ˆ    ˆ           ˆ               ˆ
                                                                                                                                                                even j : odd ]
                                                                                                                                                         [ j : (even)
                                                                                                                                             O=             O            0                            i : even ,            (93)
ϱT1 T2 T̄2 T̄1 =                                                        MXT1 T2 X̄ ′ T̄ ′ T̄ ′ M∗T ′ T ′ X ′ T̄2 T̄1 X̄ ,
                       X̄,X X ′ ,X̄ ′       T2′ ,T̄2′       T1′ ,T̄1′
                                                                                          2   1       1   2                                                      0    O(odd)                           i : odd
                                                                                                                 (91)
                       ˆ    ˆ           ˆ               ˆ                                                                        where ‘i : even (odd)’ means the row index i such that Φi
ϱT̄′ 2 T̄1 T1 T2   =                                                    MXT ′ T ′ X̄ ′ T̄2 T̄1 M∗T1 T2 X ′ T̄ ′ T̄ ′ X̄ ,        becomes Grassmann-even (odd). Similarly, ‘j : even (odd)’
                                                                           1 2
                       X̄,X T̄1′ ,T1′   T̄2′ ,T2′
                                                                                                             2 1
                                                        X ′ ,X̄ ′                                                                means the column index j such that Ψj becomes Grassmann-
                                                                                                                 (92)            even (odd). Using this basis, the SVD of O gives us
                                                                                                                                                            ∑
where the ordering of the Grassmann measures are defined                                                                                              Oij =      Uik σkl Vlj† ,       (94)
so that the coefficient tensors of ϱT1 T2 T̄2 T̄1 and ϱT̄′ 2 T̄1 T1 T2 are                                                                                                       k,l
given by ρ(t1 t2 )(t̃1 t̃2 ) and ρ(′t ′ t ′ )(t̃ ′ t̃ ′ ) in equations (80) and (81),
                                     1 2        1 2
respectively. Note that ϱT1 T2 T̄2 T̄1 and ϱT̄′ 2 T̄1 T1 T2 are Grassmann-                                                       where U and V † are unitary matrices such that
even since MXT1 T2 X̄ ′ T̄2′ T̄1′ is Grassmann-even. As we have
demonstrated in section 2.5, the SVD of these Grassmann
                                                                                                                                                                even k : odd ]
                                                                                                                                                         [ k : (even)
                                                                                                                                             U=             U           0                             i : even ,            (95)
reduced density matrices results in the SVD of their coefficient
                                                                                                                                                                 0    U(odd)                           i : odd
tensors as in equations (82) and (83). The coefficient tensors
P and Q in equations (87) and (88) are then determined sim-
ilarly with the normal HOTRG. Since resulting P and Q are                                                                                                      even j : odd ]
                                                                                                                                                        [ j :†(even)
Grassmann-even, no extra sign factor appears in calculating                                                                                V† =           V             0                               l : even ,          (96)
the right-hand side of equation (86). The Grassmann tensor                                                                                                     0     V†(odd)                             l : odd
Tn ′ in equation (86)        [∏gives the] coarse-grained Grassmann
tensor network gTr n ′ ∈Λ ′ Tn ′ with the same ordering of                                                                       and σkl = σk δkl with the singular value σ k and Kronecker’s
Grassmann measures in equation (72), as is evident from                                                                          delta δ kl in the form of
equation (89). Therefore, we can iterate the coarse-graining
transformation as in figure 7 even for the Grassmann tensor                                                                                                    even l : odd ]
                                                                                                                                                        [ l : (even)
network, without any extra difficulty.                                                                                                        σ=           σ            0                            k : even .             (97)
    We have reviewed two types of TRG algorithms so far                                                                                                         0    σ (odd)                         k : odd
and both of them employ the local approximation based on
                                                                                                                                 Recalling equations (32) and (69), we can associate k and
the (HO)SVD to define the coarse-graining transformations.
                                                                                                                                 l with new auxiliary Grassmann variables via Ξk and Ξ̄l . In
By constructing a coarse-graining transformation that includes
                                                                                                                                 equations (95) and (97), ‘k : even (odd)’ should be understood
not only the local tensor(s) but also the effects of the surround-
                                                                                                                                 as Ξk is of the Grassmann-even (odd) parity. In the same way,
ing tensors (they are usually referred to as the environment),
                                                                                                                                 ‘l : even (odd)’ in equations (96) and (97) should be understood
we can construct a more accurate transformation. This kind of
                                                                                                                                 as Ξ̄l is Grassmann-even (odd). Although we have two kinds
improvement is called the second RG [29, 32, 144]. One of
                                                                                                                                 of Grassmann variables Ξk and Ξ̄l , Kronecker’s delta δ kl in σ kl
the other ways to improve the TRG algorithms is to remove
                                                                                                                                 enforces them to be of the same Grassmann parity. Therefore,
the redundant loop structure in the tensor network [28]. This
                                                                                                                                 we can easily define the new Grassmann parity functions cor-
can be achieved by the TNR [20, 145]. We will see several
                                                                                                                                 responding to the new auxiliary Grassmann variables.
TNR-type algorithms in section 4.
                                                                                                                                      In the practical TRG computations, the Grassmann parity
    In addition to the TRG approach, many other algorithms
                                                                                                                                 functions play a crucial role in restoring the Grassmann cal-
perform approximate tensor contractions, and they are used
                                                                                                                                 culus on your computer. In the case of the Levin–Nave TRG
according to their purpose and cost performance. For example,
                                                                                                                                 explained in section 2.5, the data that should be kept on the
corner transfer matrix RG [26] and time-evolving block
                                                                                                                                 computer are the coefficient tensor Tn;xtx ′ t ′ and Grassmann
decimation [146] are widely used in the tensor network com-
                                                                                                                                 parity functions for each subscript. Here, we repeat the dis-
putations based on the Hamiltonian formalism. See [130] as a
                                                                                                                                 cussion from equation (73) to equation (76). The two kinds
recent review.
                                                                                                                                 of SVD for the Grassmann tensors in equations (73) and (74)
                                                                                                                                 correspond to
2.7. Practical remarks
                                                                                                                                                                                           D∑
                                                                                                                                                                                            LNTRG

Here, we see how to carry out the SVD of Grassmann                                                                                                                   Tn;xtx ′ t ′ ≃                 An;xta Bn;ax ′ t ′ ,    (98)
tensors in practical computations. We again consider the                                                                                                                                    a=1
Grassmann-even tensor OΦΨ in equation (67) as an example.                                                                                                                                  D∑
                                                                                                                                                                                            LNTRG
                                                                                                                                             fx ′ ( fx +ft )+fx ft
Since OΦΨ is Grassmann-even, the corresponding coeffi-                                                                                (−1)                           T   n;xtx ′ t ′   ≃            Cn;xt ′ a Dn;axt ′ ,    (99)
cient tensor Oij takes a non-zero value only when Φi Ψj is                                                                                                                                  a=1


                                                                                                                            16
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                                            Topical Review



where A, B, C and D are defined in the similar way to                                                                           In the case of the finite-temperature system, we need to impose
equations (63) and (64). Notice that the extra sign factor                                                                      the anti-periodic boundary condition for the temporal direc-
in equation (99) comes from the rearrangement of auxiliary                                                                      tion, which can be done by using the Grassmann parity func-
Grassmann variables. Using the block-diagonalized basis, one                                                                    tion f t via
can immediately find the corresponding Grassmann parity for                                                                                         [      ] ∑
                                                                                                                                                                       f f +f (N)
the new index a. Equation (75), or the diagram in figure 6, is                                                                                 gTr Tn(N) =        (−1) x t t Tn;xtxt .     (107)
now ready to be computed. We begin with contracting A and                                                                                                      x,t
D, that gives us an intermediate Grassmann tensor,
                        ˆ                                                                                                       Note that the extra sign factor (−1)ft is a result of replacing T t
                                                                                                                                                 (N)
        (AD)T1 TX̄T2 =       An;X1 T1 T Dn+1̂;X̄X̄1 T2 , (100)                                                                  with (−T) t in Tn;XTX̄T̄ .
                                                X̄1 ,X1

which is described by the corresponding coefficient tensors as                                                                  3. Examples of numerical calculations
                       ∑                        f    ( ft1 +ft )+fx1′ fx ′
    (AD)t1 tx ′ t2 =             (−1) x1                                       An;x1 t1 t δx1 x1′ Dn+1̂;x ′ x ′ t2 .            3.1. Wilson–Majorana fermions
                                                                                                                 1
                       x1 ,x1′
                                                                                                               (101)            Since the two-dimensional classical Ising model is one of the
                                                                                                                                simplest models for critical phenomena with spontaneous Z2
The contraction between C and B is                                                                                              symmetry breaking, it has been used as a benchmark to valid-
                       ˆ                                                                                                        ate various TRG methods. Here, we begin with considering a
      (CB)T̄1 XT̄T̄2 =   Cn+2̂;X2 T̄1 X Bn+1̂+2̂;T̄X̄2 T̄2 ,                                                   (102)            fermionic model that is fully equivalent to the Ising model to
                                         X̄2 ,X2                                                                                see the validity of the Grassmann TRG.
                                                                                                                                   A simple example of classical action quadratic in
which corresponds with                                                                                                          Grassmann variables is given by the action of Wilson–
                                                (           )
                     ∑                   fx 2       ft ′ +fx +fx ′ ft ′
                                                                                                                                Majorana fermions:
(CB)t ′ xt ′ t ′ =             (−1)                  1                 2   Cn+2̂;x2 t ′ x δx2 x2′ Bn+1̂+2̂;t ′ x ′ t ′ .
          1    2                                                                       1                          2 2                                                                 
                     x2 ,x2′
                                                                                                                                       1∑                   ∑               ∑
                                                                                                                                                             2               2
                                                                                                                                                                          1
                                                                                                               (103)              S=         η̄ (n) mη +       γµ ∂µS −        ∂µ ∂µ∗  η (n)
                                                                                                                                       2 n                                2
                                                                                                                                                            µ=1             µ=1
Finally, the contraction between (AD) and (CB) gives a new                                                                                                                                
                                                                                                                                         1∑                    ∑               ∑
                                                                                                                                                                 2              2
Grassmann tensor in equation (75);                                                                                                                                           1
                                                                                                                                       +          χ̄ (n) mχ +      γµ ∂µS −        ∂µ ∂µ∗  χ (n)
                    ˆ  ˆ                                                                                                                 2 n                                 2
                                                                                                                                                               µ=1             µ=1
     Tn ′ ;XTX̄T̄ =         (AD)T1 TX̄T2 (CB)T̄1 XT̄T̄2 , (104)                                                                                       (                                   )
                                                                                                                                         ∑                             1          1
                                                                                                                                              η̄ (n) γ1 ∂1S − γ2 ∂2S − ∂1 ∂1∗ + ∂2 ∂2∗ χ (n)
                               T̄1 ,T1      T̄2 ,T2
                                                                                                                                       +
                                                                                                                                           n
                                                                                                                                                                       2          2
which is restored by
                                                                                                                                                                                             (108)
                                                                                                    (               )
                               fx ( ft +fx ′ )
                                                    ∑∑                         ft1 ( ft +fx ′ )+ft ′ ft ′ +fx +ft ′
    Tn ′ ;xtx ′ t ′ = (−1)                                             (−1)                       2     1
                                                                                                                                with two-component Majorana spinors η ≡ (η1 , η2 )T and χ ≡
                                                     t1 ,t1′ t2 ,t2′
                                                                                                                                (χ1 , χ2 )T . The forward, the backward, and the symmetric dif-
                   × (AD)t1 tx ′ t2 δt1 t1′ δt2 t2′ (CB)t ′ xt ′ t ′ ,                                         (105)            ference operators are defined by ∂, ∂ ∗ , and ∂ S = (∂ + ∂ ∗ )/2,
                                                                               1      2
                                                                                                                                respectively.
where the sign factor inside the summations originates from                                                                        An equivalence between the lattice action for Wilson–
reordering auxiliary Grassmann variables in (AD) and (CB).                                                                      Majorana fermions and the Ising model has been shown
The sign factor outside the summations comes from reordering                                                                    for two-dimensional honeycomb and square lattices for
auxiliary Grassmann variables after the summations6 .                                                                           any choice of the periodic/anti-periodic boundary condi-
    The Grassmann tensor trace is also described by the coef-                                                                   tions [147]. The masses of the fermions are functions of the
                                                      (N )
ficient tensor and parity functions. Suppose Tn;XTX̄T̄ be a                                                                     ‘reverse temperature’ κ:
Grassmann tensor obtained by N times of Levin–Nave TRG
or HOTRG transformation. The path integral defined on a lat-                                                                             2 (√           )             2 (√           )
                                                     (N)
                                                                                                                                 mη =         2 − 1 − κ , mχ = −            2 + 1 + κ . (109)
tice with 2N sites is approximately given by gTr[Tn ] and one                                                                            κ                            κ
finds                                                                                                                                                                     √
                     [     ] ∑                                                                                                  The critical point of the system κc = 2 − 1 is associated
                                       f f (N)                                                                                  with vanishing mη , and this value is related to the critical
                gTr Tn(N) =       (−1) x t Tn;xtxt .       (106)
                                                              x,t                                                               point of the two-dimensional Ising model on a square lattice
                                                                                                                                βc = tanh−1 κc . Even in such a model where there are several
                                                                                                                                species of fermions, one can follow the prescription given in
                                                                                                                                the previous section to derive a tensor network representation.
6One can see that the algorithm should result in the Z2 -symmetric Levin–                                                          One can see the equivalence to the Ising model from the
Nave TRG when we set all the Grassmann parity functions to zero.                                                                specific heat of this model in figure 8. Although we omit to

                                                                                                                           17
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                         Topical Review




Figure 8. Specific heats of the Majorana–Wilson fermion system on
several sizes of lattice. The specific heat is defined as the second
derivative of free energy that is calculated by the Grassmann TRG
with DLNTRG = 64.
                                                                            Figure 9. Adapted from [117]. Relative errors of the free energy.
                                                                            Reprinted (figure) with permission from [117], Copyright (2014) by
show a detailed finite-size scaling analysis, one can see the               the American Physical Society.
logarithmic
     √        divergence of the specific heat at the critical point
κc = 2 − 1.
   In the latter section, we will show the ‘renormalization
flow’ of this model given by the plain and the improved coarse-
graining methods.

3.2. The Schwinger model

The Schwinger model, which is two-dimensional QED, is the
most famous toy model for four-dimensional QCD. This is
because the Schwinger model has similar properties to QCD,
such as fermion confinement and chiral symmetry break-
ing. A θ term can also be introduced in the massive model,
which undergoes a phase transition at θ = π. The non-trivial
effect of θ term in QCD is related to the strong CP prob-
lem. In addition, the massless model is exactly solvable and
the mass perturbation is a useful technique to investigate the
model with sufficiently small finite mass. In the HEP com-                  Figure 10. Adapted from [117]. Lee–Yang zeros in the complex
munity, Shimizu and Kuramashi did a series of pioneering                    κ-plane. Reprinted (figure) with permission from [117], Copyright
works for the two-dimensional Schwinger model with Wilson                   (2014) by the American Physical Society.
fermions [117–119].
    In [117], where the Grassmann TRG was firstly applied to                increasing bond dimension (see figure 9). They defined the rel-
a lattice gauge theory, the Schwinger model                                 ative error by using the result obtained by a large bond dimen-
                                                                            sion DLNTRG = 128, while they fixed their bond dimension for
        1 ∑∑
                  2
 S=−              ψ̄ (n) {(1 − γµ ) Uµ (n) ψ (n + µ̂)                       the analyses to DLNTRG = 96.
        2 n                                                                    In the paper, finite-size scaling analyses is taken place for
              µ=1
                                          }                                 peak heights of the chiral susceptibility and the partition func-
      + (1 + γµ ) U†µ (n − µ̂) ψ (n − µ̂)                                   tion zeros in the complex κ-plane. Using the Lee–Yang zeros,
         1 ∑                      ∑       (           (     )
                                                                            they did some reliable fittings to find out the critical exponents
      +        ψ̄ (n) ψ (n) − β      cos A1 (n) + A2 n + 1̂
        2κ n                                                                (see figure 10). Note that such an investigation in the complex
                                   n
           (      )           )                                             parameter plane fully utilizes the absence of the sign problem
      −A1 n + 2̂ − A2 (n) ,                             (110)               for the TRG approach. Both in the strong coupling limit and
                                                                            for finite couplings, the critical exponents are shown to be the
where A is the phase of U(1) link variable, was analyzed for                same as those of the two-dimensional Ising model. Moreover,
some choices of the reverse coupling constant β: the strong                 the obtained critical mass (hopping parameter) is consistent
coupling limit β = 0 and finite couplings β = 5, 10.                        with a previous work with the eight vertex model [148].
   They confirmed their formulation and code by checking                       After this work, an analysis of the same model with the
convergences of relative errors of free energy along with                   presence of θ term that impressed the efficiency of the TRG

                                                                       18
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                  Topical Review



approach on the community follows [118]. In that work,                    4.2. Removal of CDL from network
they studied the critical behavior on θ = π line, and in con-
                                                                          The TNR approaches mentioned at the beginning of this
clusion, on the line θ = π, no phase transitions occur for
                                                                          section attempt to remove the CDLs from the network. Some
κ > κc , a second order phase transition that belongs to the
                                                                          of the authors of this review article has checked the efficiency
two-dimensional Ising universality class occurs at κ = κc ≃
                                                                          of the loop-TNR and the gilt-TNR for the Majorana–Wilson
0.2415, and first-order phase transitions occur for κ < κc . This
                                                                          fermion system seen in section 3.1 and for the two-flavor
is exactly the expected result.
                                                                          Gross–Neveu model [120]. While the loop-TNR comes first
    The Schwinger model with the topological term, where the
                                                                          in the chronological order, we show here how the gilt-TNR
staggered discretization for the fermions is applied, was stud-
                                                                          removes the CDL from a network for graphical simplicity.
ied by a different group [149]. In their paper, a special prop-
                                                                              In the gilt-TNR, recursive optimization steps are taken
erty of the model was mentioned: a tensor network represent-
                                                                          place to remove a CDL on a plaquette. To illustrate this, one
ation can be constructed without having Grassmann variables
                                                                          can consider an SVD of the plaquette (see figure 13), which
on the network. Indeed this property was found in the context
                                                                          can be seen as an environment of a link. An important point
of Monte Carlo simulation in [150, 151].
                                                                          here is that, by considering the link to be open7 , the internal
                                                                          (i.e. the CDL) loop is seen to be not enclosed in the plaquette
4. Improved TRG methods for fermions                                      and connected to and only to the open link. With this logic,
                                                                          the CDL is captured by the unitary matrix (U in the figure)
The Grassmann TRG methods introduced in sections 2.5                      through the SVD. After the decomposition, the unitary matrix
and 2.6 show remarkable performances as seen for some spe-                U is replaced by another one according to the ‘environment
cific models in the latter sections. However, the key com-                spectrum’ S so that the CDL loop will be truncated down.
ponent of the RG methods was the (HO)SVD that gives the                       While the gilt-TNR is based on local replacements of tensor
best approximation for local tensors [152] rather than the                legs, the loop-TNR consists of an entanglement filtering gauge
whole network. As we mentioned in the introduction, there                 transformation of tensors and an optimization step that minim-
are possible issues to describe the RG flows near critical                izes a bit global cost function compared to the usual TRG. One
points [27, 28, 99]. This is an important motivation to seek              can see how is the CDL eliminated by the gauge transforma-
refined RG techniques. For tensor networks, the quality of the            tion in the original paper [100], and see [28], where the notion
RG methods is often related to the ‘entanglement’ of the sys-             of ‘entanglement filtering’ firstly appeared.
tem. In [20], an improved renormalization technique where                     In [120], the performance of the Grassmann version of
(dis)entanglers are introduced to the network and are variation-          loop-TNR and gilt-TNR was inspected by computing the free
ally tuned so that the short-range correlation is removed at              energy and the determination of Fisher’s zeros. In the follow-
each RG step was introduced. Also, some new and less com-                 ing, we review the RG flow of the singular values that are
putationally demanding approaches like the loop-TNR [100]                 believed to store the information of the system.
and the gilt-TNR [101] were developed. Such approaches                        The RG flow obtained by the plain Grassmann TRG and the
are called TNR. Basically, the TNR approaches require some                loop-TNR are shown in figures 14 and 158 . The vertical and
optimization steps on the network, but one can easily adapt               the horizontal axes represent the normalized singular value and
them to the Grassmann tensor networks since the treatment                 how many iterations were taken place before that, respectively.
of the Grassmann variables is exact and factored out from the             Under these iterative coarse-graining algorithms, the space-
bosonic part of the network. Note that, later in this section, we         time volume of the system grows rapidly along with the num-
also show the Grassmann version of the bond-weighted TRG                  ber of iterations; schematically it grows twice at each coarse-
(BTRG), which represents a different notion than the TNRs.                graining step. In figure 14 one cannot clearly distinguish the
                                                                          three panels where the reverse temperature is set to 0.9999κc ,
                                                                          κc , and 1.0001κc . This is on account of a contamination by
4.1. CDL structure on tensor network                                      short-range information that is difficult to properly remove by
                                                                          the normal Grassmann TRG algorithm. On the other hand, the
The difficulty in approaching critical points with the TRG has            loop-TNR shows distinguishable fixed point structures at off-
been illustrated by the CDL picture [28]. As a typical example,           critical points such as 0.9999κc , 1.0001κc . Also, at the critical-
one can easily show that a toy tensor network that consists of            ity κc , a scale-invariant structure is observed that clearly shows
the CDL tensor                                                            the superiority of the improved renormalization algorithm.
                                                                          Note that these behaviors of the singular values are qualit-
                      CDL                                                 atively the same as those observed for the two-dimensional
                     Tijkl = Λi1 l2 Λj1 i2 Λk1 j2 Λl1 k2 ,   (111)
                                                                          Ising model; the equivalence between the Ising model and the
                                                                          Majorana–Wilson fermion system can be seen in a sense like
where Λ can be assumed to be a diagonal matrix for simpli-
city, is a fixed point of the TRG [28] (see figures 11 and 12).           7 Of course, the open link is to be closed after all.
This means the usual TRG leaves short-range correlations in               8 In the reference, the result of the gilt-TNR is also shown; however, we omit
the network under each blocking step and this causes a loss of            to show it since the characteristic behavior is quite similar to that of the loop-
accuracy.                                                                 TNR.

                                                                     19
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                         Topical Review




                                             Figure 11. CDL tensor and decompositions.


                                                                         the computational cost [153]. Their idea is to introduce some
                                                                         weight on the edge of the tensor network and construct the
                                                                         coarse-graining transformation including these weights. The
                                                                         algorithm is referred to as the BTRG, which can be regarded
                                                                         as a generalization of the Levin–Nave TRG. The advantage
                                                                         of the BTRG over the Levin–Nave TRG is demonstrated by
                                                                         benchmarking with the two-dimensional Ising model in [153].
                                                                         The authors show that the BTRG outperforms the Levin–Nave
                                                                         TRG and the HOTRG at the same bond dimension. Both
Figure 12. Coarse-graining of CDL tensor network leads to a CDL          for the BTRG and the Levin–Nave TRG, their computational
tensor network.                                                          times are proportional to O(D5 ) and their memory footprints
                                                                         are O(D3 ). In the two-dimensional HOTRG, the computa-
                                                                         tional time scales with O(D7 ), and the scaling of the memory
                                                                         cost is O(D4 ). Therefore, the BTRG has shown the best per-
                                                                         formance among these three algorithms. Moreover, the authors
                                                                         numerically demonstrate that non-trivial fixed point tensors
                                                                         can be constructed in the thermodynamic limit by the BTRG.
                                                                         Recently, the BTRG has been applied to investigate the phase
                                                                         structure of the CP(1) model with a topological θ term [154].
                                                                            We review the BTRG algorithm for the two-dimensional
                                                                         square tensor network which is generated by a four-leg fun-
                                                                         damental tensor. The algorithm is schematically explained in
                                                                         figure 16. The key step in the BTRG is the low-rank approx-
                                                                         imation of the four-leg tensor based on the SVD with a hyper-
                                                                         parameter such as

Figure 13. SVD of plaquette. After this decomposition CDL loop is                             ∑
                                                                                              D
                                                                                                             (1−k)/2 k (1−k)/2 †
inside U.                                                                           Tabcd ≃            Uabi σi      σi σi     Vicd ,      (112)
                                                                                              i =1


this. For the growth of the entanglement entropy and relation-           where k ∈ R denotes the hyperparameter. If we set k = 0,
ship to the Calabrese–Cardy formula, we refer the reader to              equation (112) exactly corresponds to the tensor decom-
the discussion in [120].                                                 position employed in the Levin–Nave TRG as shown in
                                                                         equations (63) and (64). With k ̸= 0, we obtain an extra factor
                                                                         σik , which is regarded as a weight on the bond i. The authors
4.3. Bond-weighting technique                                            in [153] have given a stationary condition equation,
TRG algorithm can be improved by removing the short-range                                     [              ]4 [ ]4
                                                                                                   (1−k)/2
correlations represented by the loop entanglement. Several                                        σi             σik = σi ,               (113)
algorithms are proposed to remove these short-range cor-
relations and achieve much higher accuracy than the ori-                 that determines the optimal choice of the hyperparameter k.
ginal TRG algorithm even at the criticality as demonstrated              Equation (113) enforces the singular value spectrum to be
above. At the same time, however, these algorithms usu-                  invariant under the sequential coarse-graining procedure in
ally require more computational cost than the original TRG.              the BTRG, assuming that the local four-leg tensors and the
Recently, Adachi, Okubo, and Todo have proposed a new idea               bond weights converge after sufficiently many times coarse-
to improve the accuracy of TRG algorithms without increasing             graining transformations and the matrices U and V do not

                                                                    20
J. Phys.: Condens. Matter 36 (2024) 343002                                                                               Topical Review




Figure 14. Normalized singular values produced by Grassmann TRG at κ = 0.9999κc (left), κ = κc (middle), and κ = 1.0001κc (right).
The bond dimension is set to 64.




Figure 15. Normalized singular values produced by Grassmann loop-TNR at κ = 0.9999κc (left), κ = κc (middle), and κ = 1.0001κc
(right). The bond dimension is set to 16.


affect the spectrum. According to equation (113), the optimal            shows that the bond-weighting method improves the accuracy
value of the hyperparameter for the two-dimensional square               of the Grassmann TRG. Figure 17 shows the relative error of
tensor network should be k = −1/2, which has been numer-                 the free energy for the two-dimensional free massless Wilson
ically confirmed in [153]. We note that the application of the           fermion. With the fixed bond dimension, the bond-weighting
bond-weighing technique to other TRG algorithms has been                 method always achieves higher accuracy compared with the
discussed in [155].                                                      normal TRG. Notice that both computations require the same
   It should be emphasized that this derivation is completely            computational complexity at the same bond dimension. A
independent of the details of the lattice theory. Actually, [126]        sample implementation of the BTRG for the Gross–Neveu


                                                                    21
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                       Topical Review




Figure 16. Schematic illustration of the BTRG algorithm. Background dotted lines denote a real-space square lattice. (A) Initial tensor
network with bond weights on the lattice. (B) SVD introduces three-leg tensors and new bond weights. (C) New tensor network with bond
weights by contracting four types of three-leg tensors and four bond weights.


model with Wilson fermions at finite density is shown in [128],
whose web documentation is also provided9 . Using the code,
one can reproduce the result shown in figure 17.

4.4. Multilayered tensor network formulations for Nf -flavor
fermions

There is no difficulty in expressing the path integral of the lat-
tice fermion system as the Grassmann tensor network with
N f flavors. In practice, however, the size of the resulting
Grassmann tensor scales exponentially for N f and a O(eNf )
computational memory is required in the numerical compu-
tations. This issue has been started to be addressed recently by
Akiyama [103] and also by Yosprakob et al [127].
    Akiyama [103] has employed the matrix product decom-
position (MPD) to introduce a virtual direction so that each fla-
                                                                          Figure 17. Relative error of the free energy for the two-dimensional
vor degree of freedom is assigned to the different layers ortho-          free massless Wilson fermion. The hyperparameter is set at
gonal to the virtual direction. MPD is a common idea in the               k = −1/2.
tensor network methods such as the MPS and matrix product
operator [156, 157]. Akiyama [103] has particularly utilized
a canonical form of the MPD proposed in [158]. Thanks to
the MPD, the memory cost for each local Grassmann tensor                  is reduced from O(eNf ) to O(Nf ), and the technique has been
                                                                          benchmarked with the two-dimensional Gross–Neveu model
                                                                          at finite density with the Nf = 2, 3 Wilson fermions. Although
9   https://github.com/akiyama-es/Grassmann-BTRG.                         a naive formulation provides the two-dimensional Grassmann


                                                                     22
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                    Topical Review




Figure 18. Adapted from [103]. Multilayered Grassmann tensor network representation for the path integral of the two-dimensional Nf = 3
Gross–Neveu model with Wilson fermions. (A) Original Grassmann tensor network representation defined on the square lattice. (B)
Three-layered Grassmann tensor network. (C) MPD rewrites the original fundamental tensor (yellow) with three fundamental tensors
(green) and two kinds of singular values (purple). Reproduced from [103]. CC BY 4.0.


tensor network whose bond dimension equals 4Nf , the MPD                computations in [129], whose web documentation is also
alternatively gives the (2 + 1)-dimensional network, N f sites          provided10 . The Schwinger model is implemented as an
along the virtual direction, with the bond dimension four               example.
without any approximation. The schematic picture is shown in
figure 18. Based on the latter representation, the Silver Blaze         5. Relativistic models with fermion interactions
phenomenon in the pressure and number density is reproduced
with a relatively small bond dimension.                                 5.1. Gross–Neveu model
   Yosprakob et al [127] has proposed a compression scheme
for the initial tensor representation of lattice gauge theories         The Gross–Neveu model is a well-known toy model for QCD.
with the N f -flavor fermions. They have introduced replicas of         Since it shares several important features of the QCD such as
the original gauge field, which allows one to separate the local        asymptotic freedom and a dynamical mass generation mech-
Grassmann tensor into multiple layers associated with the fer-          anism via symmetry breaking, the model is a good test bed
mion flavor. Based on this description, each layer is individu-         for new computational methods. Here, we consider the model
ally compressed by the isometry insertion; In the case of the           defined with Wilson fermions. The lattice action reads
ZK gauge theory with the N f -flavor Wilson fermions, the ori-
                                                                                   1 ∑ ∑ { µδν,2 ( f)
                                                                                                Nf
ginal tensor whose size is 164 K10 is converted into the com-                  S=−              e    ψ̄ (n) (r1 − γν ) ψ ( f) (n + ν̂)
pressed one with D4 K2 elements, where D is the bond dimen-                        2 n,ν
                                                                                         f =1
sion introduced by the isometries. Even with D = 8, the differ-                                                                 }
                                                                                      −µδν,2 ( f)
ence between the resulting ln Z at finite gauge coupling and                       +e         ψ̄ (n + ν̂) (r1 + γν ) ψ ( f) (n)
chemical potential is suppressed less than O(10−15 ). Since                          ∑∑
gauge fields are replicated by the Kronecker deltas, they are                      +           (m + 2r) ψ̄ ( f) (n) ψ ( f) (n)
                                                                                            n   f =1
maximally entangled along the flavor direction. Therefore, the                                                   2
tensor contractions along the flavor direction are carried out
                                                                                         g2σ ∑ ∑ ( f)
before the two-dimensional spacetime coarse-graining is per-                           −         ψ̄ (n) ψ ( f) (n)
formed. Using this technique, the chiral susceptibility of the                           2Nf n
                                                                                                         f
two-dimensional infinite-coupling Z2 , Z4 , and U(1) gauge the-                                                          2
ories with the Nf = 1, 2 Wilson has been computed. The crit-                              2 ∑ ∑
                                                                                        g     
ical hopping parameters have been determined for each case.                            − π      ψ̄ ( f) (n) iγ5 ψ ( f) (n) ,        (114)
                                                                                        2Nf n
                                                                                                         f
The pressure and number density as functions of the chemical
potential have also been provided in the case of Z2 gauge the-
                                                                        where ψ ( f ) (n) and ψ̄ ( f ) (n) are the Wilson fermions with the
ory with Nf = 1, 2, 4 at finite gauge coupling, where the Silver
                                                                        flavor index f. g2σ and g2π are the four-fermi coupling constants
Blaze phenomenon has been successfully captured.
   In addition, Yosprakob has been recently develop-
ing a Python package for the Grassmann tensor network                   10   https://ayosprakob.github.io/grassmanntn/.


                                                                   23
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                     Topical Review




Figure 19. Adapted from [103]. Pressure (left) and number density (right) of the Nf = 2 Gross–Neveu model with Wilson fermions. D
denotes the bond dimension for spacetime indices and χ does the bond dimension along the virtual direction in the multilayered network.
m = 1 and g2 = 10. Reproduced from [103]. CC BY 4.0.

                                                                                  ∑{               [           (       )
and m, µ, and r represent mass, chemical potential, and the                  S=     η1 (n) γ ψ̄ (n) eµ U1 (n) ψ n + 1̂
Wilson parameter, respectively.                                                    n
                                                                                           (      ) (       )]
    Takeda and Yoshimura studied the model with Nf = 1 at                         − eµ U†1 n − 1̂ ψ n − 1̂
finite density in [102], which is the first application of the                                   [        (      )     (      ) (       )]
Grassmann TRG to the finite-density system. They employed                         + η2 (n) ψ̄ (n) U2 (n) ψ n + 2̂ − U†2 n − 2̂ ψ n − 2̂
the Levin–Nave TRG algorithm with the bond dimension up                                            }
                                                                                  +2mψ̄ (n) ψ (n) ,                                   (115)
to DLNTRG = 64 to compute the fermion number density and
its susceptibility as functions of µ. They successfully observed
                                                                           where m and µ denote the mass and chemical potential.
that the number density saturated to one with sufficiently large
                                                                           The lattice site is labelled by n = (n1 , n2 ) with the tem-
chemical potential. In addition, they considered the model on
                                                                           poral coordinate n1 and spatial coordinate n2 . The action has
an anisotropic finite lattice with (N1 , N2 ) = (64, 32), (96, 32),
                                                                           an anisotropy factor γ in the temporal hopping terms. The
where N 1 and N 2 denote the spatial and temporal lattice sizes
                                                                           link variables Uν (n) take their value on SU(3) and ψ(n)
respectively, and found that there were two peaks in the sus-
                                                                           and ψ̄(n) are the staggered quark fields described by the
ceptibility under the presence of the finite four-fermi coupling.
                                                                           three-component Grassmann variables. The staggered sign
Throughout their analysis, they pointed out that the finite bond
                                                                           function is defined via η1 (n) = 1 and η2 (n) = (−1)n1 . The
dimension effect could be enhanced not only at critical points
                                                                           (anti-)periodic boundary condition is assumed in the spatial
but also in crossover regions.
                                                                           (temporal) direction. The authors integrate out the SU(3) link
    Recently, Akiyama has investigated the model with Nf =
                                                                           variables exactly as shown in [159, 160] and the path integ-
2, 3 at finite density [103]. The pressure and number density
                                                                           ral is rewritten by the mesonic and baryonic contributions,
on a square lattice were computed as functions of chemical
                                                                           following [161], characterized by their occupation numbers.
potential using two methods, the bond-weighting method and
                                                                           Note that this kind of dual formulation has also been applied
multilayered formulation as shown in figures 19 and 20. In the
                                                                           in [162] to investigate the infinite-coupling U(N) gauge the-
zero-temperature limit, we observe the Silver Blaze phenom-
                                                                           ories with staggered fermions in three and four dimensions
ena, where the thermodynamic quantities show no depend-
                                                                           with the variant of HOTRG. Effectively reducing the number
ence on µ as long as µ is smaller than the mass of the lightest
                                                                           of configurations with vanishing contribution in the path integ-
excitation. The results obtained by the two methods are con-
                                                                           ral, they obtain the Grassmann tensor network representation
sistent and the number density saturates to two and three for
                                                                           where the bosonic tensor is defined by a four-leg tensor and
Nf = 2, 3, respectively. They noted that the finite bond dimen-
                                                                           each subscript is of dimension six. The resulting Grassmann
sion effect could be enhanced in the Silver-Blaze regime in the
                                                                           tensor network is approximately contracted by the HOTRG,
multilayered formulation.
                                                                           where the authors have applied an improved method to derive
                                                                           the projectors in the algorithm [163].
                                                                               The authors calculated the thermodynamic potential ln Z/V
5.2. QCD in the infinite-coupling limit
                                                                           as a function of µ and m on 22 and 42 lattices and confirmed
Bloch and Lohmayer made the first study of QCD in                          that the exact results were reproduced. Although the conver-
the infinite-coupling limit with the Grassmann TRG                         gence for the bond dimension DHOTRG becomes slower with
approach [104]. The lattice action with the staggered quarks               the smaller m, it has been shown that the accuracy of the
is given by                                                                calculations is sufficient to take the chiral limit: a quadratic




                                                                      24
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                   Topical Review




Figure 20. Adapted from [103]. Pressure (left) and number density (right) of the Nf = 3 Gross–Neveu model with Wilson fermions. D
denotes the bond dimension for spacetime indices and χ does the bond dimension along the virtual direction in the multilayered network.
m = 1 and g2 = 10. Reproduced from [103]. CC BY 4.0.


fit in 1/DHOTRG successfully extrapolates ln Z/V at m = 0 on                          with η1 (n) = 1. m, g0 , and µ represent the mass, four-
a 10242 lattice to DHOTRG → ∞ using the results on 40 ⩽                               fermi coupling constant, and chemical potential. When
DHOTRG ⩽ 128.                                                                         m = 0, equation (116) is invariant under the continuous
    The authors then computed the chiral condensate ⟨ψ̄ψ⟩                             transformation,
at vanishing chemical potential. They extrapolated ⟨ψ̄ψ⟩
at finite mass to DHOTRG → ∞ before taking the infinite-                                 χ (n) → eiαη5 (n) χ (n) ,    χ̄ (n) → χ̄ (n) eiαη5 (n) .   (117)
volume limit. For m ⩽ 0.005, the condensate was nicely fit-                           This global symmetry is regarded as the chiral symmetry for
ted by limV→∞ ⟨ψ̄ψ⟩ = amb , with a = 2.77 and b = 0.0414.                             the staggered NJL model.
Therefore, it is confirmed that the chiral symmetry is not                                Decomposing the hopping terms and four-fermi interac-
dynamically broken in infinite-coupling QCD with staggered                            tion term, the Grassmann tensor network representation for
quarks in two dimensions. They also studied the number dens-                          the path integral is available. Akiyama [103] introduces the
ity and chiral condensate at finite chemical potential with                           eight kinds of local tensors to describe the Grassmann tensor
m = 0.1 and DHOTRG = 64. In the infinite-volume limit, they                           network. This is because the staggered theory defined in
have found that a first-order phase transition takes place at                         equation (116) a little bit breaks the translational invariance
µc ≃ 0.3508. With µ > µc , they have observed that the num-                           on the lattice. Note that the site dependence appearing in the
ber density saturates to three and the chiral symmetry is                             local tensor is characterized just by ην (n) so that the resulting
restored. Throughout the study, they employed the stabil-                             Grassmann tensor network has a periodic structure. See refer-
ized finite-difference method developed in [164] to calculate                         ences [103, 165] for the detailed derivation of the Grassmann
the chiral condensate and number density by the numerical                             tensor network representation.
differentiation.                                                                          The authors in [103] developed the anisotropic TRG
                                                                                      (ATRG) algorithm [166] for fermions. The ATRG allows
5.3. NJL model                                                                                    ( 2d+1 ) contract d-dimensional
                                                                                      us to approximately                        ( +1 tensor
                                                                                                                                         )     networks
                                                                                      with the O DATRG     complexity and O DdATRG         memory cost.
Akiyama, Kuramashi, Yamashita, and Yoshimura made the                                                                                   ( 4d−1 )
                                                                                                      ( 2d be compared
                                                                                      These costs should        )           with the O DHOTRG      com-
first application of the TRG approach to the four-dimensional                         plexity and O DHOTRG memory cost in the HOTRG. This
lattice fermions [105]. They investigated the chiral phase                            drastic cost reduction is a result of an additional approximation
transition in the NJL model at finite density. Since this model                       for fundamental tensors. The further cost reduction technique
is an effective field theory of the QCD, the efficiency of the                        for the ATRG has been provided by Oba in [167]. In addition
TRG approach for the NJL model should be addressed from                               to the ATRG, several other algorithms have been proposed for
the viewpoint of the future application of the TRG method                             the higher-dimensional systems [168, 169]. As a validation of
toward the QCD at finite density. The model is defined with                           the ATRG for fermionic models, they first considered the NJL
the staggered fermions by the following action,                                       model in the heavy-dense limit, where m → ∞ and µ → ∞
                                                                                      keeping the ratio of eµ /m fixed. In this limit, the model can
      1 ∑∑            [                                                      ]
           4
 S=            ην (n) eµδν,4 χ̄ (n) χ (n + ν̂) − e−µδν,4 χ̄ (n + ν̂) χ (n)            be solved analytically [170]. The number density and fermion
      2 n                                                                             condensate as functions of µ in the thermodynamic and van-
          ν=1
         ∑                    ∑∑
      +m    χ̄ (n) χ (n) − g0          χ̄ (n) χ (n) χ̄ (n + ν̂) χ (n + ν̂) ,          ishing temperature limits were calculated by the ATRG with
            n                     n   ν                                               the bond dimension DATRG ⩽ 30 and the results showed a good
                                                                        (116)         agreement with the analytic ones.
                                                                                          The authors then studied the chiral phase transition in the
where χ(n) and χ̄(n) are the staggered fermions described                             cold and dense regime characterized by µ/T = O(103 ). Due
by the single-component Grassmann variables. ην (n)                                   to the sign problem, such a cold and dense regime is inaccess-
is the staggered sign function ην (n) = (−1)n1 +···+nν−1                              ible with the standard Monte Carlo simulation. They enlarged

                                                                                 25
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                           Topical Review



the bond dimension up to DATRG = 55 in the strong coupling
regime, where they found a clear signal of the first-order trans-
ition in the chiral condensate, and the chiral symmetry was
restored with µ > µc . This is exactly the predicted result by
the mean-field theory [171] and the functional RG [172]. The
authors also computed the pressure and number density which
are fundamental ingredients in the equation of state. All the
thermodynamic quantities obtained by the ATRG have shown
that the model undergoes the first-order phase transition in the
cold and dense region.

5.4. N = 1 Wess–Zumino model

The interacting two-dimensional N = 1 Wess–Zumino model,
                                                                                Figure 21. Adapted form [74]. The partition function of the N = 1
a simple supersymmetric model, shows a vanishing partition                      free Wess–Zumino model as a function of m on V = 2 × 2 lattice.
function (Witten index) [173]; i.e. This model suffers from a                   Reprinted (figure) with permission from [74], Copyright (2022) by
serious sign problem as in the case of other generic supersym-                  the American Physical Society.
metric models. (See [174] for a review.)
   The Euclidean continuum action of the model is defined by
     ˆ        {                                                    }               Even though they showed results only for the non-
             1          1          1 (                  )                       interacting case, their construction of the tensor network does
S=        2
         d x   (∂µ ϕ)2 + W ′ (ϕ)2 + ψ̄ γµ ∂µ + W ′ ′ (ϕ) ψ             ,
             2          2          2                                            not depend on the shape of the superpotential, so that further
                                                                 (118)          studies of the interacting Wess–Zumino model are awaited.

where ϕ and ψ are a one-component real scalar field and a two-
component Majorana spinor field, respectively. The superpo-                     5.5. Non-abelian lattice gauge theories coupled to fermions
tential W (ϕ) is a function of ϕ and is the source of the Yukawa-               Non-abelian gauge theories have huge internal degrees of free-
and ϕn -interactions.                                                           dom, and this fact prevents one from building a tensor net-
   The Majorana condition for ψ is given by                                     work representation in a non-expensive way. Indeed, tensor
                                                                                network studies for gauge theories are limited to abelian cases
                             ψ̄ = −ψ T C−1                       (119)          when considering coupling to fermions. Recently a subgroup
                                                                                of the authors reported a way to construct a tensor network
with the charge conjugation matrix C,                                           representation of non-abelian gauge theories with a reason-
                                                                                able numerical complexity [107]. In their work, parameter-
     CT = −C,          C† = C−1 ,            C−1 γµ C = −γµT .   (120)
                                                                                ized group elements are discretized via a Gaussian quadrat-
                                                                                ure method, and a further reduction of the degrees of freedom
The continuum action above can be shown to be invariant
                                                                                is achieved by applying the higher-order orthogonal iteration
under the supersymmetry transformation
                                                                                (HOOI) algorithm [175, 176] to plaquette tensors. Moreover,
                                                                                for the fermion part, they adopt the reduced staggered formu-
                     δϕ = ϵ̄ψ,                                   (121)
                                                                                lation [177] to completely eliminate the redundancies for the
                     δψ = (γµ ∂µ ϕ − W ′ (ϕ)) ϵ,                 (122)          fermion part.
                                                                                    They numerically showed the accuracy of their approx-
where ϵ is a two-component Grassmann number that satisfies                      imation by checking the partition function and the average
the Majorana condition (119).                                                   plaquette for the SU(2) case. Surprisingly, the convergence of
    Figure 21 shows the partition function of the free N =                      the HOOI algorithm is so fast, and the loss of accuracy by the
1 Wess–Zumino model, whose superpotential is defined by                         truncation is shown to be quite mild. In other words, the main
W(ϕ) = (1/2)mϕ2 with the mass parameter m, on V = 2 × 2                         source of the error is the quadrature method with a few num-
lattice. (Note that the sign problem occurs even in this free                   ber of Gaussian nodes; they adopted up to 5 for the quadrature
case.) The periodic boundary conditions are imposed in all dir-                 at each angle. It would be interesting to see how such a drastic
ections for both fermions and bosons. For this specific case, the               approximation is tolerable for more complicated cases.
analytical solution can be shown to be 1 for the m > 0 region
shown in the figure. The Grassmann TRG results show better
agreement with the exact solution at larger masses although                     6. The Hubbard model
smaller masses seem to be difficult. As for the difficulty in the
small mass region, the authors of [106] concluded that it is due                The Hubbard model is one of the most fundamental lattice fer-
to the lack of damping factor in the local Boltzmann weight; in                 mion models describing the itinerant spin-1/2 electrons via the
other words, the reason is that the Gauss–Hermite quadrature                    repulsive Coulomb interaction. The model is defined by the
applied to the scalar boson part does not converge.                             following Hamiltonian,

                                                                           26
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                                                    Topical Review


              ∑∑(                                   )    ∑
   H = −t                    c†i,s cj,s + c†j,s ci,s + U   ni,↑ ni,↓ .   (123)        the spin-1/2 electrons, ⟨N⟩ = 2|Λ| with the total lattice sites
              ⟨ij⟩ s=↑,↓                               i∈Λ                            |Λ|. When the filling is a single electron per site, that is ⟨N⟩ =
                                                                                      |Λ|, the situation is referred to as half-filling. At half-filling,
The first term, the tight-binding Hamiltonian, represents                             the system acquires the so-called particle-hole symmetry. To
the kinetic energy of electrons and the second one shows                              see this, let us consider the transformation such that
the repulsive Coulomb potential with U ⩾ 0. ⟨ij⟩ labels the
nearest-neighbor sites in the d-dimensional hypercubic lat-                                                           ci,s → ηi c†i,s ,                              (129)
tice Λ. The interesting aspect of the Hubbard model origin-
                                                                                      where ηi = ±1 on the even and odd lattice sites, respectively.
ates from the fact that each term is diagonalizable in a differ-
                                                                                      The Hamiltonian in equation (123) is then transformed as
ent space; the kinetic term is diagonalized in the momentum
space but the Coulomb potential term is diagonalized in the                                                   H → H + U (|Λ| − N) ,                                  (130)
real space. Exact solutions of the Hubbard model are available
only in the cases of d = 1 [178] and d = ∞ [179, 180], despite                        which means that equation (129) is the symmetry if the system
its simplicity. In the condensed matter community, the model                          is at half-filing ⟨N⟩ = |Λ|. From the viewpoint of the grand
in general dimensions and on various lattice geometries have                          canonical Hamiltonian,
been extensively studied with mean-field approaches, field-
                                                                                               H − µN → H + U (|Λ| − N) + µN − 2µ|Λ|,                                (131)
theoretical ways, and numerical methods. We recommend for
interested readers to see several recent reviews [108–114] and                        and equation (129) does describe the symmetry setting µ =
references therein.                                                                   U/2. For a detailed explanation of particle-hole symmetry, see
    Let us review the Hubbard model very briefly. The cre-                            [112], for instance.
ation and annihilation operators c†i,s and ci,s satisfy the anti-                        Let us now formulate the Hubbard model within the path-
commutation relation,                                                                 integral formalism. The path-integral representation of the
         {                   }                                                        grand partition function is
           ci,s , c†i ′ ,s ′ = δii ′ δss ′ , {ci,s , ci ′ ,s ′ } = 0,                      ˆ              [ ˆ                  {                    (          )
         {                   }                                                                 [      ]               β            ∑                    ∂
                                                                                                dψ̄ dψ exp −                                               − µ ψ ( n, τ )
           c†i,s , c†i ′ ,s ′ = 0.                                    (124)           Z=
                                                                                                                  0
                                                                                                                          dτ           ψ̄ (n, τ )
                                                                                                                                                        ∂τ
                                                                                                                                n∈Λ

The number operator is defined via                                                              ∑∑
                                                                                                 d
                                                                                                   (                                                                 )
                                                                                           −t            ψ̄ (n + σ̂, τ ) ψ (n, τ ) + ψ̄ (n, τ ) ψ (n + σ̂, τ )
                             ∑                                                                   n σ=1
                        N=      ni,s                                     (125)                                             }]
                                                                                             U(                     )2
                                              i,s                                          +    ψ̄ (n, τ ) ψ (n, τ )            ,                                    (132)
                                                                                             2
with ni,s = c†i,s ci,s . The spin operator at the site i is given by                  where the two-component Grassmann variables ψ and ψ̄ have
        (x)    (y )   (z)
Si = (Si , Si , Si ), where                                                           been defined via
                                                                                                 [           ]
                                                                                                                             [                         ]
                            (ν)       1 ∑ † (ν)                                                   ψ (n, τ )
                                                                                      ψ (n, τ ) = ↑            , ψ̄ (n, τ ) = ψ̄↑ (n, τ ) , ψ̄↓ (n, τ ) .
                        Si        =      c σ ′ ci,s ′ .                  (126)                    ψ↓ (n, τ )
                                      2 ′ i,s ss
                                        s,s                                                                                                       (133)
σ (ν) (ν = x, y, z) is the Pauli matrices. The total spin operator                    The imaginary time is parameterized by τ and β is the inverse
is                                                                                    temperature. The Grassmann fields obey the anti-periodic
                                 ∑                                                    boundary condition along the imaginary time direction. When
                             S=      Si .                   (127)                     the imaginary time direction is discretized by β = Nτ ϵ, we can
                                               i                                      identify
The Hamiltonian in equation (123) has the global U(2) =                                          ∑ [             ψ (n + τ̂ ) − ψ (n)
U(1) × SU(2) symmetry. With R ∈ U(2), equation (123) is                                       S=        ϵ ψ̄ (n)
                                                                                                      ′
                                                                                                                          ϵ
                                                                                                   n∈Λ
invariant under the transformation,
                             ∑                                                                          ∑
                                                                                                        d
                                                                                                          (                                        )
                      ci,s →    Rs,s ′ ci,s ′ .   (128)                                            −t         ψ̄ (n + σ̂) ψ (n) + ψ̄ (n) ψ (n + σ̂)
                                         s′                                                             σ=1
                                                                                                                                                         ]
                                                                                                       U(             )2
This symmetry results in the global charge conservation and                                        +      ψ̄ (n) ψ (n) − µψ̄ (n) ψ (n)                               (134)
                                                                                                       2
the spin isotropy. The global charge conservation originates
from the global U(1) symmetry and the total particle number                           as the action of the Hubbard model on the (d + 1)-dimensional
N is a good quantum number. The spin isotropy is a result of                          anisotropic lattice Λ ′ . The unit vector along the imagin-
the global SU(2) symmetry and S(z) and S2 are good quantum                            ary time direction is denoted as τ̂ in equation (134). The
numbers.                                                                              Grassmann fields ψ(n) and ψ̄(n) on Λ ′ are defined in the same
   The grand canonical Hamiltonian is given by H − µN,                                manner with equation (133). With µ = U/2, the system is at
where µ is the chemical potential. Since the model describes                          half-filling as before.

                                                                                 27
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                            Topical Review




Figure 22. Adapted from [115]. Electron number density ⟨n⟩ at
(U, t) = (4, 1) with DHOTRG = 80 and ϵ = 10−4 . N σ and N τ denote
the spatial and temporal lattice sizes, respectively. Reproduced from
[115]. CC BY 4.0.


6.1. (1+1)-dimensional model

Akiyama and Kuramashi provided a benchmark study of the
TRG method for the (1 + 1)-dimensional Hubbard model in
[115]. They also derived the Grassmann tensor network rep-
resentation of the path integral defined by equation (134)
in general spatial dimensions. In equation (134), the time
slice ϵ should be ϵ ≪ 1 to reproduce the grand partition
function in equation (132). In other words, the contribution
from the temporal hopping terms (O(1)) becomes much larger
                                                                             Figure 23. Adapted from [115]. (top) Electron number density ⟨n⟩
than that from the spatial ones (O(ϵ)) in the original action.
                                                                             at β = Nτ ϵ = 1677.7216 with ϵ = 10−4 and DHOTRG = 80. The fit
Consequently, the initial Grassmann tensor network becomes                   ansatz in equation (135) results in µc (DHOTRG ) = 2.698(1) and
extremely anisotropic in the temporal direction. The authors in              ν = 0.51(2). (bottom) µc (DHOTRG ) as a function of 1/DHOTRG . The
[115] applied the HOTRG algorithm along the temporal dir-                    solid line shows the fit result of equation (136) and the dotted curve
ection in advance before the spacetime coarse-graining took                  gives the fit result of equation (137). Reproduced from [115].
place to investigate the ground state.                                       CC BY 4.0.
    Throughout their study, the number density ⟨n⟩ is computed
as a function of the chemical potential µ on the three points
                                                                             6.2. (2+1)-dimensional model
(U, t) = (4, 0), (0, 1), (4, 1). At finite U, the Mott plateau ⟨n⟩ =
1 is reproduced as shown in figure 22. At (U, t) = (4, 1), they              Akiyama, Kuramashi, and Yamashita investigated the metal–
provided the numerical fit of the number density via                         insulator transition in the (2+1)-dimensional Hubbard
                                                                             model [116]. They used the ATRG algorithm to investigate
                 ⟨n⟩ = A + B|µ − µc (DHOTRG ) |ν .             (135)         the ground state and applied the same strategy with [115]; the
                                                                             spacetime coarse-graining after the imaginary-time evolution.
This fit gives the pseudo critical point µc (DHOTRG ) and the
                                                                             As a validation of the numerical strategy, the number density
exponent ν. With the bond dimension DHOTRG ∈ [60, 80],
                                                                             at (U, t) = (8, 0) is computed. The TRG result is consistent
the resulting ν is consistent with the exact value ν = 0.5.
                                                                             with the exact one and the Mott plateau is reproduced around
Extrapolating µc (DHOTRG ) to the limit DHOTRG → ∞, they
                                                                             µ = 4 as expected. The number density is also computed at
used the two fitting ansatz,
                                                                             (U, t) = (80, 1), (8, 1), (2, 1). With (U, t) = (80, 1), the TRG
                  µc (DHOTRG ) = µc + aD−1
                                        HOTRG ,                (136)         gives the smooth number density as a function of the chemical
                                                                             potential and µc /U ̸= 1 in contrast to (U, t) = (8, 0). These
                  µc (DHOTRG ) =     µc + bD−c
                                            HOTRG ,            (137)         results imply that the TRG calculation captures the spatial
                                                                             hopping effects even with the large repulsion parameter U. At
where the latter fit was employed to estimate uncertainty in the
                                                                             (U, t) = (8, 1), (2, 1), the transition point µc is determined by
fitting assumption. These numerical fits are shown in figure 23.
                                                                             the global fit using the quadratic function,
The extrapolation has given µc = 2.642(05)(13) which is in
agreement with the exact value µc = 2.643 · · · based on the
                                                                                                                                         2
Bethe ansatz.                                                                 ⟨n⟩ = 1 + a (µ − µc (DATRG )) + b (µ − µc (DATRG )) , (138)

                                                                        28
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                   Topical Review



with                                                                   Award Number DE-SC0019139 and DE-SC0010113. This
                                                                       research used resources of the Syracuse University HTC
                    µc (DATRG ) = µc + cD−1
                                         ATRG .          (139)         Campus Grid and NSF award ACI-1341006 and the National
                                                                       Energy Research Scientific Computing Center (NERSC), a
DATRG is varied as 56, 64, 72, 80. Their estimates are µc =            U.S. Department of Energy Office of Science User Facility
6.43(4) at (U, t) = (8, 1) and µc = 1.30(6) at (U, t) = (2, 1).        located at Lawrence Berkeley National Laboratory, operated
Based on these estimates, the authors expect that |µc − U/2|           under Contract No. DE-AC02-05CH11231 using NERSC
vanishes only at U = 0, that is the metal–insulator transition         awards HEP-ERCAP0020659 and HEP-ERCAP0023235. S
could take place with any finite repulsion.                            A acknowledges the support from the Endowed Project for
                                                                       Quantum Software Research and Education, the University
7. Conclusions                                                         of Tokyo (https://qsw.phys.s.u-tokyo.ac.jp/) and JSPS
                                                                       KAKENHI Grant Number JP23K13096.
We have reviewed the two formulations to derive the
Grassmann tensor network representations for fermionic path
                                                                       ORCID iDs
integrals. We have shown that both formulations can result in
the same ordinary tensor (bosonic tensor or coefficient tensor)        Shinichiro Akiyama  https://orcid.org/0000-0003-1415-
by properly ordering the auxiliary Grassmann variables. Exact          5620
contractions are defined as the integration of these auxiliary         Yannick Meurice  https://orcid.org/0000-0002-0995-9694
Grassmann variables. These formulations immediately allow              Ryo Sakai  https://orcid.org/0000-0002-0064-2041
us to extend any TRG algorithms for fermions. By introducing
the Grassmann parity functions, we can immediately extend
the algorithms, such as the Levin–Nave TRG and HOTRG, for              References
fermionic path integrals under explicit correspondence with
the original ones. In particular, the geometric representation           [1]   Kadanoff L P 1966 Phys. Phys. Fiz. 2 263–72
                                                                         [2]   Wilson K G and Kogut J B 1974 Phys. Rep. 12 75–199
and the connectivities remain identical.                                 [3]   Wilson K G 1975 Rev. Mod. Phys. 47 773
   These Grassmann TRG algorithms have been applied to                   [4]   Cardy J 1996 Scaling and Renormalization in Statistical
various lattice theories including the relativistic models not                    Physics (Cambridge Lecture Notes in Physics)
only in two but also in four dimensions and the Hubbard model                     (Cambridge University Press)
at finite density. Some of these models include examples                 [5]   Migdal A A 1975 Sov. Phys. JETP 42 743
                                                                         [6]   Kadanoff L P 1976 Ann. Phys. 100 359–94
where the Monte Carlo method is extremely difficult to apply             [7]   Dyson F J 1969 Commun. Math. Phys. 12 91–107
due to the sign problem.                                                 [8]   Baker G A 1972 Phys. Rev. B 5 2622–33
   We have also reviewed several TNR algorithms, the                     [9]   Meurice Y 2007 J. Phys. A: Math. Theor. 40 R39
bond-weighting method, and multilayered Grassmann tensor                [10]   Berges J, Tetradis N and Wetterich C 2002 Phys. Rep.
networks for N f -flavor fermions. Although some of these                         363 223–386
                                                                        [11]   Bervillier C, Juttner A and Litim D F 2007 Nucl. Phys. B
improved TRG methods were originally proposed for spin                            783 213–26
systems, recent numerical calculations have shown that these            [12]   Bervillier C 2013 Nucl. Phys. B 876 587–604
methods are also efficient for fermionic systems. Research on           [13]   White S R 1992 Phys. Rev. Lett. 69 2863–6
such improved algorithms has continued to progress in recent            [14]   Schollwöck U 2005 Rev. Mod. Phys. 77 259–315
years; a new algorithm of loop-TNR [181], combinations of               [15]   Schollwöck U 2011 Ann. Phys., NY 326 96–192
                                                                        [16]   Vidal G 2007 Phys. Rev. Lett. 99 220405
the Monte Carlo method and TRG [182–184], and applica-                  [17]   Cirac J I and Verstraete F 2009 J. Phys. A: Math. Theor.
tion of the machine-learning techniques to the TRG [185–                          42 504004
187]. The extension of these novel improved methods to fer-             [18]   Schollwöck U 2011 Phil. Trans. R. Soc. A 369 2643–61
mionic systems is considered important in examining whether             [19]   Orús R 2014 Ann. Phys. 349 117–58
these methods are also valid for more general physical systems          [20]   Evenbly G and Vidal G 2015 Phys. Rev. Lett. 115 180405
                                                                        [21]   Silvi P, Tschirsich F, Gerster M, Jünemann J, Jaschke D,
including fermions.                                                               Rizzi M and Montangero S 2019 SciPost Phys. Lect. Notes
                                                                                  81
                                                                        [22]   Haegeman J and Verstraete F 2017 Ann. Rev. Condens.
Data availability statement                                                       Matter Phys. 8 355–406
                                                                        [23]   Montangero S 2018 Introduction to Tensor Network
Elements of the data that support the findings of this study                      Methods: Numerical Simulations of low-Dimensional
are openly available at the following URL/DOI: https://github.                    Many-Body Quantum Systems (Springer)
com/akiyama-es/Grassmann-BTRG.                                          [24]   Ran S J, Tirrito E, Peng C, Chen X, Tagliacozzo L, Su G and
                                                                                  Lewenstein M 2020 Tensor Network Contractions:
                                                                                  Methods and Applications to Quantum Many-Body
Acknowledgment                                                                    Systems (Springer)
                                                                        [25]   Cirac J I, Perez-Garcia D, Schuch N and Verstraete F 2021
                                                                                  Rev. Mod. Phys. 93 045003
We thank M Hite and other members of the QuLAT                          [26]   Nishino T and Okunishi K 1996 J. Phys. Soc. Japan 65 891
Collaboration for valuable discussions. Y M was suppor-                 [27]   Levin M and Nave C P 2007 Phys. Rev. Lett. 99 120601
ted in part by the U.S. Department of Energy (DOE) under                [28]   Gu Z C and Wen X G 2009 Phys. Rev. B 80 155131

                                                                  29
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                  Topical Review



 [29] Xie Z Y, Jiang H C, Chen Q N, Weng Z Y and Xiang T 2009             [68] Yu J F, Xie Z Y, Meurice Y, Liu Y, Denbleyker A, Zou H,
         Phys. Rev. Lett. 103 160601                                             Qin M P, Chen J and Xiang T 2014 Phys. Rev. E 89 013308
 [30] Gu Z C, Verstraete F and Wen X G 2010 arXiv:1004.2563               [69] Denbleyker A, Liu Y, Meurice Y, Qin M P, Xiang T, Xie Z Y,
 [31] Gu Z C 2013 Phys. Rev. B 88 115139                                         Yu J F and Zou H 2014 Phys. Rev. D 89 016008
 [32] Xie Z Y, Chen J, Qin M P, Zhu J W, Yang L P and Xiang T             [70] Zou H, Liu Y, Lai C Y, Unmuth-Yockey J, Bazavov A,
         2012 Phys. Rev. B 86 045139                                             Xie Z Y, Xiang T, Chandrasekharan S, Tsai S W and
 [33] Efrati E, Wang Z, Kolan A and Kadanoff L P 2014 Rev. Mod.                  Meurice Y 2014 Phys. Rev. A 90 063603
         Phys. 86 647–67                                                  [71] Akiyama S, Jha R G and Unmuth-Yockey J 2023
 [34] Meurice Y 2013 Phys. Rev. B 87 064422                                      arXiv:2312.11649
 [35] Lyu X, Xu R G and Kawashima N 2021 Phys. Rev. Res.                  [72] Meurice Y 2019 Phys. Rev. D 100 014506
         3 023048                                                         [73] Meurice Y 2020 Phys. Rev. D 102 014506
 [36] Ueda A and Oshikawa M 2021 Phys. Rev. B 104 165132                  [74] Meurice Y, Sakai R and Unmuth-Yockey J 2022 Rev. Mod.
 [37] Ueda A and Oshikawa M 2022 Phys. Rev. E 106 014104                         Phys. 94 025005
 [38] Ueda A and Oshikawa M 2023 Phys. Rev. B 108 024413                  [75] Barthel T, Pineda C and Eisert J 2009 Phys. Rev. A
 [39] Huang C Y, Chan S H, Kao Y J and Chen P 2023 Phys. Rev.                    80 042333
         B 107 205123                                                     [76] Corboz P, Evenbly G, Verstraete F and Vidal G 2010 Phys.
 [40] Guo W and Wei T C 2023 arXiv:2305.09899                                    Rev. A 81 010303
 [41] Ueda A and Yamazaki M 2023 arXiv:2307.02523                         [77] Nishino T, Hieida Y, Okunishi K, Maeshima N, Akutsu Y
 [42] Meurice Y, Perry R and Tsai S W 2011 Phil. Trans. R. Soc. A                and Gendiar A 2001 Prog. Theor. Phys. 105 409–17
         369 2602–11                                                      [78] Gendiar A, Maeshima N and Nishino T 2003 Prog. Theor.
 [43] Byrnes T, Sriganesh P, Bursill R J and Hamer C J 2002 Phys.                Phys. 110 691–9
         Rev. D 66 013002                                                 [79] Verstraete F and Cirac J I 2004 arXiv:cond-mat/0407066
 [44] Bañuls M, Cichy K, Jansen K and Cirac J 2013 J. High               [80] Kraus C V, Schuch N, Verstraete F and Cirac J I 2010 Phys.
         Energy Phys. JHEP11(2013)158                                            Rev. A 81 052338
 [45] Buyens B, Haegeman J, Van Acoleyen K, Verschelde H and              [81] Zohar E and Burrello M 2016 New J. Phys. 18 043008
         Verstraete F 2014 Phys. Rev. Lett. 113 091601                    [82] Zohar E and Cirac J I 2018 Phys. Rev. D 97 034510
 [46] Buyens B, Haegeman J, Verschelde H, Verstraete F and Van            [83] Zapp K and Orús R 2017 Phys. Rev. D 95 114508
         Acoleyen K 2016 Phys. Rev. X 6 041040                            [84] Felser T, Silvi P, Collura M and Montangero S 2020 Phys.
 [47] Funcke L, Jansen K and Kühn S 2020 Phys. Rev. D                            Rev. X 10 041040
         101 054507                                                       [85] Magnifico G, Felser T, Silvi P and Montangero S 2021 Nat.
 [48] Dempsey R, Klebanov I R, Pufu S S and Zan B 2022 Phys.                     Commun. 12 3600
         Rev. Res. 4 043133                                               [86] Felser T, Notarnicola S and Montangero S 2021 Phys. Rev.
 [49] Okuda T 2023 Phys. Rev. D 107 054506                                       Lett. 126 170603
 [50] Honda M, Itou E and Tanizaki Y 2022 J. High Energy Phys.            [87] Emonts P, Bañuls M C, Cirac I and Zohar E 2020 Phys. Rev.
         JHEP11(2022)141                                                         D 102 074501
 [51] Itou E, Matsumoto A and Tanizaki Y 2023 J. High Energy              [88] Robaina D, Bañuls M C and Cirac J I 2021 Phys. Rev. Lett.
         Phys. JHEP11(2023)231                                                   126 050401
 [52] Kühn S, Cirac J I and Bañuls M C 2015 J. High Energy Phys.         [89] Emonts P, Kelman A, Borla U, Moroz S, Gazit S and Zohar E
         JHEP07(2015)130                                                         2023 Phys. Rev. D 107 014505
 [53] Bañuls M C, Cichy K, Cirac J I, Jansen K and Kühn S 2017           [90] Bender J, Emonts P and Cirac J I 2023 Phys. Rev. Res.
         Phys. Rev. X 7 041046                                                   5 043128
 [54] Hayata T, Hidaka Y and Nishimura K 2023                             [91] Cataldi G, Magnifico G, Silvi P and Montangero S 2023
         arXiv:2311.11643                                                        arXiv:2307.09396
 [55] Liu H, Bhattacharya T, Chandrasekharan S and Gupta R                [92] Emonts P and Zohar E 2023 Phys. Rev. D 108 014514
         2023 arXiv:2312.17734                                            [93] Emonts P and Zohar E 2020 SciPost Phys. Lect. Notes 12 1
 [56] Bruckmann F, Jansen K and Kühn S 2019 Phys. Rev. D                  [94] Troyer M and Wiese U J 2005 Phys. Rev. Lett. 94 170201
         99 074501                                                        [95] Prokof’ev N and Svistunov B 2007 Phys. Rev. Lett.
 [57] Tagliacozzo L, Celi A and Lewenstein M 2014 Phys. Rev. X                   99 250201
         4 041024                                                         [96] Chandrasekharan S 2010 Phys. Rev. D 82 025007
 [58] Silvi P, Rico E, Dalmonte M, Tschirsich F and Montangero S          [97] Dornheim T 2019 Phys. Rev. E 100 023307
         2017 Quantum 1 9                                                 [98] Dornheim T 2021 J. Phys. A: Math. Theor. 54 335001
 [59] Silvi P, Sauer Y, Tschirsich F and Montangero S 2019 Phys.          [99] Vanhecke B, Haegeman J, Van Acoleyen K, Vanderstraeten L
         Rev. D 100 074512                                                       and Verstraete F 2019 Phys. Rev. Lett. 123 250604
 [60] Rico E, Pichler T, Dalmonte M, Zoller P and Montangero S           [100] Yang S, Gu Z C and Wen X G 2017 Phys. Rev. Lett.
         2014 Phys. Rev. Lett. 112 201601                                        118 110504
 [61] Pichler T, Dalmonte M, Rico E, Zoller P and Montangero S           [101] Hauru M, Delcamp C and Mizera S 2018 Phys. Rev. B
         2016 Phys. Rev. X 6 011023                                              97 045111
 [62] Zohar E and Burrello M 2015 Phys. Rev. D 91 054506                 [102] Takeda S and Yoshimura Y 2015 Prog. Theor. Exp. Phys.
 [63] Bañuls M C, Cichy K, Cirac J I, Jansen K and Kühn S 2018                  2015 043B01
         PoS (LATTICE 2018) p 022                                        [103] Akiyama S 2023 Phys. Rev. D 108 034514
 [64] Bañuls M C and Cichy K 2020 Rep. Prog. Phys. 83 024401            [104] Bloch J and Lohmayer R 2023 Nucl. Phys. B 986 116032
 [65] Balian R, Drouffe J M and Itzykson C 1975 Phys. Rev. D             [105] Akiyama S, Kuramashi Y, Yamashita T and Yoshimura Y
         11 2104–19                                                              2021 J. High Energy Phys. JHEP01(2021)121
 [66] Itzykson C and Drouffe J 1992 Statistical Field Theory             [106] Kadoh D, Kuramashi Y, Nakamura Y, Sakai R, Takeda S and
         (Cambridge Monographs on Mathematical Physics)                          Yoshimura Y 2018 J. High Energy Phys.
         (Cambridge University Press)                                            JHEP03(2018)141
 [67] Liu Y, Meurice Y, Qin M P, Unmuth-Yockey J, Xiang T,               [107] Asaduzzaman M, Catterall S, Meurice Y, Sakai R and
         Xie Z Y, Yu J F and Zou H 2013 Phys. Rev. D                             Toga G C 2023 arXiv:2312.16167
         88 056005                                                       [108] Dagotto E 1994 Rev. Mod. Phys. 66 763–840

                                                                    30
J. Phys.: Condens. Matter 36 (2024) 343002                                                                                Topical Review



[109] Georges A, Kotliar G, Krauth W and Rozenberg M J 1996             [149] Butt N, Catterall S, Meurice Y, Sakai R and
         Rev. Mod. Phys. 68 13–125                                               Unmuth-Yockey J 2020 Phys. Rev. D 101 094509
[110] Tasaki H 1998 J. Phys.: Condens. Matter 10 4353                   [150] Gattringer C, Kloiber T and Sazonov V 2015 Nucl. Phys. B
[111] LeBlanc J P F et al (Simons Collaboration on the                           897 732–48
         Many-Electron Problem) 2015 Phys. Rev. X                       [151] Göschl D, Gattringer C, Lehmann A and Weis C 2017 Nucl.
         5 041041                                                                Phys. B 924 63–85
[112] Arovas D P, Berg E, Kivelson S A and Raghu S 2022 Annu.           [152] Eckart C and Young G 1936 Psychometrika 1 211–8
         Rev. Condens. Matter Phys. 13 239–74                           [153] Adachi D, Okubo T and Todo S 2022 Phys. Rev. B
[113] Qin M, Schäfer T, Andergassen S, Corboz P and Gull E 2022                  105 L060402
         Annu. Rev. Condens. Matter Phys. 13 275–302                    [154] Nakayama K, Funcke L, Jansen K, Kao Y J and Kühn S 2022
[114] Ostmeyer J 2023 PoS (LATTICE 2022) p 230                                   Phys. Rev. D 105 054507
[115] Akiyama S and Kuramashi Y 2021 Phys. Rev. D                       [155] Adachi D 2020 High-accuracy tensor renormalization group
         104 014504                                                              algorithms and their applications PhD Thesis The
[116] Akiyama S, Kuramashi Y and Yamashita T 2022 Prog.                          University of Tokyo
         Theor. Exp. Phys. 2022 023I01                                  [156] Östlund S and Rommer S 1995 Phys. Rev. Lett. 75 3537–40
[117] Shimizu Y and Kuramashi Y 2014 Phys. Rev. D 90 014508             [157] Dukelsky J, Martín-Delgado M A, Nishino T and Sierra G
[118] Shimizu Y and Kuramashi Y 2014 Phys. Rev. D 90 074503                      1998 Europhys. Lett. 43 457–62
[119] Shimizu Y and Kuramashi Y 2018 Phys. Rev. D 97 034502             [158] Vidal G 2003 Phys. Rev. Lett. 91 147902
[120] Asaduzzaman M, Catterall S, Meurice Y, Sakai R and                [159] Rossi P and Wolff U 1984 Nucl. Phys. B 248 105–22
         Toga G C 2023 J. High Energy Phys. JHEP01(2023)024             [160] Karsch F and Mutter K H 1989 Nucl. Phys. B 313 541–59
[121] Sakai R, Takeda S and Yoshimura Y 2017 Prog. Theor. Exp.          [161] Fromm M 2010 Lattice QCD at strong coupling PhD Thesis
         Phys. 2017 063B07                                                       Zurich, ETH
[122] Yoshimura Y, Kuramashi Y, Nakamura Y, Takeda S and                [162] Milde P, Bloch J and Lohmayer R 2022 PoS (LATTICE 2021)
         Sakai R 2018 Phys. Rev. D 97 054511                                     p 462
[123] Meurice Y 2018 PoS (LATTICE 2018) p 231                           [163] Bloch J, Lohmayer R, Meister M and Nunhofer M 2023
[124] Bao C 2019 Loop optimization of tensor network                             Nucl. Phys. B 987 116107
         renormalization: algorithms and applications PhD Thesis        [164] Bloch J, Jha R G, Lohmayer R and Meister M 2021 Phys.
         U. Waterloo (main)                                                      Rev. D 104 094517
[125] Akiyama S and Kadoh D 2021 J. High Energy Phys.                   [165] Akiyama S 2022 Tensor renormalization group approach to
         JHEP10(2021)188                                                         higher-dimensional lattice field theories PhD Thesis
[126] Akiyama S 2022 J. High Energy Phys. JHEP11(2022)030                        University of Tsukuba
[127] Yosprakob A, Nishimura J and Okunishi K 2023 J. High              [166] Adachi D, Okubo T and Todo S 2020 Phys. Rev. B
         Energy Phys. JHEP11(2023)187                                            102 054432
[128] Akiyama S 2023 arXiv:2311.17691                                   [167] Oba H 2020 Prog. Theor. Exp. Phys. 2020 013B02
[129] Yosprakob A 2023 arXiv:2309.07557                                 [168] Kadoh D and Nakayama K 2019 arXiv:1912.02414
[130] Okunishi K, Nishino T and Ueda H 2022 J. Phys. Soc. Japan         [169] Nakayama K 2023 arXiv:2307.14191
         91 062001                                                      [170] Bender I, Hashimoto T, Karsch F, Linke V, Nakamura A,
[131] Akiyama S, Kuramashi Y and Yoshimura Y 2022 PoS                            Plewnia M, Stamatescu I O and Wetzel W 1992 Nucl.
         (LATTICE 2021) p 530                                                    Phys. B 26 323–5
[132] Kadoh D 2022 PoS (LATTICE 2021) p 633                             [171] Buballa M 2005 Phys. Rep. 407 205–376
[133] Kuramashi Y and Yoshimura Y 2020 J. High Energy Phys.             [172] Aoki K I, Kumamoto S I and Yamada M 2018 Nucl. Phys. B
         JHEP04(2020)089                                                         931 105–31
[134] Fukuma M, Kadoh D and Matsumoto N 2021 Prog. Theor.               [173] Witten E 1982 Nucl. Phys. B 202 253
         Exp. Phys. 2021 123B03                                         [174] Catterall S, Kaplan D B and Unsal M 2009 Phys. Rep.
[135] Hirasawa M, Matsumoto A, Nishimura J and Yosprakob A                       484 71–130
         2021 J. High Energy Phys. JHEP12(2021)011                      [175] De Lathauwer L, De Moor B and Vandewalle J 2000 SIAM J.
[136] Dittrich B, Mizera S and Steinhaus S 2016 New J. Phys.                     Matrix Anal. Appl. 21 1253–78
         18 053009                                                      [176] De Lathauwer L, De Moor B and Vandewalle J 2000 SIAM J.
[137] Kuwahara T and Tsuchiya A 2022 Prog. Theor. Exp. Phys.                     Matrix Anal. Appl. 21 1324–42
         2022 093B02                                                    [177] van den Doel C and Smit J 1983 Nucl. Phys. B 228 122–44
[138] Luo X and Kuramashi Y 2023 Phys. Rev. D 107 094509                [178] Lieb E H and Wu F Y 1968 Phys. Rev. Lett. 20 1445–8
[139] Akiyama S and Kuramashi Y 2022 J. High Energy Phys.               [179] Metzner W and Vollhardt D 1989 Phys. Rev. Lett. 62 324–7
         JHEP05(2022)102                                                [180] Müller-Hartmann E 1989 Z. Phys. B: Condens. Matter
[140] Morita S, Igarashi R, Zhao H H and Kawashima N 2018                        74 507–12
         Phys. Rev. E 97 033310                                         [181] Homma K and Kawashima N 2023 arXiv:2306.17479
[141] Wang L and Verstraete F 2011 arXiv:1110.4362                      [182] Ferris A J 2015 arXiv:1507.00767
[142] Corboz P, Rice T M and Troyer M 2014 Phys. Rev. Lett.             [183] Huggins W, Freeman C D, Stoudenmire M, Tubman N M
         113 046402                                                              and Whaley K B 2017 arXiv:1710.03757
[143] Iino S, Morita S and Kawashima N 2019 Phys. Rev. B                [184] Arai E, Ohki H, Takeda S and Tomii M 2023 Phys. Rev. D
         100 035449                                                              107 114515
[144] Zhao H H, Xie Z Y, Chen Q N, Wei Z C, Cai J W and                 [185] Liao H J, Liu J G, Wang L and Xiang T 2019 Phys. Rev. X
         Xiang T 2010 Phys. Rev. B 81 174411                                     9 031041
[145] Evenbly G 2017 Phys. Rev. B 95 045117                             [186] Chen B B, Gao Y, Guo Y B, Liu Y, Zhao H H, Liao H J,
[146] Orús R and Vidal G 2008 Phys. Rev. B 78 155117                            Wang L, Xiang T, Li W and Xie Z Y 2020 Phys. Rev. B
[147] Wolff U 2020 Nucl. Phys. B 955 115061                                      101 220409
[148] Gausterer H and Lang C B 1995 Nucl. Phys. B                       [187] Jha R G and Samlodia A 2024 Comput. Phys. Commun.
         455 785–95                                                              294 108941




                                                                   31
