Tensor network simulation of the Kitaev-Heisenberg model at finite temperature

TL;DR

This paper employs tensor network methods to study the finite-temperature properties of the Kitaev-Heisenberg model, estimating critical temperatures and benchmarking against Monte Carlo results, thus advancing tensor network applications in complex quantum systems.

Contribution

It introduces an efficient tensor network approach for finite-temperature simulation of the Kitaev-Heisenberg model on a hexagonal lattice, including a novel lattice mapping and benchmarking against Monte Carlo data.

Findings

Critical temperatures are an order of magnitude lower than coupling constants.

The tensor network method accurately captures spin ordering crossover.

The approach is validated against Monte Carlo results for the pure Kitaev model.

Abstract

We investigate the Kitaev-Heisenberg (KH) model at finite temperature using the exact environment full update (eeFU), introduced in Phys. Rev. B 99, 035115 (2019), which represents the purification of a thermal density matrix on an infinite hexagonal lattice by an infinite projected entangled pair state (iPEPS). We show that thanks to a dynamical mapping from a hexagonal to a rhombic lattice, the eeFU on the hexagonal lattice is as efficient as the simple full update (FU) algorithm. Critical temperatures for coupling constants in the stripy and the antiferromagnetic phase are estimated. They are an order of magnitude less than the couplings in the Hamiltonian. By a duality transformation, these results can be mapped to, respectively, the ferromagnetic and zigzag phases. For the special case of the pure Kitaev model, which is tractable by quantum Monte-Carlo but the most challenging for tensor networks, the algorithm is benchmarked against the Monte-Carlo results. It recovers accurately the crossover to spin ordering and qualitatively the one to flux ordering.

Tables

Table 1. Table 1: The critical temperature T c subscript 𝑇 𝑐 T_{c} and the exponent 1 / β ~ ​ δ 1 ~ 𝛽 𝛿 1/\tilde{\beta}\delta at ϕ = − 63 ∘ italic-ϕ superscript 63 \phi=-63^{\circ} obtained by fitting T ∗ ​ ( h ) superscript 𝑇 ℎ T^{*}(h) with the scaling ansatz ( 12 ), see Fig. 6 . Here we gather results obtained by the two methods of estimating T ∗ ​ ( h ) superscript 𝑇 ℎ T^{*}(h) , using to that end location of peaks either in m ′ ​ ( t , h ) superscript 𝑚 ′ 𝑡 ℎ m^{\prime}(t,h) or C V ​ ( t , h ) subscript 𝐶 𝑉 𝑡 ℎ C_{V}(t,h) . The errors are statistical uncertainties of the fits for a confidence level 95 % percent 95 95\% . The results for different choices of χ 𝜒 \chi and d ​ β 𝑑 𝛽 d\beta suggests that the modest χ = 3 ​ D 𝜒 3 𝐷 \chi=3D and d ​ β = 0.01 𝑑 𝛽 0.01 d\beta=0.01 are sufficient to obtain converged estimates of T c subscript 𝑇 𝑐 T_{c} and 1 / β ~ ​ δ 1 ~ 𝛽 𝛿 1/\tilde{\beta}\delta . Furthermore, both methods of T ∗ superscript 𝑇 T^{*} -estimation give mutually consistent results. Combining all the results we finally estimate T c = 0.057 ​ ( 2 ) subscript 𝑇 𝑐 0.057 2 T_{c}=0.057(2) , 1 / β ​ δ = 0.46 ​ ( 5 ) 1 𝛽 𝛿 0.46 5 1/\beta\delta=0.46(5) . To obtain the final estimates we use the weighted average, taking the error as the statistical error of the average for a confidence level 99.7 % percent 99.7 99.7\% .

