Universal excitonic superexchange in spin-orbit-coupled Mott insulators

We point out the universal presence of the excitonic superexchange in spin-orbit-coupled Mott insulators. It is observed that, the restriction to the lowest spin-orbit-entangled"$J$"states may sometimes be insufficient to characterize the microscopic physics, and the virtual excitonic processes via the upper"$J$"states provide an important correction to the superexchange. We illustrate this excitonic superexchange from a two-dimensional $5d$ iridate Sr$_2$IrO$_4$ and explain its physical consequences such as the orbital-like coupling to the external magnetic flux and the nonlinear magnetic susceptibility. The universal presence of the excitonic superexchange in other spin-orbit-coupled Mott insulators such as $3d$ Co-based Kitaev magnets and even $f$ electron rare-earth magnets is further discussed.

orbitals split into the t 2g and the e g manifolds, with the former having lower energy. In the case of a large crystal field, the energy gap is large, leading to a low-spin d 5 configuration. There is one hole in the t 2g manifold which consists of orbitals d yz , d zx , and d xy . The t 2g manifold has an effective orbital angular momentum l eff = 1 [37]. Because there is a single hole in the atomic ground state, we will use the hole representation in this paper. The atomic Hamiltonian in the t 2g subspace at site i is where c † imσ (c imσ ) is the hole creation (annihilation) operator with orbital m = 1, 2, 3 and spin σ =↑, ↓ at site i. The first term is the Hubbard interaction with strength U > 0. The prime of the summation means that the term with (mσ) = (m σ ) is excluded. Note that a hole creation (annihilation) operator is an electron annihilation (creation) operator. Thus, the electron occupation number n (e) imσ = 1 − n imσ in the hole representation. The Hund's interactions are relatively small and hence neglected. The last term is the spin-orbit interaction with strength λ > 0. In the hole representation, the sign before λ is negative. The orbital index m = 1, 2, 3 corresponds to the three t 2g orbitals d yz , d zx , and d xy , respectively. In this index ordering, the orbital angular momentum matrix elements (l m ) m m = i mm m , in which l 1,2,3 refers to l x,y,z , respectively. The sign is opposite to the matrix elements of the genuine l = 1 orbitals [37]. The spin angular momentum s = σ/2, where σ = (σ x , σ y , σ z ) is the vector of Pauli matrices.
The eigenstates of the SOC term are two j eff = 1/2 states with energy −λ, and four j eff = 3/2 states with energy λ/2. In the atomic limit, the hole lies on the two-fold degenerate j eff = 1/2 level, making the ion an effective spin-1/2 system. As long as the hoppings between different sites remain relatively weak as compared with the on-site Coulomb repulsion, the system remains a spin-orbit-coupled Mott insulator.
Exchange interaction.-The interactions between the effective spins has their origin in the virtual hopping processes. When the hopping terms between different sites are taken into account, there are virtual intermediate states that contain more than one hole on the same ion site, at the cost of Coulomb repulsion energy. The intermediate states are subject to the Pauli exclusion principle, which prevents two holes with the same spin from occupying the same state. Therefore, the resulting energy correction will be dependent on the spin configurations of the initial and finial states. The overall effect is the development of an effective interaction between neighboring spins, which is called exchange interaction.
In general, the hopping processes can be described by the Hamiltonian in which h i j = h † ji is the matrix for the hopping from site j to site i, and where the Latin indices m and n refer to the ground states of the unperturbed Hamiltonian H 0 , and the Greek indices α and β refer to the excited states. ω 0α is the difference between the unperturbed energies of the ground state 0 and the excited state α.
We begin by a simple model shown in Fig. 1(a), in which the d 5 octahedrally coordinated transition metals arrange in a two-dimensional square lattice, the same as the CuO 2 plane in cuprates.
By symmetry considerations, the hopping matrices for 2 → 1 and 3 → 1 in Fig. 1(a) have the form respectively. Their subscripts hint that the hoppings are similar to the corresponding dd-type Slater-Koster parameters, although usually the ligand-mediated indirect hopping is the dominant mechanism. All the other hopping matrices up to second nearest neighbor can be derived by symmetry. The hopping terms for further neighbors are neglected. Traditionally, the excited states are limited to the ones that have two holes lying on the same j eff = 1/2 manifold, similar to the one-band Hubbard model [38]. With this restriction, only the second-order perturbation term in Eq. (3) contributes to the nearest-neighbor spin-spin interaction. By the definition of the spin operator at site i as (S i ) ss = c † is σ ss c is /2 where s and s refer to the two j eff = 1/2 states, the effective Hamiltonian for the spin interaction between site 1 and site 2 is derived as This is similar to the result for the one-band Hubbard model. If t π = t δ = t, then the coefficient 4t 2 /U is the expected antiferromagnetic Heisenberg exchange interaction parameter.
Virtual excitonic process.-In spin-orbit-coupled Mott insulator the j eff = 1/2 level is not well separated from the j eff = 3/2 level, as can be seen in the typical spin-orbitcoupled Mott insulators Sr 2 IrO 4 [39], Na 2 IrO 3 [40], and α-RuCl 3 [12]. The existence of the nearby j eff = 3/2 level will have an impact on the spin exchange interaction.
After taking the j eff = 3/2 level into consideration, the intermediate excited states represented by the indices α and β in Eq. (3) not only contain those in which two holes occupy the j eff = 1/2 level, but also those in which one hole occupies the j eff = 1/2 level and the other hole occupies the j eff = 3/2 level 1 . The process involving the latter type of excited states is called the virtual excitonic process, because it is as if one hole is excited from the j eff = 1/2 level to the j eff = 3/2 level, leaving an electron behind, and the virtual pair of the electron and the excited hole forms a virtual exciton.
Allowing the virtual excitonic process in the perturbation expansion (3), it is found that the correction to the two-spin exchange interaction comes from the third order term, which involves the participation of a third ion, as schematically shown in Fig. 2.
The involvement of a third ion has significant consequences when an out-of-plane magnetic field is applied. There will be additional Peierls phases in the hopping matrices. The accumulated Peierls phase for a closed loop should equal to 2π times the magnetic flux through the loop area divided by the flux quantum hc/e. Assuming that the accumulated Peierls phase for a plaquette is φ, it is convenient to use the gauge indicated in Fig. 1(b). Note that the second order interaction (5) is independent of the magnetic field because the areas of the hopping paths 1 → 2 → 1 and 2 → 1 → 2 are zero.
For a square lattice up to second-nearest-neighbor hopping, Fig. 1(c) shows the four hopping paths at the third perturbation order that contribute to the exchange interaction between site 1 and site 2. The resulting correction term is found to be In sharp contrast to Eq. (5), the correction term is dependent on the magnetic field. This is because the area of the third order hopping processes like 1 → 2 → 3 → 1 is nonzero, leading to a nonvanishing magnetic flux. Such field dependence is fundamentally different from the Zeeman interaction, because it is essentially coupled to the magnetic flux, which is an orbital effect. The above analysis can be easily generalized to different lattice models, so the effect should be ubiquitous in various spin-orbit-coupled Mott insulators.
At the third order perturbation, there are also the ring exchange interactions among each group of three neighboring spins. For the spins 1, 2, and 3 in Fig. 1(a), the ring exchange interaction is Because the mixed product of three spins breaks the timereversal symmetry, the 3-spin ring exchange exists only when a magnetic field is applied. Hence, it must be field-dependent. The result (7) does not involve the virtual excitonic process. The hole only hops between the j eff = 1/2 states of the three ions. The corrections by the virtual excitonic process comes at the fourth perturbation order, which engages a fourth ion to provide the j eff = 3/2 states. We expect that the correction terms have a field dependence of sin φ, because the loop area is that of a plaquette, and the product sin φ S 1 · (S 2 × S 3 ) is time-reversal invariant.
Application to Sr 2 IrO 4 .-Having demonstrated the basic concept of the virtual excitonic process and the results on a simple model, we take a look at its effects in real material.
To demonstrate it, we use Sr 2 IrO 4 as an example. Sr 2 IrO 4 is a well-known spin-orbit-coupled insulator [3,38,41]. The electron configuration of the Ir 4+ ion is 5d 5 , the same as our assumption before. It consists of corner sharing IrO 6 octahedra layers with intervening Sr layers, see Fig. 3(a). The crystal structure of Sr 2 IrO 4 is the Ca 2 MnO 4 type with space group I4 1 /acd [42]. Compared to the above simple model, the complication is that the IrO 6 octahedra rotate alternately from perfect alignment, as shown in Fig. 3(b).
The TB model without SOC is still given by the form of Eq. (2). But we note that here the d orbitals in the subscripts are defined with respect to the local coordinate system aligned with the IrO 6 octahedron at the corresponding site, as shown by the two sets of local coordinate systems x y z in Fig. 3(b). On the contrary, the spin directions are defined with respect to the global xyz axes in Fig. 3(b), because in the end we would like to have a spin interaction Hamiltonian expressed in the global coordinate system.
Considering the crystal symmetry, the hole hopping matrices between the nearest-neighbor pair 12 and the secondnearest-neighbor pair 13 in Fig. 3(b) are respectively. The overall minus signs are to emphasize that the hopping matrices for holes differ from those for electrons by a sign, while the latter is the direct result of the first-principles calculations. We note that the zero matrix elements in h 12 and h 13 are not strict because of the interlayer interactions. However, their magnitudes are around 0.1 meV, much smaller than the other matrix elements. Thus, the smallness of these matrix elements reflects the weakness of the interlayer interactions. Other hopping matrices can be inferred by the crystal symmetry. The hopping parameters t 1,2,3,4 and t 1,2,3,4 , the SOC strength λ, and the Hubbard interaction parameter U are listed in Table I. By applying the perturbation theory, the spin interaction between the pair 12 in Fig. 3

