Long-range spin-orbital order in the spin-orbital SU(2)$\times$SU(2)$\times$U(1) model

Introduction. Symmetry may be one of the fundamental concepts involved in physics. Of all these symmetries SU(2) is ubiquitous in condensed-matter physics because the spin operators are its generators. On the other hand, the success of larger symmetry groups such as SU(3) in particle physics raised a natural question whether high-symmetry groups are relevant for condensed-matter physics. It was proposed that high symmetries can be realized in spin-orbit-coupled compounds [1, 2]. This has inspired a variety of theoretical studies on spin-orbit-coupled insulators [5, 3, 4, 7, 6, 10, 9, 8, 11]. At the meanwhile, the experimental development of measuring the orbital degrees of freedom in transition-metal oxides [12, 13] has greatly boosted the investigation on orbital physics in materials with strong spin-orbit coupling and crystal-field splitting. A minimal high-symmetry model involving both the spin and orbital degrees of freedom may be the SU(4) symmetric Kugel-Khomskii model [1, 2, 10, 14, 16, 15, 17], where orbital degrees of freedom are represented as pseudospin and coupled with the spin degree of freedom on each bond. It is argued that such a model is relevant to the observed spin liquid states [18, 19, 20] in $\rm{LiNiO_{2}}$ [2], $\rm{Ba_{3}CuSb_{2}O_{9}}$ [10], and $\alpha$-$\rm{ZrCl_{3}}$ [14].

On the other hand, Yamada, Oshikawa, and Jackeli (YOJ) recently proposed [16] another spin-orbital model to describe some $d^{1}$ materials [21, 22, 23] on the triangular lattice. In addition to the geometry frustration, there is a spin-orbital frustration in this model, and it is argued [16] that such a model may host a spin-orbital liquid ground state. However, there are no systematic studies so far, which motivates us to investigate this model.

In the YOJ model, there are two kinds of terms in the Hamiltonian, as plotted in Fig. 1. The lattice sites are represented by filled circles. The nearest-neighbor sites are connected by blue solid or dashed lines, which represent the two interaction terms, respectively. The solid lines represent the term

$$H_{ij}=\left(\bm{S}_{i}\cdot\bm{S}_{j}+\frac{1}{4}\right)\left(\bm{T}_{i}\cdot\bm{T}_{j}+\frac{1}{4}\right),$$

(1)

where $\bm{S_{i}}$ and $\bm{T_{i}}$ are both pseudospin-1/2 operators defined on lattice site $i$, corresponding to effective spin and orbital degrees of freedom, respectively. The dashed lines are the term given by

$$H^{\prime}_{ij}=\left(\bm{S}_{i}\cdot\bm{S}_{j}+\frac{1}{4}\right)\left(T^{z}_{i}T^{z}_{j}-T^{x}_{i}T^{x}_{j}-T^{y}_{i}T^{y}_{j}+\frac{1}{4}\right).\\$$

(2)

The Hamiltonian $\mathcal{H}$ of the YOJ model is then the summation of all these terms.

Symmetry and method. Although Eq. (2) breaks the orbital SU(2) symmetry to U(1), it turns out that this model has a SU(2)$\times$SU(2)$\times$U(1) symmetry. To see this, one can construct [25] the generators $\bm{X}_{i}$ and $\bm{Y}_{i}$ as follows,

$$X^{\alpha}_{i}\equiv S^{\alpha}_{i}\left(\frac{1}{2}+T^{z}_{i}\right),\quad Y^{\alpha}_{i}\equiv S^{\alpha}_{i}\left(\frac{1}{2}-T^{z}_{i}\right)$$

(3)

where $\alpha=x,y,z$. It can be easily verified that

	$\displaystyle[X^{\alpha}_{i},X^{\beta}_{j}]$	$\displaystyle=i\varepsilon_{\alpha\beta\gamma}X^{\gamma}_{i}\delta_{ij},$
	$\displaystyle[Y^{\alpha}_{i},Y^{\beta}_{j}]$	$\displaystyle=i\varepsilon_{\alpha\beta\gamma}Y^{\gamma}_{i}\delta_{ij},$
	$\displaystyle[X^{\alpha}_{i},Y^{\beta}_{j}]$	$\displaystyle=0,$
	$\displaystyle[T^{z}_{i},X^{\alpha}_{j}]$	$\displaystyle=0,$
	$\displaystyle[T^{z}_{i},Y^{\beta}_{j}]$	$\displaystyle=0$		(4)

