Exact Rokhsar-Kivelson points hosting first-order phase transitions in $S=1$ compass models

Introduction

Exactly solvable models are often helpful in understanding otherwise mysterious properties of interacting many-body systems. Familiar examples include the AKLT chain [1], the Kitaev honeycomb model [2], and quantum dimer models (QDMs) at their Rokhsar–Kivelson (RK) points [3–7]. The celebrated Kitaev model [2] is but one example of a more general class of models called compass models, in which the coupling on a bond depends on its orientation and selects a different spin component along each lattice direction [8–11]. In several compass models, this structure supports subsystem symmetries that constrain ordering and dynamics [11–17]. This can partition the Hilbert space into many dynamically disconnected sectors, with further disconnected parts possible within a sector [18–21]. Such conservation laws tied to the lattice geometry, can impose mobility constraints reminiscent of fracton systems [22, 23].

The RK points alluded to earlier were originally introduced in QDMs, where the ground states are equal-amplitude superpositions of dimer coverings connected by local resonance moves [3, 4]. Later it was realized that one can construct similar RK points in other systems like quantum vertex and coloring models, establishing a broader connection between exact quantum ground states and classical statistical mechanics [24, 25]. At an exact RK point the negative of the Hamiltonian acts as the generator of a continuous-time classical Markov process on each connected component of a graph. Consequently, each component hosts an equal-amplitude ground state, and its excitation energies coincide, up to an overall rate scale, with the relaxation rates of the corresponding classical process [24, 25].

Our work has its origins in Ref. [26] which constructed exponentially many exact ground states of a bilinear-biquadratic $S = 1$ Kitaev chain using spin fractionalization. This work of Ref. [26] raises the natural question of whether analogous exactly solvable points exist in higher-dimensional $S = 1$ models. We address this question here, and in doing so, bring together RK physics and subsystem-symmetry induced Hilbert-space decomposition in $S = 1$ compass models on the square and kagome lattices as concrete examples. The kagome compass model is known more popularly as the kagome Kitaev model. On top of the standard bond-dependent bilinear interactions in such models we also include biquadratic couplings of the same spin component on each bond, with their relative strengths parametrized by an angle $\theta$. Our approach avoids spin fractionalization and, subject to geometric compatibility conditions specified later, applies to $S = 1$ models without restricting the lattice coordination number. It not only yields exact ground-state wavefunctions but also makes a statement about the entire many-body spectrum of these models by showing the presence of exact RK points. It also applies to the $S = 1$ Kitaev chain.

In quantum magnets, RK constructions include models with effective ring-exchange dynamics in constrained low-energy manifolds [27, 28], as well as spin Hamiltonians designed to support exact RVB or vertex-model ground states [30] [29]. These approaches typically involve interactions acting on more than two sites or rely on a constrained low-energy description. Here we establish a direct route to RK physics in two-dimensional $S = 1$ compass models having only nearest-neighbour two-site bilinear and biquadratic interactions. On compatible geometries, a product of sublattice-dependent on-site unitaries brings the full Hamiltonian into an adjacency-degree representation on a configuration graph in the Cartesian basis. This representation holds for arbitrary $\theta$ and reveals exact RK points at $\theta= \pi/4$ and $3\pi/4$. The construction requires neither a low-energy projection nor additional multispin interactions and directly relates the ground-state manifolds to configuration-space connectivity.

In these compass models on compatible geometries the ground-state manifold is known exactly throughout $\pi/4 \le\theta\le3\pi/4$. In the interior of this interval, both models have two product states in the ground state manifold. Each is annihilated by every term in the Hamiltonian, whose positive semidefiniteness in this interval guarantees that these states minimize the energy. At both endpoints ($\theta= \pi/4$ and $\theta= 3\pi/4$) on the two lattices, a product of on-site unitary operators maps the Hamiltonian to the Laplacian of a configuration-graph $\mathcal{G}$ [31, 32] whose vertices represent product configurations in the $S = 1$ Cartesian basis. Each connected component of $\mathcal{G}$ supports one state in the ground-state manifold, an equal-amplitude superposition over the configurations in the connected component in the transformed basis, with configuration-dependent phases in the original basis [31, 32].

Both the kagome and square $S = 1$ compass models have $m = O(L)$ independent and mutually commuting conserved $\mathbb{Z}_2$ line charges and hence $2^m$ sectors [9, 11, 15]. Thus, the subsystem symmetries enforce at least $2^m$ connected components and hence a corresponding ground-state degeneracy exponential in the linear system size at the RK points. These ground state wavefunctions have a natural PEPS representation [33, 34, 35, 36] of bond dimension two on both the lattices. These endpoints are therefore exact RK points, analogous to those in QDMs, with the exponentially many sectors labelled by the line-charge values. The RK endpoints also host first-order transitions bounding a phase with long-range order on both the lattices.

Model

