Fermionic algebraic quantum spin liquid in an octa-kagome frustrated antiferromagnet

TL;DR

This paper studies the ground state and finite-temperature properties of a spin-1/2 Heisenberg antiferromagnet on an octa-kagome lattice, revealing a fermionic algebraic quantum spin liquid phase at the isotropic point.

Contribution

It demonstrates the existence of a fermionic algebraic quantum spin liquid in an octa-kagome antiferromagnet using tensor network methods, identifying a quantum phase transition and characterizing different phases.

Findings

Identifies a quantum phase transition at J_d/J_t=0.6.

Finds a gapless spin liquid with power-law correlations at the isotropic point.

Shows linear specific heat and constant susceptibility at low temperature.

Abstract

We investigate the ground state and finite-temperature properties of the spin-1/2 Heisenberg antiferromagnet on an infinite octa-kagome lattice by utilizing state-of-the-art tensor network-based numerical methods. It is shown that the ground state has a vanishing local magnetization and possesses a $1/2$-magnetization plateau with up-down-up-up spin configuration. A quantum phase transition at the critical coupling ratio $J_{d}/J_{t}=0.6$ is found. When $0<J_{d}/J_{t}<0.6$, the system is in a valence bond state, where an obvious zero-magnetization plateau is observed, implying a gapful spin excitation; when $J_{d}/J_{t}>0.6$, the system exhibits a gapless excitation, in which the dimer-dimer correlation is found decaying in a power law, while the spin-spin and chiral-chiral correlation functions decay exponentially. At the isotropic point ($J_{d}/J_{t}=1$), we unveil that at low temperature ($T$) the specific heat depends linearly on $T$, and the susceptibility tends to a constant for $T\rightarrow 0$, giving rise to a Wilson ratio around unity, implying that the system under interest is a fermionic algebraic quantum spin liquid.

Equations

$$\text{H}=J_d\sum _{<ij>}{\text{S}}_i\cdot {\text{S}}_j+J_t\sum _{<lm>}{\text{S}}_l\cdot {\text{S}}_m-h\sum _i{S}_i^z,$$

$$|\Psi\rangle=\sum_{\{s\}}\sum_{\mu,\nu,\xi=1}^{D_{c}}\Bigg{(}\prod_{k\in\it TN}[\lambda_{1}]_{\mu\mu}[\lambda_{2}]_{\nu\nu}[\lambda_{3}]_{\xi\xi}A_{\mu\nu\xi}^{s_{k,1}s_{k,2}s_{k,3}}B_{\mu\nu\xi}^{s_{k,4}}|s_{k,1}s_{k,2}s_{k,3}s_{k,4}\rangle\Bigg{)},$$

$$|{\Psi }_{gs}⟩=\underset{\beta \to \infty }{\lim }\frac{e^{-\beta \text{H}}|\Psi ⟩}{\parallel e^{-\beta \text{H}}|\Psi ⟩\parallel },$$

$$\sum _{\mu ,{\mu }^{\prime }=1}^{D_c}{\delta }_{\mu {\mu }^{\prime }}[{\overline{\lambda }}_1{]}_{\mu \mu }[{\overline{\lambda }}_1{]}_{{\mu }^{\prime }{\mu }^{\prime }}{\overline{W}}_{\mu {\mu }^{\prime },\nu {\nu }^{\prime }}=\eta {\delta }_{\nu {\nu }^{\prime }},$$

$$\sum _{\nu ,{\nu }^{\prime }=1}^{D_c}{\overline{W}}_{\mu {\mu }^{\prime },\nu {\nu }^{\prime }}[{\overline{\lambda }}_1{]}_{\nu \nu }[{\overline{\lambda }}_1{]}_{{\nu }^{\prime }{\nu }^{\prime }}{\delta }_{\nu {\nu }^{\prime }}=\eta {\delta }_{\mu {\mu }^{\prime }}.$$

$$\sum _{g_{n\ne \epsilon }}{T}_{g_1{g}_2{g}_3{g}_4{g}_5{g}_6}^{cell}\prod _{n\ne \epsilon }{v}_{g_n}^n\propto v_{g_{\epsilon }}^{\epsilon }.$$

Full text

I Introduction

