Analytical Solutions of Symmetric Isotropic Spin Clusters Using Spin and Point Group Projectors

Abstract

Spin models like the Heisenberg Hamiltonian effectively describe the interactions of open-shell transition-metal ions on a lattice and can account for various properties of magnetic solids and molecules. Numerical methods are usually required to find exact or approximate eigenstates, but for small clusters with spatial symmetry, analytical solutions exist, and a few Heisenberg systems have been solved in closed form. This paper presents a simple, generally applicable approach to analytically solve isotropic spin clusters, based on adapting the basis to both total spin and point group symmetry to factor the Hamiltonian matrix into sufficiently small blocks. We demonstrate applications to small rings and polyhedra, some of which are straightforward to solve by successive spin-coupling for Heisenberg terms only; additional interactions, such as biquadratic exchange or multi-center terms necessitate symmetry adaptation.

1. Introduction

Spin Hamiltonians are used to model the properties of exchange-coupled magnetic solids and molecules, particularly those of the first transition-metal series [1], such as magnetic susceptibilities and heat capacities or magnetic-resonance and neutron-scattering spectra. Based on matrix diagonalization, all physical quantities within the framework of the model can be calculated exactly. While the Hilbert space grows exponentially with the number of centers, the system size limit for numerically exact or quasi-exact calculations [2] can be increased by taking advantage of symmetries [3,4,5,6,7,8,9,10,11]. However, for simulating large systems, various approximations are essential [12].

On the opposite side of the range, small clusters allow closed-form solutions, especially if the Hamiltonian is isotropic and invariant with respect to spin permutations [8] corresponding to spatial symmetries of the molecule or cluster. In addition to their pedagogical value, analytical solutions yield insights that might be obscured or unavailable in numerical data. For example, they can facilitate parametric plots of spectra and provide exact expressions for phase boundaries in parameter space to map precise quantum-phase diagrams, which would only be approximated in numerical computations. Analytical insights can also assist in aligning the spin model with experimental observations, potentially reducing the need for numerical explorations.

Some Heisenberg systems are trivially integrable—solved without diagonalization or explicit adaptation to spatial symmetries—by Kambe’s coupling method [13], which relies on successively forming subsystem spins in a hierarchical manner to ultimately produce total spin multiplets. However, not all analytically solvable cases can be subjected to this approach, as it requires certain conditions on the coupling topology [14]. Moreover, as explained in Section 3, Kambe’s method ceases to be applicable when the model is extended to include additional isotropic terms, such as biquadratic exchange or multi-site interactions.

Here we find closed-form solutions for small isotropic clusters by exploiting spin and point group (PG) symmetries to factorize the Hamiltonian into blocks, with each block corresponding to a specific irreducible representation (irrep) of SU(2) (quantum number S) and the point group (irrep label $\Gamma$). If the size of a block does not exceed $4\times 4$, then closed-form solutions for eigenvalues (energies) are guaranteed to exist and have already been obtained for some Heisenberg clusters [15,16,17,18]. Specifically, for $s=\frac{1}{2}$ rings with $N=5,6,7$ sites [15,16], simultaneous adaptation to the z-component of spin (magnetization quantum number M) and the cyclic point group $C_N$ was sufficient to obtain subspaces of manageable sizes. For $s=\frac{1}{2}$, $N=8$ or $s=1$, $N=5$, adaptation to total spin (S and M) and $C_N$ was achieved by a recursive technique designed for rings [17]. Finally, the ground state of the antiferromagnetic $s=\frac{1}{2}$ Heisenberg icosahedron was derived in a similar manner as presented here, but without explaining the method [18]. Our purpose is thus to provide a clear and easy-to-follow procedure to diagonalize any isotropic cluster that allows analytical solutions.

In the upcoming Section 2, we briefly revisit the symmetries of isotropic spin models and discuss various existing approaches for partitioning the Hamiltonian. We then explain our strategy of setting up a generalized (non-orthogonal) eigenvalue problem within a selected subspace by applying Löwdin’s spin projector and a PG projector to random states in an uncoupled basis. Our priority is to provide a procedure that is as simple as possible to implement, rather than being the most computationally efficient. Additional practical advice on implementing PG symmetry to ensure the paper is self-contained is given in the Appendix A.

For small rings and polyhedra with various local spin values s, Section 3 tabulates the dimensions of subspaces, provides selected energy expressions and boundary conditions in parameter space, and explores the ground state as a function of independent parameters. General isotropic spin Hamiltonians can involve a multitude of free parameters [19] making extensive tabulations of spectra or derived properties impractical. Our results are not directly intended to provide new insights into any exchange-coupled cluster but should assist in verifying independent implementations of the analytical diagonalization process, which could include additional terms. However, similar systems to some of those we address here, such as square or tetrahedral configurations of spin centers (see, e.g., [1,20] and references cited therein), do exist, and analytical approaches could be useful for analyzing their properties.

2. Theory

Symmetries of isotropic Hamiltonians. In the Heisenberg model, pairwise interactions are parametrized by coupling constants $J_{ij}$, Equation (1),

$${\widehat{H}}_J={\displaystyle \sum _{i<j}{J}_{ij}{\widehat{s}}_i}·{\widehat{s}}_j,$$

(1)

where ${\widehat{s}}_i=({\widehat{s}}_{x,i},{\widehat{s}}_{y,i},{\widehat{s}}_{z,i})$ is the local spin vector of site i. Biquadratic exchange is another isotropic term, Equation (2),

$${\widehat{H}}_K={\displaystyle \sum _{i<j}{K}_{ij}}{({\widehat{s}}_i·{\widehat{s}}_j)}^2$$

(2)

Note that ${({\widehat{s}}_i·{\widehat{s}}_j)}^2$ is a linear combination of scalar couplings of local spin operators of spherical tensor rank 1 (a Heisenberg-type contribution) and rank 2 [19]. The construction of rank-2 operators requires $s>\frac{1}{2}$, and therefore ${\widehat{H}}_J$ is the only isotropic pairwise interaction for $s=\frac{1}{2}$. However, for $N\ge 4$, multi-center terms occur, see Section 3.

