We explore the magnetic properties of isolated *a* − *b* planes of trinuclear organometallic crystals, Mo_{3}S_{7}(dmit)_{3}, in which an interplay of strong electronic correlations and spin molecular-orbital coupling (SMOC) occurs. The magnetic properties can be described by a XXZ+120^{0}, *S* = 1 Heisenberg model on a honeycomb lattice with single-spin anisotropy, *D*, which depends strongly on SMOC. Based on *ab initio* estimates of SMOC in Mo_{3}S_{7}(dmit)_{3} crystals, we predict that the honeycomb layers of Mo_{3}S_{7}(dmit)_{3} are Néel ordered. However, in materials with a greater degree of magnetic frustration, Néel order can give way to a large-*D* phase.

## I. INTRODUCTION

Strong spin orbit coupling (SOC) in weakly interacting systems can lead to unconventional insulating states such as topological insulators. Strongly correlated phases such as the topological Mott insulator can emerge from the interplay of strong Coulomb repulsion and SOC which may be realized in Ir-based transition metal oxides.^{1} For instance, in Na_{2}IrO_{3} and Li_{2}IrO_{3} materials, SOC removes the orbital degeneracy of 5*d* electrons leading to *S* = 1/2 pseudospins which interact through anisotropic and quantum compass exchange interactions on a honeycomb lattice. These materials are potential realizations of the Heisenberg-Kitaev model, closely related to the Kitaev model which sustains a spin liquid state.^{2}

Multinuclear coordinated organometallic complexes are strongly correlated systems in which spin-orbit coupling (SOC) can be relevant. The multinuclear complex, Mo_{3}S_{7}(dmit)_{3}, consists of honeycomb networks of Mo_{3}S_{7}(dmit)_{3} molecules stacked on top of each other along the *c*-direction of the crystal. Since the Mo_{3}S_{7}(dmit)_{3} molecules can be described by three Wannier orbitals,^{3} their packing on the honeycomb lattices within the layers lead to decorated honeycomb lattices, as shown in Fig. 1. Based on perturbative expansions including Coulomb repulsion and SMOC, the effective spin exchange model which describes the magnetic properties of these layers is a *S* = 1 XXZ+120^{0} quantum compass model on the honeycomb lattice.^{4,5}

In the present paper, we discuss the role played by the single-spin anisotropy induced by SMOC in the magnetic properties of the honeycomb layers of Mo_{3}S_{7}(dmit)_{3} crystals shown in Fig. 1. We predict that isolated honeycomb layers of Mo_{3}S_{7}(dmit)_{3} crystals are Néel ordered but increasing the magnetic frustration of the lattice can drive the system into a large-*D* phase.

## II. SPIN ORBITAL MOLECULAR COUPLING IN TRINUCLEAR ORGANOMETALLIC COMPLEXES

### A. Model for isolated trimers

The isolated triangular clusters can be described by a Hubbard-Heisenberg model in the presence of the spin molecular-orbit coupling (SMOC):

where:

with *t*_{c} the hopping between hybrid metal-ligand orbitals at sites in the triangular clusters and $ai\sigma \u2020$ creates an electron at the *i*th Wannier orbital with spin *σ*.

An important question is how SOC affects the low energy states described by these Wannier orbitals. Generically the one-electron spin-orbit coupling takes the form *H*_{SOC} = *∑*_{α}*K*_{α} ⋅*s*_{α}, where *s*_{α} is the spin of the *α*^{th} electron and *K*_{α} is a pseudovectorial operator that acts only on the spatial part of the wavefunction. For example, in the Pauli approximation,

where ** p** is the single electron momentum operator and

*ϕ*(

**) is the potential through which the electron moves.**

*r*Spherical symmetry implies that for an isolated atom ** K** ∝

**, the angular momentum operator. Thus, in**

*L**ab initio*calculations, it is common to write the SOC as a linear superposition of atomic contributions, for example, in the mean-field Breit-Pauli approximation

^{7,8}

where $ZAeff$ is the effective charge of nucleus *A* which accounts for the screening effects of the rest of the electrons on the nuclear potential, *l*_{αA} = (*r*_{α} − *R*_{A}) × *p*_{α}, is the angular momentum of the *α*^{th} electron relative to the nucleus A at *R*_{A}. This atomic approach does not provide any insight into how SOC acts on the low-energy states described by the Wannier orbitals and requires detailed first principles calculations for each new material studied. A more chemically intuitive description should consider the coupling between the electronic spin and the molecular orbital degrees of freedom.

Symmetry dictates the general form of the SMOC.^{6} For cyclic molecules with *C*_{N} symmetry the coupling is^{6}

