Higher order tensor renormalization group for relativistic fermion systems

We apply the higher order tensor renormalization group to two and three dimensional relativistic fermion systems on the lattice. In order to perform a coarse-graining of tensor networks including Grassmann variables, we introduce Grassmann higher order tensor renormalization group. We test the validity of the new algorithm by comparing its results with those of exact or previous methods.


Introduction
Monte Carlo simulations of lattice gauge theories have been shown successful as a nonperturbative numerical approach since the celebrated formulation by Wilson [1] and the first simulation by Creutz [2]. Such a stochastic approach, however, is generally suffering from the sign problem when the Boltzmann weight is complex. For example, finite quark density systems, θ term included systems, or chiral gauge theories are not easily accessible due to the sign problem. To avoid the problem, one may rely on other methods using a deterministic algorithm, say the tensor renormalization group (TRG) [3].
The basic procedure of the TRG is as follows. A target quantity which the TRG can compute is a partition function in the path integral representation. For that purpose, first of all, one has to rewrite a partition function in terms of a tensor network. In general, this step can be done by expanding the Boltzmann weight with some proper expansion scheme depending on the type of physical degrees of freedom (non-compact boson fields, compact boson fields, or fermion fields). And then after carrying out the integration of the original fields, it turns out that tensors live at lattice sites and the partition function is expressed as a contraction of them. The number of terms in the summation is exponentially large as a function of the system size. To avoid such an expensive computational cost, one may rely on a coarse-graining of the tensor network. In this step, the singular value decomposition is used to reduce the d.o.f of the system while preserving important information. By repeating local blocking transformations, one can approximately compute a value of the partition function. Actually, thanks to the coarse-graining algorithm, the cost now is proportional to the logarithm of system size.
Since the original idea of TRG was introduced by Levin and Nave [3], it has been applied to many two dimensional models associated with elementary particle physics: the lattice φ 4 model [4], the lattice Schwinger model [5,6], and the finite density lattice Thirring model [7].
In the studies of the relativistic fermion systems, the Grassmann TRG (GTRG) proposed by Gu et al. [8,9] was used. On the other hand, although the original idea of TRG was limited to two dimensional systems, a new coarse-graining method suited for any higher dimensional system was proposed in Ref. [10]. The new method uses the higher order singular value decomposition; thus it is called higher order tensor renormalization group (HOTRG). The authors in Ref. [11] analyzed the three dimensional Potts models by using the HOTRG and obtained the critical temperature and exponents precisely. The HOTRG was also applied in two dimensional systems: the XY model [12], the O(3) model [13], and the CP(1) model [14].
It is natural to extend the idea of TRG to higher dimensional systems and more complicated systems. Our final target is lattice QCD which is a four dimensional relativistic field 2 theory of quarks and gluons. The tensor network representation of non-Abelian gauge theories was already attempted in Ref. [15], and the HOTRG for boson systems was already done as mentioned above. Therefore, the final missing piece to achieve the goal is to treat fermions in higher dimensional systems. In this paper, we formulate HOTRG for fermions applicable to any dimensional system, and we call this Grassmann higher order tensor renormalization group (GHOTRG). As a concrete example, we will provide some details of GHOTRG in a two dimensional system 1 .
This paper is organized as follows. In Sec. 2 we introduce GHOTRG for a two dimensional fermion system with a focus on the Grassmann part in the tensor network. In Sec. 3 we show numerical results and compare them with exact or previous ones. A summary and an outlook for the future work are given in Sec. 4.

Grassmann higher order tensor renormalization group
In this section, after briefly reviewing the original HOTRG, we explain GHOTRG in two dimensions in detail. The three dimensional version will be obtained straightforwardly. Finally we discuss how to treat the anti-periodic boundary conditions for fermion fields.
Lattice units a = 1 are assumed in the following.