(b) is found to be
The analytical expressions for the coefficients are given in the SM, and their numerical values are listed in Table I. The terms involving J , D , and K come from the virtual excitonic processes. The phase φ = eBl 2 / , where B is the out-of-plane magnetic field, and l = 3.88Å is the distance between nearestneighbor Ir atoms. The order of magnitude φ ∼ 1 corresponds to a magnetic field B ∼ 4.4 × 10 3 T. Nonlinear susceptibility.-The positive coefficient before S z 1 S z 2 leads to the easy-plane anisotropy. The in-plane magnetic structure is canted antiferromagnetic because of the competition between the Heisenberg and the Dzyaloshinskii-Moriya interaction (the first two terms of Eq. (9)). With the application of a magnetic field B perpendicular to the IrO 2 layers, the canted antiferromagnetic structure develops an outof-plane component. We first neglect the virtual excitonic process corrections, then the classical energy per site is where S B and S B are the spins on the two sublattices, with the same z component, and µ B is the Bohr magneton. We have used the fact that the Landé g-factor of the j eff = 1/2 hole state is −2. By minimizing the energy, the tilting angle of the spins away from the xy plane can be obtained, and the perpendicular magnetic susceptibility is found to be , and the parameters of the spin interaction between site 1 and site 2 in Fig. 3(b) (bottom), in meV.
in which N is the site density. This result (11) has the same property as the mean-field result for the Néel antiferromagnetism that below the Néel temperature it depends neither on the temperature nor on the magnetic field. Now we take into account the correction term from the virtual excitonic process.
It amounts to the replacements J → J − J cos(φ/2), D → D + D cos(φ/2), and K → K + K cos(φ/2). With Taylor expansion up to O(φ 2 ), we find the magnetization where is linear susceptibility, and is the nonlinear susceptibility [55,56]. Its existence is the consequence of the virtual excitonic process. This clarifies one important source of the nonlinear susceptibility in these systems.
Discussion.-Here we discuss the applicability of the excitonic superexchange. We have shown that, the virtual excitonic process appears in the high order perturbation theory of the multiple-band Hubbard model. A smaller spin-orbitcoupling and/or a weak Hubbard interaction could enhance its contribution. Apparently, in many of these 4d/5d magnets, the spin-orbit coupling is not really the dominant energy scale, and the systems also behave like weak Mott insulators. The 4d magnet α-RuCl 3 is a good example of this kind. This material shows a remarkable thermal Hall transport result in external magnetic fields, which is likely to be related to the gapped Kitaev spin liquid [57][58][59]. More recently, it has been found that its thermal conductivity oscillates when an external magnetic field is applied [60]. It is believed that the oscillation originates from the existence of the spinon Fermi surface of some kind. The orbital quantization of the spinons requires the coupling of the spinons with an internal U(1) gauge field [61]. Such orbital effect is obviously not controlled via the Zeeman term. The virtual excitonic process and/or the strong charge fluctuation could provide the microscopic means to lock the internal U(1) gauge field with the external magnetic flux, and thereby generating the oscillatory behaviors.
For the Co-based Kitaev materials [22][23][24][25][26], the Hubbard interaction is large, but the spin-orbit coupling is weak. Thus, the excitonic superexchange may still be important. For the rare-earth magnets, the crystal field enters as an important energy scale. If the crystal field energy separation between the ground states and the excited states is not large enough compared to the exchange energy scale, the excited states would be involved into the low-temperature magnetic physics. This upper-branch magnetism has been invoked for the pyrochlore magnet Tb 2 Ti 2 O 7 [62][63][64][65].
The concept could also be applied to systems such as metal organic frameworks, molecular magnets, and Moiré systems, which have large lattice constants. Another example is the spin-liquid candidate material 1T -TaS 2 [66][67][68][69][70], which enters into a charge-density-wave (CDW) phase at low temperature. In the CDW phase, each supercell has a localized unpaired electron which provides the spin. The neighboring spins are separated at a distance l larger than ten angstroms. Moreover, the orbitals d xy , d x 2 −y 2 , and d z 2 have similar energies [71], which facilitates the excitonic process. The required magnetic field B for the phase φ ∼ eBl 2 / to be at the order of O(1) drops to hundreds of Tesla, making the effect of excitonic process more prominent.
Acknowledgments.-The first-principles calculations for this work were performed on TianHe-2. We are thankful for the support from the National Supercomputing Cen  [1,2] are constructed, and the tight-binding (TB) model is established. The firstprinciples calculations are based on density functional theory [3,4] and plane-wave pseudopotential method, as implemented in the Quantum ESPRESSO package [5,6]. The low temperature experimental crystal structure [7] is adopted in the calculations. The energy cutoff of the plane wave basis set is 120 Ry. The Brillouin zone is sampled by a 6 × 6 × 6 k-grid. The Schlipf-Gygi norm-conserving pseudopotentials [8] produced by the ONCVPSP code [9] are used. The ML-WFs and the TB model are constructed by the RESPACK program [10]. The initial guesses of the Wannier functions (WFs) are the projections of the KS states near the Fermi level to the local t 2g orbitals of each Ir atom. The spatial spread of the WFs are then minimized.  The comparison of the band structures calculated by the TB model and DFT without spin-orbit coupling (SOC) is presented in Fig. S2(a). The high symmetry paths in the Brillouin zone shown in Fig. S1 are generated by the SeeK-path tool [11,12]. It can be seen that the tight binding model reproduces the DFT results perfectly. The nearest-neighbor and next-nearest-neighbor hopping parameters are listed in Table I of the main text. We then add to the obtained TB model a local SOC term (the second term of Eq. (1) in the main text), and compare its band structure with the DFT result including SOC. The SOC strength λ is fitted by the least square method to be λ = 0.36 eV. The corresponding band structures are shown in Fig. S2(b). The details near the Fermi level (E = 0) are well reproduced by the TB model with the local SOC term.
To evaluate the Hubbard parameter U, constrained random phase approximation [13] is adopted, as implemented in the RESPACK program [10]. The Hilbert space is divided into the subspace of the t 2g orbitals in which the MLWFs are constructed, and the rest subspace (r subspace). U is taken to be the orbital-averaged on-site static Coulomb interaction partially screened by the polarization processes involving the r subspace. It is found that U = 2.232 eV.

