Numerical Renormalization Approach to Two-Dimensional Quantum Antiferromagnets with Valence-Bond-Solid Type Ground State

We study the ground-state properties of the two-dimensional quantum spin systems having the valence-bond-solid (VBS) type ground states. The ``product-of-tensors'' form of the ground-state wavefunction of the system is utilized to associate it with an equivalent classical lattice statistical model which can be treated by the transfer-matrix method. For diagonalization of the transfer matrix, we employ the product-wavefunction renormalization group method which is a variant of the density-matrix renormalization group method. We obtain the correlation length and the sublattice magnetization accurately. For the anisotropically ``deformed'' S=3/2 VBS model on the honeycomb lattice, we find that the correlation length as a function of the deformation parameter behaves very much alike as that in the S=3/2 VBS chain.


I. INTRODUCTION
There are many applications of the density matrix renormalization group (DMRG) [1,2] which was originally applied to the one-dimensional (1D) quantum systems [1,2]. Due to its remarkable success, the DMRG has now become one of the standard methods for studying 1D quantum models, two-dimensional (2D) classical models [3] and (1 + 1)-dimensional classical non-equilibrium models [4,5].
As was pointed out in Refs. [6,7], DMRG is a variational method under the matrix-product-form ansatz (MPFA) for trial wavefunctions whose usage dates back to the work of Kramers and Wannier [8][9][10]. This point of view leads to some non-trivial reformulations of the DMRG: the direct variational approach, [6,7,11,12] the product-wavefunction renormalization group (PWFRG), [13] the corner-transfer-matrix renormalization group (CTMRG). [9,10] Further, in this view, the success of the DMRG implies the unexpected accuracy of the MPFA wavefunctions.
Having seen the success of the DMRG for 1D quantum and 2D classical systems, it is natural and important to explore its higher-dimensional generalizations. For 2D quantum case, the ladder approach, [14][15][16] can be regarded as such, but it is essentially the one-dimensional algorithm. For 3D classical case, as an extension of the CTMRG, [9,10] the "corner-tensor" approach has recently been proposed, [17] which is the first among the truly higher-dimensional algorithms, although it may not be the only possibility.
In our view, the most promising one in generalizing the DMRG to higher dimensions seems to be the one which is based on a generalization of the MPFA. We should remark that, we have already known a one: a wavefunction which is a product of local tensors, typical example being the valence-bond-solid (VBS) state [18,19] for a class of 2D quantum antiferromagnets. It has been well known that the wavefunction of the 1D VBS state can be expressed as a product of 2 × 2 matrices, implying that the MPFA is exact for the system; it has also been known that the wavefunction of the 2D (or higher dimensional) VBS state can be expressed as a product of "tensors" (generalized objects of matrices). Let us simply call such form of wavefunction a tensors-product-form (TPF) wavefunction. A natural generalization of the MPFA is, then, the tensors-product-form ansatz (TPFA). Accordingly, we can formulate the TPFA-variational method which can be thought of a "higher-dimensional DMRG". The aim of the present article is to make the first step in this TPFA-variational approach to higher-dimensional quantum systems. We focus on the property of the VBStype state in 2D, and for this purpose, we develop a reliable numerical method at the same time.
The models we consider in the present article is the 2D quantum antiferromagnets which have the VBS-type ground states; [18][19][20] the TPFA is exact for these models. Let |VBS be the unnormalized ground-state vector whose wavefunction is exactly given by a product of local tensors. A main part of our problem is to evaluate the expectation value of a given observable O For 1D VBS-type models, due to the matrix-productform structure of |VBS , the RHS of (1.1) can be interpreted [18,19] as a thermal average in a 1D classical statistical-mechanical model; the transfer-matrix method allows us to evaluate the RHS of (1.1) exactly. For 2D VBS-type models, similar interpretation as a 2D classical statistical-mechanical problem is straightforward due to the TPF structure of |VBS . Based on this interpretation, Niggemann et al. made the classical Monte-Carlo study of the anisotropically generalized S = 3/2 VBS model on the honeycomb lattice. The approach we take in the present work is, in a sense, a direct generalization of the 1D case; we treat the associated 2D classical statistical-mechanical problem by the transfermatrix method. What is essential in our approach is that, for diagonalization of the transfer matrix, we em-ploy the DMRG (to be precise, PWFRG, in the present article) allowing us to make highly reliable, close-to-exact evaluation of (1.1). The use of the DMRG also has an important implication, in the light of the TPFA-variational method: 2D TPFA-variational calculation is reduced to 1D TPFA-variational method, namely the DMRG. Accordingly, we can, in principle, formulate the "nested" TPFA-variational approach where D-dimensional TPFAvariational calculation is reduced to (D − 1)-dimensional one, which in turn is reduced to (D − 2)-dimensional one, and so on. This paper is organized as follows. In Section II, we give explicit form of local tensors for the 2D VBS-type models. For actual calculations, we consider the S = 2 isotropic VBS model on the square lattice, and the "deformed" S = 3/2 VBS-type model on the honeycomb lattice proposed by Niggemann et al. In Section III, we explain our calculation method. We extend the PWFRG so that it can handle the asymmetric transfer matrices which we encounter in treating the 2D VBS-type models. In Section IV, calculated results are given. The last section (Section V) is devoted to summary and conclusion.

