Emergent Haldane phase in the $S=1$ bilinear-biquadratic Heisenberg model on the square lattice

TL;DR

This study uses advanced tensor network simulations to discover an unexpected Haldane phase in a 2D spin model, challenging prior assumptions of direct magnetic transitions.

Contribution

It reveals an emergent Haldane phase in the square lattice $S=1$ bilinear-biquadratic Heisenberg model, showing a new intermediate phase not previously identified.

Findings

Discovery of an emergent Haldane phase between magnetic orders

The phase preserves SU(2) and translational symmetry

Contradicts earlier studies suggesting direct transitions

Abstract

Infinite projected entangled pair states simulations of the $S=1$ bilinear-biquadratic Heisenberg model on the square lattice reveal an emergent Haldane phase in between the previously predicted antiferromagnetic and 3-sublattice 120$^\circ$ magnetically ordered phases. This intermediate phase preserves SU(2) spin and translational symmetry but breaks lattice rotational symmetry, and it can be adiabatically connected to the Haldane phase of decoupled $S=1$ chains. Our results contradict previous studies which found a direct transition between the two magnetically ordered states.

Equations

$$\mathcal{H}=\sum _{⟨i,j⟩}\left[\cos (\theta ){\mathbf{S}}_{\mathbf{i}}\cdot {\mathbf{S}}_{\mathbf{j}}+\sin (\theta ){\left({\mathbf{S}}_{\mathbf{i}}\cdot {\mathbf{S}}_{\mathbf{j}}\right)}^2\right],$$

Full text

The search for novel states of matter in quantum many-body systems is one of the most active areas in condensed matter physics. A fascinating example is the ground state of the $S=1$ antiferromagnetic Heisenberg chain which, unlike the $S=1/2$ chain, exhibits an energy gap, exponentially decaying spin-spin correlations, and gapless edge excitations in case of open boundaries. Thanks to Haldane’s pioneering work and conjecture that such a gapped state emerges in integer Heisenberg spin chains in general Haldane (1983a, b), this phase has been named after him.

The Haldane phase also extends to related $S=1$ models, such as weakly-coupled $S=1$ Heisenberg chains Sakai and Takahashi (1989); Koga and Kawakami (2000); Kim and Birgeneau (2000); Matsumoto et al. (2001); Wierschem and Sengupta (2014) which are realized in several quasi-1D materials Meyer et al. (1982); Renard et al. (1987); Monfort et al. (1996); Honda et al. (1997, 1998); Uchiyama et al. (1999), or the $S=1$ bilinear-biquadratic Heisenberg (BBH) chain with Hamiltonian

[TABLE]

for $\theta$ in between $-\pi/4$ and $\pi/4$. More recently, the Haldane phase has been understood as a simple example of a symmetry protected topological (SPT) phase Gu and Wen (2009); Pollmann et al. (2010, 2012); Chen et al. (2013).

In the present work we focus on the BBH model in two dimensions, which has gained much interest in recent years Papanicolaou (1988); Harada and Kawashima (2002); Tsunetsugu and Arikawa (2006); Läuchli et al. (2006); Tsunetsugu and Arikawa (2007); Bhattacharjee et al. (2006); Tóth et al. (2012); Bieri et al. (2012); Oitmaa and Hamer (2013); firstly, due to its possible connection to the triangular lattice compounds NiGa2S4 Nakatsuji et al. (2005) and Ba3NiSb2O9 COB ; Cheng et al. (2011); Bieri et al. (2012); Fak et al. (2017), and secondly, for $\theta=\pi/4$, the model is equivalent to the SU(3) Heisenberg model which can be experimentally realized using ultra-cold fermionic atoms in optical lattices Wu et al. (2003); Honerkamp and Hofstetter (2004); Cazalilla et al. (2009); Gorshkov et al. (2010). The latter has been shown to exhibit 3-sublattice order on the square and triangular lattices Tóth et al. (2010); Bauer et al. (2012), and an important question concerns the stability of this phase away from the SU(3) symmetric point. Previous studies on the square lattice based on linear flavor-wave theory Tóth et al. (2012), exact diagonalization Tóth et al. (2012), and series expansion Oitmaa and Hamer (2013) predicted a direct transition between the AF and the 3-sublattice phase for $\theta\approx 0.2\pi$. However, the accurate study of this parameter regime remains very challenging because Quantum Monte Carlo suffers from the negative sign problem.

