Ground-state properties of the spin-1/2 antiferromagnetic Heisenberg model on the triangular lattice: A variational study based on entangled-plaquette states

We study, on the basis of the general entangled-plaquette variational ansatz, the ground-state properties of the spin-1/2 antiferromagnetic Heisenberg model on the triangular lattice. Our numerical estimates are in good agreement with available exact results and comparable, for large system sizes, to those computed via the best alternative numerical approaches, or by means of variational schemes based on specific (i.e., incorporating problem dependent terms) trial wave functions. The extrapolation to the thermodynamic limit of our results for lattices comprising up to N=324 spins yields an upper bound of the ground-state energy per site (in units of the exchange coupling) of $-0.5458(2)$ [$-0.4074(1)$ for the XX model], while the estimated infinite-lattice order parameter is $0.3178(5)$ (i.e., approximately 64% of the classical value).


Introduction
In the last decades quantum Monte Carlo (QMC) [1] has emerged as one of the most powerful numerical tools to study the ground-state (GS) properties of quantum spin systems on a lattice. It has been successfully applied, for example, to determine quantities of crucial interest such as the energy and the order parameter of the two dimensional (2D) antiferromagnetic Heisenberg model (AHM) on the square lattice [2] and, more generally, can be considered the method of choice to investigate unfrustrated (bosonic) models in any dimension.
The situation is more complicated when the geometry of the system is frustrated as in the triangular or kagomé lattice. In this case QMC techniques based on imaginary time projection, which are "exact" for unfrustrated systems, suffer from a sign-problem and approximations or alternative techniques have to be employed. The density matrix renormalization group (DMRG) [3] method, although not affected by sign instability, yields very accurate results only for one dimensional (1D) spin systems, being instead hardly applicable to general 2D problems due to the unfavorable scaling of the computer resources needed with the system size [4].
Projected-entangled pair states (PEPS) [5], the natural extension to higher spatial dimensions of the matrix product states variational family, which constitutes the foundation of DMRG, can efficiently approximate GS's of local Hamiltonians [6], and have been successfully employed to simulate frustrated and bipartite quantum spin systems [7,8,9,10,11]. Their applicability, however, is restricted by the disadvantageous scaling of the computational cost to 2D models with open boundary conditions. Different approaches, which have been proposed, and applied in 2D [12,13,14,15,16,17,18,19], suffer from similar limitations and seem hardly usable for large systems [possibly three dimensional (3D) ones] and periodic boundary conditions (PBC). The idea of combining DMRG and Monte Carlo [20,21] has opened the opportunity to design new variational algorithms for Tensor-Network based ansatze [22,23]. The applicability of one of them, based on string-bond states, has been demonstrated even in the case of 3D systems [22].
In a recent paper we have introduced a class of variational states, namely entangled-plaquette states (EPS), which provides an accurate ansatz to simulate the GS properties of quantum spin models [24]. We have shown [24] that the EPS trial wave function (WF) allows to describe efficiently the GS of lattice Hamiltonians taking advantage of a simple variational Monte Carlo algorithm. Extremely accurate results have been obtained by us [24] and other investigators [25,26] for several frustrated (fermionic) and unfrustrated models.
In this paper we study, as a paradigmatic example of frustrated system, the spin-1 2 AHM on the triangular lattice. Previous works, carried out by means of both QMC approaches with a fixed-node approximation and variational calculations based on different classes of states, which explicitly incorporate antiferromagnetic 2π 3 Nèel order, found values of the energy per site extrapolated to the thermodynamic limit between ∼ −0.53 and ∼ −0.5470 [27,28,29,30,31,32]. Evidence of long-range order has been generally found with important differences, however, in the estimate of the order parameter. Our work, is motivated, therefore, not only by the interesting physics which can emerge from the GS of a quantum antiferromagnet but also by methodological aspects: it is in fact important, in order to asses its applicability to generic quantum spin models, to benchmark the general EPS ansatz against other variational choices specifically designed to tackle a particular problem.
We carried out computer simulations on systems of size up to N = 324 spins (previous applications of EPS on frustrated systems have been performed on much smaller system sizes, up to N = 64) [24]. We found for N = 36, in the case of the XX model, an error on our estimated energy relative to the exact one [33] as small as 6 × 10 −4 and an energy per site extrapolated to the thermodynamic limit of -0.4074 (1). At the Heisenberg point we obtained E EP S 0 (N = 36) = −0.55420(4), which differs from the exact value by ∼ 1%, and E EP S 0 (∞) = −0.5458(2), with a value for the order parameter of 0.3178(5).

