Estimation of biquadratic and bicubic Heisenberg effective couplings from multiorbital Hubbard models

We studied a multi-orbital Hubbard model at half-filling for two and three orbitals per site on a two-site cluster via full exact diagonalization, in a wide range for the onsite repulsion U, from weak to strong coupling, and multiple ratios of the Hund coupling J H to U. The hopping matrix elements among the orbitals were also varied extensively. At intermediate and large U, we mapped the results into a Heisenberg model. For two orbitals per site, the mapping is into a S = 1 Heisenberg model where by symmetry both nearest-neighbor (S i ⋅ S j ) and (Si⋅Sj)2 are allowed, with respective couplings J 1 and J 2. For the case of three orbitals per site, the mapping is into a S = 3/2 Heisenberg model with (S i ⋅ S j ), (Si⋅Sj)2 , and (Si⋅Sj)3 terms, and respective couplings J 1, J 2, and J 3. The strength of these coupling constants in the Heisenberg models depend on the U, J H, and hopping amplitudes of the underlying Hubbard model. Our study provides a first crude estimate to establish bounds on how large the ratios J 2/J 1 and J 3/J 1 can be. We show that those ratios appear rather limited and, as a qualitative guidance, we conclude that J 2/J 1 is less than 0.4 and J 3/J 1 is less than 0.2, establishing bounds on effective models for strongly correlated Hubbard systems. Moreover, the intermediate Hubbard U regime was found to be the most promising to enhance J 2/J 1 and J 3/J 1.


I. INTRODUCTION
The study of the one-dimensional spin-one (S = 1) Heisenberg chain by Haldane [1], with only nearestneighbor spin-spin interactions (called here "quadratic" interactions), and the prediction, and subsequent confirmation, of a spin liquid gapped ground state with protected edge states, was seminal for the start of the field of topological materials.The Haldane chain has been physically realized in several materials, such as CsNiCl 3 [2], AgVP 2 S 6 [3], NENP [4], and Y 2 BaNiO 5 [5], and recently theory predicted that doping of the fermionic two-orbital Hubbard version of the idealized Haldane chain may lead to hole pairing and eventual superconductivity [6,7].Earlier related work employing t − J model approximations also predicted superconductivity with doping although strongly competing with ferromagnetism [8].
While the solution of the Heisenberg S = 1 chain by Haldane was mathematically elegant, intuition was provided later by Affleck et al. [9] when they solved exactly an extension of the original quadratic Hamiltonian by adding "biquadratic" terms.In this exactly solvable point, the magnitude of the ratio biquadratic to quadratic couplings is β = 1/3 [9].At this special point, the model has properties qualitatively similar to those of the Haldane chain, with a unique spin-gapped ground state, exponentially decaying spin-spin correlations, and S = 1/2 spins at the edges when open boundary conditions are used.
The primary goal of the investigations reported in this publication is to study whether the more realistic electronic two-orbital Hubbard model realization of the Haldane chain recently introduced [6] can at large and/or intermediate Hubbard U and Hund J H couplings reach the biquadratic/quadratic ratio β = 1/3 when fermionic vs. pure spin Hamiltonian models are compared at low energies.Specifically, here we solve exactly the two-site problem of the fermionic model and represent the lowest energy states using the generalized Heisenberg quadraticbiquadratic model in a vast region of parameter space, including varying the elements of the hopping matrix.Our conclusion is that it is indeed possible to reach the Affleck et al. point by suitably selecting the values of U and J H . On the other hand, for the Bethe-Ansatz solvable case β = 1 we conclude that it cannot be reached using the fermionic system we studied.Our efforts were extended to the three-orbital per site Hubbard models as well, allowing us to deduce upper bounds for the biquadratic and triquadratic Heisenberg couplings emerging at large Hubbard U and low energy.

