Variational Monte Carlo method for fermionic models combined with tensor networks and applications to the hole-doped two-dimensional Hubbard model

TL;DR

This paper introduces a combined variational Monte Carlo and tensor network approach for fermionic models, demonstrating improved accuracy and revealing complex ground states in the doped Hubbard model.

Contribution

It develops a novel wave function combining pair product states, tensor networks, and symmetry projections, enhancing fermionic model simulations.

Findings

Significant accuracy improvement over individual methods.

Observation of stripe order coexisting with weak d-wave superconductivity.

Identification of metastable superconducting states at low doping.

Abstract

The conventional tensor-network states employ real-space product states as reference wave functions. Here, we propose a many-variable variational Monte Carlo (mVMC) method combined with tensor networks by taking advantages of both to study fermionic models. The variational wave function is composed of a pair product wave function operated by real space correlation factors and tensor networks. Moreover, we can apply quantum number projections, such as spin, momentum and lattice symmetry projections, to recover the symmetry of the wave function to further improve the accuracy. We benchmark our method for one- and two-dimensional Hubbard models, which show significant improvement over the results obtained individually either by mVMC or by tensor network. We have applied the present method to hole doped Hubbard model on the square lattice, which indicates the stripe charge/spin order coexisting with a weak $d$-wave superconducting order in the ground state for the doping concentration less than 0.3, where the stripe oscillation period gets longer with increasing hole concentration. The charge homogeneous and highly superconducting state also exists as a metastable excited state for the doping concentration less than 0.25.

Tables

Table 1. Table 1: Ground state energies per site for 2D Hubbard model with U = 4.0 𝑈 4.0 U=4.0 at half filling for 8 × 8 8 8 8\times 8 and 16 × 16 16 16 16\times 16 square lattices. The energies obtained from the first Lanczos step applied to the combined mVMC and FTTN, the variance extrapolation, and the QMC MPQin2016 , which are plotted in Fig.  5 , are listed. The parentheses denote the error bars in the last digit.

	VMC+FTTN+Lanczos	Var. Ext.	QMCMPQin2016
$8\times 8$	$-0.8636\left(2\right)$	$-0.8641\left(2\right)$	$-0.8642\left(2\right)$
$16\times 16$	$-0.8581\left(2\right)$	$-0.860\left(1\right)$	$-0.8605\left(2\right)$

Table 2. Table 2: Ground state energies per site of 16 × 16 16 16 16\times 16 2D square lattice Hubbard model at U = 10 𝑈 10 U=10 with the periodic-antiperiodic boundary condition for various doping concentrations ( δ ) 𝛿 (\delta) . The energies obtained from the first Lanczos step applied to the combined mVMC and FTTN, and the variance extrapolation with stripe and homogeneous initial states which are plotted in Fig.   6 , are listed. The parentheses denote the error bars in the last digit.

$\delta$	VMC+FTTN+Lanczos	Var. Ext.	homo. VMC+FTTN+Lanczos	homo. Var. Ext.
$0.00$	$-0.43018\left(6\right)$	$-0.4361\left(1\right)$
$0.05$	$-0.51782\left(8\right)$	$-0.525\left(5\right)$
$0.06$	$-0.54963\left(7\right)$	$-0.556\left(2\right)$	$-0.54497\left(6\right)$	$-0.5537\left(6\right)$
$0.08$	$-0.5822\left(2\right)$	$-0.591\left(3\right)$
$0.11$	$-0.6442\left(4\right)$	$-0.6538\left(9\right)$	$-0.6309\left(1\right)$	$-0.644\left(1\right)$
$0.14$	$-0.6988\left(3\right)$	$-0.711\left(1\right)$
$0.19$	$-0.7822\left(2\right)$	$-0.796\left(3\right)$
$0.22$	$-0.8343\left(9\right)$	$-0.8466\left(6\right)$	$-0.8318\left(2\right)$	$-0.843\left(1\right)$
$0.25$	$-0.8829\left(4\right)$	$-0.897\left(1\right)$
$0.28$	$-0.9235\left(2\right)$	$-0.938\left(2\right)$

Equations

$$H=-t\sum _{⟨i,j⟩,\sigma }\left(c_{i\sigma }^{†}{c}_{j\sigma }+h.c.\right)+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow },$$

$$|{\varphi }_{\mathrm{ref}}⟩=\mathcal{P}{\mathcal{L}}^{\mathrm{Space}}{\mathcal{L}}^S{\mathcal{L}}^K|{\varphi }_{\mathrm{pair}}⟩,$$

$$\mathcal{P}={\mathcal{P}}_{\mathrm{G}}{\mathcal{P}}_{\mathrm{J}}{\mathcal{P}}_{\mathrm{d}-\mathrm{h}}^{\mathrm{ex}},$$

$${\mathcal{P}}_{\mathrm{G}}=\exp \left(-\sum _i{g}_i{n}_{i\uparrow }{n}_{i\downarrow }\right)$$

$${\mathcal{P}}_{\mathrm{J}}=\exp \left[\frac{1}{2}\sum _{i\ne j}{v}_{ij}\left(n_{i\uparrow }+n_{i\downarrow }\right)\left(n_{j\uparrow }+n_{j\downarrow }\right)\right],$$

$${\mathcal{P}}_{\mathrm{d}-\mathrm{h}}^{\mathrm{ex}}=\exp \left(-\sum _{m=0}^4\sum _l{\alpha }_{\left(m\right)}^{\left(l\right)}\sum _i{\xi }_{i\left(m\right)}^{\left(l\right)}\right),$$

$$|{\varphi }_{\mathrm{pair}}⟩={\left(\sum _{i,j=1}^{N_s}\sum _{\sigma {\sigma }^{\prime }}{f}_{i\sigma j{\sigma }^{\prime }}{c}_{i\sigma }^{†}{c}_{j{\sigma }^{\prime }}^{†}\right)}^{N_e/2}|0⟩,$$

$$|\Psi ⟩$$

$$\mathcal{M}|q_1,q_2,\dots ,q_{N_s}⟩$$

$$|{\psi }_n⟩=\left(1+\sum _{n=1}^N{\alpha }_n{H}^n\right)|\Psi ⟩,$$

$$⟨H⟩=\frac{{\sum }_x⟨\Psi |H|x⟩⟨x|\Psi ⟩}{⟨\Psi |\Psi ⟩}=\sum _x\rho \left(x\right)E\left(x\right)$$

$$E\left(x\right)=\sum _{x^{\prime }}\frac{⟨\Psi |x^{\prime }⟩}{⟨\Psi |x⟩}⟨x^{\prime }\left|H\right|x⟩,\rho \left(x\right)=\frac{⟨\Psi |x⟩⟨x|\Psi ⟩}{⟨\Psi |\Psi ⟩},$$