Quantum spin liquid (QSL), QSL also known as quantum disorder or quantum paramagnet, has received considerable attention since it was proposed to describe a possible magnetic disordered state in interacting spin systems even at temperature down to zero. It is intuitive that in two-dimensional (2D) quantum spin models highly geometric frustration and low coordination number usually lead to strong quantum fluctuations, capable of destroying the semi-classical long range orders in the ground state, and thereby inclining to generate a QSL. QMin2D In the past decades, there have been extensive numerical simulations Kagome1-nums ; Kagome2-nums ; Kagome3-nums ; Kagome4-nums ; Kagome5-nums and experimental efforts Kagome1-exp ; Kagome2-exp showing that the spin-$1/2$ Heisenberg antiferromagnetic model (HAFM) on kagomé lattice is the most promising QSL candidate. However, since the intractability of the quantum frustrated system, some unsettled issues are still remaining in hot debate, e.g., whether the ground state of kagomé HAFM is a gapped $\mathbb{Z}_{2}$ spin liquid or a gapless Dirac QSL.

Recently, a series of new layered compounds $BiOCu_{2}(XO_{3})(SO_{4})(OH)\cdot H_{2}O$, $(X=Te,Se)$ were discovered SKM ; OKL . The 2D framework built by magnetic $Cu^{2+}$ ions in these compounds shows an extremely unique lattice (see Fig. 1). Such a lattice, we dub it as octa-kagomé lattice (OKL), does not belong to any of the 2D uniform Archimedean lattices, which has not been considered before. OKL can be regarded as a variant of the standard kagomé lattice by inserting a dimer between the corner sharing triangles along one direction, which can also be viewed as corner and edge sharing octagons. Owing to strong geometric frustrations and lower coordination numbers in OKL, the spin-1/2 HAFM on OKL could be a long-sought QSL candidate more promising and intriguing than on kagomé lattice.

Motivated by the newly synthesized layered compounds $BiOCu_{2}(XO_{3})(SO_{4})(OH)\cdot H_{2}O$, $(X=Te,Se)$, we shall study systematically, for the first time, the ground state and thermodynamic properties of the spin-$1/2$ HAFM on OKL using state-of-the-art tensor network (TN)-based numerical methods. Our results show that the system under investigation possesses a ferminonic algebraic QSL phase. This paper is organized as follows: In Sec. II, the model and TN-based simulating methods are described in detail. In Sec. III, by calculating the local magnetization, we shall show that the ground state of this system is magnetically disordered. In Sec. IV, the spatial dependence of spin-spin, dimer-dimer and chiral-chiral correlation functions of the system under interest in the ground state will be given. In Sec. V, the magnetic curves and the phase diagram in the ground state are presented. In Sec. VI, the temperature dependence of specific heat and susceptibility will be discussed. Finally, we give a conclusion in Sec. VII.

II Model and Methods

The Hamiltonian under interest reads

[TABLE]

where $\textbf{S}_{i}$ is the spin operator on the $i$th site, $J_{d}$ ($J_{t}$) is the coupling constant between nearest neighbor spins standing inside the dimer (triangle), as indicated in Fig. 1, and $h$ is the magnetic field. We set $J_{t}=1$ as energy scale.

It is usually challenging to simulate quantum many-body systems. Due to strong correlations and quantum fluctuations, most traditional methods fail to capture their novel properties. For example, mean-field theories underestimate long range fluctuations that may be critically important to the exotic many-body phenomena; quantum Monte Carlo suffers from the notorious sign-problem QMCsign when calculating frustrated spin models as well as fermionic models away from the half filling; etc.

In this work, we use state-of-the-art TN algorithms to explore the spin-1/2 Heisenberg antiferromagnet on OKL. TN scheme is free from negative-sign problem, and has been demonstrated to be a powerful numerical tool not only in strongly correlated quantum systems, QMBSDMRGPEPS ; QMBS1 ; QMBS2 but also in statistical physics, statistic1-PEPS ; statistic2-TRG ; statistic3-HOSRG ; statistic4 quantum information QuantumInformation1 ; QuantumInformation2 ; QuantumInformation3 and so on. The central task in such kind of algorithms is to compute the TN contraction, QuantumInformation2 ; TNComp i.e. to sum over all shared bonds in TN. However, except some special cases, Tree ; MERA ; exactcontrac1 ; exactcontrac2 the contraction of the TN with a regular geometry (e.g. square or honeycomb) has been shown to be NP-hard. QuantumInformation2