are satisfied. Moreover, the following seven operators are conserved quantities with respect to the Hamiltonian,

$$X^{\alpha}\equiv\sum_{i}X^{\alpha}_{i},\quad Y^{\alpha}\equiv\sum_{i}Y^{\alpha}_{i},\quad T^{z}\equiv\sum_{i}T^{z}_{i}.$$

(5)

Knowing the commutation relation of these operators, we have found seven generators of a Lie group which can be written as SU(2)$\times$SU(2)$\times$U(1), and we can use these generators to define the corresponding order parameters as usually done for spin operators [24].

To study the ground-state properties of this model, we employ the tensor-network method [26, 27, 28, 29], which is a class of numerical methods based on the tensor-network state representation of the targeted quantum state and the network contraction techniques arising from the idea of renormalization group. It is free of the sign problem, can study the thermodynamic limit directly under the help of translational invariance, and has been successfully applied to study strongly-correlated electron systems [31, 30], frustrated spin systems [33, 34, 32, 24], statistical models [35, 36, 37], topological order [38, 39, 40, 41, 42], quantum field theory [43, 44, 45], machine learning [46, 47], and even quantum circuit simulation [48], etc. In this work, we use the infinite projected entangled simplex state (PESS) ansatz [49, 50] to represent the trial wavefunction, and employ the corner transfer-matrix renormalization group (CTMRG) method [51, 52, 31] combined with the nested tensor network technique [53] to estimate the physical quantities, such as the ground-state energy and order parameters, etc. Unexpectedly, we find that in the ground state two SU(2) symmetries are broken. As we have shown, the generators of two SU(2) groups include both the spin and orbital operators, which suggests that the corresponding orders are different from the conventional magnetic order [24]. Hereafter, following the literature [14, 16] we call it a spin-orbital order.

More specifically, the PESS employed in this work is a generalization of the projected entangled pair state ansatz [27], and has been successfully applied to the highly frustrated antiferromagnetic Heisenberg model on Kagome lattice [49, 32, 54] and triangular lattice [24]. Similar to that in Ref. [24], in this work, the PESS wavefunction is represented as follows,

	$\displaystyle|\Psi\rangle=\sum_{\{\sigma,\sigma^{\prime}\}}\mathrm{Tr}$	$\displaystyle\left(...S^{(\mu\nu)}_{i_{\mu\nu}j_{\mu\nu}k_{\mu\nu}}A^{(\lambda\omega)}_{i_{\lambda\omega}j_{\lambda\omega}k_{\lambda\omega}}[\sigma_{\lambda\omega}\sigma^{\prime}_{\lambda\omega}]...\right)$
		$\displaystyle|...\sigma_{\lambda\omega}...\rangle|...\sigma^{\prime}_{\lambda\omega}...\rangle$		(6)

as illustrated in Fig. 1. Here $(\mu,\nu)$ denotes the coordinates of the upward triangles, at the center of which a rank-3 simplex tensor $S$ is introduced to characterize the entanglement in that triangle. $(\lambda,\omega)$ denote the coordinates of lattice sites, where a rank-5 projection tensor $A$ is defined, with three virtual indices labeled as $i,j,k$ and two physical indices labeled as $\sigma$ (spin) and $\sigma^{\prime}$ (orbital). Every two virtual indices associated with the same bond take the same values. “$\mathrm{Tr}$” is over all the repeated virtual indices and $\sum$ is over all the physical indices.

The bond dimension $D$, which is the highest value that the virtual indices can take, is an important parameter in tensor network states. Generally, the larger $D$ is, the more accurate the obtained representation is, but the heavier the computational cost is at the meanwhile. Therefore, in practical calculations, one has to make a good balance between accuracy and cost. In this work, using the efficient nested tensor network technique [53], we have pushed $D$ to 18 and obtained reliable features of the ground state.