where $cm\sigma \u2020=i|m|\u2211jaj\sigma \u2020ei2\pi j/N/N$. For a C_{3} symmetric molecule, such as Mo_{3}S_{7}(dmit)_{3}, this simplifies to

where $L\alpha x,y,z$ are the Cartesian components of the one-electron molecular orbital “angular momentum” operator for the *α*^{th} electron, *λ*_{xy}, is the transverse SMOC and *λ*_{z} the longitudinal SMOC.

The SMOC hamiltonian (6) gives a far clearer understanding of the relevant low-energy physics than the atomistic description (4). This is because molecular orbitals, not atomic orbitals, are the natural basis for the discussion of molecular physics. As molecular orbitals for cyclic molecules are Bloch states of the Wanniers, *H*_{SMOC} describes the coupling of orbital currents running around the plane of the molecule to the electron’s spin. The low-energy model constructed from Wannier spinors from four component relativistic band structure calculations^{9} of Mo_{3}S_{7}(dmit)_{3} yields single molecule SMOC exactly as predicted by Eq. (6).

Finally, the Hubbard-Heisenberg contribution describing the Coulomb and exchange energy in each triangular cluster is:

where *U* is the onsite Hubbard interaction, *J*_{F} is an intracluster exchange interaction, and $ni\sigma =ai\sigma \u2020ai\sigma $ the number operator. The *direct* ferromagnetic exchange, *J*_{F} < 0, plays a crucial role in generating magnetic anisotropies.^{4}

The nearest-neighbor triangular clusters are connected through the hopping amplitude, *t*:

where ⟨*lm*⟩ denotes two nearest-neighbor triangular clusters and *i* a cluster site in a given molecule.

### B. Effective spin exchange model

By straightforward diagonalization of the full hamiltonian (1), we have found that in the strong coupling limit, *U* ≫ *t*, isolated triangular clusters with four electrons effectively behave as pseudospin-1 localized moments.^{4} The many-body states of the cluster can be classified according to the *z*-component of the total angular momentum, *J*_{z} = *L*_{z} + *S*_{z}. The ground state of the isolated cluster with no SMOC is a triplet. As shown in Fig.2, SMOC splits the lowest energy triplet into a non-degenerate singlet (*j* = 0) which is the ground state and a doublet (*j* = ±1), where *j* denotes the eigenstates of *J*_{z}. Since we have an even number of electrons in the cluster, Kramers theorem does not apply and non-degenerate states are possible. SMOC can be effectively accounted for as a single-spin anisotropy contribution at each cluster. Hence, the effective spin model for the *m*-th Mo_{3}S_{7}(dmit)_{3} molecule in the crystal is just:

where $Srm$ describes the effective pseudospin-1 localized at a triangular cluster. The single-spin anisotropy, *D*, depends strongly on SMOC as shown in Fig. 3. Using perturbation theory^{4,5} we have obtained an effective spin exchange model describing the magnetic exchange coupling between the pseudospin-1 in the *a* − *b* planes of Mo_{3}S_{7}(dmit)_{3}. This model is an XXZ+120^{0} honeycomb quantum compass model^{4} leading to:

where *J* = (*J*_{xx} + *J*_{yy})/2, Δ = *J*_{zz}/*J* and *Q* = (*J*_{xx} − *J*_{yy})/2 and *ϕ*_{i} = 2*π*(*i* − 1)/3 where *i* labels the three bonds around each ▽ shown in Fig. 1. Thus, we see that the second term (proportional to *J*) is simply the XXZ model and the third term (proportional to *Q*) is the honeycomb 120° compass model.^{10}

The exchange couplings and single-spin anisotropy entering model (10) have been obtained numerically from perturbation theory^{4} as well as analytically using a canonical transformation.^{5} In Fig. 3 we show the dependence of the model (10) parameters: *J*, *Q*, Δ and *D* with *λ*_{xy} for the longitudinal SMOC, *λ*_{z} = *λ*_{xy}/2, appropriate for Mo_{3}S_{7}(dmit)_{3} crystals. With no SMOC, the exchange couplings are isotropic: *J* ≈ 0.04*t*_{c}, *Q* = 0, Δ = 1, *J*_{xz} = 0, as it should. As SMOC increases, spin exchange anisotropy, *i. e.*, *Q* ≠ 0, Δ ≠ 1 (or equivalently *J*_{xx} ≠ *J*_{yy} ≠ *J*_{zz}) increases. On the other hand, Fig. 3 shows how the single-spin anisotropy is rapidly enhanced by SMOC.

## III. LARGE-*D* PHASE INDUCED BY SMOC

One can see from the results of Fig. 3 that, for sufficiently large SMOC, the single-spin anisotropy overcomes the exchange couplings, *D* > *J*, Δ*J* = *J*_{zz}. We can expect that in the limit, *D* ≫ *J*, *J*_{zz}, the model is dominated by single-spin anisotropy and the ground state of the model is just the tensor product of $Sr\u2113z=0$ at each lattice site. Hence, it is important to first obtain an estimate of the critical *D*_{c} at which a transition to a large-*D* phase occurs and the magnitude of SMOC needed to reach such critical *D*. In order to do so, we explore below the magnetic properties of a slightly simplified version of the full model (10) in the limit of weak SMOC. This limit is relevant to the *a* − *b* planes of Mo_{3}S_{7}(dmit)_{3}, for which *ab initio* calculations^{9} give $\lambda xyMo3S7(dmit)3=0.042tc$ with *λ*_{z} = *λ*_{xy}/2. We can see from Fig. 3 that for *λ*_{xy} ≲ 0.5 exchange coupling anisotropies are negligible which justifies using an isotropic Heisenberg model for describing the magnetic properties of the *a* − *b* planes of the crystal.

We have recently explored the Néel to large-*D* transition of the isotropic version of model (10) using a *SU*(3) description of the spins which introduces three Schwinger bosons to properly account for the three projections: $Siz=0,\xb11$ of the spin-1 at each lattice site.^{11} We find a transition from the Néel ordered phase to a large-*D* phase for *D*_{c} ≈ 3*J*. From this condition and the dependence of *D* with *λ*_{xy} we can obtain a realistic estimate of SMOC needed for the large-*D* phase to occur in isolated *a* − *b* layers of Mo_{3}S_{7}(dmit)_{3}. In Fig. 3 we show the critical SMOC $(\lambda xy/tc)critic\u22480.38$ at which the *D*_{c} = 3*J* condition is satisfied. This critical value is about nine times larger than SMOC in Mo_{3}S_{7}(dmit)_{3} crystals. Hence, we conclude that isolated *a* − *b* layers of Mo_{3}S_{7}(dmit)_{3} should be in a magnetically ordered Néel state.

## IV. MAGNETIC FRUSTRATION EFFECTS

Magnetic frustration plays an important role in inducing disordered spin liquid states. We analyze the effect of a next-nearest neighbor exchange coupling, *J*′, on the magnetic properties of model (10) in the limit of weak SMOC *i. e.* with isotropic couplings. Hence, the model we study is a *S* = 1 *J* − *J*′ Heisenberg on the honeycomb lattice with single-spin anisotropy:

We first discuss the phase diagram of this model with no single-spin anisotropy, *D* = 0, based on recent SU(2) Schwinger boson mean-field theory (SBMFT) calculations.^{11} For *J*′ = 0, the ground state of the model is a Néel state, as expected since the lattice is bipartite. A direct transition from Néel order to spiral order occurs for *J*′ ≈ 0.24.

We finally discuss the effect of *J*′ on the Néel to large-*D* transition of model (11) based on SU(3) SBMFT calculations of model (11). In Fig. 4 we show the critical *D* at which a spin gap opens up signaling the formation of the paramagnetic large-*D* phase. We find that the critical *D*_{c} = 3*J* obtained for *J*′ = 0 of the unfrustrated model, is strongly suppressed by *J*′. Hence, magnetic frustration favors the large-*D* phase. The kink around *J*′ = 0.12*J* in Fig. 4 is related to the expected poor performance of the SU(3) SBMFT approach as *D*_{c} is suppressed below *J*, *D*_{c} < *J*.^{11}

## V. CONCLUSIONS

The magnetic properties of the *a* − *b* planes of Mo_{3}S_{7}(dmit)_{3} crystals can be modeled through a *S* = 1, XXZ+120^{0} quantum compass model with single-spin anisotropy.^{4} For the weak SMOC present in Mo_{3}S_{7}(dmit)_{3} crystals,^{9} quantum compass couplings and single-spin anisotropy play a minor role so that the effective model in this regime is just the isotropic *S* = 1 Heisenberg model on the honeycomb lattice. Hence, we predict that isolated honeycomb layers of Mo_{3}S_{7}(dmit)_{3} are Néel ordered. However, we find that magnetic frustration induced by a next-nearest neighbor exchange coupling, *J*′, favors the large-*D* phase.

## ACKNOWLEDGMENTS

J. M. acknowledges financial support from (MAT2015-66128-R) MINECO/FEDER, UE, and A.R. from the French National Research Agency (Contract ORGANISO, ANR-15-CE09- 0017). Work at the University of Queensland was supported by the Australian Research Council (DP160100060).