Generally, there are two ways to deal with the TN simulations: renormalization statistic2-TRG ; iTEBD ; statistic1-PEPS ; PEPS1 ; PEPS2 ; iPEPS ; TERG ; CTMRG1 ; CTMRG2 ; SRG1 ; SRG2 ; statistic3-HOSRG and encoding ODTNS ; NCD ; AOP schemes. The former follows a contraction-and-truncation scheme, while the latter encodes the TN contraction into a local self-consistent problem. Specifically, the renormalization scheme originates from Wilson’s numerical renormalization group method, NRG ; NRG-Kondo which solves successfully the Kondo Kondo problem. Then, the density matrix renormalization group (DMRG) DMRG1 ; DMRG2 was proposed by White, where the boundary condition (especially in 1D) is better considered with entanglement. For 2D systems, algorithms based on tensor renormalization group and the infinite projected entangled pair state (iPEPS) iPEPS were proposed. The degrees of freedom is coarse-grained in such a way that when the tensor is invariant under renormalization, it represents approximately an infinite system.

The encoding scheme follows an opposite way known as the “mean-field” idea that considers well the entanglement with the help of TN. The “mean-field” idea is incredibly important in numerical physics, which gives birth of the great density functional theory DFTreview1 ; DFTreview2 and ab-initio scheme which has been widely used in both physics and chemistry. To better deal with the strong correlations in many-body physics, the dynamic mean-field theory DMFT1 ; DMFT2 ; DMFT3 ; DMFT4 ; DMFT5 and density matrix embedding theory DMET1 ; DMET2 ; DMET3 were also proposed. By combining “mean-field” idea with TN and multi-linear algebra, the ab-initio optimization principle was proposed, AOP where an infinite TN is equivalently transformed into a local tensor embedded in an entanglement bath.

We here employ three kinds of TN-based algorithms, namely cluster update cluster1 ; cluster2 ; cluster3 and full update schemes iPEPS ; full-update2 ; full-update3 of the iPEPS PEPS1 ; PEPS2 ; iPEPS (a contraction-and-truncation scheme) and network contractor dynamics NCD (NCD) approach (an encoding scheme) to investigate our model for mutually validating the results obtained by each scheme. Consequently, the calculated results are consistent with each other, which manifests itself the reliability of our simulations.

II.1 Tensor-network state ansatz

We start from a TN state ansatz, as shown in Fig. 2, to describe the states at zero temperature. Such a TN state is composed of two non-equivalent tensors $A$ and $B$, and three different diagonal matrices $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$. $A$ and $B$ (each of which contains three virtual bonds that carry the entanglement of the state) are located on the two inequivalent triangles of the OKL, respectively. The physical degrees of freedom of the three spins in triangle A are put on tensor $A$, where the dimension of the physical space is 8. In this way, the dimension of the physical bond of tensor $B$ is 2, which is the Hilbert space of the spin on the right side with the $J_{d}$ coupling. Mathematically, such a TN state is written as

[TABLE]

where $k$ refers to the $k$-th unit cell of the whole lattice with $\textbf{a}_{1}$ and $\textbf{a}_{2}$ basis vectors [see Fig. 2 (b)]. To get the ground state, the imaginary time evolution is implemented to minimize the energy of the PEPS by

[TABLE]

where $\beta=1/{k_{B}T}$.

It is impossible to calculate Eq. (3) exactly in the thermodynamic limit, since the dimension of H increases exponentially with the number of lattice sites. Here, we use the Trotter-Suzuki decomposition to implement the evolution on the TN state. By splitting Hamiltonian into two parts, one has $H_{a}=\sum_{k}H_{left-trangle}^{[k]}$ and H_{b}=\sum_{k}\bigg{(}H_{dimer}^{[k]}+H_{right-trangle}^{[k]}\bigg{)}, and the first-order Trotter-Suzuki decomposition can be used to approximate the evolution operator, i.e., e^{-\beta H}\approx\big{(}e^{\tau H_{a}}e^{\tau H_{b}}\big{)}^{N}+\mathcal{O}(\tau^{2}), with $\beta=N\tau$. The approximation becomes accurate when the Trotter step $\tau$ approaches zero. In practical calculations, we decrease $\tau$ gradually from $1\times 10^{-1}$ to $1\times 10^{-5}$ so that the Trotter error becomes negligible.

