# Spontaneous symmetry breaking in a generalized orbital compass model

###### Abstract

We introduce a generalized two-dimensional orbital compass model, which
interpolates continuously from the classical Ising model to the orbital
compass model with frustrated quantum interactions, and investigate
it using the multiscale entanglement renormalization ansatz (MERA).
The results demonstrate that increasing frustration of exchange
interactions triggers a second order quantum phase transition to a
degenerate symmetry broken state which minimizes one of the
interactions in the orbital compass model. Using boson expansion
within the spin-wave theory we unravel the physical mechanism of
the symmetry breaking transition as promoted by weak quantum
fluctuations and explain why this transition occurs only surprisingly
close to the maximally frustrated interactions of the orbital compass
model. The spin waves remain gapful at the critical point, and both
the boson expansion and MERA do not find any algebraically decaying
spin-spin correlations in the critical ground state.

Published in: Physical Review B 82, 104416 (2010).

###### pacs:

75.10.Jm, 03.65.Ud, 03.67.Hk, 64.70.Tg## I Introduction

The orbital compass model (OCM) is physically motivated by the orbital interactions which arise for strongly correlated electrons in transition metal oxides with partly filled degenerate orbitals and lead to rich and still poorly understood quantum models. In these systems the orbital degrees of freedom play a crucial role in determining collective states such as coexisting magnetic and orbital order, as for instance in the colossal magnetoresistance manganitesTok00 or in the vanadate perovskites.Hor08 The orbital interactions are typically intrinsically frustrated and may strongly enhance quantum fluctuations, leading to disordered states. Fei97 While realistic orbital interactions are somewhat complex, a paradigm of intrinsic frustration is best realized in the OCM,Kho03 ; Nus04 ; Mis04 ; Dor05 ; Dou05 with the pseudospin couplings intertwined with the orientation of interacting bonds. Its two-dimensional (2D) version on a honeycomb lattice, Jac09 realized in layered iron oxides,Nag07 is equivalent to the Kitaev model.Kit06

Although conceptually quite simple, the OCM has an interdisciplinary character as it plays an important role in a variety of contexts beyond the correlated transition metal oxides, such as: (i) the implementation of protected qubits for quantum computations in Josephson lattice arrays,Dou05 (ii) topological quantum order,Nus09 or (iii) polar molecules in optical lattices and systems of trapped ions.Mil07 Numerical studies Dor05 suggested that when anisotropic interactions are varied through the isotropic point of the 2D OCM, the ground state is not an orbital liquid type but instead a first order quantum phase transition (QPT) occurs between two different types of Ising-type order dictated by one or the other interaction. Recently the existence of this transition, similar to the one which occurs in the exact solution of the one-dimensional OCM,Brz07 was confirmed using projected entangled-pair state algorithm.Oru09 This implies that the symmetry is spontaneously broken at the compass point, and the spin order follows one of the two equivalent frustrated interactions.

Knowing that the ground states of the 2D Ising model and the 2D OCM are quite different, we introduce a generalized OCM which interpolates between these two limiting cases. Using this model we will investigate: (i) the physical consequences of gradually increasing frustration in a 2D system, (ii) where a QPT occurs from the Ising ground state to the degenerate ground state of the OCM, and, finally, (iii) the order and the physical mechanism of this QPT. As increasing frustration of the orbital interactions introduces entangled states, the present problem provides a unique opportunity to use the recently developed multiscale entanglement renormalization ansatzVid07 ; Vid08 (MERA) in order to find reliable answers to the above questions. As we show below, the QPT in the generalized OCM occurs only surprisingly close to the maximally frustrated interactions in the OCM. We also explain the physical origin of this behavior using an analytic approach based on the spin-wave theory.