Model and notation in two dimensional systems
The Lagrangian density of the free Wilson fermion 2 in two dimensions is given by where n = (n 1 , n 2 ) is the lattice coordinate, the fermion fields ψ,ψ have two spinor components, and m and µ denote the mass and the chemical potential. Following the prescription described in Refs. [5,7], one can obtain tensor network representation of the partition 1 The extension to four dimensional systems is straightforward although both memory and computational costs are extremely demanding. 2 Inclusion of interaction terms is straightforward as seen in [5] and [7] for the Schwinger model and the Thirring model, respectively.
The tensor T , an elementary building block of the tensor network, is defined as (6) · η n+1,1 η n,1 x n,1 η n,2 η n+1,2 x n,2 ξ n+2,1 ξ n,1 t n,1 ξ n,2 ξ n+2,2 t n,2 is called Grassmann part. Note that the indices of tensor in Eq. (4) should be read as t n =(t n,1 , t n,2 ), reflecting the fact that the original Grassmann variable ψ has two components. Each index runs from 0 to 1; thus the tensor has 2 2×4 elements at the initial stage. Due to the Grassmann nature, the element of the tensor takes a non-trivial value only when the indices satisfy 2 i=1 x n,i + t n,i + x n−1,i + t n−2,i mod 2 = 0, for all n, otherwise the element is zero. In the following we discuss the renormalization procedure for the bosonic tensor and the Grassmann part separately.

Normal HOTRG procedure for the bosonic tensor
In this subsection, let us see the renormalization of bosonic tensors, especially focusing on the coarse-graining along the1-direction. This is a brief review of the normal HOTRG [10].
First we consider the new tensor M by contracting two bosonic tensors placed next to each other along the1-direction (see Fig. 1) where the integrated indices are defined as Then one defines M + and M − as with Next we apply the eigenvalue decomposition to M + and obtain a unitary matrix U + and the eigenvalues λ + : where the new index t n * ,b will be regarded as the second component of the index of the new tensor T new,b as it will be clear soon (see Eq. (18)). Similarly one obtains U − and λ − from M − , and then we define ǫ + and ǫ − as with a given D cut one can choose. ǫ ± represent an amount of truncation error and are used to select a unitary matrix which maintains better precision. That is, if ǫ + < ǫ − , then U + is adopted and vice versa. Now one obtains the new tensor by using the selected unitary matrix (denoted by U) and restricting the new indices 1 ≤ t n * ,b , t n * −2,b ≤ D cut , Actually, when dealing with the Grassmann variables, this part should be modified in order to incorporate the sign factors originated from the coarse-graining of them. See Eq. (31) for details. Fig. 1 Coarse-graining along the1-direction.
The indices of the new tensor now contain two kinds of coordinates n and n * . At this point, we make clear the relationship between n and n * as shown in Fig. 2. Following the relationship, the indices are changed as follows: In the last step, n * has been renamed n since all coordinates are labeled by n * . To remember that the lattice spacing of the1-direction is doubled, however, the new unit vector defined as1 * =1 +1 will be used.