the method	$D$	$\chi$	d$\beta$	$T_c$	$1/\tilde{\beta }\delta$
$m^{\prime }$ peaks	$6$	$42$	$0.01$	$0.058\left(2\right)$	$0.46\left(4\right)$
$m^{\prime }$ peaks	$7$	$21$	$0.01$	$0.059\left(14\right)$	$0.53\left(36\right)$
$m^{\prime }$ peaks	$8$	$24$	$0.01$	$0.057\left(9\right)$	$0.46\left(18\right)$
$m^{\prime }$ peaks	$6$	$24$	$0.005$	$0.059\left(2\right)$	$0.48\left(4\right)$
$m^{\prime }$ peaks	$6$	$24$	$0.02$	$0.058\left(4\right)$	$0.45\left(9\right)$
$m^{\prime }$ peaks	$6$	$18$	$0.01$	$0.058\left(8\right)$	$0.46\left(16\right)$
$m^{\prime }$ peaks	$6$	$30$	$0.01$	$0.058\left(4\right)$	$0.45\left(9\right)$
$C_V$ peaks	$6$	$42$	$0.01$	$0.055\left(3\right)$	$0.48\left(8\right)$
$C_V$ peaks	$7$	$21$	$0.01$	$0.055\left(11\right)$	$0.54\left(36\right)$
$C_V$ peaks	$8$	$24$	$0.01$	$0.055\left(11\right)$	$0.48\left(30\right)$

Table 2. Table 2: The exponents ν / β ~ ​ δ 𝜈 ~ 𝛽 𝛿 \nu/\tilde{\beta}\delta and ( β ~ − 1 ) / β ~ ​ δ ~ 𝛽 1 ~ 𝛽 𝛿 (\tilde{\beta}-1)/\tilde{\beta}\delta obtained by fitting ξ ​ ( T ∗ , h ) 𝜉 superscript 𝑇 ℎ \xi(T^{*},h) and m ′ ​ ( T ∗ , h ) superscript 𝑚 ′ superscript 𝑇 ℎ m^{\prime}(T^{*},h) in Fig. 7 .

$D$	the exponent
$6$	$\nu /\tilde{\beta }\delta =0.25\left(1\right)$
$7$	$\nu /\tilde{\beta }\delta =0.33\left(5\right)$
$8$	$\nu /\tilde{\beta }\delta =0.27\left(2\right)$
$6$	$\left(\tilde{\beta }-1\right)/\tilde{\beta }\delta =0.423\left(8\right)$
$7$	$\left(\tilde{\beta }-1\right)/\tilde{\beta }\delta =0.40\left(2\right)$
$8$	$\left(\tilde{\beta }-1\right)/\tilde{\beta }\delta =0.38\left(2\right)$

Table 3. Table 3: The critical temperatures T c subscript 𝑇 𝑐 T_{c} and the critical exponents 1 / β ~ ​ δ 1 ~ 𝛽 𝛿 1/\tilde{\beta}\delta obtained for ϕ = − 17 ∘ italic-ϕ superscript 17 \phi=-17^{\circ} in Figs. 9 . The T c subscript 𝑇 𝑐 T_{c} estimates agree within the uncertainties.

method	$D$	$T_c$	$1/\tilde{\beta }\delta$
$m^{\prime }$ peaks	$6$	$0.071\left(4\right)$	$0.31\left(3\right)$
$C_V$ peaks	$6$	$0.066\left(4\right)$	$0.35\left(4\right)$
$m^{\prime }$ peaks	$7$	$0.072\left(6\right)$	$0.29\left(4\right)$
$C_V$ peaks	$7$	$0.070\left(4\right)$	$0.38\left(4\right)$
$m^{\prime }$ peaks	$8$	$0.073\left(4\right)$	$0.32\left(2\right)$
$C_V$ peaks	$8$	$0.070\left(4\right)$	$0.38\left(4\right)$

Equations

$$H_{\mathrm{KH}}=\sum _{<i,j>}{H}_{i,j}^{(\gamma )},$$

$$H_{i,j}^{(\gamma )}=J{\mathbf{S}}_i\cdot {\mathbf{S}}_j+K{S}_i^{\gamma }{S}_j^{\gamma }$$

$$J=A\text{cos}\varphi ,K=2A\text{sin}\varphi .$$

$$\tan \tilde{\varphi }=-\tan \varphi -1,$$

$$H=H_{\mathrm{KH}}-\sum _i{h}_i{S}_i^z,$$

$${\mathrm{Tr}}_a|\psi (\beta )⟩⟨\psi (\beta )|\propto e^{-\beta H}.$$

$$e^{-d\beta H/2}\approx$$

$$G_x(d\beta /2)G_y(d\beta /2)G_z(d\beta )G_y(d\beta /2)G_x(d\beta /2),$$

$$m(t,h)$$

$$m^{\prime }(t,h)$$

$$\xi (t,h)$$

$$C_V(t,h)$$

$$T^{\ast }(h)=T_c+x^{\ast }{h}^{1/\tilde{\beta }\delta }.$$

