Ground-state properties of quantum many-body systems: entangled-plaquette states and variational Monte Carlo

We propose a new ansatz for the ground-state wave function of quantum many-body systems on a lattice. The key idea is to cover the lattice with plaquettes and obtain a state whose configurational weights can be optimized by means of a variational Monte Carlo algorithm. Such a scheme applies to any dimension, without any ‘sign’ instability. We show results for various two-dimensional spin models (including frustrated ones). A detailed comparison with available exact results, as well as with variational methods based on different ansatze, is offered. In particular, our numerical estimates are in quite good agreement with exact ones for unfrustrated systems, and compare favorably to other methods for frustrated ones.


I. INTRODUCTION
The study of the ground-state (GS) properties of quantum many-body systems is one of the most challenging tasks of theoretical physics.Exact results can be obtained numerically only for systems of relatively small size (i.e., few particles); this limitation is particularly severe, e.g., when studying phase transitions, wherein the emergence of long-range order can only be established by carrying out an extrapolation of the estimates to the thermodynamic limit.Two main numerical techniques, namely Density Matrix Renormalization Group (DMRG) 1 and Quantum Monte Carlo (QMC) 2 have successfully been employed to investigate quantum spin models on a lattice.These computational approaches, however, find their optimal applicability under different specific constraints.DMRG yields extremely accurate results in one dimension (1D) even for very large systems, but fails in describing the properties of the quantum GS in higher dimension due to the unfavorable scaling with the system size of the computational resources needed. 3MC, on the other hand, is the method of choice for quantum systems obeying Bose statistics in any dimension, but suffers from the notorious "sign problem" for Fermi systems.
Generalizations of the variational family of Matrix-Product States (MPS) underlying the DMRG have recently been investigated to go beyond the discussed limitations of DMRG itself and QMC.Specifically the most natural extension of MPS is given by Projected-Entangled Pair States (PEPS) 4 which efficiently approximate ground states of local Hamiltonians 5 and have been used to simulate, otherwise intractable, 2D quantum systems. 6,7Other variational families of states have also been proposed and tested in 2D. 8,9,10,11,12,13,14Despite the promising results obtained so far, it seems very hard to use those methods in 3D [or even for systems with periodic boundary conditions (PBC)] due to the unfavorable scaling of the required computer resources.
A new possibility has recently emerged to combine the main advantages of DMRG and Monte Carlo in order to build new algorithms. 15,16Some of them 15,17 can be used in 2D and it seems that the one based on String-Bond States 15 may be used even for some 3D systems. 18n this paper we introduce a new numerical technique that combines the strengths of QMC with an extension of PEPS to simulate lattices models, overcoming some of their limitations.We also test this technique with non-trivial models, and compare it with other techniques, including PEPS.Specifically, we propose a new class of states called Entangled-Plaquette States (EPS) which allow an accurate characterization of the GS of quantum spin systems by means of a simple Variational Monte Carlo (VMC) algorithm (see Ref. 19 for a general review).
The basic idea underlying EPS can be schematically described as follows: Assume to cut a lattice in several sub-blocks (e.g.small PEPS) and extract for each sub-block the GS wave function.The wave function of the original system, expressed as the product of the subblock wave functions yields reasonable (up to corrections which scale with the sub-block boundaries) estimates of GS energy and short-range (of the order of the sub-block size) correlations.These estimates dramatically improve if the sub-block size is increased and, more importantly, if overlapping sub-blocks (i.e., entangled plaquettes) are employed.The latter is the crucial point which allows, accounting for its correlated nature, a description of the quantum GS much more accurate than that obtainable with a simple non-overlapping-plaquette product state (i.e., in a mean-field fashion). 20A GS wave function whose weights are the product of variational parameters in one to one correspondence to the spin configuration of each entangled plaquette, appears, consequently, the natural choice for a variational ansatz.Therefore, our numerical approach is based on a family of states, namely EPS, which share many analogies with PEPS, and takes advantage of Monte Carlo sampling to estimate physical observables of interest.It can be applied to systems of any spatial dimensionality and, as a pure variational method (i.e., not involving imaginary time projection), is sign problem free.
We test our numerical protocol on a variety of quantum spin models on a square lattice comprising as many as N = L×L = 400 sites.Our energy estimates are in excellent agreement with those ("exact" in practice) obtained by QMC for a system of lattice hard core bosons.In the presence of nearest-neighbor repulsion at the Heisenberg point, we find an extrapolated (to infinite lattice size) value of the energy per site which differs from the QMC result 21 by less than 2×10 −3 , and is more accurate than VMC estimates obtained with a Jastrow wave function. 22n the case of a frustrated antiferromagnet (i.e., the socalled J 1 − J 2 model), for which a sign problem exists in QMC, our energy estimates (whose error relative to the exact ones is less than 1.5×10 −2 ) compare favorably with those obtained with PEPS, or fixed-node Green Function Monte Carlo (GFMC). 23