In this paper we show, using state-of-the-art tensor network simulations, that in between the AF and the 3-sublattice phase an intermediate quantum paramagnetic phase emerges which preserves translational and SU(2) spin symmetry, but breaks lattice rotational symmetry (see Fig. 1). We identify this intermediate phase as the Haldane phase by showing that it can be adiabatically connected to the Haldane phase of decoupled $S=1$ chains. This result at first appears surprising in view of the fact that for $\theta=0$ already a small interchain coupling $J_{y}>J^{c}_{y}=0.0436$ Matsumoto et al. (2001) is sufficient to destabilize the Haldane phase. However, we show that with increasing $\theta$ the critical interchain coupling $J^{c}_{y}(\theta)$ separating the Haldane phase from the AF phase dramatically increases, and eventually reaches the isotropic 2D limit.

Method – Our results have been obtained using infinite projected entangled-pair states Verstraete and Cirac (2004); Jordan et al. (2008); Corboz et al. (2010) (iPEPS, also called tensor product state Nishino et al. (2001); Nishio et al. (2004)) — a variational tensor-network ansatz to efficiently represent ground states in two dimensions directly in the thermodynamic limit. It can be seen as a natural generalization of matrix product states (the underlying ansatz of the powerful density-matrix renormalization group method White (1992)) to two dimensions. This approach has already been applied successfully to a variety of challenging problems in the field of frustrated magnetism and strongly correlated electrons, see e.g. Refs. Corboz et al. (2011a); Wang et al. (2016); Zhao et al. (2012); Corboz et al. (2012); Xie et al. (2014); Corboz and Mila (2013); Gu et al. (2013); Matsuda et al. (2013); Osorio Iregui et al. (2014); Corboz et al. (2014); Corboz and Mila (2014); Picot and Poilblanc (2015); Picot et al. (2016); Corboz (2016); Nataf et al. (2016); Liao et al. (2016); Zheng et al. (2016) and references therein.

The iPEPS ansatz consists of a network of order-5 tensors on a square lattice, with one tensor per lattice site. Each tensor has a physical index carrying the local Hilbert space of a site and four auxiliary indices which connect to the four nearest-neighboring tensors. Each auxiliary index goes over $D$ elements, called the bond dimension, which controls the accuracy of the ansatz. For translational invariant states an ansatz with the same tensor on each lattice site can be used, however, if translational symmetry in the ground state is broken, a larger unit cell of tensors is required. For example, to reproduce an antiferromagnetic state two different tensors (one for each sublattice) is needed, whereas the 3-sublattice 120*∘* ordered state shown in Fig. 1 requires a unit cell with 3 different tensors. In practice we run simulations using different unit cells to find out which structure yields the lowest energy state.

For more details on the method we refer to Refs. Corboz et al. (2010); Corboz and Mila (2013); Phien et al. (2015). For the experts we note that the optimization of the tensors has been done via an imaginary time evolution with the full update optimization Corboz et al. (2010) (or fast-full update Phien et al. (2015)), except for the simulations of the anisotropic model where we used the computationally cheaper simple-update optimization Vidal (2003); Jiang et al. (2008). The contraction of the infinite tensor network is done by a variant Corboz et al. (2011b, 2014) of the corner-transfer matrix method Nishino and Okunishi (1996); Orús and Vidal (2009). We also exploited the U(1) symmetry Singh et al. (2011); Bauer et al. (2011) to increase the efficiency (except in the 3-sublattice phase which breaks U(1) symmetry).

AF and 3-sublattice phases – We first discuss the well-established limits of the phase diagram in the range $\theta\in[0,\pi/4]$. For $\theta=0$ the model reduces to the $S=1$ Heisenberg model where the ground state exhibits antiferromagnetic order. Unlike for $\theta>0$, Quantum Monte Carlo has no negative sign problem in this limit, and the sublattice magnetization $m=\sqrt{\langle S_{x}\rangle^{2}+\langle S_{y}\rangle^{2}+\langle S_{z}\rangle^{2}}$ in the thermodynamic limit has been accurately determined: $m=0.805(2)$ Matsumoto et al. (2001). Our iPEPS result extrapolated to the infinite $D$ limit (see Fig. 2), $m=0.802(7)$, is in agreement with this value. For $\theta=\pi/4$ the model is equivalent to the $SU(3)$ Heisenberg model (with the fundamental representation on each site), for which a 3-sublattice ordered state has been predicted by several methods Tóth et al. (2010); Bauer et al. (2012), including previous iPEPS simulations. For $\theta$ slightly below $\pi/4$ this order corresponds to a $120^{\circ}$ order formed by the spins on 3 sublattices (see Fig. 1).