$$e^{-\tau H}\approx 1-\tau H.$$

$$\left(1-\tau H\right)|\overline{\Psi }⟩\approx |\overline{\Psi }⟩+\sum _{k=1}^{n_p}{\gamma }_k|{\Psi }^k⟩,$$

$$|{\Psi }^k⟩$$

$$f\left({\gamma }_k\right)$$

$$\frac{\partial f\left({\gamma }_k\right)}{\partial {\gamma }_k}=0.$$

$$\sum _{l=1}^{n_p}⟨{\Psi }^k|{\Psi }^l⟩{\gamma }_l=-\tau ⟨\overline{\Psi }|H|{\Psi }^k⟩.$$

$${\tilde{\alpha }}_k={\alpha }_k+{\gamma }_k.$$

$${\tau }_k=-\Delta t\left(i\right)\eta \left(i,k\right),$$

$$\Delta t\left(i\right)={\xi }^i,$$

$$P_d\left(\mathbf{r}\right)$$

$${\Delta }_d\left({\mathbf{r}}_i\right)=\frac{1}{\sqrt{2}}\sum _{\mathbf{r}}{f}_{d_{x^2-y^2}}\left(c_{{\mathbf{r}}_i\uparrow }{c}_{{\mathbf{r}}_i+\mathbf{r}\downarrow }-c_{{\mathbf{r}}_i\downarrow }{c}_{{\mathbf{r}}_i+\mathbf{r}\uparrow }\right),$$

$$f_{d_{x^2-y^2}}\left(\mathbf{r}\right)={\delta }_{r_y,0}\left({\delta }_{r_x,1}+{\delta }_{r_x,-1}\right)-{\delta }_{r_x,0}\left({\delta }_{r_y,1}+{\delta }_{r_y,-1}\right),$$

$$S_s\left(\mathbf{k}\right)=\frac{1}{3N_s}\sum _{\mathbf{r},{\mathbf{r}}^{\prime }}⟨s_{\mathbf{r}}^z{s}_{\mathbf{r}+{\mathbf{r}}^{\prime }}^z⟩e^{i\mathbf{k}\cdot {\mathbf{r}}^{\prime }},$$

$$S_c\left(\mathbf{k}\right)=\frac{1}{N_s}\sum _{\mathbf{r},{\mathbf{r}}^{\prime }}⟨n_{\mathbf{r}}{n}_{\mathbf{r}+{\mathbf{r}}^{\prime }}⟩e^{i\mathbf{k}\cdot {\mathbf{r}}^{\prime }},$$

$$S^z\left(\mathbf{r}\right)=⟨n_{\mathbf{r}\uparrow }-n_{\mathbf{r}\downarrow }⟩,$$

$$n\left(\mathbf{r}\right)=⟨n_{\mathbf{r}\uparrow }+n_{\mathbf{r}\downarrow }⟩,$$

$$\mathcal{M}|q_1,q_2,\dots ,q_{N_s}⟩$$

Full text

I Introduction

Finding the ground state of strongly correlated electron systems is one of most challenging problems in condensed matter physics. Since exact solutions only exist in few systems, deeper understanding of ground state properties strongly relies on efficient and accurate numerical algorithms. For example, one can employ exact diagonalization (ED) to find the wave function accurately, but it is only applicable for very small size systems. The density matrix renormalization group (DMRG)(DMRG1992, ) is very efficient and accurate for one dimensional systems, but it becomes inefficient for two and higher dimensional systems. The quantum Monte Carlo methods(QMC2001, ) suffer from the sign problem in general fermionic and geometrically frustrated systems.

In the past years, the tensor network algorithms have been widely developed(Niggemann1997, ; PTP2001Nishino, ; PEPS2004, ; MERA, ; TRG2007, ; PRL2007Vidal, ; SimpleUpdate2008, ; TEFR, ; SRG2009, ; SingleLayer2011, ; HOTRG2012, ; PESS2014, ; TNR, ), which are shown to be promising numerical tools. One of the simplest tensor network state is the matrix product state (MPS), which is the variational wave function of the DMRG method(Ostlund1995, ; DMRGpbc2004, ). A natural generalization of the MPS to two dimensions is the projected entangled pair state (PEPS)(PEPS2004, ), which satisfies the area law of entanglement entropy(PEPS2004, ). Besides PEPS, various types of tensor network states have been introduced, such as multi-scale entanglement renormalization ansatz (MERA)(MERA, ), tree tensor network states(TTN2009, ) and projected entangled simplex states(PESS2014, ). These tensor network wave functions are usually expressed in a real space basis which may become inefficient to capture the large amount of entanglement in itinerant fermionic systems. For example, the free fermion model, which can be exactly represented as a product state in momentum space, is extremely difficult to accurately treat by infinite fermionic PEPS algorithm(fPEPS2010, ).

The variational Monte Carlo (VMC)(Ceperley1977, ) method can be applied to study relatively large system sizes, and there is no sign problem in studies of fermionic and frustrated systems. However, the result is subject to be biased depending on the form of variational wave functions. In the region where various competing phases have very closed energy, it is very challenging to determine the correct ground state.

Sorella has developed stochastic reconfiguration (SR)(SR2001, ) method to stably optimize large number of parameters, which makes it possible to extend the variational wave function to substantially alleviate the bias. With the combination of pair-product wavefunctions and the correlation factors such as Gutzwiller Gutzwiller1963 , and Jastrow Jastrow1955 factors, as well as the quantum number projections, thousands and ten thousands of the variational parameters have been optimized, which has enabled accurate estimates of the competing ground states in terms of the comparisons with available exact results Edegger2007 ; mVMC2008 ; Becca2011 ; Misawa2014PRB . The applications have achieved fruitful and reproducible comparisons with the experimental results, for instance for the iron-based superconductors, if the method is applied to the ab initio effective HamiltoniansMisawa2014NatComm . However, how to further systematically remove the bias in the variational wave functions and reach better accuracy is still left as a challenge.

The Monte Carlo sampling techniques have been proposed as a prescription to reduce the computational cost of tensor network contraction in the MPS(MPSVMC2007, ), PEPS(PEPSVMC2011, ) and MERA(MERAVMC2012, ). Another advantage of employing the Monte Carlo sampling into tensor network methods is that it is possible to choose various types of suitable reference basis beyond the real space basis. In the VMC study of correlator product states(CPS2009, ), a Pfaffian pairing wave function has been used(Neuscamman2012, ). In the study of one dimensional fermionic system, the free fermion Slater determinant has been used as the reference wave function of the MPS(MPPSVMC2012, ), which achieves higher accuracy than the conventional MPS method with the same bond dimension. Ref. TPPSVMC2015, has generalized this idea to two dimensional systems, which employs the PEPS with various kinds of reference wave functions, such as Jastrow-type, free fermion, $d$-wave BCS and spin density wave states, in order to choose a suitable reference wave function that captures the key features of the systems.