By considering the translation invariance, we can adopt the local operation instead of evolving the whole system, and optimize the environment around the local tensors. Incidentally, for finite-temperature thermal states, the imaginary-time evolution of the density operator can be implemented similarly.

II.2 Cluster update

We choose a hexagon consisting of six tensors as the environment for cluster update, as depicted in Fig. 3 (a). The cluster tensors are transformed into a super-orthogonal form ODTNS in order to approximate the global environment optimally. Taking Fig. 3 (b) as an example, we build a double-layer structure of the cluster tensor and contract all physical indices and virtual bonds on the bra and ket layers except the bonds connected by $\lambda_{1}$. For convenience, the shaded part of Fig. 3 (b) is represented by $W$. The super-orthogonalization is much like the canonicalization for an infinite 1D lattice iTEBD-canon . The update of $\overline{\lambda}_{1}$ and $\overline{W}$ leads to the conditions

[TABLE]

[TABLE]

Fig. 3 (c) is the graphical representation of Eqs. (4) and (5). The update of $W$ is actually acting on $A$ and $B$ along the $\lambda_{1}$ direction, where $A$ and $B$ are renewed to $\overline{A}$ and $\overline{B}$. Operations on the other two directions are similar. We iterate this procedure until the cluster satisfies simultaneously the orthogonality conditions in all three directions. Then, the environment of the cluster can be best approximated by the converged diagonal matrices $\overline{\lambda}_{1}$, $\overline{\lambda}_{2}$ and $\overline{\lambda}_{3}$.

Then we permute the physical indices from A to B to evolve the interactions on the B triangles. This operation will increase the bond dimensions, and a truncation is needed. Taking Fig. 3 (d) as an example, we leave one physical index and the corresponding virtual bonds of A open and others contracted in the cluster. We use $M$ to denote the intermediate reduced density matrix, where the dimension of $M$ is $2D_{c}\times 2D_{c}$. Moreover, $M$ is a Hermitian matrix because of the double-layer structure. Then, we decompose $M$ using the SVD and only keep the basis corresponding to the $D_{c}$ dominant singular values. This procedure is shown in Fig. 4 (a), where $U_{1}$ is the unitary matrix given by the SVD holding the spared physical index of $M$, and $\overset{=}{\lambda}_{1}$ is the square root of the singular spectrum after truncation. $U_{2}$, $U_{3}$, $\overset{=}{\lambda}_{2}$ and $\overset{=}{\lambda}_{3}$ are obtained in the similar way.

Finally, we change the position of all three physical indices from $\overline{A}$ into $\overline{B}$, as depicted in Figs. 4 (b) and (c). In such a way, the evolutions given by the interactions of the triangles A and B are implemented in turn, where the geometry and the bond dimensions are kept unchanged.

II.3 Full update

Unlike the cluster update scheme, the full update scheme needs to contract all the 2D TN in order to truncate and obtain physical quantities. There are two widely used ways to simulate the whole environment, namely the iTEBD iTEBD and the corner transfer matrix renormalization group (CTMRG) CTMRG1 ; CTMRG2 . Here, we use the iTEBD in our calculation, where the TN is contracted to a matrix product state (MPS) on its boundary. Full update can achieve a higher accuracy than local optimization methods, but the computation cost is significantly large. We set the bond dimension of the MPS in iTEBD as $D_{c}^{2}$ to balance the accuracy and cost. The permutation of physical indices and variational optimization of the truncation matrices follow Refs. iPEPS, ; cluster3, ; full-update2, ; full-update3, .

II.4 Network Contractor Dynamics

NCD was first proposed to solve the TN contractions in the calculation of partition function of 2D quantum models. We adopt the NCD algorithm to optimize the environment around the local tensors in ground state simulation. Different from the renormalization, NCD follows a TN encoding strategy AOP , and the TN structure is also different from that in Fig. 2. The specific cell tensor is shown in Fig. 5 (a), which contains $6$ nonequivalent tensors with $A_{1}$, $B_{1}$, $A_{2}$ and $B_{2}$ located on triangles, and $P_{1}$, $P_{2}$ on dimers. When mapping onto the OKL, there are eight inequivalent lattice sites in the cell tensor of NCD, twice as large as the cell tensors of cluster update and full update schemes. Physical indices are on $P_{1}$, $A_{2}$ and $B_{2}$, indicating that the Hamiltonian splits into two parts, where the first part contains interactions between spins sitting on the dimer denoted by $P_{1}$ and triangles denoted by $A_{2}$ and $B_{2}$, and the rest of interactions are included in the second part of the Hamiltonian. Imaginary time evolution is applied to minimize the ground state energy, and consequently, NCD procedure plays the role of super-orthogonalization to approximate the whole TN contraction. As explained in Ref. NCD, , the whole TN contraction is simplified to a local contraction of a tensor cluster $T^{cell}$ with six contractors ${v^{n}},n\in[1,2,3,4,5,6]$. $T^{cell}$ is a six-order tensor obtained by contracting the physical indices and connected virtual bonds of the double-layer cell tensors, as represented in Fig. 5 (b), where we use bold black lines to indicate a fat index that contains double virtual bonds in one of the six directions. While each contractor $v^{n}$ is a vector with the same dimension as the $n$th index of $T^{cell}$. $T^{cell}$ and ${v^{n}}$ need to meet the following self-consistent relation