Each lattice site hosts an $S = 1$ spin, and we consider both open and periodic boundary conditions (OBC and PBC). The square lattice has a single-site unit cell, while the kagome lattice has a three-site unit cell defined as sites of an upward-pointing triangle. On the square lattice, we allow arbitrary integer dimensions $L_x, L_y \ge2$ for OBC and restrict both $L_x$ and $L_y$ to even values for PBC. On the kagome lattice, we allow arbitrary positive integer values of $L_x$ and $L_y$, both greater than 1, for both OBC and PBC ($L_x$ and $L_y$ are taken along the two primitive directions of the kagome lattice). Each bond $b$ of the two lattices is assigned an orientation-dependent label $\alpha_b \in\{x,y,z\}$ and it couples the corresponding component of the $S = 1$ spins located at its endpoints $i_b$ and $j_b$ through bilinear and biquadratic interactions:

where the sum runs over all bonds $b$ of the lattice and $\theta$ parametrizes the relative strengths of the bilinear and biquadratic terms. The square lattice contains only $x$- and $y$-type bonds, whereas the kagome lattice contains all three bond types, as illustrated in Figure 1.

Two transformations relate different regions of the parameter space in Eq. (1). First, $\hat{H}(\theta+ \pi) = -\hat{H}(\theta)$, so a shift by $\pi$ reverses the spectrum. Second, a unitary $\hat{V}$ composed of on-site $\pi$ rotations reverses the sign of every bilinear term while leaving the biquadratic terms invariant, giving $\hat{V}^{\dagger}\hat{H}(\theta)\hat{V} = \hat{H}(\pi- \theta)$. For the bond assignments in Figure 1, $\hat{V} = \prod_{i\in A}e^{i\pi\hat{S}^x_i}$ on the square lattice and $\hat{V} = \prod_{i\in A}e^{i\pi\hat{S}^x_i}\prod_{j\in B}e^{i\pi\hat{S}^y_j}$ on the kagome lattice, where $A$ and $B$ denote the sublattices. The full spectrum can therefore be reconstructed from $\theta\in[0,\pi/2]$. Within $\theta\in[0,\pi]$, the ground-state phase diagram is symmetric about $\theta= \pi/2$; in particular, the points $\theta= \pi/4$ and $3\pi/4$ are unitarily equivalent.

Cartesian basis and bond kinematics

We work throughout in the $S = 1$ time-reversal-invariant Cartesian basis $\{|x\rangle, |y\rangle, |z\rangle\}$ [37, 38], in which

In this basis, each bilinear bond operator maps a Cartesian product configuration of the spins at its endpoints to a single other product configuration or annihilates them, while each biquadratic bond operator is diagonal. Since $\hat{S}^{\alpha}|\alpha\rangle= 0$, a bond of type $\alpha$ annihilates any configuration with $|\alpha\rangle$ at either endpoint. We call such a bond blocked. The bond acts nontrivially only in the four-dimensional subspace spanned by configurations with neither endpoint in $|\alpha\rangle$ (we call this the active subspace). For a cyclic permutation $(\alpha,\beta,\gamma)$ of $(x,y,z)$, this active subspace is $\{| \beta\gamma\rangle, |\gamma\beta\rangle, |\beta\gamma\rangle, |\gamma\gamma\rangle\}$. Using Eq. (2), the bilinear terms of $\hat{H}$ generate two types of moves within this subspace:

exchange: $|\beta\gamma\rangle\leftrightarrow|\gamma\beta\rangle$, amplitude $+\cos\theta$. pairing: $|\beta\beta\rangle\leftrightarrow|\gamma\gamma\rangle$, amplitude $-\cos\theta$. {#eq:3}

The biquadratic interaction is proportional to the identity on the active subspace, contributing a diagonal energy $\sin\theta$ to each configuration for every active bond present in it.

Two consequences will be used in what follows. First, the single-bond eigenvalues are $\sin\theta\pm\cos\theta= \sqrt{2}\sin(\theta\pm\pi/4)$, each twofold degenerate, in the active subspace, and zero in the blocked subspace. Second, each active bond permits exactly one move: exchange for $|\beta\gamma\rangle$ or $|\gamma\beta\rangle$, and pairing for $|\beta\beta\rangle$ or $|\gamma\gamma\rangle$ in the Cartesian basis. This construction defines a weighted configuration graph $\mathcal{G}$ with $3^N$ vertices, one vertex for each product state in the Cartesian basis, $|s\rangle= \bigotimes_{i=1}^{N}|s_i\rangle$, where $|s_i\rangle\in\{|x_i\rangle,|y_i\rangle,|z_i\rangle\}$. The edges carry weights $\pm\cos(\theta)$, and each vertex has a weight $\sin(\theta)D(s)$ due to the biquadratic term, where $D(s)$ denotes the number of active bonds in $|s\rangle$.

Line charge induced Hilbert space decomposition

The Hamiltonian in Eq. (1) admits a family of conserved $\mathbb{Z}_2$ charges supported on straight lines, whose number scales with the linear system size [9, 11, 15]. Let $\ell$ be a maximal straight chain of bonds of a single type, and let $\tau_i$ denote the common type of the bonds incident on site $i\in\ell$ that lie outside $\ell$. The corresponding line charge is

Each rotation leaves the spin component on bonds outside $\ell$ unchanged, while reversing the component on bonds along $\ell$. Since both endpoints of every bond along $\ell$ are rotated, all bilinear and biquadratic interactions remain invariant.

On the square lattice, $\tau_i$ is constant along each line, yielding $\hat{Q}_{c}=\prod_j e^{i\pi\hat{S}_{(c,j)}^{x}}$ for columns and $\hat{Q}_{r}=\prod_i e^{i\pi\hat{S}_{(r,i)}^{y}}$ for rows [9, 11]. On the kagome lattice, $\tau_i$ alternates between the two sublattices along each line: an $x$-line ($A$ and $B$ sublattices alternate along it) carries rotations $e^{i\pi\hat{S}^{y}}$ on $A$ and $e^{i\pi\hat{S}^{z}}$ on $B$ sublattices, with the $y$- and $z$-line charges obtained by cyclic permutation [15]. The line charges satisfy $\hat{Q}_{\ell}^{2}=1$ and are diagonal in the Cartesian basis. Their eigenvalues $q_{\ell}=\pm1$ therefore partition the Hilbert space into dynamically disconnected sectors. An $L_x\times L_y$ square lattice has $L_x+L_y$ independent line charges and hence $2^{L_x+L_y}$ sectors. On the kagome lattice, the three line families provide $L_x+L_y+\gcd(L_x,L_y)$ independent charges under PBC and $2L_x+2L_y-1$ under OBC, yielding $2^{L_x+L_y+\gcd(L_x,L_y)}$ and $2^{2L_x+2L_y-1}$ sectors, respectively. The number of sectors is thus exponential in the linear system size on both the lattices.

On both the lattices the conserved line charges divide the configuration graph $\mathcal{G}$ into disconnected parts labelled by the line charge eigenvalues. Its vertices are the $3^N$ product configurations in the Cartesian basis, and its edges are the allowed moves. Each connected component lies within a line-charge sector, but any given sector can host more than a single such component. The line charges therefore give a lower bound on the number of connected components present in $\mathcal{G}$. Fully blocked configurations form isolated vertices in $\mathcal{G}$.

An exactly solvable ordered phase

For $\pi/4 < \theta< 3\pi/4$, on the active subspace, both single-bond eigenvalues $\sin\theta\pm\cos\theta$ are strictly positive. The Hamiltonian is therefore a sum of positive-semidefinite bond terms, and its zero-energy ground-state manifold is the intersection of their kernels. Each bond kernel is spanned by Cartesian-basis configurations in which that bond is blocked. Consequently, the zero-energy manifold is spanned by fully blocked configurations, and such configurations also span the ground-state space, when they exist.

On the square lattice, the two fully blocked configurations are the checkerboard arrangement of $\lvert x\rangle$ and $\lvert y\rangle$ states on the $A$ and $B$ sublattices, respectively, and its partner obtained by interchanging $\lvert x\rangle$ and $\lvert y\rangle$. On the kagome lattice, one fully blocked configuration places $\lvert x\rangle$, $\lvert z\rangle$, and $\lvert y\rangle$ on the $A$, $B$, and $C$ sublattices, respectively; its partner follows from a single application of the cyclic permutation $\lvert x\rangle\to\lvert y\rangle\to\lvert z\rangle\to\lvert x\rangle$ (More than one application leads to a state which is not blocked and hence does not lie in the ground state manifold). A schematic of these frozen configurations is given in Figure 2. Throughout this interval, both models have a two-dimensional ground-state manifold spanned by product states with long-range order. The above statements concerning frozen configurations being the ground state for $\pi/4 < \theta< 3\pi/4$ does not apply to the square lattice with PBC when either $L_x$ or $L_y$ is odd.

The RK point: as a graph laplacian

At $\theta= \pi/4$, the Hamiltonian becomes a sum of bond projectors, $\hat{H}(\pi/4) = \sqrt{2}\sum_b \hat{P}_b$, where

Its ground-state manifold is therefore the common kernel of these projectors, which also contains the two respective frozen states on each lattice introduced in § 1.4. At this point, however, the lower active-subspace eigenvalue of a bond $\sin\theta- \cos\theta$ also vanishes, enlarging each bond’s kernel and hence allowing for additional states in the ground-state manifold beyond the frozen configurations.

To make it convenient to identify the enlarged ground state manifold we perform a similarity transformation on the Hamiltonian in Eq. (1). The unitary operator is defined as,

diagonal in the Cartesian basis with all entries having unit-modulus. It is a product of on-site unitaries. The choice of $\hat{U}$ will differ between the two lattices and for each lattice it will further depend on the underlying sublattice as well. On the square lattice, we take

$g_i(x) = g_i(y) = 1$ on both $A$ and $B$ sublattices, with $g_i(z) = +i$ for $i \in A$ and $g_i(z) = -i$ for $i \in B$. On the kagome lattice the three sublattices $A$, $B$ and $C$ carry phases of:

The sublattice convention on both the lattices for the corresponding directions of chosen bond couplings is shown in Figure 1. On performing the similarity transformation on $\hat{H}$ we get,

with $g(s) = \prod_i g_i(s_i)$. It is shown in the appendix that using this similarity transformation with $\hat{U}$ composed of on-site unitaries specified earlier for the two lattices respectively, reverses the sign of every exchange amplitude and preserves that of every pairing amplitude of Eq. (3), so that both classes of moves due to the bilinear terms present in Eq. (1) (in the active subspace) acquire amplitude $-\cos(\theta)$. Combined with the fact that each active bond contributes one move in the configuration-graph of Cartesian-basis product vertices (due to the bilinear term) and one unit of diagonal amplitude equal to $\sin(\theta)$ (due to the biquadratic term), gives the central identity of this work,

where $\widetilde{H}(\mathcal{G},\theta)$ is the adjacency-degree representation of the transformed Hamiltonian on the configuration graph $\mathcal{G}$ with $3^N$ vertices. $A(\mathcal{G})$ is its adjacency matrix and $D(\mathcal{G}) = \operatorname{diag}[D(s)]$ is its degree matrix.

At $\theta= \pi/4 \equiv\theta_c$, and by unitary equivalence at $\theta= 3\pi/4$, the models realize an exact RK point, with $\widetilde{H}(\mathcal{G},\pi/4) = \mathcal{L}(\mathcal{G})/\sqrt{2}$, the graph Laplacian on $\mathcal{G}$ [31, 32]. Thus, $-\mathcal{L}(\mathcal{G})$ generates a continuous time random walk on the Cartesian-basis configurations in the transformed basis. The graph Laplacian at $\theta_c$ is positive semidefinite, and its kernel contains one equal-amplitude state per connected component of $\mathcal{G}$, which also span the ground state manifold [31, 32]. In the original basis, these ground states acquire the configuration-dependent phases specified by $\hat{U}$. Since every line-charge sector is non-empty and distinct sectors are dynamically disconnected, this implies at least $2^{O(L)}$ connected components in $\mathcal{G}$ on both the lattices and hence the ground state degeneracy is at least $2^{O(L)}$ in both cases. The frozen configurations are isolated vertices in $\mathcal{G}$.

On either lattice, a ground state in a particular line-charge sector can be constructed by projecting the product state $|P\rangle= \hat{U}\otimes_i[(|x_i\rangle+ |y_i\rangle+ |z_i\rangle)/\sqrt{3}]$ onto the corresponding line-charge eigenvalues ($\hat{U}$ is different for the two lattices, though both are composed of on-site unitaries, diagonal in the Cartesian basis). For line charge eigenvalues $\vec{q}=\{q_{\ell}\}$, this gives $\lvert\Psi^{\vec{q}}\rangle\propto\prod_{\ell}\frac{1}{2}(1+q_{\ell}\hat{Q}_{\ell})\lvert P\rangle$, where $\ell$ runs over all conserved $\mathbb{Z}_{2}$ line charges. When the sector consists of a single connected component, this projection yields its unique ground state. In sectors containing multiple such components, it yields a superposition of their ground states. These sector-projected states admit a PEPS representation [33, 35, 36]. Each site carries a tensor with one physical index of dimension three, corresponding to the three states, and four virtual indices of dimension two. The virtual indices form two pairs, one for each of the two charge lines intersecting at the site, and carry the binary variables introduced by expanding the corresponding $\mathbb{Z}_{2}$ charge projectors. Distinct sublattices have distinct tensors. A detailed discussion of the PEPS representation on both the lattices is given in the appendix.

Let $n_{\lambda}^{\mu}(s)$ denote the number of sites on sublattice $\mu$ carrying state $\lambda$ in configuration $\lvert s\rangle$, and let $n_z(s)$ be the total number of $z$ states present in $\lvert s\rangle$. The normalized ground state associated with a connected component $\mathcal{F}$ of $\mathcal{G}$ is

on the square and kagome lattices, respectively. The $(s)$ dependence of the exponents is suppressed in Eq. (8).

On both lattices the configuration amplitudes are pure phases, admitting a formal representation as square-root Boltzmann factors at imaginary inverse temperature. A similar construction in Ref. [39] projects three-color product ground states of the spin-1/2 kagome antiferromagnet onto fixed total-$S^z$ sectors. Though in that case these states can overlap, if they are in the same $S^z$ sector, and become linearly dependent, and their span need not exhaust the ground-state manifold. Here, the configuration-graph Laplacian mapping in the orthonormal Cartesian basis establishes exactly one ground state per connected component, with states from distinct components mutually orthogonal. Determining the full ground-state manifold therefore reduces to characterizing the connected components of the configuration-graph.

The sublattice-dependent unitary transformation used for the square lattice also applies to the analogous $S=1$ Kitaev chain with bilinear and biquadratic couplings [37, 38] [26], whose sites alternate between the $A$ and $B$ sublattices. In this case, the conserved quantities are bond charges that span the chain, taking the form $e^{i\pi\hat{S}_{i}^{\beta_i}}e^{i\pi\hat{S}_{i+1}^{\beta_{i+1}}}$ on a bond $b$ of type $\alpha$. A length $L$ periodic chain has $L$ conserved $\mathbb{Z}_{2}$ bond charges, defining $2^{L}$ charge sectors, with ground-state wavefunctions of the same form as those for the square lattice, given in Eq. (7). Just as the square- and kagome-lattice ground states admit PEPS representations with bond dimension two, the chain ground states admit an MPS representation with bond dimension two, presented in the Appendix. Using spin fractionalization, Ref. [26] obtained related ground-state wavefunctions for the $S = 1$ Kitaev chain.

Finally, using a variational bound, we can show that the RK endpoints host first-order phase transitions. In the product state $\lvert P\rangle$, each bond is active with probability $4/9$, giving the biquadratic bond operator of $\hat{H}$ in Eq. (1) an expectation value of $4/9$. Since $\hat{H}(\pi/4)\lvert P\rangle= 0$, the corresponding bilinear bond operator’s expectation value is $-4/9$. Thus,

for $0 \le\theta< \pi/4$, where $N_b$ is the total number of bonds. By contrast, $E_0(\theta)=0$ throughout $\pi/4 \le\theta\le3\pi/4$. The bound forces a nonzero left-hand slope at $\theta_c=\pi/4$ of the energy density $e_0=E_0/N_b$, whereas its right-hand slope vanishes, establishing a first-order transition. Unitary equivalence implies an analogous transition at $\theta=3\pi/4$. At either endpoint, the ground-state degeneracy is at least exponential in the linear system size and collapses to two upon entering the intervening frozen phase with long range order.

Exact Diagonalization

To test the analytic ground-state wavefunctions in Eqs. (7) and (8) and the predicted first-order transition at $\theta=\pi/4$, we perform exact diagonalization (ED) of $\hat{H}$ in Eq. (1) on a $2 \times4$ square lattice with PBC and $2 \times2$ kagome lattices with OBC and PBC. For these clusters, each line-charge sector contains one non-frozen connected component; though the two frozen configurations form additional isolated components within a sector. Figure 3 shows the overlap of the predicted ground-states in Eq. (7) and Eq. (8) with the ground state obtained by ED for all the sectors. In the sectors containing the frozen states, the frozen states are projected out of the numerical zero-energy subspace obtained by ED before comparison with the analytic state supported on the non-frozen component of that sector. Figure 3 shows that the overlap is unity in every sector for all three cases.

Figure 4 shows the ground-state energy density $e_0 \equiv E_0/N_b$ (where $E_0$ is the ground state energy and $N_b$ is the number of bonds present in the lattice) as a function of $\theta/\pi$ for the square and kagome lattices (here we also include square-lattice OBC data). For $\theta< \pi/4$, $e_0 < 0$ and the ground state is nondegenerate. For $\pi/4 < \theta< 3\pi/4$, $e_0 = 0$ and the ground-state manifold is spanned by two frozen states on each lattice mentioned in § 1.4, and shown in the schematic in Figure 2. The insets show a discontinuity in $de_0/d\theta$ at $\theta=\pi/4$, consistent with the first-order transition predicted in Sec. VI.

Discussion

Our results establish a connection between subsystem-symmetry induced Hilbert-space decomposition [21] and RK physics [3, 24, 25] in $S=1$ compass models on square and kagome lattices [9–11, 15], having both bilinear and biquadratic interactions, with their relative strengths parametrized by an angle $\theta$. Using a transformation built of on-site unitary operators, the Hamiltonian in Eq. (1) on both the lattices takes the adjacency-degree representation given in Eq. (6) on a configuration graph $\mathcal{G}$, whose vertices represent product states in the $S=1$ Cartesian basis. In both cases the ground state manifold is known exactly in an extended region of the parameter space. At the two exactly solvable endpoints of this interval, the Hamiltonian transforms into the Laplacian of $\mathcal{G}$ in the Cartesian basis [31, 32, 37, 38]. Ground-state counting and correspondence wavefunction construction thereby reduce to a graph-connectivity problem: each connected component of the configuration graph supports one state in the ground-state manifold, an equal-amplitude superposition over the states in the connected component in the transformed basis, with additional configuration-dependent phases in the original basis [31, 32]. Both lattices possess $m=O(L)$ conserved $\mathbb{Z}_2$ line charges, yielding $2^m$ symmetry sectors. This gives a lower bound on the number of connected components in the configuration-graph and hence on the ground state degeneracy at the RK points. The states in the ground state have a PEPS representation of bond dimension two on both the lattices.

Between the two RK endpoints, the ground-state manifold is spanned by two frozen configurations, yielding an ordered phase on both lattices. The discontinuity in the first derivative of the ground-state energy density establishes a first-order transition at each endpoint. The RK points thus lie at the boundaries of an ordered phase, reminiscent of the role of RK points in QDMs [4–7].

The adjacency-degree representation of the Hamiltonian in Eq. (6) also makes the analogy with QDMs explicit [3–7]. Active bonds play the role of flippable plaquettes: in the transformed Cartesian basis, the bilinear interaction generates transitions between configurations, while the biquadratic interaction assigns a diagonal energy to each active bond. These therefore correspond to the kinetic and potential-energy terms of a QDM, respectively. In two-dimensional QDMs on a $L \times L$ torus, integer winding numbers label $O(L^2)$ sectors on bipartite lattices, whereas winding parities label finitely many sectors on non-bipartite lattices. By contrast, the compass models considered here have $m = O(L)$ independent and mutually commuting conserved $\mathbb{Z}_2$ line charges labeling $2^m$ sectors. In either setting, additional fragments within a symmetry sector, including frozen configurations, can occur. This analogy motivates studying additional spin interactions that modify the ground-state configuration ensemble or its amplitudes, and investigating whether the resulting models can support Coulomb correlations and an emergent gauge description. Determining whether such behavior is confined to fine-tuned points or can persist throughout an extended phase remains an open question.

At the RK endpoints, our construction realizes the equal-weight case of the stochastic-matrix-form framework of Ref. [25], applied separately to each connected fragment. The transformed Hamiltonian is a graph Laplacian, so its negative directly generates a classical Markov process with a uniform stationary distribution within each fragment. Within each connected component, its relaxation rates coincide with the quantum excitation energies [24, 25]. Thus, the microscopic spin Hamiltonian acquires a direct stochastic interpretation through an on-site phase transformation alone, without any additional reweighting of configuration amplitudes. Finding exact excited states is another direction to be explored in the future.

The Laplacian representation at the RK endpoints also identifies continuous families of local Hamiltonians with exactly the same ground-state manifold. Specifically, replacing the uniform bond couplings by arbitrary strictly positive coefficients gives

where $\sigma= +1$ and $\sigma= -1$ preserve the ground-state manifolds at $\theta= \pi/4$ and $3\pi/4$, respectively. This freedom includes directional anisotropy and arbitrary positive bond disorder as long as at each bond, the biquadratic coupling equals the bilinear coupling at $\theta= \pi/4$ and its negative at $\theta= 3\pi/4$. Each Hamiltonian remains a positive sum of the same bond projectors, so their common kernel is unchanged. Equivalently, after the appropriate on-site unitary transformation, these modifications change the edge weights of the configuration-graph Laplacian while retaining its connected components and the uniform stationary distribution within each component. Further local extensions can make the transition strengths depend on nearby state occupations, provided the forward and reverse rates remain equal and strictly positive. Such terms introduce interactions involving additional sites while retaining the original two-spin moves and their configuration-dependent ground-state phases. Allowing rates to vanish can modify the ground state degeneracy as this could modify the number of connected components of the configuration graph.

The Cartesian-basis construction is independent of spatial dimension and works at fixed $S = 1$ without imposing a relation between spin and coordination number. The present construction requires compatible on-site phases that transform the Hamiltonian into an adjacency-degree representation, given in Eq. (6), on a configuration-graph. As long as at most two types of couplings meet at a site, the given construction holds, provided the subgraphs of every coupling are bipartite (for example the $S = 1$ Kitaev chain, in fact exactly the same choice of $\hat{U}$ used for the square lattice compass model works for the chain as well). In particular this construction to identify an appropriate $\hat{U}$ fails when all three bond types meet at a site, as in the $S = 1$ honeycomb Kitaev model, because the local phase constraints are inconsistent. In that case a more complicated $\hat{U}$ might exist, finding which will be a direction for the future.

Characterizing the phase of the $S = 1$ compass models on both the lattices outside the extended solvable interval remains another open direction for future work. In particular at $\theta= 0$ Eq. (6) gives $\hat{U}^{\dagger}\hat{H}(0)\hat{U} = -A$, where $A$ is the configuration-graph adjacency matrix. Thus, the question of finding the ground state energy and wavefunction at $\theta= 0$ might be inferred from studying the structure of $\mathcal{G}$ which will also be focus of future work.

More broadly, these models illustrate how RK physics can arise directly in microscopic spin Hamiltonians with only nearest-neighbour two-site interactions. A product of on-site unitaries exposes the exact adjacency-degree representation in Eq. (6) for arbitrary $\theta$, with the RK points emerging as special cases where the Hamiltonian is unitarily equivalent to a graph Laplacian. Identifying other spin models that admit such local transformations offers a route to further exact RK constructions and a possible graph-theoretic approach to their properties away from the solvable points.

Acknowledgments

This work began after the author read the arXiv preprint of Ref. [26] and explored possible extensions, including through discussions with Claude Opus 4.8. In particular, Claude Opus 4.8 provided insights into the properties of graph Laplacians and their implications for the ground-state manifold. Claude Opus 4.8 and ChatGPT Sol 5.6 assisted in writing the exact diagonalization code. These models, together with ChatGPT Sol 6.1, were also used at various stages of the project to assist with literature searches and polishing the language of the manuscript. The author thanks Ayush Sharma for collaboration on related work [38], which also motivated the present study; Kedar Damle for his guidance and encouragement; and Sanjay Moudgalya for useful comments on an earlier draft. This work was supported by the Department of Atomic Energy, Government of India, under Project Identification No. RTI-4012. The author takes full responsibility for all the results and claims.

References

Appendix

Construction of the unitary transformation

We derive the expression of $\hat{U}$ in Eq. (4) for both square and kagome lattices at once. Let $(\alpha,\beta,\gamma)$ be cyclic and consider a bond of type $\alpha$ joining sites $i$ and $j$. We will work in the $S=1$ Cartesian basis. For a product state $\bigotimes_k \lvert\phi_k\rangle$ with $\lvert\phi_k\rangle=a_k\lvert x\rangle+b_k\lvert y\rangle+c_k\lvert z\rangle$, write $u_k=\langle\beta\vert\phi_k\rangle$ and $v_k=\langle\gamma\vert\phi_k\rangle$ ($k$ runs over all the lattice sites. $a_k$, $b_k$ and $c_k$ are all non-zero and have unit modulus).

To make sure that $\hat{H}$ takes the adjacency-degree form over the configuration graph given in Eq. (6), we want to make sure that the amplitude of the two types of moves due to the bilinear terms of $\hat{H}$ on every bond on both the lattices given in Eq. (3) are exactly the same. Thus, we want both the pairing move and the exchange move to have the same sign. To enforce this we choose the transformation $\hat{U}$ such that, after the transformation, the exchange terms get an extra minus sign while the pairing terms are left invariant. Hence, we require antisymmetry under the exchange moves and equal amplitudes under the pairing moves. Using $v_k$ and $u_k$ in Eq. (3) we find,

With the ratio $r_k^\alpha\equiv v_k/u_k$ these equations become $r_j^\alpha=-r_i^\alpha$ and $r_i^\alpha r_j^\alpha=1$, whence

Thus, $\hat{H}$ takes the adjacency-degree representation on the configuration graph in the basis of product Cartesian states, if on every bond of type $\alpha$, the associated ratios $r_i^\alpha$ and $r_j^\alpha$ equal $\pm i$ and take opposite values at the two endpoints of $\alpha$. The three ratios at a site, $r_i^x=c_i/b_i$, $r_i^y=a_i/c_i$ and $r_i^z=b_i/a_i$, are not independent but obey

On the square lattice only $x$- and $y$-type bonds occur and meet at a site, so only $r^x$ and $r^y$ are constrained. Taking $a_{(m,n)}=b_{(m,n)}=1$ and $c_{(m,n)}=i(-1)^{m+n}$ gives $r^x=(r^y)^{-1}=i(-1)^{m+n}$, which reverses sign under $m\to m+1$ and under $n\to n+1$ alike, satisfying Eq. (12) on every bond (here we have used two indices $(m,n)$ to specify the coordinate of a site to make it convenient to specify the sublattice dependent phase). The result depends on the site only through the parity of $m+n$, giving the two-sublattice dependent phases quoted for the square lattice below Eq. (4) in Sec. VI.

On the kagome lattice the $A$-$B$, $B$-$C$ and $C$-$A$ bonds are of type $x$, $z$ and $y$, so each sublattice is touched by only two of the three bond types: $A$ by $x$ and $y$, $B$ by $x$ and $z$, $C$ by $z$ and $y$. Since every bond connects distinct sublattices, the alternation required by Eq. (12) is realized by a uniform assignment on each sublattice. Choosing $r_A^x=r_B^z=r_C^y=i$ fixes $r_B^x=r_C^z=r_A^y=-i$, and Eq. (13) then forces the remaining, unconstrained ratio at each site to be, $r_A^z=r_B^y=r_C^x=1$. Fixing the overall phase of each single-site state yields Eq. (5).

Verification of the Hamiltonian's adjacency-degree representation

Write $g(s)=\prod_i g_i(s_i)$ in Eq. (5), where $i$ runs over all the sites of the lattice. For an exchange move due to a bilinear term of $\hat{H}$, on a bond of type $\alpha$ joining sites $i$ and $j$, the two configurations $\lvert s\rangle$ and $\lvert s'\rangle$ differ by $(s_i,s_j)=(\beta,\gamma)\to(\gamma,\beta)$, so $(g(s))^*g(s')=(g_i(\beta)g_j(\gamma))^*g_i(\gamma)g_j(\beta)=r_i^\alpha r_j^\alpha=-1$ by Eq. (12). For a pairing move, $(\beta,\beta)\to(\gamma,\gamma)$ and $(g(s))^*g(s')=(r_i^\alpha r_j^\alpha)=+1$, due to Eq. (12). Since the bare contributions are $+\cos\theta$ and $-\cos\theta$ respectively [Eq. (3)], both transformed off-diagonal elements due to the bilinear terms of $\hat{H}$ equal $-\cos\theta$. The diagonal element due to the biquadratic term of $\hat{H}$ is $\sin\theta D(s)$, where $D(s)$ is the number of active bonds in $|s\rangle$ and each active bond contributes exactly one off-diagonal partner, so the degree of $s$ in $\mathcal{G}$ is also $D(s)$ and Eq. (6) follows.

Bond dimension two PEPS

In this appendix we give a PEPS representation of the ground-state wavefunctions of $\hat{H}$ at the RK point. At that point a ground state in a line-charge sector $\vec{q}$ on either the kagome or the square lattice (or a superposition of states in the ground state manifold when the sector contains more than one connected component) is obtained by projecting $|P\rangle$ onto that sector:

Here $\ell$ runs over the conserved $\mathbb{Z}_2$ line charges $\hat{Q}_{\ell}$, $q_{\ell}\in\{-1,1\}$ specifies the desired eigenvalue of each charge, and

The sum runs over all $3^N$ Cartesian-basis product configurations $s=(s_1,\ldots,s_N)$, and $g(s)=\prod_i g_i(s_i)$ contains the on-site phases specified in Sec. VI of the main-text for the two lattices.

For a fixed configuration $s$, the numerator of each line-charge projector can be expanded as

Here $Q_{\ell}(s)$ is the eigenvalue of $\hat{Q}_{\ell}$ on $|s\rangle$, $\tau_{\ell i}$ denotes the rotation component used by line $\ell$ at site $i$, and we have used

Each site of either lattice lies at the intersection of two $\mathbb{Z}_2$ charge lines. We introduce two virtual indices $a_j,b_j\in\{0,1\}$ along each incident line $j=1,2$. Writing $\tau_1$ and $\tau_2$ for the rotation components of these two lines at site $i$, we define the local tensor

The tensor has one physical index of dimension three and four virtual indices, each of dimension two. Contracting neighboring virtual indices along the lattice bonds, together with the Kronecker deltas in Eq. (18), propagates one common binary variable along each line.

For a closed line $\ell$, the desired charge $q_{\ell}$ is implemented by inserting

on any one virtual bond of that line. For an open line, the same factor is obtained by contracting the virtual indices at its two ends with $(1,1)$ and $(1,q_{\ell})$, respectively.

The resulting tensor contraction for a fixed physical configuration $s$ is

Here $\Omega_{\vec{q}}$ is the set of Cartesian-basis product configurations satisfying $Q_{\ell}(s)=q_{\ell}$ for every line $\ell$. The indicator function $\mathbf{1}_{\Omega_{\vec{q}}}(s)$ equals one when $s\in\Omega_{\vec{q}}$ and zero otherwise.

If $N_{\vec{q}}=|\Omega_{\vec{q}}|$ is the number of configurations in $\Omega_{\vec{q}}$ then multiplying the contraction in Eq. (20) by $(2^m\sqrt{N_{\vec{q}}})^{-1}$ gives the normalized wavefunction amplitude of $|s\rangle$. Summing these amplitudes over all configurations reproduces the normalized projected state in Eq. (14):

For the square lattice, take horizontal bonds of type $x$ and vertical bonds of type $y$. The row and column charges use rotations $e^{i\pi S^y}$ and $e^{i\pi S^x}$, respectively. With sites labelled by integers $(m,n)$, a compatible phase choice is

The local tensor is therefore

where $l,r,u,d\in\{0,1\}$ label the left, right, upper, and lower virtual legs.

For the kagome lattice, choose the bond convention $AB:x$, $AC:y$, and $BC:z$. Compatible on-site phases are

Writing $u_{\alpha},v_{\alpha}$ for the two virtual legs along a line of bond type $\alpha$, the tensors at a site $i$ on the respective sublattices are

Here $A$, $B$, and $C$ denote the three kagome sublattices shown in Figure 1(b).

MPS for the $S=1$ Kitaev chain

Consider a periodic even length $S=1$ Kitaev chain with $N$ sites, labelled $i=0,\ldots,N-1$. Bond $i$ connects sites $i$ and $i+1$ and has type $\alpha_i=x$ for even $i$ and $\alpha_i=y$ for odd $i$. We denote the complementary component in $\{x,y\}$ by $\beta_i$. The conserved bond charges are

where $\vec{q}$ labels the sector specified by the bond-charge eigenvalues $\{q_i\}$. In the Cartesian basis, define $w_\beta(s)$ through

The sublattice-dependent on-site phases are

where $s=(s_0,\ldots,s_{N-1})$ denotes a Cartesian-basis configuration.

For each physical index $s_i\in\{x,y,z\}$, introduce a $2\times2$ matrix $A_i^{s_i}$ with virtual indices $u,v\in\{-1,+1\}$:

Contracting the right virtual index at site $i$ with the left virtual index at site $i+1$ imposes $q_iw_{\beta_i}(s_i)=w_{\beta_i}(s_{i+1})$, or equivalently $w_{\beta_i}(s_i)w_{\beta_i}(s_{i+1})=q_i$. Consequently,

The amplitude therefore equals $g(s)$ for configurations in the sector $\Omega_{\vec{q}}$ and vanishes otherwise. For a sector, the normalized state consequently has the exact bond-dimension-two MPS representation