In this paper, we employ a combination of tensor network and reference wave function which consists of a Pfaffian pairing wave function and real space correlation factors, and the tensor networks can be regarded as diagonal correlation projectors, which act in the same way as that in Refs. MPPSVMC2012, and TPPSVMC2015, . In order to provide more flexible representation, we optimize all the variational parameters both in the part of the VMC and the tensor network simultaneously. Moreover, we can apply quantum number projections, such as those to restore the spin, momentum and lattice symmetries to further improve the accuracy. As a result, highly accurate ground state wave functions are obtained within a computationally tractable size of tensor bond dimension beyond the accuracy of each method if applied separately. Moreover, the accuracy can be continuously improved with the increase of the tensor bond dimension, thus providing a systematic way of removing the bias in the VMC. The combination with the VMC is particularly powerful if the tensor network method suffers from the entanglement entropy remaining beyond the area law as in the case of the itinerant fermion systems.

In the latter part of this paper, we show a fruitful application to a strongly correlated system: Hubbard model on the square lattice. Although the relevance of the Hubbard model for the mechanism of the high $T_{\rm c}$ superconductivity is an open issue, the ground states of the hole-doped two-dimensional Hubbard model has been extensively studied and debated for decades, because it is a simplest model of the cuprate superconductors and it captures several important experimental aspects. However, it was also suggested that strong competitions exist among different states including $d$-wave superconductivity and various charge inhomogeneous states such as stripe order and phase separation. Accurate determination of the phase diagram in the plane of the carrier density and the electron correlation strength is still an open question. The present method opens a possibility of studying the model at the highest accuracy among that ever studied. We show the numerical results of the $d$-wave superconducting correlations and various spin and charge correlations that indeed reveal the severe competition of the these two types of orders and clarifies how they are compromised in the best estimates of the ground states.

This paper is organized as follows. In Sec. II, we introduce our variational wave function, and describe how to optimize large number of parameters in variational Monte Carlo methods. In Sec. III, we present our benchmark results for the Hubbard model in one and two dimensions. In Sec.IV, we show results obtained by applying the present method to the hole doped Hubbard model to understand the interplay between the charge/stripe order and the $d$-wave superconductivity. Finally, we summarize the methods and results with the future scope in Sec. V.

II Numerical methods

II.1 Model

Although the method presented here is applicable to general fermionic systems on lattices, to make a presentation clearly understandable, we keep in mind the Hubbard model with the hopping amplitude $t$ between the nearest neighbor sites $\left\langle i,j\right\rangle$ and the onsite Coulomb repulsion $U$ defined by

[TABLE]

where $c_{i\sigma}^{\dagger}$ and $c_{j\sigma}$ are the creation and annihilation operators of electron with spin $\sigma$ at the $i$-th and $j$-th site, respectively, and $n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ is the number operator. Here, we take the energy unit $t=1$. We mainly consider the model on the $L\times L$ square lattice.

II.2 Variational wave function ansatz

The purpose of our work is to provide a flexible variational wave function which can be applied to efficiently capture the key features of the systems. We apply the wave function employed in the mVMC(mVMC2008, ) as the reference wave function of tensor network states, which is expressed as,

[TABLE]

where $\mathcal{P}$ is the product of real-space correlation factors

[TABLE]

in which $\mathcal{P}_{{\rm G}}$ is the Gutzwiller factor(Gutzwiller1963, ) that punishes (enhances) the double occupation of electrons on the same site defined as

[TABLE]

to take into account the local correlation effects, $\mathcal{P}_{J}$ is the Jastrow factor(Jastrow1955, ) which accounts for long-ranged density correlations through two-body operators defined as

[TABLE]

and $\mathcal{P}_{d-h}^{{\rm ex}}$ is the doublon-holon correlation factor(DH1982, ; DH1990, )

[TABLE]

where $\alpha_{\left(m\right)}^{\left(l\right)}$ are variational parameters, and $\xi_{i\left(m\right)}^{\left(l\right)}=1$ when a doublon (holon) exists at the site $i$ and $m$ holons (doublons) surround at the $l$-th nearest neighbor sites and otherwise, $\xi_{i\left(m\right)}^{\left(l\right)}=0$. The range of $l$ is taken short, normally up to the second nearest neighbors.

The two-body part $\phi_{{\rm pair}}$ is expressed in the real space representation

[TABLE]

where $f_{i\sigma j\sigma^{\prime}}$ denotes the variational parameters, and $N_{s}$, $N_{e}$ are the number of sites and electrons, respectively. Although the spins $\sigma$ and $\sigma^{\prime}$ are taken to be opposite in this paper to impose the constraint of the singlet pair, it can be easily generalized to arbitrary spins for better accuracy. The number of variational parameters are at most $N_{s}$, $N_{s}^{2}$, and $N_{s}^{2}$ for $g_{i},v_{ij}$ and $f_{ij}$, respectively. To save the computational cost, these numbers are sometimes reduced by imposing the sublattice translational invariance. In the case of the sublattice size $N_{{\rm sub}}$, the numbers are reduced to $N_{s}N_{{\rm sub}}$ for $v_{ij}$ and $f_{ij}$, while $g_{i}$ is taken site independent in the following study.

The two-body part $\phi_{{\rm pair}}$ may contain symmetry broken phases for better accuracy. On top of it, the space group, spin and momentum quantum number projection operators $\mathcal{L}^{{\rm Space}}$, $\mathcal{L}^{S}$ and $\mathcal{L}^{K}$ can be employed to recover the space group, spin SU(2) symmetry and lattice translational symmetry of the wave function to further improve the accuracy, because these symmetries are preserved for finite-size systems in most of the Hamiltonians we study. In this paper, we focus on the case of the singlet $S=0$, and the total momentum zero $K=0$ because they are satisfied in most of the ground state of models. We impose the space group symmetry later. We also note that by imposing these projections, one can also study the lowest excited state of the specified quantum number.

Combining with the tensor network algorithm, the variational wave function we are going to optimize is

[TABLE]

where $\mathcal{M}$ is the tensor network, and the physical index $q_{m}=1,2,\ldots,d$ represents $d$ local states at site $m$. Since the tensor network may break the lattice translational and rotational symmetry explicitly, the momentum projection $\mathcal{L}^{K=0}$ and space-group symmetry projection (for instance, the $\mathcal{L}^{C_{4}}$ to restore the $C_{4}$ rotational symmetry in the case of the square lattice) improve the wave function by recovering the symmetries. We apply the same quantum number projections to the Pfaffian pair state, because this preconditioning further improve the state.

