Quantum liquid crystals in the finite-field K$\Gamma$ model for $\alpha$-RuCl₃

Introduction. — The Kitaev model Kitaev (2006) is a prominent example of exactly solvable models with a quantum spin liquid (QSL) property, which is characterized by various topological signatures O’Brien et al. (2016); Yamada (2021). However, it is known that real materials cannot fully be described by this fine-tuned model and include additional interactions Wang et al. (2020), which potentially change the topological property of the ground state. Recently the K$\Gamma$ model Rau et al. (2014); Gohlke et al. (2018); Catuneanu et al. (2018); Liu and Normand (2018); Samarakoon et al. (2018); Gordon et al. (2019); Chern et al. (2020); Lee et al. (2020); Yamada et al. (2020); Gohlke et al. (2020); Buessen and Kim (2021); Luo et al. (2021); Zhang et al. (2021) has been investigated intensively in connection with a Kitaev material $\alpha$-RuCl₃, and various kinds of topological transitions were expected theoretically. Whereas various numerical methods have been used, most studies treat itinerant Majorana fermions of the Kitaev model indirectly, and the role of Majorana fermions is unclear.

Among many possible realizations of the proximate Kitaev model Jackeli and Khaliullin (2009); Plumb et al. (2014); Yamada et al. (2017a, b); Liu and Khaliullin (2018); Sano et al. (2018); Jang et al. (2019), $\alpha$-RuCl₃ is an outstanding material because a half-integer quantization of the thermal Hall conductivity has been observed Kasahara et al. (2018); Yokoi et al. (2021). At higher field, the breaking of the threefold rotation symmetry has been observed from the field angle dependence of the heat capacity at the same time as the disappearance of the thermal Hall effect Tanaka et al. (2020). The coincidence is attributed to the first-order transition from Kitaev’s $B$ phase with a Chern number 1 to Kitaev’s $A$ phase with a Chern number 0, which occurs directly maintaining an energy gap Takahashi et al. (2021). A theoretical characterization is necessary for the high-field nematic phase based on a microscopic model like the K$\Gamma$ model. While some studies Gohlke et al. (2018); Lee et al. (2020) find the nematic phase within the K$\Gamma$ model, the property of the associated nematic phase transition from the Kitaev spin liquid is not well discussed Gohlke et al. (2020).

Another interesting phase with a broken threefold rotation symmetry is a Kekulé phase. Though both a nematic phase and a Kekulé phase are $Z_{2}$ spin liquids with a toric code topological order, some properties are different. The most significant difference is the presence of a Majorana zero mode (MZM) in the $Z_{3}$ vortex for a Kekulé phase Jackiw and Rossi (1981); Yang et al. (2019). Here a $Z_{3}$ vortex means a crossing point of boundaries of three different domains of the $Z_{3}$-broken Kekulé phases. MZMs can potentially be braided by the domain wall motion, which leads to a nontrivial phase. Thus, by controlling such topological defects in this liquid crystalline phase, the protected quantum computation in the MZM states is possible.

Both of these phases are relevant to the K$\Gamma$ model even in the mean-field level, and thus we investigate these phases precisely. From now on we call the nematic phase nematic Kitaev spin liquid (NKSL) and the Kekulé phase Kekulé-Kitaev spin liquid (KKSL).

From the mean-field approximation and the exact diagonalization, we discover a vast region of NKSL and KKSL with a broken threefold rotation symmetry. Both of these phases can be regarded as quantum spin analogues of liquid crystals Kivelson et al. (1998). These phases are enabled to be discovered by the methods which treat itinerant Majorana fermions directly.

In this Letter, we study a phase diagram of the finite-field K$\Gamma$ model using the third-order perturbation. The effective model is an interacting Majorana model with four-body and six-body interactions Bravyi et al. (2010); Vijay et al. (2015); Affleck et al. (2017); Li and Franz (2018); Kamiya et al. (2018); Wamer and Affleck (2018); Rahmani et al. (2019); Rahmani and Franz (2019); Li et al. (2019); Tummuru et al. (2021). Based on the mean-field solution and the exact diagonalization, we find the phase diagram is rich enough to predict various exotic spin liquids like NKSL and KKSL.

Third-order perturbation. — We begin with the following Kitaev-$\Gamma$ or K$\Gamma$ model. These two terms are known to be dominant in $\alpha$-RuCl₃.