II. METHODOLOGY
Consider a collection of N spins 1 2 arranged on a L × L square lattice with N = L 2 (without loss of generality, here and in the following we will refer to this specific case).Provided a trial state |ψ = n W (n)|n where |n = |n 1 , n 2 , . . ., N and n i = ±1 ∀ i = 1, . . ., N , the energy expectation value on the given state is: where and W (n) = W * (n) (real weights are assumed for simplicity).According to the variational principle, E is an upper bound of the GS energy which can be evaluated by minimizing Eq. 1 with respect to the weights.At this point we have to make an ansatz for the wave function: imagine to cover the lattice with plaquettes (say one of dimension l 1 × l 2 for each site) and assign a coefficient C nP P to all the possible 2 l1×l2 spin configurations of any single plaquette.Given a global spin configuration |n , its weight can be expressed as follows: where C nP P depends only on the spin state of the n P sites belonging to the P th plaquette.With this choice the analytic expression of the derivative of Eq. 1 with respect to C nP P is: The multidimensional summations in Eq. 1 and 4 can exactly be evaluated only for small N .For large systems (N > ∼ 30) , one has to employ the Monte Carlo method.Specifically, the energy, as well as its derivatives, can be estimated from the same sample.Moreover, the only quantity depending on the plaquette coefficient with respect to which the derivative is taken is . By using Eq. 3 it turns out that D P (n P ) is simply equal to 1/C nP P (easy and fast to compute).Similarly to other works, 15,16 , the steps of the basic variational algorithm adopted to calculate the GS energy are: i) Start from a randomly chosen initial configuration; ii) generate a large set of new configurations by flipping one or more spins via the Metropolis algorithm; 24 iii) evaluate the energy and its gradient vector; iv) update all the C nP P s of a small step against the gradient direction; v) iterate from ii) until convergence of the energy is reached.It is worth mentioning that expectation values of physical observables other than the energy can be evaluated according to Eq. 1 and 2 when H is replaced by a generic operator O.
For a single-spin flip the acceptance probability is given by: where the flip is proposed for the j th spin and the products run over all the plaquettes which include such a spin.In a typical calculation, we start with 2 × 2 plaquettes; once the energy has converged, the size of the plaquettes is increased to improve the estimate.The number of coefficients which need be stored in memory for a spin-1/2 system is N × 2 l1×l2 , and can be reduced taking into account problem-dependent symmetries.For each optimization step, a few thousands updates are necessary to get a rough estimate of the energy and efficiently move the coefficients along the gradient.At the later stages of the simulation, to reach the optimal energy value, an important role is played by the gradient step which has to be carefully tuned.An example of how the error of the GS energy relative to the GFMC result decreases as a function of the plaquettes size is illustrated in Fig. 1.
The relative error is already small (less than 1%) for 2 × 2 plaquettes and reaches a value < 10 −3 when 4 × 4 plaquettes are used.Numerical data refer to a system of lattice bosons (at half filling) which interact via an infinite on-site (hard core) repulsion on a 10 × 10 lattice with PBC.Since the total magnetization along z is a good quantum number, we performed the calculation in the canonical ensemble (i.e., in the S z = 0 sector).Consequently we chose to update the configuration by FIG.1: (color online) Dependence on the plaquettes size of the error in the GS energy (computed with the method illustrated in this work) relative to the GFMC result for a system of hard core bosons at half filling on a 10 × 10 square lattice.PBC are assumed.The dashed line is only a guide to the eye.
flipping pairs of spins j and k for which S z j = −S z k .The expression of the acceptance probability for the pair update is a straightforward generalization of Eq. 5.