[TABLE]

Eq. (6) includes six self-consistent equations for $\varepsilon\in[1,2,3,4,5,6]$, which should be satisfied simultaneously. Analytically, the six contractors solved from Eq. (6) are precisely those given by rank-1 decomposition of $T^{cell}$. The rank-1 tensor, which we call a “defect”, is given by a direct product of the six contractors. The graphic representation of “defect” is shown in Fig. 5 (c). The “defect” is actually the first-order approximation of $T^{cell}$. If one substitutes the minimal number of $T^{cell}$’s with the “defects” so that no loop appears, then the original TN will become a tree framework. Thanks to the self-consistent conditions, there is no need to compute the whole contraction of such a tree, and only are the local contraction of $T^{cell}$ and the contractors ${v_{n}}$ required. In this sense, the physical quantities calculated from the “defective” TN can be viewed as a mean-field approximation of the exact one. In addition, we can introduce more loops into the defected TN to achieve a higher accuracy, but the computing cost increases inevitably.

III Disordered ground state

Let us now study the ground state properties of the spin-1/2 HAFM on the infinite OKL for isotropic point ($J_{d}=1$). To testify the reliability of our calculations, we compare the ground state energy $e_{g}$ (per site) obtained by different methods including NCD, cluster and full updates of iPEPS (Fig. 6). We found that for large bond dimension $D_{c}>4$, all schemes give consistent results, showing the reliability of our calculations. A power law dependence is found, with which the energy of infinite $D_{c}$ is $e_{g}=-0.4524$ by extrapolation, which is lower than $-0.4386(5)$, the extrapolated ground state energy of the spin-1/2 HAFM on kagomé lattice given by DMRG. Kagome3-nums ; Kagome4-nums

In Fig. 7, we present the local magnetization $\langle S_{i}^{\alpha}\rangle$ on each nonequivalent site in a unit cell for $h=0$. Small values of $\langle S_{i}^{z}\rangle$ and $\langle S_{i}^{x}\rangle$ caused by the truncation error can be observed. By increasing $D_{c}$, $\langle S_{i}^{z}\rangle$ and $\langle S_{i}^{x}\rangle$ decay rapidly, and the data are fitted (the dashed lines in Fig. 7) with the function $f(D_{c})=p(1/D_{c})^{q}$, where $p$ and $q$ are fitting parameters. Thus, the extrapolation of local magnetic moments gives a zero magnetization in $D_{c}\rightarrow\infty$ limit. In particular, as the average values of $\langle S_{i}^{z}\rangle$ and $\langle S_{i}^{x}\rangle$ only fluctuate in the vicinity of zero with the increase of $D_{c}$, we may use a linear fitting for the average magnetic moments that gives negligible intercepts about $10^{-5}$. The absence of local magnetic moments strongly suggest that it does not have conventional magnetic orders in the ground state, i.e., no traditional SO(3) symmetry is broken.

IV Spin-spin, dimer-dimer and chiral-chiral correlation functions

In Fig. 8, we present the spatial dependence of several correlation functions in the ground state for the system under interest with $J_{d}=1$. The spin-spin correlation function $|\langle S_{i}^{z}S_{j}^{z}\rangle|$ along the horizontal axis is found to decay exponentially, satisfying $f(|i-j|)=\alpha\exp(-|i-j|/\xi)$ with $\alpha=0.084$ and the correlation length $\xi=0.16$, which shows that the spin-spin correlation of this system is short-ranged and the ground state is magnetically disordered.