Quantum many-body systems exhibit several interesting collective phenomena. Recent progress in developing efficient numerical methods to study quantum systems on a lattice is remarkable and allowed to investigate complex many-body phenomena, including QPTs.Sac99 An important step here was the discovery of density matrix renormalization group, Whi92 a very powerful numerical method that can be applied to one-dimensional strongly correlated fermionic and bosonic systems.Hal06 This idea played a fundamental role in developing entanglement renormalization Vid07 ; Vid08 to study quantum spin systems on a 2D lattice. Crucial in this approach is the removal of short-range entanglement by unitary transformations called disentanglers. It generates a real-space renormalization group transformation implemented in the MERA, which was recently successfully employed to investigate several quantum spin models,Cin08 ; EveIsing ; Eve09 ; Gio ; Dav08 and interacting fermion systems.Eis09 ; Vid10 So far, the very promising MERA has been applied inter alia to the 2D quantum Ising model,Cin08 ; EveIsing and to the Heisenberg model on a kagome lattice,Kagom but other possible applications and the optimal geometries for performing sequentially disentanglement and isometry transformation were also discussed.Eve09

The paper is organized as follows. In Sec. II we introduce the generalized OCM and state the problem of the existence and nature of the QPT. Next we present the MERA algorithm in Sec. III used to investigate frustrated interactions in the OCM. Numerical results obtained using the MERA are presented in Sec. IV. In order to explain the physical mechanism of the QPT found in the OCM we performed the boson expansion within the spin-wave theory, as described in Sec. V. The details of this expansion are presented in the Appendix. The paper is summarized in Sec. VI, where the main conclusions of the present work are also given.

## Ii Generalized compass model

In this paper, we investigate the nature and position of the QPT when the OCM point is approached in a different way from that studied before,Dor05 ; Oru09 namely when frustration of interactions along two nonequivalent directions gradually increases. Therefore, we introduce a 2D generalized OCM with ferro-like interactions notefo on a square lattice in plane (we assume the exchange constant ),

(1) |

The interactions occur between nearest neighbors and are balanced along both lattice directions and . Here labels lattice sites, with () increasing along () axes, and are linear combinations of Pauli matrices describing interactions for spins:

(2) | |||||

(3) |

The interactions in Eq. (1) include the classical Ising model at for operators and become gradually more frustrated with increasing angle — they interpolate between the Ising model (at ) and the isotropic OCM (at ), see Fig. 1. The latter case is equivalent to the 2D OCM with standard interactions and along the and directions Kho03 ; Nus04 ; Mis04 ; Dou05 ; Dor05 by a straightforward unitary transformation. The model (1) includes also as a special case the 2D orbital model for electrons at ,notefo describing, for instance, the orbital part of the superexchange interactions in the ferromagnetic planes of LaMnO.Fei99

Since the isotropic model has the same interaction strength for the bonds along both and axis, it is symmetric under transformation , and the issue of the QPT between different ground states of the anisotropic compass modelOru09 does not arise. On one hand, this symmetry is obeyed by the classical Ising ground state, while on the other hand, in the ground state of the OCM this symmetry is spontaneously broken (and the ground state is degenerate). Therefore, an intriguing question concerning the ground state of the model (1) is whether it has the same high symmetry as the Ising model in a broad range of , or the symmetry is soon spontaneously broken when increases, i.e., there are degenerate ground states with lower symmetries, also for the orbital model, see Fig. 1(b). This question has been addressed by investigating the energy contributions along two equivalent lattice directions and by applying the MERA.

## Iii MERA Algorithm

### iii.1 Calculation method

In order to obtain the ground state, we use a translationally invariant MERA on infinite lattice.Eve09 The MERA is a tensor network with infinite number of layers of disentanglers and isometries. By translational invariance, all isometries (disentanglers) in a given layer are the same. Since every layer represents a coarse-graining renormalization group transformation, shown in Fig. 2 for the 9-to-1 geometry and described in more detail in Ref. Eve09, (see their Fig. 7), we assume that after a finite number of such transformations a fixed point of the renormalization group is reached (either trivial or non-trivial) and from that time on the following transformations are the same. In other words, at the bottom of the tensor network there is a finite number of non-universal layers whose tensors are different in general, but above certain level all layers are the same. The bottom layers describe non-universal short range correlations, and the universal layers above this level describe universal properties of the fixed point. The number of the non-universal bottom layers is one of the parameters of the infinite-lattice MERA. We have verified that it is enough to keep up to three non-universal layers, depending on how close the critical point is.