Isotropy means invariance with respect to spin rotations [group SU(2)], due to commutation of the Hamiltonian with all components of the total spin $\widehat{S}={\displaystyle {\sum }_i{\widehat{s}}_i}$, $[\widehat{H},{\widehat{S}}_{\alpha }]=0$ $(\alpha =x,y,z)$. Each level of an isotropic $\widehat{H}$ is a multiplet encompassing $2S+1$ states with z-projections (${\widehat{S}}_z$ eigenvalues) ranging from $M=-S$ to $M=+S$. All states of a single multiplet have an ${\widehat{S}}^2$ eigenvalue of $S(S+1)$. The basis can be spin-adapted by successive coupling, and the Hamiltonian matrix in a subspace with definite S is computed based on irreducible tensor techniques, as explained in detail elsewhere [1,21]. In the frame of exact diagonalization of the Heisenberg model, the resulting reduction in matrix sizes is highly useful, and computational packages make such calculations accessible for studying various magnetic properties [22].

In addition to spin symmetry, the Heisenberg model, and indeed any isotropic spin model, is symmetric under permutations of sites according to the spatial symmetries of the cluster [8]. This spin-permutational symmetry (SPS) is often referred to as point group (PG) symmetry. However, it is important to note that not all distinct PG symmetries of the electronic Hamiltonian are necessarily reflected in the isotropic spin model. For example, a planar hexanuclear cluster belonging to the molecular point group $D_{6h}$ could be represented as an $N=6$ Heisenberg ring, with the latter model exhibiting only $D_6$ SPS. The full group $D_{6h}=D_6\otimes {\sigma }_h$ would pertain to an anisotropic spin Hamiltonian (not considered here) that more completely represents the physics by including the consequences of spin–orbit coupling [23]; some of the group operations would then represent combinations of spin permutations and spin rotations [24,25,26].

Combining total spin (${\widehat{S}}^2$ and ${\widehat{S}}_z$) with PG is significantly more complex than using either of these two symmetries separately. It is usually impossible to successively couple individual sites into larger subsystems and ultimately into a total spin multiplet in a way that is compatible with the full point group, making demanding transformations between different coupling schemes unavoidable [8] (which are still manageable under specific circumstances [10,11]). Consequently, the application of PG symmetry is frequently limited either to a compatible subgroup [8] or—far more commonly—the full PG symmetry is utilized only in conjunction with ${\widehat{S}}_z$ (instead of ${\widehat{S}}^2$ and ${\widehat{S}}_z$) by working in an uncoupled basis $\left|m_1,\dots ,m_N\right.〉$ of definite local z-projections, $M={\displaystyle {\sum }_i{m}_i}$ [8,9,27,28]. For a concise practical explanation of the latter strategy, see Ref. [9]. We briefly mention that an alternative technique for complete adjustment to full spin and PG symmetry relies on concepts from valence bond theory but has not been widely adopted [29]. Finally, for a unitary and symmetric group approach for spin-1/2 systems, see, Ref. [30].

In contrast, we combine a PG projection operator (see below) with Löwdin’s projector [31] for full symmetry adaptation, with the aim of sufficiently reducing the dimensions of Hamiltonian blocks to enable analytical diagonalization. When used on a random state with definite M, Löwdin’s projector, Equation (3),

$${\widehat{P}}_S={\displaystyle \prod _{l\ne S}\frac{{\widehat{S}}^2-l(l+1)}{S(S+1)-l(l+1)}},$$

(3)

affords a pure-spin state $\left|S,M\right.〉$; all other contributions $(l\ne S)$ are eliminated. A similar approach (also in conjunction with spatial symmetry) has occasionally been applied in numerical calculations, e.g., in Lanczos exact diagonalization for triangular-lattice cluster models [6] but was apparently not yet employed to obtain analytical solutions.

The PG projector ${\widehat{P}}_{\lambda }^{\Gamma }$ for irrep $\Gamma$ and component $\lambda$ (the latter must be specified for multi-dimensional irreps, $d_{\Gamma }>1)$ is defined in Equation (4),

$${\widehat{P}}_{\lambda }^{\Gamma }=\frac{d_{\Gamma }}{h}{\displaystyle \sum _{g=1}^h{[D_{\lambda \lambda }^{\Gamma }(g)]}^{\ast }\widehat{G}(g)},$$

(4)

where h is the order of the group (the total number of elements g), $D_{\lambda \lambda }^{\Gamma }(g)$ is a diagonal entry of the irrep matrix $D^{\Gamma }(g)$, and $\widehat{G}(g)$ is the respective symmetry operation in spin space; the asterisk (*) denotes complex conjugation. Technical details on the practical construction of the PG projector are provided in Appendix A.

Generalized eigenvalue problem. Symmetry projectors are idempotent, ${\widehat{P}}^2=\widehat{P}$, and self-adjoint, ${\widehat{P}}^{†}=\widehat{P}$, where the dagger ${(}^{†})$ denotes the Hermitian adjoint (complex conjugate transpose). Hence, their eigenvalues can only be 0 or 1. The dimension d of the respective subspace is given by the trace, $d=\mathrm{Tr}(\widehat{P})$. We calculate the Hamiltonian and overlap matrices, h and s (the latter is not to be mixed up with a spin vector), in a space $R=(r_1,\dots ,r_d)$ of state vectors $r_i$ comprising small random integers in a symbolic representation,

$$h=R^{†}H{P}_S{P}_{\lambda }^{\Gamma }R,$$

(5)

$$s=R^{†}{P}_S{P}_{\lambda }^{\Gamma }R,$$

(6)

where H, $P_S$ and $P_{\lambda }^{\Gamma }$ are the Hamiltonian, spin and PG projector representations, respectively, in the uncoupled basis $\left|m_1,\dots ,m_N\right.〉$ for the selected magnetization M and $d(S,M,\Gamma ,\lambda )=\mathrm{Tr}(P_S{P}_{\lambda }^{\Gamma })$. The generalized eigenvalue problem,

$$hv=Esv,$$

(7)

