Approximate eigenvalue determination of geometrically frustrated magnetic molecules

Geometrically frustrated magnetic molecules have attracted a lot of interest in the field of molecular magnetism as well as frustrated Heisenberg antiferromagnets. In this article we demonstrate how an approximate diagonalization scheme can be used in order to obtain thermodynamic and spectroscopic information about frustrated magnetic molecules. To this end we theoretically investigate an antiferromagnetically coupled spin system with cuboctahedral structure modeled by an isotropic Heisenberg Hamiltonian.


Introduction
The complete understanding of small magnetic systems such as magnetic molecules is compulsorily connected to the knowledge of their energy spectra. From the energy spectra all spectroscopic, dynamic, and thermodynamic properties of the spin systems can be obtained. Unfortunately, an exact calculation of the spectrum is often restricted due to the huge dimension of the Hilbert space even if one works within the most simple isotropic Heisenberg model. The dimension grows for a system of N spins with spin quantum number s exponentially and is (2s + 1) N .
In order to get insight into the properties of large magnetic molecules one can access several numerical methods which have been developed in the past. Of course, the ultimate method of choice would be an exact numerical diagonalization yielding the complete energy spectrum. In recent years there has been enormous progress on extending the range of applicability of the exact numerical diagonalization of the Heisenberg model. To this end the use of spin-rotational symmetry [ 1,2] in combination with point-group symmetries [ 3,4,5,6] can be of great advantage with respect to a reduction of computational requirements, i.e. a need of hardware resources and computation time. Apart from the exact numerical diagonalization technique the magnetism of magnetic molecules can be very well investigated using complementary methods such as Density Matrix Renormalization Group (DMRG) [ 7,8,9], Lanczos [ 10], or Quantum Monte Carlo (QMC) [ 11,12,13] techniques. Nevertheless, also these methods suffer from theoretical limitations, QMC for instance in systems with geometric frustration.
Currently magnetic molecules which exhibit geometric frustration are of special interest due to the richness of physical phenomena like plateaus and jumps of the magnetization for varying field as well as special features of their spectra such as low-lying singlets [ 14,15,16,17]. In this respect a lot of insight has been obtained by investigating molecular representations of archimedean-type spin systems [ 18,19], i.e. systems in which participating spins occupy the vertices of an archimedean solid. Such representations have already been synthesized several years ago and exist for example as {Cu 12  In this context the molecular compound {Mo 72 Fe 30 } is probably one of the most investigated magnetic molecules. However, a theoretical explanation of its special physical properties was so far given mostly by considering purely classical models [ 23,24] or directly related quantum mechanical counterparts like the rotational-band model [ 25,26].
In this paper we want to show how the approximate diagonalization technique which has been developed and applied to unfrustrated, i.e. bipartite, magnetic molecules in Ref. [ 6] can be used in order to determine the energy spectrum of geometrically frustrated magnetic molecules. The idea of this technique is to diagonalize the full Hamiltonian in a reduced basis set. The basis set itself is an eigenbasis of the rotational-band Hamiltonian. Such an ansatz was already used by Oliver Waldmann [ 27] in order to interpret inelastic neutron spectra of {Mo 72 Fe 30 } [ 18]. We will demonstrate that in contrast to bipartite systems for a frustrated spin systems the rotationalband states of all non-trivially different sublattice colorings have to be taken into account in order to achieve a reliable convergence of energy levels. Throughout this paper we use a spin system with cuboctahedral structure and spin quantum numbers s = 1 and s = 3/2 as well as an antiferromagnetic coupling as an archetypical example of a frustrated magnetic spin system. This paper is organized as follows. In Sec. 2 the general theoretical description of the system within the Heisenberg model as well as general remarks on the use of irreducible tensor operators and point-group symmetries are given. In Sec. 3 the theoretical basics of the approximate diagonalization within the isotropic Heisenberg model are briefly reviewed and specified for the cuboctahedral spin system. The convergence behaviour of the approximate diagonalization is displayed and deeply discussed as well as the specific heat and zero-field magnetic susceptibility for a cuboctahedron with s = 3/2. Furthermore, a numerically based finding of an approximate selection rule is reported. This paper closes with a Summary in Sec. 4.

Theoretical method
In order to model the physics of antiferromagnetic molecules it has been shown that an isotropic Heisenberg Hamiltonian with an additional Zeeman-term and an antiferromagnetic nearest-neighbor coupling provides the dominant terms. Such a Hamiltonian looks like The indices of the sum are running over all pairs < i, j > of interacting spins i and j. The first part consisting of the sum over single spin operators s ∼ (i) at sites i interacting with the coupling strength J ij < 0 refers to the Heisenberg exchange whereas the second part -the Zeeman-termcouples the total spin S ∼ to an external magnetic field B. Taking without loss of generality J ij = J for interacting spins the Heisenberg part assumes the following form where each coupling is counted only once. Since due to SU (2) symmetry the commutation relations H ∼ Heisenberg , S ∼ = 0 hold, a common eigenbasis {|ν } of H ∼ Heisenberg , S ∼ 2 and S ∼ z can be found and the effect of an external magnetic field B = B · e z can be included later, i.e.
Here E ν denotes the energy eigenvalues, |ν the eigenstates, and M ν denotes the corresponding magnetic quantum number. For the matrix representation of the Heisenberg Hamiltonian (2) the irreducible tensor operator technique is used [ 1,2,6]. The Heisenberg Hamiltonian is expressed in terms of irreducible tensor operators and its matrix elements are evaluated using the Wigner-Eckart-theorem Equation (4) states that a matrix element of the q-th component of an irreducible tensor operator T ∼ (k) of rank k is given by the reduced matrix element α S||T ∼ (k) ||α ′ S ′ and a factor containing a Wigner-3J symbol [ 28]. The basis in which the Hamilton matrix is set up is of the form |α S M . α refers to a set of intermediate quantum numbers given by addition rules when coupling single spins s ∼ (i) to a total spin S ∼ with spin quantum number S. The underlying spin-coupling scheme directly influences the form of the irreducible tensor operator T ∼ (k) and further on the successive calculation process of the reduced matrix elements in Eq. (4).
By using irreducible tensor operators and the Wigner-Eckart-theorem it is possible to drastically reduce the dimensionality of the problem, i.e. of the Hamilton matrices which have to be diagonalized numerically. The total Hilbert space H can be decomposed into subspaces H(S, M = S). Such a decomposition results in block factorizing the Hamilton matrix where each block can be labelled by the total spin quantum number S.
Additionally point-group symmetries lead to a further reduction of the matrices. Apart from that, states are labelled by the irreducible representation they are belonging to. This is very helpful in understanding the physics of the system. The incorporation of these symmetries results in symmetrized basis states which are constructed by the projection operator [ 29] Here l n denotes the dimension of the n-th irreducible representation of the point-group G which is of order h. G ∼ (R) refers to the symmetry operations of G and χ (n) (R) denotes its character with respect to n. The effect of the symmetry operation G ∼ (R) on basis states of the form |α S M is discussed in Refs. [ 3,5,6].
In this section we follow a method of calculating approximate eigenstates and eigenvalues of the Heisenberg Hamiltonian in Eq. (2) and determine approximately the energy spectrum of a spin system with cuboctahedral structure. In this system 12 spins of spin quantum number s occupy the vertices of a cuboctahedron interacting along its edges. The geometrical structure of the cuboctahedron is shown in Fig. 1. A detailed description of the approximation scheme is given in Ref. [ 6] where it was successfully applied to determine approximate energy spectra and thermodynamic properties of spin rings, i.e. bipartite systems. This approximation rests on the idea to diagonalize the full Hamiltonian in a reduced basis set. The basis set itself is an eigenbasis of the rotational-band Hamiltonian which from the point of view of perturbation theory can be understood as an approximation to the full Hamiltonian. For bipartite systems this (zeroth order) approximation is already very good [ 25].
The Heisenberg Hamiltonian can be decomposed into two parts like where H ∼ RB is the rotational-band Hamiltonian [ 25,26,31] and H ∼ ′ is an operator containing the remaining terms. The rotational-band Hamiltonian which is an effective quantum mechanical Hamiltonian based on classical assumptions looks like Here N denotes the number of spins within the system and N s the number of sublattices the classical ground state of the system is composed of. The prefactor −DJ/(2N ) can be seen as the effective coupling strength between the total spin S ∼ and the sublattice spins S ∼ n which arise from coupling all single spins s ∼ (i) belonging to the n-th sublattice. The value of D directly depends on the system. It is chosen such that the energy of the ferromagnetic state of the system is matched; for the cuboctahedron it is D = 6.
The full Heisenberg Hamiltonian is now diagonalized within a reduced set {|φ i }, i = 1, . . . , n red , of basis states of H ∼ RB yielding approximate eigenstates and eigenvalues of H ∼ . The set of approximate basis states is energetically ordered. Before discussing the results of our approximate diagonalization we would like to characterize the eigenbasis of H ∼ RB . Figure 2 shows the low-lying part of the energy spectra of H ∼ RB for a cuboctahedron with s = 1 and s = 3/2 given by the rotational-band model. They exclusively consist of parabolas -so-called rotational bands. The eigenvalues E RB (S 1 , S 2 , S 3 , S) of H ∼ RB depend only on the spin quantum numbers of the sublattice spins S n , with n = 1, 2, 3, and of the total spin S. Corresponding eigenstates are trivial and analytically given as |α S 1 S 2 S 3 S M . The additional quantum number α refers to a set of intermediate spin quantum numbers which appears when coupling single spins s ∼ i of the same sublattice to the corresponding sublattice spins S ∼ i and further on coupling the sublattice spins to the total spin S ∼ . With regard to the set of intermediate spin quantum numbers α the number of all possible ways of constructing states characterized by fixed values of the sublattice and total spin quantum numbers determine the degeneracy of energy levels in the spectrum of the rotational-band Hamiltonian.
In the case of s = 1 the lowest band in Fig. 2 is given by states |α S 1 S 2 S 3 S M with sublattice spin quantum numbers S 1 = S 2 = S 3 = S max = 4 · 1 = 4 while the second band is given by a deviation of one sublattice spin of 1, i.e. |α (S max − 1) S max S max S M and permutations thereof. The other bands can then be constructed by introducing additional deviations of the sublattice spin quantum numbers. The energy spectrum with s = 3/2 shown in Fig. 2 can be constructed accordingly.
Additionally several energy levels in Fig. 2 are colored. This coloring refers to so-called superbands. A super-band consists of those rotational bands for which the sum of sublattice spin quantum numbers is the same. In contrast to bipartite systems (see Ref. [ 6]) the spectrum of the rotationalband Hamiltonian for the cuboctahedron is much denser at low energies and only the first three super-bands are well separated from the others.
The classical ground state of the system plays the key role within the approximate diagonalization. The better the quantum mechanical system can be approximated by a classical picture of the ground state the more effectively the approximate diagonalization works, i.e. the faster the approximate energy eigenvalues converge towards the exact values. From a purely classical point of view the ground state of a cuboctahedron exhibits a three-sublattice structure and is infinitely degenerate since there coexist coplanar and non-coplanar vector orientations [ 32]. Here it should be emphasized that within the approximate diagonalization, i.e. within a quantum mechanical treatment, coplanarity of the classical ground state does not play a role, but the coloring of the classical ground state does, since there is a direct impact on the set of intermediate spin quantum numbers α in the rotational-band states |α S 1 S 2 S 3 S M which are taken into account as basis states in the approximate diagonalization.  The classical ground state of the cuboctahedron exhibits 24 colorings of the spins which can be -by group theoretical considerations -decomposed into two families that are invariant under operation of the full point-group symmetry O h [ 19]. Following Ref. [ 19] these families will be denoted as Γ(q = 0)-family and M -family. It has also been shown that those irreducible representations of O h which form the Γ(q = 0)-and M -families are found in the low-lying part of the spectrum of the quantum cuboctahedron with half-integer spins s = 1/2, 3/2, 5/2 [ 19]. Figure 3 shows the colorings of the different classical ground state families of the cuboctahedron. In a classical picture each color refers to a sublattice with all spins pointing into the same direction. The angle between classical spins belonging to different sublattices is 120 • .
In order to calculate the approximate spectrum of the system one is now left with the construction of basis states of the form |α S 1 S 2 S 3 S M , i.e. quasi-classical states. Therefore spins belonging to the same sublattice have to be coupled to yield the total sublattice spins S ∼ 1 , S ∼ 2 and S ∼ 3 . Afterwards these total sublattice spins are coupled to the total spin S ∼ . The underlying coupling scheme is given by a classical ground state, i.e. a coloring from Fig. 3 and incorporated in the quantum number α. In the following the resulting basis states will be labelled with respect to the classical reference state. To this end one has to distinguish between basis states of the form The notation of the set of intermediate quantum numbers and the subscript of the states now directly point to the classical ground state colorings, i.e. the underlying coupling scheme. It is important to note, that each of these four basis sets spans the same Hilbert space H(S, M ).
In the case of the M -family three different colorings exist. When diagonalizing the Hamiltonian approximately while additionally using point-group symmetries one has to restrict to those symmetry groups where the symmetry operations do not alter the sublattice structure. A symmetry operation has no impact on the sublattice structure if the corresponding spin permutation results in recoloring of spins where all spins of a given sublattice maintain the same color, i.e. subscript (r: red, y: yellow, b: blue). For example an operation which leads to a cyclic permutation of colorings like has no impact on the sublattice structure. While the sublattice structure of the Γ(q = 0)-family is left invariant under all symmetry operations of O h , the sublattice structure (i.e. coloring) of the M -family is not. In order to restore the sublattice invariance of the classical ground state belonging to the M -family the basis states for the approximate diagonalization can be constructed from a usually overcomplete set M of basis states by an orthonormalization procedure. This set consists of rotational-band eigenstates from each coloring of the M -family, i.e. is given by while the index i is taking the values 1, . . . , n red . Here n red reflects the overall number of states contained in the incorporated rotational bands which is independent of the choice of the coupling scheme.
Since the underlying coupling scheme is different regarding the chosen coloring these states have to be converted into an -in general -arbitrary reference scheme A before calculations can be performed. A transition between states belonging to different coupling schemes can be calculated using general recoupling coefficients [ 33,34]. For the transition of a basis state |µ 1 S M M1 referring to a coloring M 1 of the M -family into the scheme A one yields |α S M A = α,S1,S2,S3 where the summation indices indicate that the summation is running over all valid combinations of values for the intermediate spin quantum numbers α as well as for the sublattice spin quantum numbers. The transitions between states of colorings M 2 and M 3 into the scheme A can be obtained analogously.   ∼ Heisenberg using only a reduced basis set of states from the lowest rotational band in the system with s = 3/2. Obviously the classification of the quasi-classical states is independent of the single spin quantum number s. Solely, the energies are changed when approximately diagonalizing the Hamiltonian of a system with different single spin quantum numbers. Based on a classical ground state belonging to the Γ(q = 0)-family the states belong to the irreducible representations A 1g (S = 0) and A 2g , E g (S=1). Taking the M -family as a starting point for the construction of basis states the states belong to A 1g , E g (S=0) and A 1g , E g , T 2u (S=1) where the T 2u -state has apart from its intrinsic degeneration, i.e. the dimension of the irreducible representation T 2u , a twofold multiplicity.
The difference in the geometric properties of the lowest states of the chosen ground-state family, expressed by their decomposition into different irreducible representations, should directly lead to a convergence behaviour which depends on the choice of the underlying coloring. Since the approximate basis set contains energetically ordered states (from lowest to higher), low-lying states are expected to converge faster with growing number of basis states used for the approximate diagonalization [ 6].

Convergence of individual colorings
In Fig. 5 the approximate eigenvalues of a cuboctahedron with s = 1 and s = 3/2 in dependence on the number of incorporated bands within the S = 0 subspaces are shown. As underlying classical ground states for each system the Γ(q = 0)-family as well as one member of the M -family have been chosen. In every subfigure the last column refers to results from a complete diagonalization within H(S = 0). Since the Γ(q = 0)-family is invariant under O h , the full point-group symmetry of a cuboctahedron (O h ) was used in order to classify the states. As mentioned before the individual members of the M -family are not invariant under O h but under D 2 . Thus, in this case the D 2 point-group symmetry was applied.
Looking at Fig. 5 it becomes apparent that the convergence is directly dependent on the choice of the underlying classical ground states. In agreement with theoretical expectations states which correspond to low-lying eigenstates of the rotational-band model converge faster than higher-lying states. Additionally, by comparing the spectra of s = 1 and s = 3/2 it can be seen that the ground state energy converges more rapidly with increasing single spin quantum number s. Nevertheless, apart from some states which exhibit a quite regular and fast convergence behaviour the convergence is -overall -rather poor. Especially the so-called low-lying singlet 1   case of half-integer spins s below the first triplet contributions, converges very slowly. This means that it is poorly approximated by one or a few eigenstates of the rotational-band Hamiltonian, but instead is a superposition of very many basis states. One significant property which becomes evident when tracing single energy levels is that the convergence is stepwise -at least in the low-lying part of the spectra. Regions in which the ground state is considerably lowered are exemplarily marked with dashed boxes in the case of a cuboctahedron s = 3/2 assuming a classical reference state of the Γ(q = 0)-family. The numbers below the boxes refer to the number of magnons existing in the states of those rotational bands which are incorporated within the marked regions. Obviously the ground state energy is lowered whenever super-bands with an even number of magnons are incorporated in the approximate diagonalization. The observation of a stepwise convergence leads to a helpful additional approximation given by an approximate selection rule which is discussed below (see Sec. 3.4). It should be mentioned here that qualitatively the same convergence behaviour can be observed within subspaces of S > 0. However, additional information cannot be extracted from graphical visualizations of the convergence in these subspaces, thus they will not be presented here.

Convergence of combined colorings
In order to improve the convergence of the energy levels in comparison to the behaviour shown in Fig. 5 the decomposition of the quasi-classical states in Fig. 4 directly leads to the starting point. Since the quasi-classical states, that result from different colorings, can be classified according to irreducible representations of 2A 1g and E g , it is a straightforward task to use linear combinations of rotational-band eigenstates of the Γ(q = 0)-family as well as of the M -family as basis states for the approximate diagonalization. The set of basis states results from an extension of M in Eq. (9) to and subsequent orthonormalization of the incorporated states. In Fig. 6 the energy difference between selected low-lying energy levels and the ground state of the system is displayed for a cuboctahedron s = 1 within S = 0 subspace. The difference is shown in dependence on the number of incorporated rotational bands. The last column refers in both subfigures to the exact values, taken from a complete diagonalization. The dotted red line indicates the energy difference to the E g -state.
In the first case (left subfigure) only rotational-band states of the Γ(q = 0)-family are taken into account for the approximate diagonalization. In the second case (right subfigure) linear combinations of states of all four colorings are used. The first noticeable difference is that when taking basis states of all colorings into account the low-lying levels start and remain in close proximity to their final (true) value. This is especially obvious when comparing the convergence of the lowest A 1g -state in both subfigures. The second difference is that when using linear combinations of states of all four colorings the convergence is smoother and more rapidly. This can be traced back to the fact that the different colorings contribute low-lying states of different irreducible representations which in the other families would only be available as high-lying states. Therefore, it is advantageous to use fewer bands but basis states of all classical colorings for the approximate diagonalization.
A short explanation might be appropriate in order to understand why the energy differences sometimes increase when taking more states into account. This is due to the fact that the ground state in such cases converges more rapidly than the excited states. Although all approximate eigenvalues improve on an absolute scale the differences get worse for a moment.

Results
In this section we like to present how good thermodynamic observables can be approximated. Figure 7 shows the specific heat C(T, B) (left) and the magnetic susceptibility dM/dB (right) for zero field B = 0 of a cuboctahedron s = 3/2. The energy spectrum of this system can be completely calculated using irreducible tensor operator technique in combination with a D 2 point-group symmetry [ 35]. In Fig. 7 the exactly calculated specific heat and the susceptibility are compared with results from an approximate diagonalization using only states of the lowest rotational band in each coloring in order to set up the basis. We also show how these observables look like if evaluated with only the lowest rotational band (L-band) of H ∼ RB (red colored energy levels in Fig. 2).  As one can see the low-temperature specific heat is very sensitive to the structure of lowlying levels. In the rotational-band model the lowest band (L-band) is gaped and describes a rotor. Therefore, the specific heat is suppressed at small temperatures and at higher temperatures approaches 3/2 k B . The approximate diagonalization, although only taking the lowest (L-band) states of all four coloring into account, achieves already an improvement for low-temperatures. The little Schottky-peak is the result of low-lying singlets and rearranging triplets. At higher temperatures this approximation displays the same behavior as the pure L-band which is expected since only L-band like states are incorporated.
In contrast to the specific heat which directly reflects the density of states within a certain energy interval the magnetic susceptibility reflects the density of magnetic states. To that end it is not surprising that the magnetic susceptibilities calculated from the approximate spectra do not differ considerably. Since the exact spectrum also possesses a rotational-band like behaviour at the lower edge the contributions from the approximate diagonalization and from the L-band reproduce the steep rise of the susceptibility for low temperatures very well. Obviously, the exact thermodynamical properties cannot be reproduced properly in the case of increasing temperature because the approximate spectra only contain a fraction of the energy levels of the full Hilbert space H.

Approximate selection rule
As mentioned before (Sec. 3.1) the stepwise convergence behaviour leads to an approximate selection rule. Using this selection rule the computational effort when setting up the Hamilton matrices can be further reduced. In Fig. 5 it was shown that when diagonalizing in the Γ(q = 0)family, regions can be marked in which single energy levels are affected by taking into account additional rotational-band states. The same can be done for the M -family. As it was already mentioned, the marked regions coincide with the incorporation of rotational-band states belonging to the same super-band.
Looking at Fig. 5 the approximate energy eigenvalues of the ground state are lowered whenever states |γ S 1 S 2 S 3 S M are additionally included which belong to a super-band with an even number of magnons. Since the convergence behaviour is state-sensitive the energy of the lowest states belonging to the 3-dimensional irreducible representation T 1u , for example, is considerably lowered whenever super-bands with an odd number of magnons are included. This relation between the incorporation of super-bands and the immediate affection on the energies of certain approximate eigenstates can also be observed if the classical ground state belongs to the M -family.
As a result of the aforementioned observations a simple rule can be conjectured for low-lying energy levels |S 1,b S 2,b S 3,b S M . The matrix elements S 1,a S 2,a S 3,a S M |H ∼ |S 1,b S 2,b S 3,b S M connecting these states with other states |S 1,a S 2,a S 3,a S M are one magnitude (or more) bigger than other matrix elements if |n a − n b | mod 2 = 0 , where n a = 3 i=1 S i,max − 3 i=1 S i,a and n b (similarly evaluated), denote the number of magnons of the super-bands the states are belonging to. The Hamilton matrices approximately split up according to Eq. (12). Thus, simultaneously computation time is saved and the dimensionality of the problem is reduced.

Summary
In this paper an approximate diagonalization scheme was proposed and used in order to determine the energy spectra of geometrically frustrated spin systems. As an example the cuboctahedron was discussed. It was shown that the convergence is clearly dependent on the coloring of the underlying classical ground state, i.e on the coupling scheme. Therefore, the approximation can be improved by using an adapted set of basis states that originates from rotational-band states of all possible sublattice colorings. Furthermore, it was shown that the approximate diagonalization is useful in order to study the low-temperature thermodynamics of geometrically frustrated systems, even in its simplest form when only states from one rotational band within each coloring are taken into account.
Although the necessary calculations are rather involved, especially for huge frustrated systems, the approximate diagonalization can be a valuable method. Numerically the evaluation of recoupling coefficients constitutes the strongest challenge, i.e. calculations are rather limited by runtime than by available memory. However, recent developments show that highly parallelized program code or public resource computing can help to overcome this barrier.