A. Two-Dimensional Valence-Bond-Solid Wavefunction and Local Tensors
Construction of the higher dimensional VBS state has already been made in the original paper. [18] We shall give a brief explanation of the TPF structure of the general VBS state.

B. Vertex-Model Interpretation
We can regard the tensor A(σ|{η}) as the Boltzmann weight of a statistical-mechanical model where both the spin variables {σ i } ("vertex spins") and the bond variables {η i α } are fluctuating variables, which we may call spin-vertex model. Then the unnormalized VBS wavefunction (2.6) is the partition function of the spin-vertex model with fixed vertex-spin configuration To obtain the quantum-mechanical expectation (1.1), we should calculate the norm VBS|VBS which can also be regarded as a partition function of a lattice model, namely, the vertex model. Assuming that the local tensors {A i (σ i |η i )} are real, we can write the vertex weight for the M -valent site as The sign factor in (2.6) can also be taken into account by introducing the modified vertex weights, for example, In this vertex model, the double index (η α , ξ α ) (α = 1, 2, . . .) plays the role of the bond variable in the ordinary vertex model. Since each component index takes two values ±1/2, the model is a 4-state vertex model. The quantum-mechanical expectation (1.1) can then be regarded as a statistical-mechanical average in the vertex model.

C. Square and Honeycomb Lattices
We treat the VBS-state-associated vertex model by the transfer matrix method. For ease of construction of the row-to-row transfer matrix, we restrict the analysis to the two cases in the present paper: S = 2 VBS model on the square lattice and the S = 3/2 VBS model on the honeycomb lattice. In the former, the row-to-row transfer matrix can be constructed in the conventional way. [21] In the latter, we can map the model on the honeycomb lattice to the one on the square lattice [20] as shown in Fig. 2. Hence either case can be treated as a vertex models on the square lattice.
fig-2 From the general formula (2.4), we can explicitly write down non-zero components of the local tensor. For S = 3/2 honeycomb-lattice VBS model, we have and for S = 2 square-lattice VBS model, we have A(2|1/2, 1/2, 1/2, 1/2) = 1, (2.14) Other non-zero components are obtained by using the permutational symmetry of A(σ|{η}) in {η} and the property (2.5). The VBS states constructed from these tensors are exact ground states of the Hamiltonian where the local nearest-neighbor (nn) -type pair Hamiltonian h ij is the projection operator into the local spin-2S space, P 2S ij , (S = 3/2 for honeycomb lattice, and S = 2 for square lattice): [18,19] P 2S ij |VBS = 0, for all nn site pair ij. In the original VBS model, the local Hamiltonian P 2S ij is a function of the inner product s i · s j , hence, is isotropic. [18,19] Anisotropic generalization of the VBS models and their associate ground-state vector has been known. [22,23] Similar generalization has also been known for the honeycomb VBS model. [20] The generalized honeycomb VBS model contains a "deformation" parameter which corresponds to the xxz-type anisotropy of the Hamiltonian. The associated wavefunction is also of the TPFA form with a deformed local tensor A(σ|{η}) a where a is the deformation parameter. Explicit form of A(σ|{η}) a is simply given by that of the isotropic case (2.13) with substitution By Ψ VBS(a) (σ 1 , σ 2 , . . . , σ N ), we denote the TPFA wavefunction made from {A i (σ i |{η i }) a }, and by |VBS(a) , corresponding state vector. Since |VBS(−a) is essentially equivalent to |VBS(a) (via a S z -diagonal unitary transformation), we treat a > 0 in the present article.
(We have confirmed that calculated results for −a are identical to those for a.) As has been shown in Ref. [20] there actually exists Hamiltonian with nn-type interactions whose ground state is precisely |VBS(a) (× normalization factor).

