Numerical renormalization approach to two-dimensional quantum antiferromagnets with valence-bond-solid type ground state

We study the ground-state properties of 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 as that in the S = 3/2 VBS chain.


Introduction
There are many applications of the density matrix renormalization group (DMRG) [1,2] which was originally applied to 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 nonequilibrium models [4,5].
As was pointed out in [6,7], DMRG is a variational method under the matrix-productform 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]. 7.2 Furthermore, in this view, the success of the DMRG implies the unexpected accuracy of the MPFA wavefunctions.
Having observed the success of the DMRG for 1D quantum and 2D classical systems, it is natural and important to explore its higher-dimensional generalizations. For the 2D quantum case, the ladder approach [14,15,16], can be regarded as such, but it is essentially the onedimensional algorithm. For the 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 algorithm in generalizing the DMRG to higher dimensions seems to be that which is based on a generalization of the MPFA. We should remark that we are already familiar with a wavefunction which is a product of local tensors, a typical example being the valence-bond-solid (VBS) state [18,19] for a class of 2D quantum antiferromagnets. It is 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 is also 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 a wavefunction a tensor-product-form (TPF) wavefunction. A natural generalization of the MPFA is, then, the tensor-product-form ansatz (TPFA). Accordingly, we can formulate the TPFA-variational method which can be thought of a higher-dimensional DMRG. The aim of this paper is to take the first step in this TPFA-variational approach to higherdimensional quantum systems. We focus on the property of the VBS-type state in 2D and, for this purpose, we develop a reliable numerical method at the same time.
The models we consider in this paper are 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-product-form structure of |VBS , the RHS of (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) exactly. For 2D VBS-type models, a 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 transfer-matrix method. What is essential in our approach is that, for diagonalization of the transfer matrix, we employ the DMRG (to be precise, PWFRG, in this paper), allowing us to make highly reliable, close-to-exact evaluation of (1). The use of the DMRG also has an important implication, in the light of the TPFA-variational method: the 2D TPFA-variational calculation is reduced to a 1D TPFA-variational method, namely the DMRG. Accordingly, we can, in principle, formulate the 'nested' TPFA-variational approach where a D-dimensional TPFA-variational calculation is reduced to a (D − 1)-dimensional one, which in turn is reduced to a (D − 2)-dimensional one, and so on. This paper is organized as follows. In section 2 we give an 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 3 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 4, calculated results are given. Section 5 is devoted to a summary and conclusion.