where we assume $K<0$, and $\Gamma$ to be a real number, and $S_{i}$ and $\sigma_{i}$ are spin-1/2 and Pauli operators defined on the $i$th site, respectively. $\langle ij\rangle\in\alpha\beta(\gamma)$ means that a nearest-neighbor (NN) bond $\langle ij\rangle$ belongs to $\gamma$-bonds as shown in Fig. 1(a), and the $\alpha\beta$-plane is perpendicular to the $\gamma$-direction. We note that our study is based on the original Kitaev model with an additional interaction, not the so-called Kekulé-Kitaev model Kamfor et al. (2010); Quinn et al. (2015); Mirmojarabian et al. (2020). The Kekulé-Kitaev model explicitly breaks the $Z_{3}$ symmetry, while our model spontaneously breaks this symmetry in the Kekulé phase.

As is usually the case, the nontrivial contribution begins from the third order in $\Gamma/K$. The third-order perturbation leads to the following effective Hamiltonian.

where $\langle ij\rangle\langle kl\rangle\langle mn\rangle\in\hexagon$ means that $\langle ij\rangle$, $\langle kl\rangle$, and $\langle mn\rangle$ are two possible patterns for $z$-, $x$-, and $y$-bonds, respectively, inside each hexagon. The factor $0.21\sim 4\times 0.227^{2}$ comes from the intermediate vortex pattern Yamada (2020); Takahashi et al. (2021), as shown in Fig. 1(b).

We then move on to the effective Majorana interacting model.

where $t\propto|K|+O(\Gamma^{2}/|K|)$, $\langle ijkl\rangle$ represents armchair-type interactions shown in Fig. 1(c), and $\langle ijklmn\rangle$ represents six-body interactions shown in Fig. 1(d). The site ordering of $\langle ijkl\rangle$ and $\langle ijklmn\rangle$ is counterclockwise, as shown in Figs. 1(c)-(d).

The same process applies to the case with a magnetic field. We can consider the Hamiltonian $H=H_{0}+H_{1}$ with

where $\vec{h}=(h^{x},h^{y},h^{z})^{t}$ is an applied magnetic field. By tracking intermediate vortex states, we would finally get

where $t\propto|K|+O(\Gamma^{2}/|K|)+O(h^{2}/|K|)$, $\langle ijkl\rangle\in\alpha\beta(\gamma)$ means that a bond $\langle jk\rangle$ belongs to $\gamma$-bonds with the $\alpha\beta$-plane perpendicular to the $\gamma$-direction, $\langle\!\langle ij\rangle\!\rangle$ represents a next-nearest-neighbor (NNN) bond, and $\langle\!\langle ijkl\rangle\!\rangle$ represents a $\Yup$-shaped, or $\Ydown$-shaped interaction for sites $ijkl$, where $ij$, $jk$, and $il$ are connected by $x$-, $y$-, and $z$-bonds, respectively, as shown in Fig 1(e). $\langle\!\langle\!\langle ij\rangle\!\rangle\!\rangle\in\alpha\beta(\gamma)$ means the next-next-nearest-neighbor (NNNN) bond which is a diagonal of hexagons in the $\gamma$-direction, where $i$ is even and $j$ is odd as usual.

From now on, we set $h^{x}=h^{y}=h^{z}=h/\sqrt{3}$ for simplicity. The absolute value of $t$ is undecidable, so we set $t=1$ as a simple normalization. We note that two-body terms are not antisymmetrized, so the hopping amplitude is $1/8$ in the Hermitian form.

The $tqs$ model. — Let us begin with the zero-field model. First, we extend the model to the following form:

with $q$ and $s$ being real parameters. When $q=s=6\Gamma^{3}/(0.21|K|^{2})$, this $tqs$ model becomes the original (K$\Gamma$) model with $h=0$.

The operator

can be defined for each hexagon plaquette $p=\langle ijklmn\rangle$, and plays an important role as they commute with each other. Indeed, the $tqs$ model is integrable in the limit $|s|\to\infty$. Eigenstates are constructed by assigning $V_{p}=\pm 1$ for each plaquette.