Coarse-graining of the Grassmann part
In this subsection, we discuss the coarse-graining of the Grassmann part. Our objective here is to coarse-grain the Grassmann parts of the two tensors T xntnx n−1 t n−2 and T x n+1 t n+1 xnt n+1−2 , which are located next to each other along the1-direction, and obtain a new Grassmann part with new Grassmann variables.
6 Fig. 2 Correspondence between n and n * .
For that purpose, first of all, we collect their Grassmann parts where the additional sign factors arise when breaking the hopping factors for ξ. x n,1( x n,1 +x n,2) + t n,1 +t n+1,1 +t n,1 t n,2 +t n+1,1 t n+1,2 · dη shown in Fig. 3, one can execute the integration of the old Grassmann variable ξ: x n,1( x n,1 +x n,2) +t n,2( t n,1 +t n,2) +t n+1,2 t n+1,1 +t n+1,2 · dη where the hopping factors for the new Grassmann variable have been shifted n * → n * +2. Note that new constraints including the new indices hold: This constraint is a consequence of Eq. (10) and Eq. (26). From a practical point of view, it is convenient to explicitly multiply the factor to the new tensor.
One can include the phase factors when contracting M and U consist of purely bosonic T in Eq. (11) or include them when making M. We observe that the latter shows a slightly better approximation. In Sec. 3, we show the results in the latter case. For later use, we summarize the new bosonic tensor including the constraint in Eq. (30), T new where the index t in the LHS should be understood as t n * = (t n * ,f , t n * ,b ) and t n * −2 = (t n * −2,f , t n * −2,b ) while x n and x n−1 * are the same as before in Eq. (7).
The ordering of the Grassmann variables and the index assignment are consistent with that of the initial tensor in Eq. (4).
The Grassmann part of the above form is similar to that of the initial tensor in Eq. (4).
The second half of the coarse-graining along the2-direction can be done as follows. First, let us interchange the1-and the2-axis and obtain a flipped tensor By applying a similar coarse-graining procedure for the1-direction to this tensor, one can carry out the coarse-graining for the2-direction. The coarse-grained tensor is then obtained by T xntnx n−1 * t n−2 * = T xntnx n−1 * t n−2 * dη where the all indices of the tensor have the structure x n = (x n,f , x n,b ). The sign factor and the constraint have been already included in the bosonic tensor. Now, the scale of the1direction is equal to that of the2-direction. We define an "iteration" of coarse-graining which consists of the coarse-graining for both directions. Then Eq. (36) is recognized as the coarsegrained tensor after the first iteration as well as the starting tensor of the second iteration. Actually, it turns out that from the second iteration, the structure of tensor in Eq. (36) does not change. The actual coarse-graining procedure from the second iteration can be done by replacing the second components of the indices to zero in the first iteration. For example, for the1-direction one only has to replace x n,2 → 0, t n,2 → 0, for all n (37) in the Grassmann part.
The GHOTRG algorithm presented here can be immediately extended to any dimensional system. We do not show any details of a coarse-graining for higher dimensional systems, but we present an initial tensor for the three dimensional free Wilson fermion with two spinor components in App. A.

Anti-periodic boundary conditions
If anti-periodic boundary conditions are imposed, the sign factor (−1) arises when a fermion line passes the boundary for the2-direction. Therefore, one can realize anti-periodic boundary conditions on the tensor network by inserting a boundary tensor B into each link along the2-direction as shown in Fig. 4, This insertion affects only the tensors around the boundary. Even when there is a boundary, the coarse-graining along the2-direction is simple. As shown in Fig. 5, there is no difficulty for the coarse-graining as long as the number of lattice sites for this direction is even. The coarse-graining along the1-direction is also straightforward although one needs to care about the additional sign factor in B. In the coarse-graining step involving the boundary shown in Fig. 6, as a result of the insertion of B, a sign factor (−1) t n−2,1 +t n−2,2 +t n+1−2,1 +t n+1−2,2 (40) appears. By using the new index t n * −2,f in Eq. (26), the factor is rewritten as (−1) t n−2,1 +t n−2,2 +t n+1−2,1 +t n+1−2,2 = (−1) t n−2,1 +t n−2,2 +t n+1−2,1 +t n+1−2,2 mod 2 = (−1) t n * −2,f .
12 Therefore we define a new boundary tensor B new for the next coarse-graining, From the second iteration, the structure of boundary tensor does not change. In this way, the boundary effect is not involved in the coarse-graining procedure for the tensor T in the bulk and appears only in the contraction procedure to compute the partition function.

Fig. 5
Coarse-graining along the2-direction in the presence of the boundary. There is no difficulty for the coarse-graining as long as the number of sites for the2-direction is even.

Fig. 6
Coarse-graining along the1-direction. For the coarse-graining involving the boundary, the additional sign factors originated from boundary tensors should be taken care. 13 3 Numerical results

