# Emergent Quadrupolar Order in the Spin-$1/2$ Kitaev-Heisenberg Model

Manodip Routh  S.N. Bose National Centre for Basic Sciences, Kolkata 700098,
India.  Sayan Ghosh  S.N. Bose National Centre for Basic Sciences, Kolkata
700098, India.  Jeroen van den Brink  Institute for Theoretical Solid State
Physics, IFW Dresden, 01069 Dresden, Germany.  Würzburg-Dresden Cluster of
Excellence ct.qmat, Germany Department of Physics, Technical University
Dresden, 01069 Dresden, Germany.  Satoshi Nishimoto  s.nishimoto@ifw-
dresden.de Institute for Theoretical Solid State Physics, IFW Dresden, 01069
Dresden, Germany.  Department of Physics, Technical University Dresden, 01069
Dresden, Germany.  Manoranjan Kumar<EMAIL_ADDRESS>S.N. Bose
National Centre for Basic Sciences, Kolkata 700098, India.

###### Abstract

Motivated by the largely unexplored domain of multi-polar ordered spin states
in the Kitaev-Heisenberg (KH) systems we investigate the ground state dynamics
of the spin-$\frac{1}{2}$ KH model, focusing on quadrupolar (QP) order in
2-leg ladder and two-dimensional honeycomb lattice geometries. Employing exact
diagonalization and density-matrix renormalization group methods, we analyze
the QP order parameter and correlation functions. Our findings reveal a robust
QP order across a wide range of the phase diagram, influenced by the interplay
between Heisenberg and Kitaev interactions. Notably, we observe an enhancement
of QP order near Kitaev quantum spin liquid (QSL) phases, despite the absence
of long-range spin-spin correlations. This highlights a complex relationship
between QP order and QSLs, offering new insights into quantum magnetism in
low-dimensional systems. Our findings provide a rational explanation for the
observed nonlinear magnetic susceptibility in $\alpha$-RuCl3.

††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

Introduction— Quantum spin liquids (QSLs) represent a quantum phase of matter
that challenges traditional understanding by lacking magnetic long-range order
(LRO), yet they may exhibit global topological order Balents (2010). These
elusive states emerge in frustrated magnetic systems, where competing spin
exchange interactions not only foster the QSL state but also lead to complex
entangled states such as dimer formations Chhajlany  _et al._ (2007),
fractional excitations like spinons Enderle  _et al._ (2010); Schlappa  _et
al._ (2018), and Majorana excitations Do  _et al._ (2017a); Koga  _et al._
(2021); Wulferding  _et al._ (2020); Hashimoto  _et al._ (2022) in the ground
state (gs). A key development in understanding QSLs was the Kitaev model,
which proposed spin-$1/2$ particles on a hexagonal lattice with direction-
dependent exchange interactions, introducing frustration into the system
Kitaev (2006). The gs of this model is a QSL, characterized by Majorana
fermions and gauge vortices as elementary excitations Hermanns  _et al._
(2015); Le Hur  _et al._ (2017); Mandal  _et al._ (2011); Baskaran  _et al._
(2007, 2008); Baskaran (2023).

The experimental pursuit of the Kitaev model has prominently featured Iridium-
based materials such as A2IrO3 (A=Na, Li), where Ir atoms form a honeycomb
structure Singh  _et al._ (2012); Hwan Chun  _et al._ (2015). In these
materials, the octahedral coordination of ligands induces a crystal field,
splitting the Ir $d$-orbitals into $t_{2g}$ and $e_{g}$ levels. Strong spin-
orbit coupling further splits the $t_{2g}$ orbitals, allowing the two highest
states of each to mimic effective spin-$1/2$ degrees of freedom. This results
in orbital-dependent, directionally anisotropic spin exchanges, i.e., so-
called Kitaev interactions Trebst and Hickey (2022a). However, the Kitaev-
Heisenberg (KH) model, incorporating residual nearest-neighbor Heisenberg
interactions, offers a more realistic description of the magnetic properties
of these materials Chaloupka  _et al._ (2010). Similar Hamiltonians have been
proposed for materials like $\alpha$-RuCl3 Suzuki  _et al._ (2021); Do  _et
al._ (2017b); Banerjee  _et al._ (2017); Trebst and Hickey (2022b); Winter
_et al._ (2016); Janssen  _et al._ (2016), $\beta$-Li2IrO3 Katukuri  _et al._
(2016); Takayama  _et al._ (2015); Trebst and Hickey (2022b); Winter  _et al._
(2016); Janssen  _et al._ (2016), and $\gamma$-Li2IrO3 Tsirlin and Gegenwart
(2022); Trebst and Hickey (2022b).