The EPS wave function
In this section we describe the entangled-plaquette ansatz offering also some detail concerning the procedure adopted to minimize the energy and compute GS physical observables. Given a collection of N spins 1 2 arranged on a lattice, a generic trial WF can be expressed as: where the W (n)'s are complex weights, n = |n 1 , n 2 , . . . , n N and n i = ±1 ∀ i = 1, 2, . . . , N is the eigenvalue of σ z i . The variational principle ensures that the energy expectation value on a trial state is a rigorous upper bound of the GS energy, therefore, such a quantity can be evaluated by minimizing the trial energy with respect to the weights. To achieve this task one has to specify the form of W (n), which corresponds to the choice of a particular variational ansatz. A desirable ansatz is one able to capture the GS physics of a system by means of few variational parameters, which can be efficiently optimized on a computer. The entangled-plaquette ansatz is based on the following idea: i) Cover the systems with M plaquettes in such a way that l sites labelled by n 1,P , n 2,P , . . . , n l,P = n P belong to the P th plaquette.
ii) Express the weight of a global spin configuration as a product of variational coefficients C n1 1 , C n2 2 , . . . , C nM M in biunivocal correspondence to the particular spin (e.g., along the z axis) configuration of the plaquettes. Hence: In order to gain additional insight into the EPS ansatz let us discuss few exemplificative cases: a) Each site is a plaquette: the wave function is of the mean field form (i.e., correlations are neglected). b) Plaquettes comprise more than one site and each site belongs to a single plaquette: with this choice it is possible to obtain reasonable estimates of the GS energy (up to corrections which scale with the plaquette boundaries) and short-range correlations (of the order of the plaquette size). c) Overlapping (i.e., entangled) plaquettes are used: this ansatz clearly produces, accounting for the entangled nature of a quantum GS, estimates of the energy and long-range correlations much more accurate than those computable with a WF whose weights are the product of non overlapping plaquettes [34,35] as in a) and b).
A variety of variational ansatze, as illustrated in [25], have a straigthforward representation in terms of EPS. This is, for example, the case of the WF introduced by Huse and Elser [27] for the AHM on the triangular lattice which contains up to three-spin correlations and a phase factor that depends on the system geometry. The accuracy of this variational choice can promptly be improved via a general (i.e., with no specific phase factor) EPS WF based on larger or differently shaped plaquettes, hence, including additional correlations. In our calculations we start from the trial state described in a) (i.e., N = M and l = 1) where the weights and the initial spin configuration in the computational basis are randomly selected; we estimate, via Monte Carlo sampling, the expectation value of the energy as well as its derivatives with respect to the C nP P 's; the value of these coefficients is then optimized by updating all of them of a small step against the gradient direction. Once convergence is reached the size of each plaquette is increased, naturally introducing correlations in the ansatz, and the procedure iterated. Finally, obtained the optimal GS WF, the mean value of observables other than the energy can be evaluated by means of the Monte Carlo method as well. An exhaustive discussion on the particular suitablility of EPS to efficiently perform the minimization procedure previously outlined is reported in [24]. It is important to mention that the entangled-plaquette ansatz, as well as the proposed numerical protocol, is applicable to systems in any spatial dimension; moreover our method, as any purely variational (i.e., not based on imaginary time projection techniques) one, is sign-problem free allowing for the characterization of both unfrtustrated and frustrated systems.

