1. Introduction
Strongly correlated electron systems with multiple internal degrees of freedom provide a fertile ground for realizing exotic quantum phases beyond the paradigms of conventional electronic orders [1–3]. In addition to the charge and spin degrees of freedom, orbital, layer, and valley indices often play an essential role in shaping the low-energy physics of a wide range of materials, including transition-metal compounds [4–6], moiré superlattices [7–10], and ultracold atomic systems [11, 12]. In such multi-component settings, these internal degrees of freedom are not merely passive labels but actively participate in collective phenomena, giving rise to intertwined electronic orders and novel forms of quantum entanglement.
A paradigmatic framework for exploring such physics is provided by Kugel–Khomskii type models [13, 14], where spin and orbital degrees of freedom are coupled through exchange interactions generated by virtual charge fluctuations. Traditionally, most studies of these models focus on phases characterized by decoupled spin and orbital orders, such as antiferromagnetic or ferromagnetic order accompanied by ferro- or antiferro-orbital order [15–17]. Even in systems with enlarged symmetries, such as SU(4)-symmetric models [18], the emphasis has largely been placed on symmetry enhancement, and exotic quantum liquid behavior [19–22]. By contrast, the structure and consequences of local spin-orbital entanglement itself have received comparatively limited attention [23, 24].
This gap in understanding invites the exploration of systems with additional, controllable internal degrees of freedom. Bilayer and multilayer materials offer a natural extension in this regard: the layer index acts as a synthetic quantum degree of freedom, which can compete or cooperate with spin, orbital, and charge dynamics. This perspective allows one to treat bilayer systems within a unified framework of multi-component correlated electrons, where the interplay between multiple internal degrees of freedom can be systematically studied and potentially engineered.
Propelled by this perspective, bilayer and multilayer quantum materials have recently emerged as a central theme in condensed matter physics, driven by advances ranging from twisted moiré superlattices [25–28] to the newly discovered high-Tc superconductors in bilayer nickelates [29–31]. These platforms are particularly intriguing because the layer index introduces a synthetic and tunable degree of freedom that can compete or cooperate with spin and charge dynamics. While much of the current research has been focused on band engineering and Fermiology [32–39], the interplay between the interlayer geometry and the orbital character of electrons offers a distinct and less explored route toward unconventional correlated states. A compelling example is found in magic-angle twisted bilayer graphene (MATBG) near quarter filling, which manifests as a ferromagnetic anomalous Hall effect insulator [40, 41]. Intriguingly, experimental evidence suggests that the charge gap in MATBG can develop prior to the onset of long-range magnetic order [42, 43], pointing toward a robust Hubbard–Mott mechanism where internal degrees of freedom—such as valley and layer pseudospins—play a decisive role [44, 45].
Another realization of correlated bilayer system is provided by the recently discovered bilayer nickelate superconductors, such as La3Ni2O7. These materials consist of NiO bilayers separated by insulating spacer layers and exhibit high-temperature superconductivity under pressure. Owing to the multi-orbital nature of the Ni 3d manifold and the pronounced structural anisotropy, bilayer nickelates naturally host strong electronic correlations together with substantial interlayer coupling effects. First-principles and spectroscopic measurements have highlighted the crucial role of orbital-dependent hybridization, with the ${d}_{{{z}}^{2}}$ orbital experiencing significant interlayer bonding–antibonding splitting, while the hopping between the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals remains largely within each layer [46]. This makes bilayer nickelates a promising platform for exploring unconventional superconductivity and spin–layer–orbital entanglement phenomena beyond the paradigms established in cuprates [47–85].
More generally, in transition-metal systems with active eg orbitals, the bilayer geometry naturally leads to a phenomenon of strongly orbital-selective interlayer hybridization. Due to the directional nature of d-orbitals, electronic coupling along the stacking direction is highly anisotropic: the orbitals extending vertically (e.g. the ${d}_{{{z}}^{2}}$ orbitals) feel a robust interlayer overlap, whereas the planar orbitals (e.g. the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals) between upper and lower layers remain effectively decoupled [39, 46, 47]. This intrinsic energy hierarchy not only reshapes the band structure, but also acts as an orbital filter that can dynamically quench specific orbital sectors at low energies. As a result, the system enters a unique regime where the surviving layer degree of freedom intertwines with spin, setting the stage for exotic composite entanglement. Beyond bilayer geometries, the concept of orbital-selective physics is also pertinent to other multi-orbital systems, such as the ruthenates (e.g. Sr2RuO4) [86, 87]. Although the specific orbital content differs (involving t2g rather than eg sectors), the essential physics governing the interplay between orbital differentiation and strong correlations shares significant conceptual similarities [88, 89].
In this work, we capitalize on this geometric mechanism to derive a low-energy theory for a quarter-hole-filled two-orbital bilayer Hubbard model. By explicitly treating the strong interlayer hybridization of ${d}_{{{z}}^{2}}$ orbitals within a molecular-orbital basis, we project out the high-energy spin singlet excitations and construct an effective anisotropic Kugel–Khomskii Hamiltonian describing the coupled dynamics of electron spin and emergent layer pseudospin degrees of freedom associated with the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals in the upper and lower layers. Through a combination of Weiss mean-field theory and generalized flavor-wave theory, we uncover a rich phase diagram. Most notably, we identify a novel spin-layer-entangled (SLE) phase in the strong spin-layer coupling regime. In this state, conventional long-range magnetic and layer orders are simultaneously melted, giving way to a hidden composite order characterized by maximal local entanglement between spin and layer sectors. We further show that this phase arises from the spontaneous breaking of an emergent O(4) symmetry, manifesting in a unique excitation spectrum with three gapless Goldstone modes.
2. Model and Hamiltonian
2.1. Bilayer two-orbital Hubbard model
We consider a bilayer two-orbital Hubbard model that captures the essential low-energy physics of transition-metal bilayer systems with strong orbital-selective interlayer hybridization. Here we take the two-orbital bilayer Hubbard model for the nickelate superconductor La3Ni2O7 as a representative example [57, 58].
In La3Ni2O7, the t2g orbitals are fully occupied and each layer hosts two eg orbitals, ${d}_{{{x}}^{2}-{{y}}^{2}}$ and ${d}_{{{z}}^{2}}$, with the average electron filling corresponding to three electrons per bilayer rung, e.g. a quarter-hole-filled configuration [29, 46]. The kinetic properties of these eg electrons are strictly governed by the highly spatially anisotropic nature of the atomic orbitals: the wave function of ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals is extended within the crystal plane, with lobes pointing towards the in-plane oxygen ligands. This geometry facilitates strong intralayer hybridization. Meanwhile, the interlayer overlap between ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals is negligible due to the lack of vertical extension of the wave function. In stark contrast, the ${d}_{{{z}}^{2}}$ orbitals feature lobe-shaped electron densities elongated along the c-axis, enabling a distinct interlayer hopping channel that is characteristic of ${d}_{{{z}}^{2}}$ symmetry [39, 90, 91].
Based on the orbital-selective physics, we start from the bilayer within the two eg orbital sector. The Hamiltonian reads
$\begin{eqnarray}H={H}_{{\rm{TB}}}+{H}_{{\rm{int}}}.\end{eqnarray}$
Here, HTB is a minimal tight-binding Hamiltonian (as sketched in figure 1(a)): $\begin{eqnarray}{H}_{{\rm{TB}}}=\displaystyle \sum _{ij\sigma \eta }{t}_{\parallel }^{{xx}}{d}_{i{x}\sigma }^{\eta \dagger }{d}_{j{x}\sigma }^{\eta }+\displaystyle \sum _{i\sigma }{t}_{\perp }^{{zz}}{d}_{i{z}\sigma }^{{\rm{T}}\dagger }{d}_{i{z}\sigma }^{{\rm{B}}},\end{eqnarray}$
where ${d}_{i\alpha \sigma }^{\eta \dagger }\,({d}_{i\alpha \sigma }^{\eta })$ creates (annihilates) an electron in orbital α (α = x, z denotes the two eg orbitals, ${d}_{{{x}}^{2}-{{y}}^{2}}$ and ${d}_{{{z}}^{2}}$, respectively) with spin σ at site i of layer η (η = T, B denotes the top and bottom layer, respectively). The hopping amplitude ${t}_{\parallel }^{{xx}}$ is the nearest-neighbor intralayer hopping of the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals, while ${t}_{\perp }^{{zz}}$ captures the strong interlayer hybridization between vertically aligned ${d}_{{{z}}^{2}}$ orbitals. Direct interlayer hopping of the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals is neglected due to their planar orbital character. The presence of strong interlayer hopping between the ${d}_{{{z}}^{2}}$ orbitals t⊥ has a significant effect on the single-site spectrum, and will be discussed in the next subsection.
Figure 1. (a) Schematic illustration of the minimal bilayer two-orbital tight-binding model where ${t}_{\perp }^{{zz}}$ and ${t}_{\parallel }^{{xx}}$ denote the orbital dependent hopping parameters. (b) Sketch of the crystal splitting of eg orbitals and the formation of the bonding–antibonding molecular orbital (MO) states between Ni z2 orbitals in the top and bottom layers. (c) Four degenerate ground-state configurations ∣Sz, τz⟩ labeled with the spin and layer quantum numbers ${S}^{{\rm{z}}}=\pm \frac{1}{2}$ and ${\tau }^{{\rm{z}}}=\pm \frac{1}{2}$ in the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbital subspace. |
The on-site term Hint consists of two contributions, the crystal field splitting and the local electron–electron interactions, which are taken in the Kanamori form [92]:
$\begin{eqnarray}\begin{array}{rcl}{H}_{{\rm{int}}} & = & \displaystyle \sum _{i\alpha \eta }{\epsilon }_{\alpha }{n}_{i\alpha }^{\eta }+U\displaystyle \sum _{i,\alpha ,\eta }{n}_{i\alpha \uparrow }^{\eta }{n}_{i\alpha \downarrow }^{\eta }\\ & & +\displaystyle \sum _{i,\alpha \lt \beta ,\sigma ,\eta }\{{U}^{{\prime} }{n}_{i\alpha \sigma }^{\eta }{n}_{i\beta \bar{\sigma }}^{\eta }+({U}^{{\prime} }-{J}_{{\rm{H}}}){n}_{i\alpha \sigma }^{\eta }{n}_{i\beta \sigma }^{\eta }\\ & & -{J}_{{\rm{H}}}({d}_{i\alpha \sigma }^{\eta \dagger }{d}_{i\alpha \bar{\sigma }}^{\eta }{d}_{i\beta \bar{\sigma }}^{\eta \dagger }{d}_{i\beta \sigma }^{\eta }+{d}_{i\alpha \sigma }^{\eta \dagger }{d}_{i\alpha \bar{\sigma }}^{\eta \dagger }{d}_{i\beta \sigma }^{\eta }{d}_{i\beta \bar{\sigma }}^{\eta })\},\end{array}\end{eqnarray}$
where ${n}_{i\alpha \sigma }^{\eta }={d}_{i\alpha \sigma }^{\eta \dagger }{d}_{i\alpha \sigma }^{\eta }$. Here U, ${U}^{{\prime} }$ and JH, represent the intra- and inter-orbital Coulomb repulsion and the Hund's coupling, respectively, satisfying ${U}^{{\prime} }=U-2{J}_{{\rm{H}}}$. We assume a crystal field splitting εx > εz, consistent with first-principles results [46] for bilayer nickelates (figure 1 (b)).2.2. Molecular orbital basis in the strong coupling limit
A unique feature of this La3Ni2O7 bilayer system is the significant orbital-selective interlayer hybridization along with fractional electron occupation per Ni ion. While the interlayer coupling between ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals is negligible, the strong interlayer hopping of the ${d}_{{{z}}^{2}}$ orbitals t⊥ leads to the formation of bonding–antibonding molecular orbitals (MOs), illustrated in figure 1(b). To capture this effect, we work with the MO basis of ${d}_{{{z}}^{2}}$ orbitals
$\begin{eqnarray}{d}_{i{\rm{z}}\sigma }^{{\rm{b}}({\rm{a}})}=\frac{1}{\sqrt{2}}({d}_{i{\rm{z}}\sigma }^{{\rm{T}}}\pm {d}_{i+{\delta }_{z}{\rm{z}}\sigma }^{{\rm{B}}}),\end{eqnarray}$
where the index b(a) corresponds to the bonding (antibonding) MO. In the MO basis, the on-site energy of ${d}_{{{z}}^{2}}$ orbital is renormalized to ${\epsilon }_{{z}}\mp {t}_{\perp }^{{zz}}$. This energy separation is the crucial mechanism that stabilizes the low-energy manifold.The Hamiltonians of equations (2 ) and (3 ) are recast in the MO basis as $H={H}_{{\rm{TB}}}^{{\rm{MO}}}+{H}_{{\rm{int}}}^{{\rm{MO}}}$, where the interaction terms are transformed accordingly. We focus on the quarter-hole filling regime, corresponding to a total occupancy of n = 3 electrons per rung. Theoretical calculations consistently suggest that the ground-state configuration relevant to La3Ni2O7 is a low-spin state [29, 38, 39, 46] where the bonding ${d}_{{{z}}^{2}}$ orbital is doubly occupied by a spin-singlet and the degenerate ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals are quarterly filled, as illustrated in figure 1(b). Note that this ground-state configuration is four-fold degenerate, reflecting the fact that the remaining electron can occupy the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbital on either the top or bottom layer with either spin orientation. To capture the layer degree of freedom of the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals, we introduce a layer pseudospin operator τ, where ${\tau }^{{z}}=\pm \frac{1}{2}$ characterizes the occupation on the top and the bottom layers, respectively. The four-fold low-spin configurations are therefore labelled by the ∣Sz, τz⟩, where ${S}^{{z}}=\pm \frac{1}{2}$ and ${\tau }^{{z}}=\pm \frac{1}{2}$ denote the total electron spin and layer pseudospin, respectively. The schematic illustration is shown in figure 1(c).
2.3. Schrieffer–Wolff transformation and second perturbation
Given the strong coupling limit where the in-plane hopping amplitude is significantly smaller than the interaction scale (e.g. ${t}_{\parallel }^{{xx}}\ll U$), we treat the in-plane kinetic term ${H}^{{\prime} }$ as a perturbation to the local Hamiltonian H0. In this regime, the ground state manifold is determined primarily by local interactions, and the effects of itinerancy enter only through higher order virtual processes, justifying a controlled strong-coupling approach. To capture the resulting low-energy physics within the degenerate ground-state manifold, we derive an effective Hamiltonian via a Schrieffer–Wolff transformation, Heff = eSHe−S, which eliminates high-energy charge excitations and projects the dynamics onto the low-energy subspace ${ \mathcal S }$ [93].
The leading contribution to the effective interaction arises from second-order virtual hopping processes [94], ${H}_{\,\rm{eff}\,}\approx \frac{1}{2}[S,{H}^{{\prime} }]$. Starting from a product of local ground-state configurations ∣Sz, τz⟩i ⨂ ∣Sz, τz⟩j, the hopping term ${t}_{\parallel }^{{xx}}$ transfers an electron in the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbital to a neighboring site, generating a high-energy intermediate state. The large energy cost of this intermediate state originates not only from the Coulomb repulsion U but also from the suppression of Hund's coupling: specifically, the rigid spin singlet formed by the electrons in the bonding ${d}_{{{z}}^{2}}$ orbital prevents the itinerant ${d}_{{{x}}^{2}-{{y}}^{2}}$ electron from aligning ferromagnetically with the ${d}_{{{z}}^{2}}$ electrons, thereby suppressing the Hund's energy gain.
We classify these virtual fluctuations into two channels based on the layer configuration as shown in figure 2. In Type-I processes (${\tau }_{i}^{{z}}={\tau }_{j}^{{z}}$), electrons reside in the same layer. The virtual hopping creates a double occupancy that strictly conserves the layer index τz while mediating spin exchange. In Type-II processes (${\tau }_{i}^{{z}}\ne {\tau }_{j}^{{z}}$), electrons reside in opposite layers. Here, the return hop allows for mixing between different basis states, mediating coupled fluctuations of both spin and layer pseudospin.
Figure 2. Schematic illustration of the major virtual hopping processes in the second-order perturbation expansion. For clarity, the bonding ${d}_{{{z}}^{2}}$ orbitals are omitted as they remain fully occupied throughout the process. (a) Type-I: Processes involving electrons in the same layer. (b) Type-II: Processes involving electrons in different layers. |
Integrating these processes yields an effective anisotropic Kugel–Khomskii type Hamiltonian:
$\begin{eqnarray}\begin{array}{rcl}{H}_{\,\rm{eff}\,} & = & \displaystyle \sum _{ij}{J}_{{\rm{s}}}{\hat{{\boldsymbol{S}}}}_{{\boldsymbol{i}}}\cdot {\hat{{\boldsymbol{S}}}}_{{\boldsymbol{j}}}+\displaystyle \sum _{\alpha }{K}_{{\rm{s}}}^{\alpha \alpha }(\frac{1}{4}-{\hat{{\boldsymbol{S}}}}_{{\boldsymbol{i}}}\cdot {\hat{{\boldsymbol{S}}}}_{{\boldsymbol{j}}}){\hat{\tau }}_{i}^{\alpha }{\hat{\tau }}_{j}^{\alpha }\\ & & +{K}_{{\rm{t}}}^{\alpha \alpha }({\hat{{\boldsymbol{S}}}}_{{\boldsymbol{i}}}\cdot {\hat{{\boldsymbol{S}}}}_{{\boldsymbol{j}}}+\frac{3}{4}){\hat{\tau }}_{i}^{\alpha }{\hat{\tau }}_{j}^{\alpha }.\end{array}\end{eqnarray}$
Due to the complex multiplet structure of the intermediate states, analytical expressions for the exchange constants are unwieldy. Instead, we can determine the relationships between effective parameters $\{{K}_{{\rm{s}}}^{\alpha \alpha },{K}_{{\rm{t}}}^{\alpha \alpha },{J}_{{\rm{s}}}\}$ and U, JH, ${t}_{\perp }^{{zz}}$, εα by numerically diagonalizing the associated two-site Hubbard Hamiltonian.It is important to note that τz is not a strictly conserved quantum number when inter-orbital fluctuations are taken into account. As a result, the system does not preserve a U(1) rotational symmetry in the layer sector. Consequently, the transverse exchange couplings are anisotropic, i.e. ${K}_{{\rm{s}}/{\rm{t}}}^{{xx}}\ne {K}_{{\rm{s}}/{\rm{t}}}^{{yy}}$. Nevertheless, since the τz violations are suppressed by the small amplitudes of the mixed layer components, this anisotropy is quantitatively negligible. That is, the system retains an approximate U(1) symmetry with Kxx ≈ Kyy. In summary, the effective model, in the Kxx = Kyy limit, has an O(3)spin × O(2)layer symmetry, where the layer O(2)layer symmetry comes from the combination of two distinct symmetries: the continuous U(1) symmetry in the transverse channel ${ \mathcal G }(\theta )\equiv \exp [{\rm{i}}{\sum }_{i}{\hat{\tau }}_{i}^{z}\theta ]$, and a discrete ${{\mathbb{Z}}}_{2}$ symmetry associated with the exchange between top and bottom layers, i.e. ${ \mathcal I }\equiv {\prod }_{i}{\hat{\tau }}_{i}^{{x}}$. The spontaneous breaking of these symmetries gives rise to different ordered phases, as will be discussed in the next section.
3. Ground-state phase diagram of the effective Kugel–Khomskii model
3.1. Weiss mean-field theory
In this section, we investigate the ground-state phase diagram of the effective anisotropic Kugel–Khomskii model of equation (5 ) within the Weiss mean-field theory. This formalism is equivalent to adopting a site factorized trial wavefunction
$\begin{eqnarray}| {\rm{\Psi }}\rangle ={\otimes }_{i}| {\boldsymbol{d}}{\rangle }_{i}.\end{eqnarray}$
Here ∣d⟩i is the local wave function at site i and can be expressed as a coherent superposition of the basis states ∣Sz, τz⟩i, $\begin{eqnarray}\begin{array}{rc}| {\boldsymbol{d}}{\rangle }_{i} & =\displaystyle \sum _{{S}^{{z}},{\tau }^{{z}}=\pm \frac{1}{2}}{d}_{i,{S}^{{z}},{\tau }^{{z}}}| {S}^{{z}},{\tau }^{{z}}{\rangle }_{i},\end{array}\end{eqnarray}$
where ${d}_{i,{S}^{{z}},{\tau }^{{z}}}$ are complex variational coefficients. The ground state can then be determined by variationally minimizing the energy E = ⟨$\Psi$∣H∣$\Psi$⟩/⟨$\Psi$∣$\Psi$⟩. While this trial wave function does not capture the quantum entanglement between different sites, this treatment retains the on-site entanglement between spin and layer pseudospin degrees of freedom. As a result, it is capable of describing exotic SLE states, as will be discussed later.To explore the generic features of the model without being restricted to a specific material La3Ni2O7, we treat the interaction strengths in equation (5 ) as independent parameters. Owing to the high dimensionality of the parameter space, our investigation focuses on a representative cross-section that captures the essential physics.
We fix the exchange scales in the pure spin and spin-layer triplet channels to Js = 0.2, ${K}_{{\rm{t}}}^{{zz}}=1.5$, and ${K}_{{\rm{t}}}^{{xx}}={K}_{{\rm{t}}}^{{yy}}=1.2$, and then systematically investigate the influence of the coupling ${K}_{{\rm{s}}}^{\alpha \alpha }$ on the system's ground state. For simplicity, we assume ${K}_{{\rm{s}}}^{{xx}}={K}_{{\rm{s}}}^{{yy}}={{\rm{\Delta }}}_{{\rm{s}}}{K}_{{\rm{s}}}^{{zz}}$. The phase diagram is shown in figure 3, which exhibits four distinct phases: (1) spin-ferromagnetic and layer staggered (FM-LS), (2) spin-antiferromagnetic and layer staggered (AFM-LS), (3) spin-antiferromagnetic and interlayer coherent (AFM-LC), and (4) SLE phase. Configurations of the first three phases are illustrated in panels (b)–(d) of figure 3.
Figure 3. (a) Ground-state phase diagram of the effective Kugel–Khomskii model in equation ( |
3.2. Layer staggered phases
We first consider the regime dominated by the Ising pseudospin anisotropy where ∣Δs∣ is small. Throughout this regime, the system spontaneously breaks the pseudospin ${{\mathbb{Z}}}_{2}$ symmetry $\langle {\hat{\tau }}^{{z}}\rangle \ne 0$, forming layer-staggered electron occupation where electrons alternately populate the top and bottom layers on adjacent sites, i.e. $\langle {\hat{\tau }}_{i}^{{z}}{\hat{\tau }}_{j}^{{z}}\rangle =-\frac{1}{4}$. This pattern manifests as a checkerboard charge order in the two layers. In addition, the spin sectors can exhibit either ferromagnetic or antiferromagnetic orders, depending on the sign of the coupling strength ${K}_{{\rm{s}}}^{{zz}}$, as shown in figures 3(b) and (c).
Within the layer-staggered phases, the ground state energy can be written as:
$\begin{eqnarray}{E}_{\,\rm{LS}\,}=\frac{4{J}_{{\rm{s}}}-{K}_{{\rm{t}}}^{{zz}}+{K}_{{\rm{s}}}^{{zz}}}{4}\langle {\hat{{\boldsymbol{S}}}}_{{\boldsymbol{i}}}\cdot {\hat{{\boldsymbol{S}}}}_{{\boldsymbol{j}}}\rangle -\left(\frac{{K}_{{\rm{s}}}^{{zz}}}{16}+\frac{3{K}_{{\rm{t}}}^{{zz}}}{16}\right).\end{eqnarray}$
A transition from the FM to the AFM spin configuration is driven by the competition between the effective exchange interactions ${K}_{{\rm{s}}}^{{zz}}$ and ${K}_{{\rm{t}}}^{{zz}}$. In particular, the AFM order becomes energetically favorable when ${K}_{{\rm{s}}}^{{zz}}+4{J}_{{\rm{s}}}\gt {K}_{{\rm{t}}}^{{zz}}$.3.3. Inter-layer coherent phase
In another limit where Δs is large and the pseudospin exhibits XY anisotropy, the system tends to break the pseudospin U(1) symmetry and establish off-diagonal correlations in the pseudospin sector: $\langle {\hat{\tau }}^{\pm }\rangle \ne 0$. This signals the emergence of spontaneous inter-layer coherence, in which electrons form coherent quantum superpositions of occupying the top and bottom layers, see figure 3(d). In the mean time, spins are ordered antiferromagnetically to synergistically minimize the total energy. The ground state energy per bond is:
$\begin{eqnarray}{E}_{\,\rm{c}\,}=-\frac{{J}_{{\rm{s}}}}{4}-\frac{1}{8}({K}_{{\rm{t}}}^{{xx}}+{K}_{{\rm{s}}}^{{xx}}).\end{eqnarray}$
Physically, the vanishing of $\langle {\hat{\tau }}^{{z}}\rangle $ implies the melting of the spatially staggered layer occupation. The emergence of in-plane layer pseudospin order corresponds to a state with spontaneous inter-layer coherence, where electrons form a quantum superposition between the top and bottom layers with a specific relative phase. In contrast to the layer staggered phase, this coherent state preserves layer symmetry and therefore does not induce any electric polarization or static charge order in real space. Moreover, since the inter-layer coherent order resides entirely in the off-diagonal pseudospin channel, it does not gap out the charge sector. As a result, no charge excitation gap associated with layer polarization is generated in this phase, as shown in the next section.3.4. Spin-layer-entangled
Interestingly, upon decreasing the anisotropy Δs, the system enters an intermediate regime where the spin and layer degrees of freedom become strongly entangled. This phase breaks both the inversion ${ \mathcal I }$ and time reversal symmetries. However, such symmetry breaking pattern cannot be characterized by the spins or pseudospins order alone: In this phase, $\langle {\hat{{\boldsymbol{S}}}}_{i}\rangle =0$ and $\langle {\hat{{\boldsymbol{\tau }}}}_{i}\rangle =0$. Instead, it is necessary to introduce a composite spin–layer order parameter ${\hat{{ \mathcal Q }}}_{i}^{\beta \alpha }\equiv {\hat{S}}_{i}^{\beta }{\hat{\tau }}_{i}^{\alpha }$, with α, β ∈ {x, y, z}, to describe the symmetry breaking of this phase. The composite spin–layer order can be well understood by rewriting the Hamiltonian equation (5 ) in terms of these operators:
$\begin{eqnarray*}\begin{array}{rcl}{H}_{\,\rm{eff}\,} & = & \displaystyle \sum _{ij}{J}_{{\rm{s}}}\hat{{{\boldsymbol{S}}}_{i}}\cdot \hat{{{\boldsymbol{S}}}_{j}}+\displaystyle \sum _{\alpha }\left(\frac{{K}_{{\rm{s}}}^{\alpha \alpha }}{4}+\frac{3{K}_{{\rm{t}}}^{\alpha \alpha }}{4}\right){\hat{\tau }}_{i}^{\alpha }{\hat{\tau }}_{j}^{\alpha }\\ & & +\displaystyle \sum _{\alpha \beta }({K}_{{\rm{t}}}^{\alpha \alpha }-{K}_{{\rm{s}}}^{\alpha \alpha }){\hat{{ \mathcal Q }}}_{i}^{\beta \alpha }{\hat{{ \mathcal Q }}}_{j}^{\beta \alpha }.\end{array}\end{eqnarray*}$
In the SLE regime of the phase diagram, the ${\hat{{ \mathcal Q }}}_{i}{\hat{{ \mathcal Q }}}_{j}$ term dominates. As a result, a SLE order with $\langle {\hat{{ \mathcal Q }}}_{i}^{\beta \alpha }\rangle \ne 0$ is established in the ground state. Meanwhile, the spin and pseudospin are disordered, $\langle {\hat{{\boldsymbol{S}}}}_{i}\rangle =0$ and $\langle {\hat{{\boldsymbol{\tau }}}}_{i}\rangle =0$, as depicted in figures B.1(a), (b). Despite this, all nine components of ${{ \mathcal Q }}_{i}^{\beta \alpha }$ are generally non-zero. The relation that $\langle {\hat{{ \mathcal Q }}}_{i}^{\beta \alpha }\rangle =\langle {\hat{S}}_{i}^{\beta }{\hat{\tau }}_{i}^{\alpha }\rangle \ne \langle {\hat{S}}_{i}^{\beta }\rangle \langle {\hat{\tau }}_{i}^{\alpha }\rangle =0$ signals the emergence of a nontrivial composite ordered phase with a maximal spin-layer entanglement. In fact, we find that the on-site ground-state wavefunction obtained by the Weiss mean-field theory satisfies the following form $\begin{eqnarray}\begin{array}{rcl}| {\boldsymbol{d}}{\rangle }_{i} & = & a| +\frac{1}{2},+\frac{1}{2}{\rangle }_{i}+{a}^{* }| -\frac{1}{2},-\frac{1}{2}{\rangle }_{i}\\ & & +{\rm{i}}b| +\frac{1}{2},-\frac{1}{2}{\rangle }_{i}+{\rm{i}}{b}^{* }| -\frac{1}{2},+\frac{1}{2}{\rangle }_{i}\end{array}\end{eqnarray}$
for the A sublattice and $\begin{eqnarray}\begin{array}{rcl}| {\boldsymbol{d}}{\rangle }_{i} & = & {\rm{i}}a| +\frac{1}{2},+\frac{1}{2}{\rangle }_{i}-{\rm{i}}{a}^{* }| -\frac{1}{2},-\frac{1}{2}{\rangle }_{i}\\ & & +b| +\frac{1}{2},-\frac{1}{2}{\rangle }_{i}-{\rm{i}}{({\rm{i}}b)}^{* }| -\frac{1}{2},+\frac{1}{2}{\rangle }_{i}\end{array}\end{eqnarray}$
for the B sublattice, where a and b are complex numbers satisfying $| a{| }^{2}+| b{| }^{2}=\frac{1}{2}$. This indicates that one can alternatively define order parameter for the SLE phase as a four-component real scalar φ ≡ (ℜa, ℑa, ℜb, ℑb)T. Hence the order parameter φ lives in emergent ${{S}}^{3}/{{\mathbb{Z}}}_{2}$ manifold, as they can be arbitrary chosen on a four-sphere with radius $\sqrt{1/2}$, and that φ and −φ corresponds to the same state. Also, the mean-field ground state possesses an emergent O(4) symmetry, although the Hamiltonian only has a smaller O(3)spin × O(2)layer symmetry.From the mean-field wave function, we observe that the composite order exhibits a characteristic spatial texture depending on the layer pseudospin index α. In particular, the transverse components display an antiferromagnetic pattern (figure B.1(c)), satisfying ${{ \mathcal Q }}_{i}^{\beta \alpha }=-{{ \mathcal Q }}_{i+\delta }^{\beta \alpha }$ for α ∈ {x, y}, whereas the longitudinal component shows a ferromagnetic distribution, ${{ \mathcal Q }}_{i}^{\beta {\rm{z}}}={{ \mathcal Q }}_{i+\delta }^{\beta {\rm{z}}}$, as illustrated in figure B.1(d).
Despite strong fluctuations in the individual spin and pseudospin degreees of freedom, we find that their relative orientation is rigidly locked. This property can be made clear by introducing a composite pseudospin operator
$\begin{eqnarray}{{\boldsymbol{ \mathcal J }}}_{i}={\hat{{\boldsymbol{S}}}}_{i}-4{{ \mathcal Q }}_{i}{\hat{{\boldsymbol{\tau }}}}_{i}.\end{eqnarray}$
One can verify that ${{ \mathcal J }}_{i}$ satisfies the angular momentum algebra $[{{ \mathcal J }}_{i}^{\alpha },{{ \mathcal J }}_{i}^{\beta }]={\rm{i}}{\epsilon }^{\alpha \beta \gamma }{{ \mathcal J }}_{i}^{\gamma }$. Further, it is straightforward to check that ${{\boldsymbol{ \mathcal J }}}_{i}^{2}| {\boldsymbol{d}}{\rangle }_{i}=0$, which indicates that the fluctuating spin Si and the transformed pseudospin degrees of freedom $-4{{ \mathcal Q }}_{i}{\hat{{\boldsymbol{\tau }}}}_{i}$ are locked antiferromagnetically. The system thus suppresses relative spin–layer fluctuations to minimize the strong coupling energy, giving rise to a hidden composite order that is not detectable by conventional probes of dipolar magnetism.4. Excitation spectrum of the effective Kugel–Khomskii model
To better understand the nature of the different phases in the phase diagram, here we discuss the excitation spectrum within each phase presented in the ground-state phase diagram of figure 3(a). We employ a generalized spin-wave theory to describe quasiparticle excitations of this model [95]. Compared to ordinary spin-wave theories where spins and pseudospins are independently treated via the Holstein–Primakoff transformation in each sector, the generalized spin-wave theory here correctly captures the on-site spin-pseudospin entanglement. The details of the formalism are presented in appendix B . The calculated excitation dispersions for the four phases are presented in figure 4.
Figure 4. The excitation dispersion along the high symmetry lines in different phases. (a) FM-LS phase with ${K}_{{\rm{s}}}^{{zz}}=0.5$, Δs = 0.3; (b) AFM-LS phase with ${K}_{{\rm{s}}}^{{zz}}=1.0$, Δs = 0.3; (c) AFM-LC phase with ${K}_{{\rm{s}}}^{{zz}}=1.0$, Δs = 1.5; (d) SLE phase with ${K}_{{\rm{s}}}^{{zz}}=1.0$, Δs = −2. Here, the red and green lines denote the spin and layer excitations, respectively. The black lines relate to the composite spin-layer excitations. |
For both the FM-LS and AFM-LS phases, the excitation spectra exhibit clear separation of energy scales, indicating different origins of excitations. The low-energy gapless excitations (red line in figures 4(a), (b)) are assigned as spin excitations that arise from the spontaneous breaking of the continuous O(3) spin-rotational symmetry. Moreover, the spin nature of low-energy excitations is also manifested in their dispersions: in the FM-LS phase where the spins develop ferromagnetic order, the Goldstone mode has a quadratic dispersion, ω ∝ k2, as shown in figure 4(a); in contrast, in the AFM-LS phase, the antiferromagnetic spin order gives rise to two degenerate linearly dispersing Goldstone modes, ω ∝ k, as shown in figure 4(b). Meanwhile, the high-energy excitations are fully gapped (green lines in figures 4(a), (b)) and well separated from the low-energy spin excitations. These modes are predominantly associated with the layer pseudospin degrees of freedom, and the presence of a finite gap is consistent with the breaking of the discrete pseudospin ${{\mathbb{Z}}}_{2}$ symmetry in the layer-staggered phases.
For the AFM-LC phase, continuous symmetries in both the spin and pseudospin sectors are spontaneously broken. Specifically, in the spin sector the antiferromagnetic order ⟨S⟩ ≠ 0 breaks the O(3) spin-rotational symmetry down to O(2), leading to two branches of linearly dispersive Goldstone modes as gapless spin fluctuations. Meanwhile, the in-plane $\langle {\hat{\tau }}^{+}\rangle \ne 0$ spontaneously breaks the U(1) layer-phase symmetry and gives rise to a third linear Goldstone mode as pseudospin fluctuations, as shown in figure 4(c). The two Goldstone modes have different velocities in general, as they are associated with two independent symmetries.
In addition to these well-separated branches, we identify a nearly flat band, marked by the black lines in figures 4(a)–(c). This mode cannot be classified as a purely spin or purely layer excitation. Instead, it corresponds to the composite spin and layer entangled excitation. Remarkably, as system parameters are tuned, this nearly flat mode gradually develops dispersion and softens. When the gap of this hybrid excitation closes, the system enters the SLE phase, indicating that this mode plays a central role in driving the transition.
The resulting excitation spectrum in the SLE phase (figure 4(d)) highlights the distinctive character of the SLE ground state. In the FM-LS, AFM-LS, and the AFM-LC phases, the low-energy excitations can be unambiguously classified according to their quantum numbers. Notably, even in the AFM-LC phase, where both spin and pseudospin sectors are gapless, the corresponding modes remain decoupled: each Goldstone mode originates from the spontaneous breaking of either spin-rotation symmetry or the layer-phase symmetry, and can therefore be identified as a purely spin or purely pseudospin excitation.
In stark contrast, the low-energy excitations in the SLE phase are intrinsically hybridized. Owing to the spin–orbital locking in the ground state, independent spin and pseudospin fluctuations are suppressed, and the relevant low-energy degrees of freedom correspond to collective rotations of the entangled composite order parameter. As a result, the gapless Goldstone modes cannot be uniquely classified as purely spin or purely pseudospin excitations. Instead, they are associated with rigid rotations of ${{ \mathcal Q }}^{\alpha \beta }$ within the emergent ${{S}}^{3}/{{\mathbb{Z}}}_{2}$ manifold, reflecting the locked nature of spin and layer dynamics. The spontaneous symmetry breaking follows the pattern O(4) → O(3), where the residual O(3) corresponds to the simultaneous rotation of locked spin and pseudospin moments. This symmetry breaking naturally accounts for the three gapless branches observed in the spectrum, as shown in figure 4(d). Note that the spin excitations (red lines) are twofold degenerate.
Although the quasiparticle dispersions obtained above fully characterize the elementary excitations of the effective model, their experimental visibility depends on how they couple to specific probes. Inelastic neutron scattering couples directly to the magnetic moment and therefore predominantly detects fluctuations in the spin sector, making the gapless spin-wave modes observable in this channel. In contrast, the pseudospin modes correspond mainly to interlayer charge redistribution and carry little magnetic spectral weight, so they are essentially invisible to neutron scattering. Such layer-resolved charge and orbital fluctuations are instead more naturally accessed by spectroscopic probes sensitive to charge and orbital dynamics, for example resonant inelastic x-ray scattering (RIXS) [96]. To make this distinction quantitative and establish a direct link to measurable response functions, we compute the dynamical structure factors (DSF) associated with both spin and layer pseudospin operators. For a local operator ${\hat{O}}_{i}$ (taken to be either the spin operator ${\hat{{\boldsymbol{S}}}}_{i}$ or the layer pseudospin operator ${\hat{{\boldsymbol{\tau }}}}_{i}$), the DSF is defined as
$\begin{eqnarray}\begin{array}{rcl}{S}_{O}({\boldsymbol{k}},\omega ) & = & \displaystyle \frac{1}{{L}^{2}}\displaystyle \sum _{i,j}{{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot ({{\boldsymbol{r}}}_{i}-{{\boldsymbol{r}}}_{j})}\\ & & \times \int {\rm{d}}t{{\rm{e}}}^{-{\rm{i}}\omega t}\langle {\hat{O}}_{i}(t){\hat{O}}_{j}(0)\rangle .\end{array}\end{eqnarray}$
The DSF computed for both spin and layer pseudospin operators in the four different phases are presented in figure 5.
Figure 5. Spin and layer dynamical structure factors along the high-symmetry lines for the same parameters as in figure 4 in different phases. (a), (b) FM-LS phase, (c), (d) AFM-LS phase, (e), (f) AFM-LC phase, and (g), (h) SLE phase. In each pair, the upper (lower) panel shows the spin (layer) dynamical structure factor. |
5. Discussions and conclusions
In this work, we have theoretically investigated the low-energy physics of a quarter-hole-filled bilayer Hubbard model, particularly focusing on the role of orbital-selective interlayer hybridization. By explicitly treating the strong vertical bonding of ${d}_{{{z}}^{2}}$ orbitals within a MO basis and projecting out high-energy states, we derived an effective anisotropic Kugel–Khomskii Hamiltonian that describes the coupled dynamics of the electron spin and an emergent layer pseudospin associated with electrons in the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals. The resulting effective model is highly anisotropic in the pseudospin space and exhibits a nontrivial interplay between spin exchange and spin–layer coupled interactions.
Using a combination of Weiss mean-field theory that fully retains on-site quantum correlations between spin and layer degrees of freedom, we established the ground state phase diagram characterized by four distinct quantum phases. In addition to conventional magnetically ordered phases accompanied by layer-staggered electron occupation, we identified a layer-coherent antiferromagnetic phase characterized by spontaneous interlayer quantum coherence without static layer polarization. More remarkably, we uncovered a SLE phase that does not exhibit any conventional dipolar order in either the spin or layer sectors. Instead, this phase is characterized by a hidden composite order parameter formed by the bilinear operators ${\hat{S}}^{\beta }{\hat{\tau }}^{\alpha }$, signaling maximal local entanglement between the two internal degrees of freedom.
We showed that the SLE phase originates from a strong and cooperative coupling between the spin and layer degrees of freedom, leading to the spontaneous breaking of an emergent O(4) symmetry. This symmetry breaking gives rise to a distinctive excitation spectrum featuring the entangled gapless Goldstone modes, which correspond to collective rotations of a rigid spin–layer-locked order parameter rather than independent fluctuations of the individual sectors.
Experimentally, the spin and layer pseudospin excitations can be distinguished using inelastic neutron scattering and RIXS, respectively. Furthermore, the static checkerboard charge order in the FM-LS and AFM-LS phases offers unique signatures: the induced electric dipoles generate characteristic Terahertz absorption peaks, and the spatial modulation is directly observable via STM. In contrast, the SLE and AFM-LC phases preserve layer symmetry and lack these static charge features, allowing for clear differentiation.
Recent experiments on the bilayer nickelate La3Ni2O7 have reported signatures of possibly intertwined spin-density-wave (SDW) and charge-density-wave (CDW) orders [97–100]. At ambient pressure, the low-temperature phase of La3Ni2O7 adopts the Amam space group, and the resulting spin or charge pattern cannot be stabilized within the unit cell of the minimal Kugel–Khomskii model we studied in this work, which considers only nearest-neighbor hopping processes. However, we expect the effective Kugel–Khomskii model is able to describe the intertwined SDW and CDW order when additional effects, such as longer-range hopping, residual itinerancy, and electron–lattice coupling [101–105], are taken into account.
Moreover, it is worth noting that the strong coupling between spin and layer degrees of freedom revealed in our study may provide a natural setting for the coexistence or mutual reinforcement of spin and charge ordering tendencies. In particular, ordering in the layer (pseudospin) sector necessarily involves a modulation of electronic occupation between the two layers, which can be viewed as a form of interlayer charge density modulation. When combined with magnetic ordering in the spin sector, such layer-selective charge redistribution may naturally accompany or enhance SDW instabilities, leading to intertwined spin and charge ordering phenomena. From this perspective, certain CDW signatures observed experimentally may not originate from a conventional Peierls-type instability driven purely by Fermi surface nesting, but could instead reflect charge modulations tied to spin-layer correlations and orbital-selective interlayer hybridization. While the minimal Kugel–Khomskii model considered here does not incorporate the lattice anisotropy and spatially modulated interactions required to stabilize SDW or CDW order at the observed wave vectors, it suggests that spin and charge degrees of freedom are intrinsically coupled in bilayer nickelates. In more realistic settings that include structural distortions and longer-range interactions, this intrinsic coupling may facilitate the emergence of coexisting or intertwined SDW and CDW orders.
Beyond La3Ni2O7, the mechanism identified in this work is expected to be applicable to a broader class of correlated systems. In particular, moiré bilayers and cold-atom optical lattices offer highly tunable platforms in which layer, orbital, or valley indices play roles analogous to the layer pseudospin considered here. In these settings, interlayer tunneling, interaction anisotropy, and filling can be controlled independently, providing promising opportunities to engineer and probe spin–layer entangled phases in a controlled manner.
Finally, we comment on the potential influence of strong quantum fluctuations beyond the mean-field level. While our current analysis identifies the ordered SLE phase arising from the emergent O(4) symmetry, the high symmetry and the intrinsic frustration between spin and layer sectors place the system in a regime susceptible to intense quantum fluctuations. In the strictly two-dimensional limit, these fluctuations could potentially melt the long-range composite order, leading to the realization of an exotic quantum spin-orbital liquid. To rigorously explore this possibility and determine the true quantum ground state, sophisticated numerical approaches capable of capturing long-range entanglement are required. In particular, tensor network methods, such as Projected Entangled Pair States, would be powerful tools for future investigations to address the competition between the SLE order and potential liquid states.
Appendix A Schrieffer–Wolff transformation and second perturbation
To derive the effective model describing the exchange interactions between these local moments, we employ the Schrieffer–Wolff transformation. The effective Hamiltonian is generated via a unitary transformation Heff = esHe−s, which can be expanded as:
$\begin{eqnarray}\begin{array}{rcl}{H}_{\,\rm{eff}\,} & = & {{\rm{e}}}^{s}H{{\rm{e}}}^{-s}={H}_{0}+[S,\,{H}_{0}]\\ & & +{H}^{{\prime} }+[S,\,{H}^{{\prime} }]\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{r}+\frac{1}{2}[S,\,[S,\,{H}_{0}]]+...\end{array}\end{eqnarray}$
This formalism allows us to systematically eliminate the high-energy degrees of freedom and construct an effective Hamiltonian Heff that acts solely within the low-energy subspace ${ \mathcal S }$. The unitary operator S is determined by the condition that the first-order terms vanish, i.e. ${H}^{{\prime} }+[{H}_{0},\,S]=0$. By solving this commutation relation, the operator S can be expressed in the general form:
$\begin{eqnarray}\begin{array}{r}S=\displaystyle \sum _{ij,\sigma {\sigma }^{{\prime} },\eta {\eta }^{{\prime} }}{A}_{ij,\sigma {\sigma }^{{\prime} }}^{\eta {\eta }^{{\prime} }}{d}_{i\sigma }^{\eta \dagger }{d}_{j{\sigma }^{{\prime} }}^{{\eta }^{{\prime} }},\end{array}\end{eqnarray}$
where the coefficients ${A}_{ij,\sigma {\sigma }^{{\prime} }}^{\eta {\eta }^{{\prime} }}$ are determined by the matrix elements of the perturbation $H^{\prime} $ between the eigenstates of the unperturbed Hamiltonian H0. Specifically, the tensor A is given by: $\begin{eqnarray}\begin{array}{r}{A}_{mn}=\frac{\langle m| H^{\prime} | n\rangle }{{E}_{m}-{E}_{n}},\end{array}\end{eqnarray}$
where ∣m⟩ and ∣n⟩ denote the local eigenstates of H0 with corresponding eigenvalues Em and En, belonging to the low-energy and high-energy subspaces, respectively. In our bilayer two-orbital system, H0 incorporates the robust interlayer ${t}_{\perp }^{{zz}}$ between ${d}_{{{z}}^{2}}$ orbitals (forming bonding and antibonding states), and the full multiorbital interactions including the Hund's coupling and pair-hopping terms. Due to the intricate mixing of these degrees of freedom, the basis states ∣m⟩ are non-trivial superpositions of the original orbital configurations. Consequently, these eigenstates and their associated energy levels cannot be obtained analytically and are instead determined through exact numerical diagonalization of H0.To evaluate the explicit form of the second-order terms, it is convenient to introduce the projection operators P and Q = 1 − P, which project onto the low-energy manifold and the excited states, respectively. Since the hopping amplitude t∥ is much smaller than the energy gap ΔE to the excited states, we can treat ${H}^{{\prime} }$ as a perturbation. The effective Hamiltonian up to the second order is given:
$\begin{eqnarray}{H}_{\,\rm{eff}\,}=P{H}^{{\prime} }P-P{H}^{{\prime} }Q\frac{1}{{H}_{0}-{E}_{0}}Q{H}^{{\prime} }P.\end{eqnarray}$
Since the hopping term ${H}^{{\prime} }$ changes the local particle number, it has no matrix element within the fixed-filling ground state manifold. Consequently, the leading contribution arises from the second-order term, which describes virtual hopping processes.
There are two major virtual hopping processes shown in figure 2. In the Type-I process, the electrons at sites i and j reside in the same layer (τi,z = τj,z) but possess antiparallel spins. The hopping event generates a high-energy intermediate state with on-site double occupancy. Subsequently, the system returns to the low-energy subspace via a second in-plane hopping. Throughout this process, the pseudospin configuration remains invariant (τz is conserved). However, the final spin arrangement has two possibilities: the spins can either recover their original configuration or undergo a spin exchange (spin-flip), as illustrated in figure 2(a).
In the Type-II process, the electrons at sites i and j initially occupy different layers (τi,z ≠ τj,z), while their spin orientation can be arbitrary. The initial hopping generates an intermediate state where the acceptor site accommodates electrons in both the top and bottom layers. During the return hop, either of the two electrons on the doubly occupied site can hop back to the donor site. This leads to multiple possible outcomes: the system may return to its original configuration, or undergo a spin flip, a pseudospin flip, as illustrated in figure 3(b). For this process, we must address the origin of the high excitation energy of the intermediate state compared to the initial configuration. The primary reason lies in the formation of a spin singlet within the doubly occupied bonding ${d}_{{{z}}^{2}}$ orbital. When the ${d}_{{{x}}^{2}-{{y}}^{2}}$ orbitals are one in each layer, the presence of this rigid ${d}_{{z}^{2}}$ singlet suppresses the effective inter-orbital Hund's coupling. Unlike in a high-spin atomic limit where electrons align to lower the energy, the z2 singlet prevents ferromagnetic alignment with the ${d}_{{{x}}^{2}-{{y}}^{2}}$ electrons. Consequently, the system is unable to gain the Hund's correlation energy, resulting in a significantly higher total energy.
Appendix B The order parameters in different phases
The order parameters in different phases are shown in figure B.1. Here we define the averaged local moments in the spin and pseudospin sectors as $\langle {O}_{s}\rangle =\frac{1}{{L}^{2}}{\sum }_{i,\alpha }{\langle {S}_{i}^{\alpha }\rangle }^{2}$ and $\langle {O}_{\tau }\rangle =\frac{1}{{L}^{2}}{\sum }_{i,\alpha }{\langle {\tau }_{i}^{\alpha }\rangle }^{2}$, respectively. In the FM-LS and AFM-LS phases, both spin and pseudospin moments are finite, as shown in figures B.1(a), (b). In the SLE phase, both moments are strongly suppressed, indicating the absence of conventional dipolar order in either sector. Figures B.1(c), (d) shows the momentum-space distribution of spin-layer correlations ${S}_{\mu }({\boldsymbol{k}})=\tfrac{1}{{L}^{2}}{\sum }_{i}{{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot {{\boldsymbol{r}}}_{i}}{\sum }_{\alpha ,\beta }{{ \mathcal Q }}_{i}^{\beta \alpha }$, where β ∈ {x, y, z}. The transverse component (μ = ⊥) showa a Bragg peak at k = (π, π), and the longitudinal component (μ = ∥) shows a peak at k = (0, 0).
Figure B.1. (a), (b) Averaged local moments in the spin and pseudospin sectors along the dashed line in figure 3(a), defined as $\langle {O}_{s}\rangle =\frac{1}{{L}^{2}}{\sum }_{i,\alpha }{\langle {S}_{i}^{\alpha }\rangle }^{2}$ and $\langle {O}_{\tau }\rangle =\frac{1}{{L}^{2}}{\sum }_{i,\alpha }{\langle {\tau }_{i}^{\alpha }\rangle }^{2}$, respectively. (c), (d) Momentum-space distribution of spin-layer correlations ${S}_{\mu }({\boldsymbol{k}})=\tfrac{1}{{L}^{2}}{\sum }_{i}{{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot {{\boldsymbol{r}}}_{i}}{\sum }_{\alpha ,\beta }{{ \mathcal Q }}_{i}^{\beta \alpha }$, where β ∈ {x, y, z}. (c) shows the transverse component (μ = ⊥) summing over α ∈ {x, y}, and (d) shows the longitudinal component (μ = ∥) with α = z. |
Appendix C SU(4) flavor wave theory
Given that the local Hilbert space has dimension Dl = 4, we introduce four Schwinger bosons with annihilation (creation) operators ${b}_{m,i}\,({b}_{m,i}^{(\dagger )})$, with the flavor index m ∈ {0, 1, 2, 3}.
The eigenstates of ${S}_{i}^{z}$ are expressed in terms of these bosons as ${b}_{m,i}^{\dagger }| {\rm{\varnothing }}\rangle =| m{\rangle }_{i}$. The local constraint ${\sum }_{m=0}^{3}{b}_{m,i}^{\dagger }{b}_{m,i}=1$ projects the bosonic operators onto the physical Hilbert space. To describe the fluctuations above the ordered ground state obtained from our variational calculation, we perform a site-dependent unitary transformation to align the local quantization axis with the classical ground state direction. We introduce a new set of bosonic operators ${\tilde{b}}_{n,i}$ (n ∈ {0, 1, 2, 3}), which are related to the original operators bm,i via a local unitary matrix ${{ \mathcal U }}_{i}$:
$\begin{eqnarray}{b}_{m,i}=\displaystyle \sum _{n=0}^{3}{({{ \mathcal U }}_{i})}_{mn}{\tilde{b}}_{n,i}.\end{eqnarray}$
We then treat the system within the Holstein–Primakoff approximation. Assuming the ground state is macroscopically occupied, we replace the operator ${\tilde{b}}_{0,i}$ with a classical variable by condensing the boson in the n = 0 channel:
$\begin{eqnarray}{\tilde{b}}_{0,i}={\tilde{b}}_{0,i}^{\dagger }\approx \sqrt{1-{\sum }_{n=1}^{3}{\tilde{b}}_{n,i}^{\dagger }{\tilde{b}}_{n,i}}.\end{eqnarray}$
By substituting these expressions back into the original Hamiltonian and retaining terms up to quadratic order, we obtain the spin-wave Hamiltonian:
$\begin{eqnarray}{H}_{\,\rm{SW}\,}=\frac{1}{2}\displaystyle \sum _{{\boldsymbol{k}}}{{\boldsymbol{\Psi }}}_{{\boldsymbol{k}}}^{\dagger }{ \mathcal H }({\boldsymbol{k}}){{\boldsymbol{\Psi }}}_{{\boldsymbol{k}}}+\,\rm{const}\,{,}\end{eqnarray}$
where $\Psi$k is the Nambu spinor containing the Fourier-transformed boson operators.This quadratic Hamiltonian is diagonalized by a Bogoliubov transformation, yielding the single-particle dispersions Ekα:
$\begin{eqnarray}{H}_{\,\rm{SW}\,}=\displaystyle \sum _{k,\alpha }{E}_{k\alpha }\left({a}_{k\alpha }^{\dagger }{a}_{k\alpha }+\frac{1}{2}\right)+{E}_{0}.\end{eqnarray}$