The ground-state wavefunction is obtained by imaginary-time evolution. In order to make the calculation more efficient, the wavefunction is updated by the simple update algorithm [55, 56]. Though for a given $D$, simple update might be less accurate than the full update [57] or direct variational calculation [58], it can produce wavefunction with a much larger $D$ and the numerical accuracy can thus be remedied properly, as exemplified in Ref. [32]. Moreover, to avoid bias and reduce the Trotter error, we start from a wavefunction randomly generated on complex field, and gradually reduce the Trotter step $\tau$ from $0.5$ to $10^{-3}$, which turns out to be sufficiently small to estimate the physical quantities.

After mapping the obtained infinite PESS wavefunction to a square network [59], we measure the physical observables through the CTMRG method developed for an arbitrary unit cell on the square lattice [31]. The central idea of this method is to represent the surrounding environment of a local tensor in terms of corner matrices and edge tensors approximately. A key parameter therein is the bond dimension $\chi$ of the environment tensors, which controls both the accuracy and cost during the measurement. In principle, our results should also depend on $\chi$. We will show later such dependence is very weak and thus it can be neglected when the $\chi$ is large enough, i.e., $\chi\geq{D^{2}}$. In our calculation, the maximal $\chi$ is no less than $D^{2}$ to ensure a reliable result.

Ground-state energy. In the tensor-network simulations, the cluster size should be compatible with the unit cell if the ground state has a long-range order. However, since the ground state is unknown a priori, we need to try different cluster sizes to determine the unit cell. A correct cluster size should have the lowest energy in correspondence to the ground state. For this purpose, we compare the energy obtained from the PESS ansatz with various clusters. We have checked four different clusters compatible with the Hamiltonian symmetry, i.e., $2\times 2$, $3\times 3$, $4\times 2$, $4\times 4$, which are all illustrated in Fig. 1. The energy $E$ obtained from the ansatz with these four clusters is shown in Fig. 2 for comparison. It shows that the energy obtained on the $2\times 2$ and $3\times 3$ clusters is obviously higher than that on the $4\times 2$ and $4\times 4$ clusters. Moreover, the energy is the same on the latter two clusters, suggesting that the $4\times 2$ is the minimal cluster compatible with the ground state. Therefore, in all the rest calculations, we focus on the ansatz on the $4\times 2$ cluster.

In Fig. 2, we extrapolate the energy per bond $E$ as a function of $1/D$ to the large-$D$ limit by the formula

$\displaystyle E(1/D)=E_{g}+a_{0}(1/D)+a_{1}(1/D)^{2},$

(7)

where $E_{g}$, $a_{0}$ and $a_{1}$ are the fitting parameters. We finally end with the ground-state energy per bond $E_{g}=-0.0793(5)$. It may serve as a benchmark for future works.

Spin-orbital order. It is argued that the quantum fluctuation is strong [60, 61, 14] in high-symmetry models and thus spin-orbital liquid is favored therein. This argument is supported by several works [10, 15, 17] on the SU(4) spin-orbital model on various lattices. In this context, the main concern on this model is whether its symmetries are broken or not, in particular, whether the SU(2)$\times$SU(2)$\times$U(1) symmetry is broken or not. To check this, we calculate the corresponding local order parameters, which are defined as the expectation values of the generators in the ground state. The magnitude of the spin-orbital order at site $i$ is defined by $m_{i}^{X}=\sqrt{\sum_{\alpha}{\langle X_{i}^{\alpha}\rangle}^{2}}$ and $m_{i}^{Y}=\sqrt{\sum_{\alpha}{\langle Y_{i}^{\alpha}\rangle}^{2}}$ corresponding to two SU(2) groups. Their average in one unit cell is then defined by $m^{X/Y}=\sum_{i}{m_{i}^{X/Y}}/N$ with $N=4\times{2}=8$.

In Fig. 3, we show $m^{X/Y}$ as a function of $1/D$ with $D$ ranging from 4 to 18. One may notice that $m^{X}$ and $m^{Y}$ coincide very well. This can be understood as follows. Let us check the two SU(2) groups. Although their generators $X_{i}^{\alpha}$ and $Y_{i}^{\beta}$ commute, we can see that $X_{i}^{\alpha}\leftrightarrow{Y_{i}^{\alpha}}$ under the transformation $T_{i}^{z}\rightarrow{-T_{i}^{z}}$. In addition, $m^{X/Y}$ decreases monotonically as $D$ increases. We try to fit such a curve by various polynomial functions of $1/D$, which gives $m^{X/Y}=0.085(10)$ in the large-$D$ limit. Obviously, the finite magnitudes tell us that two SU(2) symmetries are spontaneously broken. Moreover, the magnitudes are about half of that of the spin-1/2 anfiferromagnetic Heisenberg model (see, for example, Table I in Ref. [24]). This is consistent with the argument that quantum fluctuation is enhanced [60, 61, 14] in high-symmetry models. In addition, $|\langle T^{z}\rangle|=\sum_{i}{|\langle T_{i}^{z}\rangle|}/N$ is always zero within our error bar, suggesting that the U(1) symmetry is not broken.