The chiral-chiral correlation function is defined by $|\langle C_{i}C_{j}\rangle|=|\langle[\textbf{S}_{i_{1}}\cdot(\textbf{S}_{i_{2}}\times\textbf{S}_{i_{3}})][\textbf{S}_{j_{1}}\cdot(\textbf{S}_{j_{2}}\times\textbf{S}_{j_{3}})]\rangle-\langle\textbf{S}_{i_{1}}\cdot(\textbf{S}_{i_{2}}\times\textbf{S}_{i_{3}})\rangle\langle\textbf{S}_{j_{1}}\cdot(\textbf{S}_{j_{2}}\times\textbf{S}_{j_{3}})\rangle|$, where the lattice sites $i$ and $j$ belong to the left-triangles along the horizontal direction. It is found that the chiral-chiral correlation function also decays exponentially with $\alpha=0.0054$ and $\xi=0.11$, revealing the absence of a long-range spin chiral order.

The dimer-dimer correlation function, which is defined by $|\langle D_{i}D_{j}\rangle|=|\langle(S_{i}^{z}S_{i+1}^{z})(S_{j}^{z}S_{j+1}^{z})\rangle-\langle S_{i}^{z}S_{i+1}^{z}\rangle\langle S_{j}^{z}S_{j+1}^{z}\rangle|$ for the $i$-th and $j$-th dimers, is disclosed to exhibit a power-law decay as of the form $|\langle D_{i}D_{j}\rangle|\sim 1/|i-j|^{\eta}$ with $\eta=17.96$ [Fig. 8(b)]. This fact signatures possible existence of an algebraic QSL in this system.

Here it is interesting to ask if the correlations along the vertical axis behave the same as those along the horizontal axis. To answer this question, we also calculated the three correlation functions in the vertical direction. The results show that the behaviors in this direction are different, as shown in Fig.  8. It is seen that all three correlations along the vertical axis decay exponentially, fitted by the function $f(|i-j|)=\alpha\exp(-|i-j|/\xi)$, with $|\langle S_{i}^{z}S_{j}^{z}\rangle|$ fitted with $\alpha=0.034$ and $\xi=0.94$, $|\langle C_{i}C_{j}\rangle|$ fitted with $\alpha=0.00075$ and $\xi=0.47$, and $|\langle D_{i}D_{j}\rangle|$ fitted with $\alpha=0.00038$ and $\xi=0.68$. This is owing to the nonequivalent lattice structure along the two axes. It is the introduction of $J_{d}$ that causes the lattice essentially distinct from a combination of decoupled zig-zag spin chains, of which the ground state is a valence bond state (VBS) with two-fold degeneracy and a finite magnetic excitation gap. sawtooth1 ; sawtooth2 A strong $J_{d}$ coupling (especially at the isotropic point) is crucial in the critical phase, which will be discussed later.

We would like to mention that the nature of the correlations presented here is similar to the case with a resonating valence bond (RVB) wave function constructed on a square lattice, where an exponentially decaying spin-spin correlation and a power-law decaying dimer-dimer correlation were also observed. NNRVB square ; J1-J2 square

V Magnetic curves and phase diagram in ground state

The magnetization per site $M$ as a function of magnetic field $h$ in the spin-1/2 HAFM on OKL is presented in Fig. 9. One may observe that in magnetic curves [Fig. 9 (a)], for $J_{d}=0.2$, three plateaux with $M=0,\frac{1}{4}$ and $\frac{1}{2}$ are observed, while for $J_{d}=0.8$ and $1$, apart from the two plateaux with $M=\frac{1}{4}$ and $\frac{1}{2}$, no $M=0$ plateau is found. These results imply that depending on $J_{d}$, there may be two phases in the system, one phase with a zero-magnetization plateau and the other phase without. As the width of the $M=0$ plateau gives the gap from the singlet ground state to the first triplet excited state, we find that in the phase with small $J_{d}$ the spin excitation is gapful, while in the other phase with large $J_{d}$ it is gapless. For a closer inspection, we calculated the cases with small $J_{d}$ under weaker magnetic fields, as given in Fig. 9 (b). The results demonstrate that the spin gap decreases with increasing $J_{d}$, suggesting that there must be a critical point $J_{d}^{c}$, at which a quantum phase transition (QPT) happens: for $J_{d}<J_{d}^{c}$ the ground state is in a gapped phase, and for $J_{d}>J_{d}^{c}$ it is in a gapless phase.