The remaining task is to choose the appropriate tensor network $\mathcal{M}$ in Eq. (8). In one dimension, it is natural to choose $\mathcal{M}$ as an MPS

[TABLE]

because it holds the lattice translational symmetry and the computational cost is as low as $O\left(D^{3}\right)$ for the periodic boundary MPS, where $A^{m}[q_{m}]$ are $D\times D$ matrices. The operation Tr is to trace out all the matrix indices.

In two dimensional systems, various types of tensor network states may be employed. Among them we select and employ to meet several requirements. The first requirement is that the network structure should keep the lattice symmetry as much as possible. Therefore, the MPS will not be considered. The second requirement is that the tensor network exactly reproduces the ground state in the limit of infinite bond dimensions with nonzero reference basis, so that the accuracy is improved systematically. Therefore, the string bond state(SBS2008, ) is ruled out, because it does not guarantee the covering of the whole Hilbert space even when the infinite bond dimensions are taken. The third is that we prefer a tensor network in which the contraction can be done without truncation, so that the variational principle will not be broken during the optimization. Therefore, the PEPS is not employed, since the computational cost grows exponentially if the truncation is excluded.

In order to satisfy the above requirements, we build up the tensor network based on the idea of tree tensor network (TTN) as described in Appendix A. The TTN is a two dimensional tensor network which can be contracted exactly in polynomial time. If we employ the Monte Carlo sampling on the real space configuration, the physical indices of the leaf tensors are fixed, so every leaf tensor becomes a vector. Therefore, we can start from the contraction of the vector at the leaf tensor and then continue the contraction of rank 3 tensors at the higher hierarchical levels, of which the computational cost scales as $O\left(N_{s}D^{3}\right)$, where $D$ is the dimension of the virtual indices.

In the TTN for the $L\times L$ lattice, the number of bonds of tensors connecting any 2 sites is at most $O\left(\log\left(L\right)\right)$, which has a potentiality to capture the long range entanglement efficiently. However, the number of bonds of tensors connecting 2 neighboring sites can be as large as $O\left(\log\left(L\right)\right)$, if the neighboring sites belong to a different large block, which poses a limitation. For instance some of nearby sites belonging to different blocks are sparsely connected via internal-node tensors at a high hierarchical level.

We modify the standard network structure of the TTN to improve the efficiency while keeping the computational cost in the same order. The aim is to keep the neighboring sites as neighbors also in the tensor network. For this purpose, we propose a fat tree tensor network (FTTN), which contains redundant physical indices in the leaf node tensors. Consider an $L\times L$ square lattice system in which the local Hilbert space dimension of each site is $d$. As shown in Fig. 1 for $4\times 4$ lattice, the FTTN is composed of a set of tensors $t_{j,i}$, where $j=1,2,\cdots,R+1\left(R=\log_{2}N_{s}\right)$, and $i=1,2,\ldots,N_{s}/2^{j-1}$. The FTTN is connected as a binary tree structure, which can be expressed as

[TABLE]

where $t_{1,i}$ is the leaf tensor (meaning the end tensor of the tree structure) containing $4$ physical indices $q_{i},q_{a_{i}},q_{b_{i}},q_{c_{i}}$ of four sites on the plaquette with site $i$ on the top left vertex and $1$ virtual index, and $t_{2,i}$, $t_{3,i}$, $\ldots$, $t_{R+1,1}$ are internal-node tensors which only contain virtual indices. Every $2$ neighboring leaf tensors share 2 common physical indices. In the TTN, the number of leaf tensors connecting 2 neighboring sites is always 1, while it is 2 in the FTTN. The FTTN can reproduce the area law, because the TTN can be regarded as a subset of the FTTN with every leaf tensor contains only one physical index, and the TTN can reproduce the area lawTTN2009 . In the Monte Carlo sampling, the physical degrees of freedom are fixed in every sample, so the tensor network contraction in the FTTN will have the same computational complexity as in the TTN. One can also increase the number of sites in every leaf tensor to make further improvement.

It should be pointed out that regardless of the TTN or FTTN, the lattice translational symmetry and rotational symmetry are broken. (Nevertheless, that breaking is weaker than the MPS). To restore the symmetry, we operate the quantum number projection by the summation over the spatially translated and rotated tensor networks.

Besides tensor networks, one can also include backflow correlationsTocchio2008 ; Tocchio2011 ; Ido to further improve the correlation effect of the variational wave function. The backflow correlations can be implemented in the pair-product wave functionIdo , and in this case the additional computational cost of calculating kinetic energy arises because the pairing amplitude with backflow correlations is dependent on a real space configuration of electrons, which scales as $O\left(\gamma N_{s}^{2}\right)$. In the simplest consideration of the nearest neighbor sites backflow correlations in two dimensions, the prefactor $\gamma$ is about $4000$. As a result, the numerical cost of employment of backflow correlations is demanding when calculation system size as large as $16\times 16$. Therefore, we do not implement backflow correlations in our calculation.

To further improve the accuracy of the VMC calculation, one can apply the Lanczos methodVMCLanczos1993 . After obtained optimized wave function, we extend the wave function by multiplying the Hamiltonian as

[TABLE]

where $\left(1+\sum_{n=1}^{N}\alpha_{n}H^{n}\right)$ can be regarded as projection operator with variational parameters $\alpha_{n}$. One can minimize the energy by choosing appropriate $\alpha_{n}$. In principle, the accuracy can be systematically improved by increasing $N$, but the computational cost grows exponentially. Therefore, we employ Lanczos method up to the first step in our calculation.

II.3 variational Monte Carlo

We calculate the ground state by the variational Monte Carlo method with the variational wave function Eq. (8) provided in the last section. The expectation values of the Hamiltonian can be calculated by Monte Carlo samplings of the real space electron configurations $x$,

[TABLE]

where

[TABLE]

and $\rho$ is the weight in the importance sampling.

To find the ground state wave function, we optimize the variational parameters by the SR method(SR2001, ) in the variational Monte Carlo. The SR method starts from an approximate power method by the imaginary time evolution operatorIdo

[TABLE]

When $\tau$ is sufficiently small, it is reasonable to approximate $\left(1-\tau H\right)\left|\overline{\Psi}\right\rangle$ as a linear combination of the current wave function and its first derivatives,

[TABLE]

where $\left|\overline{\Psi}\right\rangle=\left|\Psi\right\rangle/\sqrt{\left\langle\Psi|\Psi\right\rangle}$ and $\left|\Psi^{k}\right\rangle$ is the derivative of normalized wave function with respect to a variational parameter $\alpha_{k}$