$$m^{\prime }(t^{\ast },h)\propto h^{(\tilde{\beta }-1)/\tilde{\beta }\delta },$$

$$\xi (t^{\ast },h)\propto h^{-\nu /\tilde{\beta }\delta }.$$

$${\tilde{H}}_{i,j}^{\gamma }=\tilde{A}\left[\cos (\tilde{\varphi }){\mathbf{S}}_i\cdot {\mathbf{S}}_j+2\sin (\tilde{\varphi })S_i^{(\gamma )}{S}_j^{(\gamma )}\right].$$

$$T_c=0.057(2),1/\tilde{\beta }\delta =0.46(5),$$

$$T_c=0.091(4),1/\tilde{\beta }\delta =0.46(5).$$

$$\nu /\tilde{\beta }\delta =0.41(1).$$

$$(\tilde{\beta }-1)/\tilde{\beta }\delta =0.26(1).$$

$$T_c=0.072(3),1/\tilde{\beta }\delta =0.35(2).$$

$$T_c=0.062(3),1/\tilde{\beta }\delta =0.35(2).$$

Full text

I Introduction

Weakly entangled quantum states can be efficiently represented by tensor networks Verstraete et al. (2008); Orús (2014): either a 1D matrix product state (MPS) Fannes et al. (1992), its 2D generalization to a projected entangled pair state (PEPS) Verstraete and Cirac (2004), or a multi-scale entanglement renormalization ansatz (MERA) Vidal (2007, 2008); Evenbly and Vidal (2014a, b). The MPS is a compact representation of ground states of 1D gapped local Hamiltonians Verstraete et al. (2008); Hastings (2007); Schuch et al. (2008) and purifications of their thermal states Barthel (2017). It is also the ansatz underlying the powerful density matrix renormalization group (DMRG) White (1992, 1993); Schollwöck (2005); Schöllwock (2011). Analogously, the 2D PEPS is expected to represent ground states of 2D gapped local Hamiltonians Verstraete et al. (2008); Orús (2014) and their thermal states Wolf et al. (2008); Molnar et al. (2015), though representability of area-law states in general was shown to have its limitations Ge and Eisert (2016). Tensor networks evade the sign problem plaguing the quantum Monte Carlo, hence they can deal with fermionic systems Corboz et al. (2010a); Pineda et al. (2010); Corboz and Vidal (2009); Barthel et al. (2009); Gu et al. (2010), as was shown for both finite Kraus et al. (2010) and infinite PEPS Corboz et al. (2010b, 2011).

The PEPS was proposed originally for ground states of finite systems Verstraete and Cirac (2004); Murg et al. (2007) generalizing earlier attempts to construct trial wave-functions for specific models Nishio et al. (2004). Efficient numerical methods for infinite PEPS (iPEPS) Jordan et al. (2008); Jiang et al. (2008); Gu et al. (2008); Orús and Vidal (2009) promoted it to a versatile tool for strongly correlated systems in 2D. Examples of their potential include a solution of the long standing magnetization plateaus problem in the highly frustrated compound $\textrm{SrCu}_{2}(\textrm{BO}_{3})_{2}$ Matsuda et al. (2013); Corboz and Mila (2014), demonstration of the striped nature of the ground state of the doped 2D Hubbard model Zheng et al. (2017), and new evidence supporting gapless spin liquid (SL) in the kagome Heisenberg antiferromagnet Liao et al. (2017). Recent developments in iPEPS optimization Phien et al. (2015); Corboz (2016a); Vanderstraeten et al. (2016), contraction Fishman et al. (2018); Xie et al. (2017), energy extrapolations Corboz (2016b), and universality class estimation Corboz et al. (2018); Rader and Läuchli (2018); Rams et al. (2018) open possibility of applying it to even more difficult problems, including simulation of thermal states Czarnik et al. (2012); Czarnik and Dziarmaga (2014, 2015a); Czarnik et al. (2016a); Czarnik and Dziarmaga (2015b); Czarnik et al. (2016b, 2017); Dai et al. (2017); Czarnik et al. (2019); Czarnik and Corboz (2019); Kshetrimayum et al. (2019), mixed states of open systems Kshetrimayum et al. (2017); Czarnik et al. (2019), exited states Vanderstraeten et al. (2015), or real time evolution of 2D quantum states Czarnik et al. (2019); Hubig and Cirac (2019) .