The gs phase diagram of the honeycomb-lattice KH model has been extensively
studied, revealing four magnetic LRO phases (Néel, zigzag (ZZ), stripy (ST),
and ferromagnetic (FM)) and two QSL phases; FM Kitaev (FK) and
antiferromagnetic Kitaev (AFK) QSLs by tuning the Heisenberg and Kitaev
interactions Chaloupka  _et al._ (2010); Gotfryd  _et al._ (2017); Schaffer
_et al._ (2012); Feng  _et al._ (2007). However, the potential for a quantum
phase involving multi-polar spin states, such as a quadrupolar (QP) or spin-
nematic state, in the KH model remains largely unexplored. From an
experimental perspective, third-order positive susceptibility probing a
higher-order correlation such as QP has been observed at zero magnetic field
in $\alpha$-RuCl3 Holleis  _et al._ (2021). Yet, the intricate mechanisms
behind these observations are not fully elucidated. In addition to these
experimental findings, the present study is further motivated by the inclusion
of a staggered QP operator in the Hamiltonian of the KH model (see below).

The QP phase, established in liquid crystal systems Kelker (1973, 1988);
Mermin and Wagner (1966); Alexander  _et al._ (2012), manifests in magnetic
systems as spin nematics, where magnetic quadrupole moments create
orientational order without magnetic LRO Läuchli  _et al._ (2006); Tsunetsugu
and Arikawa (2006); Lacroix  _et al._ (2011). Initially identified in spin-1
systems with biquadratic exchange Blume and Hsieh (1969), this order has been
explored in various Heisenberg spin-1 models Hu  _et al._ (2019); Podolsky and
Demler (2005); Tanaka and Hotta (2020). Extending this concept to spin-$1/2$
systems has spurred research, particularly in $J_{1}$$-$$J_{2}$ chain systems
under strong magnetic fields Hikihara  _et al._ (2008); Parvej and Kumar
(2017); Sudan  _et al._ (2009). In the context of Kitaev systems, a spin-
nematic phase in the spin-$1/2$ Kitaev-Ising model has been predicted, arising
from the interplay between Kitaev QSL and magnetic LRO, but lacking
topological order Nasu  _et al._ (2017). Additionally, four-body interactions
among Majorana fermions in the Kitaev QSL could induce a topological nematic
phase transition, evolving from the chiral QSL phase to the toric code phase
Takahashi  _et al._ (2021).

In this letter, we explore the gs dynamics of the spin-$\frac{1}{2}$ KH model,
focusing on the QP order. We investigate two geometries: the 2-leg ladder and
the two-dimensional (2D) honeycomb lattice, using exact diagonalization (ED)
and density-matrix renormalization group (DMRG) methods. Our primary
contribution is a detailed exploration of the QP order parameter and
correlation functions within these systems. We present a novel gs phase
diagram highlighting the QP order as a function of the ratio between
Heisenberg and Kitaev interactions. A striking aspect of our findings is the
robust presence and stability of QP order across a wide range of the phase
diagram, with a notable enhancement near the Kitaev QSL phases. This
enhancement occurs despite the absence of longer-range spin-spin correlations,
adding a new dimension to our understanding of the interplay between QP order
and QSLs.

Quadrupolar operator— The QP state emerges as a quantum ordered phase
characterized by collective bimagnon excitations. This intriguing state, while
exhibiting a preference for the orientation of paired magnons, notably does
not break time-reversal symmetry. Its detection hinges on the analysis of a
symmetric and traceless rank-2 tensor operator Penc and Läuchli (2011):

$\displaystyle\hat{Q}_{ij}^{\alpha\beta}=S_{i}^{\alpha}S_{j}^{\beta}+S_{j}^{\alpha}S_{i}^{\beta}-\frac{2}{3}(S_{i}.S_{j})\delta_{\alpha\beta}$
(1)

where $\alpha$ and $\beta$ represent the Cartesian coordinates such as $x$,
$y$, and $z$. As explained below, the Hamiltonian of the KH model can be
reformulated to explicitly include a specific component of Eq. (1):

$\displaystyle\hat{Q}^{x^{2}-y^{2}}_{ij}=\frac{1}{2}(S^{+}_{i,j}S^{+}_{i+1,j}+S^{-}_{i,j}S^{-}_{i+1,j}).$
(2)

Figure 1: Lattice structures used in our calculations: (a) KH ladder or
brickwall lattice which is obtained by deforming the honeycomb-lattice KH
model. (b) 2-leg KH ladder. (c) 24-site honeycomb-lattice KH cluster with
periodic boundary conditions.

