Overcoming the Sign Problem at Finite Temperature: Quantum Tensor Network for the Orbital $e_g$ Model on an Infinite Square Lattice

TL;DR

This paper demonstrates that tensor network methods can effectively overcome the quantum Monte Carlo sign problem at finite temperature for the orbital $e_g$ model, accurately capturing phase transition properties and critical exponents.

Contribution

It introduces a tensor network approach to study a sign-problematic 2D quantum model at finite temperature, providing precise estimates of critical temperature and exponents.

Findings

Confirmed finite order parameter below $T_c$

Estimated critical temperature $T_c$ accurately

Critical exponents within 1% of 2D Ising universality

Abstract

The variational tensor network renormalization approach to two-dimensional (2D) quantum systems at finite temperature is applied for the first time to a model suffering the notorious quantum Monte Carlo sign problem --- the orbital $e_g$ model with spatially highly anisotropic orbital interactions. Coarse-graining of the tensor network along the inverse temperature $\beta$ yields a numerically tractable 2D tensor network representing the Gibbs state. Its bond dimension $D$ --- limiting the amount of entanglement --- is a natural refinement parameter. Increasing $D$ we obtain a converged order parameter and its linear susceptibility close to the critical point. They confirm the existence of finite order parameter below the critical temperature $T_c$, provide a numerically exact estimate of~$T_c$, and give the critical exponents within $1\%$ of the 2D Ising universality class.

Figures

Tables

Table 1. Table 1: Fitted T c subscript 𝑇 𝑐 T_{c} and γ 𝛾 \gamma obtained for different symmetry breaking field values δ ​ h 𝛿 ℎ \delta h with D = 8 𝐷 8 D=8 and M = 72 𝑀 72 M=72 . Here data for 0.3566 < T < 0.3677 0.3566 𝑇 0.3677 0.3566<T<0.3677 were used. Changes of the fitted T c subscript 𝑇 𝑐 T_{c} and γ 𝛾 \gamma with decreasing δ ​ h ≤ 10 − 6 𝛿 ℎ superscript 10 6 \delta h\leq 10^{-6} are negligible as compared to their dependence on D 𝐷 D or range of data used to fit T c subscript 𝑇 𝑐 T_{c} and  γ 𝛾 \gamma .

	$\delta h$	$T_c$	$\gamma$
	${10}^{-6}$	0.356631	1.7324
	${10}^{-7}$	0.356633	1.7317
	${10}^{-8}$	0.356633	1.7317

Equations

$$H=-J\sum _j\sum _{\alpha =a,b}{\tau }_j^{\alpha }{\tau }_{j+e_{\alpha }}^{\alpha }.$$

$${\tau }_j^a=\frac{1}{4}\left(-{\sigma }_j^z+\sqrt{3}{\sigma }_j^x\right),{\tau }_j^b=\frac{1}{4}\left(-{\sigma }_j^z-\sqrt{3}{\sigma }_j^x\right).$$

$$m(T)\equiv ⟨{\sigma }_j^x⟩.$$

$$m(T)\propto (T_c-T{)}^{\beta },$$

$$\chi (T)={\frac{dm}{dh}|}_{h=0}.$$

$$\chi (T)\propto (T-T_c{)}^{\gamma },$$

$$T_c=0.3566\pm 0.0001,$$

$$T_c^{\mathrm{Ising}}=\frac{1}{2\log (1+\sqrt{2})}\approx \mathrm{0.567296.}$$

$$m(0)=\mathrm{0.993.}$$

$${\beta }_{\mathrm{Ising}}=\frac{1}{8},{\gamma }_{\mathrm{Ising}}=\frac{7}{4}.$$

$$\chi =\frac{\delta m}{\delta h},$$

Full text

[Published in: Physical Review B 96, 014420 (2017)]

I Introduction

Frustration in quantum spin systems occurs by competing exchange interactions and often leads to disordered spin liquids Bal10 ; Lucile . This is in contrast to Ising spins on a square lattice where periodically distributed partial frustration in form of exchange interactions with different signs does not suppress a phase transition at finite temperature $T_{c}$ Lon80 , while complete frustration gives a disordered classical phase Vil77 . Frustration may also be generated by a different mechanism — when Ising-like interactions for different pseudospin components compete on a square lattice in the two-dimensional (2D) compass model Nus05 ; Dou05 ; Dor05 ; Tro10 or on the honeycomb lattice in the Kitaev model Kit06 . While the short-range spin liquid is realized in the Kitaev model Bas07 , the pseudospin nematic order stabilizes below $T_{c}$ in the 2D compass model Wen10 ; Cza16 . In such cases entanglement plays an important role Ami08 and advanced methods of quantum many-body theory have to be applied.