where v is an eigenvector, has real eigenvalues E (energy levels). When solving Equation (7) with symbolic computer algebra packages like Mathematica (Eigenvalues[h,s]) or the MATLAB symbolic toolbox (eig(h, s)), the energy expressions in spaces with $d=3$ or $d=4$ are usually very long, even when explicitly assuming symbolic Hamiltonian parameters to be real, and the built-in functions for algebraic simplification might not always produce significantly shorter forms. We observed that it is sometimes possible to obtain more concise results by a similarity transformation of h, Equation (8),

$$\stackrel{⌣}{h}=L^{-1}h{(L^{-1})}^{†},$$

(8)

based on the Cholesky decomposition, $s=L{L}^{†}$; $\stackrel{⌣}{h}$ has the same eigenvalues as the original problem, but, in contrast to $s^{-1}h$, it is Hermitian. Still, most solutions of cubic $(d=3)$ or quartic $(d=4)$ polynomial equations are impracticably lengthy functions of the parameters. Therefore, Section 3 presents only a few illustrative and reasonably concise results for brevity.

Additional considerations. The only prerequisites for following our recipe are symbolic representations of the ${\widehat{s}}_i$ and the generators C for site permutations with their corresponding irrep matrices $D^{\Gamma }$. Instead of constructing the ${\widehat{s}}_i$, one can directly compute the scalar products ${\widehat{s}}_i·{\widehat{s}}_j$ of all relevant pairs in a magnetization subspace and build the projectors and all model terms considered in this paper (Heisenberg, biquadratic and four-center terms) from them (typically, $M=0$ or $M=\frac{1}{2}$, which encompasses all multiplets). This avoids working in the full Hilbert space throughout.

A symbolic calculation of $P_S$ (with a significant fraction of non-zero elements) can become a bottleneck. Therefore, one may choose to first build the $(M,\Gamma ,\lambda )$ basis. A direct full diagonalization of $P_{\lambda }^{\Gamma }$, kee** only the eigenvectors with eigenvalue 1, is not always feasible with symbolic computer algebra, but is indeed not required, because the $(M,\Gamma ,\lambda )$ space can be generated by scanning the rows of $P_{\lambda }^{\Gamma }$ and selecting only the first column with a non-zero entry, discarding all other columns of $P_{\lambda }^{\Gamma }$ that have a non-zero entry in the given row. (The described construction of a $(M,\Gamma ,\lambda )$ subspace does not explicitly require the $P_{\lambda }^{\Gamma }$ matrices and could instead be achieved as detailed in Ref. [9]. However, we believe that our current method, which applies (at least conceptually) a combined PG and spin projector to random states, offers more pedagogical clarity and would be slightly simpler to implement. Forming $P_{\lambda }^{\Gamma }$ in a symbolic representation usually does not pose a significant computational cost for systems that are small enough to have closed-form solutions.) The idea behind this procedure is that each uncoupled state $\left|m_1,..,m_N\right.〉$ appears in at most one distinct state in the $(M,\Gamma ,\lambda )$ basis. (As noted in Ref. [25], this is not true for all multidimensional irreps in all point groups if one does not separate components $\lambda$ but instead uses a simplified projector, ${\widehat{P}}^{\Gamma }=(d_{\Gamma }/h){\displaystyle {\sum }_g{\chi }^{\ast }(g)\widehat{G}(g)}$, based on characters, $\chi (g)=\mathrm{Tr}[D^{\Gamma }(g)]$, that summarily includes all components $\lambda$ of a multidimensional irrep $\Gamma$.) The thus selected columns of $P_{\lambda }^{\Gamma }$ are subsequently normalized and collected in a rectangular matrix $p_{\lambda }^{\Gamma }$, which represents a complete orthonormal set in the $(M,\Gamma ,\lambda )$ space, ${(p_{\lambda }^{\Gamma })}^{†}{p}_{\lambda }^{\Gamma }=1$. Spin and Hamiltonian matrices are transformed accordingly, ${\tilde{S}}^2={(p_{\lambda }^{\Gamma })}^{†}{S}^2{p}_{\lambda }^{\Gamma }$, and $\tilde{H}={(p_{\lambda }^{\Gamma })}^{†}H{p}_{\lambda }^{\Gamma }$, and ${\tilde{P}}_S$ is constructed from ${\tilde{S}}^2$ instead of $S^2.$ Finally, $\tilde{H}$ and ${\tilde{P}}_S$ are applied to a set of random states to set up the generalized eigenvalue problem in $(S,M,\Gamma ,\lambda )$.

As long as $d(S,M,\Gamma ,\lambda )\le 4$ for all S, it may be possible to directly diagonalize $\tilde{H}$, even when $d(M,\Gamma ,\lambda )>4$. On the other hand, if spin adaptation is necessary, an alternative to forming ${\tilde{P}}_S$ is to diagonalize ${\tilde{S}}^2$ and to then transform $\tilde{H}$ into the space with the desired S. This method parallels how Schumann solved the Hubbard model on a square [32]. (Schumann additionally used so-called pseudospin symmetry, where applicable. This symmetry of bipartite Hubbard lattices does not exist in spin-only models. He adapted a basis with definite particle number and spin magnetization first to pseudospin, then to total spin, and lastly to PG symmetry.) However, there is no guarantee that ${\tilde{S}}^2$ can be diagonalized in symbolic form (at least not within a practical time frame), and this has indeed turned out to be impossible in some cases, like $M=0$ and $\Gamma =E_g$ $(d=12)$ in the $s=1$ octahedron. Therefore, the diagonalization of ${\tilde{S}}^2$ is not a universal alternative to using Löwdin’s projector.

We performed symmetry adaptation and diagonalization using a custom-written MATLAB program. Further analyses, such as deriving analytical conditions for phase boundaries as a function of free parameters, which also form the basis for creating the phase diagrams, were carried out in Mathematica.

3. Results