[TABLE]

To determine the coefficients $\gamma_{k}$, we need to minimize the cost function

[TABLE]

which can be obtained by computing the derivative with respect to $\gamma_{k}$ and Let

[TABLE]

Then, $\gamma_{k}$ can be found by solving the following coupled linear equation,

[TABLE]

Once we obtain the coefficients $\gamma_{k}$, the variational parameters are updated as follows

[TABLE]

Then, we repeat these steps until the energy converges.

To avoid local minimum, we gradually decrease the step width $\tau$ and randomize the update of every parameter

[TABLE]

where $\Delta t\left(i\right)$ is a gradually reduced function of SR step number $i$

[TABLE]

and $\eta\left(i,k\right)$ is a random number in the interval of $\left(r\left(i\right),1\right)$, with $r\left(i\right)=1-\left(1-r\left(0\right)\right)\nu^{i}$ to gradually reduce the randomness by selecting $\nu\in\left(0,1\right)$.

In order to optimize large number of parameters, one can solve the SR equation iteratively by conjugate gradient (CG) methodNeuscamman2012 . The detailed implementation is described in Appendix B. The computational cost for solving the SR equation by CG scales as $O\left(n_{s}n_{p}n_{\rm iter}\right)$ instead of $O\left(n_{s}n_{p}^{2}+n_{p}^{3}\right)$ needed for the solution of Eq.(19) by the simple matrix inversion, where $n_{s}$ is the number of samples, $n_{p}$ is the number of variational parameters and $n_{\rm iter}$ is the number of iterations in CG method. For comparison, the computational cost for the energy calculation scales as $O\left(n_{s}\left(N_{s}^{3}+\log_{2}\left(N_{s}\right)N_{s}D^{3}\right)\right)$, where $N_{s}^{3}$ comes from the calculation of Pfaffian, and $\log_{2}\left(N_{s}\right)N_{s}D^{3}$ comes from the contraction of tensor networks. Assume that we employ the full sublattice, namely, $N_{{\rm sub}}=N_{s}$, in the variational wave function, the number of parameters scales as $n_{p}\sim O\left(N_{s}^{2}+N_{s}D^{3}\right)$, where $N_{s}^{2}$ comes from the Jastrow factor and $\phi_{{\rm pair}}$, and $N_{s}D^{3}$ comes from the tensor network. Typically, we use $n_{\rm iter}\backsim O\left(10^{3}\right)$ in solving the SR equation by CG, so the dominant contribution of computational cost comes from solving SR equation. In the calculation of $16\times 16$ square lattice, the computational cost of the VMC part (with the CG) and the tensor network part are comparable when $D=8$.

III Benchmark results

To benchmark our method, we test the Hubbard model defined in Eq.(1).

III.1 1D Hubbard model

We first calculate the 1D Hubbard model with the periodic boundary condition (PBC). Figure 2 shows the $D$ dependence of the relative error of the ground state energy with respect to the exact result. The conventional MPS improves the accuracy about two orders of magnitude from $D=2$ up to $D=12$. However, even with $D=12$, the error is nearly $10^{-3}$. On the other hand, the conventional mVMC method ($D=0$) already achieves the error slightly less than $10^{-3}$ and the combination with the MPS substantially improves the accuracy with nearly the two orders of magnitude smaller error $\sim 10^{-5}$ for $D=12$. By combining with the mVMC, the accuracy of MPS is enhanced nearly two orders of magnitude.

III.2 2D Hubbard model

In the following, we calculate the 2D Hubbard model with the PBC. It is necessary to rely on the Monte Carlo sampling to combine the VMC and the tensor network procedure. When the Monte Carlo sampling is introduced to the tensor network part, the initial (reference) wave function $|\phi_{\rm ref}\rangle$ in the 2D case is required to be refined in advance as in the case of the VMC wavefunction. If we employ a simple state for the reference wavefunction such as $\sum_{i}|x_{i}\rangle$, which represents the equal-weight linear combination of all the real space basis function as employed in the conventional tensor network methods, the statistical error from the Monte Carlo sampling causes numerical instabilities. In other words, the tensor network calculation is made possible only by combining with the VMC if the Monte Carlo sampling is employed to reduce the computational cost of the contraction in the tensor network. Therefore, in the 2D case, we will not show the comparison with the solely FTTN result.

Figure 3 shows the $D$ dependence of the relative error of the ground state energy with respect to the exact result. We find that both of the combination of the TTN and FTTN with mVMC improve the accuracy of conventional mVMC results. In particular, the combination with the FTTN shows more significant improvement than that with the TTN. With the increase of the FTTN bond dimension $D$ up to 16, the accuracy of the mVMC is improved by nearly $1$ order of magnitude. Therefore, the combination of the FTTN and the mVMC provides a systematic method to improve the accuracy of each method applied separately.

In Fig. 4, we plot the energy as a function of the variance $\Delta_{{\rm var}}=\left(\left\langle H^{2}\right\rangle-\left\langle H\right\rangle^{2}\right)/\left\langle H\right\rangle^{2}$. Since the energy is linearly proportional to $\Delta_{{\rm var}}$ for sufficiently small varianceKashima2000 ; Kashima2001 ; SR2001 , we can perform the linear fitting to extrapolate to the energy of zero variance, so that more accurate ground state energy can be obtained. Figure 4 shows that the extrapolated energy agrees well with the exact result within the error bar of the linear fitting, which indicates the order of the relative error as small as $\sim 10^{-4}$. Though the accuracy is substantially improved, note that the strict variational principle satisfied before the extrapolation is not hold after the extrapolation, because of the possible extrapolation error.

We now perform the calculation for larger system sizes. Figure 5 shows the ground state energies for 2D Hubbard model with $U=4.0$ at half filling on $4\times 4$, $8\times 8$ and $16\times 16$ square lattices with the periodic-antiperiodic boundary conditions. We show comparisons among the mVMC, mVMC combined with FTTN, the first Lanczos step applied to the combined mVMC and FTTN, the variance extrapolation, and the QMC results. Our calculated energies agree well with the QMC results, which is expected to be practically exact. The relative error with respect to the QMC energiesMPQin2016 is about $0.3\%$ on the $16\times 16$ lattice size after the first-step Lanczos operation, and the extrapolated energies agree well with the QMC result within the error bar of the linear fitting. The first Lanczos step applied to the combined mVMC and FTTN, the variance extrapolation, and the QMC energies for the $8\times 8$ and $16\times 16$ square lattices are listed in Table 1.

IV doped Hubbard model on square lattice

