Gradient optimization of fermionic projected entangled pair states on directed lattices

TL;DR

This paper extends a stochastic gradient optimization method to fermionic PEPS, introducing a fermi arrow notation, and demonstrates its effectiveness in simulating complex fermionic models with improved accuracy and scalability.

Contribution

The authors develop a fermionic PEPS gradient optimization method with a new fermi arrow notation, enabling efficient simulation of fermionic systems with larger bond dimensions.

Findings

Gradient optimization improves simple update results.

Larger bond dimensions are needed for convergence.

Method offers lower scaling than direct contraction methods.

Abstract

The recently developed stochastic gradient method combined with Monte Carlo sampling techniques [PRB {\bf 95}, 195154 (2017)] offers a low scaling and accurate method to optimize the projected entangled pair states (PEPS). We extended this method to the fermionic PEPS (fPEPS). To simplify the implementation, we introduce a fermi arrow notation to specify the order of the fermion operators in the virtual entangled EPR pairs. By defining some local operation rules associated with the fermi arrows, one can implement fPEPS algorithms very similar to that of standard PEPS. We benchmark the method for the interacting spinless fermion models, and the t-J models. The numerical calculations show that the gradient optimization greatly improves the results of simple update method. Furthermore, much larger virtual bond dimensions ($D$) and truncation dimensions ($D_c$) than those of boson and spin systems are necessary to converge the results. The method therefore offer a powerful tool to simulate fermion systems because it has much lower scaling than the direct contraction methods.

Tables

Table 1. Table 1: The ground state energies of the interacting spinless fermion models on a 4 × 4 4 4 4\times 4 lattice. The ground state energies obtained from the simple update (SU) method, and the gradient optimization (GO) method are compared with those from exact diagonalization method. The relative errors between the GO results and exact results are in the order of 10 -5 .

$V$	SU	GO	Exact	relative error
0.1	-0.66590	-0.67124	-0.67125	1$\times$10-5
0.8	-0.59255	-0.59309	-0.59312	5$\times$10-5
2	-0.48136	-0.50643	-0.50646	5$\times$10-5

Table 2. Table 2: Compare the ground state energies of the free fermion model of SU ( D 𝐷 D =8), and GO ( D 𝐷 D =6, 8) with the exact results. For the 4 × \times 4 lattice, the SU result is numerically instable for D 𝐷 D =8.

	SU(D=8)	GO(D=6)			GO(D=8)
Size	Energy	Energy	$D_c$	relative error	Energy	$D_c$	relative error	exact
4$\times$4	-	-0.68398	16	4$\times$10-5	-0.68401	32	1$\times$10-5	-0.68402
6$\times$6	-0.67721	-0.73269	24	5$\times$10-4	-0.73305	52	5$\times$10-5	-0.73309
8$\times$8	-0.74763	-0.75414	40	1$\times$10-3	-0.75492	84	2$\times$10-4	-0.75510
10$\times$10	-0.75387	-0.76619	55	1$\times$10-3	-0.76705	110	5$\times$10-4	-0.76748
12$\times$12	-	-0.77094	80	5$\times$10-3	-	-	-	-0.77538

Table 3. Table 3: Compare the ground state energies of t-J model with hole filling n ¯ h subscript ¯ 𝑛 ℎ \bar{n}_{h} =0.125 calculated by SU and GO methods. A virtual bond dimension D = 12 𝐷 12 D=12 is used.

size	SU	GO
4$\times$4	-0.55108	-0.56420
4$\times$8	-0.57994	-0.59055
6$\times$8	-0.59431	-0.60349
8$\times$8	-0.60849	-0.61184
8$\times$10	-0.61068	-0.61738
8$\times$12	-0.61707	-0.62164
12$\times$12	-0.62307	-0.62973
$L\to \infty$	-0.66757	-0.67008

Equations

$${\hat{I}}_h(i,j)|0⟩=\sum _{k=0}^{D-1}|k{⟩}_{(i,j{)}_r}|k{⟩}_{(i,j+1{)}_l},$$

$$|{\Phi }_0⟩={\Pi }_{i=1}^{L_1-1}{\Pi }_{j=1}^{L_2-1}{\hat{I}}_h(i,j){\hat{I}}_v(i,j)|0⟩.$$

$$\hat{P}[i,j]$$

$$|{\Phi }_{\mathrm{fPEPS}}⟩={\Pi }_{i,j}\hat{P}[i,j]|{\Phi }_0⟩.$$

$$\tilde{p}({\beta }_3,{\beta }_4)$$

$$\sum _{k=0}^{D-1}{a}_A^{†k}{a}_B^{†k}|0⟩=\sum _{k=0}^{D-1}(-1{)}^{\tilde{p}(k)}{a}_B^{†k}{a}_A^{†k}|0⟩.$$

$$\sum _{\beta ,{\beta }_A,{\beta }_B}{A}_{{\beta }_A}{a}_A^{{\beta }_A}{B}_{{\beta }_B}{a}_B^{{\beta }_B}{a}_B^{†\beta }{a}_A^{†\beta }|0⟩=\sum _{\beta }{A}_{\beta }{B}_{\beta }|0⟩$$

$$\sum _{\beta ,{\beta }_A,{\beta }_B}{A}_{{\beta }_A}{a}_A^{{\beta }_A}{B}_{{\beta }_B}{a}_B^{{\beta }_B}{a}_A^{†\beta }{a}_B^{†\beta }|0⟩=\sum _{\beta }\tilde{p}(\beta )A_{\beta }{B}_{\beta }|0⟩$$