We focused on two rings (symmetric triangle and square) and three polyhedra (tetrahedron, octahedron and cube). Incidentally, except for the cube, these specific systems are trivially integrable within the Heisenberg model (see below), but they require matrix diagonalization when other isotropic terms are included. Tables of wave functions, energies or other properties [24] are usually based on implicit assumptions about which independent parameters are negligible. The systematic construction of all possible terms, whose number rises quickly with N and s, was described in Ref. [19]. For instance, the isotropic Hamiltonian for a group of four $s=\frac{1}{2}$ sites has 9 independent parameters, whereas ten sites permit 8523 parameters [19] (this number, which corresponds to the independent ways of coupling ten rank-1 operators to form a scalar, would be lowered by spatial symmetries), although most of these would be negligibly small in practice. In our analysis, we primarily consider nearest-neighbor (NN) Heisenberg exchange and additionally consider biquadratic exchange or four-center terms. Our tables therefore do not aspire to be useful for analyzing all specific cases but should allow others to effectively verify independent implementations. To apply our method, only the spin and point group symmetry need to be present, and any additional spin Hamiltonian terms that have these symmetries can be just as easily integrated into the framework. All energies are reported in units of the uniform NN coupling constant J, which is chosen to be antiferromagnetic, that is, we set $J=1$.

Triangle. An $s=\frac{1}{2}$ triangle with three different coupling constants is the smallest system that necessitates matrix diagonalization, because there are two $S=\frac{1}{2}$ levels. On the other hand, an isosceles triangle has exchange symmetry, $[\widehat{H},{\widehat{P}}_{12}]=0$; with $J_{13}=J_{23}$, the square of the pair spin ${\widehat{s}}_{12}={\widehat{s}}_1+{\widehat{s}}_2$ is a good quantum number, $[{\widehat{H}}_J,{\widehat{s}}_{12}^2]=0$. This pair spin is then coupled with ${\widehat{s}}_3$ to obtain a total spin multiplet,

$$\begin{array}{l}{\widehat{H}}_J=J_{12}{\widehat{s}}_1·{\widehat{s}}_2+J_{13}({\widehat{s}}_1+{\widehat{s}}_2)·{\widehat{s}}_3=\\ \frac{J_{12}}{2}[{\widehat{s}}_{12}^2-2s(s+1)]+\frac{J_{13}}{2}[{\widehat{S}}^2-{\widehat{s}}_{12}^2-s(s+1)]\end{array},$$

(9)

yielding the spectrum of Equation (10):

$$E_J=\frac{J_{12}}{2}[s_{12}(s_{12}+1)-2s(s+1)]+\frac{J_{13}}{2}[S(S+1)-s_{12}(s_{12}+1)-s(s+1)].$$

(10)

This is the simplest example of Kambe’s method [33]. In the symmetric triangle, all three couplings are equal, thus Equation (9) becomes ${\widehat{H}}_J=\frac{J}{2}[{\widehat{S}}^2-3s(s+1)]$, meaning that all multiplets with the same S are degenerate [1].

However, when including biquadratic exchange in the isosceles triangle $(K_{13}=K_{23}\ne 0)$, ${\widehat{P}}_{12}$ remains a symmetry, $[{\widehat{H}}_K,{\widehat{P}}_{12}]=0$, but ${\widehat{s}}_{12}^2$ is no longer a good quantum number, because $[{({\widehat{s}}_1·{\widehat{s}}_3)}^2+{({\widehat{s}}_2·{\widehat{s}}_3)}^2,{\widehat{s}}_{12}^2]\ne 0$, rendering Kambe’s method inapplicable.

Table 1. Dimensions of combined spin and PG subspaces $(S,\Gamma )$ in symmetric triangles (point group $D_3$) for various values s (first column). Each cell represents the number of multiplets for ascending total spin values S ^a.

s	${\mathbf{A}}_1$	${\mathbf{A}}_2$	E
1/2	0 1 0	0 0 0	1 0 0
1	0 1 0 1	1 0 0 0	0 1 1 0
3/2	0 1 1 0 1	0 1 0 0 0	1 1 1 1 0
2	1 0 1 1 1 0 1	0 1 0 1 0 0 0	0 1 2 1 1 1 0
5/2	0 1 1 1 1 1 0 1	0 1 1 0 1 0 0 0	1 1 2 2 1 1 1 0
3	0 1 0 2 1 1 1 1 0 1	1 0 1 1 1 0 1 0 0 0	0 1 2 2 2 2 1 1 1 0
7/2	0 1 1 1 2 1 1 1 1 0 1	0 1 1 1 1 1 0 1 0 0 0	1 1 2 3 2 2 2 1 1 1 0
4	1 0 1 1 2 1 2 1 1 1 1 0 1	0 1 0 2 1 1 1 1 0 1 0 0 0	0 1 2 2 3 3 2 2 2 1 1 1 0
9/2	0 1 1 1 2 2 1 2 1 1 1 1 0 1	0 1 1 1 2 1 1 1 1 0 1 0 0 0	1 1 2 3 3 3 3 2 2 2 1 1 1 0
5	0 1 0 2 1 2 2 2 1 2 1 1 1 1 0 1	1 0 1 1 2 1 2 1 1 1 1 0 1 0 0 0	0 1 2 2 3 4 3 3 3 2 2 2 1 1 1 0
11/2	0 1 1 1 2 2 2 2 2 1 2 1 1 1 1 0 1	0 1 1 1 2 2 1 2 1 1 1 1 0 1 0 0 0	1 1 2 3 3 4 4 3 3 3 2 2 2 1 1 1 0
6	1 0 1 1 2 1 3 2 2 2 2 1 2 1 1 1 1 0 1	0 1 0 2 1 2 2 2 1 2 1 1 1 1 0 1 0 0 0	0 1 2 2 3 4 4 4 4 3 3 3 2 2 2 1 1 1 0
13/2	0 1 1 1 2 2 2 3 2 2 2 2 1 2 1 1 1 1 0 1	0 1 1 1 2 2 2 2 2 1 2 1 1 1 1 0 1 0 0 0	1 1 2 3 3 4 5 4 4 4 3 3 3 2 2 2 1 1 1 0

^a For example, for $s=\frac{1}{2}$, there are 0, 1 and 0 multiplets with PG label ${\mathrm{A}}_1$ for S = 0, 1 and 2, respectively.