In parallel with iPEPS, progress was made in simulating systems on cylinders of finite circumference with DMRG. This method of high numerical stability is routinely used to investigate 2D ground states Zheng et al. (2017); Cincio and Vidal (2013) and recently was applied also to thermal states on a cylinder Bruognolo et al. (2017); Chen et al. (2018a, 2019); Li et al. (2019) but the exponential growth of the bond dimension limits the circumference to a few lattice sites. Among alternative approaches are methods of direct contraction and renormalization of a 3D tensor network representing a 2D thermal density matrix Li et al. (2011); Xie et al. (2012); Ran et al. (2012, 2013, 2018); Peng et al. (2017); Chen et al. (2018b); Ran et al. (2019).

The Kitaev model Kitaev (2006) is an exactly solvable pseudospin-$1/2$ system on a hexagonal lattice with Ising-like couplings of $\gamma=x,y,z$ components of nearest-neighbor pseudospins with strength $K_{\gamma}$ along $\gamma$-bonds, see Fig. 1. In any of its three A-phases, when one of the three couplings dominates, the model reduces to the effective toric code Hamiltonian Kitaev (2003). On the other hand, in the B-phase, where the three couplings are of similar strength, the ground state of the highly frustrated model is a critical quantum spin liquid (SL). In the SL a magnetic field opens a finite energy gap protecting a chiral topological order. Its non-abelian anyonic excitations can be employed to perform universal topological quantum computation Nayak et al. (2008). This motivates intensive search for a robust physical implementation of the model.

In spin systems their $SU(2)$ symmetry constrains the interaction to be of the Heisenberg type. In order to break the symmetry – and introduce the bond-anisotropy at the same time – the spins can be mixed with orbital degrees of freedom, as originally argued for iridium oxides Chaloupka et al. (2010). The resulting Kitaev-Heisenberg (KH) model was considered as a minimal model to study stability of the Kitaev spin liquid phase in materials like $\textrm{Na}_{2}\textrm{IrO}_{3}$, $\alpha\textrm{-Li}_{2}\textrm{IrO}_{3}$, $\textrm{Li}_{2}\textrm{RhO}_{3}$, and $\alpha\textrm{-RuCl}_{3}$, though recent results suggest that more general extensions of the KH model are requiredWinter et al. (2017); Rusnačko et al. (2019). In this paper we follow Refs. Chaloupka et al., 2010; Reuther et al., 2011; Chaloupka et al., 2013; Price and Perkins, 2012; Nasu et al., 2015 and investigate the basic KH model at finite temperature as a first step towards its more realistic extensions.

The paper is organized as follows. In section II we review the KH model and its phase diagram at zero temperature. In section III we briefly outline the tensor network method to simulate thermal states of a quantum Hamiltonian. Here we introduce the dynamical mapping from a hexagonal to rhombic lattice that renders the exact environment full update (eeFU) as efficient as the simplified full update (FU) algorithm. In section IV a scaling theory is discussed that is necessary to extrapolate results obtained with a finite symmetry-breaking bias to zero bias field. In section V we present results for the striped and antiferromagnetic phases of the KH model. In section V.3 we benchmark our method against quantum Monte-Carlo resultsNasu et al. (2015) in the pure Kitaev model with a quantum spin-liquid ground state. We conclude in section VI.

II Model

The model Chaloupka et al. (2010, 2013); Gotfryd et al. (2017) is a sum of nearest-neighbor terms on a hexagonal lattice,

[TABLE]

where

[TABLE]

depends on bond direction $\gamma=x,y,z$, see Fig. 1. Here $\mathbf{S}=(S^{x},S^{y},S^{z})=\frac{1}{2}(\sigma_{x},\sigma_{y},\sigma_{z})$ are spin-$1/2$ operators defined by Pauli matrices. The coupling constants are parameterized by an angle $\phi$:

[TABLE]

Here $A>0$ is a constant.

The zero-temperature phase diagram is shown in Fig. 2. It was obtained by a variety of methods Chaloupka et al. (2010, 2013); Gotfryd et al. (2017) and corroborated by iPEPS Osorio Iregui et al. (2014). Pairs of angles, $\phi$ and $\tilde{\phi}$ satisfying

[TABLE]

