Effects of anisotropy in spin molecular-orbital coupling on effective spin models of trinuclear organometallic complexes

TL;DR

This paper investigates how anisotropy in spin molecular-orbital coupling influences effective spin models in layered trinuclear organometallic complexes, revealing phase transitions and anisotropic magnetic behaviors.

Contribution

It introduces a detailed effective spin model incorporating anisotropic SMOC effects and predicts phase transitions in these complex materials.

Findings

Effective XXZ + 120° honeycomb quantum compass model derived

Strong anisotropies lead to a transition from Haldane to D-phase

External magnetic fields induce a transition from Haldane to magnetic order

Abstract

We consider layered decorated honeycomb lattices at two-thirds filling, as realized in some trinuclear organometallic complexes. Localized $S=1$ moments with a single-spin anisotropy emerge from the interplay of Coulomb repulsion and spin molecular-orbit coupling (SMOC). Magnetic anisotropies with bond dependent exchange couplings occur in the honeycomb layers when the direct intracluster exchange and the spin molecular-orbital coupling are both present. We find that the effective spin exchange model within the layers is an XXZ + 120$^\circ$ honeycomb quantum compass model. The intrinsic non-spherical symmetry of the multinuclear complexes leads to very different transverse and longitudinal spin molecular-orbital couplings, which greatly enhances the single-spin and exchange coupling anisotropies. The interlayer coupling is described by a XXZ model with anisotropic biquadratic terms. As the correlation strength increases the systems becomes increasingly one-dimensional. Thus, if the ratio of SMOC to the interlayer hopping is small this stabilizes the Haldane phase. However, as the ratio increases there is a quantum phase transition to the topologically trivial `$D$-phase'. We also predict a quantum phase transition from a Haldane phase to a magnetically ordered phase at sufficiently strong external magnetic fields.

Figures

Tables

Table 1. Table 1: List of parameters entering our microscopic model for Mo 3 S 7 (dmit) 3 . The exchange couplings of our derived effective spin exchange model ( 37 ) using the actual DFT parameters jacko2017 obtained for the crystal are also tabulated. The exchange couplings are isotropic so α 𝛼 \alpha can be x , y , z 𝑥 𝑦 𝑧 x,y,z . Parameters of the effective model that are smaller than 10 − 4 superscript 10 4 10^{-4} are not included. All energy units are in eV.

tc	t	tz	${\lambda }_{xy}$	${\lambda }_z$	Jab	Jc
0.06	0.047	0.041	0.0025	0.005	0.0024	0.01296

Equations

$$H_{\text{eff}}$$

$$H=H_0+H_{SMOC}+H_{U-J_F}.$$

$$H_0=-t_c\sum _{⟨ij⟩\sigma }\left(a_{i\sigma }^{†}{a}_{j\sigma }+H.c.\right),$$

$$H_{SOC}=\mathbf{K}\cdot \mathbf{S},$$

$$\mathbf{K}=\frac{\hslash }{4m^2{c}^2}[\mathbf{p}\times \nabla V(\mathbf{r})].$$

$$H_{SMOC}$$

$$H_{U-J_F}$$

$$H_0=\sum _{k\sigma }{ϵ}_k{b}_{k\sigma }^{†}{b}_{k\sigma },$$

$$b_{k\sigma }^{†}=\frac{1}{\sqrt{3}}\sum _{j=1}^3{e}^{ik\varphi (j-1)}{a}_{j\sigma }^{†},$$

$$c_{k\sigma }^{†}={\text{sgn}}^k(-k)\frac{1}{\sqrt{3}}\sum _{j=1}^3{e}^{ik\varphi (j-1)}{a}_{j\sigma }^{†}.$$

$$c_{0,\sigma }^{†}=b_{0\sigma }^{†}$$

$$c_{1\sigma }^{†}=-b_{1,\sigma }^{†}$$

$$c_{-1\sigma }^{†}=b_{-1,\sigma }^{†}.$$

$$H_0$$

$$L^{+}$$

$$L^{-}$$

$$L_z$$

$$H_{SMOC}=\frac{{\lambda }_{xy}}{\sqrt{2}}(b_{0\downarrow }^{†}{b}_{-1\uparrow }-b_{1\downarrow }^{†}{b}_{0\uparrow }-b_{0\uparrow }^{†}{b}_{1\downarrow }+b_{-1\uparrow }^{†}{b}_{0\downarrow })+\frac{{\lambda }_z}{2}(b_{1\uparrow }^{†}{b}_{1\uparrow }-b_{1\downarrow }^{†}{b}_{1\downarrow }-b_{-1\uparrow }^{†}{b}_{-1\uparrow }+b_{-1\downarrow }^{†}{b}_{-1\downarrow }).$$

$$H$$

$$H_{U-J_F}$$

$$H_m^{\text{eff}}=D({\mathcal{S}}_{{\mathbf{r}}_m}^z{)}^2.$$

$$H=H_{\ell }+H_m+H_{kin},$$

$$H_{kin}^{dumbbell}=-t\sum _{\sigma }\left(a_{\ell 1\sigma }^{†}{a}_{m1\sigma }+a_{m1\sigma }^{†}{a}_{\ell 1\sigma }\right),$$

$$H_{kin}^{dumbbell}=-\frac{t}{3}\sum _{k_1{k}_2\sigma }\left(b_{\ell k_1\sigma }^{†}{b}_{m{k}_2\sigma }+b_{m{k}_2\sigma }^{†}{b}_{\ell k_1\sigma }\right),$$

$$H_{kin}^{tube}=-t_z\sum _{i\sigma }\left(a_{\ell i\sigma }^{†}{a}_{mi\sigma }+a_{mi\sigma }^{†}{a}_{\ell i\sigma }\right).$$

$$H_{kin}^{tube}=-t_z\sum _{k\sigma }\left(b_{\ell k\sigma }^{†}{b}_{mk\sigma }+b_{mk\sigma }^{†}{b}_{\ell k\sigma }\right),$$

$$H_{eff}^{(2),conf}$$

$$|m_0⟩$$

$$H_{eff}^{(2),conf}$$

$$\overline{H}(\epsilon )$$

$$P_{\mu }H{P}_{\nu }(1-{\delta }_{\mu \nu })+i{P}_{\mu }H{P}_{\mu }(P_{\mu }S{P}_{\nu })-i(P_{\mu }S{P}_{\nu })P_{\nu }H{P}_{\nu }=0.$$

$$i{P}_{\mu }H{P}_{\nu }$$

Full text

I Introduction

The interplay of strong Coulomb interaction and spin-orbit coupling (SOC) can lead to emergent quantum phases balents2014 and new phenomena which remain poorly understood. The conventional Mott transition can be strongly affected by SOC leading to a topological Mott insulator with bulk charge gap but fractionalized surface states carrying spin but no charge pesin2010 . Such states may be realized in Ir-based transition metal oxides such as Sr2IrO4. In contrast to conventional Mott insulators, the spin exchange interactions arising in Mott insulators with SOC are typically anisotropic with quantum compass nussinov2015 instead of the conventional Heisenberg interactions. A possible realization of a quantum compass model on a hexagonal lattice, i.e., a Heisenberg-Kitaev model,jackeli2009 ; perkins2014a ; perkins2014b may be found in Na2IrO3 and Li2IrO3 materials in which SOC removes the orbital degeneracy of the 5d electrons leading to effective $S=1/2$ pseudospins. Interestingly, the Kitaev model is exactly solvable: it sustains a spin liquid ground state whose low energy excitations are Majorana fermions kitaev2006 . In other iridates with strong SOC such as Sr2IrO4, an antisymmetric Dzyaloshinski-Moriya (DM) interaction arises associated with the lack of an inversion symmetry center.

There are several strongly correlated molecular materials in which spin-orbit coupling is relevant including, metal-organic frameworks,oshikawa2016 layered organic salts,winter2017 ; elise2017 and multinuclear coordinated organometallic complexes.hoffman ; khosla2017 ; jacko2017 ; merino2016 ; powell2016 The elementary building blocks of multinuclear complexes are molecular clusters containing transition metal ions whose d-orbitals are hybridized with molecular orbitals where each of the hybrids is typically described by a single Wannier orbitaljacko2015 . The coupling of the spin with the electron currents around the Wannier orbitals describing each molecule gives rise to a spin molecular-orbital coupling (SMOC) khosla2017 ; jacko2017 .

A typical multinuclear complex is Mo3S7(dmit)3. Here the honeycomb networks of Mo3S7(dmit)3 molecules are stacked on top of each other along the $c$-direction of the crystal. Mo3S7(dmit)3 molecules can be described by three Wannier orbitals jacko2015 , and their packing on a honeycomb lattice within the layers leads to a decorated honeycomb lattices, as shown in Fig. 1. The electronic and magnetic properties of the decorated honeycomb lattice are interesting both in the weakly and strongly interacting limit. At weak coupling, a tight-binding model on such a lattice leads to topological insulating phases when SOC is turned on which display the quantum spin hall effectruegg2010 as predicted in graphenekane2005 . At strong coupling, the exact ground state of the Kitaev model on the decorated honeycomb lattice kivelson2007 , is a chiral spin liquid. Therefore, it is interesting to find possible realizations of the decorated honeycomb lattice in actual materials to probe such rich physics. Furthermore, similar models arise naturally in a number of other organicelise2017 and organometallic materialshoffman ; Silveira ; Wang1 ; Wang2 ; Wang3 and inorganic compounds with decorated lattices.Sheckelton ; Bao ; Chen

Organometallic complexes have intrinsic structural properties which make them different to transition metal oxides. A crucial difference comes from the fact that isolated molecular clusters break the spherical symmetry present in isolated transition metal ions. While the total angular momentum of the ion is conserved, it is only the component perpendicular to the molecular plane that is conserved in cyclic molecular clusters. Hence, in these systems, anisotropies are intrinsic to the molecules constituting the material, whereas in transition metal oxides anisotropies can only be achieved via the environment surrounding the ions in the crystal. This suggests that anisotropic spin exchange interactions may be easily generated in organometallic complexes due their intrinsic structure. These anisotropies may be further enhanced by the anisotropic SMOC typically found in these systems. jacko2017 ; powell2016 SMOC is an emergent coupling between electron currents circulating around the cyclic molecules and the electron spin. Also by tuning the relative orientation between molecules in the crystal a Dzyaloshinskii-Moriya interaction can be generatedpowell2016 . All the above suggests that these materials are ideal playgrounds for the realization of quantum compass models.nussinov2015

Recently merino2016 ; powell2016 we derived an effective super-exchange Hamiltonian that captures the magnetic properties of trinuclear coordinated complexes at strong coupling. The onsite Coulomb repulsion, $U$, leads to $S=1$ moments localized at each triangular cluster whence SMOC, $\lambda$, induces a single-spin anisotropy, $D$. The $S=1$ moments behave as weakly coupled chains due to the decorated lattice structure of trinuclear organometallic complexes. merino2016 ; powell2016 The lattice structure is such that three hopping amplitudes connect two nearest-neighbor molecules along the $c$-direction while only one hopping amplitude connects nearest-neighbor molecules in the $a$-$b$ planes, cf. Fig. 2. As $U$ is increased exchange of electrons between nearest-neighbor molecules in the $a$-$b$ plane is suppressed as compared to exchange between molecules along the $c$-direction. This leads to a quasi-one-dimensional effective spin exchange model of $S=1$ localized moments which is in the Haldane phase.powell2016

Here we extend our previous work, which focused on Mo3S7(dmit)3, by studying the more general problem of trinuclear organometallic complexes with strong correlations and strong anisotropic SMOC. After introducing our general combined analytical and numerical approach to extract exchange coupling parameters in these systems we show how anisotropy in SMOC plays a crucial role in determining the level of anisotropy of the effective spin Hamiltonian. We show that the effective spin exchange Hamiltonian for two-thirds filled trinuclear coordination crystals is

[TABLE]

where ${\mathcal{S}}^{\alpha}_{\bm{r}}$ is the $\alpha$th component ($\alpha=x,y,z$) of the pseudospin-one operator at position $\bm{r}$, $\bm{r}_{\ell}$ is the position of site $\ell$, $\bm{\delta}_{z}=(0,0,c)$, c is the interlayer spacing, $j=1,2,3$ labels the nearest neighbour bonds as marked in Fig. 1, $\phi_{j}={2\pi(j-1)/3}$, $\bm{\delta}_{j}=(\sin\phi_{j},\cos\phi_{j},0)a_{g}$ is the vector, of length $a_{g}$, connecting one sublattice to its three nearest neighbours in the plane, and $\sum_{\ell\in\bigtriangledown}$ indicates that the sum runs over only the sublattice of triangles that point down in Fig. 1.

For large $U$, the magnitude of the antiferromagnetic exchange coupling between nearest neighbor clusters in the $c$-direction, $J^{c}$, is much larger than the exchange coupling between nearest-neighbor clusters in the $a$-$b$ plane, $J^{ab}$, we conclude that the magnetic properties of two-thirds filled trinuclear coordination crystals can be effectively described by $S=1$ XXZ chains with a local single-spin anisotropy, $D^{*}$ and anisotropic biquadratic terms, $P_{\alpha\beta}$. We explore the effect of anisotropic SMOC, $\lambda_{xy}\neq\lambda_{z}$, finding that the largest anisotropic spin exchange couplings and single-spin anisotropies emerge when $\lambda_{xy}/\lambda_{z}<1$, which is the relevant parameter regime for Mo3S7(dmit)3.

For Mo3S7(dmit)3 ab initio estimates of SMOCjacko2017 indicate that $\lambda_{xy}\approx 0.042t_{c}$, and $\lambda_{xy}\approx\lambda_{z}/2$. This, suggests that single-spin anisotropies are smaller than the exchange coupling along the $c$-direction, $D^{*}<J^{c}$, so that Mo3S7(dmit)3 is in the Haldane phase rather than in the topologically trivial ‘$D$-phase’, i.e., the tensor product of the $j=0$ singlets (where $j$ is the $z$-component of the total angular momentum) at each cluster, which is expected for $D^{*}>J^{c}$. In spite of the small SMOC values found in Mo3S7(dmit)3 (see Fig. 3), the chemical flexibility of molecular crystals can significantly enhance $\lambda_{xy}$ and $\lambda_{z}$, and suppress $t_{z}$. Together this could drive other related systems into the $D$-phase and enhance anisotropies in the exchange interactions.

In Fig. 3 we show how the critical SMOC, $\lambda^{critical}_{xy}$, at which the transition from the Haldane to the $D$-phase occurs i. e. when $D^{*}(\lambda_{xy})\sim J^{c}$, is strongly suppressed by reducing $t_{z}$ and/or by a ferromagnetic intracluster exchange, $-J_{F}$. Variations in the SMOC anisotropy (not shown) can also significantly vary $D^{*}$ [see Appendix B and particularly, Eq. (51)]. On the other hand, increasing $U$ by, say, a factor of two does not change $\lambda^{critical}_{xy}$ since $J^{c}$ is moderately influenced by $U$ when $U\rightarrow\infty$. Intracluster charge fluctuations not captured by our spin model but present in the original Hubbard model are found to strongly suppress the spin gap nourse2016 . For the microscopic parameters found from density functional theory (DFT) jacko2015 ; jacko2017 for Mo3S7(dmit)3 the transition line is given by $D^{*}\sim 0.066J^{c}$. The charge fluctuation effect suppresses $\lambda^{critical}_{xy}$ even further becoming comparable to the SMOC in Mo3S7(dmit)3 crystals. Hence, even though SMOC is small in Mo3S7(dmit)3 it may be possible to drive it from the Haldane to the $D$-phase by modifying crystal parameters, in particular, by suppressing $t_{z}$. This may be achieved by applying negative uniaxial pressure along the $c$-direction of the crystal which increases the interlayer distance. Alternatively, an expansion along the $c$-direction can be achieved by applying uniaxial (positive) pressure on the $a-b$ directions through the Poisson effect. However, this procedure can lead to changes in the in-plane arrangement of the molecules distorting the physics of the honeycomb lattice discussed here.

We analyze the possible magnetic anisotropies arising in the decorated honeycomb lattice of Fig. 1, which can be realized by isolating the $a$-$b$ planes of trinuclear clusters. More specifically, we analyze the role played by the interplay of Coulomb repulsion, intracluster exchange and SMOC in producing anisotropic exchange couplings. We study the role played by SMOC anisotropy, $\lambda_{xy}\neq\lambda_{z}$, which is generically the case in these systems and has not been considered in previously. We find that the effective exchange couplings within the $a$-$b$ planes are anisotropic only when both SMOC and intracluster exchange, $J_{F}$, are present. These magnetic anisotropies lead to a spin-one XXZ + 120∘ degree honeycomb quantum compass model with single spin anisotropy. In the limit of $J_{F}\rightarrow 0$, our effective spin exchange model reduces to the conventional isotropic $S=1$ antiferromagnetic Heisenberg model on a honeycomb lattice.

We predict that under a sufficiently large external magnetic field, the Haldane phase can be destroyed giving way to a three-dimensional ordered magnet. This occurs at a critical magnetic field, $h_{c}\sim\Delta_{s}$, where $\Delta_{s}$ is the zero-field Haldane gap of the $S=1$ chain.

The present paper is organized as follows. In Section II we introduce the minimal strongly correlated model for describing the electronic properties of isolated triangular molecules in the presence of SMOC. The physics of a single molecule described by this model is discussed in the Appendix A. In Section III we analyze the electronic structure of two coupled trimers arranged as two nearest-neighbor molecules in the $a$-$b$ plane and also as two nearest-neighbor molecules along the $c$-direction. The energy level spectra of two coupled trimers is obtained exactly and compared to second order perturbation theory. In Section IV the combination of the numerical perturbative approach with an analytical canonical transformation (see also Appendix B), used to extract the exchange interactions between the nearest neighbor pseudospins, is detailed. In Section V, we discuss the qualitative phase diagram expected for the quasi-one-dimensional spin model arising from our approach. Finally, in Section VI, we conclude providing an outlook of our work.

II Model of isolated trimers in the presence of SMOC

Here we condisder crystals formed of triangular tri-nuclear molecules. In order to understand the effects of SMOC on the electronic and magnetic properties of these systems we first discuss the relevant model for isolated triangular clusters. The simplest strongly correlated model is a Hubbard model on a triangle in the presence of SMOC janani2014a ; janani2014b ; jacko2015 :

[TABLE]

In general all operators should also have a molecular label but this is suppressed throughout the current Section as we deal only with a single complex.

The tight-binding part reads

[TABLE]

where $t_{c}$ is the hopping between the hybrid metal-ligand orbitals at nearest-neighbor sites in the cluster and $a^{\dagger}_{i\sigma}$ creates an electron at the $i$th Wannier orbital with spin $\sigma$.

The general SOC contribution iskhosla2017

[TABLE]

where ${\bf S}$ is the electron spin and ${\bf K}$ is a pseudovector operator,

[TABLE]

We project onto a basis of one Wannier orbital per site of the model illustrated in Fig. 1. The two spin states of the Wannier orbital are a Kramers pair thus this projection removes all non-trivial effects of the atomic SOC. For example, in Mo3S7(dmit)3 the Mo atoms are in a $C_{1}$ environment. Consider an atom with ${\bf L}\cdot{\bf S}$ atomic SOC in a $C_{1}$ environment with time reversal symmetry. The most general coupling between two states is ${\bf B^{*}}\cdot{\bf S}+C^{*}\mathbb{1}$ (the most general $2\times 2$ Hamiltonian). However, we require that these states remain degenerate to maintain time-reversal symmetry. Thus only the $C^{*}$ term remains, providing a constant energy shift as the only effect of atomic SOC in the subspace of the Krammers pair. Note that the ${\bf B^{*}}$ term is a Zeeman splitting term; if we had projected out more states, this term could be non-zero. It has been argued that this is relevant to some transition metal oxides khaliullin05 ; kim2012 where this projection induces an effective anisotropy on the atomic SOC. Thus the only SOC term possible in our model is the direct coupling of the spin to currents running around the plane of the molecule (SMOC), which has no analogue in transition metal oxides.

For $C_{3}$ symmetric molecules it can be shownkhosla2017 the SMOC is

[TABLE]

where $L$ is the molecular orbital angular momentum of electrons in the cluster, $\lambda_{xy}$ describes the transverse SMOC while $\lambda_{z}$ describes the longitudinal contribution.

Finally, the Hubbard-Heisenberg term reads

[TABLE]

where $U$ is the onsite Hubbard interaction, $J_{F}$ is an intracluster exchange interaction, and $n_{i\sigma}=a^{\dagger}_{i\sigma}a_{i\sigma}$ the number operator. The direct exchange, $J_{F}$, between electrons at nearest-neighbor sites is generically non-zero and favors ferromagnetic tendencies, i.e., it is expected to be negative, $J_{F}<0$. We will see below that, even if it is much smaller than the direct Coulomb interaction, $J_{F}$ plays a crucial role in generating magnetic anisotropies. It plays a similar role as the Hunds coupling in transition metal oxidesjackeli2009 , which also generates magnetic exchange anisotropies between spins in the lattice.

The non-interacting part (3) can be readily diagonalized:

[TABLE]

using Bloch operators:

[TABLE]

with $\phi={2\pi/3}$. $k=0,\pm 1$ correspond to the allowed $0,\pm{\frac{2\pi}{3}}$ momenta in the first Brillouin zone of the triangular cluster with energies $\epsilon_{0}=-2t_{c}$ and $\epsilon_{1}=\epsilon_{-1}=t_{c}$.

The SMOC contribution to $H$ is most naturally described using ‘Condon-Shortley’ states which are eigenstates of the $z$-component of the angular momentum, $L_{z}$, of the clusterkhosla2017 ; powell2015 ,

[TABLE]

More explicitly we have

[TABLE]

Note that as Bloch’s theorem applies to the cluster, the $z$-component of angular momentum is defined up to $3n$ (in units of $\phi$) with $n$ an integer, i.e., Bloch states with momentum $k^{\prime}$ satisfying $k=k^{\prime}\pm 3n$ are equivalent to the $k=0,\pm 1$ states.

Hence, the tight-binding part of the Hamiltonian $H_{0}$ can be expressed either in terms of the Condon-Shortley or Bloch operators as

[TABLE]

Similarly, from the expressions of the angular momentum in terms of the Bloch states:

[TABLE]

the SMOC contribution to the Hamiltonian of the isolated cluster reads:

[TABLE]

We may also express $H$ in the site basis, ${|i\sigma\rangle}$, using the transformation of Eq. (9) which leads to:

[TABLE]

with $A={(e^{i\phi}-1)\over 6}$, $B={i\over 3}\sin(\phi)$, and $\sigma=\pm 1$. It is evident from the above Hamiltonian that SMOC can be understood as a spin-dependent hopping between nearest-neighbor sites of the trimers.

Four-component relativistic ab initio calculations jacko2017 for Mo3S7(dmit)3 have found anisotropic SMOC: $\lambda_{xy}\approx\lambda_{z}/2>0$, cf. Table 1. Below we will fix $t_{c}>0$ as the unit of energy and explore different values of SMOC and different $\lambda_{xy}/\lambda_{z}$ ratios. Note that the electronic properties of the model are invariant under the particle-hole transformation $a^{\dagger}_{i}\rightarrow h_{i},a_{i}\rightarrow h^{\dagger}_{i}$, where $h^{\dagger}_{i}$ and $h_{i}$ are hole operators together with the transformation $t_{c}\rightarrow-t_{c}$, $\lambda_{xy}\rightarrow-\lambda_{xy}$, $\lambda_{z}\rightarrow-\lambda_{z}$. The onsite Coulomb repulsion within each Wannier orbital, $U$, is comparable to or even larger than the bandwidth of the relevant Mo3S7(dmit)3 bands crossing the Fermi energy. We will assume $U=10t_{c}$ as a reasonable estimate. Since the Mo3S7(dmit)3 crystal is at $2/3$-filling there are $N=4$ electrons per triangular cluster in the crystal. In order to fully characterize the electronic structure of two coupled clusters through perturbation theory techniques we have analyzed triangular clusters with $N=3,4,5$ electrons and the parameters $t_{c},\lambda_{xy},\lambda_{z}>0$, relevant to Mo3S7(dmit)3 crystals. Through the particle-hole transformation we can also obtain the electronic structure of triangular clusters with $N=1$ ($N=2$) electrons from the $N=5$ ($N=4$) solutions by switching the sign of $\lambda_{xy},\lambda_{z},t_{c}$.

Since $J_{z}=L_{z}+S_{z}$ is a conserved quantity: $[J_{z},H]=0$, it is convenient to use the $(k,\sigma)$ representation instead of the site representation to classify the basis states according to their quantum number: $j=k+\sigma$. We have already expressed $H_{0}+H_{SMOC}$ in the $(k,\sigma)$ basis through Eq. (12)-(14). The Hubbard-Heisenberg contribution is expressed in the $(k,\sigma)$ basis as

[TABLE]

For the triangular clusters studied here ${1\over 3}\sum_{k,k^{\prime},k\neq k^{\prime}}(U-J_{F}\cos((k-k^{\prime})\phi)-J_{F})n_{k\uparrow}n_{k^{\prime}\downarrow}=(U/3-J_{F}/6)\sum_{k,k^{\prime},k\neq k^{\prime}}n_{k\uparrow}n_{k^{\prime}\downarrow}$. Note that while for the Hubbard-Heisenberg model the effective Coulomb repulsion between electrons is different for electrons in different orbitals, in a pure Hubbard model ($J_{F}=0$), all Coulomb interactions are equal to $U/3$. This has been shown to be important for finding spin exchange anisotropies in the context of transition metal oxidesaharony ; perkins2014a .

Hence, the full Hamiltonian can be explicitly expressed in the $(k,\sigma)$ basis using the expressions for $H_{0}$, $H_{SMOC}$ and $H_{U-J_{F}}$ in Eq. (12), (14) and (16), respectively.

In Appendix A we present results for the electronic structure of trimers with $N=3,4,5$ electrons expressed in this basis. From this analysis, we conclude that isolated trimers with $N=4$ electrons in the presence of SMOC effectively behave as pseudospin-one localized moments. In Fig. 4 we show that under SMOC the lowest energy triplet splits into a non-degenerate singlet ($j=0$) and a doublet ($j=\pm 1$), where $j$ is the $z$-component of total angular momentum. Higher energy excitations are doublets or non-degenerate under SMOC. Note that since we have an even number of electrons in the cluster, Kramers theorem does not apply and non-degenerate states are possible. Hence, SMOC induces a single-spin anisotropy at each cluster so that the effective spin model for $N=4$ electrons in the $m$-th Mo3S7(dmit)3 molecule in the crystal reads:

[TABLE]

As shown in Fig. 4 the overall energy level structure of the cluster i.e., level splittings and degeneracies remain unaffected by anisotropies in SMOC, $\lambda_{xy}\neq\lambda_{z}$ and/or intracluster exchange $J_{F}\neq 0$. However, the absolute value of $D$ is strongly enhanced when $\lambda_{xy}/\lambda_{z}<1$ as shown in Fig. 5. This is directly relevant to Mo3S7(dmit)3 crystals in which $\lambda_{xy}/\lambda_{z}\approx 1/2$.

III Two coupled triangular clusters

We now consider two triangular coupled clusters. We analyze the electronic structure of two nearest neighbor triangular clusters as arranged in Mo3S7(dmit)3 crystals and shown in Fig. 2. In Fig. 2(a) we show two nearest-neighbor clusters in the $a$-$b$ plane, whereas in 2(b) we show two nearest-neighbor clusters along the $c$-direction. The molecules in the “dumbbell” configuration of Fig. 2(a) are related by inversion symmetry as in Mo3S7(dmit)3. Molecules in the “tube” configuration of Fig. 2(b) are related by a rigid translation along the $c$-axis but no inversion symmetry is present. We first report exact results for the energy level structure. This gives key information about the type of spin exchange acting between the effective pseudospins localized at each trimer. These exact results are also used to benchmark perturbation theory calculations discussed in Section III.2.

III.1 Electronic structure

Consider a model of two trimers, $\ell$ and $m$ coupled by $H_{kin}$:

[TABLE]

where $H_{\ell}$ is the Hubbard-Heisenberg model of an isolated trimer, $\ell$, in the presence of SMOC as introduced previously, Eq. (15) in Sec. IIB. The hopping between two neighbor clusters is described through, $H_{kin}$.

As shown in Fig. 2a, in the coplanar dumbbell arrangement, there is only one hopping amplitude connecting the trimers, so $H_{kin}$ reads

[TABLE]

which connects, say, site 1 of the $\ell$-cluster with site 1 of the $m$-cluster. Here $a^{(\dagger)}_{mi\sigma}$ annihilates (creates) an electron with spin $\sigma$ in the $i$th Wannier orbital on molecule $m$. The kinetic energy contains off-diagonal hopping matrix elements in the Bloch basis:

[TABLE]

showing that the orbital momentum is not conserved in this case due to the breaking of trigonal symmetry.

In the tube arrangement, Fig. 2b the three vertices of the two clusters are connected by a hopping, $t_{z}$, and $H_{kin}$, reads

[TABLE]

As the tubes respect the trigonal symmetry of the isolated trimers, the angular momentum about the $C_{3}$ axis is conserved. Hence, the kinetic energy between two trimers in the tube arrangement is diagonal when expressed in the Bloch basis:

[TABLE]

where $k=0,\pm 1$, are the allowed momenta at each trimer. of isolated trimers.

We have exactly diagonalized model (18) for two coupled triangular clusters in the presence of SMOC. We consider the case in which each cluster is filled with $N=4$ electrons which is the relevant case for Mo3S7(dmit)3 crystals. In Fig. 6(a) and (b) we show the dependence of the eigenenergies, $E_{i}$, on $\lambda=\lambda_{xy}=\lambda_{z}$ (isotropic SMOC) for $U=10t_{c}$, $t=0.785t_{c}$ and $J_{F}=0$ in the dumbbell (a) and tube (b) arrangements. For $\lambda=0$ we find that the eigenspectrum of the coupled trimers consists of a ground state non-degenerate singlet, a triplet and a pentuplet. This is the eigenspectrum expected for an isotropic antiferromagnetic exchange interaction between two localized $S=1$ moments janani2014b . As $\lambda$ is increased the energy levels are split partially removing $\lambda=0$ degeneracies. The ground state of the coupled trimers is found to be non-degenerate for any value of $\lambda$.

In Figs. 6(c) and (d) we show the dependence of $E_{i}$ on $t$ for fixed SMOC, $\lambda=0.25t$ and $\lambda=t$. In both cases the eigenenergies depend quadratically on $t$, $E_{i}\propto t^{2}$ up to large values of $t/t_{c}\sim 1$ indicating that second order perturbation theory ($O(t^{2})$) is reliable. Below, we will further analyze the accuracy of the $O(t^{2})$ calculation for the model parameters that are relevant to Mo3S7(dmit)3 crystals.

In order to understand these spectra, it is important to understand the symmetries of the models. This can be a little subtle when SMOC is included.

In the absence of SMOC the dumbbell model is $D_{2h}$ symmetric as it also contains three mutually perpendicular two-fold rotation axes (cf. Fig. 2a). If two molecules, $\ell$ and $m$ are related to one another by inversion symmetry then the pseudovectorial nature of angular momenta requires that the SMOC is equal on both molecules: $\lambda_{\ell,xy}=\lambda_{m,xy}$ and $\lambda_{\ell,z}=\lambda_{m,z}$. On the other hand if two molecules are related by a $\pi$-rotation about, say, the $z$-axis this yields $\lambda_{\ell,z}=\lambda_{m,z}$, but $\lambda_{\ell,xy}=-\lambda_{m,xy}$. This leads to significant changes in the effective interactions between the molecular spins, which we have discussed elsewhere.merino2016 ; powell2016 Thus the case $\lambda_{\ell,xy}=\lambda_{m,xy}$ and $\lambda_{\ell,z}=\lambda_{m,z}$, which we consider here, lowers the symmetry to $C_{i}$ (triclinic).

In the absence of SMOC the tube model is D3h symmetric. This is lowered to C3v in the presence of SMOC, which can be understood as follows. In our model $\lambda_{\ell,xy}=\lambda_{m,xy}$ and $\lambda_{\ell,z}=\lambda_{m,z}$. Under a mirror reflection with respect to a plane perpendicular to the $z$-axis passing through the middle of the tube, i.e., a $\sigma_{h}$ operation, there is a change in sign of the transverse SMOC contribution: $\lambda_{\ell,xy}=-\lambda_{m,xy}$, which would be inconsistent with our model, except for $\lambda_{\ell,xy}=0$. We find that our model is symmetric under $C_{3}$ rotations and only has three $\sigma_{v}$ reflection planes. Hence, we conclude that the point group symmetry for the tube in the presence of SMOC is C3v.

In both the dumbbell (Fig. 6a) and tube (Fig. 6b) configurations the $\lambda=0$ triplet is split into a singlet and a doublet while the pentuplet is split into two doublets and a singlet. The energy levels are found to depend quadratically on $\lambda$: $E_{i}\propto\lambda^{2}$, indicating the absence of the linear DM antisymmetric exchange. In both cases this is expected on symmetry grounds. For the dumbbell this is straightforward, since there is an inversion center at the midpoint between the two triangular clusterspowell2016 ; DM2 . For the tube the C3 rotation symmetry implies that ${\bm{D}}^{DM}\parallel z$-axis (Moriya’s rule 5; Ref. DM2, ) and the $\sigma_{v}$ reflection symmetry implies that ${\bm{D}}^{DM}\parallel xy$-plane (Moriya’s rule 3). Both conditions taken together lead to ${\bm{D}}^{DM}={\bm{0}}$, and there is no DM coupling between the two spins in the tube arrangement.

In Mo3S7(dmit)3 the symmetry of the tube is lowered from $C_{3v}$ to $C_{3}$ by small intermolecular interactions neglected in the current model jacko2017 . This allows for a non-zero DM coupling parallel to the $C_{3}$ axis, which points along the crystallographic c-axis powell2016 .

The level degeneracies for both pairs of coupled clusters (Fig. 6) are those expected for an isotropic antiferromagnetic Heisenberg model with a trigonal single ion anisotropy described by Eq. (17), which we have seen arises for non-zero $\lambda$. This is expected for the tube, as in $C_{3v}$ symmetry there are two-fold degenerate states corresponding to the $E$ irreducible representation.

However, the $C_{i}$ symmetry of the dumbbell configuration admits only one-dimensional irreducible representations. Thus, one expects the level degeneracies associated with the trigonal symmetry to be fully lifted in the presence of SMOC. We will denote these level splittings as triclinic splittings. The absence of such triclinic splittings for $J_{F}=0$ in the dumbbell arrangement therefore indicates a hidden symmetry in the model. This is broken for $J_{F}\neq 0$. To quantify the degree of hidden symmetry breaking we plot the difference in energy between the second and third eigenstates, $E_{3}-E_{2}$ in Fig. 7. For $J_{F}=0$ no level splitting is present for any $\lambda_{xy}/\lambda_{z}$ ratio. However, a triclinic splitting arises as $-J_{F}$ is increased, saturating at sufficiently large $-J_{F}$. The largest splittings are found when SMOC is anisotropic, particularly when $\lambda_{xy}/\lambda_{z}>1$.

Thus, it is apparent that hidden symmetry is related to Coulomb matrix and is present in the absence of direct exchange interaction. For $J_{F}=0$ the symmetric and antisymmetric spin exchange tensors are proportional, but this is lifted for $J_{F}\neq 0$. This hidden symmetry plays a similar role in controlling the anisotropy of effective spin models of transition metal oxides.sea

We stress that the C3 rotation symmetry of the tube conformation forbids trigonal level splittings, even for $J_{F}\neq 0$. Consistent with this expectation, no triclinic level splittings are observed in our calculations for the tube configuration.

III.2 Second order perturbation theory in the intercluster hopping

In order to derive a low energy effective Hamiltonian for the two coupled clusters we now perform perturbation theory calculations to $O(t_{conf}^{2})$, where $conf=dumbbell$, $tube$ and $t_{dumbbell}=t$ and $t_{tube}=t_{z}$. The effective Hamiltonian for two coupled clusters with $N$ electrons in each cluster is given by

[TABLE]

where $E_{0}(N,j_{iz})$, is the energy of the isolated trimer, $i$, with $j_{iz}=0,\pm 1$ with $N$ electrons ($N=4$ in the case of interest here), with corresponding eigenstate $|N,j_{iz}\rangle$. In the expression above we are implicitly assuming that the ground state of isolated uncoupled trimers is three-fold degenerate even for non-zero SMOC. From a comparison to exact results and the canonical transformation, discussed below, we find that this approximation is very accurate for the parameter regime analyzed. The $\{|m_{0}\rangle\}$ are the complete set of virtual excitations in which an electron is transferred from one cluster to the other and may be written as

[TABLE]

where $A_{\gamma_{i}}(N\pm 1,\mu_{i})=\langle N\pm 1,\mu_{i}|N\pm 1,\gamma_{i}\rangle$, $\gamma_{i}$ denotes the excitations and $\mu_{i}$ runs over the the Hilbert state configurations with $N\pm 1$ electrons on trimer $i=\ell,m$.

Introducing these states in Eq. (LABEL:eq:hameff) we find for a given configuration of the coupled clusters

[TABLE]

where the excitation energies are $\Delta\epsilon(N-1,\gamma_{\ell};N+1,\gamma_{m})=2E_{0}(N)-(E_{\gamma_{\ell}}(N-1)+E_{\gamma_{m}}(N+1))$.

It is important to test the reliability of the present second order perturbative calculation for the values of the inter cluster hopping amplitudes relevant to Mo3S7(dmit)3 crystals. We have checked the accuracy of the second order perturbation theory calculations by comparing the nine lowest energy eigenstates with the exact eigenspectrum in our previous workmerino2016 . From Fig. 3 of [merino2016, ] we concluded that the second order, $O(t^{2})$, calculation is very accurate in the dumbbell arrangement with $U=10t_{c}$, even for the large inter-molecular hopping amplitude, $t=0.785t_{c}$ relevant to Mo3S7(dmit)3 crystals.

In the tube arrangement, comparable accuracies can only be achieved at larger $U$. The poorer accuracy at intermediate $U$ in the tube configuration is due the stronger charge fluctuations in this configurationnourse2016 ; janani2014a . In the tube particles can be exchanged between the two clusters through $\sim t_{z}^{2}/t_{c}$ processes without paying energy cost $\sim U$.powell2016 In contrast, in the dumbbell case, since particles can only be exchanged through the single hopping connecting the two vertices there is always an energy cost $\sim U$ inherent to the exchange process $\sim 4t^{2}/U$. In spite of this, at sufficiently large values of $U$ we find that the second order perturbation theory is sufficiently accurate for both the dumbbell and tube arrangements even for the large values of $t=0.785t_{c}$ and $t_{z}=0.683t_{c}$ extracted from DFT for Mo3S7(dmit)3.jacko2015 ; jacko2017

IV Effective magnetic spin exchange model

In order to determine the analytical form of the pseudospin exchange Hamiltonian, we have performed a canonical transformation. Analytical expressions of the pseudospin model valid to $O(\lambda^{2})$ and $O(t^{2})$, are obtained assuming a $t$-$J$ model for the triangular clusters, specified in Appendix B. By equating the matrix elements of the effective pseudospin exchange Hamiltonian obtained from the canonical transformation to the matrix elements of $H_{eff}^{(2)}$ evaluated in the low energy subspace, $\{|j_{\ell},j_{m}\rangle\}$, with $j_{\ell},j_{m}=0,\pm 1$, we are able to extract the parameters entering the pseudospin exchange model.

IV.1 Canonical transformation for a nearly degenerate low-energy subspace

Consider an arbitrary Hamiltonian, $H=H_{0}+H_{1}$ where $H_{0}=\sum_{\nu}P_{\nu}HP_{\nu}$, $H_{1}=\sum_{\mu\neq\nu}P_{\nu}HP_{\mu}$, and $P_{\nu}$ is a projector onto the $\nu$th subspace. Now define $H(\varepsilon)=H_{0}+\varepsilon H_{1}$. Let

[TABLE]

We choose $S$ so that the linear term vanishes, i.e., such that $iH_{1}=\left[H_{0},S\right].$ This implies that

[TABLE]

because $P_{\mu}P_{\nu}=P_{\mu}\delta_{\mu\nu}$ and $\sum_{\mu}P_{\mu}=1$. For $\mu=\nu$ this yields $P_{\mu}SP_{\mu}=\gamma P_{\mu}$ for $\gamma\in\mathbb{C}$. While, for $\mu\neq\nu$ we find

[TABLE]

If we choose the projectors such that they project onto strictly degenerate subspaces then

[TABLE]

Therefore, keeping only second order $O(\varepsilon^{2})$ terms, we find that

[TABLE]

Finally, we find the effective low-energy Hamiltonian by projecting onto the low-energy subspace, henceforth denoted $\cal L$. Here it is convenient to associate all of the subspaces with the states chosen so that the low energy subspace is diagonal, i.e., $P_{\mu}HP_{\nu}=0$ if $\mu\neq\nu$ and both $\mu$ and $\nu\in\cal L$. (This is always possible provided we can solve the problem restricted purely to $\cal L$, as in elementary degenerate perturbation theory.) We then find that

[TABLE]

where $P_{\cal L}=\sum_{\mu\in\cal L}P_{\mu}.$ In the case that $\cal L$ is strictly degenerate this reduces to the standard result. In the case where there is a small spread of energies in $\cal L$ and these are treated as a single subspace, as in the derivation of the $t$-$J$ model, a similar result holds but is approximate because the replacement of $P_{\mu}HP_{\mu}$ by its expectation value in Eq. (29) is no longer exact. We note that this is precisely the approximation made in Eq. (LABEL:eq:hameff) where we neglected the single-ion splitting of the ground state triplet in the denominator.

The effective Hamiltonian derived from this canonical transformation describing the coupling between two isolated nearest-neighbor trimers, $\ell$ and $m$, in the tube arrangement of Fig. 2b is

[TABLE]

where $J^{c}_{\alpha\beta}$ is diagonal and $J^{c}_{xx}=J^{c}_{yy}\neq J^{c}_{zz}$, and the anisotropic biquadratic couplings, $P_{\alpha\beta}=P_{\beta\alpha}$, obey $P_{xx}=P_{yy}=P_{xy}$ and $P_{zx}=P_{zy}=(P_{zz}+P_{xx})/2$. Both numerically and analytically we find $P_{xx}\ll P_{zz}$, indeed we find numerically that $P_{xx}$ is negligibly small and thus do not discuss it further below. $D^{c}=D+\Delta D^{c}$ is the single-spin anisotropy including corrections, $\Delta D^{c}$, due to hopping processes between the clusters. The perturbative expressions for these parameters are given in Appendix B. Thus, one can recast the bilinear terms of $H^{c}_{\ell m}$ in the familiar XXZ form. Doing so, one finds that the Hamiltonian for a single chain is

[TABLE]

where $J^{c}=J_{xx}^{c}$ and $\Delta^{c}=J_{zz}^{c}/J_{xx}^{c}$.

For two isolated nearest-neighbor trimers in the dumbbell arrangement with the $t$ bond connecting the two sites labeled ‘1’ (cf. Figs. 1 and 2a), the exchange Hamiltonian is

[TABLE]

$D^{ab}=D+\Delta D^{ab}$ is the single-spin anisotropy including corrections, $\Delta D^{ab}$, due to hopping processes between the clusters and is plotted in Figs. 8 and 9. We find that $\Delta D^{ab}$ is very small so that $D^{ab}\sim D$.

To derive the effective Hamiltonian for the full crystal we know need to note that we have, so far, only considered the $t$-bonds between Wannier orbitals labeled ‘1’, cf. Figs. 1 and 2, and Eq. 19. Rather than repeating the derivation for ‘2’ and ‘3’ bonds we can simply use the $C_{3}$ symmetry of the molecules and note that the $\bm{\mathcal{S}}_{\bm{r}_{m}}$ operators transform as vectors under rotation. Hence we can replace

[TABLE]

in Eq. (34), where $j$ labels the bond, as shown in Fig. 1.

Firstly, one finds that the $K_{\pm\pm}$ and $K_{z\pm}$ terms vanish in the full crystal due to cancellation among the contributions from the three nearest-neighbor bonds. Transforming the other terms, one can rewrite that Hamiltonian as

[TABLE]

where $J^{ab}=(J_{xx}^{ab}+J_{yy}^{ab})/2$, $\Delta^{ab}=J_{zz}^{ab}/J^{ab}$ and $Q=(J_{xx}^{ab}-J_{yy}^{ab})/2$. The perturbative expressions for these parameters are given in Appendix B. Thus, we see that the second term (proportional to $J^{ab}$) is simply the XXZ model and the third term (proportional to $Q$) is the honeycomb 120∘ compass model.nussinov2015

Finally, combining the results obtained above we obtain the full effective spin exchange model for the crystal, which reads:

[TABLE]

where $D^{*}=D+\Delta D^{c}+\Delta D^{ab}$. This expression neglects ‘three molecule’ terms analogous to the ‘three site’ terms neglected in the usual formulation of the $t$-$J$ model.Chao ; Harris We will see below that $J_{xz}^{ab}$ is extremely small. On neglecting this term one finds that the effective Hamiltonian is given by Eq. (I).

The parameters governing the spin exchange between molecules $\ell$ and $m$ in our spin exchange Hamiltonian, $H_{\text{eff}}$, are obtained by comparing the canonical transformation with our numerical second order perturbation theory

[TABLE]

recall $H^{(2),conf}_{eff}$ is defined in Eq. (25). The above equations are solved for a given set of parameters: $U$, $J_{F}$, $t_{c}$, $t$, $t_{z}$, $\lambda_{xy}$, and $\lambda_{z}$ entering our original microscopic model (2).

IV.2 Anisotropic exchange in the $ab$-plane

We have explored anisotropies arising in the exchange couplings of the effective exchange model, Eq. (34) for the two clusters coupled as in Fig. 2(a). Since the non-pseudospin-conserving $K_{\alpha\beta}$ terms exactly cancel in the crystal they will not be discussed any further. We find that when $J_{F}=0$ the exchange coupling tensor is diagonal and isotropic, $J^{ab}_{\alpha\beta}=J^{ab}\delta_{\alpha\beta}$. This is consistent with our previous results (see Fig. 4(a) of Ref. [merino2016, ]) and the lack of triclinic splittings observed in the energy level spectrum for two clusters in the dumbbell configuration shown in Fig. 7.

As shown in Fig. 8 anisotropic exchange couplings $J^{ab}_{xx}\neq J^{ab}_{yy}\neq J^{ab}_{zz}$ arise when $J_{F}\neq 0$, which are consistent with the triclinic splittings found in the exact level spectrum of Fig. 7. Also we find off-diagonal exchange couplings: $J_{xz}^{ab}\neq 0$ to all orders of SMOC consistent with the analytical expression for $J_{xz}^{ab}$ derived in our previous work merino2016 valid to O($\lambda_{xy}^{2},\lambda_{z}^{2}$). However, we typically find small values of $J_{xz}^{ab}\sim-0.00044t_{c}~{}(\lambda_{xy}=\lambda_{z}/2)$ and $J_{xz}^{ab}\sim-0.0003t_{c}~{}(\lambda_{xy}=\lambda_{z})$ and so this parameter is not displayed in Fig. 8. Therefore, to an excellent approximation, the in-plane Hamiltonian is an XXZ + 120∘ honeycomb model with single ion anisotropy.

Comparing the results shown in Fig. 8 for different $\lambda_{xy}/\lambda_{z}$ ratios, we observe that the anisotropies in the exchange couplings are enhanced for $\lambda_{xy}/\lambda_{z}\neq 1$. In fact, larger anisotropies are found to occur for $\lambda_{xy}=\lambda_{z}/2$, which is the parameter regime relevant to Mo3S7(dmit)3 crystals.jacko2017 Also note from Fig. 8 the strong dependence of the magnitude of $D$ on the SMOC anisotropy. The single-spin anisotropy increases rapidly with SMOC, becoming equal to the exchange couplings, $D\sim J^{ab}$ at $\lambda\approx 0.45t_{c}$ ($\lambda=\lambda_{xy}=\lambda_{z}$), at $\lambda_{xy}\approx 0.22t_{c}$ ($\lambda_{xy}=\lambda_{z}/2$) and at $\lambda_{xy}\approx 1.045t_{c}$ ($\lambda_{xy}=2\lambda_{z}$). At sufficiently large $D\gtrsim J^{ab}$ we expect the $D$-phase, i.e., a tensor product of $j=0$ states located at each cluster of the crystal. Hence, a $D$-phase is favored by anisotropic SMOC with $\lambda_{xy}<\lambda_{z}$.

In Fig. 8 we also show results for an antiferromagnetic exchange coupling inside the cluster, $J_{F}>0$. This could arise in, say, Mo3S7(dmit)3 due to superexchange via the sulphur atoms in the core. We find similar spin exchange anisotropies for both ferromagnetic and antiferromagnetic $J_{F}$. In the antiferromagnetic case we find that $D^{ab}$ becomes negative for sufficiently large SMOC and $\lambda_{xy}=2\lambda_{z}$, consistent with our perturbative results for the $t$-$J$ model [cf. Eqs. (51), (53a), and (55a)]. This signifies a switch of the ground state of the isolated cluster from the $j=0$ singlet to the $j=\pm 1$ doublet. In contrast, in the ferromagnetic cases, $J_{F}<0$, we have explored a large parameter set and we always find $D^{ab}>0$.

In order to understand the effect of exchange couplings with SMOC anisotropy, we show in Fig. 9 exchange couplings, $J^{ab}_{\alpha\alpha}$, and $D$ in two extreme cases: $\lambda_{xy}=0$ and $\lambda_{z}=0$ with $J_{F}=-0.3t_{c}$. The $J^{ab}_{\alpha\alpha}$ are suppressed (enhanced) with SMOC for $\lambda_{xy}=0$ ($\lambda_{z}=0$), consistent with the results shown in Fig. 8. Only when $\lambda_{xy}$ is turned on, does one find that the transverse couplings become different i.e., $J_{xx}\neq J_{yy}$. Furthermore, the single-spin anisotropy is much more strongly enhanced by $\lambda_{xy}$ than by $\lambda_{z}$ (by more than an order of magnitude), consistent with the analytical expressions [see Eqs. (51), (53a), and (55a)].

IV.3 Anisotropies in the exchange interactions along the $c$-direction

The exchange couplings between two neighboring clusters in the $c$-direction are shown in Fig. 10. We find a diagonal exchange tensor: $J^{c}_{\alpha\beta}=J^{c}_{\alpha\alpha}\delta_{\alpha\beta}$, with $J^{c}_{xx}=J^{c}_{yy}\neq J^{c}_{zz}$ for any $J_{F}\neq 0$ and $\lambda_{xy}/\lambda_{z}$ ratio. The higher symmetry than for a pair of molecules in the $ab$ plane is due to the $C_{3}$ rotational symmetry of the tube dimer, as discussed above.

The largest anisotropies with $J^{c}_{xx}=J^{c}_{yy}>J^{c}_{zz}$ are seen in the case of anisotropic SMOC with $\lambda_{xy}=\lambda_{z}/2$ as shown in Fig. 10(b). The only non-negligible biquadratic exchange terms, $P_{zz}>P_{zx}$, increase rapidly with $\lambda_{xy}$ starting to saturate around $\lambda/t_{c}\sim 1-1.5$. The single-spin anisotropy equals the exchange coupling, $D=J^{c}$, at $\lambda_{xy}=0.65t_{c}$ for $\lambda_{xy}/\lambda_{z}=1/2$ and at $\lambda_{xy}=1.457t_{c}$ for $\lambda_{xy}/\lambda_{z}=1$, while for $\lambda_{xy}/\lambda_{z}=2$ there is no critical $\lambda_{xy}$ at which $D\sim J^{c}$ within the parameter range explored. Hence, anisotropic SMOC with $\lambda_{xy}<\lambda_{z}$ again favors the $D$-phase as in the dumbbell arrangement.

Finally, in Fig. 11 we compare the dependence of the exchange couplings on $J_{F}$ for $\lambda_{xy}=\lambda_{z}=1$. The couplings in the $ab$ plane, $J^{ab}_{\alpha\alpha}$ are suppressed and become gradually anisotropic, $J^{ab}_{xx}\neq J^{ab}_{yy}\neq J^{ab}_{zz}$ as $J_{F}$ increases. This is in contrast to the exchange couplings in the $c$-direction which do not display larger anisotropies but rather $J^{c}_{xx}=J^{c}_{yy}\neq J^{c}_{zz}$ for any $J_{F}$.

V Discussion of properties of the quasi-one-dimensional pseudospin-one model

Our analysis shows that the magnetic properties of layered decorated honeycomb lattice model at strong coupling, $U\gg t_{c},t,t_{z},\lambda_{xy},\lambda_{z}$, are captured by model (37) with the exchange couplings obtained from our combined approach described above. On comparing $J^{ab}$ in Fig. 8 with $J^{c}$ in Fig. 10 we find that $J^{c}\sim 5J^{ab}$ for $U=10t_{c}$. This is related to the fact that two clusters in the tube arrangement are connected by three hoppings so that they can exchange electrons without paying an energy costmerino2016 ; powell2016 $\sim U$. This mechanism is generic to decorated lattices and not specific to the model considered here.elise2017 In contrast, neighboring clusters in the dumbbell arrangement pay energy, $U$, since they can only exchange particles through a single hopping connecting them. Hence, $J^{ab}$ is strongly suppressed by $U$ in contrast to $J^{c}$, leading to an increase of the $J^{c}/J^{ab}$ ratio. Hence, at large $U$ the system becomes quasi-one-dimensional consisting on a set of weakly coupled pseudo spin-one antiferromagnetic chains.

An isotropic version of model, (IV.1) i.e., $J^{c}_{\alpha\beta}=J^{c}\delta_{\alpha\beta}$ $P_{\alpha\beta}=P\delta_{\alpha\beta}$ and $D^{c}=0$ is just the bilinear-biquadratic model: $H=J^{c}\sum_{\ell}\bm{\mathcal{S}}_{\bm{r}_{\ell}}\cdot\bm{\mathcal{S}}_{\bm{r}_{\ell}+\bm{\delta}_{z}}+P\sum_{\ell}(\bm{\mathcal{S}}_{\bm{r}_{\ell}}\cdot\bm{\mathcal{S}}_{\bm{r}_{\ell}+\bm{\delta}_{z}})^{2}$, which becomes the Affleck- Kennedy-Lieb-Tasaki (AKLT) model for $P/J^{c}=1/3$. The AKLT model can be solved exactly and has the valence bond solid ground state and is in the Haldane phaseAKLT .

We finally note that the next-nearest-neighbor exchange couplings between clusters in the $c$-direction can be neglected since recent estimates merino2016 suggest that they are about 20 times smaller than the nearest neighbor exchange coupling. This is because the small parameter in the perturbation theory is $t_{z}/3$ so fourth order terms (such as next-nearest-neighbor exchange couplings) must be at least an order of magnitude smaller than second order terms (such as nearest-neighbor exchange coupling).

V.1 One-dimensional antiferromagnetic $S=1$ Heisenberg chains

When no interchain coupling is present, $J^{ab}=0$, and $D^{*}<J^{c}$, the system consists on a set of uncoupled one-dimensional $S=1$ antiferromagnetic chains that are in the Haldane phase. The Haldane phase is characterized by exponentially decaying spin correlations associated with dmrg the Haldane spin gap $\Delta_{s}=0.4107(3)J^{c}$ to the lowest triplet state and string order. It is a symmetry-protected topological phase with nonlocal string order and fractionalized edge states haldane1983a ; haldane1983b ; gomezsantos1989 . Topological protection can arise from either (i) the dihedral group of $\pi$-rotations around the $x$ and $y$ axis, (ii) time-reversal symmetry or (iii) reflection through a plane perpendicular to the chain (or bond-center inversion symmetry, which is equivalent in one-dimension).pollmann2012 In the underlying fermionic model, charge fluctuations imply that topological protection can only come from reflection symmetry with respect to a plane perpendicular to the $c$-axis at the midpoint of a bondnourse2016 .

On the other hand, when $D^{*}\gg J^{c}$, the ground state is adiabatically connected to a trivial state consisting on the tensor product of the ${\mathcal{S}}^{z}_{\bm{r}_{\ell}}=0$ at each cluster. The lowest energy excitations of the $D$-phase which reside in the ${\mathcal{S}}^{z}_{\bm{r}_{\ell}}=\pm 1$ sector, are gapped and consist of pairs of excitons and antiexcitons which can be bound. Numerical studies oitmaa2009 ; normand2011 ; langari2013 ; tzeng2008 have established that in the pure spin model the quantum critical point separating the $D$-phase and Haldane phase occurs at $D^{*}/J^{c}\sim 0.96-0.971$. It has been found that in a pure spin model such as the one discussed here, a quantum phase transition between the Haldane phase and the topologically trivial $D$-phase is signalled by the change in sign of an inversion-symmetry-based order parameter langari2013 which is a non-local topological order parameter. Hence, a transition from a Haldane phase to a $D$-phase occurs when increasing SMOC until $D^{*}\sim J^{c}$.

From our analysis of Fig. 10 (b), which is the relevant SMOC ratio to Mo3S7(dmit)3, (assuming $J_{F}=-0.3t_{c}$), we predict a transition from the Haldane to the $D$-phase at $\lambda_{xy}\sim 0.65t_{c}$. Ab initio estimates of SMOCjacko2017 in Mo3S7(dmit)3 find that $\lambda_{xy}=\lambda_{z}/2=0.042t_{c}$, which would naively mean that the single-spin anisotropy is too small, $D^{*}\ll J^{c}$, to induce a $D$-phase in the crystal. By moving to suitable materials containing heavier elements khosla2017 , SMOC can be increased by, at most, a factor of $4-5$ leading to $\lambda_{xy}\approx 0.2t_{c}\ll\lambda^{critical}_{xy}$ which means that the system is still in the Haldane phase. However, the critical $\lambda_{xy}$ for the transition can be reduced by suppressing $t_{z}$ and increasing $-J_{F}$ as shown in Fig. 3. Also in the underlying fermionic model (neglecting SMOC), the Haldane gap is suppressed by more than an order of magnitude by charge fluctuations.nourse2016 More specifically, charge fluctuations renormalize the critical condition $D^{*}\sim J^{c}$ to $D^{*}\sim 0.066J^{c}$ for the parameters relevant to Mo3S7(dmit)3. This leads to a smaller $\lambda^{critical}_{xy}$ as shown in Fig. 3. The above discussion indicates that a series of materials related to Mo3S7(dmit)3 with slight variations in model parameters could easily effectively span the phase Haldane–to–$D$-phase transition. Furthermore, a material on the D-phase side of the transition could be driven into the Haldane phase by uniaxial pressure along the c-axis. In particular, our results above suggest that the critical ratio $D^{*}/J^{c}$ could be exceeded by moving to suitable materials containing heavier elements.khosla2017 Furthermore, one expects that the interlayer hopping $t_{z}$ will be extremely sensitive to chemical details. As $J^{c}\sim t_{z}^{2}$ structures with increased interlayer separation will strongly favor the $D$-phase.

V.2 Effect of the interchain couplings

When the quantum pseudospin-one chains are coupled through a sufficiently strong interchain coupling, $J^{ab}$, the Haldane phase becomes unstable to 3D magnetic order. In previous numerical studies of weakly coupled $S=1$ antiferromagnetic Heisenberg chains (with $D^{*}=0$), it was estimatedsengupta2014 that the critical value for the transition from the Haldane to the ordered 3D magnet occurs around $J^{ab}/J^{c}\geq(0.08-0.11)z\sim 0.3$, where the coordination number $z=3$ for the honeycomb lattice. Since we find that $J^{ab}/J^{c}\lesssim 0.2$, we expect that the ground state of our model is in the Haldane phase when $D^{*}=0$. This critical ratio, $J^{ab}/J^{c}$, for the onset of 3D magnetic order is suppressed by $D^{*}$ as shownsakai1990 by mean-field treatments of the interchain coupling, $J^{ab}$.

V.3 Effect of an external magnetic field

An external magnetic field suppresses the 1D quantum fluctuations and the Haldane gap, $\Delta_{s}$, closesaffleck1990 at $h_{c}\sim\Delta_{s}$, whence a transition to a 3D ordered magnet occurs. A quantum critical region with a V-shape emerges around $h_{c}$ in the temperature versus magnetic field, $T$-$h$, phase diagramsakai1990 ; giamarchi2007 ; bera2015 . The temperature, $T\sim J^{ab}$, sets the energy scale at which 3D quantum criticality for $T<J^{ab}$ crosses over to 1D behavior for $T>J^{ab}$. Similarly the three-dimensional magnetically ordered phase found for $h>h_{c}$ and $T=0$ crosses over to a gapless Tomonaga Luttinger Liquid (TLL) at temperatures $T>J^{ab}$. We note that, strictly speaking, the TLL behavior should only occurcapponi2016 in the range $J^{ab}<T<J^{c}$, since at too large temperatures, $T\gg J^{c}$, classical behavior sets in. In the presence of a nonzero and small $D^{*}$, with $D^{*}\ll J^{c}$, the lowest triplet state is split into a $j=\pm 1$ doublet with energy $\Delta_{\pm}$ above the ground state and a $j=0$ singlet at energy $\Delta_{0}$ with $\Delta_{\pm}<\Delta_{0}$. Hence, under an applied magnetic field, $\Delta_{\pm}$, is suppressed and the transition from the Haldane phase to the 3D ordered phase occurs around $h_{c}=\Delta_{\pm}<\Delta_{s}$. Apart from the downward shift of $h_{c}$, we can expect, qualitatively, a similar $T$-$h$ phase diagram as in the case with no single-spin anisotropy, $D^{*}=0$.

VI Conclusions

We have analyzed the magnetic properties of the trinuclear organometallic materials, such as Mo3S7(dmit)3. These materials are potential candidates for realizing compass interactions in their layers. In order to explore such possibilities we have derived an effective magnetic model describing the magnetic interactions between the pseudospin-one at each molecular cluster arising from strong Coulomb repulsion, lattice structure and SMOC. In spite of the crystals being nearly isotropic, we find that the exchange coupling between nearest-neighbor pseudospins along the $c$-direction is much larger than between pseudospins within the hexagonal $a$-$b$ planes. Hence, the spin exchange model for these crystals is effectively quasi-one-dimensional. Magnetic anisotropies are found to arise under the simultaneous effect of spin orbit coupling and intra-cluster exchange interaction. These anisotropies are further enhanced by SMOC anisotropy, particularly when $\lambda_{xy}<\lambda_{z}$, which is naturally present in organometallics. Our analysis suggests that Mo3S7(dmit)3 is most probably in the Haldane phase since the efffective model consists of weakly coupled $S=1$ antiferromagnetic chains in the presence of small single-spin anisotropy induced by SMOC. However, by increasing the interlayer distances through changes in the chemistry of the material, increasing the anisotropy of magnitude of the SMOC it should be possible to effectively drive it into to the $D$-phase. A larger SMOC should be realised in complexes containing heavier metals.khosla2017

The Haldane phase is strongly sensitive to an external magnetic field. Under applied magnetic fields larger than the Haldane gap, $h>h_{c}\sim\Delta_{s}$, the Haldane phase is destroyed and a three-dimensional magnet may be stabilized. We have estimated this critical field, $h_{c}$, based on our present analysis using DFT parametersjacko2017 for Mo3S7(dmit)3 (Table 1) with an onsite $U=10t_{c}$ and $J_{F}=-0.3t_{c}$. Using these parameters we extract $J^{c}=0.0126$ eV from our Fig. 10(b) which leads to a critical magnetic field $h_{c}\sim\Delta_{s}\sim 41.4$ T assuming the Haldane spin gap, $\Delta_{s}=0.414J^{c}\approx 0.09t_{c}$, in the pure Haldane chain. However, recent DMRG calculations on Hubbard tubesnourse2016 have shown that charge fluctuations strongly suppress the spin gap when decreasing $U$. For the parameter range considered here, we would find: $\Delta_{s}\sim 0.006t_{c}$, implying experimentally accessible critical fields: $h_{c}\sim 3$ T. A V-shaped quantum critical region in the $T-h$ phase diagram separating the Haldane phase from the three-dimensional magnetically ordered phase should then emerge as observed in inorganic Haldane chain materials.bera2015

Exfoliation or growth of a monolayer of trinuclear complexes arranged as in the ab-planes of Mo3S7(dmit)3, would lead to the realization of a decorated hexagonal lattice which is known to contain rich physics. We have found that at large $U$ and no SMOC, the magnetic interactions between the pseudospin-one is that of a conventional nearest-neighbor antiferromagnetic Heisenberg model on an hexagonal lattice.footnote The ground state of this model is a pure Néel antiferromagnet. However, if crystal parameters are tuned so that magnetic exchange anisotropies are enhanced, disordered spin liquid phases sheng2015 may be achieved. For instance, if the relative orientation between the molecules in the crystal is modified so that inversion symmetry within the planes is broken, a DM interaction arisespowell2016 which competes with the magnetic order,cepas2008 which can lead to interesting spin liquid phasesmessio2017 . All this illustrates how isolated layers of trinuclear organometallic complexes are ideal playgrounds to explore the quantum many-body phases realized in a decorated honeycomb lattice.

Acknowledgements.

J. M. acknowledges financial support from (Grant No. MAT2015-66128-R) MINECO/FEDER, Unión Europea.. Work at the University of Queensland was supported by the Australian Research Council (FT13010016 and DP160100060) and by computational resources provided by the Australian Government through Raijin under the National Computational Merit Allocation Scheme..

Appendix A Electronic structure of isolated triangular clusters

Here, we provide the details of the electronic structure of isolated clusters with different numbers of electrons.

A.1 Isolated triangular cluster with five electrons

We start studying isolated trimers with $N=1$ electrons. This is due to its intrinsic importance and due to the fact that the electronic structure of trimers with $N=5$ electrons and $t_{c},\lambda_{xy},\lambda_{z}>0$, relevant to Mo3S7(dmit)3 can be obtained from the $N=1$ case by a particle-hole transformation switching the sign of the parameters: $t_{c}\rightarrow-t_{c},\lambda_{xy}\rightarrow-\lambda_{xy},\lambda_{z}\rightarrow-\lambda_{z}$ apart from a rigid energy shift.

For only one electron in the cluster, $N=1$, the Hamiltonian is just $H=H_{0}+H_{SMOC}$. Since $[H,J_{z}]=0$, where $J_{z}=L_{z}+S_{z}$, then the projection of the total momentum along the $z$-axis is a good quantum number. In the following we denote the basis states for a fixed number of particles, $N$, as $|N;j,n\rangle$ where $j=k+\sigma$ and $n$ numbers the different possible configurations for each $j$-sector. Hence, in this case the possible basis states are

[TABLE]

The eigenenergies, $E_{n}(N;j)$ of the Hamiltonian are

[TABLE]

Hence the level spectra for $N=1$ consists of three doublets with the energies given above. The ground state of the system with one electron, $N=1$, is a doublet with energy, $E_{0}$. Time-reversal invariance of the Hamiltonian, $[T,H]=0$, and Kramers theorem ensures that all states should have a minimum degeneracy of two since the cluster has an odd number of electrons. Note that the level spectra of the triangular cluster with $N=5$ electrons (one hole) would be the same as (LABEL:eq:spect) but with the signs reversed: $t_{c}\rightarrow-t_{c}$, $\lambda_{xy}\rightarrow-\lambda_{xy},\lambda_{z}\rightarrow-\lambda_{z}$ and with an upward rigid shift of all energies by $+2U$.

To make contact with previous work on transition metal oxides it is illustrative to consider our model Hamiltonian: $H=H_{0}+H_{SMOC}+H_{U-J_{F}}$, with $H_{0}$, $H_{SMOC}$ and $H_{U-J_{F}}$ expressed in the $(k,\sigma)$ basis as given by Eq. (12), (14) and (16), respectively. For $U,J_{F}=0$, this model is reminiscent of a model previously considered jackeli2009 ; perkins2014a ; perkins2014b for Ir4+ ions in $A_{2}$IrO3 ($A$=Na,Li) compounds. In these systems, five electrons occupy the lowest $t_{2g}$ manifold of the Ir ions which is well separated from the high energy $e_{g}$ doublet. The low energy effective model for the hole in the $t_{2g}$ manifold of the *isolated * Ir-ions includes a trigonal crystal field resulting from the surrounding oxygen octahedra and a large SOC contributionperkins2014b : $H=\Delta(L^{z})^{2}+\lambda{\bf L}\cdot{\bf S}$, with $\Delta>0$.

Through the particle-hole transformation discussed above, the three-fold degenerate $t_{2g}$ manifold of the isolated Ir ion with one hole is equivalent to our model of the isolated molecule with one electron, $N=1$, with the signs of $\lambda=\lambda_{xy}=\lambda_{z}$ and $t_{c}$ reversed. Full rotational symmetry is only recovered for $t_{c}\rightarrow 0$ in our model when $U,J_{F}=0$. In that case, $[H_{0}+H_{SMOC},L]=0$, so that the total angular momentum, $L$, is a good quantum number, as it should. In this situation, we find that isotropic SOC ($\lambda_{xy}=\lambda_{z}=\lambda$) splits the $(2L+1)(2S+1)=6$ manifold ($L=1,S=1/2$) into a $j=1/2$ doublet with energy $E_{0}(1;j=1/2)=-\lambda$ and a $j=3/2$ quadruplet with energy $E_{1}(1;j=3/2)={\lambda\over 2}$. This situation corresponds to removing the crystal field acting on the $d$-orbital manifold in transition metal oxides.

A.2 Isolated triangular clusters with four electrons

The basis states with $N=N_{\uparrow}+N_{\downarrow}=4$ electrons includes states with total spin $S_{z}=0$, ($N_{\uparrow}=2,N_{\downarrow}=2$), $S_{z}=1$ ($N_{\uparrow}=3,N_{\downarrow}=1$) and $S_{z}=-1$ ($N_{\uparrow}=1,N_{\downarrow}=3$). Noting that basis states with total momentum $k^{\prime}$ are equivalent to $k$ if they satisfy $k=k^{\prime}\pm 3n$, we find that the basis states can be classified according to three possible values: $j=0,\pm 1$. Since the Hamiltonian does not mix sates with different $j$, the original $15\times 15$ matrix can be expressed in block diagonal form consisting of $5\times 5$ matrices corresponding to $j=0,\pm 1$. We now explicitly show the classification of the $(k,\sigma)$ basis states according to $j=0,\pm 1$ and the analytical diagonalization of the matrices corresponding to each of the $j$-sectors. We keep the $|N;j,n\rangle$ classification of the basis states.

A.2.1 $j=0$ sector

The three possible configurations with $k=\sigma=0$ are

[TABLE]

[TABLE]

There is only one configuration for either $k=-1,\sigma=1$,

[TABLE]

or $k=1,\sigma=-1$

[TABLE]

Hence, the $j=0$ Hamiltonian reduces to a $5\times 5$ matrix:

$H(4;j=0)=\begin{bmatrix}-2t_{c}+{4U\over 3}-{7J_{F}\over 6}-\lambda_{z}&{U-J_{F}/2\over 3}&-{U+J_{F}\over 3}&0&0\\ {U-J_{F}/2\over 3}&-2t_{c}+{4U\over 3}-{7J_{F}\over 6}+\lambda_{z}&-{U+J_{F}\over 3}&-{\lambda_{xy}\over\sqrt{2}}&{\lambda_{xy}\over\sqrt{2}}\\ -{U+J_{F}\over 3}&-{U+J_{F}\over 3}&4t_{c}+{4U\over 3}-{5J_{F}\over 3}&{\lambda_{xy}\over\sqrt{2}}&-{\lambda_{xy}\over\sqrt{2}}\\ 0&-{\lambda_{xy}\over\sqrt{2}}&{\lambda_{xy}\over\sqrt{2}}&U+t_{c}-J_{F}+{\lambda_{z}\over 2}&0\\ 0&{\lambda_{xy}\over\sqrt{2}}&-{\lambda_{xy}\over\sqrt{2}}&0&U+t_{c}-J_{F}+{\lambda_{z}\over 2}\\ \end{bmatrix}$

A.2.2 $j=-1$ sector

We work in the basis

[TABLE]

The first three states have $k=-1=2$, $\sigma=0$, the fourth has $k=0$, $\sigma=-1$ and the fifth has $k=1=-2$, $\sigma=1$. The $j=-1$ Hamiltonian is

[TABLE]

A.2.3 $j=+1$ sector

It is convenient to take the basis states as the time-reversed analogues of the $j=-1$ sector:

[TABLE]

[TABLE]

Thus one immediately sees that $H(4;j=+1)=H(4;j=-1)$. Hence, there is a double degeneracy of the eigenvalues $E_{i}(4;j=+1)=E_{i}(4;j=-1)$.

For $\lambda=0$, the ground state is three-fold degenerate corresponding to the $S=1$ triplet combination of the two unpaired spins in the cluster. These lowest three degenerate states correspond to $j=0,\pm 1$. From the above analysis we conclude that isolated clusters with four electrons can be described through the effective Hamiltonian given in Eq. (17) where $D$ is an increasing function of SMOC as discussed in the main text.

A.3 Isolated triangular clusters with three electrons

The basis for $N=3$ electrons consists of 20 configurations: 18 configurations with $S_{z}=1/2$ ($N_{\uparrow}=2,N_{\downarrow}=1$) or $S_{z}=-1/2$ ($N_{\uparrow}=1,N_{\downarrow}=2$) and 2 configurations with $S_{z}=3/2$ ($N_{\uparrow}=3,N_{\downarrow}=0$) or $S_{z}=-3/2$ ($N_{\uparrow}=0,N_{\downarrow}=3$). The only allowed $j$ values for the cluster with $N=3$ electrons are $j=\pm{1\over 2},+{3\over 2}$ with the largest ($8\times 8$) matrix corresponding to $j=+{3\over 2}$. The $j=-{3\over 2}$ sector is not given here since the configurations are just the same as the ones in the $j=+{3\over 2}$ sector.

A.3.1 $j=+3/2$

The configurations with $j=3/2$ are

[TABLE]

Yielding the $8\times 8$ Hamiltonian matrix

[TABLE]

A.3.2 $j=\pm 1/2$

We take the basis

[TABLE]

and analogously for $j=-1/2$. The $6\times 6$ Hamiltonian matrix reads

$H(3;j=+1/2)=\begin{bmatrix}-3t_{c}+{2U\over 3}-{5J_{F}\over 6}-{\lambda_{z}\over 2}&{U+J_{F}\over 3}&{U+J_{F}\over 3}&0&0&{\lambda_{xy}\over\sqrt{2}}\\ {U+J_{F}\over 3}&3t_{c}+{2U\over 3}-{5J_{F}\over 6}+{\lambda_{z}\over 2}&-{U+J_{F}\over 3}&0&-{\lambda_{xy}\over\sqrt{2}}&0\\ {U+J_{F}\over 3}&-{U+J_{F}\over 3}&{2U\over 3}-{5J_{F}\over 6}&0&0&0\\ 0&0&0&{2U\over 3}-{J_{F}\over 3}&{U-J_{F}/2\over 3}&-{U-J_{F}/2\over 3}\\ 0&-{\lambda_{xy}\over\sqrt{2}}&0&{U-J_{F}/2\over 3}&{2U\over 3}-J_{F}/3+\lambda_{z}&{U-J_{F}/2\over 3}\\ {\lambda_{xy}\over\sqrt{2}}&0&0&-{U-J_{F}/2\over 3}&{U-J_{F}/2\over 3}&{2U\over 3}-{J_{F}\over 3}-\lambda_{z}\\ \end{bmatrix}$ Due to Kramers theorem the eigenstates, $E_{n}(3;j=1/2)=E_{n}(3;j=-1/2)$ and the energy levels for $E_{n}(3;j=3/2)$ are at least doubly degenerate. With no SMOC present, $E_{n}(3;j=\pm 1/2)=E_{n}(3;j=3/2)$ and the eigenstates are four-fold degenerate. However, when SMOC is present $E_{n}(3;j=\pm 1/2)\neq E_{n}(3;j=3/2)$ and the four-fold degeneracy is broken leading to two-fold degenerate levels.

Appendix B Expression for effective spin models from the canonical transformation of the $t$-$J$ model

In this appendix we model the $\ell$th trinuclear complex by the three site $t$-$J$ model, i.e.,

[TABLE]

where $h^{\dagger}_{\ell i\sigma}=a_{\ell i\sigma}$ creates an hole with spin $\sigma$ in the $i$th Wannier orbital and $P_{0}$ projects out states that contain empty sites. Note that it is important to retain the ‘three site’ terms here, as we will need to consider states far from half-filling. For a single molecule the effective low-energy model, retaining only the three lowest energy states is given by Eq. (17) with

[TABLE]

The $t$-$J$ model of the interlayer coupling between neighbouring molecules $\ell$ and $m$ is

[TABLE]

where now three are no three site terms because of the topology of underlying tight-binding model [cf. Eq. (21) and Fig. 2b]. Performing the canonical transformation described in section IV.1 and retaining quadratic terms in $t_{z}$, linear terms in $J_{z}$ (as $J_{z}$ is already quadratic in $t_{z}$) and quadratic terms in the SMOC (i.e., up to order $\lambda_{z}^{2}$, $\lambda_{xy}^{2}$, or $\lambda_{xy}\lambda_{z}$) yields an effective Hamiltonian described by Eq. (IV.1) with

[TABLE]

The $t$-$J$ model of the in-plane coupling between molecules $\ell$ and $m$ along a ‘1-bond’ (cf. Fig. 1) is

[TABLE]

again the three site terms vanish because of the underlying tight-binding model [Eq. (19)]. Performing the canonical transformation, adding in the 2- and 3-bonds, as described in section IV.1, and retaining quadratic terms in $t_{g}$, linear terms in $J_{g}$ and quadratic terms in the SMOC yields an effective Hamiltonian described by Eq. (36) with

[TABLE]