III. CALCULATION METHOD
In the previous section, we explained how to translate physical quantities from the quantum-mechanical frame (Eq. (1.1) ) into the classical statistical-mechanical one (Eq. (2.11) and/or Eq. (2.12)).
In this section, to deal with the classical statisticalmechanical frame, we introduce a transfer matrix. So now we focus on a way of renormalization of an asymmetric transfer matrix by using the "classical PWFRG" [13]. Our transfer matrix is composed of 4-state vertex weights.
We denote this 4-state vertex weight by E p q (a, b), which is depicted by Fig. 3. We regard this Let us explain how to obtain the expectation value of any one-point operator A i at site i (by the PWFRG): where ψ represents the ground state of the VBS model (for instance, |ψ = |VBS(a) ). We represent the numerator of Eq. (3.1) by transfer matrices T and T A , whose definitions will be shown later (Eq. (3.3) and Eq. (3.4)): where N is a linear size. T and T A in this equation are defined in the following (for short , we represent where 2M is another linear size (along this direction, we do not impose a periodic boundary condition). And E A is composed of a quantum-mechanical operator A which is sandwiched between classical weights (Eq. (2.4) and/or Eq. (2.9)). Note that E ≡ E id (id: a local identity operator). Using these notations, where R is composed of column right eigenvectors, and D is diagonal matrix whose diagonal elements {λ j } (these λ 1 , λ 2 , . . . are in descending order of the absolute value) is the set of eigenvalues of the transfer matrix T and L is consist of row left eigenvectors. Here we should keep in mind that T is, generally, an asymmetric matrix, upon whose left and right eigenvectors L and R we impose the normalization condition LR = 1.
where the subscript "1" represents the eigenstate with the largest eigenvalue of T in absolute value and we denote its left eigenvector and right eigenvector by l| and |r respectively. The denominator of Eq. (3.1) (the "partition function") is obtained by To get l|, T A , |r and λ 1 in this equation, we use the PWFRG method (infinite method). As we stated earlier, this transfer matrix is asymmetric. So we decide how to renormalize the asymmetric transfer matrix using the PWFRG. This is because there are several ways to treat an asymmetric transfer matrix. In this paper, we treat it in the following way which is different from our previous treatment. [4] In our previous treatment, we adopted asymmetric density matrices in the DMRG. Unlike this, in this paper, our treatment is equivalent to treating symmetric density matrices in the PWFRG (namely, we use the singular value decomposition (SVD) for left and right eigenvectors of T ). It is easier to perform the renormalization with the SVD than renormalization with diagonalization of asymmetric matrices, because the former is free from difficulties of manipulating complex numbers. In compensation for this difficulty, we only have to introduce identity operators as will be shown later (see Eq. (3.11)). These identity operators is very important for the correct renormalization of an asymmetric transfer matrix.
First of all, we present the outline of our prescription in this paper before explaining details.
Our prescription in this paper is: P-1 Multiply left eigenvector and right eigenvector obtained in the previous iteration by transfer matrix (and by an identity operator) several times to get the improved left eigenvector and right eigenvector.
P-2 Perform SVD for the improved left and right eigenvectors, respectively to get "projectors".
P-3 Select important states and renormalize transfer matrices, identity operators and one-point operators. P-4 Get improved "projectors" using the "recursion relation" [13] and, from these "projectors", construct new left and right eigenvectors.
P-5 Evaluate one-point functions such as A i by using Eq. (3.9). Check the convergence of both left and right wavefunctions and the value of the one-point function. If convergence is true then we calculate 2-point function or evaluate λ 2 to calculate the correlation length (see section III B).
The followings are the details in the above prescription: (Hereafter, we use Greek indices to represent the block state.) First, a new notation is introduced. We renormalize bare T into renormalized transfer matrix T where, in this equation, T L and T R are renormalized as Eq. (3.15). Eq. (3.10) is expressed according to Nishino's diagram [3] in Fig. 4. fig-4 Second, we give caveats in our prescription below: As for P-1: To concretely explain what is mentioned in P-1, we take the case of improving right eigenvector |r imp . It is important that whenever we multiply by a transfer matrix, we simultaneously must multiply by left and right renormalized identity operators id L and id R . Assuming that k is some integer, where (id L ) T is the transpose of id L etc. More concretely, rewriting this equation in its component (we treat the case k = 1 for notation simplicity), where [id L ] α µ and [id R ] β ν are left and right renormalized identity operators (see Eq. (3.16)).
As for P-2: SVD is performed as: And the left identity operator: id L is renormalized as At the very first iteration (when we treat a 4-site chain of vertices T = EEEE), we set This is an initial condition for the recursion relation Eq. (3.16). Next, we prepare an object TA L . This is for Eq. (3.21). In the same way as Eq. (3.10), we renormalize bare T A into T A as follows. (3.18) where TA L is made of As for P-4: For example, U R 's recursion relation is: (3.20) As for P-5: We show how to obtain the expectation value of any onepoint operator A i in the next way.

B. Correlation Length
In what follows, we explain two methods of calculating correlation lengths. For concreteness, we explain the case of the spin-3/2 honeycomb lattice model. It is also easy to deal with the spin-2 square lattice model.

Calculus by observation of dumping correlators
We set j − i ≡ r(> 0) and assume N ≫ r ≫ 1 in the sequel. And we locate each operator A in the position "b" in Fig. 2.
So, we have like Eq.(3.7) From this equation, it follows that we measure correlation length in the direction depicted in Fig. 5. fig-5 In off-critical region, we expect We bear in mind that in Eq. (3.22), ζ and τ includes ζ = 1 and τ = 1 (let us recall that the subscript "1" represents the eigenstate with the largest eigenvalue of T in absolute value). So, in the region A i = 0, we easily evaluate correlation length from Eq. (3.22) and Eq. (3.23) as

Evaluation by using the power method
Another method by which correlation lengths are evaluated is the following. By making use of Eq. (3.5), we can further transform Eq.(3.22) into (3.25) where a constantC A is 3.26) and the definition of ξ pm in this equation is 27) and l 2 | and |r 2 are the left and right eigenvectors of the eigenvalue λ 2 , respectively. From Eqs. (3.25), (3.26) and (3.27), we evaluate λ 2 to obtain the correlation length ξ −1 pm . We use the power method to evaluate λ 2 that is, we perform the power method in the space which is orthogonal to the ground state vectors l| and |r , as follows.
First, we prepare the initial vectors l(1) | and |r (1) . And then we get the improved left eigenvector l(n) | and the improved right eigenvector |r (n) . which are orthogonal to l| and |r . Namely, and |r (n) ≡ T |r (n−1) − |r l|T |r (n−1) l|r . (3.29) We should note that our l(n) | and |r (n) do not satisfy Eq.(3.6). So, slightly modifying Eq. (3.8), we calculate 2nd largest eigenvalue λ 2 of T in the following way.
It should be noted that we save memory resource in a computer by using the power method (That is we need only O(m 2 ) memory, not O(m 4 ), where m is a number of retained base in the PWFRG). So, if necessary, we can perform the above calculation with a larger number m.

