A cluster-based mean-field, perturbative and coupled-cluster theory description of strongly correlated systems
Athanasios Papastathopoulos-Katsaros, Carlos A. Jiménez-Hoyos, Thomas M. Henderson, Gustavo E. Scuseria
AA cluster-based mean-field, perturbative and coupled-cluster theorydescription of strongly correlated systems
Athanasios Papastathopoulos-Katsaros, a) Carlos A. Jiménez-Hoyos, Thomas M. Henderson,
1, 3 and Gustavo E.Scuseria
1, 3 Department of Chemistry, Rice University, Houston, Texas 77005, USA Department of Chemistry, Wesleyan University, Middletown, Connecticut 06459, USA Department of Physics and Astronomy, Rice University, Houston, Texas 77005, USA (Dated: 19 February 2021)
We introduce cluster-based mean-field, perturbation and coupled-cluster theories to describe the ground state ofstrongly-correlated spin systems. In cluster mean-field, the ground state wavefunction is written as a simple tensorproduct of optimized cluster states. The cluster-language and the mean-field nature of the ansatz allows for a straight-forward improvement based on perturbation theory and coupled-cluster, to account for inter-cluster correlations. Wepresent benchmark calculations on the 2D square J − J Heisenberg model, using cluster mean-field, second-orderperturbation theory and coupled-cluster. We also present an extrapolation scheme that allows us to compute ther-modynamic limit energies very accurately. Our results indicate that, even with relatively small clusters, the correlatedmethods can provide an accurate description of the Heisenberg model in the regimes considered. Some ways to improvethe results presented in this work are discussed.
I. INTRODUCTION
Spin lattices and more specifically, Heisenberg models, areof significant chemical importance. For example, iron-sulfurclusters relevant to nitrogen fixation or photosynthesis, suchas ferredoxins, have been treated according to the Heisenbergmodel. Single molecule magnets have possible applicationsin quantum computers as the smallest practical unit for mag-netic memory. These molecules are usually metal clusters,and the magnetic coupling between the spins of the metal ionscan also be described by a Heisenberg Hamiltonian. Lastly,electrides, conjugated hydrocarbons and a few superconduc-tors have some of their features modelled after Heisenberg ex-change interactions.
A common feature of all those systemsis that they are strongly-correlated.Despite tremendous effort and progress, the accurate andefficient description of the ground state of strongly-correlatedsystems represents an open problem in quantum chemistry.The defining feature of strongly-correlated systems is theirmulti-reference nature, which makes single-reference meth-ods inadequate. Accordingly, approaches based on compos-ite particles have been proposed for treating these systems,e.g. resonating valence bond as a ground state candidate forhigh- T c superconductors, suggested by Anderson. In this work, which is a continuation of the work in Ref. 8,we use composite many-spin cluster states to describe theground state of strongly-correlated spin lattices. More specif-ically, we divide the lattice into clusters, each containing apredefined number of sites which in this work are chosen us-ing proximity in the lattice. These clusters can have any shapeand any size, although for 2D systems, compact shapes pro-vide better results. The cluster states are a subset of all theavailable many-particle states (most often all the states of the S z = a) Electronic mail: [email protected] as a product of cluster states. It is important to note that themany-particle state in each cluster is determined in the pres-ence of other clusters. The resulting cluster mean-field (cMF)state is variational. The optimization provides not only theoptimal cMF state, but also a renormalized Hamiltonian ex-pressed in term of cluster states. Traditional many-body ap-proaches can then be used, on this renormalized Hamiltonian,to account for the missing inter-cluster correlations.In related work, Isaev, Ortiz, and Dukelsky considered,in their hierarchical mean-field (HMF) approach, a similaransatz to ours for the same Hamiltonian. Our approach dif-fers from that used in Ref. 9 in not requiring the individ-ual clusters to share the same ground state. That is, theground state of each cluster is optimized independently allow-ing for (translational and spin) symmetry-broken solutions.In addition, we here consider two common approaches fromquantum chemistry, Rayleigh-Schrödinger perturbation the-ory (RS-PT) to second-order and coupled-cluster (CC) asa means to obtain a correlated approach defined in terms ofclusters. Our coupled-cluster approach is inspired by Li’s block-correlated coupled-cluster method. In that work, how-ever, the ground state of each cluster was not optimized inthe presence of other clusters as we do here. Block-correlatedcoupled cluster has been used with high success in quantumchemistry to describe strongly-correlated molecular systemsusing a complete active-space reference state. In addition,Mayhall went beyond cMF by implementing a selectiveconfiguration interaction framework and had very accurate re-sults for fermionic systems. Lastly, a cluster product stateis also connected with tensor network (TN) techniques thathave been gaining popularity for treating strongly-correlatedsystems. For more information on this connection and theadvantages of cMF in contrast to other more sophisticated ap-proaches, we refer the reader to Ref. 8.There are many positive aspects that cMF posseses, whichwe would like to point out. First of all, straightforward sym-metry breaking ( S symmetry) can partially account for inter-cluster correlations. Moreover, the fact that the cluster mean- FIG. 1. Nearest-neighbor ( J ) and next-nearest neighbor ( J ) inter-actions. field state constitutes the ground state of a mean-field Hamil-tonian of which the full set of eigenstates can be easily con-structed, allows us to exploit already developed techniquesfrom the single-reference family of methods, in order to ac-count for the missing inter-cluster correlations.Our objective with this work is two-fold. First, we presentthe cMF formalism and provide details of the RS-PT (cPT2)and the coupled-cluster formulation we use (cCCSD). Oursecond objective is to apply these techniques to a sim-ple strongly-correlated system: the Heisenberg model in a2D square lattice. It is important to underline that the1D case is exactly solvable, so we will focus our atten-tion on the 2D model. This model has received numer-ous studies in the past two decades, using various meth-ods such as exact diagonalization, coupled-cluster, density-matrix renormalization group (DMRG), matrix-product or tensor-network based algorithms, resonatingvalence bond (RVB) and quantum Monte Carlo (QMC). We compare our results with calculations from Ref. 20, whichwe use as a reference. Those results are not exact in the ther-modynamic limit, because they are extrapolated from differ-ent 2D shapes (more information in Ref. 20), but we considerthem useful for a semi-quantitative comparison. Our resultsshow that cPT2 and cCCSD significantly improve upon cMFand can provide an accurate description of the ground state ofthe square J1-J2 Heisenberg lattice.The rest of the article is organized as follows. In Sec. II wepresent the formalism behind cMF, cPT2 and cCCSD. SectionIII provides some practical computational details regardingthe calculations presented in this work. In Sec. IV we presentthe results of cMF, cPT2 and cCCSD calculations for the 2Dsquare J − J Heisenberg model. A brief discussion follow-ing the results is presented in Sec. V, as well as some ideas asto how to improve the approaches presented here. Lastly, Sec.VI is dedicated to some general conclusions.
II. FORMALISMA. Heisenberg model
In this work, we focus our attention on the J - J Heisen-berg model in two-dimensions on a square grid. The Heisen-berg model describes a collection of spins in a lattice (of finite
FIG. 2. Néel (left) and collinear (right) antiferromagnetic phases ofthe square J − J Heisenberg model. In between, there is a non-magnetic phase whose exact form is debated. size L) interacting through the Hamiltonianˆ H = J ∑ (cid:104) i j (cid:105) (cid:126) S i · (cid:126) S j + J ∑ (cid:104)(cid:104) i j (cid:105)(cid:105) (cid:126) S i · (cid:126) S j (1)where (cid:126) S i is the spin- operator on site i , J and J are thenearest-neighbor and the next-nearest neighbor coupling co-efficients respectively (see Fig. 1), and the notation (cid:104) i j (cid:105) im-plies interaction among nearest-neighbors, while (cid:104)(cid:104) i j (cid:105)(cid:105) im-plies interaction among next-nearest neighbors. In the fol-lowing, we confine ourselves to the antiferromagnetic (AFM)case J , J > (cid:46) J / J (cid:46) . J . In J / J (cid:38) .
6, the ground state displays anAFM phase with collinear long-range order character due tothe dominance of the next-nearest-neighbor coupling J (seeFig. 2). In the regime 0 . (cid:46) J / J (cid:46) .
6, the Néel and thecollinear orders compete. The nature of this intermediateground state is still a much debated issue, as arethe type of the phase transitions and the transition points. Be-cause cMF has already been used to shed light on these topics(see Ref. 39) we will not focus our attention on these issues,but we need to mention there is a second-order transition fromthe Néel to the paramagnetic phase, whereas there is a first-order transition from the paramagnetic to the collinear anti-ferromagnetic phase.
B. Cluster Mean-Field
Our formalism of cluster mean-field (cMF) is based onRef. 8, but in this paper we confine ourselves to the spincase, which consists of a subset of configurations found infermionic systems. For more details we refer the reader tothat work, but below we present the general framework.Let the lattice states be grouped, according to some crite-rion (such as proximity in real space), into clusters of size l , l , ... l n , where n is the number of such clusters. Formally,the Hilbert space of each cluster is of size 2 l i , as in each site wecan either have a spin-up or a spin-down. We choose to workwith eigenstates of S z in each cluster, thereby reducing the ef-fective dimension of the Hilbert space. Similarly to Ref. 27,the Hilbert space of the full system is simply given by the ten-sor product of the Hilbert spaces of all clusters.A second-quantized formulation in terms of cluster productstates can also be established. Let A † I , c ( A I , c ) create (annihilate)the I-th state in cluster c. This I-th state is a linear combinationof many-spin basis states (possibly mixing states with differ-ent S z ) constructed as products of the states in the cluster. Weformally write | I (cid:105) c = A † I , c |−(cid:105) c (2)where | I (cid:105) c is the vacuum state in cluster c, a useful abstractconstruct.Each cluster product state is formally built as | I (cid:105) | J (cid:105) ... | Z (cid:105) n ≡ A † I , A † J , ... A † Z , n |−(cid:105) (3)where |−(cid:105) is a vacuum state for the full system. In this work,we consider a cluster product (mean-field) state as a varia-tional ansatz for the ground state wavefunction. That is, theansatz | Φ (cid:105) for the ground state is given by | Φ (cid:105) = | (cid:105) | (cid:105) ... | (cid:105) n (4)where the 0 label indicates the ground state of each clusterin the presence of the other clusters. The lowest energy cMFstate is obtained by a variational minimization scheme, as out-lined in Ref. 8.Defining excited configurations is straightforward. We canwrite them as | Φ Ii (cid:105) = | (cid:105) ... | I (cid:105) i ... | (cid:105) n , (5) (cid:12)(cid:12) Φ Ii ; J j (cid:11) = | (cid:105) ... | I (cid:105) i ... | J (cid:105) j ... | (cid:105) n , (6) (cid:12)(cid:12) Φ Ii ; J j ; Kk (cid:11) = | (cid:105) ... | I (cid:105) i ... | J (cid:105) j ... | K (cid:105) k ... | (cid:105) n (7)for singly-, doubly-, and triply-excited clusters.Before proceeding further, let us comment on the nature ofthe cluster product states considered in this work. We indi-cated above that the ground state of each cluster is expressedas a linear combination of the many-spin basis states in it. Theexpansion over those states can be restricted and for our pur-poses, we confine ourselves to cases with a specific S z . This isdone in order for the cluster product state | Φ (cid:105) to be an eigen-function of S z and to reduce the dimension of the ground statevector in each cluster. Lastly, we choose to break S in eachcluster to develop long-range antiferromagnetic ordering. C. Matrix elements and cMF optimization
The evaluation of the matrix elements is again similar toRef. 8, but an important difference that we need to point outis that the spin cluster Hamiltonian has 1- and 2-cluster ele-ments but not 3- or 4-cluster elements. This difference arises because of the form of the Hamiltonian, which does not havemore than 2-spin interactions. This significantly simplifies theprocedure for both the cMF optimization, as well as the cPT2and cCC extensions mentioned later. For details regarding thecMF optimization, we refer the reader to Ref. 8. We note thatthe zero-th order cluster Hamiltonian for a given number offermions is given byˆ H c = ∑ pr ∈ c (cid:104) p | ˆ t | r (cid:105) a † p a r + ∑ pqrs (cid:104) pq | ˆ v | rs (cid:105) a † p a † q a s a r + ∑ pr ∈ c a † p a r ∑ c (cid:48) (cid:54) = c ∑ qs ∈ c (cid:48) ρ c (cid:48) sq ( (cid:104) pq | ˆ v | rs (cid:105) − (cid:104) pq | ˆ v | sr (cid:105) ) (8)where ρ c (cid:48) sq is the one-particle density matrix in cluster c (cid:48) . Wechoose to perform the optimization self-consistently, in orderto minimize the energy. The formula above was given in itsfermionic form, since we expect most readers will be morefamiliar with fermionic Hamiltonians than spin Hamiltonians.Lastly, we remind the reader that this is a spin system so thereis no orbital optimization to be considered.We think one final comment contrasting cMF and standarddiagonalization techniques is necessary. Let us focus on thecase of a finite lattice with periodic boundary conditions. Thescaling of cMF with respect to cluster size is similar to thatof exact diagonalization performed on a full lattice of thesame size as the cluster. cMF, however, has a larger prefactorsince the equations need to be solved self-consistently and theansatz explicitly breaks the translational (and, possibly, spin)symmetry of the lattice. Thus, we can reach cluster sizes com-parable to those achievable with exact diagonalization results(though of course cMF can have many clusters). In this work,the largest cMF calculation reported used a 6x6 cluster, andthe length of the eigenvector in each tile is ∼ . Note alsothat we can use less costly approximate methods in each clus-ter instead to reduce the scaling of cMF. Also note that cMFallows us to find different solutions for different regimes. D. Perturbation theory
Once again, the theory is very similar to Ref. 8, but sig-nificantly simplified. Generally, in RS-PT, the second-ordercorrection to the ground state energy is evaluated as E ( ) = ∑ µ (cid:54) = | ˆ V µ | ε − ε µ (9)where ˆ V = ˆ H − ˆ H and V µ = (cid:104) Φ | ˆ V | µ (cid:105) . Here, µ labels theeigenstates of ˆ H and ε µ are the corresponding eigenvalues.The excited states framework was explained in section B. Asmentioned earlier, the evaluation of the matrix elements is eas-ier compared to the fermionic case, because there are no 3-and 4-cluster interactions, so the cost of cPT2 is O ( n M ) ,where M is the number of excited states in each cluster. Wewould like also to remind the reader that we use all the statesin the Hilbert space of each cluster (although in practice forcPT2 we only need states with S z = m , m + , m −
1, where m is the S z value used in cMF), even though in the cMF optimiza-tion only one S z sector of the Hilbert space was considered. Inaddition, due to the structure of the Hamiltonian, we only havetwo types of 2-tile interactions ( S z -preserving and S z mixing),as opposed to 5 (see Ref. 8) in the fermionic case. E. Coupled-cluster theory
In this paper, our aim is to go beyond perturbation theoryto a coupled-cluster framework. Our work bears connection toRef. 12, but it is important to note that in that work, the groundstate of each cluster was neither optimized in the presence ofother clusters nor tested on the same model.The cluster coupled-cluster (cCC) expansion of the ground-state wave function can be written in an intermediate normal-ized form as follows: | Ψ (cid:105) = e ˆ T | Φ (cid:105) (10)where ˆ T is the cluster operator. In this work, we focus ourattention to singles and doubles, thereforeˆ T = ˆ T + ˆ T (11)and in the cluster language, ˆ T and ˆ T operators can be writtenas ˆ T = ∑ i ∑ I ( i ) t I i B + I i B − i (12)ˆ T = ∑ i ∑ j (cid:54) = i ∑ I ( i ) ∑ J ( j ) t I i J j B + I i B − i B + J j B − j (13)where the coefficients t I i and t I i J j are the single and doubleamplitudes respectively, the operator B + I i excites cluster i tostate I and operator B − i de-excites cluster i to its ground state.By projecting onto | Φ (cid:105) and the space of singly anddoubly excited configuration functions, and by utilizingthe fact that for spin systems we have up to 2-cluster ex-citations, the cCCSD equations naturally truncate and become E cCCSD = (cid:104) Φ | H | Φ (cid:105) + (cid:104) Φ | H | T Φ (cid:105) + (14a) (cid:104) Φ | H (cid:12)(cid:12)(cid:12)(cid:12) ( T + T ) Φ (cid:29) E cCCSD t I i = (cid:104) Φ I i | H | Φ (cid:105) + (cid:104) Φ I i | H | T Φ (cid:105) + (14b) (cid:104) Φ I i | H (cid:12)(cid:12)(cid:12)(cid:12) ( T + T ) Φ (cid:29) + (cid:104) Φ I i | H (cid:12)(cid:12)(cid:12)(cid:12) ( T T + T ) Φ (cid:29) E cCCSD ( t I i J j + t I i t J j ) = (cid:10) Φ I i ; J j (cid:12)(cid:12) H | Φ (cid:105) + (cid:10) Φ I i ; J j (cid:12)(cid:12) H | T Φ (cid:105) + (14c) (cid:10) Φ I i ; J j (cid:12)(cid:12) H (cid:12)(cid:12)(cid:12)(cid:12) ( T + T ) Φ (cid:29) + (cid:10) Φ I i ; J j (cid:12)(cid:12) H (cid:12)(cid:12)(cid:12)(cid:12) ( T T + T ) Φ (cid:29) + (cid:10) Φ I i ; J j (cid:12)(cid:12) H (cid:12)(cid:12)(cid:12)(cid:12) ( T + T T + T ) Φ (cid:29) where eq. 14a is the equation for the energy, and 14b and 14care for the single and double amplitudes respectively. If weinsert eq. 12 and eq.13, we get a set of nonlinear equationsand like in traditional coupled-cluster, we solve the equationsself-consistently.It should be emphasized that in a cCCSD calculation thecomputation of the last term in the right-hand side of eq.14cis the most time consuming step, which makes the cCCSDmethod computationally a O ( n ) procedure, where n is thenumber of clusters.Lastly, similar to the conventional truncated CC expansion,the truncated cCC expansion is also size extensive. For moreinformation on the proof, we encourage the reader to checkRef. 12. III. COMPUTATIONAL DETAILS
The cMF, cPT2 and cCCSD calculations presented in thiswork were carried out with a locally prepared code. In allthe calculations with even number of spins in each cluster,we use the same number of up and down spins, whereas inall the cases with odd number of spins, we use +/- numberof up and down spins in order to be able to construct Néeland collinear antiferromagnetic phases respectively. In thiswork, we choose to work with eigenstates of S z in each clusterfor computational convenience. At this point, it is importantto note that the S z of each cluster is not a symmetry of theHamiltonian, so this choice is a constraint on the cMF. Thefull relevant S z sector of Hilbert space within each cluster wasused in constructing the cluster ground state | (cid:105) c . For smallcluster sizes, the ground state in each cluster was found bya standard diagonalization of the local cluster Hamiltonian.For larger cluster sizes, a Davidson algorithm was used tosolve for the ground state. In cPT2 and cCCSD calculations,we loop over all the relevant excited cluster states (althoughin principle the number of states can be truncated, which wasdone in Ref. 8). For solving the cCCSD equations, the tradi-tional CC iteration scheme was used. IV. RESULTS
In this section we present results of cMF, cPT2 and cCCSDcalculations on the 2D Heisenberg model. We start by pro-viding the basic idea in Sec. IV A, where we dive into somedetails regarding the optimization of cMF states and the wayin which other results are presented. In Sec. IV B we showand compare the results from cPT2 to cMF calculations. InSec. IV C, we show and compare our cCCSD results to thecMF and cPT2 ones. These results are compared to accuratenumerical estimates from Ref. 20.
A. Basic idea
In this section we discuss most of the aspects regardingthe optimization of cMF states. In this way, we hope thatthe results presented in subsequent sections will become moretransparent to the reader. We consider rectangular Heisenbergperiodic lattices with J , J ≥
0. As a first step, we choose thecorresponding tiling scheme (i.e., defining how the spins aregrouped into clusters). The optimized state in the cluster isexpressed as a linear combination of all the possible (combi-natorial number) resulting configurations and it is optimizedby a self-consistent diagonalization of the appropriate clusterHamiltonian.It is important to mention that we are interested in the ther-modynamic limit properties (very large system sizes). This isachieved by having both the cluster size l → ∞ and the systemsize L → ∞ . For this work, however, only a limited numberof sizes was tried, because the exact diagonalization in largeclusters becomes very expensive (recall that we require notonly the ground state in each cluster but also, in principle, allexcited states). However, we tried an extrapolation schemeand we computed results to the thermodynamic limit ( l → ∞ ).Lastly, we define the thermodynamic limit system size corre-sponding to a specific cluster size ( L → ∞ for a fixed l ) as thethe smallest system size for which the results do not changeif we increase it further. It is rigorously shown ‡ that for theHeisenberg Hamiltonian and clusters of even dimensions, thethermodynamic limit system size of cMF is twice as large asthe cluster size in each dimension of the rectangle, and forcPT2 and cCCSD, it is three times larger. B. cMF results
We start by considering the cMF results. It is importantto note that some of these results are already published inRef. 39. All calculations in this section were performed in pe-riodic square lattices. Only uniform tiling schemes were con-sidered; clusters were rectangles of l lattice sites, each filledwith l spins. All the results are assumed to be at the thermo-dynamic limit L → ∞ for a fixed l , as explained in the pre-vious section. We note that broken-symmetry cMF solutionscan be achieved, that is, a non-zero magnetization develops ‡ This is based on the fact that the cMF solution is uniform (same wavefunctionon every cluster) for the specific model for clusters of even dimensions. Also,cPT2 and cCCSD only correlate neighbor tiles; there are no Hamiltonian ma-trix elements between separated tiles. Therefore, it follows that cPT2 andcCCSD yield the thermodynamic limit energy with a 3x3 supracluster lattice. M agne t i z a t i on o f a c en t e r s i t e J /J cMF-2x2cMF-2x4cMF-2x6cMF-2x8cMF-4x4 FIG. 3. The magnetizations of a center site for different cluster sizesare shown. The Néel (on the left) and the collinear (on the right)antiferromagnetic phases can easily be observed as there is strongmagnetization. In between, there is a non-magnetic phase, which isdepicted by the zero magnetization. on each lattice site. Regarding the collinear case and the rect-angular clusters, the collinear correlations are observed in thesmaller dimension. This was chosen because its energy waslower than the case with collinear correlations in the larger di-mension. We present in Fig. 3 the magnetizations of a centersite (for the 2x2 case every site is equivalent) in cMF calcu-lations for different values of l and J / J . For readers notfamiliar with spin systems, magnetizations are equivalent tospin densities. For recognizing the critical points, we com-pared the paramagnetic and the antiferromagnetic solutionsand saw whether their energy difference changes sign. Dif-ferent cluster sizes correspond to different critical points, soan extrapolation scheme is needed to compute the phase tran-sition points of the thermodynamic limit ( L → ∞ and l → ∞ ).This was performed in Ref. 39 and it was shown that cMF wascomparable to other methods in the field.In Figs. 4-6 we demonstrate the energy per site obtainedin cMF calculations for different values of l and J / J in thethermodynamic limit ( L → ∞ ). Fig. 4 describes the conver-gence of 2xN to the 2x ∞ limit. It is interesting that in thiscase the second critical point is at ∼ L we notice that there is no sharptransition to a paramagnetic regime. Lastly, Fig. 5 describesthe convergence of 4xN to the 4x ∞ limit. In this case, thevalues predicted for the second transition point are the largest.We emphasize the significance of the choice of the cluster,which can be pointed in two different cases. First, the 2x2cluster results are significantly different from all the other re-sults, which suggests that at the cMF level, 2x2 clusters arenot sufficiently large to capture all the physics of the system(luckily this is not the case for the correlated methods, as weshall see later). Second, comparing the 2x8 to the 4x4 clusters,we notice that each phase can be better studied by a differentcluster shape, although the total number of sites is the same. -0.65-0.6-0.55-0.5-0.45-0.4 0 0.2 0.4 0.6 0.8 1 E ne r g y pe r s i t e ( E / N ) J /J cMF-2x2cMF-2x4cMF-2x6cMF-2x8cMF-2x10cMF-2x12Extrapol. cMF-2xinf FIG. 4. Energies per site at the cMF level for different clusterschemes are shown, with every cluster having one dimension equalto 2. cMF-2x2 stands for cluster mean-field with 2x2 tiles, etc. Thepurpose of this figure is to show the convergence of 2xN to the 2x ∞ limit (black circles), extrapolated as described in Sec. III B. -0.65-0.6-0.55-0.5-0.45-0.4 0 0.2 0.4 0.6 0.8 1 E ne r g y pe r s i t e ( E / N ) J /J cMF-4x4cMF-4x6cMF-4x7 FIG. 5. Energies per site at the cMF level for different clusterschemes are shown, with every cluster having one dimension equalto 4. cMF-4x4 stands for cluster mean-field with 4x4 tiles, etc. Thepurpose of this figure is to describe the convergence of 4xN to the4x ∞ limit. This probably occurs because of the long-range order of thecollinear phase, which is better depicted in the 2x8 case. Oneimportant observation to point out is that the paramagnetic so-lution is exactly equivalent to the energy of a single tile withopen boundary conditions: ie, in cMF the inter-cluster inter-action vanishes as there is no magnetization along the clusterboundaries.In this work, we also went one step beyond and tried to ex-trapolate the cMF energy to the thermodynamic limit ( L → ∞ and l → ∞ ). In order to achieve that, we plotted the cMFenergy with respect to 1 / N (note that N = l / x , where x isone of the dimensions of the cluster). Figures 7 and 8 showa nearly linear behavior as the size of the cluster increases, -0.65-0.6-0.55-0.5-0.45-0.4 0 0.2 0.4 0.6 0.8 1 E ne r g y pe r s i t e ( E / N ) J /J cMF-2x2cMF-3x3cMF-4x4cMF-5x5 FIG. 6. Energies per site at the cMF level for different square clusterschemes are shown. cMF-2x2 stands for cluster mean-field with 2x2tiles, etc. The purpose of this figure is to describe the convergence ofthe NxN to the ∞ x ∞ limit. which can also be explained rigorously * . The J / J = 0 limithas been well studied with careful extrapolations carried outby several methods (see Ref. 20 and references therein). Ourresults are summarized in Tab. I. It is shown that they gen-erally agree with other methods in the field (Ref. 25, 32) andthe agreement with reference calculations (Ref. 20) is accu-rate to (0.001 - 0.002) x J . Lastly, an extrapolation of themagnetization of the central site was attempted using a 1 / N behavior and square lattices, similar to the energy. The result-ing extrapolated magnetization (for J / J =
0) in the thermo-dynamic limit (both L → ∞ and l → ∞ ) is 0.320, which is ingood agreement with 0.310 from Ref. 20. Clusters cMF (E/N)2x ∞ -0.624453x ∞ -0.638364x ∞ -0.645845x ∞ -0.650356x ∞ -0.65341 ∞ x ∞ (1) -0.66873 ∞ x ∞ (2) -0.66789Reference Ref. 20 -0.67010CCM Ref. 25 -0.66936QMC Ref. 32 -0.66944 * The difference between a 2x6, a 2x8 and a 2x10 tile lies in the addition ofextra "internal" sites in the tile, ie., the boundaries have very similar magne-tizations. If this is the case, then the error in the energy of the 2xN vs 2x ∞ should behave linearly, as observed in the plot. On the other hand, the linearnature of the NxN extrapolation with respect to 1/N can be explained becausethe error is proportional to the surface (or in this case the perimeter) of thecluster. −0.65−0.64−0.63−0.62−0.61−0.60−0.59−0.58 c M F ene r g y pe r s i t e ( E / N ) FIG. 7. Energy per site obtained in cMF calculations for J / J = N . It is shown thatthe behavior is nearly linear, so the extrapolation scheme is reliable.Accurate estimates of the energies for the 2x ∞ , 3x ∞ , 4x ∞ , 5x ∞ and6x ∞ clusters can be computed. −0.66−0.64−0.62−0.60−0.58 c M F ene r g y pe r s i t e ( E / N ) NxinfNxN
FIG. 8. Energy per site obtained in cMF calculations for J / J = N . Nx ∞ energiesrefer to the extrapolated energies depicted in Fig. 7 and NxN refer toenergies computed with square clusters. The extrapolated results at(1 / N = 0) were obtained by fitting a line using the two largest clustersizes in each series and can be compared with extrapolated referencecalculations from Ref. 20 (see Tab. I.)TABLE I. Energy per site obtained by extrapolating cMF calculationsat J / J = ∞ x ∞ (1) is the result from extrapolatingNx ∞ → ∞ x ∞ , and ∞ x ∞ (2) is the result from extrapolatingNxN → ∞ x ∞ . The extrapolated result from exact diagonalizations(reference calculations) from Ref. 20 is shown for comparison, aswell as coupled-cluster (CCM) and quantum Monte Carlo (QMC)results. The agreement with Ref. 20 is accurate to (0.001 - 0.002) x J . C. cPT2 results
We continue by considering the cPT2 results. All the com-putational parameters are the same as for the cMF results. Weremind the reader that all the possible excited states were usedto compute the cPT2 energy, except for the 4x4 case wherefewer states ¶ were used. It was ensured, however, that the en-ergy was converged to at least 4 decimals, compared to whenusing all the states. The criterion for choosing those stateswas their corresponding eigenvalues, from which we chose thelowest ones. In Fig. 9 we demonstrate the energy per site ob-tained for different values of l and J / J . Even though cPT2is not variational, there is evidence that the exact energy islower than the cPT2 energy, because of the results of Ref. 20.There are three important observations. First, increasing thecluster size is not as important for cPT2 as for cMF. Thiscan be useful for real systems, because for cCCSD, we donot have to use large clusters, whose cost can be prohibitive.We have to underline, however, that for our case, cPT2 andcCCSD should converge to the exact answer as the size of thecluster increases. Second, the energy improves very signifi-cantly even for the 2x2 case, which suggests that a large partof the inter-cluster correlations can be treated perturbatively.Third, even though the second-order critical point does notchange significantly with the cPT2 correction, the first-orderone shifts significantly to J / J ∼ .
62, compared to ∼ . L → ∞ and l → ∞ ). We have to remind the reader that itshould approach 0, because clusters of infinite size capture allthe energy at the cMF level. Similarly to the cMF analysis,the correlation energy was plotted with respect to the inverseof the cluster size, and at J / J = − . D. cCCSD results
We continue by considering the cCCSD results. All thecomputational parameters are the same as for the cMF andcPT2 results. We remind the reader that all the possible ex-cited states were used to compute the cCCSD energy. InFig. 10 we demonstrate the energy per site obtained for differ-ent values of J / J . We have also included, for comparison,the extrapolated results (to the thermodynamic limit L → ∞ and l → ∞ ) calculated from exact calculations computed inRef. 20. It is important to note that in the intermediate regime,the extrapolated results are not reliable, for reasons mentionedin Ref. 20 and that is why we do not show them. The energy ¶ S z =
0, 1250/11440 for S z = +
1, 1250/11440 for S z = − -0.75-0.7-0.65-0.6-0.55-0.5-0.45-0.4 0 0.2 0.4 0.6 0.8 1 E ne r g y pe r s i t e ( E / N ) J /J cMF-2x2cPT2-2x2cMF-2x4cPT2-2x4cMF-3x3cPT2-3x3cMF-2x6cPT2-2x6cMF-4x4cPT2-4x4 FIG. 9. Comparison of cMF and cPT2 for different cluster schemes isshown. cMF-2x2 stands for cluster mean-field with 2x2 tiles, etc andcPT2-2x2 stands for cPT2 with 2x2 tiles, etc. The circles denote cMFand the squares cPT2. It is worth noting that the cPT2 energies arevery close to each other, which suggests that for correlated methods,a minimal cluster size is sufficient. -0.7-0.65-0.6-0.55-0.5-0.45-0.4 0 0.2 0.4 0.6 0.8 1 E ne r g y pe r s i t e ( E / N ) J /J cMF-2x2cPT2-2x2cCCSD-2x2Reference (Ref. 20) FIG. 10. A comparison between cMF, cPT2 and cCCSD for the en-ergy per site. cCCSD improves the energy significantly compared tocPT2. The exact energy is taken from Ref. 20. Again, we need toemphasize, that in the intermediate regime, the extrapolated resultsare not reliable (a strange behavior is also depicted in the graph). improves compared to the cPT2 energy in the two antiferro-magnetic phases, but the correction is minimal for the para-magnetic phase. We think that this is due to the nature ofthe paramagnetic phase (see introduction for the debate of thenature of the paramagnetic phase), which suggests that thereexist correlations that cannot be captured by using just 2x2tiles and a low order coupled-cluster theory. One can arguethat the improvement over cPT2 is relatively small consider-ing the error of cPT2. This may imply that triples or quadru-ples are needed (discussed in the following section). We haveto emphasize, however, that at some point cCC should becomeexact.
V. DISCUSSION
In Sec. II, we have described the cluster mean-field ap-proach to treat strongly-correlated spin systems. A cMF stateis used as a variational ansatz for the ground state wavefunc-tion, which is guaranted by construction to provide better vari-ational estimates than HF when the size of the cluster is largerthan 1. Because of the simple cluster language, a RS-PTscheme can be easily adopted to account for the missing inter-cluster correlations. The results presented in Secs. IV B, IV C,and IV D provide evidence that a cluster-based approach can(semi)-quantitatively capture the physics of the ground stateof the 2D square J − J Heisenberg model. Due to the natureof the Heisenberg Hamiltonian, contributions to the second-order energy arise only from two-cluster spin interactions.From calculations performed in 1D (not in the manuscript),cMF, cPT2 and cCCSD are not as efficient for 2D as for 1D.A significant improvement to the ground state energy is ob-tained with cPT2 as well as cCCSD, which gives the best re-sults. We also notice that enlarging the size of the cluster inmean-field calculations is worse than performing a cPT2 or acCCSD calculation. The good quality of cPT2 and cCCSDresults suggest that the zero-th order Hamiltonian is suitablefor describing spin lattices, at least for the antiferromagneticregimes.In the rest of this section we discuss possible strategies thatcan improve the results presented in this article. The simpleststrategy, also discussed in Ref. 8, is to use the full Hilbertspace (not restricted to a given S z sector). This will give sig-nificantly more variational freedom to the cMF ansatz. Thisapproach, however, requires a Hilbert space of much largerdimension. As regards the cPT2 and the cCCSD methods, themost straightforward way to go beyond those is to use cPTnor cCCSDTQ, etc. To do so, we must truncate the numberof states used, because the computational cost will be pro-hibitive. This truncation scheme could be either based onthe local character of the clusters or can be found stochas-tically. Another possible route could be to exploit locality.For example, we can treat the interaction between nearest-neighbor clusters with cCC and with further clusters withcPT. One more advantage of the cluster-based approaches,is that even though we have used those approaches to studystrongly interacting systems, they may be used in other con-texts. More specifically, systems which can be effectively rep-resented in terms of weakly interacting fragments of strongly-correlated subsystems can be very efficiently described bycPT2 or cCCSD. Lastly, another route for correlating cMFwould be to write the ansatz as a linear combination of differ-ent cMFs of different tilings. This has been tried for dimersby Garcia-Bach and has yielded very promising results. VI. CONCLUSIONS
In this work, the optimization of the cluster mean-field statehas been carried out with the restriction that the cluster statehas well-defined S z quantum number. The restrictions are im-posed in order to preserve total S z in the full system and fa-cilitate the computation of the matrix elements. The clusterproduct state constitutes an eigenstate of a mean-field (zero-thorder) Hamiltonian, which allows us to go beyond mean-fieldin a perturbative and a coupled-cluster framework. We havepresented mean-field, second-order perturbative and coupled-cluster results of the ground state energies (and magnetiza-tions) of the square 2D Heisenberg model in the thermody-namic limit ( L → ∞ ). Also, we have presented a very accurateextrapolation scheme for thermodynamic limit ( L → ∞ and l → ∞ ) energies. In general, we observe that cPT2 energieswith small clusters are often better than cMF results with sig-nificantly larger ones, and the same applies to cCCSD. Over-all, the results of this work suggest that a cluster mean-fieldapproach can provide an excellent starting point and a path toa highly accurate, efficient description of strongly-correlatedsystems, while cPT2 and cCCSD provide an accurate quan-titative description. Several strategies to improve the mean-field description as well as correlated approaches built on topof it have been suggested. ACKNOWLEDGMENTS
This work was supported by the U.S. Department of En-ergy, Office of Basic Energy Sciences, Computational andTheoretical Chemistry Program under Award No. DE-FG02-09ER16053. CAJH is grateful for support from a start-uppackage at Wesleyan University. G.E.S. acknowledges sup-port as a Welch Foundation Chair (Grant No. C-0036).
DATA AVAILABILITY
The data that support the findings of this study are availablefrom the corresponding author upon reasonable request. S. Sharma, K. Sivalingam, F. Neese, and G. K.-L. Chan, “Low-energyspectrum of iron–sulfur clusters directly from many-particle quantum me-chanics,” Nature Chem , 927–933 (2014). D. Gatteschi, R. Sessoli, and J. Villain,
Molecular nanomagnets , Meso-scopic physics and nanotechnology No. 5 (Oxford University Press, 2006). J. Wu, T. G. Schmalz, and D. J. Klein, “An extended Heisenberg model forconjugated hydrocarbons,” The Journal of Chemical Physics , 9977–9982 (2002). J. L. Dye, “Electrides: From 1D Heisenberg Chains to 2D Pseudo-Metals,”Inorg. Chem. , 3816–3826 (1997). A. B. Van Oosten, R. Broer, and W. C. Nieuwpoort, “Heisenberg exchangein La2CuO4,” Int. J. Quantum Chem. , 241–243 (1995). I. W. Bulik, T. M. Henderson, and G. E. Scuseria, “Can Single-ReferenceCoupled Cluster Theory Describe Static Correlation?” J. Chem. TheoryComput. , 3171–3179 (2015). P. W. Anderson, “The Resonating Valence Bond State in La2CuO4 and Su-perconductivity,” Science , 1196–1198 (1987). C. A. Jiménez-Hoyos and G. E. Scuseria, “Cluster-based mean-field andperturbative description of strongly correlated fermion systems: Applica-tion to the one- and two-dimensional Hubbard model,” Phys. Rev. B ,085101 (2015). L. Isaev, G. Ortiz, and J. Dukelsky, “Hierarchical mean-field approach tothe J1–J2 heisenberg model on a square lattice,” Phys. Rev. B , 024409(2009). E. Schrödinger, “Quantisierung als Eigenwertproblem,” Ann. Phys. ,361–376 (1926). I. Shavitt and R. J. Bartlett,
Many-body methods in chemistry and physics:MBPT and coupled-cluster theory , Cambridge molecular science (Cam-bridge University Press, 2009). S. Li, “Block-correlated coupled cluster theory: The general formulationand its application to the antiferromagnetic Heisenberg model,” The Journalof Chemical Physics , 5017–5026 (2004). T. Fang, J. Shen, and S. Li, “Block correlated coupled cluster method with acomplete-active-space self-consistent-field reference function: The imple-mentation for low-lying excited states,” The Journal of Chemical Physics , 234106 (2008). J. Shen and S. Li, “Block correlated coupled cluster method with the com-plete active-space self-consistent-field reference function: Applications forlow-lying electronic excited states,” The Journal of Chemical Physics ,174101 (2009). V. Abraham and N. J. Mayhall, “Selected Configuration Interaction in a Ba-sis of Cluster State Tensor Products,” J. Chem. Theory Comput. , 6098–6113 (2020). J. I. Cirac and F. Verstraete, “Renormalization and tensor product states inspin chains and lattices,” J. Phys. A: Math. Theor. , 504004 (2009). H. Bethe, “Zur Theorie der Metalle: I. Eigenwerte und Eigenfunktionen derlinearen Atomkette,” Z. Physik , 205–226 (1931). E. Dagotto and A. Moreo, “Phase diagram of the frustrated spin-1/2 Heisen-berg antiferromagnet in 2 dimensions,” Phys. Rev. Lett. , 2148–2151(1989). H. J. Schulz and T. A. L. Ziman, “Finite-Size Scaling for the Two-Dimensional Frustrated Quantum Heisenberg Antiferromagnet,” Europhys.Lett. , 355–360 (1992). J. Richter and J. Schulenburg, “The spin-1/2 J1–J2 Heisenberg antiferro-magnet on the square lattice: Exact diagonalization for N=40 spins,” Eur.Phys. J. B , 117–124 (2010). L. Capriotti and S. Sorella, “Spontaneous Plaquette Dimerization in the J1– J2 Heisenberg Model,” Phys. Rev. Lett. , 3173–3176 (2000). M. Mambrini, A. Läuchli, D. Poilblanc, and F. Mila, “Plaquette valence-bond crystal in the frustrated Heisenberg quantum antiferromagnet on thesquare lattice,” Phys. Rev. B , 144422 (2006). H. J. Schulz, T. A. Ziman, and D. Poilblanc, “Magnetic Order and Disor-der in the Frustrated Quantum Heisenberg Antiferromagnet in Two Dimen-sions,” J. Phys. I France , 675–703 (1996). D. Schmalfuß, R. Darradi, J. Richter, J. Schulenburg, and D. Ihle, “Quan-tum J 1 - J 2 Antiferromagnet on a Stacked Square Lattice: Influence of theInterlayer Coupling on the Ground-State Magnetic Ordering,” Phys. Rev.Lett. , 157201 (2006). R. Darradi, O. Derzhko, R. Zinke, J. Schulenburg, S. E. Krüger, andJ. Richter, “Ground state phases of the spin-1/2 J 1 – J 2 Heisenberg anti-ferromagnet on the square lattice: A high-order coupled cluster treatment,”Phys. Rev. B , 214415 (2008). R. F. Bishop, D. J. J. Farnell, and J. B. Parkinson, “Phase transitions in thespin-half J 1 - J 2 model,” Phys. Rev. B , 6394–6402 (1998). H.-C. Jiang, H. Yao, and L. Balents, “Spin liquid ground state of the spin-1 2 square J 1 - J 2 Heisenberg model,” Phys. Rev. B , 024424 (2012). V. Murg, F. Verstraete, and J. I. Cirac, “Exploring frustrated spin systemsusing projected entangled pair states,” Phys. Rev. B , 195119 (2009). J.-F. Yu and Y.-J. Kao, “Spin- 1 2 J 1 - J 2 Heisenberg antiferromagnet on asquare lattice: A plaquette renormalized tensor network study,” Phys. Rev.B , 094407 (2012). L. Wang, D. Poilblanc, Z.-C. Gu, X.-G. Wen, and F. Verstraete, “Construct-ing a Gapless Spin-Liquid State for the Spin- 1 / 2 J 1 - J 2 HeisenbergModel on a Square Lattice,” Phys. Rev. Lett. , 037202 (2013). L. Capriotti, F. Becca, A. Parola, and S. Sorella, “Resonating Valence BondWave Functions for Strongly Frustrated Spin Systems,” Phys. Rev. Lett. ,097201 (2001). A. W. Sandvik, “Finite-size scaling of the ground-state parameters ofthe two-dimensional Heisenberg model,” Phys. Rev. B , 11678–11690(1997). M. P. Gelfand, “Series investigations of magnetically disordered groundstates in two-dimensional frustrated quantum antiferromagnets,” Phys. Rev.B , 8206–8213 (1990). M. E. Zhitomirsky and K. Ueda, “Valence-bond crystal phase of a frus-trated spin-1/2 square-lattice antiferromagnet,” Phys. Rev. B , 9007–9010(1996). K. Takano, Y. Kito, Y. ¯Ono, and K. Sano, “Nonlinear sigma Model Methodfor the J 1 - J 2 Heisenberg Model: Disordered Ground State with PlaquetteSymmetry,” Phys. Rev. Lett. , 197202 (2003). V. Lante and A. Parola, “Ising phase in the J1 - J2 Heisenberg model,” Phys.Rev. B , 094427 (2006). W.-J. Hu, F. Becca, A. Parola, and S. Sorella, “Direct evidence for a gaplessZ 2 spin liquid by frustrating Néel antiferromagnetism,” Phys. Rev. B ,060402 (2013). T. Li, F. Becca, W. Hu, and S. Sorella, “Gapped spin-liquid phase in theJ 1 - J 2 Heisenberg model by a bosonic resonating valence-bond ansatz,”Phys. Rev. B , 075111 (2012). Y.-Z. Ren, N.-H. Tong, and X.-C. Xie, “Cluster mean-field theory study ofJ -J Heisenberg model on a square lattice,” J. Phys.: Condens. Matter , 115601 (2014). E. R. Davidson, “The iterative calculation of a few of the lowest eigenvaluesand corresponding eigenvectors of large real-symmetric matrices,” Journalof Computational Physics , 87–94 (1975). G. E. Scuseria, T. J. Lee, and H. F. Schaefer, “Accelerating the convergenceof the coupled-cluster approach,” Chemical Physics Letters , 236–239(1986). M. Garcia-Bach, “Long-range spin-pairing order and spin defects in quan-tum spin- ladders,” Eur. Phys. J. B14