Intermediate phase – In Refs. Tóth et al. (2012); Oitmaa and Hamer (2013) a direct transition between the AF state and the 3-sublattice state has been predicted to occur around $\theta\approx 0.2\pi$ based on exact diagonalization, linear flavor-wave theory and series expansion. We first attempt to reproduce this result with iPEPS by determining the critical value $\theta_{c}$ for which the energies of the two states — distinguished by different unit cells — intersect. To do so, we initialize a simuluation from deep within the AF (3-sublattice) phase, and slowly increase (decrease) $\theta$ up to the point where the energies cross. The resulting critical value $\theta_{c}\approx 0.21\pi$ (for $D=10$) is close to the previous prediction. However, from a systematic analysis of the AF order parameter shown in Fig. 2 we find that the AF order actually vanishes long before $\theta_{c}$, i.e. that the AF phase is only stable up to $\theta=0.189(2)\pi$ COI . This indicates the presence of an intermediate non-magnetic phase in between the AF- and the 3-sublattice phase.

We next explore the region around $\theta=0.2\pi$ in more detail. When starting from different random initial states with a 2-sublattice ansatz, we observe a competition between a weakly magnetized state and a non-magnetized state which breaks lattice rotational symmetry but preserves SU(2) and translational symmetry CVB . This non-magnetized state can also be found by restricting the simulation to a 1-site unit cell CU ; it exhibits the lowest variational energy for large $D$.

The rotational symmetry breaking manifests itself in energy differences in $x$ and $y-$direction, as illustrated by the different thicknesses of the bonds in the middle panel of Fig. 1, reminiscent of coupled 1D chains. Since the ground state of the BBH chain for $\theta\in[0,\pi/4]$ is in the Haldane phase, the question naturally arises whether the intermediate 2D phase could possibly be adiabatically connected to the Haldane phase by continuously decreasing the $y$-coupling to zero. A first hint that this picture is correct comes from the observation that when initializing the iPEPS as a product of chains in the 1D Haldane phase the simulation converges to the same non-magnetized state as with randomly initialized tensors. In order to confirm this picture we study the stability of the Haldane phase in the anisotropic BBH model in the following.

