Permutation Symmetry in Spin Adapted Many-Body Wave Functions

In the domain of exchange-coupled PNTM clusters, local emergent symmetries exist which can be exploited to greatly increase the sparsity of the CI eigensolutions of such systems. Sparsity of the CI secular problem is revealed by exploring the site permutation space within spin-adapted many-body bases, and highly compressed wave functions may arise by finding optimal site orderings. However, the factorial cost of searching through the permutation space remains a bottleneck for clusters with a large number of metal centers. In this work, we explore ways to reduce the factorial scaling, by combining permutation and point group symmetry arguments, and using commutation relations between cumulative partial spin and the Hamiltonian operators, (cid:20)(cid:16) ˆ S ( n ) (cid:17) 2 , ˆ H (cid:21) . Certain site orderings lead to commuting operators, from which more sparse wave functions arise. Two graphical strategies will be discussed, one to rapidly evaluate the commutators of interest, and one in the form of a tree search algorithm to predict how many and which distinct site permutations are to be analyzed, eliminating redundancies in the permutation space. Particularly interesting is the case of the singlet spin states for which an additional reversal symmetry can be utilized to further reduce the distinct site permutations.

In the domain of exchange-coupled PNTM clusters, local emergent symmetries exist which can be exploited to greatly increase the sparsity of the CI eigensolutions of such systems.Sparsity of the CI secular problem is revealed by exploring the site permutation space within spin-adapted manybody bases, and highly compressed wave functions may arise by finding optimal site orderings.However, the factorial cost of searching through the permutation space remains a bottleneck for clusters with a large number of metal centers.In this work, we explore ways to reduce the factorial scaling, by combining permutation and point group symmetry arguments, and using commutation relations between cumulative partial spin and the Hamiltonian operators, Ŝ(n) 2 , Ĥ .Certain site orderings lead to commuting operators, from which more sparse wave functions arise.Two graphical strategies will be discussed, one to rapidly evaluate the commutators of interest, and one in the form of a tree search algorithm to predict how many and which distinct site permutations are to be analyzed, eliminating redundancies in the permutation space.Particularly interesting is the case of the singlet spin states for which an additional reversal symmetry can be utilized to further reduce the distinct site permutations.
Understanding the ground-and excited-state electronic structures of polynuclear transition metal (PNTM) clusters and their magnetic and catalytic properties represent important challenges in modern electronic structure theory. 1 These compounds, abundant in nature, are the active sites of metalloproteins and make key enzymatic reactions possible.2][13][14][15][16] Biomimetic PNTM clusters also exist, [17][18][19][20] such as the Co(II) 3 Er(III)(OR) 4 cubane relevant in the context of the artificial water splitting process. 21,22NTM clusters consist of transition metal ions (lanthanide or actinide ions are also possible), at variable oxidation states, whose valence d (or f) orbitals are partially occupied.In general, the unpaired electrons at each magnetic site obey the first Hund's rule and feature parallel spins.However, strong ligand field effects may cause some of the valence electrons to pair, effectively reducing the formal local spin.[25] The local spins couple across the magnetic sites, leading to a plethora of energetically low electronic states. 26In this context, electronic states can be distinguished into collinear states, with spins across the sites parallel or anti-parallel aligned, and noncollinear states, with spin across magnetic centers interacting at an angle (in the simple angular momentum vector model).
In the context of exchange-coupled PNTM complexes, including Fe(III) 4 S 4 and Co(II) 3 Er(III)(OR) 4 , algorithms have been advised within DMRG to unitarily transform the active molecular orbitals (MO) into one-electron bases that are more conducive to fast convergence with respect to the M value. 176,177Similarly, within the spin-adapted FCIQMC algorithm based on the graphical unitary group approach (GUGA), 118,[178][179][180][181][182][183][184] MO transformations have been identified, that increase the sparsity of the CI Hamiltonian and the corresponding ground-and excited-state eigenvectors, dramatically favoring the convergence of FCIQMC with respect to the walker numbers. 26,86,118,119,170,185,186It has numerically been shown that MO transformations that are opti-mal for DMRG differ from the ones that are optimal for GUGA-FCIQMC. 185,186In the case of FCIQMC, large active space calculations, correlating up to 44 electrons in 32 orbitals for the Fe(III) 4 S 4 cubane, CAS(44,32), 119,170 and up to 56 electrons in 56 orbitals for a Co(II) 3 Er(III)(OR) 4 cubane, CAS(56,56), 186 have been performed at a relatively lower walker number for specific orbital orderings.
Orbital (site) orderings can be identified that lead to sparse CI Hamiltonian matrices with a unique (quasi) block diagonal structure, when expressed in spin-adapted bases.Precisely this block diagonal structure minimizes the mixing of CSFs, from which the sparse wave functions arise.Also, in virtue of the blocking it has been possible to selectively target excited states. 170sing general model Heisenberg Hamiltonian operators [187][188][189][190][191][192] for clusters of various dimensions, shapes, and local spin value at each magnetic site (in general S local > 1/2), Li Manni and coworkers 193 have identified the rationale behind the wave function sparsity and the blocking of the CI Hamiltonian matrix in the commutation relations between cumulative partial spin operators, where M runs over the first n magnetic sites, and the Hamiltonian operator, Ĥ .The choice of the Heisenberg model Hamiltonian is two-fold.First, Heisenberg models describe the magnetic interactions in exchange coupled PNTM clusters with an high degree of accuracy, and any valuable observation in terms of sparsity can be promptly transferred from the Heisenberg model to the full ab initio Hamiltonian. 26,118,170,185,193Moreover, site reordering is relevant only for open-shell electronic configurations, and its effects in terms of sparsity become more important as the number of unpaired electron (seniority) increases; from this point of view the Heisenberg model is most suited as the CI space consists exclusively of unpaired electrons.Specific site orderings may lead to vanishing commutators Commuting operators are associated with compatible observables, which can be simultaneously determined.Furthermore, considering that cumulatively spin adapted bases, such as the GUGA configuration state functions (CSFs), are by construction eigenbases of cumulative partial spin operators, it follows that CSFs with different expectation value of partial spin operators may lead (when Eq. 2 is verified) to vanishing Hamiltonian matrix elements, causing the corresponding Hamiltonian matrix to exhibits a unique block diagonal structure (vide infra).How can we identify a site ordering that fulfills Eq. 2? Is such site ordering unique?Considering that the smallest relevant partial spin operator is the one of a single site, Ŝ(1) 2 , and that a CSF can be eigenfunction of such operator only if orbitals (and their electrons) residing on one magnetic center are adjacent, the grouping of orbitals on each site will always be assumed in the remainder of this document.Therefore, the permutation space to be investigated scales as N!, where N is the number of magnetic centers, and the search for optimal site orderings becomes rapidly computational demanding for clusters with an increasingly larger number of magnetic sites.Permutation and point group symmetries can be exploited to substantially reduce the search for site orderings that maximize wave function sparsity.This represents the main focus of the present work.
Site ordering not only affect the sparsity of CI Hamiltonian matrix of exchange coupled systems; in recent literature, its role has been discussed in the context of the Jordan-Wigner transformation for converting systems of spins into systems of fermions and vice versa; and in that context an extended Jordan-Wigner transformation has been proposed to circumvent this dependency. 194,195Commutators Between Partial Spin and Hamiltonian Operators The connection between site ordering, partial spin operators, commutators and blocking of the CI Hamiltonian matrix in spin adapted bases is quite remarkable, and it emerges from the wellknown cumulative spin couplings of the GUGA method, which we briefly recall below.
As in other genealogically constructed spin eigenfunctions, 196 in GUGA, the spin of each individual electron is coupled to all previous ones in a cumulative manner; therefore, each CSF is by construction an eigensolution of any cumulative partial spin operator (not necessarily over magnetic sites).In the context of exchange-coupled PNTM clusters, it is reasonable to classify CSFs on the basis of their expectation values over the cumulative partial spin operators of the magnetic sites.The following CSFs |uuu, uud, ddd⟩ , |uuu, udu, ddd⟩ , |uuu, udd, udd⟩ , If we consider two CSFs with different eigenvalues (a and if Ŝ(n) 2 , Ĥ = 0, we may write and since a ̸ = b, it must be deduced that ⟨CSF| H |CSF ′ ⟩ = 0, from which the block diagonal structure of the Hamiltonian matrix arises.In the example above, this well known property of commuting operators 197 would imply that CSFs with different ŜAB 2 expectation values will not couple over Ĥ , and the CI Hamiltonian matrix will exhibit a blocked structure over the ŜAB 2 expectation values.Thus, it is pivotal to establish commutation relations between cumulative partial spin and the Hamiltonian operators for different site orderings to understand and gain control over the structure and sparsity of the Hamiltonian matrix.
Considering that often magnetic interactions in exchangecoupled PNTM clusters are well described by Heisenberg models, and that for the Heisenberg model the evaluation of the above mentioned commutators is relatively easy, we will describe commutation relations and blocking on the basis of these models; and, in virtue of the correspondence between the two Hamiltonian operators, any insight gained over the simpler model can be transferred in first approximation to the more complex ab initio Hamiltonian.The strategy of transferring information from a simple model to a more complex one in the context of wave function compression has been already utilized in previous works. 118,170,193It is to be noted that in general exchange coupled PNTM clusters exhibit local spin S local > 1/2, and the analysis below will deal with this general case, rather specific examples with, say, S = 1/2 Heisenberg models.
Although Heisenberg models only approximately describe the full ab initio Hamiltonians of molecules, they are useful to understand the effect of orbital orderings on wave function compression.This is attributed to the fact that only u and d couplings are involved in wave function compression.Empty and doubly occupied orbitals are invariant under orbital orderings, This implies that CSFs with high seniority contribute more than low seniority CSFs to the wave function compression effects, as shown in Figure SI.1.As an edge case, zero seniority CSFs are merely single Slater determinants, that are not affected by the ordering in which orbitals are coupled.
Thus, rather than focusing on ab initio Hamiltonians, we simplify the problem by restricting our analysis to Heisenberg models, which still hold the essential features of wave function compression via orbital orderings.As illustrated with a Co(II) 3 Er(III)(OR) 4 cubane cluster in Section 5, our analysis on Heisenberg models can be easily extended to ab initio systems.
It is important to highlight that already Shavitt compared CI expansions based on different orbital transformations. 198Unfortunately, the orbital transformations he suggested and explored and the chemical systems investigated were not conducive to any major advantage in terms of compression, or reduced coupling (sparsity), and no other work followed.On the contrary, in the case of exchange coupled PNTM clusters the compression effects that follow localizations (in case of ab initio models, based on molecular orbitals) and site reorderings, are dramatic, with many orders of magnitude reduction of the non-vanishing coefficients in the CI expansion.This wave function compression is particularly appealing for methods that take advantage of the sparsity, including FCIQMC, and thus worthy investigating.
In the following we summarize two approaches to calculate the , Ĥ commutators, the conventional algebraic approach based on the Levi-Civita symbols, and a faster graphical approach.

Levi-Civita Approach
The commutator between a generalized Heisenberg Hamiltonian and a cumulative partial spin operator over the first consecutive n sites (Eq. 1) can be written as a sum of commutators between the scalar products of local spin operators, multiplied by the corresponding magnetic coupling constants, J PQ , The scalar products may further be expanded into their Cartesian components (i, j ∈ {x, y, z}), The rank reduction relation 56 may be applied to each component of Eq. 6.
Considering commutation relations of local spin operators, Ŝi P , Ŝ j Q = iδ PQ ε i jk Ŝk P , from angular momentum theory 199 , Eq. 7 suggests that non-vanishing commutators are of the form which contain three spin operators over different sites (i.e., , and the forth one, ŜX , is one of the sites on the left of the commutator, X = P or Q.Thus, ŜA • ŜA , ŜB • ŜC = 0, as the spin operator appearing twice, ŜA , only appears on the left side of the commutator.On the contrary, ŜA • ŜB , ŜA • ŜC ̸ = 0.
For convenience, for the non-vanishing commutators we introduce the ternary quantity TABC , where the prefactor −2 is conveniently chosen to absorb the overall sign in Eq. 4 as well as the prefactor of the cross product ŜA • ŜB that arises in In Eq. 9, we have employed the Levi-Civita symbol, ε i jk , which exposes permutations in a way compatible with tensor analysis.We also notice that the values of non-vanishing commutators are compatible with the definition of scalar triple product.By definition, TABC is zero unless there are three different subscripts, (ABC).Changing the order of the subscripts results in a factor of (−1) p , where p is the number of pairwise swaps of the subscripts to recover the original (ABC) order.For example, TABC = TBCA = − TBAC .
We consider the Heisenberg Hamiltonian operator of an isosceles triangle as an example where the J BC = J AC equality applies.The commutation relation is promptly recognized, as its components are of the type ŜX • ŜX , ŜP • ŜQ = 0 (two identical spin operators on the same side of the commutator).If the ABC site ordering is chosen, each CSF is an eigenfunction of ŜA + ŜB 2 .The commutator • ŜB and that Eq. 12 applies.Thus, where in the last step we have used the equality T BAC = −T ABC .
For the commutator vanishes, the corresponding Hamiltonian matrix is certainly blocked on the basis of the partial spin eigenbasis and the eigensolutions of Ĥ will be compressed.
In the ACB ordering all CSFs will be eigenfunctions of the operator ŜA + ŜC 2 .The following commutator would be of interest While in Eq. 13 the commutator vanishes, leading to the advantageous blocking of the corresponding Hamiltonian matrix, in Eq. 14 the commutator in general does not vanish, as J BC ̸ = J AB .Thus, in the latter case a non-block diagonal structure of the Hamiltonian matrix is to be expected under the ŜA + ŜC 2 metric.
As shown in Fig. 6 of Reference 193, the sparsity monotonically increases as J BC approaches the J AB value.

Graphical Approach
The observations from Section 1.1 enable a convenient graphical strategy to evaluate commutators between cumulative partial spin and the Heisenberg Hamiltonian operators, Ŝ(n) 2 , Ĥ , by relying on the triangle rules summarized in Algorithm 1.
For a given ordering of subscripts, say IJK in Eq. 15, the coupling constant, J JK , between the site K in G 2 and the site to its left in the subscripts (J in the case of IJK subscripts) has a positive sign, while the coupling constant J IK with the other site, I, has a negative sign.
From the two sets, two triangles can be formed, namely (132)  and (134), which are marked in green and blue in Figure 1a, respectively.The two triangles identify the entries in the summation in Eq. 15 (see Algorithm 2).Similarly, in calculating Ŝ1 + Ŝ2 + Ŝ3 2 , Ĥ , the two groups G 1 = {1, 2, 3} and G 2 = {4} are identified, from which three triangles can be drawn, namely (124), ( 234) and (134) (different colors in Figure 1b), corresponding to the three summands in Algorithm 3. Eq. 15 shows that a source of (quasi) vanishing commutators is the cancellation of pairs of magnetic coupling constants with the same (or close) value, when an appropriate orbital/site or-dering is chosen; thus, generally wave function compression can be achieved both for ferromagnetic and antiferromagnetic interactions, as long as a orbital/site ordering exists that leads to such cancellation.
2 Permutation Symmetry within the Hamiltonian While in the previous section a graphical procedure has been discussed to access the commutators between cumulative partial spin and the Hamiltonian operators, in this section the reduction of the dimensionality of the permutation space is discussed on the basis of the internal symmetries of the Hamiltonian operator.Two forms of Heisenberg Hamiltonian operators will be utilized, one in which the Hamiltonian operator commutes with all cumulative partial spins operator (Section 2.1), and a less strict one that commutes only with the first m cumulative spins (Section 2.2).

Hamiltonian Commuting with All Cumulative Partial Spins
A general expression for an Heisenberg Hamiltonian that commutes with all cumulative partial spins operators follows The Heisenberg Hamiltonian operator of the isosceles triangle and of the 4-site pyramidal structure (see Fig. 2) feature this form of Hamiltonian operator.Circumsphere center Fig. 2 Pyramidal structure that fulfills Eq. 16.Notice that the base of the pyramid is an isosceles triangle and the apex (vertex 4) is equally spaced from the vertices of the base (1, 2 and 3).In the middle of the pyramid the circumsphere center is also marked for a 5-site cluster that fulfills Eq. 16.
In order to obtain an isosceles triangle, given the first two sites, 1 and 2, the third site, 3, must reside in the plane dividing the base in two equal parts and perpendicular to the 12 direction (a two-dimensional degree of freedom).The 4-site polyhedron that correspond to the Hamiltonian operator of Eq. 16 is obtained by adding to the isosceles triangle a forth point residing on the line passing from the triangle's circumcenter and perpendicular to the triangle plane (a one-dimensional degree of freedom).On the basis of geometrical criteria, a 5-site system featuring an Hamiltonian as given by Eq. 16 is also possible; starting from the 4site polyhedron of Figure 2 the fifth site can only be uniquely located at the center of the circumsphere (no degree of freedom).It is clear that on the basis of geometrical considerations, a 6site cluster that fulfills the conditions of the Hamiltonian operator of Eq. 16 cannot exist, as a point in the three-dimensional space that is equally spaced from all previous 5 sites does not exists.In molecular clusters the magnetic coupling constants are not only dependent on the geometrical features of the cluster; on the contrary, environmental effects could affect their values, and ultimately magnetic symmetries may differ from the geometric ones.It is therefore relevant, for pure amusement, to investigate whether such a highly symmetric 6-site PNTM cluster could exists whose Hamiltonian features symmetries that go beyond the fundamental geometrical limits.Such a system would have an astonishingly sparse Hamiltonian matrix if expressed in a cumulatively spin-adapted basis.
For a 4-site pyramidal cluster (Fig. 2), Eq. 16 reads as The (1234) site ordering guarantees all commutators of such an Hamiltonian with the cumulative partial spin operators to vanish.It is worth mentioning that when applying the graphical approach (Algorithm 1) the (1234) ordering renders uniquely isosceles triangles, namely (123), ( 124), (134), and (234) that make (J JK −J IK ) in Eq. 15 consistently vanish.Are there site permutations equivalent to (1234) in terms of retaining vanishing commutation relations?From inspection of Eq. 17, it emerges that upon (12) permutation, all cumulative partial spin operators still commute with the Hamiltonian operator, making the (2134) site ordering identical to (1234); this equivalence arises from the commutativity of dot products, i.e., Ŝi • Ŝ j = Ŝ j • Ŝi .Any other permutation would lead to one or more non-vanishing commutators, thus compromising the block structure of the corresponding Hamiltonian matrix and the sparsity of the corresponding eigensolutions.

Less Constrained Hamiltonians
Model Hamiltonians may also exist that commute only with the first m cumulative partial spin operators The 4-site square model Hamiltonian, is an example of a cluster that fulfills Eq. 18, for the first two sites, m = 2 (see also Reference 193).In the following, two other examples will be discussed together with the permutation symmetries that can be exploited in those cases.In practical applications, even the simplest commutation relation between a local spin operator and the Hamiltonian operator, Ŝi 2 , Ĥ = 0, may lead to important sparsity, in virtue of the separation of CSFs fulfilling and violating first Hund's rule.

The Diamond Cluster
The Heisenberg Hamiltonian operator of a 4-site cluster in a rhombus shape (Figure 1) reads as where J 13 could be larger, equal or smaller than J 24 .The square cluster with J 13 = J 24 = J diag is a special case of the rhombus.As for the square cluster, the Hamiltonian operator (Eq.20) can directly be derived from Eq. 18, with m = 2 (with some relabeling and it is to be preferred over the (1234) ordering, which instead does not exhibit such a vanishing commutation relation.The permutation group of the rhombus is the D 2 dihedral group of order h = 4, which is made of the identity (E) and the three C 2 rotations around the Cartesian axes (Fig. 1), and it is to be distinguished from the D 2h point group symmetry the system belongs to.The list of N!/h = 6 symmetrically non-equivalent site permutations consists of the (1234), ( 2314), ( 2341), ( 1243), ( 1324) and (2413) orderings.These site orderings are listed in separate columns in Table 1, together with the symmetry equivalent site orderings that are obtained by applying the symmetry operations of the D 2 dihedral group.In Table 1, an S loc = 3/2 has been utilized (see Table SI.1 in the Supplementary Material † for S loc = 1/2).The level of compression (measured by the L 4 -norm of the eigen- b We use J 12 = J 23 = J 34 = J 14 = −1.789,J 13 = −1.000,and J 24 = −2.000. vectors) as a function of different site orderings, for the energetically lowest states with total spin S tot = 0, 1, 3, are also reported in Table 1, for the geometry depicted in Fig. 1.
The L 4 -norm changes substantially as a function of the site ordering.For example, for S tot = 0, the ground-state wave function reaches its maximum compression, L 4 = 1.00, for the (2413) and the (1324) orderings, indicating (within the cumulatively spin adapted basis of CSFs) a single reference wave function.On the contrary, the other four site orderings feature a less compact ground-state wave function with L 4 = 0.69.This difference is to be attributed to the fact that for the (2413) ordering Ŝ2 + Ŝ4 2 , Ĥ = 0, while for any of the other site orderings, We also notice that while for (1324) and (2413) the commutators of the Hamiltonian operator with the 2-site partial spin operators vanish, i.e., Ŝ1 + Ŝ3 2 , Ĥ = Ŝ2 + Ŝ4 2 , Ĥ = 0, the commutator with the three-site partial spin operator does not vanish, and it differs for the two site orderings, namely It is precisely this difference that makes L 4 (1324) ̸ = L 4 (2413) for S tot = 1, 2. It is to be noted that there is no a single site ordering that is best suited for all spin states.However, the (1324) and (2413) orderings are similarly best suited for all spin states, if we consider averaging the L 4 -norms over all spin states.
While identical eigenvectors are to be expected for point group symmetry equivalent site orderings, Table 1 shows that eigenvectors with identical L 4 -norm may also arise for symmetrically nonequivalent site orderings.For example, for the triplet spin state, S tot = 1, the six non-equivalent site orderings further group in four sets.This argument applies irrespectively of the actual magnetic coupling constants considered.This is an interesting observation, indicating that additional symmetries exist that further reduce the number of non-equivalent site orderings.Thus, in the search of an optimal site ordering, instead of searching the entire permutation space (N! = 24) or in virtue of dihedral symmetry considerations, searching the reduced (N!/h = 6) space, it should be possible to further restrict the search over a yet smaller space.This space can be identified by a recursive procedure that relates permutations of tuples to the dihedral symmetry of the cluster and will be introduced in Section 3.

The Kite Cluster
Squares and rhombi are special cases of the more general 4-site kite clusters (Fig. 3), that will be discussed in the present section.The dihedral group of the kite is D 1 (h = 2), which is identical to C 2 in non-standard orientation, where the C ′ 2 axis of the D 1 group is chosen along the z and not the x direction.The symmetrically non-equivalent site permutations are thus N!/h = 12, as summarized in Table 2.The (2413) site ordering is the one that leads to maximum compression, because it is the only ordering for which the commutation relation Ŝ2 + Ŝ4 2 , Ĥ = 0 is fulfilled.This result is also graphically recognized (see Fig. 3) as the (24) pair splits the kite into two isosceles triangles, (241) and (243).As already observed for the rhombus case, symmetry non-equivalent site orderings further group on the basis of the L 4 -norm values, so that for example for the S tot = 1 and S tot = 3 states the 12 symmetry non-equivalent site orderings lead to only 7 groups on the basis of the L 4 -norm value.In Section 3, the further reduction of non-equivalent site permutations will be addressed and a recursive graphical procedure discussed.Moreover, in Table 1 and Table 2 we also observe that the number of groups further reduces to only two for the S tot = 0 state.The further reduction of distinct site orderings for the singlet spin  b We use J 12 = J 14 = −1.638,J 23 = J 34 = −0.894,J 13 = −0.741,and J 24 = −1.000.

Faraday Discussions Accepted Manuscript
states is discussed in detail in Section 5.

Exploiting Permutation Symmetry via the Tree Search Algorithm
In the previous sections it has been shown that different site orderings have an impact on the sparsity of the eigensolutions, as measured by the L 4 -norm of the many-body expansion, and that commutators between cumulative partial spin and the Hamiltonian operators can be used as descriptors of such sparsity.The number of non-equivalent site permutations can trivially be reduced by the order (h) of the point (or dihedral) group symmetry of the considered cluster (N!/h).However, it has been observed that site permutations not related by point group symmetry operations may still be grouped together under the metric of the L 4 -norm.Thus, it should be possible to further reduce the space of non-equivalent site permutations by exploiting additional internal symmetries of the Hamiltonian.A recursive graphical Tree Search Algorithm is presented, that is based on the commutators between cumulative partial spin and Hamiltonian operators, and that is relevant to identify the unique non-equivalent site orderings under the L 4 metric, and thus to reduce the factorial scaling of the permutation space.
A given site ordering, s, defines an array of commutators for an N-site cluster.Two arrays, A (s) and A (t) (where s and t define two site permutations), are considered equivalent if it is possible to establish equivalence between the commutators within the two arrays, otherwise numerically if the L 4 -norms (and thus the wave function structure) of the two site orderings, s and t, are identical.The size of the permutation space can then be reduced by excluding redundant site orderings leading to equivalent commutator arrays.
The effect of the point group symmetry on the array of commutators, A (s), as a whole, can be written as where Ô represents a symmetry operation within the cluster's point group symmetry.One other symmetry we can exploit to reduce the permutation space size is the commutativity of dot products, i.e., Ŝi • Ŝ j = Ŝ j • Ŝi .As in general the two symmetries are not independent, i.e., a point group symmetry operation might already cover the dot product symmetry or vice versa, one should consider the two symmetries simultaneously to span the entire permutation space of unique site orderings.
In order to establish the A (s) = A (t) equivalence we proceed recursively by applying point group symmetry to the commutators , Ĥ , for n values increasingly larger, n = 1, 2, 3, . .., and so on.This is best done graphically, by means of a tree search algorithm.Figure 4 shows the tree search algorithm applied to the kite cluster.The tree is formed by nodes and in each node a subset of sites are marked with circles from which Ŝ(n) is built.For example, the top-left node of Fig. 4  partial spin operators are derived from the level-2 by adding circles with level-specific colors (in Fig. 4 blue circles are utilized for level-2 and green circles are utilized for the newly introduced sites in level-3).Symmetry non-equivalent nodes lead to different branches of the tree.The tree grows until the last level is reached (level-N, for an N-site cluster).The level-N nodes span the permutation space of unique site orderings.In the case of the kite cluster (Fig. 4), there are four distinct nodes in level-2, and a total of seven distinct nodes in level-3.Branching from level-3 nodes to level-4 nodes is trivial and thus not shown in the figure.The seven distinct orderings obtained from the tree search algorithm match the ones identified in Table 2.
The further reduction of distinct site orderings for singlet spin states is explained in Section 5, and there its connection to the Fig. 4 Graphical representation (tree) of the cumulative partial spin operators and their commutation relation with Ĥ for a 4-site kite cluster.The seven nodes at the bottom correspond to the seven non-equivalent site orderings identified in Table 2 on the basis of the L 4 -norm.

Faraday Discussions Accepted Manuscript
graphical tree search algorithm will also be made.

A Hairy Example: The Six-Site Cluster
In this section the tree search algorithm is applied to a larger 6site hexagonal cluster.The regular hexagon is in a D 6 dihedral group (h = 12), and, in absence of the tree search algorithm discussed in Section 3, it would be characterized by N!/h = 60 nonequivalent site permutations.By the application of the tree search algorithm we show that this number can further be reduced.The complete tree is shown in Fig. 5 At the level-2 row we may identify 3 distinct pairs of sites, (12), ( 13) and (14).Any other pair is identical to the three above either for point group symmetry considerations or for the commutativity of the dot product.Thus, the three level-2 nodes indicates that there are only 3 non-equivalent commutators of the type Ŝ(2) 2 , Ĥ .Following the tree search algorithm recursively until level-5 reveals a total of 30 symmetry non-equivalent site orderings, which represents the number of distinct permutations under the L 4 -norm metric (see Table SI.5).The tree search algorithm allows to reduce the permutation search from the initial 6! number of possible site orderings to only 30 distinct site permutations.The number of distinct site orderings reduces to only 12 when the permutation space of the singlet spin states is considered, a substantially smaller number.This further reduction will be discussed in detail in Section 5.

Higher Reversal Symmetry for Singlet Spin States
Table 1 (rhombus), Table 2 (kite) and Table SI.5 (regular hexagon) show that the number of distinct site orderings on the basis of the L 4 -norm values further reduces in the case of singlet spin states, S tot = 0, indicating that an additional symmetry exists that could easily be exploited.
For singlet spin states, the Hamiltonian matrix elements for a given site ordering, s, are identical to the ones obtained when the reverse site ordering s is utilized where to permutation of rows and columns, and so are their L 4 norms.symmetry is not present for spin states with tot > 0, as it would lead to unphysical negative cumulative spins.Therefore, while the (1234) and (2314) site orderings in Table 1 have the same L 4 norms for any spin state in virtue of Eq.21 and Eq.22, for the singlet spin states also the (1234) and (4321) site orderings have the same L 4 norms, precisely because they are the reverse orderings of each other.Similar arguments apply to the and (4231) site orderings for the rhombus, as well as all the equivalent orderings emerging for the singlet spin states of the kite example (Table 2) and the regular hexagon (Table SI.5).
The reversal symmetry for singlet spin states can also be identified graphically in the tree search algorithm, by following the tree from the bottom level to the top level and the ordering of the holes instead of the circles.For example, in the case of the kite (Fig. 4), the (2314) ordering arises from choosing the third node in level-2 and the fourth node in level-3.Reading this path in reverse ordering and following the holes (sites not marked by circles) we would obtain the (4132) ordering.These two site orderings will be in the same group on the basis of the L 4 norm for the singlet spin states.
The equivalence between reverse site orderings identified for singlet spin states explains the number of distinct site orderings found for the Co(II) 3 Er(III)(OR) 4 cubane cluster studied in Reference 186.For the Co(II) 3 Er(III)(OR) 4 cubane an exhaustive site permutation search over all 4! site orderings for the lowest singlet spin state has revealed only 3 distinct site orderings on the basis of the weight of the dominant CSF (see Table 1 of Reference 186).The system is characterized by a C 1 point group symmetry, and the magnetic interactions across the metal centers should in first approximation be described by a 6J-Heisenberg Hamiltonian; thus, no permutation symmetry can be exploited in the tree search algorithm, except the commutativity of the scalar products between the first two chosen sites, that reduces the number of non-equivalent site orderings to 4!/2.However, for the singlet spin states, site orderings and their reverse are equivalent (halving the distinct site orderings), and commutativity can also be applied to the reverse site orderings (another factor of 2 reduction).Thus, despite the low point group symmetry of the system, the permutation space is substantially reduced from 4! to only 3, on the ground of reversal symmetry and commutativity of the first two scalar products.Fig. 5 Graphical representation (tree) of the cumulative partial spin operators and their commutation relation with Ĥ for a 6-site cluster in a regular hexagonal shape.The rows identify the levels corresponding to the
The 30 nodes at the bottom level correspond to the non-equivalent site orderings identified in Table SI.5 on the basis of the L 4 -norm.
a Not allowed as it leads to |∆b| > 2

Proof of the Reversal Symmetry for Singlet States
In the cases where S loc > 1/2, we may assume that there are multiple parallel aligned electrons on each site (first Hund's rule).Thus, a general equation for the Heisenberg Hamiltonian operator reads as where p (and p ′ ) and q run over electrons on site P and Q, respectively.J H ≫ 1 has been introduced to guarantee ferromagnetic on-site couplings.When there is only one electron on each site, Eq. 24 reduces to Eq. 4.
Heisenberg Hamiltonians can be expressed in terms of the exchange two-body spin-free excitation operators using the following relation 185 Ŝp where êqppq = êpqqp = Êpq Êqp − Êpp , the one-body spin-free excitation operator Êpq = ∑ σ =↑,↓ â † pσ âqσ excites an electron of spin σ from the orbital q to the orbital p, and the operators â † pσ and âqσ are the electron creation and annihilation operators, respectively.
Inserting Eq. 25 into Eq.24 we obtain where we replace êppqq and êppp ′ p ′ with 1 because they are products of the number operators of each orbital that always accommodate a single electron.The constant terms are invariant under the ordering reversal operation and will be ignored in the remainder of this document.The ordering reversal operation converts the electron labels p and q as where f is the number of electrons on each site and N is the total number of sites, thus f N is the total number of electrons (for sake of simplicity it is assumed that all sites have the same number of electrons f ).In virtue of the reversal operation, the magnetic coupling constant J PQ in the original site ordering is equal to J P Q in the reversed ordering.Therefore, the equality represents a sufficient condition to verify Eq. 23.
In GUGA, the matrix elements over exchange type two-body operators, êpqqp , are calculated as The abscissa and the ordinate represent the k level and the corresponding cumulative b k value, respectively.In general, paths are followed from the left to the right (direct paths).One advantage of using genealogical branching diagrams over Shavitt graphs for CSFs with only unpaired electrons is that reading a CSF path from right to left (reverse path, follow the k levels in Fig. 6) is equivalent to reversing the CSF, with decreasing segments becoming increasing segments and vice versa.In any direct path, the node on the right of each segment indicates the cumulative b k value, while in any reverse path the node on the left of each segment represents the cumulative b k.When using the genealogical branching diagrams, the segments contributing to a direct coupling element, ⟨i| êpqqp | j⟩ are also the contributing segments to the inverse coupling element, ⟨ ĩ| ê p q q p | j⟩.For example, in Fig. 6 the contributing segments of ⟨i| ê7337 | j⟩ (filled in gray) are also the contributing segments of ⟨ ĩ| ê7337 | j⟩ = ⟨ ĩ| ê2662 | j⟩, (Eq.27 has been used for the orbital label inversion).Although the same segments are used for the matrix element calculations of the original and the reversed twobody operators, the segment values which may not necessarily be the same as the segment values in the k indices.This behavior is shown with an example in Table 5 and Table 6, where the ⟨i| ê3773 | j⟩ and the ⟨ ĩ| ê2662 | j⟩ are evaluated.The individual segment values for the direct and reversed matrix elements are in general not the same, however their products are and numerically proves Eq. 28.A more general proof will be offered in the following.
Instead of calculating the segment values using the original GUGA Table 3 with    ).Segment values outside of the range [3, 7] are not defined and marked with -.W 0 and V 0 values are not reported as at least one of the contributing factors vanishes.
Table 6 Evaluation of ⟨ ĩ| ê2662 | j⟩ = ⟨uduududd| ê2662 |uuduuddd⟩ Hamiltonian matrix element using tabulated value reported in Table 3. Segment values outside of the range [2, 6] are not defined and marked with -.W 0 and V 0 values are not reported as at least one of the contributing factors vanishes.
Table 3 and Table 4 show that the only difference in the contributions to W 0 for x = 0 is the sign of RL and RL segment values; however, as any Hamiltonian matrix elements always contain only one RL and one RL segment, the final products are always the same in sign and value for x = 0 terms.The tables also show coinciding segment values for (d ′ , d) = (1, 1) and (d ′ , d) = (2, 2) for the RL and the RL cases, as well as for the RL case with ∆b = 0.
Next, we prove that , and , which is the same as the ratio for W 1 .The equivalence of ratios holds for all d ′ d = 12 and d ′ d = 21, and both for RL and RL segments types (left as exercise for the motivated reader).
Thus, any RL and RL segments can always be expressed as and any segment product can always be written as simple product of internal RL segments.In the evaluation of segment products, and moving towards the proof of Eq. 28, in virtue of Eq. 30 the r(d ′ k d k , ∆b, b) can be omitted, as it takes identical values for direct (W ) and reversed (V ) segment products.
The details discussed above can be utilized to demonstrate Eq. 28 by induction.The simplest possible loop types are considered as base case of the induction (see Fig. 7).In one case (Fig. 7a The V 1 counterpart of this loop (relying on tabulated values of Table 4) lead to the same result A similar proof holds for the second loop type (Fig. 7b).
In the induction step, we insert a segment to the left of the closing segment; we may refer to the new segment as the n-segment, making the step vector of size n + 1 after insertion (see Fig. 7c and Fig. 7d).Only segment insertions of the d ′ d = 11 and d ′ d = 22 type will be considered, that retain the ∆b value of the level prior to the insertion; other segment pairs lead either to closing loops that is invalid as the node of the inserted segment is not connected to the existing segment, or |∆b| > 2 which makes the coupling element zero, equally proving W 1 ≡ V 1 and thus Eq. 28.As an example, we show the induction step for a generic loop with a W 1 (RL;  Equivalent conclusions are obtained by applying insertions to to d ′ d = 21 loop-opening segments (Fig. 7d).Thus, the induction completely proves Eq. 23 for any Heisenberg Hamiltonians within the singlet spin symmetry.

Conclusions
Within cumulatively spin adapted many-body bases, site reorderings considerably impact the sparsity of the Hamiltonian matrix and its eigenfunctions, allowing methods that approximate the full CI solutions, including the spin-adapted FCIQMC algorithm, to converge rapidly to the exact solution.Specifically, site permutations exist that lead to Hamiltonian matrices with (quasi) unique block diagonal structure.Within each block, the spin-adapted functions (often referred to as configuration state functions, CSFs) have common expectation value of the cumulative partial spin over the first n sites, ⟨ Ŝ(n) 2 ⟩.The servables).The sparsity that follows vanishing commutators is of broad impact.It applies to general systems with S local ≥ 1 2 ; it equally applies to homonuclear and heteronuclear systems; it holds for PNTM clusters whose magnetic centers have same or different local spins, as long as their interactions are dominantly of exchange type.For vanishing commutators the numerical values of the magnetic coupling constants do not have any impact on the sparsity of the Hamiltonian matrix, and sparsity persists whether J pq > 0 or J pq < 0.
Permutation and point group symmetries can be exploited such that the number of non-equivalent site orderings, on the basis of the L 4 -norm metric, is reduced as compared to the factorially growing (N!) permutation space.A tree search algorithm could be utilized to identify the distinct site orderings.The number of distinct site orderings is substantially reduced for the singlet spin states, on the basis of reversal symmetry and commutativity of scalar products.The importance of this strategy resides on its ability to reduce the space of distinct site orderings as compared to the factorial growth of the permutation space with the number of magnetic centers.
The equivalence between reverse site orderings for singlet spin states further reduces the number of distinct site orderings, as found for example for the Co(II) 3 Er(III)(OR) 4 cubane cluster studied in Reference 186.In Reference 186, only three distinct site orderings have been identified upon searching the 4! permutation space of the singlet spin sector.The Co(II) 3 Er(III)(OR) 4 system is characterized by a C 1 point group symmetry; thus, in principle, no point-group symmetry arguments can be invoked to justify the reduction of distinct site orderings.In the present work, we demonstrated the reduction of distinct site orderings on the ground of the reversal symmetry and commutativity of the scalar products.Moreover, via the tree algorithm we have been able to predict what these three orderings are.

(
seniority)are utilized as an example.In equation 3, u and d represent the cumulative spin coupling with S = 1/2 (spin-up) and S = −1/2 (spin-down), respectively.In each CSF of Eq. 3, sites A, B and C are separated by commas.The three CSFs have a common expectation value over ŜA 2 , S A (S A + 1) = 15/4.However, on the basis of the ŜAB 2 = ŜA + ŜB 2 metric, the first two CSFs are grouped together, sharing the expectation value S AB (S AB + 1) = 6, which is different for the third CSF, S AB (S AB + 1) = 2.

Fig. 1 Algorithm 1
Fig. 1 Schematic representation of a rhombus cluster and the partitioning of the triangles on the basis of Algorithm 1. Filled and empty circles represent points in G 1 and in G 2 , respectively.The inset collects the symmetry operations of the D 2 dihedral group.
represents the Ŝ1 + Ŝ3 2 partial spin operator (nodes 1 and 3 are marked with blue circles), from which the Ŝ1 + Ŝ3 2 , Ĥ commutator is evaluated.Each row of the tree is characterized by a unique n value of the Ŝ(n) 2 partial spin operator, starting with n = 2 for the top row.Level-1 nodes are not considered as Ŝ(1) 2 always commutes with any Heisenberg Hamiltonian operators, implying that any choice of level-1 is equally good.We refer to the top row entries as level-2 nodes.Due to the commutativity of dot products, ŜA • ŜB = ŜB • ŜA , we do not need to consider the relative ordering of the first two sites, which are consequently marked in the same color.Each row collects all non-equivalent cumulative partial spin operators, on the basis of point-group symmetry considerations and commutativity of dot products.For example, in the top row of Fig. 4, the node labeled as (1 − 3) is identical to (3 − 1) for commutativity of the dot product Ŝ1 • Ŝ3 = Ŝ3 • Ŝ1 , from which it also follows that Ŝ1 + Ŝ3 2 , Ĥ = Ŝ3 + Ŝ1 2 , Ĥ .Thus, only one between (1 − 3) and (3 − 1) is reported.Similarly, the (1 − 2) node is identical to (2 − 1) for commutativity, and identical to (1 − 4) and (4 − 1) for point-group symmetry considerations.The tree grows downwards to the next row (level 3).The Ŝ(3) 2
ĩ⟩ and | j⟩ refer to the reversed CSFs of |i⟩ and | j⟩, respectively, in which the orbital ordering is reversed and u and d spin couplings flipped; for example, the reversed of |i⟩ = |uduudd⟩ is | ĩ⟩ = |uuddud⟩.A proof of Eq. 23 is given in Section 5.1.Thus, the original and the reversed Hamiltonian matrices are identical, up
where Y k refers to the segment type symbol; the segment value d k = 1 if the k-th coupling of the ket CSF, | j⟩, is u (cumulatively spin coupling with s = +1/2) and d k = 2 if the coupling is d (cumulatively spin coupling with s = −1/2); the d ′ k variable refers to the segment values of the bra CSF, |i⟩; b ′ k and b k are the cumulative spin values, 2S k up to the k-level of the bra and ket functions, respectively; and ∆b k refers to the b k − b ′ k value. 200In Eq. 29, all segment values, W x (Y k ; d ′ k d k , ∆b k , b k ), over the k segments for the same x value are multiplied, and the resulting products over the two x values summed.For exchange type operators, êpqqp , only RL, RL, and RL segment type symbols are to be considered out of all possible combinations (see Tables VI of Reference 200).If k = q, the segment type is RL (tail node), if k = p, the segment type is RL (head node), and the other internal segments are of the RL type.While within the GUGA framework, CSFs and their step vectors are generally represented graphically via Shavitt's graphs, for the special Heisenberg model case (where orbitals can only be singly occupied by u and d spins) a slightly modified version of the simpler genealogical branching diagrams can be utilized.In Figure 6, two CSFs, |i⟩ = |uududdud⟩ and | j⟩ = |uuuddudd⟩, are graphically represented by a dashed and a solid line, respectively.

Fig. 6
Fig. 6 branching diagram representing the | j⟩ = |uuuddudd⟩ and |i⟩ = |uududdud⟩ CSFs.An ordering reversal operation of a CSF can be viewed as reading the same CSF in the reverse path; thus, | j⟩ = |uuduuddd⟩.The loop between the two CSFs in the range k ∈ [3, 7] is filled in gray and indicates the segments contributing to the ⟨i| ê3773 | j⟩ coupling element.
variables (b ′ , b, ∆b, d ′ , d) of the reversed CSFs and the reversed generator, say ê2662 , one could first convert Table 3 to a reversed table, where changes to the b ′ , b, ∆b, d ′ , d are applied (Table4), and then use the variables of the original CSFs and generator, say ê3773 , with segments value collected from the reversed table.
Table 4 is obtained by noting that ifd k = 1 then d k = 2 (an up-spin in the direct branch is a down-spin in the reverse branch) and b k = b k − 1; and vice versa, a d k = 2 is reversed to d k = 1, and b k = b k + 1. Moreover the RL becomes RL and vice versa.This strategy makes the evaluation of the W x values convenient, as one can consistently use the d ′ k d k , ∆b k , b k variables of the direct CSFs.Faraday Discussions Accepted Manuscript Open Access Article.Published on 23 April 2024.Downloaded on 5/15/2024 7:33:39 PM.This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.View Article Online DOI: 10.1039/D4FD00061G with an example.Let us consider the segment values,W 1 (RL; 21, 2, b) = √ (b−1)(b+1) b ) a loop-opening segment of the type d ′ k d k = 12 is immediately closed by the loop-closing d ′ k+1 d k+1 = 21 segment, and in another case (Fig. 7b) a loop-opening segment of the type d ′ k d k = 21 is followed by the loop-closing d ′ k+1 d k+1 = 12 segment.The undulating lines outside the loops represent branches that are identical in the two coupled CSFs, ⟨i| êpqqp | j⟩, outside the range of the generator.The coupling element value of the loop in Fig. 7a, upon taking advantage of Eq. 31 and omitting r(d ′ k d k , ∆b, b), is W 1 (RL; 12, −2, b − 1)W 1 (RL; 21, 0, b) 1) b(b + 2) .

Fig. 7
Fig. 7 All possible loops that can appear in Heisenberg Hamiltonian matrices, represented via genealogical branching diagrams.The simplest loops are (a) for d ′ d = 12 and (b) for d ′ d = 21 loop openings.Arbitrary length-n loops are viewed as inserting n − 2 segments of d ′ d = 11 and/or d ′ d = 22 couplings in between the loop-opening and -closing segments: (c) a loop opened by d ′ d = 12 segments and (d) a loop opened by d ′ d = 21 segments.Indices of the loop-opening segment (m) and the one before the loop-closing segment (m + n − 2) are presented for (c) and (d).The b value of the node marked in red is β (see Eq. 34).Undulating lines denote any possible couplings before and after the opening and closing segments.
unique block structure is explained by the commutation relation Ŝ(n) 2 , Ĥ ; for some site orderings such commutator vanishes, allowing the Hamiltonian and the cumulative partial spin operators to share common eigensolutions (compatible ob-Faraday Discussions Accepted Manuscript Open Access Article.Published on 23 April 2024.Downloaded on 5/15/2024 7:33:39 PM.This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.View Article Online DOI: 10.1039/D4FD00061G

Table 1
Symmetry Non-Equivalent Site Orderings and Corresponding L 4 -Norms a for a 4-Site Rhombus Cluster with S loc = 3/2 b a The largest L 4 -norms are highlighted in bold and indicate the best site orderings.

Table 2
Symmetry Non-Equivalent Site Orderings and Corresponding L 4 -Norms a for the 4-Site Kite Cluster with S loc = 3/2 b a L 4 -norm calculated from the CI eigenvectors for each site ordering.

Table 4 Table
of reversed segments as a function of the original variables before applying the reversing transformation.The segment values are denoted as