This limit ($|s|\to\infty$) is also known as the Vijay-Hsieh-Fu (VHF) surface code Vijay et al. (2015) and has a characteristic $Z_{2}$ topological order with an exact $S_{3}$ anyon symmetry. To check the existence of such a phase we did exact diagonalization of the 54-site cluster, which is shown in Fig. 2(a). This VHF surface code is indeed realized in the $tqs$ model with a large $|s|$ region, as shown in Fig. 2(b). In Fig. 2(b) the right half is mostly blue connected to the $s\to\infty$ limit, and the left half is mostly red connected to the $s\to\infty$ limit. This means that the large $|s|$ phase is rather stable with respect to the perturbation of $t$ and $q$.

It seems that even on the $q=s$ line, most of the $q>0$ region is connected to this VHF surface code phase, and from this we conclude that even in the original (K$\Gamma$) model Eq. (5) this phase is realized in the large $\Gamma/|K|$ region. Since the resolution of Fig. 2(b) is not good enough, the exact region of $\Gamma/|K|$ where this phase is realized will be discussed later.

K$\Gamma$ model with a magnetic field. — Then, we go back to the original (K$\Gamma$) model Eq. (5). With a magnetic field, it indeed becomes easy to identify the phase because we can use a Chern number to diagnose gapped topological phases. The Chern number is connected to the topological order of the spin model according to Kitaev’s 16-fold way. For example, phases with a Chern number $\pm 1$ in the Majorana model can be identified with the Ising topological order, which we call a chiral spin liquid (CSL) with a Chern number $\pm 1$ in the following.

The mean-field phase diagram is shown in Fig. 3(a) Phases are identified from the Chern number and bond order parameters. The Chern number is computed by the Fukui-Hatsugai-Suzuki method Fukui et al. (2005). The details of the calculation is included in Supplemental Material (SM) SM . The method is based on Refs. Li and Franz (2018); Takahashi et al. (2021). From the phase diagram we find the following phases:

1. 1.

   a KKSL with a Chern number $0$.

2. 2.

   Kitaev’s non-Abelian CSL with a Chern number $1$.

3. 3.

   an Abelian CSL with a Chern number $-2$.

4. 4.

   an NKSL with a Chern number $0$.

5. 5.

   a zigzag nematic phase with a Chern number $-1$.

The bond ordering patterns for a KKSL, an NKSL, and a zigzag nematic phase are shown in Figs. 1(f), (g), and (h), respectively.

Exact diagonalization. — From now on we will confirm the results of the mean-field theory by exact diagonalization of Eq. (5). The exact diagonalization is done using the 54-site cluster again. We note that this size is much larger than what has been used for the Kitaev model exact diagonalization. We take an average of expectation values for all the degenerate ground states. First, we check the operator $\{\psi,\psi^{\dagger}\}/2$, corresponding to the absolute value of the Kekulé order parameter Pujari et al. (2015) with

where the definition of $\Delta_{i}^{\prime}$ ($i=1$–$9$) is included in SM SM . In Fig. 3(b), the orange triangular region with $\Gamma<0$ corresponds to the Kekulé-ordered phase. The area of this phase becomes larger than that of the mean-field phase diagram. Surprisingly, it looks like the Kekulé order appears with an infinitesimal $\Gamma<0$ in the 54-site calculation. This implies that the Kitaev spin liquid is unstable with respect to a negative $\Gamma$ perturbatively and it is possible that the Kekulé order opens a gap for Dirac cones (see SM SM for more details).

Next we check the nematic order parameter $\{\phi,\phi^{\dagger}\}/2$ Pujari et al. (2015) with

where the definition of $\Delta_{i}$ ($i=1$–$9$) is included in SM SM . In Fig. 3(c), NKSL is identified with the small orange region around $\Gamma/|K|=0.1$. The area of this phase becomes smaller than that of the mean-field phase diagram. We now confirmed the existence of NKSL even in the exact diagonalization, but the zigzag nematic phase disappears in the same calculations. There remains a violet region which is different from the zigzag nematic phase found in the mean-field approach. We note that the order parameter is nonzero on the whole parameter space in Fig. 3(b) and (c) due to the finite-size effect.

Finally, we check the value of $V_{p}$ to identify the mysterious area which cannot be identified by the calculation of the nematic order parameter. As shown in Fig. 3(d), the value of $|V_{p}|$ reaches around 0.75, which is very close to the highest value $+1$. Combined with what is obtained from the $tqs$ model, we believe that this phase is connected to the large $|s|$ phase of the $tqs$ model. Thus, this phase appearing in the large $\Gamma/|K|$ region should be described by the VHF surface code. This region indeed has the largest value of the entanglement entropy (see Fig. 4), which reconfirms that it is a highly entangled phase which cannot be described by the mean-field picture.