Another interesting phenomenon in magnetic curves is the occurrence of $M=1/4$-plateau, which can also be called $M/M_{s}=1/2$-plateau (briefly 1/2-magnetization plateau) with $M_{s}=1/2$ the saturation magnetization per site. The frustrated Heisenberg models on lattices with triangular structures lead usually to 1/3-magnetization plateau, which has been found in, e.g., kagomé kagomeplateau1 ; kagomeplateau2 ; kagomeplateau3 and Husimi Husimi lattices. The occurrence of the 1/2-magnetization plateau in the present system is understandable, because the unit cell of the OKL contains four inequivalent lattice sites, leading to the periodicity $n$ of the ground state is 4, consistent with the condition of $n(S-M)=integer$. To explore the nature of this 1/2-plateau, we calculated the local magnetic moment at four inequivalent sites in a unit cell, and found in this plateau phase the spin configuration is of up-down-up-up (UDUU), as illustrated in the inset of Fig. 9 (a). Such a plateau is a commensurate, classical state stabilized by quantum fluctuations.

To determine accurately the quantum critical point (QCP) $J_{d}^{c}$, we calculated the spin gap as a function of $J_{d}$, as given in Fig. 10, which gives $J_{d}^{c}=0.6$. To further confirm this point, we also studied the second-order derivative of the ground state energy with respect to $J_{d}$ (the inset of Fig. 10), which reveals a sharp dip at the same point, indicating the QPT indeed appears at $J_{d}^{c}=0.6$.

By summarizing our calculated results, we present the ground state phase diagram of the spin-1/2 HAFM on OKL in the $J_{d}-h$ plane, as shown in Fig. 11. It can be seen that when $h=0$, the phase for $J_{d}<0.6$ is a VBS, as in the limit of $J_{d}\rightarrow 0$, the system approaches to an uncoupled zig-zag spin chain, whose ground state is a VBS with twofold degeneracy and a finite magnetic excitation gap. sawtooth1 ; sawtooth2 Because there is no quantum phase transition for $J_{d}<J_{d}^{c}$, the system should stay in the same VBS phase. For $J_{d}>0.6$, the system enters into a gapless QSL state, which is evidenced by the algebraically decaying dimer-dimer correlations and vanishing local magnetic moments. When $h>0$, the VBS state is gradually melted by closing the spin gap, and the system enters into a spin canted state. By increasing the magnetic field further, the system enters into the 1/2-magnetization plateau (with $M=1/4$) phase, in which the spin gap opens again, and the spin configurations are arranged in UDUU alignments. Above the UDUU phase, it enters into another spin canted phase. By increasing $h$ further, all spins are polarized. It should be remarked that all the phase boundaries in this phase diagram are obtained by observing various critical magnetic fields.

VI Thermodynamic properties

Next, we explore the thermodynamic properties of the spin-$1/2$ HAFM on OKL using the optimized decimation of tensor network state. ODTNS The free energy can be obtained by collecting all renormalization factors down to the targeted temperature. Alternatively, one can also get the physical quantities by calculating the expectation values of local operators with tensor-network thermodynamic states. Considering the precision and cost of thermal-state TN algorithms, we choose the cluster update scheme to contract the “environment” around the local inequivalent tensors. The energy as well as other thermodynamic quantities including specific heat and susceptibility are thus calculated. To keep a higher accuracy, we adopted the second-order Trotter-Suzuki decomposition and fix the Trotter slice to be 0.01 in the calculations of thermodynamic properties.

We obtain the temperature dependence of the specific heat by $C(T)={\partial f}/{\partial T}$, where $f$ is the free energy per site. Fig. 12 gives the results for $J_{d}=10^{-4}$ and $1$. It is observed that at high temperature, both go to converge, and $C(T)$ decreases down to zero with increasing temperature. But at low temperature (see the inset of Fig. 12), both cases show intrinsically distinct behaviors: the specific heat for $J_{d}=10^{-4}$ exhibits two peaks and is pretty close to the exact diagonalization (ED) result of the zig-zag spin chain with 8 triangles, which also verifies the reliability of our method. When $T\rightarrow 0$, $C(T)$ shows an exponentially decaying behavior, suggesting that there should be a finite excitation gap, being well consistent with the result in the ground state, as the system in this case is almost dimerized; for $J_{d}=1$, the specific heat exhibits a single peak, and when $T$ is very low, $C(T)$ obeys a polynomial behavior of the form $C(T)=1.741T-7.96T^{2}+22.51T^{3}-25.61T^{4}$. When $T\rightarrow 0$, $C(T)$ is linearly dependent on temperature, which indicates the existence of gapless excitations, and implies that the system is critical. It is also consistent with the preceding result that the ground state is an algebraic QSL.