on the right and left of the diagram, respectively, are related by a duality transformation Chaloupka et al. (2010). There are two self-dual points, $\phi=\tilde{\phi}=\pm 90^{\circ}$, where the model reduces to the Kitaev model. Each of these two Kitaev points is surrounded by a gapless quantum spin liquid (SL). The same duality maps the antiferromagnetic (AF) and stripy phases on the right to, respectively, the zigzag and ferromagnetic (FM) phases on the left.

For $\phi=0^{\circ}$ and $180^{\circ}$ the model reduces to the antiferromagnetic and ferromagnetic Heisenberg model, respectively. By the Mermin-Wagner theorem, its $\textrm{SU}(2)$-symmetry prevents spontaneous symmetry breaking at any $T>0$. The duality transformation maps these two points to $\phi=135^{\circ}$ and $-45^{\circ}$, respectively. Their hidden $SU(2)$ symmetry also prevents the ordering at finite $T$.

The frustrated model is not tractable by quantum Monte Carlo, except the pure Kitaev model Nasu et al. (2015). Its mean-field theory is $SU(2)$-symmetric Chaloupka et al. (2013); Gotfryd et al. (2017) suggesting no finite-$T$ ordering at any $\phi$ but a spin-wave expansion and plaquette mean-field suggest a disorder-induced-order at low temperatures stabilized by both quantum and thermal fluctuations Chaloupka et al. (2013); Reuther et al. (2011). The latter effect is confirmed by classical Monte Carlo simulations Price and Perkins (2012, 2013). The model is also tractable by a high-$T$ series expansion Singh and Oitmaa (2017).

In this work we treat the finite-$T$ KH model with a quantum tensor network for the first time. Previously we used quantum tensor networks to simulate the closely related compass and $e_{g}$ models Czarnik et al. (2016a, 2017) at finite $T$ achieving good accuracy. In order to simulate the model in neighbourhood of non-analytic critical points efficiently, we add a tiny symmetry breaking bias,

[TABLE]

with a magnitude $h=|h_{i}|$. $h_{i}=h$ in the FM phase and is staggered in the AF phase. To obtain critical properties of the Kitaev-Heisenberg model we extrapolate to $h=0$ as described in Sec. IV.

III Tensor network

In this work we apply the exact environment full update (eeFU) introduced and benchmarked in Ref. Czarnik et al., 2019. Here we just outline the algorithm emphasizing its adjustments to the KH model referring for more details to Ref. Czarnik et al., 2019. The most important development is the dynamical mapping from a hexagonal to rhombic lattice that makes the eeFU as efficient as the FU algorithm.

Thanks to the duality transformation (4), it is enough to consider the AF and FM phases only. They require only two sublattices: $A$ and $B$. We enlarge the Hilbert space by accompanying every pseudospin-$1/2$ with a pseudospin-$1/2$ ancilla. The iPEPS tensor network in Fig. 3(a) represents a thermal state’s purification $|\psi(\beta)\rangle$ in the enlarged space. Here $\beta=1/T$ is an inverse temperature. Its partial trace over the ancillas ($a$) yields the Gibbs state as a thermal density matrix:

[TABLE]

The purification is evolved in the imaginary time $\beta$ with the eeFU algorithm: $|\psi(\beta)\rangle=e^{-\beta H/2}|\psi(0)\rangle$.

The time evolution is represented by a product of $N$ small time steps $e^{-\beta H/2}=\left(e^{-d\beta H/2}\right)^{N}$, where $N=\beta/d\beta$. Each time step is subject to a second order Suzuki-Trotter decomposition Suzuki (1966, 1976); Trotter (1959):

[TABLE]

where $G_{\gamma}(d\beta)=\prod_{\langle i,j\rangle\|\gamma}e^{-\frac{1}{2}d\beta H^{(\gamma)}_{ij}}$ is a product of nearest neighbor gates over all $\gamma$-bonds. Here $H^{(z)}_{ij}$ includes also the bias in Eq. (5).

The action of $G_{\gamma}$ on one of the $\gamma$-bonds is shown in Fig. 3(b). A contraction of the “old” tensors $A,B$ with the gate $e^{-\frac{1}{2}d\beta H^{(\gamma)}_{ij}}$ becomes a contraction of new exact tensors $A^{\prime}$, $B^{\prime}$ with an enlarged $\gamma$-bond dimension $4D$. The exact contraction $A^{\prime}-B^{\prime}$ is approximated by a contraction of new tensors $A^{\prime\prime}-B^{\prime\prime}$ with the original bond dimension $D$. The new tensors are optimized to minimize the error introduced by this approximation to the whole purification.