Model— We consider a spin-$1/2$ KH model on both a ladder and a honeycomb
lattice Chaloupka  _et al._ (2010). The Hamiltonian is given by

$\displaystyle\mathcal{H}=J\sum_{i,j}\vec{S}_{i}\cdot\vec{S}_{j}+K\sum_{\gamma}\sum_{i,j}S_{i}^{\gamma}S_{j}^{\gamma}.$
(3)

Where $S_{i}^{\gamma}$ is the $\gamma$ ($=x$, $y$, or $z$) component of the
spin-$1/2$ operator $\vec{S}_{i}$ at site $i$. Here, $J$ and $K$ denote the
Heisenberg and Kitaev interactions, respectively. To facilitate analysis, we
introduce an angular parameter $\phi$ ($\in[0,2\pi]$), defining $J=\cos\phi$
and $K=\sin\phi$.

For the ladder or cylinder geometry, the Hamiltonian is adapted to include the
spin-$1/2$ operator $\vec{S}{i,j}$ at site $i$ on the $j$-th leg:

$\displaystyle\mathcal{H}_{\rm
leg}=\frac{2J+K}{4}\sum_{j=1}^{n}\sum_{i=1}^{L}(S_{i,j}^{+}S_{i+1,j}^{-}+S_{i,j}^{-}S_{i+1,j}^{+})$
$\displaystyle+\frac{K}{2}\sum_{j=1}^{n}\sum_{i=1}^{L}(-1)^{i+j}(\hat{Q}^{x^{2}-y^{2}}_{ij})+J\sum_{j=1}^{n}\sum_{i=1}^{L}(S_{i,j}^{z}S_{i+1,j}^{z})$
$\displaystyle\equiv\mathcal{H}_{\rm exc}+\mathcal{H}_{\rm Q}+\mathcal{H}_{\rm
z}.$ (4)

In this formulation, $\mathcal{H}_{\rm Q}$ signifies the QP operator. Its
dominance suggests the potential stabilization of QP order. The possibility of
QP order arising through similar mechanisms has been previously discussed in
the context of $J_{1}$-$J_{2}$ chain systems under magnetic fields Grafe  _et
al._ (2017).

Although the Kitaev model was originally proposed for the honeycomb lattice,
it is recognized that any three-coordinated lattice with Kitaev interactions
can exhibit similar intriguing properties. In this sense, the KH model on a
2-leg ladder [Fig. 1(b)], derived from a brick-wall lattice analogous to the
honeycomb lattice, aligns with the necessary geometry. This alignment
indicates that the 2-leg KH ladder can offer vital insights into the behavior
of the 2D honeycomb-lattice KH model. Notably, the gs phase diagram of the
ladder model displays significant similarities with the 2D counterpart in
terms of magnetic order and Kitaev QSL, with the key distinction being the
replacement of the Néel phase by the rung-singlet (RS) phase Agrapidis  _et
al._ (2019); Catuneanu  _et al._ (2019). Therefore, initiating our
investigation with the 2-leg KH ladder is a strategic choice. This setup not
only allows for precise numerical calculations but also facilitates reliable
extrapolations to the thermodynamic limit, essential for exploring novel
phenomena in the 2D KH model. Building upon this, we extend our investigation
to 4-leg and 6-leg cylinders [Fig. 1(a)], progressively approaching the 2D
bulk limit. This systematic increase in the number of legs serves to bridge
our understanding from the simpler ladder systems to the more complex dynamics
inherent in the full 2D KH model.

Methods— To investigate the intricate dynamics of the KH model, we employ two
complementary computational techniques. For the 2-leg ladder system [Fig.
1(b)], encompassing system sizes up to $12\times 2$, and for the 24-site
periodic boundary condition (PBC) cluster [Fig. 1(c)], we utilize the exact
diagonalization (ED) method. For larger systems, specifically extended 2-leg
ladders, as well as 4-leg and 6-leg cylinders [Fig. 1(a)], we implement the
DMRG method. In our DMRG calculations, we retain up to 5000 density-matrix
eigenstates during the renormalization process. This high number of kept
states ensures the accuracy of our results, with the largest discarded weight
being maintained at an impressively low level of approximately $2\times
10^{-5}$.

Figure 2: Results for 2-leg KH ladder. (a) Absolute values of the QP order
parameter for leg bonds with various system sizes. The inset represents the
same physical quantity but without taking absolute values. (b) Ground-state
Phase diagram for magnetic order (inner) and QP order (outer). Four long-range
ordered phases: rung-singlet (RS), zigzag (ZZ), ferromagnetic (FM), stripy
(ST) and two QSL phases: antiferromagnetic Kitaev (AFK) and ferromagnetic
Kitaev (FK). (c) Compilation of the behavior of the QP correlation functions
$|P_{\rm L1L1}(r)|$ with $N=100\times 2$ cluster for representative parameters
(see text). (d) Inter-leg and rung-rung QP correlation functions at
$\phi/\pi=0.5$, $1.2$, and $1.5$. Inset: Examples for the cases where the QP
correlation function exhibits an exponential decay.