IV. RESULTS
A. Isotropic VBS models on honeycomb and square lattices As has been mentioned in Sec. II, both the S = 3/2 VBS model on the honeycomb and the S = 2 VBS model on the square lattice can be mapped to classical vertex models on the square lattice. We made PWFRG calculation for the sublattice magnetization and the correlation length in a way as described in Sec. III.
We have confirmed that sublattice magnetization is zero (Néel order is absent) in both models; both models are in the disordered phase, in agreement with previous studies. [18,19,24] As for the correlation length ξ, we obtain We should note that, in the honeycomb case, the length unit for the above value of ξ −1 honeycomb is the lattice spacing of the mapped square lattice. Unlike the isotropic case, the anisotropic (xxz-like) generalization of honeycomb VBS model can be either in the disordered phase or the ordered (Néel) phase, depending on the "deformation parameter" a. [20] The PWFRG allows us to study the behavior of the system in much more detail, which we present in this subsection. In the actual calculation, we have observed that the PWFRG calculations are rapidly convergent in increasing m (number of retained bases), except for the cases very near the transition point. Consequently, for most values of a, very small m, say m = 12, is sufficient for the accuracy required in the present study.

Staggered magnetization
In Fig. 6, we show the PWFRG result of the sublattice magnetization (staggered magnetization) M st as a function of the parameter a. fig-6 The Néel order which exists in the large-a region disappear at a = a c ∼ = 2.54, in agreement with the Monte-Carlo result. [20] It may well be expected that the Isinglike anisotropy make the system be in the Ising universality class as regards this anisotropy-induced phase transition. To check this, we plot (M st ) 8 versus a as shown in