Starting with randomly chosen tensors, the structure is optimized layer by layer, from the top to bottom and back. In given layer , we calculate an environment of each tensor type by means of renormalized Hamiltonians and density matrices computed from other layers. The environments are aimed at updating tensors to minimize total energy. In the universal layer, this updating technique is slightly different: and are fixed points of the renormalization procedure defined by tensors in this layer. The above steps are iterated until the convergence of energy is achieved. For given , we obtain the ground state for different values of bond dimension . It turns out that in most cases it is sufficient to work with , which is the same in each layer. However, it is necessary to increase to in the neighborhood of the critical point. The number of operations and the required memory scale as and respectively. The inset in Fig. 4(a) shows the convergence of the energy of the ground state with an increasing bond dimension. Here we also present a comparison of results obtained with the alternative 5-to-1 geometry.

The algorithm is implemented in `c`

++ and optimized in order to
work on multi-processor computers. On an eight-core 2.3 GHz
processor, it takes about half an hour to update the whole tensor
network which consists of four layers of tensors with . Near
the critical point, i.e., at , the convergence
requires several thousands of iterations whereas it is significantly
faster far from . When is scanned from
to (or back), it is more efficient to use the previous
ground state as an initial state for the next discrete value of
instead of starting from a random initial state for each
value of . We have carefully verified convergence to the
ground state by scanning back and forth and comparing the
results with those obtained from random initial states for selected
values of .

### iii.2 Correlations

In order to calculate correlations, we take advantage of the special structure of the renormalization group transformation in Fig. 2. A site of the lattice that lies in the center of a decimation block (number 9 in Fig. 2) undergoes renormalization in a particularly easy manner. Since no disentangler is applied to this central site, a one-site operator at this site is mapped by the -th renormalization group transformation to a coarse-grained one-site operator

(4) |

Here is a renormalizer superoperator built out of contracted isometries only, as shown in Fig. 3(a):

(5) |

The meaning of the transformations and is given in Figs. 2 and 3(a). Thus, if we have non-universal layers at the bottom of the geometry of MERA, then renormalized one-site operators at the central sites just below the universal layer are given by [see Fig. 3(b)]:

(6) |

where denotes a physical, microscopic one-site operator at one of the central sites at the very bottom of the MERA tensor network.

To extract information on the correlations, it is convenient to write eigen-decomposition of the renormalizer in the universal layer:

(7) |

It is straightforward to verify the basic property of the spectrum of : . The ortonormality of the vectors in Eq. (5) implies that the identity operator is an eigenvector with eigenvalue . In our numerical calculations this is the only eigenvalue with modulus .

After the operator is decomposed as , a repeated action of the renormalizer in the universal layers can be written as

(8) |

A correlator between two central sites and separated by a distance in the horizontal (vertical) direction is thus given by:

(9) | |||||

(10) | |||||

(11) |

where and

(12) |

Here is a two-site reduced density matrix in a universal layer derived from as depicted in Fig. 3(c).

Correlations corresponding to the leading eigenvalue do not decay with the distance between and . They describe long range order in the operator and can be used to extract its expectation value :

(13) |

where we use the property: that holds for and the fact that which is a consequence of being an identity. Thus only a one-site operator with a non-zero coefficient has non-zero expectation value. A trivial example is the identity . Indeed, we obtain in Eq. (6), which is equivalent to , and Eq. (13) yields as expected.

## Iv Numerical results

### iv.1 Symmetry breaking transition

Information about the ground state of the OCM Eq. (1) is contained in average energy per bond and energy anisotropy :

(14) | |||||

(15) |

In the classical limit of Ising interactions and . Due to increasing frustration, the energy gradually increases for increasing angle in Eq. (1) and reaches a maximum of in the OCM, see Fig. 4(a). This increase is smooth and does not indicate the existence of a QPT.