lists the subspace dimensions for symmetric triangles up to $s=\frac{13}{2}$, which is the smallest s that does not permit obtaining the full spectrum in closed form, because there are five $(S=6,\Gamma =\mathrm{E})$ levels. Note that Griffith had already classified terms in triangles, albeit for smaller s [34]. Table 1 shows that up to $s=\frac{3}{2}$, there is at most one level in each $(S,\Gamma )$ sector, so symmetry adaptation suffices to determine the eigenfunctions, and all energies depend linearly on any parameters; phase boundaries in a two- or three-dimensional parameter space would then be straight lines or planes, respectively.

For $s=2$ and $s=\frac{5}{2}$, we exemplarily collect the full spectra as a function of the biquadratic exchange constant K in

Table 2. Full spectrum of an $s=2$ triangle as a function of biquadratic exchange (K), with $J=1$. The second column attributes numbers to all levels for easy comparison with Figure 1. Conditions on K for a specific level to be the ground state are given in the last column.

$(\mathit{S},\mathbf{\Gamma })$	#	Energy	Ground State
(0, A1)	1	$27K-\frac{9}{2}$	$-\frac{1}{30}\le K\le \frac{1}{4}$
(1, A2)	2	$40K-4$	$\varnothing$
(1, E)	3	$31K-4$	$\varnothing$
(2, A1)	4	$72K-3$	$K\le -\frac{1}{30}$
(2, E)	5 6	$\frac{93}{2}K\pm \frac{3}{2}\sqrt{193}\left|K\right|-3$	$\varnothing$ $\varnothing$
(3, A1)	7	$33K-\frac{3}{2}$	$\varnothing$
(3, A2)	8	$15K-\frac{3}{2}$	$K\ge \frac{1}{4}$
(3, E)	9	$51K-\frac{3}{2}$	$\varnothing$
(4, A1)	10	$37K+\frac{1}{2}$	$\varnothing$
(4, E)	11	$19K+\frac{1}{2}$	$\varnothing$
(5, E)	12	$24K+3$	$\varnothing$
(6, A1)	13	$48K+6$	$\varnothing$

and

Table 3. Energy spectrum and ground state conditions for an $s=\frac{5}{2}$ triangle as a function of biquadratic exchange (K), with $J=1$.

$(\mathit{S},\mathbf{\Gamma })$	#	Energy	Ground State
(1/2, E)	1	$\frac{1}{16}(975K-102)$	$-\frac{1}{46}\le K\le \frac{1}{9}$
(3/2, A1)	2	$\frac{1}{16}(867K-90)$	$\frac{1}{9}\le K\le \frac{7}{39}$
(3/2, A2)	3	$\frac{1}{16}(1155K-90)$	$\varnothing$
(3/2, E)	4	$\frac{1}{16}(1443K-90)$	$\varnothing$
(5/2, A1)	5	$\frac{1}{16}(2447K-70)$	$K\le -\frac{1}{46}$
(5/2, A2)	6	$\frac{1}{16}(1295K-70)$	$\varnothing$
(5/2, E)	7 8	$\frac{1}{16}(-70+1679K\pm 24\sqrt{19}|K|)$	$\varnothing$ $\varnothing$
(7/2, A1)	9	$\frac{1}{16}(1971K-42)$	$\varnothing$
(7/2, E)	10 11	$\frac{3}{16}(-14+401K\pm 32\sqrt{34}|K|)$	$\varnothing$ $\varnothing$
(9/2, A1)	12	$\frac{1}{16}(783K-6)$	$\varnothing$
(9/2, A2)	13	$\frac{1}{16}(399K-6)$	$K\ge \frac{7}{39}$
(9/2, E)	14	$\frac{1}{16}(1359K-6)$	$\varnothing$
(11/2, A1)	15	$\frac{1}{16}(1091K+38)$	$\varnothing$
(11/2, E)	16	$\frac{1}{16}(611K+38)$	$\varnothing$
(13/2, E)	17	$\frac{1}{16}(975K+90)$	$\varnothing$
(15/2, A1)	18	$\frac{1}{16}(1875K+150)$	$\varnothing$

, respectively. (Our Hamiltonian is not exhaustive, e.g., three-center terms, which would occur for $s>\frac{1}{2}$, are ignored.)

The ground states as a function of K and magnetic field strength B (Zeeman term, ${\widehat{H}}_B=B{\widehat{S}}_z$) are plotted in Figure 1 $(s=2)$ and Figure 2 $(s=2)$. These phase diagrams were built on simple analytical conditions, which resemble the ground state conditions in the last columns of Table 2 and Table 3, respectively, but also consider the magnetic field strength B as another independent parameter.

Square. The Heisenberg square, with sites numbered consecutively, is again integrable by Kambe’s method [1,35], because $[{\widehat{H}}_J,{\widehat{s}}_{13}^2]=[{\widehat{H}}_J,{\widehat{s}}_{24}^2]=0$, but ${\widehat{s}}_{13}^2$ or ${\widehat{s}}_{24}^2$ do not commute with the biquadratic ${\widehat{H}}_K$. In addition, there exist two independent four-center interactions in the $s=\frac{1}{2}$ square (For $s>\frac{1}{2}$, additional four-center (and three-center) interactions involve local rank-2 operators): the frequently discussed cyclic exchange,

$${\widehat{H}}_C=C\left[({\widehat{s}}_1·{\widehat{s}}_2)({\widehat{s}}_3·{\widehat{s}}_4)+({\widehat{s}}_1·{\widehat{s}}_4)({\widehat{s}}_2·{\widehat{s}}_3)-({\widehat{s}}_1·{\widehat{s}}_3)({\widehat{s}}_2·{\widehat{s}}_4)\right],$$

(11)

and the less common non-cyclic exchange [19],

$${\widehat{H}}_{\tilde{C}}=\tilde{C}\left[({\widehat{s}}_1·{\widehat{s}}_2)({\widehat{s}}_3·{\widehat{s}}_4)+({\widehat{s}}_1·{\widehat{s}}_4)({\widehat{s}}_2·{\widehat{s}}_3)+6({\widehat{s}}_1·{\widehat{s}}_3)({\widehat{s}}_2·{\widehat{s}}_4)\right].$$

(12)