In real systems pseudospin interactions concern the orbital degrees of freedom. The case of $e_{g}$ orbitals is paradigmatic here as it (i) is related to the 2D compass model Cin10 and (ii) initiated spin-orbital physics Kug82 ; Fei97 ; Tok00 ; Cor12 ; Brz15 — the well known systems with $e_{g}$ orbitals are: KCuF3 Bella ; Dei08 ; Pav08 , LaMnO3 Dag01 ; Fei99 ; Ish02 ; Ish03 ; Ole05 ; Pav10 ; Kov10 ; Sna16 , and LiNiO2 Rey01 ; Mil04 ; Rei05 . This field is very challenging due to the interplay and entanglement of spins and orbitals which leads to remarkable consequences Ole12 ; Brz12 . However, when spin order is ferromagnetic, as in the $(a,b)$ planes of KCuF3 and LaMnO3, spins disentangle and one is left with the 2D orbital $e_{g}$ model vdB99 ; Tan05 where hole propagation is possible by the coupling to orbitons vdB00 . Surprisingly, the tendency towards long-range order with such excitations is then opposite to that for spin systems Mat06 , i.e., $e_{g}$ orbital order occurs in a 2D square lattice below $T_{c}$ Ryn10 ; Wen11 , for instance in K2CuF4 kcuf4 ; Mos04 , while the role of quantum fluctuations increases with increasing dimension vdB99 ; Fei05 .

In this article we investigate a phase transition at $T_{c}$ in the 2D orbital $e_{g}$ model. A better understanding of the signatures of this phase transition provides a theoretical challenge. We present a very accurate estimate of $T_{c}$ and the critical exponents being in the 2D Ising universality class. These results could be achieved due to a remarkable recent progress in tensor networks due to the formulation of an algorithm at finite temperature using a projected entangled-pair operator (PEPO) var .

The paper is organized as follows. Sec. II gives brief overview of tensor network methods. Sec. III introduces simulated model. Sec. IV introduces 2D finite temperature tensor network method used to simulate the model. Numerical results are presented in Sec. V. Sec. VI summarizes the paper. Appendix A gives detailed description of results convergence analysis which enabled us to obtain trustworthy results for the model. Technical details of simulations are given in Appendix B. Finally Appendix C gives additional results for low temperature regime of the model.

II Tensor networks

Since the discovery of the density matrix renormalization group (DMRG) White ; Sch05 — that was later shown to optimize the matrix product state (MPS) variational ansatz Sch11 — quantum tensor networks proved to be an indispensable tool to study strongly correlated quantum systems Che15 . MPS ansatz was later generalized to a 2D projected entangled pair state (PEPS) PEPS ; CMRPEPS and supplemented with the multiscale entanglement renormalization ansatz (MERA) MERA . The networks do not suffer from the notorious sign problem fermions and in the doped case fermionic PEPS provided better variational energies for the $t$-$J$ model PEPStJ and the Hubbard model Cor15 than the best available variational Monte Carlo results. A combination of different tensor networks, supplemented with other sign-error free methods, seems to have finally settled the controversy on the ground state of the underdoped Hubbard model stripesHubbard . The networks — both MPS Yan11 ; Cin13 ; Ran15 and PEPS Poi12 ; Poi13 ; Wan13 — also made some major breakthroughs in the search for topological order. This is where, like in the $e_{g}$ model Ryn10 , geometric frustration often prohibits the traditional quantum Monte Carlo.

Thermal states of quantum Hamiltonians were explored much less than their ground states. In one dimension they can be represented by an MPS ansatz prepared with an accurate imaginary time evolution Ver04 ; WhiteT . A similar approach can be applied in 2D models Czarnik ; self , where the PEPS manifold is a compact representation for Gibbs states Molnar but the accurate evolution proved to be more challenging. Alternative direct contractions of the 3D partition function were proposed ChinaT but, due to local tensor update, they are expected to converge more slowly with increasing refinement parameter. Even a small improvement towards a full update can accelerate the convergence significantly HOSRG .

In order to avoid these problems, in the pioneering work var two of us introduced an algorithm to optimize variationally a projected entangled-pair operator (PEPO) representing the Gibbs state $e^{-\beta H}$ of a 2D lattice system ($\beta\equiv 1/T$). Its first challenging benchmark applications include the quantum compass Cza16 and Hubbard varHubbard models where it provided accuracy comparable to the best conventional methods.