Correlation length
In the transfer-matrix method, the correlation length is usually calculated from the ratio between the largest and the next-to-largest eigenvalues of the transfer matrix. This definition of the correlation length is, in our notation, ξ pm (see Eq. (3.27)), whose PWFRG result is shown in Fig.8. fig-8 The linear behavior of ξ −1 pm near a c is consistent with the expectation that the system is in the Ising-model universality-class, namely, In addition to the critical behavior at a c , a notable feature of ξ pm is that it exhibits cusp singularity at the isotropic point a = a iso = √ 3. To clarify the origin of this cusp, we calculated longitudinal correlation length ξ l (a) and the transverse correlation length ξ t (a) from the correlation functions S z i S z j and S x i S x j , as shown in Fig.9  fig-9 Totally different a-dependence between ξ −1 l (a) and ξ −1 t (a) is seen; the isotropic point a = √ 3 is the crossing point of the two curves. Since ξ pm is determined by the correlator with the largest correlation length (or, equivalently, smallest inverse correlation length), we have which explains the cusp at a = a iso = √ 3 in Fig.8. It is interesting to compare the correlation lengths in the 2D deformed VBS model with those in the 1D counterpart (deformed S = 3/2 chain ) [23], which we denote ξ which also show the crossover behavior at the isotropic point a = √ 3. It should be noted that magnitudes of these 1D correlation lengths are almost the same as those in the 2D case in a wide range of parameter values. Major differences are the absence of the critical point a c and the non-vanishing behavior of ξ (1D) t −1 in the a → 0 limit (see Fig. 10). fig-10 In Ref. [20], the "asymptotic equivalence" between the deformed honeycomb VBS model in the large-a region and the free-fermion model is pointed out, which gives fairly accurate estimation of the critical point a c . As a further check of this equivalence, we compare the correlation length ξ l with that of the free-fermion model which we denote by ξ ff = ξ ff (a). The method given in Refs. [25][26][27] allows us to obtain ξ ff (a) from the solution of Q(iω * ) = 0, where Q(φ) is given by Eq. (24) in Ref. [20] (ξ −1 ff = ω * in the disordered phase and ξ −1 ff = 2ω * in the ordered phase). In Fig.11, we see a remarkable agreement in a unexpectedly wide parameter range including the small-a region where the "asymptotic equivalence" may not hold any longer. fig-11 The cusp behavior of the correlation length implies that the isotropic VBS model may be a "singular" point in the parameter space. We should point out that a similar cusp-like behavior of the correlation length at the VBS point has been known in the β − ξ −1 curve [28] and the θ − ξ curve [29] of the S = 1 bilinear-biquadratic chain, where −β and θ relate to the relative amplitude of the biquadratic exchange term. To explore the relation between these cusps in the different parameter spaces is an interesting problem, which may help us to characterize the disordered phase in the 2D VBS-like models, in more detail: like the string order parameter [30] which characterizes the Haldane phase in the 1D case.