Quadrupolar order in the 2-leg ladder— We commence our analysis with the 2-leg
ladder system. To obtain a comprehensive phase diagram for QP order as a
function of $\phi$, we calculate the QP order parameter
$\langle\hat{Q}^{x^{2}-y^{2}}_{ij}\rangle$ using ED. Our findings, focusing on
a single bond along the leg, are depicted in Fig.2(a). We find that
$\langle\hat{Q}^{x^{2}-y^{2}}_{ij}\rangle$ exhibits finite values over a broad
range of $\phi$ ($0.48\pi\lesssim\phi\lesssim 1.56\pi$) with minimal finite-
size effects. This order parameter parameter shows significant values in the
AFK, FK, and ZZ phases of the magnetic phase diagram, while it rather
diminishes in the FM phase [Fig.2(b)]. Conversely, when $i$ and $j$ are
positioned on a rung bond, $\langle\hat{Q}^{x^{2}-y^{2}}_{ij}\rangle$
consistently equals zero.

To further substantiate these findings, using DMRG we compute leg-leg and
rung-rung QP correlation functions, defined as:

$\displaystyle P_{\rm LmLn}(r)=\langle
S_{i,m}^{+}S_{i+1,m}^{+}S_{i+r,n}^{-}S_{i+1+r,n}^{-}\rangle$ (5)

and

$\displaystyle P_{\rm RR}(r)=\langle
S_{i,m}^{+}S_{i,m+1}^{+}S_{i+r,m}^{-}S_{i+r,m+1}^{-}\rangle,$ (6)

respectively. The results, summarized in Fig.2(c,d), reveal that the intra-leg
correlation $P_{\rm L1L1}(r)$ exhibits QP LRO within the $\phi$ range where
$\langle\hat{Q}^{x^{2}-y^{2}}_{ij}\rangle$ is finite. Note that $P_{\rm
L1L1}(r)$ alternates in sign with distance. Meanwhile, near both Kitaev points
$\phi=\pm\pi/2$, the inter-leg correlation $P_{\rm L1L2}(r)$, and near the FK
phase, the rung-rung correlation $P_{\rm RR}(r)$, also exhibit QP LRO. Unlike
$P_{\rm L1L1}(r)$, both $P_{\rm L1L2}(r)$ and $P_{\rm RR}(r)$ are always
positive. Interestingly, despite the rung QP order parameter being
consistently zero across all $\phi$ values, $P_{\rm RR}(r)$ still exhibits LRO
in the FK phase. This suggests that states of QP order with broken spin-
rotational symmetry are degenerate, leading to a zero local order parameter,
yet detectable through correlation functions.

From these correlation function analyses, the $\phi$ domain for QP order can
be classified into four phases: QP1 (only $P_{\rm L1L1}(r)$ shows LRO), QP2
(both $P_{\rm L1L1}(r)$ and $P_{\rm L1L2}(r)$ show LRO), QP3 (all of $P_{\rm
L1L1}(r)$, $P_{\rm L1L2}(r)$, and $P_{\rm RR}(r)$ show LRO), and a phase with
no QP LRO. The schematic picture of each QP phase is illustrated in the
Supplementary Material sup . The correspondence with the magnetic order phases
is shown in the phase diagram in Fig.2(b). The phase diagram for QP order can
be reasonably explained by considering the spin structures of the
corresponding magnetic phases. In the FM phase ($0.81\lesssim\phi/\pi\lesssim
1.38$), the near-full magnetization results in smaller values of $P_{\rm
L1L1}(r)$ and $|\langle\hat{Q}^{x^{2}-y^{2}}_{ij}\rangle|$, and they vanish
without $\mathcal{H}_{\rm Q}$ at $\phi=\pi$. The ZZ phase
($0.53\lesssim\phi/\pi\lesssim 0.81$) exhibits Ising-type LRO in each leg with
opposite polarization between legs Agrapidis  _et al._ (2019), leading to the
anticipated development of $P_{\rm L1L2}(r)$. However, $P_{\rm RR}(r)$ does
not show LRO due to antiparallel spin orientation on the rungs. In fact, a
detailed wave-function analysis suggests the dominant contribution of double-
flip configurations in the gs (see the Supplementary Material sup ). In the RS
($-0.26\lesssim\phi/\pi\lesssim 0.48$) and ST ($1.57\lesssim\phi/\pi\lesssim
1.74$) phases, the each leg possesses a Néel-type LRO, precluding the
dominance of the $\mathcal{H}_{\rm Q}$ term. Thus, a QP order is absent.