II. PREVIOUS INVESTIGATIONS USING QUADRATIC-BIQUADRATIC S = 1 SPIN MODELS
Interest in spin Heisenberg models with spin higher than 1/2 started developing years ago in the context of finding exactly solvable Hamiltonians, in dimension one or more, to uncover disordered spin liquid ground states in antiferromagnets.Of particular interest were valence bond (VB) states, which could serve as toy models for the ideas of Anderson using S = 1/2 resonant valence bonds related to high-T c superconductivity [10].Affleck and collaborators extended the notion of VB states to spins higher than 1/2 [9], as explained in the introduction.For S = 1 adding a biquadratic nearest-neighbor term with coupling J 2 in addition to the standard (quadratic) Heisenberg interaction with coupling J 1 , they found that for β = J 2 /J 1 = 1/3 (with J 1 and J 2 both positive, thus antiferromagnetic) the ground state is exactly solvable and indeed made out of valence bonds.The same model but for the case |β| = J 2 /J 1 = 1 has been solved using the Bethe Ansatz method [11].At this special point the ground state is gapless with a power-law decay.This point, with J 1 > 0 and J 2 < 0, could separate the spin liquid gapped phase for β > −1 from a dimerized phase for β < −1.Our analysis below indicates that Affleck et al.'s case β = 1/3 indeed can be realized with a twoorbital per site electronic model at intermediate Hubbard U , but the ratio |β| = 1 is too large and will require more general electronic models.
To summarize, the isotropic S = 1 Heisenberg model with a biquadratic term has been previously studied.The phase diagram of the model in 1D was obtained via DMRG [12].These authors verified that for β = 1/3 the ground state is indeed a VB state.In addition they obtained the following phases: (i) For J 1 > 0 and J 2 = 0 the system has a non-degenerate disordered ground state with antiferromagnetic spin correlations that decay exponentially indicating the presence of a gap in the spectrum (i.e. the Haldane state); (ii) at β = 1/3 with J 1 and J 2 both positive, the system has the VB ground state with a spin gap in the spectrum, as described before (i.e. the Affleck et al. state); (iii) β = 1 with both J 1 and J 2 positive, indicates the critical point where the Hamiltonian is integrable with a gapless ground state [13]; (iv) for β = 1 with both J 1 and J 2 negative, the model is also integrable and it has a gapless ground state; (v) for J 1 = 0 and J 2 = −1 the system is in a dimerized state; this dimerized state is the ground state when J 2 < 0 and |β| ≥ 1; finally, (vi) the ground state is FM for J 1 = 0 and for β < −1 and for J 1 < 0 and |β| ≤ 1.The region between J 1 > 0 and |β| ≤ 1 has a gapped ground state [14].Several of these phases have been theoretically confirmed via a variety of approaches such as calculation of static and dynamic structure factors using recursion methods [15].
Moreover, recent efforts by one of the present coauthors and collaborators searched for spin liquids in two dimensions focusing on the SU(3) point where the strength of the quadratic and biquadratic interactions are equal, and adding further interactions [16,17].Indeed spin liquids were unveiled for the spin-only models, a conceptually interesting result.But, as already explained, it is difficult to establish which electronic fundamental multiorbital model can realize these complex spin models in the limit of large U , with the exception of the state of Affleck et al. which is reachable with the two-orbital per site model studied here.Our investigations establish limits based on basic Hubbard models on what range of β could be realized in practice.For larger values of β more complex fermionic models will be required.
In addition, it was shown that for certain values of parameters higher spin Heisenberg Hamiltonians in one dimension posses conformal invariance, property that allows an analytical determination of critical exponents [18].The integrable high-spin Heisenberg models are given by a Hamiltonian with a polynomial form in powers of nearest-neighbor Heisenberg interactions rang-ing from 1 to 2S.This was demonstrated via a mapping into the Wess-Zumino-Witten model at specific values of the Hamiltonian parameters [19,20].Various numerical studies of higher spin Heisenberg Hamiltonians were performed to understand whether the higher spin anisotropic Heisenberg Hamiltonians belong to the same universality class as the S = 1/2 isotropic model or, instead, the isotropic integrable higher spin ones [21][22][23].
The isotropic S = 3/2 Heisenberg model has not been as much explored.The isotropic and anisotropic cases were studied in order to determine whether they belong to the same universality class as the S = 1/2 Heisenberg model, which was confirmed using Lanczos and DMRG approaches [9,[21][22][23].Other authors have explored the S = 1 biquadratic model in two dimensions using DMRG in the context of high-T c superconductors finding nematic phases [16], although at robust J 2 /J 1 .
Spin 1 systems are also realized in the area of two dimensional ruthenates [24] often using three-orbital per site Hubbard models with four electrons in those three orbitals, leading to a net S = 1 per site.Rich phase diagrams were reported.But in these ruthenates S = 1 effective Hamiltonians are rarely employed.Spin 1 systems often appear also in the area of iron superconductors because Fe 2+ , with n = 6 electrons in the 3d shell, is the usual iron valence, either in planes or ladders.However, these iron-materials are considered to reside in the intermediate U region [25,26] and, again, they are not often described via purely spin systems but with multi-orbital electronic models instead [27].
Our study also has limitations.For example, the addition of a Zeeman magnetic term to the biquadratic S = 1 model was explored using DMRG [28], and a spin nematic phase was observed in a triangular lattice [29].The addition of single-ion anisotropy to the S = 1 spin Heisenberg model was studied using quantum Monte Carlo and series expansions [30], and for the model with biquadratic term [31] with density matrix renormalization group (DMRG).The addition of a next-nearest neighbor term to the S = 1 Heisenberg model with biquadratic coupling was explored as well with DMRG [32].More recently, research on this model has focused on entanglement and topological properties [33,34].Because our study relates to a single bond, we cannot distinguish between square and triangular lattices.Including more than a single bond, terms such as (S i • S j )(S i • S k ) with sites (i, j, k) belonging to the same plaquette, will also appear in the large U expansion.Note that the models described in this paragraph have either a Zeeman term, single-ion anisotropy, or next-nearest neighbor interactions.Thus, it is too early to make statements on whether these models can or cannot be realized with fermionic two-orbital Hubbard models.As a consequence, our study should be considered qualitative, but it still provides a crude but valuable estimation of how large some extra terms beyond the canonical quadratic Heisenberg interactions can be.

A. Multi-Orbital Hubbard Model
For the exact-diagonalization calculations, we work with the multi-orbital Hubbard model mentioned in [6,7] and described as follows: where c † i,γ,σ (c i,γ,σ ) creates (annihilates) an electron at site i, with orbital γ, and spin projection along the z-axis σ.The first term represents the inter-and intra-orbital hopping between only nearest-neighbor sites.General hopping matrices for the two-and three-orbitals per site cases are displayed in Eqs.( 2) and (3), respectively, and in our study we allowed for the hoppings to vary over broad ranges to search for the largest ratios of Heisenberg interactions.The second term is the standard onsite Hubbard repulsion U between spins ↑ and ↓ electrons, at the same orbital.The third term contains the onsite inter-orbital repulsion, with the usual relation U = U −2J H due to rotational invariance.The fourth term involves the Hund's coupling J H that explicitly shows the ferromagnetic character between orbitals.The last term represents the onsite inter-orbital electron-pair hopping P i,γ = c i,γ,↑ c i,γ,↓ .All these terms in the Hubbard model are canonical.
The general hopping matrices used here for the exactdiagonalization calculation of two-and three-orbitals per site on the two-site system are: The allowed high-order Heisenberg model for any spin-S system can be written generically as: Using the above equation we can write the general Hamiltonian for S = 1 spin system as: Energy (E/J1) vs. J2/J1 for the two-site S = 1 Heisenberg model.The shaded area depicts the region with ordering singlet, triplet and quintuplet in increasing order of energies, as it occurs in the more fundamental two-orbital per site Hubbard model.
We diagonalize the Hamiltonian in Eq. ( 5) for the twosite system and obtain the following three energy levels: In Fig. 1, we illustrate the plot of these energy levels vs. J 2 /J 1 .For J 2 /J 1 < 1/3, the ordering of these levels strictly follows the singlet-triplet-quintuplet sequence in increasing order of energies.This is vital as the same sequence appears in the more fundamental two-orbital per site Hubbard model in strong coupling.
Of course, when comparing these energies mentioned in Eq. ( 6) with the Hubbard results obtained from exactdiagonalization in the strong coupling regime a constant offset in energies must be included, leading generically to E a = E a + E of f where a = s, t, q and E of f is the offset energy.Based on this information and the energies provided in Eq. ( 6) one can compute the ratio J 2 /J 1 in terms of the Hubbard energies obtained from exactdiagonalization E a 's as: The above equation is used for calculating the values of J 2 /J 1 in our two-site two-orbitals per site exactdiagonalization study, in the range where the Hubbard model energies are in the expected singlet-tripletquintuplet order, starting from the singlet ground state (this assumption tends to break down only in weak coupling, already outside the range of the Heisenberg model description).
Similarly for S = 3/2 the high-order Heisenberg Hamiltonian reads: We diagonalize this Hamiltonian in Eq. ( 8) for the twosite system and obtain four energy levels: Following the same reasoning as in the case of S = 1, i.e. considering an offset energy, then E a = E a + E of f where a = s, t, q, v, and using the set of equations provided in Eq. ( 9) the analytical expression for J 2 /J 1 and J 3 /J 1 in terms of E a 's for S = 3/2 becomes and Equations ( 10) and ( 11) were used for calculating the values of J 2 /J 1 and J 3 /J 1 in our two-site three-orbitals per site exact-diagonalization study, respectively.Here, we do not include a figure like Fig. 1 for the case of S = 3/2 because it would require a three-dimensional plot of energy vs. J 2 /J 1 and J 3 /J 1 which would be difficult to visualize.For this reason, we simply have included here the relevant equations that were employed.

IV. RESULTS
In this section, we will discuss our numerical results via exact-diagonalization for the two-site system.Note that not only U and J H will be varied, but the most time-consuming portion of the calculation arises from the large number of hopping amplitude ratios that were studied (using t 11 as unit of reference).As a consequence, we have analyzed hundreds of different ratios of Hamiltonian parameters and in all cases mapped the lowenergy results into the corresponding Heisenberg models.Specifically, on average we run over 30 values of U and 12 values of J H /U , for a fixed set of hopping amplitudes.This already amounts to 360 runs.For two orbitals per site, we used 36 combinations of t 22 /t 11 and t 12 /t 11 for a total of 360×36 = 12,960 cases.For three orbitals per site, we used 196 combinations of the ratios t 22 /t 11 , t 33 /t 11 , t 12 /t 11 , t 13 /t 11 , and t 23 /t 11 for a total of 360×196=70,560 cases.Crudely estimating, the total number of cases studied is of the order of 4×10 4 .We computationally automatized the fittings, and from the vast array of numbers we isolated approximately 150 sets of data containing the largest ratios for J 2 /J 1 and J 3 /J 1 .Those special cases were plotted and visually inspected.From that set, the very small subset displayed in the figures shown in the present publication is the subset that in our judgement best represents the cases where the Heisenberg coupling ratios are the largest in absolute value, because our primary aim is to establish upper bounds on those quantities.These ratios can be positive or negative.2. J2/J1 plotted vs. U for a two-orbitals per site system via exact-diagonalization, at the various ratios JH /U indicated.The hopping parameters chosen are t22 = t11 and t12 = t21 = 0 for this example, namely the unit matrix.To the readers, the hopping matrix is presented as inset in the plot.The bandwidth for this set of hopping parameters is W = 4t11.Note that "bandwidth" is defined with regards to the tight-binding model with the hoppings used here but in the bulk limit.3. J2/J1 plotted vs. U for a two-site two-orbital per site system via exact diagonalization, at the various ratios JH /U indicated.The hopping parameters chosen are t22 = t12 = t21 = t11.This hopping matrix is presented in the plot as inset.The bandwidth for these set of hopping parameters is W = 8t11 (for definition of bandwidth see caption of Fig. 2).
Here, we present our two-site two-orbitals per site exact-diagonalization results.All the results presented below have the same low-energy order: first a singlet (with the total spin S T ot = 0) for the ground state, then a triplet (S T ot = 1) for the first excited state, and finally a quintuplet (S T ot = 2) for the second excited state.
In both Fig. 2 and Fig. 3 we first performed exact diagonalization of the multi-orbital Hubbard model defined in Eq. ( 1).The hopping parameters used for each of these figures is shown as an inset, for better visualization, and also in the caption.As already explained, for each J H /U , we identified the range of U that gives us the ordering: singlet, triplet and quintuplet for the ground-state, first excited-state, and second excited state, respectively.After the proper regime of couplings was identified, the energies of these respective states were used to calculate the ratio J 2 /J 1 using Eq.(7).
Our main result is that the largest ratio observed (in absolute value) is close to 0.4.For a wide variety of "less symmetric" hopping amplitudes, namely employing neither the unit matrix or the matrix with all elements equal, we observed that |J 2 /J 1 | is smaller than those reported in Figs. 2 and 3. Two important details to remark are: (a) the ratios shown can be both positive and negative and for this reason the two examples shown were chosen.In both cases, positive and negative, the largest magnitudes of the ratios are not too different.(b) As obvious from the figures, the largest ratios are obtained as U is reduced from very strong coupling.This makes sense because in the limit where a perturbative expansion in t 11 /U is valid, J 1 is the lowest order and J 2 the next leading order.Naturally, their ratio of coefficients scales as t 11 /U and J 2 /J 1 converges to zero as U diverges.As a consequence, we can firmly conclude that the most promising region to observe the effects of the biquadratic term is the range of U/W ∼ 1, which is the intermediate coupling regime.This region of parameter space usually contains a variety of exotic phases because here several tendencies are in close competition leading to "frustration" effects which are hidden, namely not obvious at the Hamiltonian level.U for the two-site three-orbitals per site system obtained via exact diagonalization, at the several values of JH /U indicated.The hopping parameters chosen for this particular plot are t11 = 1.0, t22 = 0.25, t33 = 0, t12 = 0.75, t23 = 0.125 and t13 = 0.They are also shown in an inset in the figure for the benefit of the readers.The bandwidth for this particular set of hopping parameters is W ≈ 5.86t11 (for definition of bandwidth see caption of Fig. 2).The color convention for the ratio JH /U , as well as the hopping matrix, is common to both panels.

B. Two-Site Three-Orbitals per site
In this subsection, we present our two-site threeorbitals per site exact diagonalization results.All the results presented below have the same energy ordering, namely singlet (S T ot = 0) for the ground state, triplet (S T ot = 1) for the first-excited state, quintuplet (S T ot = 2) for the second-excited state, and septuplet (S T ot = 3) for the third-excited state.The latter origi-nates in the three orbital per site nature of the problem, and it does not appear for two orbitals per site.Namely, the extra spin manifold occurs because the total number of electrons in the system is 6 which allows total spins 3, 2, 1, and 0, contrary to a total of 4 electrons in the previous subsection.
Unlike the two-site two-orbital per site case, here we observe that it is the "less symmetric" (as mentioned in section IV A) hopping amplitudes that gives large values of the ratios |J Qualitatively, the conclusions in Fig. 4 resemble those for the two-orbital per site case.Once again, the ratios are the largest as U/t 11 decreases from the strong coupling regime.Thus, again we conclude that the intermediate coupling region U/W ∼ 1 is the most promising to observe sizable values for J 2 and J 3 .Also, the largest values of J 2 /J 1 are similar to those of the two-orbital per site case.However, as expected from the strong coupling expansion, J 3 /J 1 is an order of magnitude smaller than J 2 /J 1 because it requires the next order in the large U expansion to develop as compared with J 2 /J 1 .
Figure 5 illustrates the dependence of the results varying slightly the hopping amplitudes.Focusing on the matrices contained in both panels, the only difference between the two cases resides in the hopping t 23 , which varies by a factor 2. However, this relatively small modification leads to a reduction in approximately a factor two in the values of J 2 /J 1 .This high sensitivity to small changes in the hoppings is surprising.Such effect manifest the most at intermediate couplings, while in strong coupling the ratios are less sensitive to small hopping modifications.
In Figure 6, we illustrate the case where the hoppings reside only along the diagonal, but one of them, i.e. t 33 , is zero.Surprisingly, in this case the fits lead to negative values for both J 2 /J 1 and J 3 /J 1 .The strength is also reduced when compared with Fig. 5.In Fig. 7, J 2 /J 1 is shown now increasing t 33 from zero, as compared with Fig. 6.Here we see that as we increase t 33 the largest value of the ratio J 2 /J 1 decreases slowly, indicating that in order to find the maximum possible value of the ratio J 2 /J 1 the hopping amplitude t 33 must be zero.Similarly, we tune other hopping amplitudes and find the best possible scenario where we achieve the largest value of J 2 /J 1 and J 3 /J 1 .U for a two-site three-orbitals per site system using exact diagonalization, at the indicated several ratios of JH /U (color convention is the same for both panels).The hopping parameters chosen for panel (a) are t11 = t22 = 1.0, t33 = 0.125 and t αβ = 0 for all α = β and for panel (b) are t11 = t22 = 1.0, t33 = 0.25 and t αβ = 0 for all α = β.These hoppings are also shown in the insets.The bandwidth for both sets of hopping parameters is W = 4t11 (for definition of bandwidth see caption of Fig. 2).

V. DISCUSSION
In our study we have focused on a two-site electronic multi-orbital Hubbard model to deduce what range of Heisenberg couplings are possible at intermediate and large values of the interaction strength U .In particular, for two orbitals per site we focused on how large the biquadratic coupling strength J 2 can become, in terms of the canonical quadratic Heisenberg superexchange interaction J 1 .We found that J 2 /J 1 can be of both signs, but in magnitude cannot exceed ∼0.4.This allows for the model used by Affleck and collaborators [9] to be realized employing electronic models.It would be interesting to investigate if this electronic model -namely selecting suitable Hubbard U , Hund coupling J H , and hoppings such that J 2 /J 1 = 1/3 -will also lead to a valence bond ground state, although likely the said electronic model will not be exactly solvable.On the other hand, the exactly solvable case J 2 /J 1 = 1 cannot be realized with the model we used.For the case of spin S = 3/2, namely using three orbitals per site, the conclusions are similar: once again J 2 /J 1 cannot exceed ∼0.4,while J 3 /J 1 is even smaller by a factor approximately two.
Our study also suggests that some spin-only models that were studied to search for quantum spin liquids must impose constraints on the parameter space explored.To realize spin liquids using electronic models the most optimal path continues being the addition of hoppings beyond nearest-neighbors to create explicit frustration.
Note that our conclusions are in agreement with calculations using a two-orbital per site Hubbard model [35], carried out perturbatively at small t/U up to fourth order (first and third order cancel; the second order gives the canonical quadratic Heisenberg model, and the fourth order is the one relevant for our discussion involving the biquadratic contribution).By this fairly different procedure, nevertheless a conclusion similar to ours was reached: the ratio J 2 /J 1 is severely limited at large U .Our results, valid at any value of U because they do not rely on perturbation theory, suggest that intermediate U is more promising than strong U , but still J 2 /J 1 cannot reach values above 0.4.The methodology proposed in Ref. [36], adding to the problem an extra orbital residing in a neighboring site, may reduce J 1 , providing a promising path to enhance the ratio J 2 /J 1 [37].Our next goal is to investigate the effect of spin-orbit coupling in these two-site multi-orbital Hubbard systems [38] and estimate the range of biquadratic-quadratic Heisenberg couplings.
0.25 0.50 0.75 1.00 1.25 FIG.2.J2/J1 plotted vs. U for a two-orbitals per site system via exact-diagonalization, at the various ratios JH /U indicated.The hopping parameters chosen are t22 = t11 and t12 = t21 = 0 for this example, namely the unit matrix.To the readers, the hopping matrix is presented as inset in the plot.The bandwidth for this set of hopping parameters is W = 4t11.Note that "bandwidth" is defined with regards to the tight-binding model with the hoppings used here but in the bulk limit.
FIG.3.J2/J1 plotted vs. U for a two-site two-orbital per site system via exact diagonalization, at the various ratios JH /U indicated.The hopping parameters chosen are t22 = t12 = t21 = t11.This hopping matrix is presented in the plot as inset.The bandwidth for these set of hopping parameters is W = 8t11 (for definition of bandwidth see caption of Fig.2).

FIG. 4 .
FIG. 4. (a)J2/J1 and (b) J3/J1 vs. U for the two-site three-orbitals per site system obtained via exact diagonalization, at the several values of JH /U indicated.The hopping parameters chosen for this particular plot are t11 = 1.0, t22 = 0.25, t33 = 0, t12 = 0.75, t23 = 0.125 and t13 = 0.They are also shown in an inset in the figure for the benefit of the readers.The bandwidth for this particular set of hopping parameters is W ≈ 5.86t11 (for definition of bandwidth see caption of Fig.2).The color convention for the ratio JH /U , as well as the hopping matrix, is common to both panels.

FIG. 6 .
FIG.6.(a) J2/J1 and (b) J3/J1 plotted vs. U for two-site three-orbitals per site system obtained via exact diagonalization, at the several values of JH /U indicated.Color convention is common to both panels.The hopping parameters chosen for the plot are t11 = t22 = 1.0, t33 = 0 and t αβ = 0 for all α = β.The bandwidth for these set of hopping parameters is W = 4t11 (for definition of bandwidth see caption of Fig.2).
t 11 t 12 t 13 t 21 t 22 t 23 t 31 t 32 t 33 where t αβ represents the nearest-neighbor hopping element from orbital α to orbital β.Due to rotational symmetry of the two-site system, t αβ = t βα .This reduces the number of hopping elements from N 2 o to N o (N o + 1)/2, where N o is the number of orbitals.B.Heisenberg Model with Higher Order Terms