Except for $s=\frac{1}{2}$, ${\widehat{H}}_C$ and ${\widehat{H}}_{\tilde{C}}$ do not commute with ${\widehat{s}}_{13}^2$ and ${\widehat{s}}_{24}^2$. In other words, biquadratic exchange and multi-center interactions generally prevent trivial spin-coupling solutions. Dimensions of the $(S,\Gamma )$ sectors are listed in

Table 4. Dimensions of the $(S,\Gamma )$ spaces in the square (point group $D_4$). For further information, see caption and footnote to Table 1.

s	${\mathbf{A}}_1$	${\mathbf{A}}_2$	${\mathbf{B}}_1$	${\mathbf{B}}_2$	E
1	2 0 2 0 1	0 1 0 0 0	1 0 1 0 0	0 1 1 1 0	0 2 1 1 0
3/2	2 0 3 1 2 0 1	0 1 1 1 0 0 0	2 0 2 0 1 0 0	0 2 1 2 1 1 0	0 3 2 3 1 1 0
2	3 0 4 1 4 1 2 0 1	0 2 1 2 1 1 0 0 0	2 0 3 1 2 0 1 0 0	0 2 2 3 2 2 1 1 0	0 4 3 5 3 3 1 1 0
5/2	3 0 5 2 5 2 4 1 2 0 1	0 2 2 3 2 2 1 1 0 0 0	3 0 4 1 4 1 2 0 1 0 0	0 3 2 4 3 4 2 2 1 1 0	0 5 4 7 5 6 3 3 1 1 0

, and the spectrum for $s=1$ as a function of biquadratic and cyclic exchange (K and C, respectively) is given in

Table 5. Spectrum for an $s=1$ square as a function of biquadratic (K) and four-center coupling (C), at $J=1$. The ground state conditions with K set to zero are defined in the last column.

$(\mathit{S},\mathbf{\Gamma })$	#	Energy	Ground State $(\mathit{K}=0)$	Ground State $(\mathit{C}=0)$
(0, A1)	1 2	$8K+\frac{3}{2}C-\frac{3}{2}\pm \frac{1}{2}\sqrt{64K^2+96KC-32K+41C^2-34C+9}$	$\varnothing$ $C\le \frac{1}{4}$	$\varnothing$ $K\le \frac{1}{2}$
(0, B1)	3	$8K+3C-1$	$\varnothing$	$\varnothing$
(1, A2)	4	$5K-2C-\frac{1}{2}$	$\varnothing$	$\varnothing$
(1, B2)	5	$9K+2C-\frac{5}{2}$	$C=\frac{1}{4}$	$\varnothing$
(1, E)	6 7	$\frac{11}{2}K+C-\frac{3}{4}\pm \frac{1}{2}\sqrt{9K^2+20KC-K+20C^2-10C+\frac{9}{4}}$	$\varnothing$ $\varnothing$	$\varnothing$ $K=\frac{1}{2}$
(2, A1)	8 9	$\frac{13}{2}K-\frac{3}{4}\pm \frac{1}{2}\sqrt{25K^2-48KC+K+32C^2+8C+\frac{9}{4}}$	$\varnothing$ $\frac{1}{4}\le C\le \frac{13}{2}$	$\varnothing$ $K\ge \frac{1}{2}$
(2, B1)	10	$5K+\frac{1}{2}$	$\varnothing$	$\varnothing$
(2, B2)	11	$4K$	$\varnothing$	$\varnothing$
(2, E)	12	$7K-2C-\frac{1}{2}$	$\varnothing$	$\varnothing$
(3, B2)	13	$4K-3C$	$C\ge \frac{13}{2}$	$\varnothing$
(3, E)	14	$4K+C+1$	$\varnothing$	$\varnothing$
(4, A1)	15	$4K+C+2$	$\varnothing$	$\varnothing$

.

The (K, C) ground state criteria are rather complicated and thus not shown here. However, one clear observation, which could not be obtained directly from numerical computations, is that, under the assumption of antiferromagnetic Heisenberg coupling, $J=1$, none of the levels 1, 3, 4, 6, 8, 10, 11, 12, 14 or 15 are the ground states for any set (K, C). Phase diagrams including a magnetic field, setting either $K=0$ or $C=0$, are shown in Figure 3.

Polyhedra. Here we shall summarily discuss a few highly symmetric polyhedra. In the tetrahedron, all possible pairs are coupled equally, like in the triangle. Therefore, ${\widehat{H}}_J$ is proportional to ${\widehat{S}}^2$, ${\widehat{H}}_J=\frac{J}{2}[{\widehat{S}}^2-4s(s+1)]$. Biquadratic exchange again lifts degeneracies because it does not commute with pair spins, $[{\widehat{H}}_K,{\widehat{s}}_{12}^2]\ne 0$, etc. For $s=\frac{1}{2}$, where a single four-center term is compatible with $T_d$ symmetry,

$${\widehat{H}}_T=T[({\widehat{s}}_1·{\widehat{s}}_2)({\widehat{s}}_3·{\widehat{s}}_4)+({\widehat{s}}_1·{\widehat{s}}_4)({\widehat{s}}_2·{\widehat{s}}_3)+({\widehat{s}}_1·{\widehat{s}}_3)({\widehat{s}}_2·{\widehat{s}}_4)],$$

(13)

the pair spins commute with ${\widehat{H}}_T$, but, as in the case of the square, the respective commutators are non-zero for $s>\frac{1}{2}$. The dimensions of symmetry subspaces are shown in

Table 6. Dimensions of the $(S,\Gamma )$ spaces in the tetrahedron (point group $T_d$).

s	${\mathbf{A}}_1$	${\mathbf{A}}_2$	E	${\mathbf{T}}_1$	${\mathbf{T}}_2$
1/2	0 0 1	0 0 0	1 0 0	0 0 0	0 1 0
1	1 0 1 0 1	0 0 0 0 0	1 0 1 0 0	0 1 0 0 0	0 1 1 1 0
3/2	1 0 1 1 1 0 1	1 0 0 0 0 0 0	1 0 2 0 1 0 0	0 1 1 1 0 0 0	0 2 1 2 1 1 0
2	1 0 2 0 2 1 1 0 1	0 0 1 0 0 0 0 0 0	2 0 2 1 2 0 1 0 0	0 2 1 2 1 1 0 0 0	0 2 2 3 2 2 1 1 0
5/2	1 0 2 1 2 1 2 1 1 0 1	1 0 1 0 1 0 0 0 0 0 0	2 0 3 1 3 1 2 0 1 0 0	0 2 2 3 2 2 1 1 0 0 0	0 3 2 4 3 4 2 2 1 1 0
3	2 0 2 1 3 1 3 1 2 1 1 0 1	1 0 1 1 1 0 1 0 0 0 0 0 0	2 0 4 1 4 2 3 1 2 0 1 0 0	0 3 2 4 3 4 2 2 1 1 0 0 0	0 3 3 5 4 5 4 4 2 2 1 1 0