In the AFK ($0.48\lesssim\phi/\pi\lesssim 0.53$) and FK
($1.38\lesssim\phi/\pi\lesssim 1.57$) phases, the competition between dipole
fluctuations ($\mathcal{H}_{\rm exc}$) and QP fluctuations ($\mathcal{H}_{\rm
Q}$) complicates gs determination. Nevertheless, the larger coefficient of
$\mathcal{H}_{\rm Q}$ supports the emergence of QP order along the leg, as
corroborated by our numerical results. Remarkably, in the FK phase, the rung
QP correlation function $P_{\rm RR}(r)$ also exhibits LRO, with a saturation
value exceeding $P_{\rm L1L1}(r)$. This is attributed to the FM Ising
character of rung interactions, conducive to bimagnon formation Shannon  _et
al._ (2006). The approximate wave function of the FK phase is detailed in the
Supplementary Material sup .

Given these insights, we hypothesize that QP order might also be present in
the 2D honeycomb KH system, akin to the 2-leg ladder scenario. To test this
hypothesis, we extend our investigation to isotropic honeycomb cluster as well
as 4-leg and 6-leg cylinders, in the parts that follow.

Figure 3: Results for honeycomb-lattice KH model. (a) Absolute values of the
QP order parameter for $x$ or $y$ bond, calculated using the 24-site PBC
cluster. For comparison, those for the 2-leg KH ladder are also plotted. The
inset represents the same physical quantity but without taking absolute
values. (b) Ground-state Phase diagram for magnetic order (inner) and QP order
(outer).

Quadrupolar order in the 2D honeycomb lattice— Analogous to the approach for
the 2-leg ladder, we first aim to grasp the overall picture of the QP order
phase diagram as a function of $\phi$ by calculating the QP order parameter
using ED on the 24-site PBC cluster.The results, illustrated in Fig. 3(a), are
compared with those from the 2-leg ladder. In the context of magnetic order
phases for the 2D honeycomb-lattice KH model, the QP order parameter is finite
in the AFK ($0.494\lesssim\phi/\pi\lesssim 0.506$), ZZ
($0.506\lesssim\phi/\pi\lesssim 0.847$), FM ($0.847\lesssim\phi/\pi\lesssim
1.452$), FK ($1.452\lesssim\phi/\pi\lesssim 1.539$), and ST
($1.539\lesssim\phi/\pi\lesssim 1.694$) phases, and it is zero only in the
Néel phase ($-0.306\lesssim\phi/\pi\lesssim 0.494$). While the results are
fundamentally similar to those for the 2-leg ladder, a notable difference is
that the QP order parameter remains finite even in the ST phase. For the
$z$-bonds, similar to the rung QP in the 2-leg ladder, the QP order parameter
is zero.

Figure 4: DMRG results for the QP correlation functions $P_{\rm L1L1}(r)$ at
representative value of $\phi$ for six quantum phases using 2-leg, 4-leg, and
6-leg cylinders.

To ascertain the persistence of QP order in the bulk limit, we compute QP
correlation functions using DMRG for 4-leg and 6-leg cylinders, in addition to
the 2-leg ladder. Focusing on the intra-leg QP correlation function, $P_{\rm
L1L1}(r)$, we observe in Fig. 4 that, consistent with the QP order parameter
results, $P_{\rm L1L1}(r)$ converges to a finite value in the long-distance
limit for all magnetic phases except the Néel phase. Importantly, this
convergence value does not diminish with an increasing number of legs,
suggesting its stability in the bulk limit. In the Néel phase, $P_{\rm
L1L1}(r)$ exhibits a power-law decay indicating no long-range QP order.

Figure 5: DMRG results for the inter-leg QP correlation functions, $P_{\rm
L1Ln}(r)$ for $n=1-3$ at magnetic ordered phases, including ZZ
($\phi/\pi=0.65$), FM ($\phi/\pi=1.2$), ST ($\phi/\pi=1.6$) and QSL phase
($\phi/\pi=0.5$), using 6-leg cylinders.