$$R_{I,J}=\sum _{i_1,i_2,\cdots ,i_N}\prod _{l=1,k}\tilde{p}(i_l)A_{I,i_N,i_{N-1},\cdots ,i_1}{B}_{i_1,i_2,\cdots ,i_N,J}$$

$$|\Psi ⟩$$

$$E$$

$$E(\mathbf{n})=\sum _{{\mathbf{n}}^{\prime }}\frac{W({\mathbf{n}}^{\prime })}{W(\mathbf{n})}⟨{\mathbf{n}}^{\prime }|H|\mathbf{n}⟩.$$

$$\frac{\partial E}{\partial T[\mathbf{i}{]}_{n,{\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4}}=2⟨\Delta [\mathbf{i}{]}_{n,{\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4}(\mathbf{n})E(\mathbf{n})⟩$$

$$-2⟨\Delta [\mathbf{i}{]}_{n,{\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4}(\mathbf{n})⟩⟨E(\mathbf{n})⟩,$$

$$\Delta [\mathbf{i}{]}_{n,{\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4}(n)=\frac{1}{W(\mathbf{n})}\frac{\partial W(\mathbf{n})}{\partial T[\mathbf{i}{]}_{n,{\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4}},$$

$$\frac{\partial W({\bf n})}{\partial T[{\bf i}]_{n,\beta_{1},\beta_{2},\beta_{3},\beta_{4}}}=\mathcal{C}\Big{[}T[1]\cdots T[{\bf i-1}]T[{\bf i+1}]\cdots T[N]\Big{]}\,.$$

$$\delta =\frac{||X⟩-|\Psi ⟩{|}^2}{⟨\Psi |\Psi ⟩},$$

$$⟨X^{i,j}|X⟩=⟨X^{i,j}|\Psi ⟩,$$

$$H=-t\sum _{⟨i,j⟩}(c_i^{†}{c}_j+H.c.)+V\sum _{⟨i,j⟩}{n}_i{n}_j,$$

$$H=-t\sum _{⟨i,j⟩,\sigma }(c_{i,\sigma }^{†}{c}_{j,\sigma }+H.c.)+J\sum _{⟨i,j⟩}({\vec{S}}_i{\vec{S}}_j-\frac{1}{4}{n}_i{n}_j)$$

$$\Delta E=E(D_c)-E(D_c^{\max }),$$

$$\hat{A}$$

$$\hat{B}$$

$$\widehat{A^{\prime }}$$

$$\widehat{B^{\prime }}$$

$$\hat{C}=\hat{U}\hat{S}\hat{V}{\hat{I}}_{\hat{U}\to \hat{S}}{\hat{I}}_{\hat{S}\to \hat{V}},$$

$$\hat{C}$$

$$\hat{U}$$

$$\hat{S}$$

$$\hat{V}$$

Full text

I introduction

Interacting quantum many-body systems pose some of the most exciting open problems in physics. Particularly, fermion systems are central to many of the most fascinating effects in condensed matter physics, such as high-temperature superconductivity,Lee et al. (2006) the fractional quantum Hall effect,Stormer et al. (1999) and Mott insulator transitions Edwards and Hewson (1968); Imada et al. (1998). The simulation of the strongly correlated fermion system plays the critical role to understand these system and is also one of the most challenging problems in condensed matter physics.

The Quantum Monte Carlo (QMC)Foulkes et al. (2001) method as one of the leading methods in studying many-body physics has achieved great success in bosonic and spin systems since its first proposed. However, except in some special cases, Li and Yao (2018) the fermion systems are extremely difficult to treat using QMC simulationsLoh et al. (1990); Troyer and Wiese (2005) because of the notorious “sign problems”.

Recently, the methods based on tensor network states (TNS), especially the projected entangled states (PEPS) Schollwöck (2011); Perez-Garcia et al. (2007); Verstraete et al. (2008); Jiang et al. (2008); Vidal (2008); Verstraete and Cirac (2004); Sfondrini et al. (2010); Verstraete et al. (2006) have shown their power on simulation of the strongly correlated many-particle systems. The PEPS is sign-problem free and has achieved great successes in studying the frustrated spin models Wang et al. (2016); Vidal (2007); M.-H. et al. (2012). The PEPS method has been extended to study fermion models (namely fPEPS) by different approaches Barthel et al. (2009); Corboz et al. (2010, 2011); Gu et al. (2010, 2013); Kraus et al. (2010). Apperently, the fPEPS are more complicated than PEPS because of the anti-commutation properties of the fermion operators. In addition, fermion systems are highly frustrated. It has been proven that the entanglements of the ground states of some fermion systems are beyond the area law Hastings (2007); Wolf (2006). Therefore, to faithfully simulate such models, it usually requires very large bond dimensions ($D$). Furthermore, it has been shown that the imaginary time evolution with simple update Jiang et al. (2008) method may have large errors because the environment effects are oversimplified. To exactly consider the environment, the traditional methods, e.g., the full update method Lubasch et al. (2014a, b), suffer from extremely high computational scaling to the bond dimensions. This problem is more serious for the fermion models when large $D$ is required. We note that the recently developed infinite PEPS (iPEPS) with full update method has achieved great success,Corboz et al. (2014, 2010) by making use of the translation symmetry, which may greatly reduce the number of independent tensors. However, not all systems have such symmetry, e.g., defects, disorders and systems with spontaneous symmetry broken, etc. In these cases, the finite PEPS method is essential.

The recently developed Monte Carlo sampling techniques for PEPS can greatly reduce the computational scalingSandvik and Vidal (2007); Schuch et al. (2008); Sfondrini et al. (2010); Bañuls et al. (2017); Dong et al. (2017); Liu et al. (2017). By combining with stochastic gradient optimization (GO) method, one can achieve great precision in obtaining the ground states. Liu et al. (2017); He et al. (2018) In this work, we extended the stochastic gradient method Liu et al. (2017); He et al. (2018) to optimize the fPEPS wave functions for fermion systems. To simplify the implementation of the fPEPS algorithms, we introduce a “Fermi arrow” notation to specify the order of the fermion operators in the entangled EPR pairs. With this notation and some local operation rules associated with the Fermi arrows, we can greatly simplify the implementation of the stochastic gradient optimization method (and other methods) for fPEPS. We implement this local operation rules for fPEPS in our recently developed TNSpack Dong et al. (2018), in which the anti-commutation properties of the fermion operators are automatically taken account of. Therefore, one can implement fPEPS algorithm very similar to that of the standard PEPS without worry too much about the details of the anti-commutation between the fermion operators.

We benchmark the stochastic gradient method for fPEPS on the interacting spin-less fermion models, and the t-J models. The numerical calculations show that the gradient optimization greatly improves the results of simple update method. Furthermore, for these models, very large virtual bond dimensions $D$ and truncation dimensions $D_{c}$ are necessary to converge the results which is the dominate difficult to simulate the fermion systems. Therefore the present method is advantageous because it has much lower scaling than the traditional direct contraction method.

II Definition of fPEPS based on directed network

The definition of the fPEPSKraus et al. (2010) on a lattice is similar to that of the standard PEPSVerstraete and Cirac (2004); Verstraete et al. (2006). Without lose of generality, we take a fermion system on a $L_{1}\times L_{2}$ square lattice as an example, where the physical dimension of each site is $d$. For each horizontal bond connecting sites $(i,j)$ and $(i,j+1)$, there is a EPR pair, i.e., a Bell type entangled state,

[TABLE]

where $|k\rangle_{(i,j)_{r}}$ and $|k\rangle_{(i,j+1)_{l}}$ are the fermion states on site $(i,j)$ and site $(i,j+1)$. States $|k\rangle$ are generated as, $|k\rangle=|k_{1}k_{2}\cdots k_{p}\rangle={a_{1}^{{\dagger}k_{1}}}{}{a_{2}^{{\dagger}k_{2}}}\cdots{a_{p}^{{\dagger}k_{p}}}|0\rangle$, where $(k_{1}k_{2}\cdots k_{p}$) is the binary representation of $k$ and $|0\rangle$ is the vacuum state. $a_{i}$s and $a^{\dagger}_{j}$s are the fermion operators that satisfy $\{a_{i},a^{\dagger}_{j}\}=\delta_{ij}$. For convenience, we denote the state $|k_{(i,j)_{r}}\rangle={a_{(i,j)_{r}}^{{\dagger}k}}|0\rangle$. Similarly, for each vertical bond connecting site $(i,j)$ and $(i+1,j)$, there is also a Bell type entangled state, (in short) $\hat{I}_{v}(i,j)|0\rangle=\sum^{D-1}_{k=0}{a_{(i,j)_{d}}^{{\dagger}k}}{a_{(i+1,j)_{u}}^{{\dagger}k}}|0\rangle$. Therefore, a standard virtual mother state of a fPEPS can be defined as,

[TABLE]

To define a quantum state in the real physical space, we project $|\Phi_{0}\rangle$ to the physical space. The projector on site $(i,j)$ is defined as:

[TABLE]

Here, $a^{{\dagger}}_{(i,j)}$ is the creation operator of the physical particle on site $(i,j)$ whereas $a_{(i,j),m}^{n}$ ($m=l,d,r,u$ and $n=0,1,\cdots,D-1$) are the annihilation operators of the state $|n\rangle_{(i,j),m}$. The fPEPS is then defined as,

[TABLE]

To make the fPEPS well defined, the state $|\Phi_{\rm fPEPS}\rangle$ should be independent of the order of the projectors up to a global phase, i.e, the parity of all elements in a projector should be the same. The parity of the element $T_{\beta,\beta_{1},\beta_{2},\beta_{3},\beta_{4}}[i,j]a^{{\dagger},\beta}_{(i,j)}a_{(i,j)_{l}}^{\beta_{1}}a_{(i,j)_{d}}^{\beta_{2}}a_{(i,j)_{r}}^{\beta_{3}}a_{(i,j)_{u}}^{\beta_{4}}$ of the projector $\hat{P}[i,j]$ is obtained by $\tilde{p}(\beta)\tilde{p}(\beta_{1})\tilde{p}(\beta_{2})\tilde{p}(\beta_{3})\tilde{p}(\beta_{4})$, where $\tilde{p}(x)$=-1, if the parity of $x$ is odd, and $\tilde{p}(x)$=+1 if the parity of $x$ is even. Therefore, the parity of all elements can be obtained by the lower indices of tensor $T_{\beta,\beta_{1},\beta_{2},\beta_{3},\beta_{4}}[i,j]$. Without lose of generality, we assume all nonzero projector elements have even parity in this paper. As a consequence, the elements with odd parity vanish, i.e., $T_{\beta,\beta_{1},\beta_{2},\beta_{3},\beta_{4}}[i,j]$=0, if $\beta$+$\beta_{1}$+$\beta_{2}$+$\beta_{3}$+$\beta_{4}$ is odd. In this definition of fPEPS, we may interchange the positions of any two projectors and EPR pairs, because they all have even parity.

One of the key issues in the fPEPS is the order of the fermion operators, including the operators in the projectors and in EPRs. We define the standard order of the fermion operators in each projector operators on the square lattice as followings, physical creation operator, left, down, right, and up virtual operators (i.e., anti-clockwise order), which is the same as the order of the lower indices in the tensor $T_{\beta,\beta_{1},\beta_{2},\beta_{3},\beta_{4}}$ (see Eq. II). When changing the order of fermion operators, a sign which is determined by the parity of the indices will appear. For example, if we exchange the two adjacent fermion operators $a_{(i,j),r}^{\beta_{3}}$ and $a_{(i,j),u}^{\beta_{4}}$, there will be an extra phase, i.e., $T_{\beta,\beta_{1},\beta_{2},\beta_{4},\beta_{3}}=\tilde{p}(\beta_{3},\beta_{4})T_{\beta,\beta_{1},\beta_{2},\beta_{3},\beta_{4}}$, where

[TABLE]

Beside the fermion operators appeared in projector $\hat{P}$, we also need to specify the operators’ order in the EPR pairs, which is not given in the tensors explicitly. In this work, we introduce a Fermi arrow notation to specify the order of the EPR pairs. For example, as shown in Fig. 1(a), the arrow points from site A to site B , and the corresponding EPR state is $\hat{I}_{A\rightarrow B}|0\rangle=\sum^{D-1}_{k=0}a_{B}^{{\dagger}k}a_{A}^{{\dagger}k}|0\rangle$, whereas in Fig. 1(b), the arrow points from B to A, and the corresponding EPR state is $\hat{I}_{B\rightarrow A}|0\rangle=\sum^{D-1}_{k=0}a_{A}^{{\dagger}k}a_{B}^{{\dagger}k}|0\rangle$. The two states can be transformed to each other as follows,

[TABLE]

We may also assign Fermi arrows to the physical indices: the Fermi arrows point into the sites for the physical creation operators, and pointing out of the sites for the annihilation operators. This definition is equivalent to insert EPR pairs between the physical operators when contracting the physical indices e.g., $\langle 0|\sum_{\beta_{A},\beta_{B}}A^{\dagger}_{\beta_{A}}a_{A}^{\beta_{A}}B_{\beta_{B}}a_{B}^{{\dagger},\beta_{B}}|0\rangle=\langle 0|\sum_{\beta_{A}}A^{\dagger}_{\beta_{A}}a_{A}^{\beta_{A}}\sum_{\beta_{B}}B_{\beta_{B}}a_{B}^{\beta_{B}}\sum_{k}a_{B}^{{\dagger}k}a_{A}^{{\dagger}k}|0\rangle$. With this definition, we can treat the physical indices and the virtual indices in the same manner, and do not need to distinguish the real fermions and virtual fermions during operations. We can now uniquely define a fPEPS on a directed lattice, as shown for example in Fig. 2, on a 4$\times$4 lattice.

By defining some calculation rules associated with Fermi arrows, we are able to perform fPEPS calculations. Contraction is one of the most important operations in PEPS algorithms. When we contract the tensors on two sites in a fPEPS, we need to consider the Fermi arrow direction. We take the two situations in Fig. 1 as an example, which gives two different contraction formula,

[TABLE]

for Fig. 1(a) and,

[TABLE]

for Fig. 1(b). The anti-commutation relation of fermions has been used. The Fermi arrow helps to distinguish the two situations of contraction in the graphical notions of Fig. 1. Using the graphic representation may greatly simplify the notation.

More generally, giving two tensors $\mathbf{A}$ and $\mathbf{B}$, connected via multi virtual bonds (EPRs), where $\{i_{1},i_{2},\cdots,i_{k},i_{k+1},\cdots,i_{N}\}$ are the joint bonds to be contracted. Assume that bonds $\{i_{1},i_{2},\cdots,i_{k}\}$ have Fermi arrows pointing from $\mathbf{B}$ to $\mathbf{A}$ , and the rest bonds $\{i_{k+1},\cdots i_{N}\}$ have Fermi arrows pointing from $\mathbf{A}$ to $\mathbf{B}$. We first reshape $\mathbf{A}$ to $A_{I,i_{N},i_{N-1},\cdots,i_{1}}$ and reshape $\mathbf{B}$ to $B_{i_{1},i_{2},\cdots,i_{N},J}$, where $\{I\}$, $\{J\}$ are the bonds that are not to be contracted in $\mathbf{A}$ and $\mathbf{B}$ respectively. For the convenience of notation, we assume the signs due to the change of bond order in the tensors according to Eq. 7 have been absorbed into the tensors, then the resulting tensor is,

[TABLE]

Other often used operations associated with Fermi arrows are given in the Appendix. We implement these operation rules for fPEPS in our recently developed TNSpack Dong et al. (2018), in which the anti-commutation properties of the fermion operators are automatically taken account of by these rules. Therefore, one can implement fPEPS algorithm very similar to that of standard PEPS without worrying too much about the details of the anti-commutation between the fermion operators.

The Fermi arrows define the fermionic order for the fPEPS. In some previous methods, Corboz et al. (2010) the EPS pairs are not explicitly used. We note that in this work, the EPR pairs are only used in the derivation of the operation rules associated with Fermi arrows. Once we have these rules, one may ignore the underlying EPR pairs, and use only Fermi arrows for all operations. In Ref. Barthel et al., 2009, the authors proposed a general fermionization procedure using so called fermionic operator circuits (FOCs), in a bra and ket notation, instead of EPR pairs. Our Fermi arrows are similar to the contraction arcs defined in Ref. Barthel et al., 2009 albeit the starting point and detailed implementations of the two methods are different.

III Stochastic gradient optimization of fPEPS

In order to find the ground state of a given Hamiltonian using fPEPS, different methods have been introduced. The leading method is the imaginary time evolution (ITE) method.Jiang et al. (2008) However, due to the high computation complexity to obtain the exact environment during the time evolution, some kinds of approximations are necessary. The simple update methodJiang et al. (2008) has been widely used, which however may have large errors because the environment is over simplified. Several methods have been developed to treat the environment more rigorously, such as the full update methodLubasch et al. (2014a, b), and the gradient methodVanderstraeten et al. (2016); Liu et al. (2017), which may significantly improve the results, but the scaling to $D$ of these methods is rather high.

We recently developed stochastic gradient optimization method for PEPS, combined with Monte Carlo sampling techniques. Liu et al. (2017); He et al. (2018) This method gives remarkable accuracy of the results which may be even better than the results of full update method at given $D$.Liu et al. (2017) The method has two advantages. First, the environments of tensors are treated rigorously, and therefore, the results are more accurate than SU and even FU methodsLiu et al. (2017). Secondly, the Monte Carlo sampling technique may greatly reduce the scaling of the method to the virtual bond dimension $D$ from $D^{10}$ to $D^{6}$ for OBC,Sandvik and Vidal (2007); Schuch et al. (2008); Sfondrini et al. (2010); Bañuls et al. (2017) which is even more crucial for fPEPS, where larger $D$ is often needed to converge the results. In this work, we extended this method to fPEPS.

The fPEPS wave functions of a many-particle system in Eq. 4 can be rewritten as,

[TABLE]

where ${\bf i}$=$(i,j)$ is the site index of the lattice, and $\mathcal{C}$ means to contract all the entangled virtual fermions according to the rules defined in Sec. II. $W({\bf n})$ is the coefficient of the physical state $|{\bf n}\rangle$= $|{n_{1},n_{2},\dots,n_{N}}\rangle$ in the particle number representation. The energy of the system can be written as,

[TABLE]

where,

[TABLE]

The total energy of the system for a given fPEPS can be evaluated via Monte Carlo sampling over the physical configurations space. Sandvik and Vidal (2007); Schuch et al. (2008); Sfondrini et al. (2010); Bañuls et al. (2017); Liu et al. (2017)

To optimize the energy function, we need the derivatives of the energy with respect to the tensor elements,

[TABLE]

where $\langle\cdots\rangle$ denotes the MC average. $\Delta[{\bf i}]_{n,\beta_{1},\beta_{2},\beta_{3},\beta_{4}}$ is defined as

[TABLE]

and the derivative of $W({\bf n})$ is

[TABLE]

The derivatives can be also evaluated by the MC samplings.Liu et al. (2017)

Once we have the energy and its gradients, we can optimize the system energy using stochastic gradient methodSandvik and Vidal (2007); Liu et al. (2017), which has been successfully applied to the standard PEPS method.

The overall algorithm for fPEPS is similar to that of PEPS. We need to contract the fPEPS tensors at given particle configuration to obtain $W({\bf n})$ and the gradients. However, contracting a fPEPS is much more complicate than contacting a standard PEPS, because of the anti-commutation relation of fermions. One must be very careful about the contraction order and underlying fermions’ order in EPR pairs. We show here that with the help of Fermi arrows and the operation rules associated with them, the contraction can be done easily as in the standard PEPS algorithms.

To obtain $W({\bf n})$, we need to contact a single layer of fPEPS with fixed particle configuration $|{\bf n}\rangle$. We adopt the boundary-MPO method, Lubasch et al. (2014a, b)where we need to find a fermionic matrix product operator (fMPO) denoted as $|X\rangle$ with bond dimension $D_{c}$ [see Fig.3(b)] to approximate the two rows of fPEPS with bond dimension $D$, denoted as $|\Psi\rangle$ [see Fig.3(a)]. To find such $|X\rangle$, we minimize

[TABLE]

which lead to the linear equation for each tensor $T^{i,j}$ on site $(i,j)$,

[TABLE]

where $|X^{i,j}\rangle$ is obtained by taking the tensor $T^{i,j}$ out of $|X\rangle$, as graphically shown in Fig. 4.

To solve the equation, we first contract the tensors on the left side of Fig. 4. We change the arrow directions from Fig. 3(b) to Fig. 5(a), i.e., all arrows are pointing into site $(i,j)$. The rule of changing the directions of Fermi arrows are given in the Appendix. As will be seen in the following text, the change of Fermi arrow directions is to take the advantages of the canonic form of fMPO.Perez-Garcia et al. (2007)

We next do QR decomposition to the tensor on the first site of $|X\rangle$, resulting in two tensors, $Q_{1}$ and $R$ as shown in Fig. 5(b). The rules for QR (and other decompositions) in the presence of Fermi arrows are also given in the Appendix. We then contract the $R$ tensor with the second tensor on the right site, and perform QR decomposition on the second site again to obtain the $Q_{2}$ tensor. We repeat this process until reach the tensor $T^{i,j}$. We contract the last $R$ tensor with the $T^{i,j}$, resulting in a new tenor $T^{i,j}_{R}$. Similarly, we perform the LQ decomposition on the right side of $|X\rangle$, starting from the last site to the site ($i$, $j$), and contract the last $L$ tensor with $T^{i,j}_{R}$ to get $T^{i,j}_{RL}$. During the LQ (QR) decompositions, new Fermi arrows have been inserted between L (Q) tensors and Q (R) tensor. After these processes, we obtain $|X\rangle$ in Fig. 5(c).

We perform the same operations for $|X^{i,j}\rangle$. After these operations, the left side of Fig. 4 become that of Fig. 6. By using the orthogonality of $\hat{Q}$ and $\hat{Q}^{\dagger}$, i.e., $\hat{Q}\hat{Q}^{\dagger}$=${\bf I}$, which is discussed in the Appendix for the fPEPS with Fermi arrows, we obtain the right side of Fig. 6. The original equation Fig. 4 become of Fig. 7, which can be solved iteratively as in standard boundary MPO method,Lubasch et al. (2014a, b) which usually converges in a few sweeps.

The contraction in Eq. 17 can be calculated in the same procedures. Once we have $W({\bf n})$, and $\Delta_{n_{i}\beta_{1}\beta_{2}\beta_{3}\beta_{4}}$, the energies and their gradients can be easily calculated.

In our calculations, we first perform ITE with simple update method to obtain a approximate ground state,Liu et al. (2017) which usually have energy errors around $10^{-2}$. We further optimize the fPEPS via gradient decent method to obtain the highly accurate ground state.

IV benchmark results

We benchmark our method for two typical fermion models on finite size square lattices, including the interacting spin-less fermions model and the t-J model. We demonstrate that our method can give highly accurate results compared with the exact results.

IV.1 The spin-less fermions model

The interacting spin-less fermions model reads,

[TABLE]

where $c_{i}^{\dagger}$ and $c_{j}$ are the creation and the annihilation operators, and $\langle i,j\rangle$ denotes the nearest neighbor pairs. The chemical potential $\mu$ is set to zero here. We set hopping parameter $t$=1, and the interaction strength $V\geq 0$. In general this model is not exactly solvable, and has been numerically investigated in Ref. de Woul and Langmann, 2010 by the mean field theory and in Ref. Corboz et al., 2010 by the iPEPS method. Both methods give similar phase diagrams. For the parameters we used, the ground state is in a uniform metallic phase when $V$ is small and moves towards the phase boundary between the uniform phase and the phase separation with the increasing of $V$. Therefore the ground state of the model is expected to have entanglement beyond area law. Wolf (2006)

We firstly calculate this model on a 4$\times$4 square lattice so we can compare the fPEPS results with those obtained from the exact diagonlization method. In the calculations, we take $D=10$, and $D_{c}=30\sim 40$. The convergence of these parameters will be discussed in details in Sec.V. The results are presented in Table 1 for various $V$. As seen from Table I, the SU method may give the results with errors around 5$\times$10*-3* when $V$ is small, but the errors increase for larger $V$. When $V$=2, the error of SU is about 10*-2*. The GO may significantly improve the ground state energies. By using the given $D$ and $D_{c}$, we are able to obtain an impressive highly accurate ground state with relative error near $\sim 10^{-5}$.

We now consider a special case of $V$=0, where the model reduces to the free fermion model. Although in this case, the model is exactly solvable, it is a challenging model for the fPEPS method because the free fermions have strong entanglement in real space $S\sim L^{d-1}{\rm log}L$ that violates the area law Wolf (2006). Especially at $\mu$=0, the Fermi surface is very large, making the problem more difficult. One may expect that to obtain the high accuracy results it requires very large $D$ and $D_{c}$. Furthermore, the required parameters $D$ and $D_{c}$ will generally increase rapidly with the increasing of the size of the system to keep the given accuracy. In Table 2, we list the ground states energies of the free fermion model on the square lattice with different sizes obtained from the SU and GO methods, compared with the exact results $E_{ex}$. We see that the relative errors of the SU are usually about $10^{-2}$ for $D$=8, but sometimes the SU method may have numerical instability in some small systems. The performance of GO is much better, and we always get stable results. Even with a small $D$=6, the relative errors are about $10^{-3}\sim 10^{-4}$, and reduce to $\sim$ $10^{-5}$ when $D$=8 is used. On the other hand, the violation of the area law is also showed in this table, that the accuracy gets lower in larger systems for a given $D$.

From the above tests, we find that the SU method sometimes may give rather accurate results ($\sim$ 10*-2* - 10*-3*), but the situation may change from case to case. On the other hand the GO always gives reliable and highly accurate results ($\sim$ 10*-5*).

IV.2 t-J model

In this section, we benchmark our method on the t-J model,

[TABLE]

where $\sigma=\uparrow,\downarrow$ is the spin index and $\vec{S}_{i}$ is the spin $1/2$ operator on site $i$. $n_{i}=\sum_{\sigma}c_{i,\sigma}^{\dagger}c_{i,\sigma}$ is the number of electrons on site $i$. In t-J model, the electron double occupancy is forbidden. The t-J model is one of the key models to understand many important physical phenomenaZhang and Rice (1988), such as high $T_{c}$ superconductivity Lee et al. (2006). Here, we calculate the model with $J/t=0.4$ and hole filling of $\bar{n}_{h}=$0.125. The U(1) symmetry is adopted to enforce the particle number conservation. But true physics of the system at this point, whether the ground state is a stripe state Hellberg and Manousakis (1999); White and Scalapino (1998); Corboz et al. (2011, 2014) or an uniform phaseSherman and Schreiber (2003); Vineet Mallik et al. (2018), is still under debate. Without doubt, the energy is one of the critical criterions to determine the ground state of the system. We calculate the ground energies of different system sizes, using $D=12$, $D_{c}=36\sim 50$, and the results are shown in Table 3 for both SU and GO methods. Again we see GO method greatly improves the energies obtained from SU method. By extrapolating the energy to thermodynamic limit, we obtain that the ground state energy -0.6701, which is lower than the value -0.6693 Corboz et al. (2014) from state of art DMRG calculations White and Scalapino (1998) with $\chi\rightarrow\infty$, and -0.6619 obtained from variation quantum Monte Carlo plus $p$-step Lanczos methods. Hu et al. (2012) More results of t-J model Don will be published in a separate paper.

V convergence of fPEPS

Fermion systems may have large entanglement that beyond the area law Wolf (2006) and therefore it may need large $D$ to represent the many-particle state. One may also expect that the $D$ and $D_{c}$ will increase with the size of the system. The speed of the increasing of $D$ and $D_{c}$ along with the size of the system indicates the efficiency of the simulation methods. It is important to understand how the fPEPS calculations converge respect to $D$ and $D_{c}$. For finite systems, we can explicitly exam what $D$ and $D_{c}$ are needed to converge the results as the size of the systems grow up. In this section, we will discuss the convergence of the parameters $D$ and $D_{c}$ respectively. We show that the behavior depends strongly on the models.

We first investigate the convergency of the ground state energies to $D_{c}$ in a given system with fixed parameter $D$. We calculate the error of energy defined as,

[TABLE]

where $E(D_{c})$ is the energy with a giving truncation parameter $D_{c}$ and $E(D_{c}^{\rm max})$ is the converged energy where the maximal $D_{c}$ is used.

Figure 8(a) depicts the results of the spin-less fermion model with $V=0$ (free fermion) and $V=2$; and different system sizes, the $4\times 4$ and $10\times 10$ lattices. We fix the bond dimension at $D=8$. We use $D_{c}^{\rm max}=48$ for the $4\times 4$ system and $D_{c}^{\rm max}=64$ for the $10\times 10$ system. We first note that $\Delta E$s approach 0 in a non-trivial way, which are not always from above (i.e., $\Delta E>$0). This means that the ground state energy is not variational to $D_{c}$, and therefore one must be very careful to extrapolate $D_{c}$ to infinite. The convergency of energy is model dependent. As shown in the figure, $D_{c}$ converge much faster for $V=2$ (correlated electrons) than for $V$=0 (free electrons). In both cases, the convergency of energy strongly depend on the size of the systems. In the cases of small sizes $L$=4, the energies converge rather fast with $D_{c}$. However, for the $10\times 10$ system, $\Delta E$ converge much slower as functions of $D_{c}$. For $V=2$, the energy is well converged at $D_{c}=26$ (about 3$D$), whereas the energy of free fermions is not well converged even at $D_{c}=48$.

We investigate the convergence of the t-J model at hole doping $\bar{n}_{h}$=0.125, and the results are shown in Fig. 8(b). In the calculations, $D=12$ is used, and the result of $D_{c}^{\rm max}$=50 is used as a reference. Interestingly, we find the ground state energies converge rather fast with $D_{c}$. The errors reduce to 3$\times 10^{-4}$ for $D_{c}$=2$D$, and the errors reduce to 1$\times 10^{-5}$ for $D_{c}$=3$D$, More importantly, unlike the interacting fermion model, $D_{c}$ is only slightly dependent on the size of the system.

The energy errors in the calculations are induced by the contraction errors due to bond dimension truncation. We further test the relationship between the convergent truncation dimension $D_{c}$ and the size of the system, i.e. we exam the minimal $D_{c}$ needed to ensure the relative contracting error $\delta<10^{-6}$ (see Eq. 18) along with the increasing of the system size. The bond dimension used here is fixed to a relatively small one $D=6$. We compare the truncation errors for the spin-less interacting electron model at $V$=0 and $V$=2. For the t-J model, we compare two situations, the hole doping $\bar{n}_{h}$=0.125, and the $\bar{n}_{h}$=0, and the later one reduces to the Heisenberg model. The results are shown in Fig. 9(a). We find that the required $D_{c}$ is almost independent of the size of the system for the Heisenberg model, and for the t-J model with hole filling $\bar{n}_{h}$=0.125. However, the required $D_{c}$ increases rapidly with the size of the system for the interacting electron model, especially for the free electrons. At $L$=12, $D_{c}$=80 is required to ensure the desired contraction accuracy for the free electron model and $D_{c}$=40 for the $V$=2 model.

We also exam the relationship between $D_{c}$ and the bond dimension $D$. In this test, we fixed the size of the system to $L$=10. The results are shown in Fig. 9(b). We see that the required $D_{c}$ increase roughly linearly with $D$ for these models. For the Heisenberg model (and even J1-J2 model) Liu et al. (2018) and the t-J model with hole filling $\bar{n}_{h}=$0.125, D_{c}$$\sim 2$D$ - 3$D$ is enough to ensure the accuracy of contraction, whereas for the interacting electron model, D_{c}$$\sim 9$D$ - 15$D$ are required to ensure the desired contraction accuracy, which becomes the major difficult to simulate these models. We note that in the standard contraction method, the bond truncation dimension $D_{c_{2}}$ for a double layer tensor network should scale as $D^{2}$, Lubasch et al. (2014a, b) making the simulation of fermions with large $D$ even more difficult.

With the convergent $D_{c}$ for each $D$, we can analyze the convergence of the energy against $D$ for a given system size. The energy of a model in the thermodynamic limit can be further extracted by finite size scaling method. The convergence of the energy against the parameter $D$ are shown in Fig. 10, where $\Delta E$ is defined as the energy differences compared to those of maxima $D$, which are $D^{\rm max}$=8 for the free fermion model and $D^{\rm max}$=10 for the interacting fermion model with $V$=2. A $D^{\rm max}$=12 is used for the t-J model with $\bar{n}_{h}$=0.125. Surprisingly, for the free fermion model ($V$=0), we have $\Delta E\approx 4\times 10^{-4}$ for $D$=7 on the 10$\times$10 lattice, as shown in Fig. 10(a), where one may expect a much larger error. For the interacting fermion with $V=2$, which is shown in Fig. 10(b), the energy is also converged to $\Delta E\approx 1\times 10^{-4}$ at $D$=8. On the other hand, for the t-J model, the energies converge rather slowly with $D$. For the $12\times 12$ system, the energy errors reduces to about $3\times 10^{-4}$ at $D=11$. The non-trivial dependence of $D$ and $D_{c}$ for different models may pose some interesting questions to understand the structure of fPEPS. We leave these problems for future studies.

VI Summary

In this work, we extend the stochastic gradient optimization method combined with Monte Carlo sampling techniques to optimize the fPEPS wave functions for fermion systems. The Monte Carlo sampling techniques may greatly reduce the scaling of the calculation, and therefore allow using larger bond dimensions ($D$) and bond truncation dimensions ($D_{c}$) in the calculations, which is important for the faithful simulations of fermion systems.

We benchmark the method on the interacting spinless fermion models, and the t-J models. The numerical calculations show that the gradient optimization may greatly enhance the accuracy of the results over the simple update method. We further investigate the converge of fPEPS calculation with respect to $D$ and $D_{c}$ for the models. The free fermion model is most challenging to simulate with fPEPS, because the $D_{c}$ increase very rapidly with $D$ and the size of the system. For t-J models, we find that large $D$s are needed to converge the results. Our method therefore offer a powerful tool to simulate fermion systems because it has much lower scaling in both computational time and memory than direct contraction methods.

VII Acknowledge

This work was funded by the National Key Research and Development Program of China (Grant No. 2016YFB0201202), the Chinese National Science Foundation (Grants No. 11774327, No. 11874343, No. 11474267), and the Strategic Priority Research Program (B) of the Chinese Academy of Sciences (Grant No. XDB01030200). China Postdoctoral Science Foundation funded project (Grant No. 2018M632529). The numerical calculations have been done on the USTC HPC facilities.

Appendix A Rules for Fermi arrows

In this Appendix, we give the rules of operations associated with the Fermi arrows in fPEPS. These rules are straightforward to prove.

A.0.1 Reversing Fermi arrows and the Hermitian conjugate

Sometimes, we need to reverse the direction of a Fermi arrow. The rule of reversing Fermi arrows is giving as follows. Suppose,

[TABLE]

are two projectors in a fPEPS that are connected by a Fermi arrow pointing from $\hat{A}$ to $\hat{B}$, as shown on the left side of Fig. 11. We may reverse the Fermi arrow, pointing from $\hat{B}$ to $\hat{A}$, and resulting in two possible (but equivalent) forms that are given on the right side of Fig. 11. It is easy to prove that,

[TABLE]

When we calculate the expectation value of a physical quantity, $\langle\Psi_{\rm fPEPS}|\hat{O}|\Psi_{\rm fPEPS}\rangle$, we need take the Hermitian conjugate of ket state $|\Psi_{\rm fPEPS}\rangle$ to get the bra state $\langle\Psi_{\rm fPEPS}|$. When taking the Hermitian conjugate of the projectors in a fPEPS, we need to (i) reserve the orders of the indices of the tensor associated with the projectors, e.g., change tensor $A_{l,d,r,u}$ to $A^{{\dagger}}_{u,r,d,l}$, as shown in Fig. 12; and (ii) reverse all the Fermi arrows associated with the projectors. Note that here the reversion of the Ferim arrows is required by the Hermitian conjugate, and no change is needed for the tensors during the process.

A.0.2 Matrix decompositions and contractions

The operations such as tensor decompositions also have close relation to the Fermi arrows. For example, in the standard PEPS, when we do SVD to a matrix $C$, we have $C$=$USV$. However, in fPEPS, two Fermi arrows should be inserted to the inner bonds after the decomposition, i.e., the Fermi arrow pointing from $\hat{U}$ to $\hat{S}$, and the one pointing from $\hat{S}$ to $\hat{V}$ as follows, and schematically shown in Fig. 13(a),

[TABLE]

where

[TABLE]

Other matrix decompositions such as LQ/QR decompositions follow the similar rules, i.e., one need to insert Fermi arrows (i.e., directed EPR pairs) between the decomposed matrices, as shown in Fig. 13(b),(c).

In standard PEPS, we often use so called canonical form of MPS in the MPO algorithmLubasch et al. (2014a, b) to contract the PEPS, taking the advantage of the orthogonality of the $Q_{i}$ tensors obtained from LQ/QR decompositions (or the $U$ and $V$ matrices from SVD decompostions), Schollwöck (2011) i.e., $Q_{i}Q_{i}^{{\dagger}}$=${\bf I}$, where ${\bf I}$ is a unit matrix. However, this relation cann’t be directly used in the fPEPS, where we need to take the Fermi arrows into consideration during the contractions. It is easily prove that only when the Fermi arrows have “consistent directions”, i.e., all Fermi arrows point from $\hat{Q}$ to $\hat{Q}^{{\dagger}}$, or from $\hat{Q}^{{\dagger}}$ to $\hat{Q}$, we can use the orthogonality condition for $\hat{Q}$ matrix. The results after contraction are Fermi arrows pointing to the right or to the left, as schematically shown in Fig. 14. If the Fermi arrows are not “consistent”, we need to rearrange the directions of the Fermi arrows first to make them “consistent”, before we can use the orthogonality condition. This is done in Sec.III, when we contract two rows of fPEPS via a MPO scheme.