Such criticality is further evidenced by the susceptibility at low temperature.The susceptibility is calculated according to $\chi(T)=[M(h+\Delta h)|_{T}-M(h)|_{T}]/\Delta h$, where $\Delta h=0.01$ is taken. The results for $J_{d}=10^{-4}$ and $1$ are presented in Fig. 13. One may see that both curves obey the Curie-Weiss law at high temperature and exhibit a sharp peak at low temperature due to antiferromagnetic interactions. Significant differences occur when $T\rightarrow 0$. As shown in the inset of Fig. 13, for $J_{d}=10^{-4}$, $\chi$ goes to zero in an exponential way, revealing the existence of a finite spin gap, while for $J_{d}=1$, $\chi$ converges to a finite constant in a polynomial of the form $\chi(T)=0.09566-0.06828T+2.438T^{2}$, being reminiscent of a Luttinger liquid behavior, and consistent again with the critical feature of the ground state.

In addition, it is quite interesting to look at the Wilson ratio (WR) $R_{w}$ for the present critical system at the isotropic point. The WR is defined by $R_{w}=(4/3)(\pi k_{B}/g\mu_{B})^{2}\chi/(C/T)$, where $\chi$ is the susceptibility, $C$ is the specific heat, $k_{B}$ is the Boltzmann constant, $g$ is the Landé factor, and $\mu_{B}$ is the Bohr magneton. For simplicity we have assumed $k_{B}=g\mu_{B}=1$. It is known that for free electron gas, $R_{w}=1$. For most QSL theories, the WR is usually less than one. QSL For the present system with $J_{d}=1$, at $T\rightarrow 0$, $\chi$ tends to a constant, and $C(T)\sim T$, which gives $R_{w}\approx 0.72$. In consideration of the fact that the linear temperature dependence of the specific heat resembles the Luttinger liquid behavior, and the WR $R_{w}$ is on the order of unity, which are analogous to the behaviors induced by fermionic quasiparticles, we conclude that the present isotropic system should be a fermionic gapless QSL.

VII Conclusion

The ground state and thermodynamic prosperities of the spin-$1/2$ HAFM on OKL have been systematically studied with the aid of powerful TN numerical simulations. We adopted three kinds of TN algorithms in calculations of the ground state energy per site, which gives $-0.4524$ by the infinite $D_{c}$ extrapolation in the thermodynamic limit, lower than $-0.4386(5)$ on kagomé lattice. The magnetic order is melted in the ground state due to strong frustration induced by corner sharing triangles. A QPT is found in this system. It is disclosed that below the QCP, the system has a finite spin gap and is in a VBS state, while above the QCP, the system is in a gapless QSL state. At the isotropic point, we uncover that the dimer-dimer correlation function decays algebraically, while the spin-spin and chiral-chiral correlation functions behavior in an exponential way. In addition, the specific heat at low temperature is shown to depend linearly on temperature, exhibiting a Luttinger liquid behavior, and the susceptibility tends to a finite constant when $T\rightarrow 0$, which indicates a gapless excitation in the system. The Wilson ratio is found to be 0.72, close to 1. All these features reveal that the isotropic spin-1/2 HAFM on OKL is a fermionic gapless QSL.

Acknowledgements.

The authors acknowledge Wei Li, Xin Yan, and Yi-Zhen Huang for useful discussions, and also appreciate Guang-Zhao Qin and Xuan-Ting Ji for kind help. In particular, we thank Ying-Ying Tang and Zhang-Zhen He for useful discussions on the compounds with OKL. This work was supported in part by the MOST of China (Grants No. 2012CB932900 and No. 2013CB933401), the NSFC (Grant No. 14474279), and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB07010100). SJR was supported by ERC ADG OSYRIS, Spanish MINECO (Severo Ochoa grant SEV-2015-0522, FOQUS grant FIS2013-46768), Catalan AGAUR SGR 874, Fundació Cellex, and EU FETPRO QUIC.