. The dimensions for $s=2$ were previously reported in Table 4 of Ref. [36].

The complete analytical spectrum for the $s=\frac{3}{2}$ tetrahedron as a function of K and T is exemplarily provided in

Table 7. Spectrum of the $s=\frac{3}{2}$ tetrahedron $(J=1)$ as a function of biquadratic and four-center coupling, K and T.

$(\mathit{S},\mathbf{\Gamma })$	#	Energy
(0, A1)	1	$\frac{327}{8}K+\frac{327}{16}T-\frac{15}{4}$
(0, A2)	2	$\frac{135}{8}K+\frac{135}{16}T-\frac{15}{4}$
(0, E)	3	$\frac{423}{8}K+\frac{423}{16}T-\frac{15}{4}$
(1, T1)	4	$\frac{179}{8}K+\frac{99}{16}T-\frac{13}{4}$
(1, T2)	5 6	$\frac{299}{8}K+\frac{227}{16}T-\frac{13}{4}\pm \sqrt{81K^2+120KT+46T^2}$
(2, A1)	7	$\frac{315}{8}K+\frac{267}{16}T-\frac{9}{4}$
(2, E)	8 9	$\frac{243}{8}K+\frac{51}{16}T-\frac{9}{4}\pm \frac{3}{2}\sqrt{36K^2+36KT+37T^2}$
(2, T1)	10	$\frac{219}{8}K+\frac{27}{16}T-\frac{9}{4}$
(2, T2)	11	$\frac{219}{8}K-\frac{69}{16}T-\frac{9}{4}$
(3, A1)	12	$\frac{303}{8}K-\frac{297}{16}T-\frac{3}{4}$
(3, T1)	13	$\frac{159}{8}K-\frac{81}{16}T-\frac{3}{4}$
(3, T2)	14 15	$\frac{279}{8}K-\frac{93}{16}T-\frac{3}{4}\pm \frac{3}{4}\sqrt{144K^2-360KT+289T^2}$
(4, A1)	16	$\frac{287}{8}K-\frac{153}{16}T+\frac{5}{4}$
(4, E)	17	$\frac{143}{8}K+\frac{63}{16}T+\frac{5}{4}$
(4, T2)	18	$\frac{191}{8}K-\frac{153}{16}T+\frac{5}{4}$
(5, T2)	19	$\frac{171}{8}K+\frac{27}{16}T+\frac{15}{4}$
(6, A1)	20	$\frac{243}{8}K+\frac{243}{16}T+\frac{27}{4}$

.

For most of the levels, the (K, T) criteria for a level to be the lowest-energy state are again quite complex. For $J=1$, none of the levels 4, 5, 7, 8, 10, 11, 14, 16, 17, 19 or 20 can be the ground states for any (K, T), and 1 is a ground state only on the line $T=-2$, $K\ge -\frac{1}{28}$. Figure 4 shows phase diagrams including a magnetic field, setting $K=0$ or $T=0$.

For the octahedron, a formulation of ${\widehat{H}}_J$ as in Equation (14),

$${\widehat{H}}_J=\frac{J}{2}({\widehat{S}}^2-{\widehat{s}}_{14}^2-{\widehat{s}}_{25}^2-{\widehat{s}}_{36}^2),$$

(14)

shows that eigenstates have definite pair spins for diametrically opposite sites. Accidental degeneracies between terms belonging to different $(S,\Gamma )$ sectors are lifted when incorporating multi-center terms [19] or biquadratic exchange.

Table 8. Dimensions of the $(S,\Gamma )$ spaces in the octahedron (point group $O_h$).

s	${\mathbf{A}}_{1\mathbf{g}}$	${\mathbf{A}}_{2\mathbf{g}}$	${\mathbf{A}}_{1\mathbf{u}}$	${\mathbf{A}}_{2\mathbf{u}}$	${\mathbf{E}}_{\mathbf{g}}$	${\mathbf{E}}_{\mathbf{u}}$	${\mathbf{T}}_{1\mathbf{g}}$	${\mathbf{T}}_{2\mathbf{g}}$	${\mathbf{T}}_{1\mathbf{u}}$	${\mathbf{T}}_{2\mathbf{u}}$
1/2	0 1 0 1	1 0 0 0	0 0 0 0	1 0 0 0	0 1 1 0	0 0 0 0	0 0 0 0	0 1 0 0	1 0 1 0	0 1 0 0
1	3 0 3 1 2 0 1	0 2 0 2 0 0 0	1 0 0 0 0 0 0	0 1 0 1 0 0 0	1 2 4 2 2 1 0	0 1 1 0 0 0 0	0 2 1 1 0 0 0	2 1 3 1 1 0 0	0 4 2 3 1 1 0	1 2 3 2 1 0 0
3/2	0 5 2 7 3 4 2 2 0 1	3 1 4 3 4 1 2 0 0 0	0 2 0 2 0 0 0 0 0 0	3 0 3 1 2 0 1 0 0 0	1 6 8 8 7 6 3 2 1 0	1 2 4 2 2 1 0 0 0 0	2 4 6 5 4 2 1 0 0 0	0 7 5 8 4 4 1 1 0 0	5 4 11 7 9 4 4 1 1 0	1 7 7 9 6 5 2 1 0 0

shows that the $s=\frac{1}{2}$ octahedron is solved directly by spin and PG adaptation, and that the full spectrum can be obtained in closed form for $s=1$ (not detailed here, due to excessively long expressions for solutions in three- and four-dimensional spaces).