However, by investigating the anisotropy Eq. (15) between and bonds, we identified an angle at which starts to grow. Although a gradual evolution of the ground state staring from might be also expected, the Ising-type state is first surprisingly robust in a broad range of angles , and the energy associated with bonds along the and axes remains the same, i.e., . Next, the symmetry between the and directions is spontaneously broken above , where a finite value of is found, and then grows rapidly with further increasing angle , i.e., large spin correlations develop along only one of the two equivalent directions and . This QPT was detected by the MERA at , see Fig. 4(b).

### iv.2 Magnetization in the ground state

To understand better the QPT at let us consider the expectation value of the spontaneous magnetization derived from the long range order in the correlation function:

(16) |

where . For the interactions in Eq. (1) one finds for any .

We found that the ground state obtained using the MERA for is characterized by and Ising-type long range order of which gradually decreases but remains rather large, , in this parameter range. The symmetry between the directions and is broken above by appearance of a nonzero component .

The value of the total magnetization

(17) |

obtained from the MERA decreases continuously from in the Ising model to in the OCM, see Fig. 4(c). Thus the reduction in the order parameter by quantum fluctuations arising from the admixture of the -th component, is here rather small, and reproduces qualitative results obtained for the orbital model within the linear orbital wave theory vdB99 . Furthermore, by a closer inspection of we have found that the derivative does not exist at .

As expected from the behavior of , the obtained symmetry breaking shown in Fig. 4 implies that the direction of spontaneous magnetization , parametrized by an orientation angle

(18) |

begins to change when increases above , see Fig. 4(d). For , the magnetization has only one component with , pointing either parallel or anti-parallel to which is half-way between and , see Fig. 1. Below the ferromagnetic ground state is doubly degenerate and the magnetization is . When increases above the magnetization begins to rotate in the plane by the non-zero angle Eq. (18) with respect to the initial magnetization below , and each of these two states splits off into two ferromagnetic states rotated by with respect to the -axis. As a result, one finds four degenerate states above , and each of them is tilted with respect to , either toward , or toward , depending on the sign of the rotation angle . In the OCM limit is approached, the magnetization angle approaches . In this limit there are four degenerate Ising-type ferromagnetic states, with magnetization either along (and ), or (and ).

Qualitatively the same results were obtained from the embedded clusters and they are also shown in Fig. 4 for comparison. While cluster is too small and the quantum fluctuations are severely underestimated, the two larger and clusters are qualitatively similar and estimate the QPT point from above, see Fig. 4. Rather slow convergence of these results toward the MERA result for and for demonstrates the importance of longer-range correlations for the correct description of the QPT at .

Altogether, these results show that the degenerate ground state of the generalized OCM consists of a manifold of states with broken symmetry. This confirms that the OCM is in the Ising universality class,Mis04 ; Dor05 with no quantum coupling between different broken symmetry Ising-type states. However, we found the large value of rather surprising and we investigated it further using spin-wave theory. These calculations are presented in the next Section.

Another surprise is the absence of any algebraically decaying spin-spin correlations in the MERA ground state at . They could arise from the subleading eigenvalues which we found to be non-zero. However, their corresponding coefficients with or in Eq. (11) are small (at most ) and they decay with increasing dimension and especially the number of non-universal layers . As a result, the only nonvanishing term in Eq. (11) is the leading one for , describing the non-decaying long range order. Notice that this observation does not exclude non-trivial short range correlations up to a distance described by the non-universal layers. We believe that when is too small, then the missing short range correlations find a way to show up in the small but non-zero universal coefficients , but these coefficients decay quickly with increasing as the short range correlations become accurately described by the increasing number of non-universal layers.

## V Spin wave expansion

Since the spin wave expansion in powers of becomes exact when the spin , we introduced a large- extension of the generalized OCM Hamiltonian Eq. (1) with rescaled spin operators: and . We consider first the classical energy per site:

(19) |

obtained using the mean-field (MF) for the ordered state of classical spins , with the magnetization direction given by Eq. (18). The classical energy has a minimum at for the entire range of . However, when the angle approaches , the minimum becomes more and more shallow, and finally disappears completely at . Thus, the classical ground state becomes very sensitive to quantum fluctuations in the vicinity of the maximally frustrated interactions in the OCM.