Phase diagram of the $K\Gamma$ model. — As shown in Fig. 4, we find the following phases.

1. 1.

   a KKSL with a Chern number $0$.

2. 2.

   Kitaev’s non-Abelian CSL with a Chern number $1$.

3. 3.

   an Abelian CSL with a Chern number $-2$.

4. 4.

   an NKSL with a Chern number $0$.

5. 5.

   a VHF surface code.

Except for the fact that the zigzag nematic phase is replaced by the VHF surface code, the phases obtained are the same as those found in the mean-field theory. We note that Kitaev’s non-Abelian CSL with a Chern number $1$ is identified because this phase is stable on the $\Gamma/|K|=0$ line Takahashi et al. (2021), while an Abelian CSL with a Chern number $-2$ is just a speculation.

Especially in the experimentally relevant region, $\Gamma/|K|\sim 0.1$, we find an NKSL phase with a Chern number 0, which potentially explains the observed topological nematic transition from Kitaev’s non-Abelian CSL to the toric code phase.

Overall, we have obtained a very rich phase diagram within the flux-free assumption, and disclosed the existence and properties of NKSL expected from experiments. Though our calculation is somewhat idealized, a more realistic phase diagram should be obtained by dealing with a $\Gamma^{\prime}$ term, for example, additionally Takikawa and Fujimoto (2019, 2020). In addition, in the future it is important to check whether an Abelian CSL with a Chern number $-2$ is connected to the one discovered previously in the K$\Gamma^{\prime}$ model Takikawa and Fujimoto (2020).

Effect of the Heisenberg term. — Here we will discuss the effect of the Heisenberg term. With this term the third-order perturbation leads to the following additional interaction.

where $\langle\!\langle\!\langle ijkl\rangle\!\rangle\!\rangle$ represents an zigzag-type interaction for sites $ijkl$, and $\langle\!\langle\!\langle ij\rangle\!\rangle\!\rangle^{\prime}$ represents an NNNN hopping which is not a diagonal of hexagons. This term promotes nematic order, and potentially enhances NKSL. Thus, the additional small Heisenberg term should be relevant to the field-induced nematic transition, and will be investigated in the future.

Discussion. — We invent a systematic approach to treat non-Kitaev interactions in a perturbative manner by fully incorporating many-body effects of Majorana interactions. We discover the emergence of a Kekulé order in the K$\Gamma$ model for $\Gamma<0$, and a nematic order for $\Gamma>0$, as well as a more exotic VHF surface code phase. We also discover an Abelian CSL with a Chern number -2, which can be detected by the sign change of the thermal Hall effect or spin Seebeck effect Takikawa et al. (2021).

In the experimentally relevant region $\Gamma/|K|\sim 0.1$ for $\alpha$-RuCl₃, the discovered nematic phase potentially explains the nematic transition observed in the high-field phase. While the position of the NKSL phase $\Gamma/|K|\sim 0.1$ is consistent the density matrix renormalization group Gohlke et al. (2018, 2020), its nature is totally different from nematic phases discovered previously Lee et al. (2020). NKSL is not a simple nematic phase but a spin liquid phase with a partial symmetry breaking. We believe that NKSL discovered in our study is not the same phase as the ones in the previous studies, and that the approximate coincidence of the nematic region $\Gamma/|K|\sim 0.1$ is only accidental. Our study endorses the topological nematic transition scenario proposed in Ref. Takahashi et al. (2021) with a more realistic model and interactions. This state is potentially detectable by nuclear magnetic resonance (NMR) or Mössbauer spectroscopy Yamada and Fujimoto (2021).

Another direction to explore is to engineer the Kekulé order in artificial systems. The sign and magnitude of the $\Gamma$ term is known to be strongly dependent on the spin-orbit coupling (SOC), so it can potentially be controlled by the proximity effect with a substrate. The substrate with heavy elements may strengthen the SOC of $\alpha$-RuCl₃, and allow the fine tuning of $\Gamma$. The signature of the Kekulé phase is the existence of an MZM at the $Z_{3}$ vortex core, which would potentially be discovered by scanning tunneling microscope (STM), leading to a possible topological quantum computation.