In order to minimize the error of the infinite purification, we need a tensor environment for the considered bond $\gamma$. To this end we treat the exact two-site contraction $A^{\prime}-B^{\prime}$ as if it were a single iPEPS tensor on the two sites, see Fig. 4. Effectively, every two nearest neighbor sites connected by a $\gamma$-bond are fused into a single site. The hexagonal lattice is replaced by a rhombic one, which can be treated as a square lattice. This way we can employ full potential of the robust square-lattice corner transfer matrix renormalization group Baxter (1978); Nishino and Okunishi (1996); Orús and Vidal (2009); Corboz et al. (2014) (CTMRG) to obtain the environment for the $\gamma$-bond, see Fig. 4(c). However, the main advantage is that every enlarged $4D$-dimensional $\gamma$-bond index is hidden inside the square-lattice composite iPEPS tensor and, hence, it does not slow down the CTMRG which is the main bottleneck of the whole algorithm.

IV Estimation of critical temperature

The evolution near a critical point is challenging Czarnik et al. (2012, 2019). In particular, finite $\chi$ limits the correlation length which can be obtained by CTMRG Nishino et al. (1996), hence a large $\chi$ is necessary to render the environment of the $\gamma$-bond accurate enough to obtain correct new tensors $A^{\prime\prime}$ and $B^{\prime\prime}$. Therefore, in Refs. Czarnik et al., 2012, 2019 a small symmetry-breaking bias $h$ was introduced to turn the transition into a smooth crossover making the correlation length finite and allowing for results well converged in $\chi$. However, in order to estimate $T_{c}$ an extrapolation back to $h=0$ was necessary. To this end, a systematic scaling theory was used Czarnik et al. (2019) yielding very accurate results for the quantum Ising model. Here we follow the same approach.

According to the scaling theory the order parameter $m$, its derivative with respect to $T$, and the correlation length $\xi$ satisfy the scaling laws:

[TABLE]

respectively. Here $t=T-T_{c}$, the prime is a derivative with respect to $t$, $f(x)$, $g(x)$ and $h(x)$ are non-universal scaling functions, while $\tilde{\beta},\delta,\nu,\alpha$ are universal critical exponents. In order to estimate $T_{c}$ we use an observation that, for a fixed $h$, the slope $m^{\prime}(t,h)$ has a peak at $t^{*}=T^{*}-T_{c}>0$. In the regime of small $h$ its position, determined by the maximum $x^{*}$ of the scaling function $f^{\prime}(x)$, should scale as

[TABLE]

Fitting numerical data for the pseudo-critical temperature, $T^{*}(h)$, with the function on the right hand side we estimate three parameters: $x^{*}$, $1/\tilde{\beta}\delta$ and, most importantly, $T_{c}$. Similarly we observe that also $C_{V}(t,h)$ has a maximum at $\tilde{t}^{*}=\tilde{T}^{*}-T_{c}>0$, which position scales the same as $T^{*}(h)$ in (12) making it also suitable to estimate $1/\tilde{\beta}\delta$ and $T_{c}$.

Furthermore, we use the behavior of $\xi(t^{*},h)$ and $m^{\prime}(t^{*},h)$ to test self-consistency of the scaling theory. We observe that $\xi(t^{*},h)$ is close to the maximal correlation length for a given bias $h$ and $m^{\prime}(t^{*},h)$ is the maximal magnetization’s slope by definition. Equations (9,10,12) imply two power laws that do not depend on the unknown $T_{c}$:

[TABLE]

[TABLE]

Therefore, they provide a reliable test whether $h$ is small enough to achieve the critical scaling regime.

V Results

We choose to study two angles: $\phi=-63^{\circ}$ and $-17^{\circ}$. The former sits midway between the zero-temperature phase boundary at $-80^{\circ}$, separating the stripy phase from the spin liquid, and the $\textrm{SU}(2)$-symmetric point at $-45^{\circ}$. Likewise, the latter sits midway between the stripy-AF phase transition and the $SU(2)$-symmetric Heisenberg point at $0^{\circ}$. This is why we expect a relatively high critical temperature at both angles. Furthermore, $\phi=-63^{\circ}$ lies near the range $J/K=-0.3...-0.1$ reported recently Winter et al. (2017) for a proximate Kitaev spin liquid material $\alpha\textrm{-RuCl}_{3}$, making it a good starting point for a future study of more realistic extensions of the minimal KH model.