This behavior of the classical ground state energy explains why small energy contributions due to quantum fluctuations may play so crucial role in the generalized OCM only in the regime of close to , where they trigger a QPT by splitting the shallow symmetric classical energy minimum at into two symmetry-broken minima at finite values — we show an example of this behavior in Fig. 5 for a particular value of . Since the quantum fluctuations induce here symmetry breaking instead of making the ground state more symmetric, this mechanism goes beyond the Landau functional paradigm.

We analyzed the effects of quantum fluctuations and the arising symmetry breaking using the Holstein-Primakoff representation of spin operators via bosons:

(20) | |||||

(21) |

Operators satisfy standard bosonic commutation relations: and . In this approach, we are looking for a critical value , above which it is energetically favorable to change the direction of magnetization from the symmetric state to a symmetry-broken state with a finite value of . We expanded the square root in Eq. (21) in powers of and obtained an expansion of Hamiltonian Eq. (1) in powers of the operators . As we applied Wick’s theorem to reduce the obtained Hamiltonian to an effective quadratic Hamiltonian, the terms proportional to the odd powers of do not contribute and are skipped below (for more details see the Appendix). When truncated at the sixth order term this expansion reads

(22) |

Here is a sum of all terms of the -th order in operators. In a similar way, and denote expansions truncated at the fourth and second order terms, respectively. We have found a posteriori that the second order expansion (noninteracting spin waves) does not suffice and higher order terms are necessary. Consequently, we consider below Hamiltonian Eq. (1) expanded up to the sixth order.

For given and , we can approximate the ground state of the boson Hamiltonian given by Eq. (22) by a Bogoliubov vacuum obtained as the ground state of the quadratic Hamiltonian obtained using the MF averaging of four- and six-boson terms. Details of this calculation can be found in the Appendix A.

First we performed separate calculations for , and for several values of spin when the -expansion given in Eq. (22) is convergent. The quadratic fails for large , where its spectrum becomes gapless and the magnetization Eq. (17) diverges. In contrast, and give only small reduction in in the entire range of , see Fig. 6(a) and 6(c). Interestingly, the Bogoliubov spectrum remains gapful at in both the fourth and sixth order expansions and, just like in the MERA, there are no algebraically decaying spin-spin correlations. The critical angle at which the symmetry-breaking QPT occurs increases toward with increasing when the quantum fluctuations become less significant. Therefore, the magnetization increases with increasing and it tends to in the classical limit .

Encouraged by these results, we also performed similar calculations for the generalized OCM Eq. (1) with , where the convergence of the -expansion becomes problematic. Unlike for , we find that the fourth order expansion is insufficient as it predicts the first order QPT [Fig. 6(d)] and does not agree qualitatively with the prediction of the MERA, see Sec. IV. Only in the sixth order one finds a qualitative agreement between the present boson expansion and the MERA, both giving the second order QPT at . A cusp in seen in Fig. 6(c) shows that even the sixth order expansion is not quite converged for . Again, the Bogoliubov spectrum remains gapful at in the sixth order expansion, with a finite gap equal , and one finds no algebraically decaying spin-spin correlations.

## Vi Conclusions

Summarizing, we found that a second order quantum phase transition in the generalized orbital compass model Eq. (1) occurs at which is surprisingly close to the compass point , i.e., only when the interactions are sufficiently strongly frustrated. There is spontaneous ferromagnetic magnetization at any angle . Below the ferromagnetic ground state is doubly degenerate with the spontaneous magnetization, either parallel or antiparallel to the average direction . None of the directions, neither nor , is preferred in this symmetric phase. In contrast, when increases above the symmetry between and becomes spontaneously broken and the ferromagnetic magnetization begins to align parallel/antiparallel to either or . The ground state is fourfold degenerate in this symmetry-broken phase. The spontaneous magnetization is close to and quantum fluctuations remain small in the whole range of .