Results
In the following, we show results obtained with EPS for the 2D AHM on the triangular lattice. Computer simulations have been performed on several system sizes comprising up to N = 324 sites. Our numerical findings are compared with exact diagonalization (ED) results available for N = 36 [33] and estimates obtained with alternative computational approaches or variational ansatze [27,28,29,30,31,32]. The AHM Hamiltonian is: where J xy = J z = 1, S α i is the spin-1 2 operator in the α direction acting on site i and the summations run over nearest-neighbour sites. The calculations were carried on in the S z T OT = N i S z i = 0 sector, being the total magnetization along z a good quantum number. Therefore, during the simulation, we update the system configuration proposing a spin-flip for couples of spins with opposite S z . As a first benchmark we estimate the GS energy of the XX model (i.e., J z = 0, J xy = 1) with PBC. Numerical results, shown in table 1, refer to different size l of the plaquettes. For l = 4 the relative error of the EPS energy of the N = 36 system with respect to the exact result is smaller than 2%, and decreases down to 6.5 × 10 −4 when l = 16. into account the statistical uncertainty, from the results obtained when N = 324. For this size, therefore, the thermodynamic limit appears fully recovered. The EPS extrapolated energy is slightly lower than that computed by extrapolating exact results up to N = 36 (i.e., E ED 0 (N = ∞) = −0.4066), illustrating the importance of studying system sizes larger than those treatable via ED techniques. It has to be mentioned that the small deviation shown by our data in figure 1 from a linear scaling, being the same regardless the size of the plaquettes, should not be due to a loss in accuracy of the EPS WF as the system size increases; more plausibly, it can be ascribed to higher order corrections in the series expansion of the energy as a function of 1 N . Moreover the energy estimates (for a fixed lattice size) do not show evidence of convergence with the plaquette sizes that we have been able to utilize in this work. This is not surprising (in fact the EPS ansatz is formally exact when the plaquette size is as big as that of the whole system) and simply means that the exact result, even if the energy differences for l > 9 are quite small, is not yet reached. As a consequence one can expect, on further increasing the plaquette size, an improvement of the GS energy. It is however remarkable, taking into account the methodological aspect of our paper, that our best variational upper bounds, obtained by using plaquettes of 16 sites, already provide extremely accurate energy estimates.
Next, we discuss our findings for the GS of the Hamiltonian (3) at the Heisenberg point (i.e., J xy = J z = 1). GS energies per site are reported in table 2. The EPS GS energy computed for N = 36 differs from the exact one [33] by ∼ 1%, while, calculations based on ordered WF's yield typically higher values of the GS energy between ∼ −0.5367 [27] and ∼ −0.5515 [31]. When N increases, the general EPS WF remains competitive with respect to accurate choices specifically adapted to the triangular lattice geometry. Particularly, our energies for N = 144 and N = 324 differ from the best (to our knowledge) variational estimates [31] by less than 0.25%. Table 2. Variational GS energy per site (in units of J xy ) computed, by using the entangled-plaquette ansatz, for the spin-1 2 AHM on a triangular lattice of size N with PBC. Error bars are in parenthesis; N plaquettes comprising 16 sites have been employed. The exact result for N = 36 is also shown for comparison [33].  (2), which is in excellent agreement with the most accurate results derived in previous works [29,31].
GS estimates of the spin-spin correlation function Corr EP S 0 (r MAX , N ) = S 0 · S rMAX EP S computed at the maximum distance (taking into account the PBC) on the lattice are reported in table 3. The magnetic order parameter of the system defined as m 2 = lim N →∞ Corr 0 (r MAX , N ) is a quantity of central importance, being directly related to the presence of long-range order. Wether or not the spin-1 2 AHM displays long-range order has been a long lasting-question, however, the positive answer, on the basis of analytical [36,37] and numerical [27,28,29,30,31,33,38,39] results, is commonly accepted. Estimates of the order parameter expressed as a fraction of its classical value range from ∼ 0.4 [29,38] to ∼ 0.72 [30]. Our computed spin-spin correlations for N = 144 and N = 324 are indistinguishable (within the statistical uncertainty). Hence we can consistently assume m EP S = (Corr EP S 0 (r MAX , N = 324))    (3) approach for the AHM on the square lattice [24], is consistent with the long-range ordered nature of the system found in previous numerical studies [29,37,38,28,30].

Conclusions
The general entangled-plaquette wave function has been employed to study the ground-state properties of the spin-1 2 antiferromagnetic XX and Heisenberg Hamiltonian on the triangular lattice. With this ansatz and a simple variational Monte Carlo algorithm for the energy minimization, we found extremely accurate energy estimates for finite systems, as well as, by extrapolating, in the thermodynamic limit. The accuracy of our results compares well with that provided by the best alternative methods or variational ansatze, which, when possible, (as in the case of the triangular lattice), are adjusted, to improve the estimates, by incorporating physical features of the particular system of interest. It is crucial to stress that the trial state utilized in this work, although not including any problem dependent term, provides a reasonable estimate even of the system order parameter; therefore, the variational class of entangled-plaquette states can be regarded as a reliable option to describe the GS of general frustrated (i.e., not treatable with standard quantum Monte Carlo methods) physical models on a lattice with no need of any a-priori knowledge.