Two-dimensional valence-bond-solid wavefunction and local tensors
Construction of the higher dimensional VBS state has already been made in an earlier paper [18]. We shall give a brief explanation of the TPF structure of the general VBS state.
Consider a spin-S operator situated at a site i on a lattice which may not necessarily be a regular one. Let {|σ } (σ = −S, −S + 1, . . . , S) be the corresponding s z -diagonal basis set. We assume that the site i is 'M -valent' with M = 2S. We then prepare M spin-1/2 operators, and regard the spin-S operator as that constructed from these component operators. Let |η α (η α = ±1/2, α = 1, . . . , M) be the s z -diagonal base of the αth component spin, and denote |{η} ≡ |η 1 , η 2 , . . . , η M ≡ |η 1 ⊗ |η 2 ⊗ · · · ⊗ |η M . ( Consider the symmetrized state where Sym η stands for the symmetrization with respect to the indices {η α }, and N η is a factor introduced for proper normalization, ({η})|({η}) = 1. Clearly, we have We then define a local tensor A(σ|{η})(= A(σ|η 1 , η 2 , . . . , η M )) by which is the building block of the VBS state (see figure 1). Note that, by definition, the tensor A(σ|η 1 , η 2 , . . . , η M ) is symmetric with respect to the indices {η α }. Furthermore, we adopt the phases of the vectors {|σ } so that the tensor has the property 7.4 Let A i (σ i |{η i }) be the local tensor associated with the site i. The VBS wavefunction Ψ VBS (σ 1 , σ 2 , . . . , σ N ) is given by where we have introduced the sign factor S ij (η i , η j ) for each connected site pair < i, j >, corresponding to formation of the valence bond [18,19]. For the valence bond formed from the αth bond at the site i and βth bond at the site j, the sign factor is explicitly given by where The sign factor can be absorbed into the definition of modified local tensor as We can introduce modified tensors where two or more indices (at αth, βth, . . ., positions) are replaced by those with bars (Ā(σ|η 1 , . . . ,η α , . . . ,η β , . . .), etc.). In terms of the modified tensors, we can express the VBS wavefunction as a product of these tensors summed over the bondvariables {η} in a form similar to (7), with the sign factor S i,j (η i , η j ) simply replaced by the Kronecker delta where we have assumed that the αth bond at site i and βth bond at site j are to be connected, and that the tensor at site j is of the modified type at the βth position.

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 the spin-vertex model. Then the unnormalized VBS wavefunction (7) is the partition function of the spin-vertex model with fixed vertex-spin configuration To obtain the quantum mechanical expectation (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 (7) can also be taken into account by introducing the modified vertex weights, for example,W  (1) can then be regarded as a statistical-mechanical average in the vertex model.

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 figure 2. Hence either case can be treated as a vertex models on the square lattice. From the general formula (5), we can explicitly write down non-zero components of the local tensor. For the S = 3/2 honeycomb-lattice VBS model, we have and for the S = 2 square-lattice VBS model, we have Other non-zero components are obtained by using the permutational symmetry of A(σ|{η}) in {η} and the property (6). The VBS states constructed from these tensors are exact ground states of the Hamiltonian where the local nearest-neighbour (nn)-type pair Hamiltonian h ij is a q p b Figure 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 α ).

'Deformed' VBS-type model on the honeycomb lattice
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. The explicit form of A(σ|{η}) a is simply given by that of the isotropic case (14) with the substitution By Ψ VBS(a) (σ 1 , σ 2 , . . . , σ N ), we denote the TPFA wavefunction made from {A i (σ i |{η i }) a }, and by |VBS(a) , the corresponding state vector. Since |VBS(−a) is essentially equivalent to |VBS(a) (via a S z -diagonal unitary transformation), we treat a > 0 in this paper. (We have confirmed that calculated results for −a are identical to those for a.) As has been shown in [20], there actually exists a Hamiltonian with nn-type interactions whose ground state is precisely |VBS(a) (× normalization factor).

Calculation method
In the previous section, we explained how to translate physical quantities from the quantum mechanical frame (equation (1)) into the classical statistical-mechanical one (equation (12) and/or equation (13)). In this section, to deal with the classical statistical-mechanical frame, we introduce a transfer matrix. So we now 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 in figure 3. We regard this E p q (a, b) as a matrix element of E p q . 7.7

Magnetization
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 equation (18) by transfer matrices T and T A , whose definitions will be shown later (equation (20) and equation (21)): where N is a linear size. T and T A in this equation are defined in the following (for short , we represent {. . . p −1 , p 0 , p 1 , . . .} by {p}).
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 (equation (5) and/or equation (10)). Note that E ≡ E id (id: a local identity operator). Using these notations, Furthermore T is decomposed into where R is composed of column right eigenvectors, and D is a diagonal matrix whose diagonal elements {λ j } (these λ 1 , λ 2 , . . . are in descending order of the absolute value) are the set of eigenvalues of the transfer matrix T and L consists 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 Using equations (19), (22) and (23), 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 equation (18) (the 'partition function') is obtained by From equations (18), (24) and (25), we get the one-point function 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. In contrast, 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 the difficulties of manipulating complex numbers. In compensation for this difficulty, we only have to introduce identity operators as will be shown later (see equation (28)). These identity operators are 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 the details.

P1 Multiply the left and right eigenvectors obtained in the previous iteration by a transfer matrix
(and by an identity operator) several times to get the improved left eigenvector and right eigenvector. P2 Perform SVD for the improved left and right eigenvectors, respectively, to get 'projectors'. P3 Select important states and renormalize transfer matrices, identity operators and one-point operators. P4 Get improved projectors using the 'recursion relation' [13] and, from these projectors, construct new left and right eigenvectors. P5 Evaluate one-point functions such as A i by using equation (26). Check the convergence of both left and right wavefunctions and the value of the one-point function. If convergence is true then we calculate the 2-point function or evaluate λ 2 to calculate the correlation length (see section 3.2). P6 Return to step P1.
We now present the details of the above prescription. (Hereafter, we use Greek indices to represent the block state.) First, we introduce a new notation. We renormalize T into the renormalized transfer matrix T where T L and T R are renormalized as equation (32). Equation (27) is expressed according to Nishino's diagram [3] in figure 4.  Second, we give caveats in our prescription below.

P1
To concretely explain the step P1, we take the case of the 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 components (we treat the case k = 1 for simplicity of notation), where [id L ] α µ and [id R ] β ν are left and right renormalized identity operators (see equation (33)). P2 SVD is performed as: and where {[σ L ] µ } µ and {[σ R ] µ } µ are sets of singular values, and U L , V L , U R and V R are orthogonal matrices.

P3
We present some examples. The new renormalized transfer matrix: T (new) L is given by 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

7.10
This is an initial condition for the recursion relation equation (33). Next, we prepare an object TA L . This is for equation (38). In the same way as equation (27), we renormalize T A into T A as follows.
where TA L is composed of

P4
As an example, U R 's recursion relation is:

P5
We show how to obtain the expectation value of any one-point operator A i as follows:

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 following. We also locate each operator A in position b in figure 2. So, as in equation (24), we have From this equation it follows that we can measure correlation length in the direction depicted in figure 5.
In the off-critical region, we expect We bear in mind that in equation (39), ζ and τ includes ζ = 1 and τ = 1 (we recall that the subscript 1 represents the eigenstate with the largest eigenvalue of T in absolute value). So, in

Evaluation by using the power method
Another method by which correlation lengths are evaluated is the following. By making use of equation (22), we can further transform equation (39) into where a constantC A is and the definition of ξ pm in this equation is and l 2 | and |r 2 are the left and right eigenvectors of the eigenvalue λ 2 , respectively. From equations (42), (43) and (44), 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) . We then obtain 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 .
We should note that our l (n) | and |r (n) do not satisfy equation (23). So, slightly modifying equation (25), we calculate the second largest eigenvalue λ 2 of T as follows: where λ 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 the number of retained bases in the PWFRG). So, if necessary, we can perform the above calculation with a larger number m.