It was not quite unexpected. Just like for the ground-state PEPS, the accuracy of the thermal PEPO is limited by its finite bond dimension $D$, i.e., the size of tensor indices connecting nearest-neighbor lattice sites. This size limits the entanglement within the ground/thermal state. However, by its very definition the Gibbs state is the mixed state that maximizes the entropy for a given average energy. Since this maximal entropy is actually the entropy of entanglement with the rest of the universe, then — thanks to the monogamy of entanglement — the Gibbs state also minimizes its internal entanglement. Among all states with the same average energy it is the one most suited to be represented by a tensor network. Encouraged by the benchmarks tests, in this work we apply the algorithm for the first time to a model that evades treatment by quantum Monte Carlo Ryn10 ; Wen11 . Numerical convergence and self-consistency alone allow us to make definitive statements on the physics of the model demonstrating the power of this method.

III The $e_{g}$ orbital model

The quantum $e_{g}$ model on an infinite square lattice is defined by the Hamiltonian

[TABLE]

Here $j$ labels lattice sites, $e_{a}(e_{b})$ are unit vectors along the $a(b)$ axis and $\tau^{\alpha}_{j}$ are orbital operators represented by Pauli matrices:

[TABLE]

The coupling in the orbital space depends on the spatial orientation of the bond. In what follows $J=1$.

At low temperature a spontaneous breaking of symmetry takes place and the system orders according to the strongest interaction $\propto\frac{3}{16}\sigma^{x}_{i}\sigma^{x}_{j}$ Cin10 . This symmetry breaking implies a finite real order parameter

[TABLE]

Unlike the 2D compass model Wen10 , the model (1) is not tractable by Monte Carlo Wen11 , but the order parameter suggests the 2D Ising universality class for the finite temperature transition which is confirmed by our simulations.

IV The algorithm at $T>0$

The algorithm was described in all technical detail elsewhere Cza16 . Its aim is to represent matrix elements of the operator $\rho=e^{-\beta H/2}$ by the 2D tensor network in Fig. 1. Here we show only a small $4\times 4$ unit of an infinite square lattice and each geometrical shape (here a green ball) represents a tensor. There is one tensor at every lattice site. Each line sticking out of the tensor represents one index. A (black) line connecting two tensors represents a tensor contraction through the connecting index. There is one bond index along every nearest neighbor bond. It has a finite bond dimension $D$. The dashed bond lines connect the $4\times 4$ unit with the rest of the lattice. The open (red) vertical indices number the orbital basis’ states. Those pointing up/down number bra/ket states. The desired 2D network in Fig. 1(f) — known as PEPO — can be contracted efficiently to obtain local expectation values. A finite $D$ is sufficient to represent Gibbs states with their limited entanglement.

On the other hand, the 2D operator $e^{-\beta H/2}\equiv U(\beta)$ can be naturally represented by a 3D network, the third dimension being the imaginary time $\beta$. The evolution is split into $N$ small time steps ($d\beta\ll 1$), $U(\beta)=U(d\beta)^{N}$. With a Suzuki-Trotter decomposition, each step can be represented by a 2D layer in Fig. 1(a). In the $e_{g}$ model, its bond indices have dimension $4$. The product of $N$ steps is the 3D network in Fig. 1(b). Here we show only three layers; the remaining $N-3$ ones are represented by the vertical dashed lines.

The 3D network is too hard to treat directly. Formally, it can be compressed to a 2D network by contracting along each vertical column first. The resulting 2D network in Fig. 1(c) arises at the price of a huge bond dimension $4^{N}$. Fortunately, we know that just a tiny $D$-dimensional subspace in the $4^{N}$ dimensions is enough to accommodate all correlations. Therefore, it is justified to insert every bond line with a $D$-dimensional projection made of two isometries. There are two independent projections along the axes $a$ and $b$, see Fig. 1(d). After the insertion, every isometry is absorbed into its nearest tensor truncating its bond index down to a tractable size $D$, see Fig. 1(e). The outcome is the desired PEPO $U(\beta)$ in Fig. 1(f), and the Gibbs state is $e^{-\beta H}=U^{\dagger}(\beta)U(\beta)$.

Now the problem is how to handle the huge isometries from $4^{N}$ to $D$. Fortunately, by a divide-and-conquer strategy, each of them can be split into a hierarchy of smaller isometries connected into a tree tensor network Cza16 . It is possible to optimize the smaller isometries one-by-one to obtain the most accurate projection available for a given $D$. The cost of the algorithm is polynomial in $D$ and only logarithmic in the number of steps $N$, allowing for $d\beta$ small enough to make the Suzuki-Trotter decomposition numerically exact at very little expense.

V Numerical Results