Now we have demonstrated that there is a stable long-range spin-orbital order in the YOJ model. Next, we will show the pattern of such an order. In Fig. 4 the spin-orbital order corresponding to the three components of $\bm{X}_{i}$ is plotted as a vector. It shows clearly that the ground state exhibits a stripe long-range order. From Fig. 4(a) we can see that along the dashed bonds, where the interaction is given in Eq. (2), the spins are ferromagnetic, and along the two solid bonds the spins are antiferromagnetic. However, we want to point out that the spins can be ferromagnetic along one of two solid bonds, i.e., (b) and (c) are also the ground-state configurations. This can be understood as follows. First, let us see the equilateral triangle in Fig. 4 with its three sides $\alpha$, $\beta$, and $\gamma$. The lines $\alpha$ and $\beta$ join at the site, say, $m$. Lines $\alpha$ and $\gamma$ join at site $n$. After the transformation $T_{m}^{x}\rightarrow-T_{m}^{x},T_{m}^{y}\rightarrow-T_{m}^{y}$, $T_{n}^{x}\rightarrow-T_{n}^{x},T_{n}^{y}\rightarrow-T_{n}^{y}$, the side $\beta$ becomes solid and side $\gamma$ becomes dashed. By a similar operation, all the bonds parallel to side $\beta$ can become solid, and all bonds parallel to side $\gamma$ can be dashed. Keep in mind that our order parameters do not depend on $T_{i}^{x}$ and $T_{i}^{y}$, which means their expectation values do not change. We then rotate the lattice counterclockwise by $2\pi/3$ around an axis perpendicular to the plane at any lattice site, and we obtain the figure in panel (b). Similarly, we can obtain the configuration in panel (c).

Convergence analysis. In order to obtain reliable data, we not only pushed the bond dimension $D$ as large as possible, but also checked the convergence of the expectation values with respect to the environment tensor dimension $\chi$ for each $D$ as in Refs. [53, 32, 24]. For example, in Fig. 5 we plot $m^{X}$ as a function of $\chi$ for several $D$s. It shows that the data converge quickly as $\chi$ increases, which demonstrates that $\chi\sim D^{2}$ is sufficiently large to produce a reliable result for all $D$s up to 18. Therefore, with $\chi\geq{D^{2}}$ our results should be reliable.

Conclusion. In summary, by using tensor-network algorithms with PESS wavefunction ansatz [49, 50], we have studied the ground state of the spin-orbit-coupled YOJ model [16] on the triangular lattice, which possesses a SU(2)$\times$SU(2)$\times$U(1) symmetry. The trial wavefunction was optimized by the imaginary-time evolution method, and the expectation values were estimated by the multi-sublattice CTMRG algorithm in combination with the nested tensor-network technique. We found that the two SU(2) symmetries are broken, leading to long-range spin-orbital orders with a stripe pattern. The origin of these orders is different from the conventional magnetic order. A careful finite bond-dimension scaling analysis gives the magnitudes of the spin-orbital orders $m^{X}=m^{Y}=0.085(10)$. The reason for $m^{X}=m^{Y}$ is also discussed. Our resuls impose a strong constraint on the microscopic Hamiltonian in searching for quantum spin-orbital liquid in real materials on the triangular lattice.

We are grateful to Hong-Hao Tu and Shaojin Qin for helpful discussions and Ken Chen for plotting Fig. 1. This work is supported by the National key R$\&$D Program of China (Grants No. 2022YFA1402704, No. 2017YFA0302900 and No. 2016YFA0300503), by the National Natural Science Foundation of China (Grants Nos. 12274187, 12274458, 12047501, 11874188, 11834005, 11774420), and by the Research Funds of Renmin University of China (Grants No. 20XNLG19).