Anisotropic model – We introduce different coupling strengths in $x$ and $y-$direction and study the phase diagram in the $(\theta,J_{y})$ plane (setting $J_{x}=1$ in the following). For $J_{y}=0$ the model simply reduces to independent $S=1$ chains, which are known to lie in the Haldane phase (for $\theta\in[0,\pi/4[$). The goal is now to determine the critical coupling $J_{y}^{c}(\theta)$ separating the Haldane phase from the AF phase (or the 3-sublattice phase for large $\theta$), for different values of $\theta$. In order to obtain an estimate of the phase transition for a fixed value of $\theta$ we initialize the iPEPS in the Haldane phase and in the AF phase (or 3-sublattice phase) respectively, run simulations for different values of $J_{y}$, and determine the value $J_{y}^{c}(\theta)$ for which the energies of the two states intersect, using a fixed bond dimension $D=10$ and the simple update optimization (see data in the supplemental material SUP ). We note that this approach provides only an approximate phase boundary, in contrast to the extrapolated full update simulations used in the isotropic case. However, it is computationally much more efficient, which becomes significant when probing the extended two-dimensional parameter space $(\theta,J_{y})$. Moreover, a comparison with Monte Carlo and full update results (see below) indicates that the relative error on the phase boundary is only a few percent, which is accurate enough for our current purpose.

For $\theta=0$ we find a critical value $J_{y}^{c}(0)=0.042$ which is close to the Quantum Monte Carlo result 0.0436 from Ref. Matsumoto et al. (2001) (see also Refs. Sakai and Takahashi (1989); Koga and Kawakami (2000); Kim and Birgeneau (2000); Wierschem and Sengupta (2014)). This value lies far away from the isotropic 2D limit $J_{y}=1$. However, we find that $J_{y}^{c}(\theta)$ monotonously increases with $\theta$ as shown in the phase diagram in Fig. 3, and that beyond $\theta_{c}=0.200\pi$ no phase transition occurs, showing that the 1D Haldane phase can indeed be adiabatically connected to the isotropic 2D limit. Finally, for $0.213\pi\leq\theta<\pi/4$ we find a finite transition value $J_{y}^{c}(\theta)$ between the Haldane and the 3-sublattice phase which decreases with increasing $\theta$.

Comparing to the full update results (Fig. 2), which predict the two transitions to be at $0.189(2)\pi$ and $0.217(4)\pi$ respectively, we see that the simple update underestimates the extent of the Haldane phase at the isotropic point. Moreover, it does so by a margin of at most $0.01\pi$ (and by much less for $\theta=0$), indicating that the continuous path that connects the intermediate 2D phase to the 1D Haldane phase persists also when taking the error margin on the phase boundary into account.

Transition from Haldane to 3-sublattice phase – We next focus again on the isotropic 2D case ($J_{y}=1$) and accurately determine the transition from the Haldane to the 3-sublattice phase by pushing the simulations up to $D=16$ (Haldane state) and $D=10$ (3-sublattice state) using the full update optimization, and compare the energies of the two states in the infinite $D$ limit. Figure 4(a) shows the energies extrapolated in the so-called truncation error $w$ (see Ref. Corboz (2016) for details). For $\theta=0.21\pi$ the state in the Haldane phase is clearly lower than the 3-sublattice state, whereas for $\theta=0.22\pi$ the opposite is true. By linear interpolation of the energies, taking into account the extrapolation error, we find a critical value of $\theta_{c}=0.217(4)\pi$. Finally, the squares in Fig. 4(b) show the difference in bond energies $\Delta E=E_{y}-E_{x}$ in $x$ and $y-$direction of the Haldane state. In the infinite $D$ limit $\Delta E$ tends to a finite value, e.g. $\Delta E=0.07(1)$ for $\theta=0.21\pi$, which shows that the rotational symmetry in the Haldane phase is indeed spontaneously broken.

Nature of the phase transitions in isotropic case – Because the Haldane and 3-sublattice phases break different translational and rotational symmetries, the corresponding phase transition is expected to be first order. This picture is confirmed by the occurrence of hysteresis around the transition point, which allows us to simulate both phases on both sides of the phase transition. Moreover, the sublattice magnetization is strictly positive throughout the 3-sublattice phase (see Fig. 4(b)) — even for $\theta=0.21\pi$ where the 3-sublattice state is no longer the lowest energy state — implying that the magnetization does not go to zero when approaching the transition from above. Since the magnetization is zero in the Haldane phase, it jumps to zero at the transition, showing that the transition is clearly of first order.

As for the AF to Haldane phase transition, the absence of a clear hysteresis in the full update simulations and the fact that the sublattice magnetization in the AF phase goes to zero as we approach the Haldane phase indicates either an unconventional second order or a weak first order phase transition. However, due to the error bars in Fig. 2(b) close to the critical point, we cannot exclude one of the two based on our data.

Conclusion – We have studied the $S=1$ BBH model on a square lattice where, in contrast to previous predictions, we found an intermediate quantum paramagnetic phase between the AF and 3-sublattice 120∘ magnetically ordered phases in the parameter range $\theta/\pi\in[0.189(2),0.217(4)]$. This intermediate phase is characterized by (1) translational symmetry, (2) an absence of magnetic and quadrupolar order, i.e. SU(2) spin symmetry is preserved, and (3) a spontaneous rotational symmetry breaking with stronger bonds in $x$ (or $y$) direction (i.e. lattice nematic order). The above features are reminiscent of the ones of weakly-coupled $S=1$ chains in the Haldane phase, which motivated us to study the anisotropic BBH model. With increasing $\theta$ we found that the critical coupling $J^{c}_{y}(\theta)$ separating the Haldane and AF phases monotonically increases up to $\theta=0.189(2)\pi$, after which no phase transition occurs as a function of $J_{y}$, i.e. the Haldane phase persists all the way up to the isotropic 2D limit. From this we identified the intermediate phase as a continuous 2D extension of the Haldane phase.

It is interesting to note that a similar situation has previously been encountered in the $J_{1}$-$J_{2}$ $S=1$ Heisenberg model on a square lattice Jiang et al. (2009), in which an intermediate Haldane phase between an AF and a stripe phase appears that also survives up to the isotropic limit. Moreover, our findings provide an additional example of a nematic quantum paramagnet which in Ref. Wang et al. (2015) was proposed to likely emerge in spin-1 systems with competing interactions and suggested to be potentially relevant to understand the nematic phase in the iron-based superconductor FeSe.

Finally, our results further highlight the potential of iPEPS as a powerful tool for challenging open problems in frustrated magnetism where Quantum Monte Carlo suffers from the negative sign problem. As a future work it will be interesting to see whether the Haldane phase can also be found in or nearby the isotropic 2D limit on the triangular lattice, potentially offering further understanding of the unusual behavior of NiGa2S4 Nakatsuji et al. (2005) and Ba3NiSb2O9 Cheng et al. (2011).