For each $T<T_{c}$ the order parameter $m$ (3) was converged in $D$ in the symmetry broken phase, see Fig. 2. For each $D$ it was fitted with a power law,

[TABLE]

see Fig. 3. Here $\beta$ is the order-parameter critical exponent (not to be confused with the inverse temperature $\beta=1/T$). For $D\geq 7$ the estimates: $0.35660<T_{c}<0.35664$, and $0.1258<\beta<0.1261$, do not depend significantly on increasing $D$. They slowly drift towards $T_{c}=0.35661$ and $\beta=0.125$, respectively. For more details see Appendix A.

In the symmetric phase above $T_{c}$, we calculated the magnetic susceptibility using the linear approximation,

[TABLE]

Here $h$ is an infinitesimal symmetry-breaking field $h\sum_{i}\tau_{i}^{x}$ added to the Hamiltonian (1). The derivative was approximated accurately by a finite difference between $h=10^{-6}$ and $h=0$. More details on $\chi(T)$ numerical calculation are given in Appendix B, see Fig. 7 and Table 1.

The susceptibility was converged in $D$ (Fig. 4) and fitted with a power law,

[TABLE]

see Fig. 3 and Appendix A. Again, for $D\geq 7$ the estimates: $0.35660<T_{c}<0.35665$, and $1.732<\gamma<1.740$, almost do not depend on increasing $D$, and drift towards $T_{c}=0.35661$ and $\gamma\simeq 1.75$. Altogether, both exponents are less than $1\%$ away from the exact $\beta=\frac{1}{8}$ [see Fig. 5(a)] and $\gamma=\frac{7}{4}$ in the 2D Ising universality class.

Remarkably, $T_{c}$ found from $m(T)$ (4) and $\chi(T)$ (6) is identical up to the four-digit precision. We propose

[TABLE]

deduced from the scatter of the data for $D\geq 7$ in Fig. 3(b) multiplied by a factor of $3$, see also Fig. 5(b). It is worthwhile to compare the above estimate (7) with the 2D Ising model Ons44 with interaction $\frac{1}{4}\sigma^{z}_{i}\sigma^{z}_{j}$,

[TABLE]

Exchange interactions in the dominating term $\frac{3}{16}\sigma^{x}_{i}\sigma^{x}_{j}$ in Eq. (1) are reduced by the factor $\frac{3}{4}$ from the 2D Ising model, so this reduction alone would give instead $T_{c}=0.75T_{c}^{\rm Ising}$. De facto, the obtained value in Eq. (7) is $T_{c}\simeq 0.6286T_{c}^{\rm Ising}$, i.e., it is further reduced by $\sim 16$% by quantum fluctuations activated at finite $T$ due to $\propto\frac{\sqrt{3}}{4}(\sigma^{x}_{i}\sigma^{z}_{j}+\sigma^{x}_{i}\sigma^{z}_{j})$ and $\propto\frac{1}{4}\sigma^{z}_{i}\sigma^{z}_{j}$ terms in Eq. (1). The order parameter (3) at $T=0$ is almost saturated as quantum fluctuations are negligible at $T\to 0$,

[TABLE]

More details on $m(0)$ simulation are given in Appendix C, see Fig. 8. The value in Eq. (9) was obtained by the present method and agrees with the ground state MERA calculations Cin10 . This shows that the quantum fluctuation effects in the $e_{g}$ orbital model (1) are very weak indeed at $T=0$ vdB99 , while at $T>0$ the fluctuations are activated and reduce significantly the value of the critical temperature down to $T_{c}\simeq 0.3566$, see Eq. (7). Indeed quantum fluctuations play a role here but are not as significant as for the 2D SU(2) symmetric Heisenberg antiferromagnet Mat06 . Yet, the entanglement between the orbital operators is here much reduced from that in the 2D compass model var and therefore such an accurate estimate of $T_{c}$ (7) is possible.

VI Summary

Being a paradigmatic frustrated system, the orbital $e_{g}$ model evades treatment by quantum Monte Carlo but it proves to be accurately tractable by our thermal tensor network. The notorious sign problem — often inescapable for quantum Monte Carlo — is not an issue for our method. Instead the relevant issue is if the entanglement in a thermal state can be accommodated within a bond dimension that is small enough to fit into a classical computer. This criterion is satisfied by the thermal state of the $e_{g}$ model and a four-digit estimate of the critical temperature and a better than $1\%$ accuracy of the critical exponents could be achieved. Since the Gibbs state is the least entangled one among all excited states with the same average energy, it is potentially the easiest target for a suitable tensor network.

Acknowledgements.