The duality transformation (4) maps the results for $\phi=-63^{\circ}$ and $-17^{\circ}$ to $\phi=-136.09^{\circ}$ in the ferromagnetic phase and $145.23^{\circ}$ in the zig-zag phase, respectively.

V.1 Stripy (ferromagnetic) phase:

$\phi=-63^{\circ}$ ($\phi=-136.09^{\circ})$

The duality transformation (4) maps $\phi=-63^{\circ}$ in the stripy phase to $\tilde{\phi}=-136.09^{\circ}$ in the ferromagnetic phase where we actually make simulations taking advantage of the fact that we need only two sublattices there. After the transformation the nearest neighbor terms in the Hamiltonian become

[TABLE]

where $\tilde{A}=-A\cos(\phi)/\cos(\tilde{\phi})$. $\tilde{A}\approx 0.63$ when we set $A=1$. The order parameter for the stripy phase equals the ferromagnetic order parameter of the transformed Hamiltonian $m=2\sqrt{\langle S^{x}\rangle^{2}+\langle S^{y}\rangle^{2}+\langle S^{z}\rangle^{2}}.$

The specific heat for different values of the small bias $h$ is shown in Fig. 5. At higher temperatures, comparable to the coupling constants in the Hamiltonian, we obtain a broad peak which does not depend on the bias. At an order of magnitude lower temperatures we can see sharp peaks whose position and magnitude depend strongly on the applied bias. This sensitivity suggests that they indicate spontaneous symmetry breaking in the direction of the bias. Below we analyze this low temperature regime in detail employing numerical data obtained with bond dimensions $D=6,7,8$ for biases in the range $6.1\times 10^{-4}\leq h\leq 8.5\times 10^{-3}$. This is where the results appear converged in $D$.

Fig. 6(a,b) shows the pseudo-critical temperatures $T^{*}(h)$ obtained from the peaks of $m^{\prime}$. The data are fitted accurately by the scaling ansatz (12), indicating a second-order phase transition. Fig. 6(c) shows the same for the pseudo-critical temperatures $\tilde{T}^{*}(h)$ obtained from the peaks of $C_{V}$. The results are again fitted accurately by the scaling ansatz (12). In Tab. 1 we collect $T_{c}$ and $1/\tilde{\beta}\delta$ fitted in Fig. 1. We find that the results for different $D$ are mutually consistent. We remark that the results obtained with $D=7,8$ are more ”noisy” than $D=6$, which may be related to stability issues in the full updatePhien et al. (2015); Hubig and Cirac (2019); Hasik and Becca (2019). Furthermore, the results obtained from the $m^{\prime}$ and $C_{V}$ peaks are also mutually consistent. Finally, for $\phi=-63^{\circ}$ we obtain

[TABLE]

see details in Tab. 1. Using the duality transformation, we obtain for $\phi=-136.09^{\circ}$

[TABLE]

As a further self-consistency check, we analyze the correlation length $\xi(t^{*},h)$ at $T^{*}$. We extract $\xi$ from the iPEPS with the precise method in Ref. Rams et al., 2018. Figure 7(a) shows a log-log plot of $\xi(t^{*},h)$ in function of $h$, which is fitted well by a linear behaviour predicted by the scaling ansatz (14). The fits for different $D$ yield close to each other values of the exponent $\nu/\tilde{\beta}\delta$, see Tab. 2. A weighted average combining the results for different $D$ yields

[TABLE]

Notice that the largest $\xi(t^{*},h)$ of $3.5$ lattice constants of the rhombic lattice obtained in Fig. 4(b) is beyond reach of the state of the art finite cluster exact diagonalization or DMRG on a cylinder. We analyze also $m^{\prime}(t^{*},h)$, see Fig. 7(b). Again we find that the results can be accurately fitted by the scaling ansatz (13). The estimates of $(\tilde{\beta}-1)/\tilde{\beta}\delta$ obtained with different $D$ are close to each other, see Tab. 2. Their combination yields

[TABLE]

V.2 Antiferromagnetic (zig-zag) phase:

$\phi=-17^{\circ}$ ($\phi=145.23^{\circ}$)