II. ANALYTICAL EXPRESSIONS FOR THE SPIN INTERACTION PARAMETERS
The analytical expressions for the spin interaction parameters of Sr 2 IrO 4 in the main text are J = 4 (t 1 + t 2 + t 3 ) 2 − 4t 2 . The parameters of the spin interaction between site 1 and site 2 in Fig. 3(b) of the main text, in meV. The nonvanishing tetragonal crystal field is taken into account.
where the parameter The corresponding numerical values are listed in Table I of the main text.

III. EFFECTS OF THE TETRAGONAL CRYSTAL FIELD
In deriving the above analytical expressions, we have assumed a perfect octahedral environment of the Ir 4+ ion. As a result, the three t 2g orbitals are degenerate. In the crystal of Sr 2 IrO 4 , their is a tetragonal crystal field which makes the energy of the d xy orbital different from that of the degenerate d xz and d yz orbitals. It gives rise to an additional term in the atomic Hamiltonian where C † iσ = [c † i,yz,σ , c † i,zx,σ , c † i,xy,σ ] as in the main text. It is found that ∆ = 159 meV. There are no simple analytical expressions for the spin interaction parameters for nonvanishing ∆. The numerical results are listed in Table S1. Compared with the corresponding parameters in the Table I of the main text, it can be seen that the tetragonal crystal field has only minor influences on the spin interactions.