V. SUMMARY AND CONCLUSION
In this paper, we have presented a new way of application of the product wavefunction renormalization group (PWFRG) which is a variant of the density-matrix renormalization group (DMRG), to two-dimensional (2D) valence-bond-solid (VBS) type models.
The exact tensors-product-form (TPF) structure of the groundstate wavefunctions have allowed us to map the systems into two-dimensional classical statistical-mechanical models, namely, the vertex models which can be treated by the transfer-matrix method. We extend the PWFRG method so that it can handle the asymmetric transfer matrix associated with the VBS-type models.
For the isotropic 2D VBS models on the honeycomb lattice (S = 3/2) and the square lattice (S = 2), we have confirmed that they are in the disordered phase without the Néel order. We have obtained accurate values of correlation length for both models. For the anisotropically generalized (or deformed) VBS-type model on the honeycomb lattice [20], we have obtained detailed parameterdependence of the Néel order and the correlation length. We have confirmed the anisotropy induced phase transition at a critical value of the deformation parameter, and that this second-order phase transition belongs to the 2D Ising-model universality class. In the disordered phase, we have found a cusp-like behavior of the correlation length. We have explained this behavior as the crossover between the longitudinal and the transverse correlation length. Further, the PWFRG-calculated longitudinal correlation length well agrees with the correlation length of the free-fermion model in a fairly wide range of the deformation parameter, which far exceeds the prediction of the "asymptotic equivalence". [20] As has been mentioned in the Introduction, the present study is also aimed at higher-dimensional generalization of the DMRG-type numerical renormalization approach, based on the TPF-ansatz on the trial wavefunctions. From this point of view, what we have done in the present paper is calculation of physical quantities under given TPF-wavefunction with fixed local tensor. What we should do next is to vary the local tensor to make the variational calculation. A "direct numerical" variational calculation can actually be done, where we sweep the variational parameter and numerically find the minimum of the energy-expectation. The PWFRG method is also useful for reliable determination of the optimal variational value of the parameter. Studies along this line is now undertaken, whose details will be published elsewhere. The mapping of the model on the honeycomb lattice to the model on the square lattice. Circles (•), which express quantum spins (S = 3/2), are connected by solid lines that form the honeycomb lattice. The small letters "a" and "b" indicate sites of the two sublattices, respectively (see Sec. III B 1).
FIG. 3. E p q (a, b) is expressed according to Nishino's diagram [3]. Each circle (•) represents the mapped classical site, whose states are expressed by "p", "q", "a" and "b". These state variables take 4 states of the double index (η i α , ξ i α ).