We thank Philippe Corboz for insightful discussions. We kindly acknowledge support by Narodowe Centrum Nauki (NCN, National Science Centre, Poland) under Projects: No. 2013/09/B/ST3/01603 (P.C. and J.D.) and No. 2016/23/B/ST3/00839 (A.M.O.). The work of P.C. on his Ph.D. thesis was supported by NCN under Project No. 2015/16/T/ST3/00502.

Appendix A Convergence of the results

The bond dimension $D$ (see Fig. 1) has to be large enough to accommodate the entanglement in the thermal state. Furthermore, an environmental bond dimension $M$ that is used in the analysis of the effective 2D tensor network depicted in Fig. 1(f) (see Ref. Cza16 for details) has to be large enough to accommodate long range correlations. In general, these requirements cannot be satisfied at the critical temperature $T_{c}$ but the phase transition can be approached from both sides close enough to fit the critical power laws. In this appendix we demonstrate that indeed we are able to approach $T_{c}$ close enough to obtain stable and converged fits.

All results presented here, which were obtained with $M=72$, are converged in $M$. Another potential source of errors are Trotter errors. They are not a significant issue for our approach as its cost scales at most logarithmically with the the inverse Trotter time step $1/d\beta$. Our results were obtained with $d\beta\leq 0.001$ and are converged in $d\beta$.

The convergence of the critical exponents, $\beta$ for the magnetization $m(T)$ and $\gamma$ for the susceptibility $\chi(T)$, is shown in Figs. 5(a) and 6(a) where we compare them with the 2D Ising model exponents,

[TABLE]

For $D\geq 7$ we see that the exponents approach the Ising values while $T_{lim}$ is approaching $T_{c}$. For $T_{lim}$ sufficiently close to $T_{c}$ they no longer depend significantly on range of $T$ depending instead primarily on $D$. In this regime all fitted exponents fall within $1\%$ of 2D Ising universality class, drifting towards $\beta_{Ising}$ or $\gamma_{\rm Ising}$ with increasing $D$. The obtained behavior of the exponents indicates the 2D Ising universality class of the transition.

The data collected in Figs. 5(b) and 6(b) demonstrate similar convergence behavior of fitted $T_{c}$ as for the exponents. For $D\geq 7$ fitted $T_{c}$ approaches $T_{c}=0.3566$ when $T_{lim}$ is approaching the critical point. For $T_{lim}$ sufficiently close to $T_{c}$ the critical point $T_{c}$ begins to depend primarily on $D$ rather than on $T_{lim}$. Reaching this regime where the fits become stable with respect to $T_{lim}$ justifies taking into account only their $D$ dependence to obtain the final $T_{c}$ estimate Eq. (7).

We remark that our estimate of $T_{c}$ is based on two independent $T_{c}$ estimates, coming either from the $\chi(T)$ or $m(T)$ fits, which agree up to five digits for the largest $D$.

Appendix B Numerical details

In our simulations we use the algorithm described in detail in Ref. Cza16 . In prticular we use corner matrix renormalization (CMR) to contract approximately tensor networks representing thermal states CMR ; CMRPEPS . To reach convergence of the observables $m$ and $\chi$ approximately $10$ iterations of the optimization loop were necessary. The isometries at the beginning of the loop were initialized by a local truncation scheme based on higher-order singular value decomposition. The CMR procedure made $\sim 1000$ iterations in the whole loop. The further away from the phase transition, the fewer CMR iterations were necessary to reach convergence.

Linear susceptibility $\chi(T)$ defined by Eq. (5) was calculated from a finite difference of the order parameter $\delta m$ corresponding to finite difference of the symmetry breaking field $\delta h=10^{-6}$:

[TABLE]

where $\delta m=m(h=\delta h)-m(h=0)$. Fig. 7 shows that $\chi(T)$ is already converged in $\delta h$ for $\delta h=10^{-6}$. More accurate benchmark of $\delta h$ convergence is given by Table 1 showing that decreasing $\delta h$ further results in changes of fitted $\gamma$ and $T_{c}$ that are negligible as compared to their dependence on $D$ or the range of $T$.

All simulations were done in Matlab with an extensive use of the Ncon procedure encon . To give an idea of the actual time and computer resources needed to generate the data, the most challenging data points nearest to the phase transition, with the largest bond dimensions $D=11$ and $M=72$, required $1-2$ days on a desktop.

Appendix C Simulation of the low temperature phase

The entanglement in the low $T$ phase is small enough to converge the curve $m(T)$ in $D$ already for $D=4$, see Fig. 8.

Thanks to a short correlation length at low temperature, the calculations are much less demanding numerically than close to the critical point. Because of that we were able to generate the data shown in Fig. 8 during one day using a laptop.