The error of the free energy and the hierarchy of eigenvalues
Here we exclusively consider the massless free Wilson fermion. Figure 7 shows the free energy as a function of the chemical potential µ on the 2 × 2 space-time lattice computed by using the GHOTRG. For the maximal D cut (= 16), it agrees with the analytical exact result up to the machine precision. On the other hand, when D cut decreases, the difference becomes larger especially in the region µ < 1.0.
For larger D cut , the error rapidly reduces as expected. In the figure, we observe that the error becomes larger for µ ≈ 0.0 while it is smaller for µ ≈ 2.0. This tendency is seen for all D cut < 16. In order to qualitatively understand this behavior, we investigate the hierarchy of eigenvalues of M ± 5 . We show them at µ = 0.0 and µ = 2.0 in Fig. 9. For µ = 0.0, the hierarchy is milder than that of µ = 2.0. Therefore we confirm the strong correlation between the error of the free energy and the hierarchy of eigenvalues; a low compression in the coarse-graining procedure causes a large error in a physical quantity.
We also investigate the error for a larger lattice volume 32 × 32 as shown in Fig. 10. The error tends to be large around µ = 1.0. The eigenvalues of M ± is shown in Fig. 11 where the hierarchy at µ = 1.0 is quite poor compared with that of µ = 2.0. This behavior is very similar to that of GTRG [7]. Figure 12 shows the fermion number density

The fermion number density of the lattice Thirring model
for the lattice Thirring model whose tensor network representation is presented in Ref. [7].
In the large µ region, the fermion number density reaches the saturation density. When switching on the interaction g = 0, the Silver Blaze like phenomenon seems to occur around the small µ region where the fermion number density is constantly zero while the onset is observed at finite µ. Such a behavior was already seen in Ref. [7] where GTRG is used. In any case, the exact massless free results are well reproduced by GHOTRG with D cut = 32 The hierarchy of eigenvalues of M ± at µ = 0.0 and µ = 2.0 with D cut = 16 for massless free system. An anisotropic coarse-graining along the1-direction is called 0.5-iteration. in a wide range of µ although the approximation gets worse around µ = 1.0; this tendency was also observed in the GTRG analysis [7].

The free energy
We apply GHOTRG to the three dimensional massless free Wilson fermion and calculate the free energy. Since it is very demanding to carry out the full D cut calculation even for the 2 × 2 × 2 lattice, we have to compromise to use an anisotropic lattice, say, 2 × 1 × 1 to check the algorithm and the code. As a check, we carry out the coarse-graining from 2 × 1 × 1 to 1 × 1 × 1 with several D cut as shown in Fig. 13. With the full D cut = 16, the GHOTRG result is consistent with the exact result up to the machine precision. Even for smaller D cut < 16, the accuracy is reasonable, and for larger D cut the free energy rapidly converges to the exact one. Note that at the smaller D cut , the free energy may be complex due to the truncation in the coarse-graining procedure; thus the real part of the free energy is shown in the figure.  Fig. 12 The fermion number density for the lattice Thirring model. The exact results for massless free case is also shown. The lattice volume is 32 × 32.  Fig. 13 The free energy of three dimensional massless free Wilson fermions as a function of µ on 2 × 1 × 1 lattice. The real part of the free energy is shown. 18 4 Summary and outlook In this paper we have formulated Grassmann higher order tensor renormalization group and applied it to the free fermion system and the lattice Thirring model. The numerical results are consistent with analytical or previous ones. Thus we conclude that GHOTRG is a correct algorithm.
In the study of two dimensional finite chemical potential systems, we observed the very poor hierarchy of eigenvalues at µ = 1.0, and this is the reason why the accuracy gets worse around this parameter region. This shows that the situation of accuracy for GHOTRG is similar to that of the GTRG.
The algorithm of GHOTRG can be straightforwardly extended to three dimensions or more. In principle, now one can deal with any dimensional fermion system by using GHOTRG. It is, however, still hard to study four dimensional fermion systems using GHOTRG owing to a huge computational cost. We hope that our formulation is a starting point towards the four dimensional lattice QCD. Apart from lattice QCD, there are some interesting models in low dimensions such as lattice chiral gauge theories and lattice SUSY which are generally suffering from the sign problem. We are now in the position to approach such models in two or three dimensions by using GHOTRG.
The bottleneck of (G)HOTRG is a tensor contraction whose cost is proportional to D 4d−1 cut with the dimensionality d. This cost could be reduced by randomizing the contraction. So far, many random algorithms have been proposed in the frameworks of MPS, PEPS, and TTN [16][17][18]. Recently, new algorithms combined with MERA are proposed [19,20]. A tensor network Monte Carlo method introduced by Ferris [21] may be applied to the HOTRG. It seems worthwhile to pursue this direction.