These results were obtained using the MERA and the mechanism of the QPT was explained within the spin-wave theory. For classical spins the minimum of energy is at one of the two symmetric states with the magnetization either parallel or antiparallel to , see Fig. 5. The minimum becomes more and more shallow as the compass point is approached. However, the quantum fluctuations are weak due to the gapful orbital wave excitations, and only very close to the above OCM point become strong enough to split the shallow minimum into two distinct minima in the vicinity of the OCM point. In this way the symmetry between the axes and is spontaneously broken. For this reason the orbital model with ferro-orbital interactions, considered in Ref. vdB99, and corresponding to a “moderate” value of [see Fig. 1(b)], orders in a symmetric (uniform) phase induced by the stronger (here ) interaction component.

Interestingly, since — unlike in the Landau paradigm — the symmetry in the present model Eq. (1) is broken rather than restored by quantum fluctuations, we do not find any algebraically decaying spin-spin correlations at the critical point found in the generalized orbital compass model Eq. (1). The spin waves also remain gapful at this point.

###### Acknowledgements.

We thank P. Horsch for insightful discussions. L.C., J.D. and A.M.O. acknowledge support by Polish Ministry of Science and Higher Education under Projects No. N202 175935, N202 124736 and N202 069639, respectively. A.M.O. was also supported by the Foundation for Polish Science (FNP). *## Appendix A Details of the spin wave calculation

In this section we present the details of the spin wave calculation. We consider the general case of an square lattice, with being odd for convenience. The results presented in Section V are obtained after taking the thermodynamic limit .

For given and , the ground state of the boson Hamiltonian Eq. (22) is approximated by a Bogoliubov vacuum obtained as the ground state of a mean-field (MF) quadratic Hamiltonian (to be derived later on). Terms , and in Eq. (22) are given by:

(23) |

(24) |

(25) |

where , and . The signs mean here that both terms, with and sign separately, must be taken into account.

To derive the quadratic approximation , we replace the boson terms in and with two-boson terms and proper averages by means of the MF approximation and Wick’s theorem. This justifies a posteriori why the expansion (22) is limited only to the terms with even number of boson operators. As an example of this approximation, consider one of the contributions to in Eq. (24): , which is replaced with a quadratic term:

(26) |

The above replacement procedure leads to six MF parameters that should satisfy self-consistency conditions. These are in fact all possible combinations of operators defined on nearest-neighbor sites that cannot be derived one from another by commutation relations and translational invariance of the lattice, i.e., , , , , , and .

The obtained Hamiltonian is diagonalized by the Fourier transformation followed by the Bogoliubov transformation. Fourier transformation which is consistent with periodic boundary conditions and has the following form:

(27) |

where is the momentum. In the sum, momentum components and take the values (for odd considered here):

(28) |

Diagonalization of is completed by the Bogoliubov transformation:

(29) |

where the modes and are normalized such that . The obtained modes are used to calculate new values of the MF parameters (). For instance, one of them reads: . Starting from random values, the above steps are iteratively applied until full convergence of all is reached, which results in satisfying the self-consistency conditions.

## References