In this section we show applications of the present method to the carrier doped Hubbard model on the square lattice to gain insight into a long-standing issue of the high-$T_{\rm c}$ superconductivity and severe competitions among various orders and fluctuations, given that the Hubbard model captures some essential physics of the high-$T_{\rm c}$ superconductivity. We show comparisons among the mVMC, mVMC combined with FTTN and the first Lanczos step applied to the combined mVMC and FTTN.

IV.1 Energy

Figure 6 shows the doping concentration dependence of the ground state energy on $16\times 16$ square lattice for $U=10.0$. The energy difference between periodic-periodic boundary and periodic-antiperiodic boundary results are negligible on $16\times 16$ lattice size (See upper inset of Fig. 6), so we have performed all the calculation with the periodic-antiperiodic boundary condition in this section. We employ the $16\times 2$ sublattice in the variational wave function to be compatible with the possible stripe orders. The result shows that the combination of the FTTN substantially lower the energy of the corresponding mVMC result. We have performed the optimization from a homogeneous $d$-wave superconducting state and a stripe order coexisting with $d$-wave superconducting order, and we find that the stripe ordered state coexisting with the weak $d$-wave superconductivity provides lower energy, while the state optimized from the homogeneous $d$-wave superconducting state stays metastable as an excited state at least for $\delta<0.25$. The energies of the first Lanczos step applied to the combined mVMC and FTTN, and the variance extrapolation, which are obtained from stripe and homogeneous initial states, are listed in Table 2.

In Fig. 7, we show the same data as the main panel of Fig. 6, but in an enlarged scale of the vertical axis after subtracting a linear function $f(\delta)$ to emphasize the difference among different methods. Note that the phase separation region determined by drawing the tangent line from the point at $\delta=0$ as shown in the dashed line for the VMC+FTTN+Lanczos data suggests that the phase separation for $0\leq\delta\lesssim 0.1$, which is narrower than the phase separation region $0\leq\delta\lesssim 0.2$ in Ref. Misawa2014PRB, , while below $\delta\sim 0.1$, the survival of the phase separation is robust. The difference of the present result from Ref. Misawa2014PRB, is ascribed to the fact that the present calculation allows the finite-period phase separation, namely the stripe order, replacing the phase separation, which is regarded as the “infinite period” charge order. The charge/spin stripe order suppresses the $d$-wave superconducting correlation as we see below.

IV.2 $d$-wave Superconducting Correlation

In Fig. 8, we show the superconducting correlations defined as

[TABLE]

where

[TABLE]

and $f_{d_{x^{2}-y^{2}}}$ is the $d_{x^{2}-y^{2}}$ superconducting pairing symmetry factor

[TABLE]

where $\mathbf{r}=\left(r_{x},r_{y}\right)$. Both of the results for the ground states with the coexisting stripe and superconductivity as well as the excited states with charge uniform superconducting states presented in Figs. 6 and 7 are shown. The both results show that the combination of the FTTN slightly suppresses the superconducting correlation on the long range part in comparison to the mVMC results.

Since the ground state we obtained has a stripe order, we show the superconducting correlation along $x$ and $y$ directions separately for the ground state in Fig. 9. We see that the combination of the FTTN suppresses the superconducting correlation both along $x$ and $y$ directions with factors two to three, and the superconducting correlation along the $y$ direction, which is the stripy direction, shows larger long range correlation than along the $x$ direction, indicating that the charge modulation suppresses the superconductivity as expected. In contrast the long-ranged part of the amplitude of the superconducting order is more than one order of magnitude larger for the charge homogeneous excited states. It should be noted that, even for the stripe direction in the hole rich region, the superconducting correlation is much lower than the case of the charge homogeneous states. Since the long-ranged part of $P_{d}$ is the square of the order parameter, the order parameter $\langle\Delta_{d}\rangle$ is more than factor 3 larger for the charge-uniform excited states. It suggests that the superconductivity can be substantially enhanced from the equilibrium ground state if one can keep the metastable charge-uniform state. Nevertheless, despite weak order, the charge-inhomogeneous ground state preserves the saturated long-ranged correlation particularly in the stripe direction ($y$ direction). In the direction crossing the stripe ($x$ direction), the correlation shows the long-ranged saturation to a smaller value with oscillation with the period of the charge stripe. Since the charge-stripe long-range order may be sensitively destroyed by the randomness such as that caused by the distribution of the dopant atoms in the real compounds of the cuprate superconductors and may be replaced with domain structures, the long-range superconducting order may further be suppressed than the values in the present ideally periodic stripe order.

Here in Fig. 10, we show the superconducting correlation by taking the origin $\mathbf{r}_{i}$ in Eq. (23) at the maximum and minimum values of the hole density (namely, columns of the stripe with smallest (largest) electron densities). It clearly shows that the superconducting correlation along the stripe direction stays large at long distance only for the hole-rich columns while that at the hole poor region is extremely small. This is because the hole-poor region is essentially Mott insulating like (see the data for the charge density discussed below). The correlation in the $x$ direction shows an oscillatory behavior and it confirms that the hole-rich superconducting strips are bridged each other essentially by the mechanism of the Josephson junction through the hole-poor strips, which is the reason why the superconducting order stays smaller than the charge homogeneous state.

IV.3 Spin and Charge Correlations

To identify the stripe order in the ground state, we show the spin structure factor

[TABLE]

where $s_{\mathbf{r}}^{z}=n_{\mathbf{r}\uparrow}-n_{\mathbf{r}\downarrow}$. We also present the charge structure factor

[TABLE]

where $n_{\mathbf{r}}=n_{\mathbf{r}\uparrow}+n_{\mathbf{r}\downarrow}$.

Figure 11 shows the spin and charge structure factor at $\delta\sim 0.11$. The peak of spin structure factor is at $\left(\frac{7\pi}{8},\pi\right)$, and the peak of charge structure factor is at $\left(\frac{\pi}{4},0\right)$, which indicates that the ground state has a stripe order with $(l_{c},l_{s})=(8,16)$, where $l_{c}$ ($l_{s}$) denotes the charge (spin) wave length in a stripe phase.

Then, we compare the stripe order obtained by different numerical methods in Fig. 12 and Fig. 13. Our calculation shows that the combined mVMC, FTTN and first Lanczos step very slightly lower the peak of the spin structure factor in comparison to the VMC+FTTN results. The difference between the VMC and VMC+FTTN results is also small.

Although the true ground state of a finite-size system must preserves the translational symmetry, in our calculated results of the stripe ordered states, the translational symmetry is explicitly broken if the momentum projection is not imposed, because the system size is fairly large. To show the spin and charge stripe patterns in the real space, we have computed the local spin density along $z$ direction

[TABLE]

and the local charge density

[TABLE]