Isotropic VBS models on honeycomb and square lattices
As has been mentioned in section 2, 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 section 3.
We have confirmed that sublattice magnetization is zero (the Néel order is absent) in both models; both models are in the disordered phase, in agreement with previous studies [18,19,24]. 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.

S = 3/2 deformed VBS model on the honeycomb lattice
Unlike the isotropic case, the anisotropic (xxz-like) generalization of the 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 behaviour 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 with 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. figure 6, we show the PWFRG result of the sublattice magnetization (staggered magnetization) M st as a function of the parameter a.

Staggered magnetization In
The Néel order which exists in the large-a region disappears at a = a c ∼ = 2.54, in agreement with the Monte Carlo result. [20] It may well be expected that the Ising-like anisotropy causes the system to be in the Ising universality class as regards this anisotropy-induced phase transition.

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 equation (44)), whose PWFRG result is shown in figure 8. The linear behaviour 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 behaviour at a c , a notable feature of ξ pm is that it exhibits a 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 figure 9 Totally different a-dependence between ξ −1 l (a) and ξ −1 t (a) is observed; 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 figure 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 behaviour at the isotropic point a = √ 3. It should be noted that the magnitudes of these 1D correlation lengths are almost the same as those in the 2D case for a wide range of parameter values. Major differences are the absence of the critical point a c and the non-vanishing behaviour of ξ (1D) t −1 in the a → 0 limit (see figure 10). In [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 a 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 [25,26,27] allows us to obtain ξ ff (a) from the solution of Q(iω * ) = 0, where Q(φ) is given by equation (24) in [20] (ξ −1 ff = ω * in the disordered phase and ξ −1 ff = 2ω * in the ordered phase). In figure 11, we see a remarkable agreement in a unexpectedly wide parameter range including the small-a region where the 'asymptotic equivalence' may no longer hold.
The cusp behaviour 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 behaviour 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, such as the string order parameter ξ D Figure 11. Inverse correlation lengths of the free-fermion model (solid curves) and honeycomb VBS model (•: ξ −1 l (a), +: ξ −1 pm (a)).
[30] which characterizes the Haldane phase in the 1D case.

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 tensorproduct-form (TPF) structure of the ground-state wavefunctions has allowed us to map the systems into two-dimensional classical statistical-mechanical models, namely 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 parameter-dependence 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 behaviour of the correlation length. We have explained this behaviour as the crossover between the longitudinal and the transverse correlation length. Furthermore, the PWFRG-calculated longitudinal correlation length agrees well 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 higherdimensional 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 to calculate physical quantities under given TPF-wavefunction with a fixed local tensor. The next step 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 7.17 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 these lines are now being undertaken, the results of which will be published elsewhere.