III. RESULTS
Estimates of the GS energy per site of a system of lattice hard core bosons, for three different lattice sizes, are presented and compared to GFMC results in Tab.I.Although the GFMC method is not purely variational (i.e., the GS wave function is projected out from a trial state via imaginary time evolution) our numerical data are in excellent agreement with GFMC ones, even for the largest system considered in this work (L = 20).As a second case study we investigate the J 1 −J 2 quantum spin Hamiltonian: where the first (second) summation runs over nearest-(next-nearest-) neighbor sites.For J 1 = 1 and J 2 = 0 the above Hamiltonian describes the antiferromagnetic Heisenberg model, which is isomorphic to a system of hard core bosons with nearest-neighbor repulsion of strength V =2t.This model has been extensively studied numerically in the past (see for instance Refs.21 and  22).GS energies per site for various lattice sizes are presented in Tab II and compared to exact calculations 25 and Stocastic Series Expansion (SSE) data.The energy per site is a monotonic increasing function of the lattice size.The simple formula 26 where β = 1.4377, 27 gives an extrapolated energy per site E ∞ = −0.6683(3),which is in good agreement with the most accurate QMC result: 21 E SSE ∞ = 0.669437(5) and lower than other extrapolated values obtained by purely variational methods.For example, on the basis of a Jastrow wave function, Trivedi and Ceperley found E JST ∞ ≃ −0.6590. 22alculations carried on with open boundary conditions (OBC) yield GS energies lower than those obtained with the general PEPS, or SBS method.For example, we get E = −0.6258(1)for L = 10 while the PEPS result is E P EP S = −0.62515and the SBS one E SBS ∼ −0.6225.
Estimates of the spin-spin correlation function computed at the maximum distance on the lattice according to the formula Corr(L/2, L/2) = S r • S r ′ where r − r ′ = (L/2, L/2) are shown in Tab III.We find an extrapolated value of the staggered magnetization defined by   26 It has to be mentioned that the discrepancy between our result and the SSE one might be due to a violation of the socalled "area law". 28y including the next-nearest-neighbor interaction (J 2 > 0) in the J 1 − J 2 Hamiltonian the system becomes frustrated and is believed to undergo a phase transition for J 2 ≃ 0.6.This model cannot be simulated by QMC due to the sign problem (arising in turn from the underlying Fermi statistics).Our method, instead, is purely variational and can be applied without the occurrence of any sign instability.With the aim of comparing our estimates with exact results (available only for N ≤ 36), we computed the GS energy, as a function of J 2 /J 1 , of the J 1 − J 2 model on a 6 × 6 square lattice.Numerical values are shown in Tab.IV; the error of our estimates relative to the exact result is shown in Fig. 2.This quantity never exceeds 1.5×10 −2 , moreover it has to be  mentioned that the GS energies computed by means of EPS compare favorably to those obtained by GFMC with the fixed-node approximation 23 and, except in a narrow region of J 2 /J 1 ∼ 0.5, where the the agreement is however remarkable, to variational results obtained with a BCS-type ansatz. 29,30lso in this case, for OBC, we test our scheme against the PEPS one.The estimated energies are lower than those computed with PEPS even for a 8 × 8 lattices (see Tab. V), where the PEPS approach performs at its best.

IV. CONCLUSIONS AND OUTLOOK
In this work we have presented Entangled-Plaquette States: an ansatz for the GS wave function of quantum many-body systems which allows the accurate estimate of physical observables by means of Variational Monte Carlo.Our approach not only gives accurate results for unfrustrated systems (where other methods are in principle "exact" ) but, most importantly, applied to systems for which QMC simulations suffer from a sign problem, yields estimates whose accuracy seems to be not obtainable with different techniques.Therefore our method appears as a very promising avenue to investigate the effects of Fermi statistics (e.g., frustration).The extension of our computational approach to 3D can be easily implemented (i.e., taking cubic plaquettes) and several possible improvements (e.g., cover the lattice with plaquettes of different shapes and sizes), are being currently investigated.

1 FIG. 2 :
FIG.2:(color online) Error (relative to the exact result) of the GS energy of the J1 −J2 model computed with the method illustrated in this work.The square lattice comprises N = 36 spins; PBC are assumed.The dashed line is only a guide to the eye.

TABLE I :
GS energy per site (in units of the nearest-neighbor hopping integral t) for a system of hard core bosons on a L×L lattice with PBC.For L = 20, 4 × 3 plaquettes have been employed, 4×4 ones in the remaining cases.GFMC estimates have been also computed and are shown for comparison.

TABLE II :
GS energy per site of the antiferromagnetic Heisenberg model on a L × L square lattice with PBC.SSE estimates and exact results are also shown for comparison.Numerical values are given in units of J1.

TABLE III :
Spin-spin correlation function (computed at the maximum distance on the lattice) of the antiferromagnetic Heisenberg model on a L × L square lattice with PBC.SSE estimates (adjusted for different factors in the definition) are also shown for comparison.Numerical values are given in units of J1.

TABLE IV :
GS energy per site (in units of J1) for the J1 − J2 model on a square lattice (with PBC) comprising 36 sites.Exact results are also shown for comparison.

TABLE V :
GS energy per site (in units of J1) for the J1 − J2 model on a square lattice (with OBC) comprising 64 sites.PEPS results are also shown for comparison.
a data provided by V.Murg