- (1) Y. Tokura and N. Nagaosa, Science 288, 462 (2000); A. Weiße and H. Fehske, New J. Phys. 6, 158 (2004); E. Dagotto, ibid. 7, 67 (2005).
- (2) P. Horsch, A. M. Oleś, L. F. Feiner, and G. Khaliullin, Phys. Rev. Lett. 100, 167205 (2008); A. M. Oleś and P. Horsch, in: Properties and Applications of Thermoelectric Materials — The Search for New Materials for Thermoelectric Devices, ed. V. Zlatic and A. C. Hewson (Springer, New York, 2009), p. 299; A. M. Oleś, Acta Phys. Pol. A 115, 36 (2009).
- (3) L. F. Feiner, A. M. Oleś and J. Zaanen, Phys. Rev. Lett. 78, 2799 (1997).
- (4) K. I. Kugel and D. I. Khomskii, Sov. Phys. Usp. 25, 231 (1982); D. I. Khomskii and M. V. Mostovoy, J. Phys. A 36, 9197 (2003).
- (5) Z. Nussinov, M. Biskup, L. Chayes, and J. van den Brink, Europhys. Lett. 67, 990 (2004).
- (6) A. Mishra, M. Ma, F.-C. Zhang, S. Guertler, L.-H. Tang, and S. Wan, Phys. Rev. Lett. 93, 207201 (2004).
- (7) J. Dorier, F. Becca, and F. Mila, Phys. Rev. B72, 024448 (2005); S. Wenzel and W. Janke, ibid. 78, 064402 (2008).
- (8) B. Douçot, M. V. Feigel’man, L. B. Ioffe, and A. S. Ioselevich, Phys. Rev. B71, 024505 (2005); Z. Nussinov and E. Fradkin, ibid. 71, 195120 (2005).
- (9) G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 102, 017205 (2009).
- (10) A. Nagano, M. Naka, J. Nasu, and S. Ishihara, Phys. Rev. Lett. 99, 217202 (2007).
- (11) A. Kitaev, Ann. Phys. (N.Y.) 321, 2 (2006).
- (12) Z. Nussinov and G. Ortiz, Ann. Phys. (N.Y.) 324, 977 (2009).
- (13) P. Milman, W. Maineult, S. Guibal, L. Guidoni, B. Douçot, L. Ioffe, and T. Coudreau, Phys. Rev. Lett. 99, 020503 (2007).
- (14) W. Brzezicki, J. Dziarmaga, and A. M. Oleś, Phys. Rev. B75, 134415 (2007); E. Eriksson and H. Johannesson, ibid. 79, 224424 (2009).
- (15) R. Orús, A. C. Doherty, and G. Vidal, Phys. Rev. Lett. 102, 077203 (2009).
- (16) G. Vidal, Phys. Rev. Lett. 99, 220405 (2007).
- (17) G. Vidal, Phys. Rev. Lett. 101, 110501 (2008).
- (18) S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, U.K., 1999).
- (19) S. R. White, Phys. Rev. Lett. 69, 2863 (1992); Phys. Rev. B48, 10345 (1993).
- (20) K. Hallberg, Adv. Phys. 55, 477 (2006).
- (21) L. Cincio, J. Dziarmaga, and M. M. Rams, Phys. Rev. Lett. 100, 240603 (2008).
- (22) G. Evenbly and G. Vidal, Phys. Rev. Lett. 102, 180406 (2009).
- (23) G. Evenbly and G. Vidal, Phys. Rev. B79, 144108 (2009).
- (24) V. Giovannetti, S. Montangero, and R. Fazio, Phys. Rev. Lett. 101, 180503 (2008); V. Giovannetti, S. Montangero, M. Rizzi, and R. Fazio, Phys. Rev. A79, 052314 (2009); S. Montangero, M. Rizzi, V. Giovannetti, and R. Fazio, Phys. Rev. B80, 113103 (2009); P. Silvi, V. Giovannetti, P. Calabrese, G. E. Santoro, and R. Fazio, J. Stat. Mech. L03001 (2010).
- (25) C. M. Dawson, J. Eisert, and T. J. Osborne, Phys. Rev. Lett. 100, 130501 (2008).
- (26) T. Barthel, C. Pineda, and J. Eisert, Phys. Rev. A80, 042333 (2009); C. Pineda, T. Barthel, and J. Eisert, Phys. Rev. A81, 050303 (2010).
- (27) G. Evenbly and G. Vidal, Phys. Rev. B81, 235102 (2010).
- (28) G. Evenbly and G. Vidal, Phys. Rev. Lett. 104, 187203 (2010).
- (29) The orbital interactions in transition metal oxides favor alternating orbital order, see Refs. Tok00, ; Hor08, — the corresponding OCM is equivalent to the present ferro-orbital model Eq. (1) by a unitary transformation.
- (30) L. F. Feiner and A. M. Oleś, Phys. Rev. B59, 3295 (1999).
- (31) J. van den Brink, P. Horsch, F. Mack, and A. M. Oleś, Phys. Rev. B59, 6795 (1999).
- (32) The comparison has been made for geometries with the same bond dimension in all layers. The performance of 5-to-1 geometry can be further improved by considering different bond dimensions on some indices.