which are shown in the color scale plot of the spin and charge configuration in the lower panel of Fig. 12 and Fig. 13, respectively. We note that the stripe order has the amplitude as large as 0.2, implying that the charge modulation extends from the Mott insulating density ($\delta\sim 0$) to $\delta\sim 0.2$. In the realistic condition with the potential randomness and long-range Coulomb interaction, this amplitude of the long range order may be weakened.

In the overdoped region, the ground state may show a different stripe order. We calculate the spin and charge structure factor at $\delta\sim 0.22$, which are shown in Fig. 14. The peak of spin structure factor is at $\left(\frac{3\pi}{4},\pi\right)$, and the peak of the charge structure factor is at $\left(\frac{\pi}{2},0\right)$, which indicate that the ground state has a stripe order with $(l_{c},l_{s})=(4,8)$. The reduction of the period with the increasing hole concentration is intuitively understood from the deceasing mean hole distance with doping and also consistent with the experimental indications of the diffuse peak in neutron scattering in the cuprates Yamada1998 ; Dai1998 .

From Figs. 15 and 16, we can see that the combination of the mVMC with the FTTN and first Lanczos step provide nearly the same structure factor. The spin and charge stripe patterns can be seen in the color scale plots of the spin and charge configurations in the lower panels of Fig. 15 and Fig. 16, respectively.

In order to show clear comparison, we pictorially depict the spin and charge orders at $\delta\sim 0.11$ and $\delta\sim 0.22$ in Fig. 17.

To emphasize the difference between the charge inhomogeneous ground state, and the charge-homogeneous and superconducting excited state, we here show the charge and spin correlations of the charge-homogeneous excited state for an example at $\delta\sim 0.11$ and $U=10$ in Fig. 18. The charge structure factor does not have an appreciable peak confirming the charge homogeneity, while the spin correlation has a peak at the commensurate wave number $(\pi,\pi)$. The coexistence of the superconductivity with the antiferromagnetic order around $\delta\sim 0.1$ was already found for the charge homogeneous state Misawa2014PRB .

V Summary and discussion

We have proposed a method, which combines the tensor network and the variational Monte Carlo method by taking advantage of the both to study fermionic lattice models. In order to perform fast contraction and preserve the lattice symmetry, we have introduced the FTTN into the variational wave function. Our calculation shows that this combined method substantially improves the accuracy in comparison to the accuracies separately achieved by the conventional VMC calculation and the tensor network.

Tensor network states usually satisfy entanglement area law which may become inefficient to capture the large amount of entanglement in itinerant fermionic systems. Recently, there exists an attempt to alleviate the limitation of the area law by the PEPS simulations with the help of the energy extrapolationPEPSextrapolate and it was applied to the issue of the competing phases in the doped Hubbard modelstripeU8 The mVMC provides a more flexible reference wave function instead of the basis of real space product states, so that the combination with the mVMC extends the power of the tensor network algorithms particularly for highly entangled correlated metals. Detailed comparisons with the single PEPS algorithm about the accuracy are left for future study.

We have applied the present method to study the ground state of the hole doped Hubbard model. The ground states show coexisting stripe-type charge and spin orders and weak $d$-wave superconducting order in the lightly doped region. The results indicate that the strong superconducting order realized by imposing spatial homogeneity becomes an excited state and the true ground state is replaced by the state with the superconducting correlation substantially weakened by the emergence of the charge/spin stripe order with the period depending on the doping concentration, which can be seen from Fig. 8. This compromised ground state shows the universal feature of the strong-coupling superconductivity that is subject to the spatial inhomogeneity including the phase separation and charge/spin order. The high-$T_{{\rm c}}$ superconductivity in the strongly correlated electron systems needs to overcome the simultaneous charge inhomogeneity that weakens the superconductivity. The excitation energy of the charge-uniform superconducting state is very small and roughly around 0.01 with small doping concentration dependence (see Fig. 7).

Ref. stripeU8, has studied Hubbard model at $U=8$ with $\frac{1}{8}$ doping by various kinds of latest numerical methods, which is consistent with the existence of stripe orders found here. Throughout our calculation on the doped Hubbard model, we fix at $U=10$, and we leave the $U$ dependence for future work to further study the stability of the stripe order.

We have studied only the lattices with the sizes $2^{n}\times 2^{n}$. Our calculation on $16\times 16$ lattice shows that the ground states below $\delta\sim 0.12$ and above $\delta\sim 0.2$ show the stripe orders with $(l_{c},l_{s})=(8,16)$ and $(l_{c},l_{s})=(4,8)$, respectively, while roughly in the region $0.12<\delta<0.2$ the energy is convex implying the phase separation. However, one can speculate that other periodicities or structures of the stripes with the period between (8,16) and (4,8) that are not compatible with this size may have slightly lower energy filling the convexity like $t$-$J$ modelStripetJ and precludes the phase separation in this region. From the energy curve in Fig. 7 the phase separation is expected as well in the region $0<\delta<0.12$, which looks more robust. Systematic studies along this line are left for future studies. In addition, analyses on different lattice structures as well as effects of the intersite Coulomb interaction are intriguing future issues. In particular the intersite interaction may substantially change the behavior of the charge modulation.

Acknowledgements.

We would like to thank Naoki Kawashima, Frank Pollmann and Ying-Jer Kao for stimulating discussion. The authors thank the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo for the facilities. This work was financially supported by Japan Society for the Promotion of Science through Program for Leading Graduate Schools (MERIT), the MEXT HPCI Strategic Programs for Innovative Research (SPIRE), the Computational Materials Science Initiative (CMSI) and Creation of New Functional Devices and High-Performance Materials to Support Next-Generation Industries (CDMSI). We thank the computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (under the project number hp130007, hp140215, hp150211, and hp160201). This work was also supported by a Grant-in-Aid for Scientific Research (No. 22104010, No. 22340090 and No. 16H06345) from Ministry of Education, Culture, Sports, Science and Technology, Japan.

Appendix A Tree tensor network state

In this appendix, we will describe the tree tensor network state. To reduce computational cost, we employ a binary tree in this work. As an example, we consider a square lattice with $N_{s}=L\times L$ sites, where each site has a local Hilbert space with dimension of $d$. As shown in Fig. 19 for $8\times 8$ lattice, the TTN is composed of a set of tensors $t_{j,i}$, where $j=1,2,\cdots,R\left(R=\log_{2}N_{s}\right)$, and $i=1,2,\ldots,N_{s}/2^{j}$. The TTN is connected as a binary tree structure, which can be expressed as

[TABLE]

where $t_{1,i}$ is the leaf tensor which directly contains 2 physical indices (2 sites) and 1 virtual index, and $t_{2,i},t_{3,i},\ldots,t_{R,1}$ are internal-node tensors which only contain virtual indices.