The decay of QP correlation functions perpendicular to the leg direction,
namely along the $z$-bond, is also of interest. We calculate the inter-leg QP
correlation function, $P_{\rm L1Ln}(r)$, for various leg distances in a
$50\times 6$ system. The calculated $P_{\rm L1Ln}(r)$ for $n=1-3$ is plotted
in Fig. 5. In the ZZ and FM phases, the saturation value of $P_{\rm L1Ln}(r)$
at long distances follows power-law decay [Fig. 5(b,c)], suggesting quasi-LRO
along the $z$-bond direction (QPxy1 phase). Near the Kitaev points (AFK and
FK), a rapid decay of inter-leg correlations is observed [Fig. 5(a)],
indicating LRO in the $xy$ zigzag direction alone (QPxy2 phase). The ST phase
presents a unique case, where the decay along the $z$-bond direction does not
follow a simple pattern [Fig. 5(d)], hinting at the potential for LRO (QP2D
phase). For all $\phi$ values, the rung-rung QP correlation function, $P_{\rm
RR}(r)$, exhibits exponential decay for both 4-leg and 6-leg cylinders. The
correspondence between QP order and magnetic order is summarized in the phase
diagram in Fig. 3(b). A notable distinction from the 2-leg ladder case is the
presence of QP order in the ST phase, likely due to the unique orientation of
ST order in this setting, although the magnitude of QP correlations in this
phase is relatively small compared to other phases.

Summary— We have explored the gs dynamics of the spin-$\frac{1}{2}$ KH model,
focusing particularly on the emergence and characteristics of QP order. Our
study spans two distinct geometries: the 2-leg ladder and the 2D honeycomb
lattice. Employing a synergistic approach combining ED and DMRG techniques, we
have meticulously analyzed the QP order parameter and correlation functions in
these systems. Our key contribution is the unveiling of a comprehensive gs
phase diagram, which distinctly maps out the QP order in relation to the
interplay between Heisenberg and Kitaev interactions. A notable discovery is
the pronounced stability of QP order across a broad range of the phase
diagram, underscoring its robustness in the KH model. Additionally, our
findings reveal a remarkable enhancement of the QP order parameter and its
correlation functions near the Kitaev QSL phases. This enhancement is
particularly intriguing as it occurs in the absence of long-range spin-spin
correlations, suggesting a nuanced relationship between QP order and QSLs. In
fact, the third-order positive susceptibility, indicative of a general
signature of QP order, has been observed in $\alpha$-RuCl3 Holleis  _et al._
(2021), where the relationship between $K$ and $J$ appears to be $-K\gg J>0$
Yadav  _et al._ (2016), suggesting a proximity to the FK QSL phase.

These insights not only advance our understanding of the KH model but also
contribute significantly to the broader discourse on quantum magnetism in low-
dimensional systems, highlighting the intricate interplay between different
types of quantum order.

Acknowledgments$-$ MK thanks Saptarshi Mandal, Krishnendu Sengupta, David
Logan, Jay D Sau, Rajiv R. P. Singh for the fruitful discussions. We thank
Ulrike Nitzsche for technical assistance. MK acknowledges support from SERB
through Grant Sanction No. CRG/2020/000754. MR acknowledges the support from
Institute for Theoretical Solid State Physics of IFW Dresden and DST-India. SG
expresses appreciation for financial support from DST-INSPIRE. This project is
funded by the German Research Foundation (DFG) via the projects A05 of the
Collaborative Research Center SFB 1143 (project-id 247310070) and through the
Wurzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum
Matter ct.qmat (EXC 2147, Project No. 39085490).