Finally, the Heisenberg cube cannot be solved by Kambe’s method, as it lacks conserved subsystem spins. The dimensions of the symmetry spaces for $s=\frac{1}{2}$, the only s value for which the cube is fully solvable, are collected in

Table 9. Dimensions of the $(S,\Gamma )$ spaces in the $s=\frac{1}{2}$ cube (point group $O_h$).

${\mathbf{A}}_{1\mathbf{g}}$	${\mathbf{A}}_{2\mathbf{g}}$	${\mathbf{A}}_{1\mathbf{u}}$	${\mathbf{A}}_{2\mathbf{u}}$	${\mathbf{E}}_{\mathbf{g}}$	${\mathbf{E}}_{\mathbf{u}}$	${\mathbf{T}}_{1\mathbf{g}}$	${\mathbf{T}}_{2\mathbf{g}}$	${\mathbf{T}}_{1\mathbf{u}}$	${\mathbf{T}}_{2\mathbf{u}}$
3 0 2 0 1	0 0 0 0 0	1 0 0 0 0	0 2 0 1 0	2 0 2 0 0	0 1 1 0 0	0 2 0 0 0	1 2 2 1 0	0 3 1 1 0	1 1 1 0 0

. As far as we know, the complete spectrum of ${\widehat{H}}_J$ was not derived previously, and we therefore present it in

Table 10. Spectrum of the $s=\frac{1}{2}$ Heisenberg cube.

$(\mathit{S},\mathbf{\Gamma })$	#	Energy
(0, A1g)	1 2 3	$\begin{array}{l}-\frac{5}{3}+\frac{2\sqrt{7}}{3}\cos \alpha \pm \frac{2\sqrt{21}}{3}\sin \alpha \\ -\frac{5}{3}-\frac{4\sqrt{7}}{3}\cos \alpha \end{array}$
(0, A1u)	4	0
(0, Eg)	5 6	$-1\pm \sqrt{2}$
(0, T2g)	7	−1
(0, T2u)	8	−2
(1, A2u)	9 10	−4 0
(1, Eu)	11	−1
(1, T1g)	12 13	$-\frac{1}{2}(1\pm \sqrt{5})$
(1, T2g)	14 15	$-\frac{1}{2}(3\pm \sqrt{5})$
(1, T1u)	16 17 18	$\begin{array}{l}-\frac{2}{3}-\frac{\sqrt{10}}{3}\cos \beta \pm \frac{\sqrt{30}}{3}\sin \beta \\ -\frac{2}{3}+\frac{2\sqrt{10}}{3}\cos \beta \end{array}$
(1, T2u)	19	0
(2, A1g)	20 21	$-1\pm \sqrt{2}$
(2, Eg)	22 23	$\frac{1}{2}(1\pm \sqrt{5})$
(2, Eu)	24	0
(2, T2g)	25 26	$\pm 1$
(2, T1u)	27	−1
(2, T2u)	28	1
(3, A2u)	29	0
(3, T2g)	30	1
(3, T1u)	31	2
(4, A1g)	32	3

. We note accidental degeneracies for $E=-1,0,+1,-1\pm \sqrt{2}$.

$$\alpha =\frac{1}{3}{\tan }^{-1}\left(\frac{3\sqrt{591}}{13}\right), \beta =\frac{1}{3}{\tan }^{-1}\left(3\sqrt{111}\right)$$

4. Summary and Conclusions

While numerical methods are commonly used to investigate spin models of magnetic solids and molecules, small clusters with spatial symmetry allow for analytical diagonalization of the Hamiltonian, and a symbolic representation of spectra or wave functions may provide deeper insights that go beyond mere numerical data. To factor the Hamiltonian into symmetry subspaces and thus enable analytical diagonalization, we provided a simple yet effective approach to adapt the basis to both total spin (using Löwdin’s projector) and point group (PG) symmetry. The construction of PG projectors was extensively discussed. Overall, our procedure for employing spin and PG symmetry to set up a generalized eigenvalue problem in a subspace is not intended to be computationally optimal but designed to be easily followed and implemented. We chose small rings and polyhedra as examples and elucidated how additional interactions (beyond the Heisenberg model) prevent trivial integrability. Our aim was not to compile exhaustive tables but rather to highlight specific results that may be useful for verifying independent implementations of the analytical diagonalization scheme.

It is worth noting that a similar procedure would also be applicable to anisotropic systems. Although a general anisotropic Hamiltonian does not conserve spin, $[\widehat{H},{\widehat{S}}^2]\ne 0$ and $[\widehat{H},{\widehat{S}}_z]\ne 0$, one can still make use of PG symmetry and apply respective projectors to random states. With anisotropy, most site permutations must be combined with spin rotations to represent symmetries [24,26], and this necessitates working with the respective double group for systems with half-integer spin. Analytical solutions for anisotropic models are more severely restricted in terms of system size, because the group is smaller (lacking spin symmetry). Lastly, the present method could also be adapted for use with the Hubbard model, whose additional pseudospin symmetry on a bipartite lattice allows to further block-diagonalize the Hamiltonian [32,37]. However, the Hubbard model has itinerant electrons, and, at half-filling, its state space is larger than that of the $s=\frac{1}{2}$ Heisenberg model with the same number of sites. Therefore, the limits on system size are stricter.

Note that we leveraged symmetry to factor H into subspaces for easier diagonalization, but a classification of multiplets in terms of S and $\Gamma$ holds qualitative value too, e.g., for deducing spectroscopic selection rules [23,38,39] or for assessing the momentum-transfer dependencies of inelastic neutron scattering intensities [23,39,40]. Symmetry classifications also aid in analyzing the mixing of multiplets by anisotropic terms [26,38].

Analytical solutions can be directly transferred when adding a further spin that has equal couplings to all other sites. The interaction of the original system with such a central site can be handled with Kambe’s method, as detailed in Ref. [8]. However, the number of non-trivial systems that allow for the closed-form solutions of their entire spectrum is naturally limited by the requirement that neither subspace dimension exceeds four. In some situations, analytical diagonalization could still be used in the smaller spaces of large S values, which are important to consider for ferromagnetic coupling or in high magnetic fields.