The representation of a wave function with a TTN on real space basis can be interpreted as real space coarse-graining transformation. Each layer of node tensors reduces the effective lattice size by a factor of 2, so the height (the number of the hierarchical levels) of the tree is $\log_{2}\left(L\times L\right)$. Differently from coarse-graining transformation, the node tensors in the TTN are not necessarily isometric.

The structure of the hierarchy of the node tensors are equivalent to that of elimination tournament play. The total number of node tensors in the binary tree is $N_{s}-1$. Since there is no loop in the TTN, exact contraction is feasible. If we employ the Monte Carlo sampling on the real space configuration, the physical indices of the leaf tensors are fixed so every leaf tensor becomes a vector. Therefore, we can start from the contraction of the vector at the leaf tensor and then continue the contraction of rank 3 tensors at the higher hierarchical levels, of which the computational cost scales as $O\left(N_{s}D^{3}\right)$, where $D$ is the dimension of the virtual indices. Since we only move one or two electrons in every Monte Carlo step, we can reuse the intermediate result from the tensor network contraction on the previous configuration. Therefore, the computational cost can be reduced to $O\left(\log_{2}\left(N_{s}\right)D^{3}\right)$.

Appendix B Stochastic reconfiguration with conjugate gradient

To solve the SR equation Eq. (19) by a direct method, one should calculate the overlap matrix $S_{kl}=\left\langle\Psi^{k}|\Psi^{l}\right\rangle$. By introducing

[TABLE]

$S_{kl}$ can be calculated as

[TABLE]

where $j$ is the sample index, and $\rho\left(j\right)$ is the weight function for the importance sampling defined as

[TABLE]

The diagonal elements of $S$ matrix is the sample variance of $\mathcal{O}$

[TABLE]

and the variance of $S_{kk}$ is given by

[TABLE]

If $\text{Var}\left(S_{kk}\right)/S_{kk}$ is larger than a threshold, we truncate the $k$-th parameter in the SR equation to stabilize the optimization.

According to Eq. (32), the time cost for explicit construction of $n_{p}\times n_{p}$ $S$ matrix is $O\left(n_{s}n_{p}^{2}\right)$, and the time cost for solving the SR equation is $O\left(n_{p}^{3}\right)$, while the memory cost is $O\left(n_{p}^{2}\right)$, where $n_{s}$ is the number of samples and $n_{p}$ is the number of variational parameters. This is both the main time and memory consuming part, when the number of variational parameters becomes large.

To reduce the cost, one can solve the SR equation iteratively by conjugate gradient (CG) method, so that the explicit construction of $S_{kl}$ is not required(Neuscamman2012, ). We only need to realize the matrix-vector multiplication $\sum_{l=1}^{n_{p}}S_{kl}\gamma_{l}$, which can be calculated by the Monte Carlo sampling as

[TABLE]

To reduce the memory cost, one can change the summation order of the variational parameters and samples as

[TABLE]

As a result, the computational cost is reduced from $O\left(n_{s}n_{p}^{2}+n_{p}^{3}\right)$ to $O\left(n_{s}n_{p}n_{\rm iter}\right)$, and the memory cost is reduced from $O\left(n_{s}n_{p}+n_{p}^{2}\right)$ to $O\left(n_{s}n_{p}\right)$. By this application of the CG method, it allows an efficient way of solving Eq. (19).

Appendix C Fast update of Pfaffian

In each Monte Carlo sample of real space configuration

[TABLE]

one can calculate $\left\langle x|\phi_{{\rm pair}}\right\rangle$ by the computation of Pfaffian of the matrix

[TABLE]

If a new sample $x^{\prime}$ is proposed by changing the position of one electron in $x$, then one can calculate ${\rm Pf}\left(B\right)=\left\langle x|\phi_{{\rm pair}}\right\rangle$ by

[TABLE]

where matrices $A$ and $B$ have same elements except for $\alpha$-th row and column. From Eq. (40), one can update Pfaffian at the cost of $O\left(M\right)$, where $M$ is the dimension of the matrix.

This fast update of Pfaffian is only possible if the inverse of $A$ is known, but the inverse need only to be directly computed once, as it can be efficiently updated when one row or one column of the matrix changes by using the Sherman-Morrison formulaShermanMorrison1949 ,

[TABLE]

where $u$ and $v$ are column vectors, where the computational cost of Eq. (41) is $O\left(M^{2}\right)$. The detailed explanation of the update technique can be found in Ref. SMorita2015, However, during the update of matrix $A$, the round-off error accumulates. As a result, the computation of matrix inverse would be inaccurate, so that the fast update of Pfaffian would contain large error.

To solve this problem, we estimate the condition number by matrix 1-norm

[TABLE]

where the 1-norm of matrix is defined as

[TABLE]

If the condition number becomes larger than threshold $\kappa\left(A\right)>\kappa_{c}$, the Pfaffian of $A$ should be computed from scratch, instead of using fast update procedure. The computational cost of 1-norm scales as $O\left(M^{2}\right)$, which is the same as the computation cost of the matrix inverse by the Sherman-Morrison formula, so that the estimate of the condition number does not increase the order of the cost in the update of the Pfaffian.

Appendix D Rescaling of the variational parameters

If the condition number of overlap matrix $S_{kl}=\left\langle\Psi^{k}|\Psi^{l}\right\rangle$ is large, the solution of the SR equation would be inaccurate. Therefore, to stabilize the optimization, we shift the diagonal elements in $S$

[TABLE]

where $\epsilon$ is a small constant.

For a real space sample $x$, the derivative of every tensor can be expressed as

[TABLE]

where $t_{\mu}^{i}$ is the $\mu$-th element of the $i$-th tensor. The ${\rm Tr}$ in the denominator is to trace out virtual indices of all tensors, and the ${\rm Tr}^{\prime}$ in the numerator is to trace out virtual indices of all tensors except $t^{i}$, and $\mu\left(\nu\right)$ is the simplified notation of all virtual indices of one tensor. If $t^{i}$ is rescaled by a factor of $\eta^{i}$

[TABLE]

the wave function is unchanged, but the derivative of the corresponding tensor has been rescaled as

[TABLE]

If inappropriate rescaling factor $\eta^{i}$ is employed, the derivative of every tensor could be significantly different in the order of magnitude. Therefore, the shift of diagonal elements in $S$ as Eq. (43) will suppress the change of parameters in tensors which have small derivatives.

To change parameters efficiently while stabilizing the optimization in the SR procedure, we determine the rescaling factor according to the amplitude of the derivative of each tensor as

[TABLE]

instead of rescaling the parameters according to the amplitude of tensor elements as

[TABLE]

As a result, the derivative of every tensor would be more or less in the same order of magnitude, so that SR method can optimize elements of every tensor efficiently.