## References

  * Balents (2010) L. Balents, Nature 464, 199 (2010).
  * Chhajlany  _et al._ (2007) R. W. Chhajlany, P. Tomczak, A. Wójcik, and J. Richter, Physical Review A 75, 032340 (2007).
  * Enderle  _et al._ (2010) M. Enderle, B. Fåk, H.-J. Mikeska, R. Kremer, A. Prokofiev, and W. Assmus, Physical review letters 104, 237207 (2010).
  * Schlappa  _et al._ (2018) J. Schlappa, U. Kumar, K. Zhou, S. Singh, M. Mourigal, V. Strocov, A. Revcolevschi, L. Patthey, H. Rønnow, S. Johnston, _et al._ , Nature Communications 9, 5394 (2018).
  * Do  _et al._ (2017a) S.-H. Do, S.-Y. Park, J. Yoshitake, J. Nasu, Y. Motome, Y. S. Kwon, D. Adroja, D. Voneshen, K. Kim, T.-H. Jang, _et al._ , arXiv preprint arXiv:1703.01081 (2017a), 10.48550/arXiv.1703.01081.
  * Koga  _et al._ (2021) A. Koga, Y. Murakami, and J. Nasu, Physical Review B 103, 214421 (2021).
  * Wulferding  _et al._ (2020) D. Wulferding, Y. Choi, S.-H. Do, C. H. Lee, P. Lemmens, C. Faugeras, Y. Gallais, and K.-Y. Choi, Nature communications 11, 1603 (2020).
  * Hashimoto  _et al._ (2022) A. Hashimoto, Y. Murakami, and A. Koga, in  _Journal of Physics: Conference Series_, Vol. 2164 (IOP Publishing, 2022) p. 012028.
  * Kitaev (2006) A. Kitaev, Annals of Physics 321, 2 (2006), january Special Issue.
  * Hermanns  _et al._ (2015) M. Hermanns, K. O’Brien, and S. Trebst, Physical Review Letters 114, 157202 (2015).
  * Le Hur  _et al._ (2017) K. Le Hur, A. Soret, and F. Yang, Physical Review B 96, 205109 (2017).
  * Mandal  _et al._ (2011) S. Mandal, S. Bhattacharjee, K. Sengupta, R. Shankar, and G. Baskaran, Physical Review B 84, 155121 (2011).
  * Baskaran  _et al._ (2007) G. Baskaran, S. Mandal, and R. Shankar, Phys. Rev. Lett. 98, 247201 (2007).
  * Baskaran  _et al._ (2008) G. Baskaran, D. Sen, and R. Shankar, Phys. Rev. B 78, 115116 (2008).
  * Baskaran (2023) G. Baskaran, arXiv preprint arXiv:2309.07119 (2023).
  * Singh  _et al._ (2012) Y. Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W. Ku, S. Trebst, and P. Gegenwart, Physical review letters 108, 127203 (2012).
  * Hwan Chun  _et al._ (2015) S. Hwan Chun, J.-W. Kim, J. Kim, H. Zheng, C. C. Stoumpos, C. Malliakas, J. Mitchell, K. Mehlawat, Y. Singh, Y. Choi, _et al._ , Nature Physics 11, 462 (2015).
  * Trebst and Hickey (2022a) S. Trebst and C. Hickey, Physics Reports 950, 1 (2022a), kitaev materials.
  * Chaloupka  _et al._ (2010) J. c. v. Chaloupka, G. Jackeli, and G. Khaliullin, Phys. Rev. Lett. 105, 027204 (2010).
  * Suzuki  _et al._ (2021) H. Suzuki, H. Liu, J. Bertinshaw, K. Ueda, H. Kim, S. Laha, D. Weber, Z. Yang, L. Wang, H. Takahashi, _et al._ , Nature communications 12, 4512 (2021).
  * Do  _et al._ (2017b) S.-H. Do, S.-Y. Park, J. Yoshitake, J. Nasu, Y. Motome, Y. S. Kwon, D. Adroja, D. Voneshen, K. Kim, T.-H. Jang, _et al._ , Nature Physics 13, 1079 (2017b).
  * Banerjee  _et al._ (2017) A. Banerjee, J. Yan, J. Knolle, C. A. Bridges, M. B. Stone, M. D. Lumsden, D. G. Mandrus, D. A. Tennant, R. Moessner, and S. E. Nagler, Science 356, 1055 (2017).
  * Trebst and Hickey (2022b) S. Trebst and C. Hickey, Physics Reports 950, 1 (2022b).
  * Winter  _et al._ (2016) S. M. Winter, Y. Li, H. O. Jeschke, and R. Valentí, Physical Review B 93, 214431 (2016).
  * Janssen  _et al._ (2016) L. Janssen, E. C. Andrade, and M. Vojta, Phys. Rev. Lett. 117, 277202 (2016).
  * Katukuri  _et al._ (2016) V. M. Katukuri, R. Yadav, L. Hozoi, S. Nishimoto, and J. van den Brink, Scientific reports 6, 29585 (2016).
  * Takayama  _et al._ (2015) T. Takayama, A. Kato, R. Dinnebier, J. Nuss, H. Kono, L. Veiga, G. Fabbris, D. Haskel, and H. Takagi, Physical review letters 114, 077202 (2015).
  * Tsirlin and Gegenwart (2022) A. A. Tsirlin and P. Gegenwart, physica status solidi (b) 259, 2100146 (2022).
  * Gotfryd  _et al._ (2017) D. Gotfryd, J. Rusnačko, K. Wohlfeld, G. Jackeli, J. Chaloupka, and A. M. Oleś, Physical Review B 95, 024426 (2017).
  * Schaffer  _et al._ (2012) R. Schaffer, S. Bhattacharjee, and Y. B. Kim, Physical Review B 86, 224417 (2012).
  * Feng  _et al._ (2007) X.-Y. Feng, G.-M. Zhang, and T. Xiang, Phys. Rev. Lett. 98, 087204 (2007).
  * Holleis  _et al._ (2021) L. Holleis, J. C. Prestigiacomo, Z. Fan, S. Nishimoto, M. Osofsky, G.-W. Chern, J. van den Brink, and B. Shivaram, npj Quantum Materials 6, 66 (2021).
  * Kelker (1973) H. Kelker, Molecular Crystals and Liquid Crystals 21, 1 (1973), https://doi.org/10.1080/15421407308083312 .
  * Kelker (1988) H. Kelker, Molecular Crystals and Liquid Crystals Incorporating Nonlinear Optics 165, 1 (1988), https://doi.org/10.1080/00268948808082195 .
  * Mermin and Wagner (1966) N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).
  * Alexander  _et al._ (2012) G. P. Alexander, B. G.-g. Chen, E. A. Matsumoto, and R. D. Kamien, Rev. Mod. Phys. 84, 497 (2012).
  * Läuchli  _et al._ (2006) A. Läuchli, F. Mila, and K. Penc, Phys. Rev. Lett. 97, 087205 (2006).
  * Tsunetsugu and Arikawa (2006) H. Tsunetsugu and M. Arikawa, Journal of the Physical Society of Japan 75, 083701 (2006), https://doi.org/10.1143/JPSJ.75.083701 .
  * Lacroix  _et al._ (2011) C. Lacroix, P. Mendels, and F. Mila,  _Introduction to frustrated magnetism: materials, experiments, theory_ , Vol. 164 (Springer Science & Business Media, 2011).
  * Blume and Hsieh (1969) M. Blume and Y. Hsieh, Journal of Applied Physics 40, 1249 (1969).
  * Hu  _et al._ (2019) W.-J. Hu, S.-S. Gong, H.-H. Lai, H. Hu, Q. Si, and A. H. Nevidomskyy, Phys. Rev. B 100, 165142 (2019).
  * Podolsky and Demler (2005) D. Podolsky and E. Demler, New Journal of Physics 7, 59 (2005).
  * Tanaka and Hotta (2020) K. Tanaka and C. Hotta, Phys. Rev. B 101, 094422 (2020).
  * Hikihara  _et al._ (2008) T. Hikihara, L. Kecke, T. Momoi, and A. Furusaki, Phys. Rev. B 78, 144404 (2008).
  * Parvej and Kumar (2017) A. Parvej and M. Kumar, Phys. Rev. B 96, 054413 (2017).
  * Sudan  _et al._ (2009) J. Sudan, A. Lüscher, and A. M. Läuchli, Phys. Rev. B 80, 140402 (2009).
  * Nasu  _et al._ (2017) J. Nasu, Y. Kato, J. Yoshitake, Y. Kamiya, and Y. Motome, Phys. Rev. Lett. 118, 137203 (2017).
  * Takahashi  _et al._ (2021) M. O. Takahashi, M. G. Yamada, D. Takikawa, T. Mizushima, and S. Fujimoto, Phys. Rev. Res. 3, 023189 (2021).
  * Penc and Läuchli (2011) K. Penc and A. M. Läuchli, “Spin nematic phases in quantum spin systems,” in  _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_, edited by C. Lacroix, P. Mendels, and F. Mila (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011) pp. 331–362.
  * Grafe  _et al._ (2017) H.-J. Grafe, S. Nishimoto, M. Iakovleva, E. Vavilova, L. Spillecke, A. Alfonsov, M.-I. Sturza, S. Wurmehl, H. Nojiri, H. Rosner, _et al._ , Scientific reports 7, 6720 (2017).
  * Agrapidis  _et al._ (2019) C. E. Agrapidis, J. van den Brink, and S. Nishimoto, Phys. Rev. B 99, 224418 (2019).
  * Catuneanu  _et al._ (2019) A. Catuneanu, E. S. Sørensen, and H.-Y. Kee, Phys. Rev. B 99, 195112 (2019).
  * (53) Supplemental materials on the following topics in relation to the main text: (I) Schematic picture of quadrupolar ordered states, (II) Wave function analysis in two-leg ladder, (III) Quadrupolar order and correlation in 2-leg ladder from ED, (IV) Finite-size effects in the quadrupolar correlation functions in 2-leg ladder.
  * Shannon  _et al._ (2006) N. Shannon, T. Momoi, and P. Sindzingre, Phys. Rev. Lett. 96, 027213 (2006).
  * Yadav  _et al._ (2016) R. Yadav, N. A. Bogdanov, V. M. Katukuri, S. Nishimoto, J. Van Den Brink, and L. Hozoi, Scientific reports 6, 37925 (2016).