The order parameter is the staggered magnetization $m=\sqrt{\langle S^{x}_{A}-S^{x}_{B}\rangle^{2}+\langle S^{y}_{A}-S^{y}_{B}\rangle^{2}+\langle S^{z}_{A}-S^{z}_{B}\rangle^{2}}$. In the following we present results for the bond dimensions $D=6,7,8$ and the time step $d\beta=0.01$. We begin with Fig. 8 showing the specific heat in function of temperature. Similarly as at $\phi=-63^{\circ}$, there is a broad maximum at temperatures comparable to the coupling constants in the Hamiltonian that, in the limit of small biases, does not depend on the applied bias. Therefore, it is not related to symmetry breaking in contrast to the sharp peaks at low temperatures that depend on the bias in a systematic way. We investigate the peaks in detail below finding results consistent with a continuous, symmetry breaking, phase transition.

In Fig. 9 we use $T^{*}$ obtained from the peaks of $m^{\prime}$ and $C_{V}$ to estimate the critical temperature and the exponents. Their best fits are collected in table 3. Averaging the results for $D=8$ we obtain:

[TABLE]

Using the duality (4) we obtain for $\phi=145.23^{\circ}$ in the zigzag phase:

[TABLE]

For a further self-consistency check, in the log-log plots in Fig. 10 we test the scaling ansatzes (13,14). We see that for $\xi(t^{*},h)$ deviations from the power law are significant. Furthermore, the range of $\xi(t^{*},h)$ is more limited here than for $\phi=-63^{\circ}$. For the maximal slope $m^{\prime}(t^{*},h)$ the deviations from the power law scaling are less significant than for $\xi(t^{*},h)$.

The results for $\phi=-17^{\circ}$ are more significantly affected by the deviations from the asymptotic scaling than the ones for $\phi=-63^{\circ}$. Nevertheless, they provide evidence that $T_{c}$ is small w. r. t. the couplings in the Hamiltonian.

V.3 Kitaev quantum spin liquid: $\phi=90^{\circ}$

The pure Kitaev model is tractable by quantum Monte-Carlo Nasu et al. (2015). At the same time its critical spin-liquid ground state makes the low temperature physics challenging for tensor networks as they may require large bond dimension $D$. Therefore, it is an ideal case to test limitations of our method by benchmarking it against the quantum Monte-Carlo results. The model was shown Nasu et al. (2015) to have two cross-over temperatures: $T_{1}<T_{2}$. Near $T_{2}$ it crosses over from a spin disordered paramagnet to a state with NN spin correlations. Near $T_{1}$ flux ordering takes place. This is where an expectation value of the plaquette flux operator, $W$, becomes non-zero before it converges to $1$ in the ground state at $T=0$.

In the absence of any phase transitions at finite temperature, there is no need to smooth the evolution by any bias field. Figure 11(a) shows specific heat in function of temperature. We can clearly see the peak at $T_{2}$ which is well converged in the bond dimension $D$. The second peak at the lower $T_{1}$ builds up with increasing $D$ in a systematic way. Location of the peaks is consistent with Ref. Nasu et al., 2015. The origin of the second peak is corroborated in Fig. 11(b) showing the average flux operator $W$ in function of temperature.

VI Conclusion

We applied the recently introduced tensor network algorithm to obtain thermal states of the Kitaev-Heisenberg model with a focus on their critical properties. As a technical advancement, we also show that the dynamical mapping from the hexagonal to the rhombic lattice makes the exact environment full update (eeFU) as efficient as the simple full update (FU) algorithm where the infinite tensor environment is delayed with respect to the imaginary time evolution. In the stripy phase at $\phi=-63^{\circ}$ we provide evidence for the second order phase transition and estimate its critical temperature at $T_{c}=0.056(4)$. Furthermore, for $\phi=-17^{\circ}$ in the antiferromagnetic phase, we estimate $T_{c}=0.076(15)$. Both critical temperatures are small w. r. t. the couplings in the Hamiltonian. By the duality transformation, these results can be mapped to, respectively, ferromagnetic and zigzag phases. Finally, we benchmark our method against quantum Monte-Carlo results in the special case of pure Kitaev model with a challenging spin-liquid ground state. We recover two crossovers for spin and flux ordering.

Acknowledgements.

We acknowledge insightful discussions with Andrzej M. Oleś on the KH model. We thank Joji Nasu for the QMC data in the Kitaev limit. This research was supported in part by the Polish Ministry of Science and Education under grant DI2015 021345 (AF) and the National Science Centre (Narodowe Centrum Nauki) under grant 2016/23/B/ST3/00830 (PC,AF,JD).