Exchange interactions, Jahn-Teller coupling, and multipole orders in pseudospin one-half $5\boldsymbol{d}^\mathbf{2}$ Mott insulators

We develop a microscopic theory of multipole interactions and orderings in 5$d^2$ transition metal ion compounds. In a cubic environment, the ground state of 5$d^2$ ions is a non-Kramers $E_g$ doublet, which is nonmagnetic but hosts quadrupole and octupole moments. We derive pseudospin one-half Hamiltonians describing various spin-orbital exchange processes between these ions. Direct overlap of the $t_{2g}$ orbitals results in bond-dependent pseudospin interactions similar to those for $e_g$ orbitals in manganites. The superexchange process via oxygen ions generates new types of pairwise interactions. In perovskites with 180$^\circ$ bonding, we find nearly equal mixture of Heisenberg and $e_g$ orbital compass couplings. The 90$^\circ$ superexchange in compounds with edge-shared octahedra is most unusual: despite highly anisotropic shapes of the $E_g$ wavefunctions, the pseudospin interactions have no bond dependence and show instead a hidden SU(2) symmetry, which equally supports quadrupole and octupole orders. We consider the $E_g$ pseudospin models on various lattices and obtain their ground state properties using analytical, classical Monte Carlo, and exact diagonalization methods. On the honeycomb lattice, we observe a duality with the extended Kitaev model, and uncover a critical point where the quadrupole and octupole states are exactly degenerate. On the triangular lattice, an exotic pseudospin state, corresponding to the coherent superposition of vortex-type quadrupole and ferri-type octupole orders, is realized due to geometrical frustration. We also consider Jahn-Teller coupling effects and lattice mediated interactions between $E_g$ pseudospins. Possible implications of the results for recent experiments on double perovskite osmates are discussed, including effects of local distortions on the pseudospin wavefunctions and interactions.


I. INTRODUCTION
As a hallmark of strong correlations, the spin-orbital multiplet structure of ions is largely preserved in transition metal (TM) compounds. At low temperatures, the spin and orbital degeneracy of these multiplet levels has to be lifted one way or another. Apart from exotic means of the entropy quenching such as formation of quantum spin and orbital liquids, this is typically done by longrange ordering of spins and orbitals, or their composites, through symmetry breaking phase transitions.
Broadly speaking, the interactions driving these phase transitions have three different microscopic origins: (a) Jahn-Teller orbital-lattice coupling, (b) Kugel-Khomskii type spin-orbital exchange, and (c) relativistic spin-orbit coupling (SOC). Depending on the multiplet structure of constituent ions and the nature of chemical bonds in a crystal, the interplay between these couplings may take various forms, resulting in rich spin-orbital physics in TM compounds.
Lifting the orbital degeneracy via a cooperative Jahn-Teller (JT) structural transition is most common in e g orbital systems like manganites. At this transition, the orbitals are (self)trapped by static lattice distortions. The JT driven orbital order is essentially independent of spins and happens well before magnetic ordering. In this picture, the low energy physics is given by "spin-only" Hamiltonians, with the exchange parameters dictated by the Goodenough-Kanamori rules [1,2].
In t 2g orbital systems with relatively weak JT coupling, the spin S and orbital L degrees of freedom are no longer separated [3]. They may instead develop joint dynamics driven by the spin-orbital exchange interactions, as well as by intraionic SOC which unifies the two sectors by forming total angular momentum J = S + L. The latter root to the "spin-orbital-entangled" physics is especially relevant to 4d and 5d electron compounds.
In the strong SOC limit, the JT orbital-lattice coupling and Kugel-Khomskii exchange interactions have to be reformulated in terms of total angular momentum J of the lowest multiplet level, as is usually done in 4f electron systems. This leads to a number of important consequences. First, the JT orbital order is "converted" into quadrupole order of J moments, involving also the spin sector which was initially "blind" to JT physics. Effective JT coupling is typically reduced, due to a partial suppression of the initial orbital degeneracy. Second, exchange interactions between effective J moments ("pseudospins") may become highly anisotropic and bond-directional; this is due to the nonspherical shape of the spin-orbit entangled wavefunctions. Third, pseudospin states may carry not only dipole or quadrupole moments, but also higher-rank multipoles 2 such as a magnetic octupole.
The physical content of pseudospin wavefunctions is decided by a filling factor n of d-orbital levels. In combination with the lattice and chemical bonding geometry in a given material, this leads to a variety of non-trivial interactions and ground states among different d n compounds. This includes a possible realization of Kitaev spin-liquids, excitonic magnetism, and multipole orders (for a recent review, see Ref. [4]).
In this paper, we focus on spin-orbital physics in compounds based on d 2 ions. The d 2 configuration with twoelectron spin S = 1 and effective orbital moment L = 1 is special, because its total angular momentum J = 2 is isomorphic to a single d-electron orbital moment, l = 2. This analogy has interesting implications for the symmetry and physical properties of d 2 ions. Namely, in a cubic environment, a J = 2 level has to split into E g doublet and T 2g triplet levels [see Fig. 1(a)], just like the d-electron l = 2 level splits into e g and t 2g orbital levels [5]. While the T 2g triplet hosts an effective angular momentum J = 1 (with a familiar relation J = −J ) [6], the non-Kramers E g doublet is similar to an e g doublet and carriers no dipole moment. This implies that d 2 ions with non-Kramers E g ground states may show high-rank multipole orders similar to rare-earth f 2 non-Kramers Γ 3 ions [7].
Experimentally, a single phase transition around 30-50 K is observed in 5d 2 double perovskite (DP) compounds [8][9][10][11]. This is very different from 5d 1 Kramers ion DPs which show two separate transitions [12][13][14][15], corresponding to quadrupole (structural) and dipole orders of J = 3/2 states [16,17]. Having a single transition is natural for pseudospin-1/2 doublet systems, and this clearly points to the E g doublet physics in 5d 2 DPs. However, the precise nature of this transition is not yet fully established. The structural changes at this transition, if any, are found to be below 0.1% [8]. While no magnetic Bragg peaks were seen in neutron diffraction data, timereversal (TR) symmetry breaking is detected by muon spin relaxation. To reconcile these observations, a ferrotype octupolar order of the E g doublets has been proposed [8,18,19].
The octupole is a third-rank magnetic multipole which carries no dipole moment, and its long-range order is observed in rare-earth compounds (see Ref. [7] for a review of multipole orders). The possibility of octupolar order in d electron systems is intriguing. It is actually quite unexpected because an E g doublet is subject to JT physics: its partners have different charge density shapes (planar and elongated), see Fig. 1(a). Therefore, a conventional quadrupole order like in e g orbital systems [2] is the most natural instability to expect in the first place. To realize the octupolar order, exchange interactions between the octupole moments must be strong enough to overcome the quadrupolar interactions contributed by the Kugel-Khomskii exchange and orbital-lattice JT couplings.
Early theoretical work [20] on d 2 DP systems with strong SOC assumed that cubic splitting of the J = 2 level ∆ c is smaller than the exchange couplings and therefore neglected it. The obtained phase diagram contains dipolar and quadrupolar ordered states. Here we develop a theory of d 2 electron systems starting from the opposite limit, i.e. when cubic splitting ∆ c is large and the E g doublet is well separated from the virtual T 2g states, as actually seen in experiment [8]. Having in mind 5d 2 materials other than DP compounds, we keep the discussion as general as possible, considering various spinorbital exchange processes typical in TM oxides. The resulting E g doublet interactions are represented in terms of pseudospin one-half Hamiltonians. In most cases, the interactions are dominated by quadrupolar couplings. In a 90 • exchange geometry however, the quadrupole and octupole channels are equally present, and effective interactions on a single bond can be written in a Heisenberg form with no preference for either of these two channels.
The resulting multipole orders of E g doublets in different lattices are considered. On a honeycomb lattice, we show that the E g pseudospin model can be mapped to the extended Kitaev model, thereby uncovering a hidden SU (2) symmetry point that separates quadrupole and octupole orders. The pseudospins on a geometrically frustrated triangular lattice show more complex phase behavior, including a coherent mixture of different rank (quadrupole and octupole) orders in the ground state. The order parameters are reduced by quantum fluctuations. In DP lattices, we find that the exchange interactions favor a quadrupole order.
We further discuss orbital-lattice coupling effects, and show that JT phonon mediated interactions cooperate with exchange interactions to support quadrupolar order. This is similar to conventional e g orbital systems. We suggest that in DP lattices, where the magnetic ions are widely separated and have no common oxygen, a dynamical Jahn-Teller effect may develop to reduce the structural distortions induced by quadrupolar order. We also consider modifications of the pseudospin wavefunctions by symmetry lowering distortions (caused by site disorder or other defects), and find that they induce a magnetic dipole moment on the E g doublet. In general, d 2 compounds represent an interesting class of materials where all three main actors -the electron exchange, orbital-lattice interaction, and relativistic SOC -play an essential role in determining the ground states and lowenergy excitations.
The paper is organized as follows: Sec. II introduces the E g doublet states and their pseudospin-1/2 description. In Section III, we derive pseudospin Hamiltonians considering different orbital exchange geometries which are typical in TM compounds. Sec. IV studies pseudospin orderings and excitations on various lattice structures. Sec. V discusses Jahn-Teller coupling and disorder effects in the context of experiments in DP compounds. Sec. VI summarizes the main results. The E g doublet wavefunctions written in the J z basis are [6]: 1 √ 2 (|2 + | − 2 ) and |0 . We regard them as pseudospin s = 1/2 states | ↑ and | ↓ , correspondingly.
To get an idea about the orbital shapes, one can represent these functions in terms of two-electron spin and orbital |S z , L z states: In the pseudospin-up state with L z = ±1, one of the electrons must occupy l z = 0 planar orbital d xy , flattening the overall charge density as shown in Fig. 1(a). However, the pseudospin-down state is dominated by an L z = 0 component, where the electrons occupy l z = 1 and l z = −1 complex orbitals ∓(d yz ± id zx )/ √ 2; thus, its charge density is elongated towards apical oxygen O z . Under cubic rotations, these wavefunctions transform in a standard way, similar to e g orbital pair, x 2 − y 2 and 3z 2 − r 2 .
Within the E g doublet, the J = 2 quadrupole operators have matrix elements ± 1 2 |O 3 | ± 1 2 = ±1 and ± 1 2 |O 2 | ∓ 1 2 = 1. Thus, the following correspondence between the pseudospin s z and s x components, and E g quadrupoles follows: s z = 1 2 O 3 and s x = 1 2 O 2 . The third component s y = 1 2 T xyz describes the octupolar moment T xyz = 1 √ 3 J x J y J z with threefold symmetry axis [111]. The projections of the octahedral x, y, z axes onto the twodimensional pseudospin (s z , s x ) plane [111] make 120 • angles between them, and the pseudospin s z axis is parallel to the octahedral z axis projection, see Fig. 1(a). This is the most natural choice, because s z is related to the O 3 quadrupole moment (3) directed along z axis. As we will see below, this also results in one-to-one correspondence between the exchange bond labels γ ∈ {x, y, z} and octahedral (x, y, z) axes. The basis rotations within the (s z , s x ) plane by φ = 2π/3 correspond to the cyclic permutations among J x , J y , J z . Finally, we note that the s z and s x operators are TR-even, while the s y octupole is TR-odd; this implies that the pairwise interactions of the type s z i s y j and s x i s y j are not allowed, unless TR symmetry is broken.
Following e g orbital pseudospin formalism [2,21], we introduce the following pseudospin combinations: Here, the pseudospin index γ = (z, x, y) also specifies the corresponding angles φ γ = (0, 2π/3, 4π/3). In essence, (τ γ ,τ γ ) play the role of (s z γ , s x γ ) operators defined in the rotated basis of pseudospin functions: Physically, τ x (τ y ) andτ x (τ y ) correspond to the quadrupolar operators of 3x 2 − r 2 (3y 2 − r 2 ) and y 2 − z 2 (z 2 − x 2 ) symmetries, respectively. The notations τ γ and τ γ are useful since one may derive the exchange Hamiltonian H (γ) for γ = z type bonds in terms of (s z , s x ) pair, and then restore H (γ) for all γ by simply replacing s z → τ γ and s x →τ γ . In perovskites with 180 • bonding, the z-type bond is parallel to the octahedral z axis; while in other cases, e.g. in a honeycomb lattice, the z-type bond is orthogonal to the octahedral z axis (a convention also used in the Kitaev model literature).
To derive pseudospin exchange interactions, one has to project Kugel-Khomskii type spin-orbital Hamiltonians -which are already known from previous worksonto the low-energy E g doublet subspace. We note that conventional e g orbital exchange interactions operate in the quadrupolar sector (s z , s x ) exclusively [2]. In contrast, we will see below that the E g "orbital" exchange may involve interactions between the octupole moments s y , as well; this is because the E g pseudospin states are spin-orbit entangled objects. Combined with the specific hopping geometry of t 2g orbitals, this results in a nontrivial structure of the E g interactions. We consider below some basic exchange processes which commonly appear in transition metal compounds.

III. PSEUDOSPIN EXCHANGE HAMILTONIANS
A. Single-orbital exchange: direct t2g orbital overlap We start with the simple case where one specific orbital is active on a given nearest-neighbor (NN) exchange bond. Two examples of single-orbital coupling are shown in Fig. 1(b): direct d xy orbital hopping on z-type bonds in the honeycomb lattice, and d xy orbital hopping in the ab plane of the DP lattice [20]. In this case, we expect that the E g exchange Hamiltonian is similar to that for e g orbitals in ferromagnetic manganites. Indeed, spinorbit E g and pure orbital e g states have the same (Γ 3 ) symmetry properties. Moreover, the Kugel-Khomskii e g exchange process also involves a single-orbital, specific to a given bond [2].
Octupolar order Quadrupolar order β α 2. (a) Shifts and splitting of the J = 0, 1, 2 levels of d 2 en considering mixing of the ground state t 2 2g configuraith t2geg states by virtue of SOC interaction. Focusing lowest levels, we find the originally five-fold degenerate states to split into Eg doublet and T2g triplet. Evalperturbatively for ⇣ ⌧ 10Dq, the splitting comes out rtional to ⇣ 2 /10Dq. (b) Tetragonal compression leads to reased repulsion of d-electrons from apical oxygens and r singles out "planar" states from Eg and T2g sets. The upolar moment hosted by the Eg doublet gets pinned ay. (c) Complex combinations of the Eg doublet states xpose the octupolar moment of "cubic" shape.
neously satisfy all the interacting bond directions. he spin-orbital entanglement, this property is inhery the pseudospin interactions. 8 reover, apart from the exchange interactions driven rtual electron hoppings, there are other contributo the orbital interactions, especially in d 1 and ses. These are mediated by coupling to the JTons, and by electrostatic multipolar interactions bed-orbitals on di↵erent sites. 9 In low-energy efe Hamiltonians, these interactions transform into pseudospin multipolar couplings, driving both structural and "spin-nematic" transitions breaking cubic symmetry in real and pseudospin spaces. The JT driven interactions are also important in d 4 systems, as they split the excited J = 1 levels and hence promote magnetic condensation. 10 In the case of d 5 systems with Kramersdegenerate J = 1/2 pseudospins, the e↵ective Hamiltonians are predominantly of exchange origin, as in usual spin-1/2 systems, although JT orbital-lattice coupling still shows up in fine details of pseudospin dynamics, in a form of pseudospin-lattice coupling. 10 In general, the low-energy pseudospin Hamiltonians may take various forms depending on the electron configuration d n and symmetry of the crystal structure. Sensitivity of orbital interactions to bonding geometry is a decisive factor shaping the form of the pseudospin Hamiltonians. We illustrate this considering spin-orbital exchange process in two di↵erent cases -when metaloxygen octahedra MeO 6 share the corners, and when they share the edges. These two cases are common in TM compounds and referred to as 180 and 90 bonding geometry, reflecting the approximate angle of the Me-O-Me bonds. For simplicity, we limit ourselves to the case of d 5 ions with wavefunctions [c.f. Fig. 3(a)]: 2g configuration with t2geg states by virtue of SOC interaction. Focusing on the lowest levels, we find the originally five-fold degenerate J = 2 states to split into Eg doublet and T2g triplet. Evaluated perturbatively for ⇣ ⌧ 10Dq, the splitting comes out proportional to ⇣ 2 /10Dq. (b) Tetragonal compression leads to an increased repulsion of d-electrons from apical oxygens and further singles out "planar" states from Eg and T2g sets. The quadrupolar moment hosted by the Eg doublet gets pinned this way. (c) Complex combinations of the Eg doublet states that expose the octupolar moment of "cubic" shape. multaneously satisfy all the interacting bond directions. Via the spin-orbital entanglement, this property is inherited by the pseudospin interactions. 8 Moreover, apart from the exchange interactions driven by virtual electron hoppings, there are other contributions to the orbital interactions, especially in d 1 and d 2 cases. These are mediated by coupling to the JTphonons, and by electrostatic multipolar interactions between d-orbitals on di↵erent sites. 9 In low-energy effective Hamiltonians, these interactions transform into pseudospin multipolar couplings, driving both structural and "spin-nematic" transitions breaking cubic symmetry in real and pseudospin spaces. The JT driven interactions are also important in d 4 systems, as they split the excited J = 1 levels and hence promote magnetic condensation. 10 In the case of d 5 systems with Kramersdegenerate J = 1/2 pseudospins, the e↵ective Hamiltonians are predominantly of exchange origin, as in usual spin-1/2 systems, although JT orbital-lattice coupling still shows up in fine details of pseudospin dynamics, in a form of pseudospin-lattice coupling. 10 In general, the low-energy pseudospin Hamiltonians may take various forms depending on the electron configuration d n and symmetry of the crystal structure. Sensitivity of orbital interactions to bonding geometry is a decisive factor shaping the form of the pseudospin Hamiltonians. We illustrate this considering spin-orbital exchange process in two di↵erent cases -when metaloxygen octahedra MeO 6 share the corners, and when they share the edges. These two cases are common in TM compounds and referred to as 180 and 90 bonding geometry, reflecting the approximate angle of the Me-O-Me bonds. For simplicity, we limit ourselves to the case of d 5 ions with wavefunctions [c.f. Fig. 3

(a)]:
|f"i = + sin # |0, "i cos # | + 1, #i, 2g configuration with t2geg states by virtue of SOC interaction. Focusing on the lowest levels, we find the originally five-fold degenerate J = 2 states to split into Eg doublet and T2g triplet. Evaluated perturbatively for ⇣ ⌧ 10Dq, the splitting comes out proportional to ⇣ 2 /10Dq. (b) Tetragonal compression leads to an increased repulsion of d-electrons from apical oxygens and further singles out "planar" states from Eg and T2g sets. The quadrupolar moment hosted by the Eg doublet gets pinned this way. (c) Complex combinations of the Eg doublet states that expose the octupolar moment of "cubic" shape. multaneously satisfy all the interacting bond directions. Via the spin-orbital entanglement, this property is inherited by the pseudospin interactions. 8 Moreover, apart from the exchange interactions driven by virtual electron hoppings, there are other contributions to the orbital interactions, especially in d 1 and d 2 cases. These are mediated by coupling to the JTphonons, and by electrostatic multipolar interactions between d-orbitals on di↵erent sites. 9 In low-energy effective Hamiltonians, these interactions transform into pseudospin multipolar couplings, driving both structural and "spin-nematic" transitions breaking cubic symmetry in real and pseudospin spaces. The JT driven interactions are also important in d 4 systems, as they split the excited J = 1 levels and hence promote magnetic condensation. 10 In the case of d 5 systems with Kramersdegenerate J = 1/2 pseudospins, the e↵ective Hamiltonians are predominantly of exchange origin, as in usual spin-1/2 systems, although JT orbital-lattice coupling still shows up in fine details of pseudospin dynamics, in a form of pseudospin-lattice coupling. 10 In general, the low-energy pseudospin Hamiltonians may take various forms depending on the electron configuration d n and symmetry of the crystal structure. Sensitivity of orbital interactions to bonding geometry is a decisive factor shaping the form of the pseudospin Hamiltonians. We illustrate this considering spin-orbital exchange process in two di↵erent cases -when metaloxygen octahedra MeO 6 share the corners, and when they share the edges. These two cases are common in TM compounds and referred to as 180 and 90 bonding geometry, reflecting the approximate angle of the Me-O-Me bonds. For simplicity, we limit ourselves to the case of d 5 ions with wavefunctions [c.f. Fig. 3(a)]: Octupolar order Quadrupolar order β α IG. 2. (a) Shifts and splitting of the J = 0, 1, 2 levels of d 2 on when considering mixing of the ground state t 2 2g configuraion with t2geg states by virtue of SOC interaction. Focusing n the lowest levels, we find the originally five-fold degenerate = 2 states to split into Eg doublet and T2g triplet. Evalated perturbatively for ⇣ ⌧ 10Dq, the splitting comes out roportional to ⇣ 2 /10Dq. (b) Tetragonal compression leads to n increased repulsion of d-electrons from apical oxygens and urther singles out "planar" states from Eg and T2g sets. The uadrupolar moment hosted by the Eg doublet gets pinned his way. (c) Complex combinations of the Eg doublet states hat expose the octupolar moment of "cubic" shape.
ultaneously satisfy all the interacting bond directions. ia the spin-orbital entanglement, this property is inherted by the pseudospin interactions. 8 Moreover, apart from the exchange interactions driven y virtual electron hoppings, there are other contribuions to the orbital interactions, especially in d 1 and 2 cases. These are mediated by coupling to the JThonons, and by electrostatic multipolar interactions beween d-orbitals on di↵erent sites. 9 In low-energy efective Hamiltonians, these interactions transform into configuration. The e↵ective angular momentum is indicated by the rotating arrow, spin by the color following the convention of Fig. 1. For both contributions, Lz and Sz sum up to Jz = +1/2. (b) Similar decomposition of the J = 0 two-hole ground state of d 4 . The total orbital angular momentum Lz and spin Sz may be combined in three ways here. The internal compensation in L and S creates a cubic-shaped object showing no spin polarization.
pseudospin multipolar couplings, driving both structural and "spin-nematic" transitions breaking cubic symmetry in real and pseudospin spaces. The JT driven interactions are also important in d 4 systems, as they split the excited J = 1 levels and hence promote magnetic condensation. 10 In the case of d 5 systems with Kramersdegenerate J = 1/2 pseudospins, the e↵ective Hamiltonians are predominantly of exchange origin, as in usual spin-1/2 systems, although JT orbital-lattice coupling still shows up in fine details of pseudospin dynamics, in a form of pseudospin-lattice coupling. 10 In general, the low-energy pseudospin Hamiltonians may take various forms depending on the electron configuration d n and symmetry of the crystal structure. Sensitivity of orbital interactions to bonding geometry is a decisive factor shaping the form of the pseudospin Hamiltonians. We illustrate this considering spin-orbital exchange process in two di↵erent cases -when metaloxygen octahedra MeO 6 share the corners, and when they share the edges. These two cases are common in TM compounds and referred to as 180 and 90 bonding geometry, reflecting the approximate angle of the Me-O-Me bonds. For simplicity, we limit ourselves to the case of d 5 ions with wavefunctions [c.f. Fig. 3(a)]: Cubic splitting c of J = 2 level, and the patial shapes of E g doublet wavefunctions. Right panel hows the pseudospin coordinate axes with respect to oxyen octahedra. (b) Direct hopping between xy orbitals in oneycomb (left) and DP (right) lattices, resulting in bondependent pseudospin ⌧ -interactions (10) between spin-orbit ntangled E g states. (c) Two-orbital superexchange via 180 e-O-Me bonding geometry. Hopping is orbital conserving: z $ xz and yz $ yz. This process results in pseudospin nteractions (13) comprising isotropic Heisenberg and bondependent compass-type couplings. (d) Two-orbital superexhange via 90 bonding geometry. Hopping interchanges the rbital labels: xz $ yz. This process leads to the interactions 15), which are anisotropic in pseudospin space but have no ond dependence.
ond is orthogonal to the octahedral z axis (a convention sed also in Kitaev model literature).
To derive pseudospin exchange interactions, one has o project Kugel-Khomskii type spin-orbital Hamiltonins -which are already khown from the previous works onto low-energy E g doublet subspace. We note that conventional e g orbital exchange interactions operate n quadrupolar sector (s z , s x ) exclusively [2]. In conrast, we will see below that the E g "orbital" exchange ay involve interactions between the octupole moments y as well; this is because the E g pseudospin states are pin-orbit entangled objects. Combined with the specific opping geometry of t 2g orbitals, this results in nontrival structure of the E g interactions. We consider below ome basic exchange processes which commonly appear n transition metal compounds.

III. PSEUDOSPIN EXCHANGE HAMILTONIANS
A. Single-orbital exchange: direct t 2g orbital overlap We start with simple case when one specific orbital is ctive on a given bond. Shown in Fig. 1(b) are two exmples of a such single-orbital exchange: via a direct d xy rbital hopping on z-type bonds of honeycomb lattice, or xy orbital hopping in the ab plane of DP lattice [19]. In his case, we expect the E g exchange Hamiltonian simlar to that for e g orbitals in ferromagnetic manganites. ndeed, spin-orbit E g and pure orbital e g states have the ame ( 3 ) symmetry properties, and Kugel-Khomskii e g xchange process also involves a single-orbital specific to given bond [2]. Neglecting Hund's coupling e↵ects in the intermediate tates, direct hopping t d (d † xy,i d xy,j + H.c.) gives the folowing exchange Hamiltonian, written in terms of spin = 1 and orbital L = 1 moments of d 2 ion [21]: isomorphic to a single d-electron orbital moment l = 2. This analogy has interesting implications for the symmetry and physical properties of d 2 ions. Namely, in cubic environment, J = 2 level has to split into E g doublet and T 2g triplet levels, see Fig. 1(a), just like the d-electron l = 2 level does split into e g and t 2g orbital levels [5]. While T 2g triplet hosts an e↵ective angular momentum e J = 1 (with a familiar relation e J = J ) [6], the non-Kramers E g doublet is similar to e g doublet and carriers no dipole moment. This implies that d 2 ions with non-Kramers E g ground states may show high-rank multipole orders similar to rare-earth f 2 non-Kramers 3 ions [7].
Experimentally, a single phase transition around 30-50 K is observed in 5d 2 double perovskite (DP) compounds [8][9][10][11]. This is very di↵erent from 5d 1 Kramers ion DPs which show two separate transitions [12][13][14], corresponding to quadrupole (structural) and dipole orders of J = 3/2 states [15,16]. Having a single transition is natural for pseudospin-1/2 doublet systems, and this clearly points to the E g doublet physics in 5d 2 DPs. The structural changes at this transition, if any, are found to be below 0.1% [8]. (For comparison, distortions about 0.1% are observed in a quadrupole ordered 5d 1 DP system [14]). The time-reversal (TR) symmetry breaking is detected by muon spin relaxation. To reconcile these observations, a ferro-type octupolar order of the E g doublets has been proposed [8,17,18].
The octupole is a third-rank magnetic multipole which carriers no dipole moment, and its long-range order is observed in rare-earth compounds (see Ref. [7] for a review of multipole orders). Possibility of octupole order in d electron systems is intriguing. It is actually quite unexpected because E g doublet is subject to JT physics: its partners have di↵erent charge density shapes (planar and elongated), see Fig. 1(a). Therefore, a conventional quadrupole order like in e g orbital systems [2] is the most natural instability to expect in a first place. To realize the octupolar order, exchange interactions between the octupole moments must be strong enough to overcome the quadrupolar interactions contributed by the Kugel-Khomskii exchange and orbital-lattice JT couplings.
Early theoretical work [19] on d 2 DP systems with strong SOC assumed that cubic splitting of J = 2 level c is smaller than the exchange couplings and therefore neglected it. The obtained phase diagram contains dipolar and quadrupolar ordered states. Here we develop a theory of d 2 electron systems starting from the opposite limit, i.e. when cubic splitting c is large and E g doublet is well separated from the virtual T 2g states, as actually seen in the experiment [8]. Having in mind the 5d 2 materials other than DP compounds, we keep discussion as general as possible, considering various spin-orbital exchange processes typical in TM oxides. The resulting E g doublet interactions are represented in terms of teraction on a single bond can be written in a Heisenberg form with no preference in either of these two channels. The multipole orders of E g doublets in di↵erent lattices are considered. On a honeycomb lattice, we show that the E g pseudospin model can be mapped to the extended Kitaev model, uncovering thereby a hidden SU(2) symmetry point that separates quadrupole and octupole orders. On a triangular and double-perovskite lattices, the exchange interactions favor a quadrupole order, with the order parameters reduced by quantum fluctuations of pseudospins.
We further discuss orbital-lattice coupling e↵ects, and show that JT phonon mediated interactions cooperate with exchange interactions to support quadrupole order. This is similar to a conventional e g orbital systems. We suggest that in DP lattices, where the magnetic ions are widely separated and have no common oxygen, a dynamical Jahn-Teller e↵ect may develop to reduce the structural distortions induced by quadrupole order. We also consider modifications of the pseudospin wavefunctions by low symmetry distortions (caused by site disorder or other defects), and find that they induce a magnetic dipole moment on E g doublet. In general, d 2 compounds represent interesting class of materials where all three main actors -the electron exchange, orbital-lattice interaction, and relativistic SOC -play an essential role in determining the ground states and low-energy excitations.
The paper is organized as follows: Sec. II introduces E g doublet states and their pseudospin-1/2 description. In Section III, we derive pseudospin Hamiltonians considering di↵erent orbital exchange geometries which are typical in TM compounds. Sec. IV studies pseudospin orderings on various lattice structures. Sec. V discusses Jahn-Teller coupling e↵ects in the context of experiments in double-perovskites. Sec. VI summarizes the main results.
In the pseudospin-up state with L z = ±1, one of the electrons must occupy l z = 0 planar orbital d xy , flattening For the x (y) bonds where the d yz (d zx ) orbital exchange is active, L z is replaced by L x (L y ). Projection of this Hamiltonian onto the E g subspace results in: with J τ = 4 9 t 2 d U , and τ γ given by Eq. (5). This interaction has the same structure as the Kugel-Khomskii e g orbital Hamiltonian, but with the reduced exchange constant due to a small fraction of the active orbital (e.g. d xy for the z bond) in the two-electron E g wavefunction.
Representative values of t d ∼ 0.1 − 0.2 eV and U ∼ 2 eV would give an energy scale of J τ ∼ 2 − 9 meV (the lower end is appropriate for DP lattice where d-ions are well separated and thus hopping t is small).
We note that there is a subtle difference between Eq. (10) and the Kugel-Khomskii e g exchange in perovskites [2]. In the latter, the active e g -orbital (say 3z 2 − r 2 on z bond) is quasi-one dimensional and bondoriented, thus enforcing pseudospins τ γ to be aligned along the interacting γ-bond directions (hence the name "pseudodipolar" or "compass" model [2]). In contrast, the τ γ quadrupoles in Eq. (10) try to avoid the bond directions; e.g. for z-type bond we have s z i s z j coupling but with the Ising s z axis being perpendicular to the z bond direction. Physically, the pseudospin orientation specifies the shape of the quadrupolar charge distribution and can be probed in the experiment. Formally however, the two models can be converted into each other by a 90 • rotation within the (s x , s z ) quadrupolar plane [i.e. replacing τ in Eq. (10) byτ ]. This point has to be kept in mind while comparing the present τ -model results with those in canonical compass model studies [23][24][25][26][27].

B. Two-orbital superexchange: 180 • bonding geometry
This case is typical for a metal-oxygen-metal (Me-O-Me) superexchange process in perovskites, see Fig. 1(c). On the z-type bond, two orbitals a = d yz and b = d zx equally contribute, and hopping is orbital-conserving: −t(a † iσ a jσ + b † iσ b jσ ). The spin-orbital Hamiltonian (at J H = 0) reads as [31]: where orbital operator for z type bond reads as Operators O (x) and O (y) for x and y bonds follow from cubic permutations among L x , L y , L z .
A projection of the above Hamiltonian onto the E g subspace gives with the exchange constant J = 2 3 t 2 U . In this equation, the Heisenberg term gives an equal coupling in quadrupolar (s x , s z ) and octupolar s y sectors. This term is not present in the Kugel-Khomskii e g orbital Hamiltonian, but is realized here because the E g orbitals have a complex internal structure, and they are made of t 2g orbitals with hopping rules different from those for e g orbitals.
The second bond-dependent τ γ term in Eq. (13) is a direct analogue of the Kugel-Khomskii e g orbital exchange. It operates only in the quadrupolar channel, thus disfavoring octupolar correlations. It is important to note that the easy axis orientations in this term exactly coincide with the bond directions, as dictated by the shapes of the active complex orbitals, e.g. ∓(d yz ± id zx )/ √ 2 orbitals having rotational symmetry around z bond (like 3z 2 −r 2 axial symmetry in e g models). Thus the E g pseudospins in the 180 • bonding geometry behave exactly as the e g orbital compasses do in cubic lattices [2], orienting themselves along the bond directions. This follows from a general observation that in case of axial symmetry, the spin-1/2 anisotropy term should have a dipole-dipole interaction form [32]. For the same reason, the compasslike τ -interaction also appears for non-Kramers doublets in f 2 electron system [33]. Due to differences between d and f orbital hopping geometries however, the isotropic term in Eq. (13) is not present in the f 2 case. On square or cubic lattices, we expect that the Hamiltonian (13) would have two-sublattice quadrupolar order, with alternating planar and elongated E g states, as selected by the anisotropic τ -term via order-from-disorder mechanism.

C. Two-orbital superexchange: 90 • bonding geometry
This process is typical for nearest-neighbor Me-O 2 -Me superexchange in delafossite derived structures with edge shared octahedra, see Fig. 1(d). On the z type bond, two orbitals a = d yz and b = d zx equally contribute again, but hopping is orbital non-conserving: −t(a † iσ b jσ + b † iσ a jσ ). The resulting spin-orbital Hamiltonian reads as in the 180 • case, see Eq. (11), but now with the modified orbital part [i.e. interchanging L x j ↔ L y j in Eq. (12)] [31]: The orbital non-conservation during the hoppings has dramatic consequences for the exchange symmetry, as observed previously in spin-orbit J = 1/2 [3,34] and J = 0 [31] systems. In the present non-Kramers E g doublet case, this results in the pseudospin Hamiltonian: U . This result is remarkable in several aspects. It has no γ-bond dependence, since the quadrupolar (s x , s z ) part is isotropic and does not change under the rotations (5) and (6), and the octupolar s y moment is not affected by C 3 rotations around [111] axis and is thus independent of γ as well. This is unlike the d 5 Kramers doublet case, where the cubic rotations affect all three components of the J = 1/2 vector, via cyclic permutations of its x, y, z components (see, e.g. Eq. (5.8) in Ref. [3]). Nevertheless, SOC results in strong exchange anisotropy: quadrupoles are ferro-correlated, while the octupolar components s y are coupled in an antiferro-fashion.
In bipartite (e.g. honeycomb) lattices, this anisotropic Hamiltonian can conveniently be converted into an AF Heisenberg form Js i · s j , by changing the sign of the s x and s z components on one of the sublattices; such hidden symmetries are common to spin-orbit pseudospin-1/2 Hamiltonians [3,35,36]. After this transformation, one observes an exact degeneracy between quadrupole and octupolar orderings, with the out-of-plane Goldstone mode representing a smooth rotation from one type order to the other one at no energy cost. Such exact degeneracy and coherent mixture of different (even/odd) rank order parameters and related gapless modes is rather unusual, but have previously been discussed in the context of t 2g orbital Hamiltonians, see Refs. [3,37] for details.

IV. PSEUDOSPIN ORDER: QUADRUPOLAR VERSUS OCTUPOLAR STATES
In this section, we discuss the phase behavior and excitations of the E g pseudospin models on different lattices.

A. Simple cubic lattice
The two-orbital 180 • -exchange Hamiltonian, Eq. (13), is applicable to perovskite lattices. As we already mentioned in that section, we expect a two-sublattice quadrupolar order in this case. This is conceptually similar to e g orbital order in 3d systems; the only difference is that the E g "orbitals" are spin-orbit coupled objects. Like in the e g orbital case, both exchange and JT couplings will contribute to the quadrupolar ordering, and they typically support each other. Formally, the Hamiltonian (13), comprising an AF Heisenberg interaction and an anisotropic compass-like terms is very similar to the model studied in Ref. [28]. Thus, its excitation spectrum should acquire a sizeable gap due to the order-by-disorder mechanism.

B. Honeycomb lattice
A honeycomb lattice is derived from the delafossite structure with edge-shared octahedra. In general, two different channels are operative in this case: direct hopping t d considered in Sec. III A, and indirect t superexchange via 90 • bonding considered in Sec. III C. It is known that for pseudospin J = 1/2 exchange in d 5 compounds, there is also a combination of these two processes (i.e. t times t d terms) resulting in the off-diagonal, socalled Γ interaction [38]. Interestingly, such a cross-term is absent in the present E g problem. So, the full Hamiltonian in honeycomb or triangular lattices is comprised of the bond-dependent τ -model (10), and the 90 • exchange J-interaction (15) which is also anisotropic, but bondindependent: Physically, both J τ and J are positive, and their ratio can be arbitrary. While the first term operates in the pure quadrupolar (s z , s x ) sector, J coupling equally supports AF octupolar and FM quadrupolar states. As noticed above, the J interaction is actually dual to the AF Heisenberg model (on bipartite lattices). A finite τterm breaks this symmetry and selects the ordering type. Since the (s z , s x ) part of the J term has a negative sign, a small admixture of positive J τ reduces the quadrupole interactions. As a result, two-sublattice staggered order of octupole moments s y is favored at J J τ . The ground state wavefunction is complex, ψ A/B = (|↑ ± i|↓ )/ √ 2, and has a cubic shape, see Fig. 2(c) of Ref. [4].
In the opposite limit J τ J, it is obvious that s y octupole order has to give way to ordering of the τquadrupoles that live in (s z , s x ) plane. The quadrupole order is TR invariant (i.e. the condensate wavefunction is real) but breaks cubic symmetry. The transition is of a spin-flop type: spins flop from the [111] direction into the honeycomb plane. In terms of the condensate wavefunction, this corresponds to the phase-jump from π/2 to 0 in the relative phase factor e iφ between | ↑ and | ↓ states.
Interestingly, the transition point and quadrupole order pattern that replaces octupole order can be obtained from symmetry considerations alone, by virtue of the duality transformations in pseudospin honeycomb models [36]. To this end, we use the explicit form of τ γ given in Eq. (5) and re-write Eq. (16) as follows: . (17) Here, λ = 2J/J τ , and the overall energy scale equal to J τ /2 is not shown. This equation has exactly the same structure as the extended Kitaev model, written in the hexagonal coordinate frame [36]. Simple re-labeling of the spin axes (x, y, z) ↔ (Y, Z, X), and a term-by-term comparison of Eq. (17) with Eq. (A1) of Ref. [36] gives the following correspondence: J XY = 1 − λ, J Z = λ, A = 1, and B = 0. (We note that B term of Ref. [36] couples in-plane and out-of-plane components of spins; for the present E g problem, finite B would imply linear quadrupole-octupole coupling which is forbidden by TR symmetry). Next, we use the relations (A2-A5) of Ref. [36] to obtain the parameters K, Γ,J, and Γ , which define the extended Kitaev model in the octahedral axes frame [38] (we useJ to avoid confusion with J in our models): So far, we have shown that Eqs. (16) and (17) correspond to the extended Kitaev model at the specific parameter set. The virtue of this mapping is that at λ = 2J/J τ = 1, we see thatJ = Γ = 0. Thus, at this point, the model is isomorphic to the K = Γ = 1 model, which in turn, is dual to the isotropic Heisenberg model, see Table I of Ref. [36]. This leads to a remarkable observation that at J τ = 2J, the highly anisotropic Hamiltonian (16) is dual to the effective FM Heisenberg model H ij = −Js i ·s j . The duality transformation involves a six-sublattice rotation matrix T 6 [36], which converts inplane FM order of effective spinss into a vortex pattern of s z and s x moments in our model. This quadrupole order is shown in Fig. 2(a) (cf. Fig. 2(e) of Ref. [36]). On the other hand, out-of-plane FM order ofs corresponds to octupolar AF order of s y moments already discussed above. Being dual to the eigenstates of a hidden FM Heisenberg model, these vortex and AF states are "fluctuation free", and low-energy excitations are magnons with a quadratic dispersion.
The exact degeneracy of these two states is lifted as soon as λ deviates from its critical value 1. From Eq. (17), we see that the corrections to the SU(2) point Hamiltonian H(λ=1) read as (1 − λ)(s z i s z j + s x i s x j − s y i s y j ). This term acts as an easy-plane or easy-axis anisotropy, selecting quadrupolar vortex order if 1 − λ > 0, and s y octupolar AF order if 1 − λ < 0. It also opens a finite gap in magnon spectra.
The symmetry-based considerations above are confirmed by numerical studies. Fig. 2(a) shows the phase diagram of Hamiltonian (17), obtained by classical theory and exact diagonalization (ED) on a C 3 symmetric 24-site cluster. The first-order transition between the vortex-quadrupole and AF-octupole phases occurs at λ = 1. Across the spin-flop transition, the inplane quadrupole moments (black arrows) flip to the out-of-plane octupole moments (red arrows). We investigated the quantum effects using linear spin wave theory (LSWT). In the quadrupole phase, the zero-point magnon energy δE depends on the vortex pattern orientation [specified by the angle ϕ in Fig. 2(a)]. The state with ϕ = π/6 has the lowest energy, see Fig. 2(b). Note that the pinning potential is extremely weak, so the spins are almost free to rotate (globally) within a quadrupolar plane. This implies the presence of lowenergy quadrupole moment fluctuations.
The calculated magnon gaps are presented in Fig. 2(c) near the spin-flop transition area (λ ∼ 1). In the planartype quadrupole phase, there are two different gaps, ∆ ⊥ and ∆ , associated with the out-of-plane and in-plane magnon modes. The out-of-plane gap is finite already within LSWT, and proportional to the deviation from the hidden FM SU(2) point: ∆ ⊥ ∝ (1 − λ). The in-plane gap, on the other hand, is zero within LSWT because the classical energy of the vortex pattern is independent of the in-plane rotation angle, ϕ. Planar anisotropy δE(ϕ) and the corresponding gap only appears beyond LSWT level, via quantum order-from-disorder mechanism, and thus is small: ∆ ∝ (1 − λ) 2 . The octupole phase (λ > 1) has uniaxial symmetry and a two-fold degenerate magnon dispersion with the gap ∆ ∝ (λ − 1). As expected, the ordered moments (dashed line) are fully saturated at the hidden FM point. Away from this point, they are reduced by quantum fluctuations (the effect is stronger in the quadrupole phase owing to the presence of a soft in-plane magnon mode).

C. Triangular lattice
The behavior of the model (16), or its equivalent (17), on a triangular lattice is of interest too. It is also relevant to double-perovskites, where the face-centered-cubic (fcc) lattice of magnetic ions can be viewed as triangular planes stacked along the [111] direction. The triangular lattice is non-bipartite and a paradigmatic example of geometrical frustration for AF Ising-type models. This is exactly the case for octupolar interactions here: see the Js y i s y j or λs y i s y j terms with positive J and λ values in Eqs. (16) and (17), respectively.
Leaving full exploration of the model for future study, we discuss now its global phase behavior on a triangular lattice. Inspection of the classical ground states, supported by classical Monte Carlo simulations suggest that there are (at least) three distinct phases, shown in Fig. 3(a), which are realized when the parameter λ = 2J/J τ varies from pure J τ limit to a dominant J regime. Stripy-quadrupole phase I at small λ is essentially the same state as found earlier in compass model [25,26]; note, however, that the spin pattern in our "anti-compass" τ model is rotated by 90 • , for the reasons discussed in Sec. III A. This state has a classical energy per site E I = −(3 − λ)s 2 (in units of J τ /2).
At large λ, the ground state is driven by the J interaction, which is of ferro-type in the quadrupolar channel, while octupole AF coupling is equally strong but frustrated. This results in simple ferro-quadrupole order (phase III), with energy E III = −3(λ − 1)s 2 . Classically, moments can freely rotate within a quadrupolar plane, but quantum effects generate in-plane anisotropy, pinning the ordered moments along the middle of the two bonds (ϕ = π 6 ). The anisotropy is rather weak: magnon zero-point energy, calculated at λ = 2, varies only by δE(ϕ) 7 × 10 −4 J τ .
In the above states I and III, the octupolar Ising interaction λs y i s y j was left "unused" because of its frustrating nature. At intermediate λ values, however, this coupling is actually larger than the quadrupolar one (as J τ and J quadrupole terms are of different sign and compete). Therefore, an intermediate state between I and III, which finds a way to resolve "triangular" frustration and activates large octupole couplings, is expected.
The pseudospin ordering pattern, whose unit cell is shown in Fig. 3(b) and further detailed in Fig. 3(c), does exactly this job. In this state, the original triangular lattice is divided into two, honeycomb and triangular sublattices. The spins on the honeycomb sublattice are canted and carry both octupole and quadrupole moments. The latter condense into a vortex pattern similar to what shown in Fig. 2(a). While the out-of-plane components form a ferro-octupole order. The honeycomb octupole moment is largely (but not fully) compensated by down-oriented octupoles residing at the middle of every hexagon. As a whole, phase II represents a coherent superposition of vortex-quadrupole and ferri-octupole orders. Such a mixture of the different rank multipoles is rather unusual. We also note that this order is noncoplanar and has a large unit cell which helps to relieve the frustrations inherent to spin-orbital models. In this sense, the case is similar to a complex behavior of spinorbit pseudospins J = 1/2 of d 5 ions on a triangular lattice [3,39,40].
As a function of spin canting angle θ, the classicial energy of the mixed state II is obtained as follows: Here, the second term originates from coupling between honeycomb lattice octupoles s y ∝ sin θ with those residing at the hexagon centers. Minimizing E II (θ) with respect to θ, we find sin θ = λ/(1 + λ). This gives a ground state energy of phase II (per site): Comparing this result with E I and E III obtained above, we find the phase transition points λ 1 and λ 2 : 1.69. (24) In between λ 1 and λ 2 , the angle θ varies from 34 • to 39 • , and the size of octupole moment on honeycomb sites 2| s y | = sin θ 0.56 − 0.63. This nearly compensates pure octupole moments from triangular sublattice, leaving rather small total octupole moment per site: 2| s y | tot = 1 3 (2 sin θ − 1) 0.04 − 0.09. Fig. 3(b) are two different vortex patterns, related to each other by in-plane rotations of spins. Classically, these states are degenerate. Quantum zeropoint energies, calculated within LSWT for these two ground states, slightly differ; the right one is lower by δE 2 × 10 −4 J τ . This result implies that in-plane magnon excitations acquire a small but finite gap, generated by the order-from-disorder mechanism. Out-ofplane excitations, corresponding to fluctuations between quadrupole and octupole sectors, are gapped out already on a classical level.

Displayed in
Overall, the E g pseudospin J τ −J model (16) on triangular lattice contains rich physics yet to be fully explored theoretically. The above results also should encourage experimental work finding and studying 5d 2 compounds with quasi-two dimensional honeycomb and triangular lattice structures.

D. Double perovskites
Now, we move to the DP lattice which motivated this study. DPs are special because the magnetic ions are widely separated from each other and thus interact weakly. This implies that the pseudospin one-half description, which assumes that the intersite interactions are less than on-site cubic splitting ∆ c , is best justified in DP compounds.
The dominant exchange channel in DPs is due to the single-orbital process considered in Sec. III A. This results in Kugel-Khomskii type Hamiltonian (10) acting in the pure quadrupole τ -channel: J τ τ iγ τ jγ . To our knowledge, the ground state of this model on fcc lattice (formed by magnetic ions in DPs) has not yet been considered. To address the behavior of the τ -model on the highly frustrated fcc lattice, we perform classical Monte-Carlo simulations on the related model, J τ n iγ n jγ , where pseudospins τ = 1/2 are replaced by classical vectors of unit length (n 2 = 1).
The simulated annealing Monte Carlo is performed for DP and, for comparison, triangular lattices. We use 1372 sites (7 × 7 × 7 unit cell) of DP and 1296 sites (36 × 36) of triangular lattice with periodic boundary conditions. Monte Carlo simulations were performed using the ALPS project library [41][42][43]. We find a collinear AF-quadrupole order at low temperatures, and the ordering pattern in DP lattice is displayed in Fig. 4(a). The moment is along the bond direction, and there are 8 anti-parallel and 4 parallel nearest-neighbors. Within the [111] planes, moments form a stripy pattern as in the triangular lattice, see phase I in Fig. 3(a). Temperature dependence of the ordered moment length n in Fig. 4(b) shows that the ordering sets in at T c 1.6J τ in DP lattice. This is about three times higher than T c /J τ in the triangular lattice, most likely due to increased dimensionality. Quantifying the quadrupole moment reduction by quantum fluctuations in fcc lattice is an interesting but challenging problem, and left for future study. It should also be noticed that T c in the actual model with quantum spin τ = 1/2 is different from the above Monte Carlo result. Roughly, an upper limit for the rescaling factor can be obtained by replacing n 2 = 1 by s(s + 1). For the present case of spin one-half, this gives an estimate of T c ∼ J τ . A rather low value of T c (despite a large coordination number 12) is presumably due to frustrations of the model on fcc lattice.
In principle, the 90 • bonding superexchange via nonmagnetic B sites is possible in DP lattice, due to the extended nature of 5d orbitals. However, this process involves many hopping steps (B i -O-B -O-B j ), and the corresponding indirect hopping t and hence J must be small. Therefore, even though the full exchange Hamiltonian in DPs is formally given by Eq. (16), the singleorbital quadrupole interaction J τ should dominate over J coupling between the octupoles. This implies that spinorbital exchange in DPs uniquely supports AF order of quadrupole moments, breaking underlying discrete point group symmetries both in real and pseudospin spaces, but preserves TR symmetry.
We should note that the above result is obtained at large cubic splitting ∆ c between the excited T 2g states and pseudospin E g doublet. The previous calculations in the other limit, i.e. neglecting cubic splitting and using full J = 2 Hilbert space instead [20] found no octupolar instability, too. Ref. [18] suggested that in the intermediate case, when the E g doublet is formed and T 2g triplet is not too high, the virtual states may generate octupolar interactions that are strong enough to overcome quadrupolar couplings. We now inspect this possibility, by considering contributions of the virtual T 2g states to the effective pseudospin E g Hamiltonian. This brings us to Jahn-Teller physics which operates not only within the ground state E g doublet, but also connects it with the T 2g triplet. We consider a linear coupling between the octahedral normal modes Q Γ and electron quadrupolar moments O Γ of symmetry Γ. For t 2g orbital systems, the Γ 3 doublet (Q 3 and Q 2 modes of 3z 2 − r 2 and x 2 − y 2 symmetries, respectively) as well as Γ 5 triplet modes (Q xy , etc) are relevant [6]. Microscopically, orbital-lattice coupling in the Γ 3 channel splits the t 2g orbital levels, while coupling to Γ 5 modes is orbital non-diagonal, e.g. Q xy distortion mixes d yz and d zx wavefunctions.
In terms of two-electron J = 2 quadrupoles, the Jahn-Teller couplings in the above two channels read as follows (summation over lattice sites i is implied): where O 3 and O 2 quadrupoles are defined in Eqs. (3) and (4), while O xy = (J x J y + J y J x )/2 √ 3, etc. Within the ground state E g doublet, Γ 3 coupling takes a form familiar from the e g orbital JT problem [2,21]: On the other hand, quadrupolar operators in H JT (Γ 5 ) have no matrix elements within the pseudospin subspace; instead, they connect the E g doublet to the excited T 2g states. This leads to a so-called second-order or "pseudo-Jahn-Teller" effect [44] which operates through the mixing of the ground and excited states. In spin-orbit coupled systems, this effect modulates the spatial shape of the pseudospin wavefunctions and generates new terms in low-energy effective Hamiltonians [45].

B. Pseudospin interactions mediated by Jahn-Teller coupling
Spatial correlations between the octahedral deformations on different sites mediate interactions between quadrupolar moments [21,46]. Typically, these interactions cooperate with the Kugel-Khomskii mechanism of orbital ordering [2]. In most TM compounds, the JT centers share common oxygens and thus stay in direct contact with each other. This leaves little room for singleion JT dynamics. In the DP lattice, however, the JT ions have no common oxygen, so they have to interact by exchanging virtual phonons. Since JT phonon modes disperse weakly, this interaction is much smaller than in perovskites. This has two important consequences: (i) cooperative JT couplings are weak enough so that the E g pseudospin description remains valid, and (ii) single-ion JT dynamics, intrinsic to non-Kramers E g states, may develop. (a) Anti-ferro correlated Q3 distortions on xy plane of the DP lattice, leading to V coupling (28) between quadrupole moments of Γ3 symmetry. (b) Anti-ferro correlated Qxy distortions on xy plane, leading to V coupling (29) in the Γ5 quadrupolar channel.
We start with pseudospin interactions mediated by JT coupling in the Γ 3 channel (27). Two ions in the xy plane couple most efficiently via Q 3 type distortions, as illustrated in Fig. 5(a). The corresponding quadrupole interaction is whose strength V = −g 2 Q 3i Q 3j ω=0 > 0 is given by non-local static susceptibility of the Q 3 modes. In DPs, antiferro-type intersite correlations Q 3i Q 3j < 0 arise due to a finite dispersion δω q of the corresponding optical phonons with energy ω 0 . This gives a rough scale of V as a small fraction ∝ δω q /ω 0 ∼ 0.1 of a single-ion JT stabilization energy E JT . For d 1 Os in DPs, Refs. [47,48] evaluated E JT ∼ 20 meV (this might be lower for d 2 ), suggesting V ∼ 1 meV, i.e. of the same order as the exchange coupling J τ for DPs. Relatively small values of E JT (and hence JT mediated coupling V ) may indicate a rather weak coupling of the diffuse 5d orbitals to the lattice. We note that considering both Q 3 and Q 2 modes would induce also s x i s x j and s z i s x j type terms; however, these are not essential in the context of possible octupole order. The main result is that the quadrupole interactions, mediated by Γ 3 type JT phonons cooperate with the exchange coupling, i.e. J τ → J τ + 4V , as in the case of the usual e g orbital problem in manganites.
The JT coupling in the Γ 5 channel (26) works differently, and it actually leads to octupolar coupling s y i s y j .
For a pair in the xy plane, we consider that Q xy type distortions shown in Fig. 5(b) are most relevant. This leads to the following interaction between Γ 5 -type quadrupole moments: with positive V = −(g ) 2 Q i xy Q j xy ω=0 . The interactions on yz and zx planes are induced by Q yz and Q zx type modes. For t 2g orbital systems, constants V and V are expected to be of the same order, but their ratio depends on material details. In d 1 Os double-perovskite, coupling to Γ 3 type modes are stronger [47] which might also be the case in d 2 osmates.
As said above, Γ 5 quadrupoles have no matrix elements within E g doublet; instead, they create transitions from the E g doublet to excited T 2g states. For example, where hard-core bosons f and T belong to pseudospin (f ↑ , f ↓ ) and triplet (T x , T y , T z ) sectors, correspondingly. Thus, the pairwise interaction V O i xy O j xy may (a) excite a pair of triplons, and (b) lead to dispersion and broadening of the T -excitons. In the present context, we are interested in the pair generation process which dynamically mixes up the ground and excited state wavefunctions, modifying thereby the pseudospin exchange Hamiltonian. Specifically, this process activates the composite operators (SL 2 z ) i and (SL 2 z ) j in the exchange Hamiltonian of Eq. (9). These operators have nondiagonal E g ↔ T 2g matrix elements [e.g., , and are thus sensitive to the admixture of JT induced triplet states in pseudospin wavefunctions. Now, we assume that the exchange J τ and quadrupole V couplings are small compared to cubic splitting ∆ c . In other words, we assume that the dispersion and broadening of triplon excitations, caused by these interactions, is smaller than ∆ c , and thus the pseudospin description is valid. This is exactly what is observed in experiment [8]. Then we proceed along the lines of Ref. [18], eliminating virtual triplon pairs perturbatively. This results in the following effective Hamiltonian (for pairs on the xy plane), which now includes the exchange as well as JT-coupling mediated interactions: For γ = x, y bonds, one has to replace s z → τ γ and s x →τ γ . Effective parameters read as: This implies that quadrupolar τ -interactions dominate over those in the octupolar s y sector. Moreover, J y /J z can be reduced by dynamical JT effects, which suppresses the interactions involving s y operators by the Ham factor [6]. The result (33) is natural for 5d ions on DP lattices, where intersite interactions J τ , V , and V between widely separated ions should be much less than the single-ion energy ∆ c . The latter is driven by large SOC for 5d electrons; we also think that the cubic splitting ∆ c is further enhanced by the dynamical JT effect stabilizing the E g doublet against the T 2g states. It would be interesting to check the latter point by quantum chemistry calculations.
The above findings suggest that a single phase transition observed in 5d 2 osmium DP oxides [8] is driven by a quadrupolar ordering of E g doublets. Comparison of the Monte-Carlo result (Sec. IV D) for the quadrupolar ordering temperature ∼ J τ + 4V (including now JT-phonon mediated coupling V ) with the experimental transition temperature 30 − 50 K gives J τ + 4V ∼ 4 meV, consistent with the above estimates of intersite couplings. Concomitant lattice distortions might be small for several reasons: first, JT coupling is already weakened by the SOC effect which partially removes orbital degeneracy; second, on-site JT vibronic dynamics [6,44] and intersite quantum fluctuations reduce the pseudospin order parameter and hence the static distortions, to the levels that are difficult to detect directly by x-ray diffraction. In general, JT coupling seems to be rather weak for 5d orbitals; indeed, quadrupole order induced distortions in 5d 1 DPs have been found to be extremely small [14] or below the resolution limit [15]. However, quadrupolar order should lead to changes in phonon spectra that should be well visible in Raman and optical data. Possible signatures of the dynamical JT effect, e.g. transitions between vibronic levels [49] are of special interest. Also, a quadrupolar order can be probed by nuclear magnetic and quadrupole resonance experiments. Further experiments, especially on single crystal samples, are necessary to identify the nature of the "hidden" order parameter in DP 5d 2 compounds.

C. Coupling to local distortions: Induced magnetic moments
While a quadrupole order of JT active E g doublets sounds natural, this picture cannot explain TR symmetry breaking observed in 5d 2 DPs [9][10][11]. One possible explanation is to attribute this effect to magnetic moments induced by defects (e.g., B ↔ B site disorder). Near defects, non-cubic crystal fields can modify the E g doublet wavefunctions, or even completely destroy the pseudospin description, thus recovering the J dipole moments at least partially. Nuclear magnetic and quadrupole resonance lineshapes may quantify such magnetic state inhomogeneities. In fact, the signatures of spin disorder and freezing are rather common in d 2 DPs [9][10][11]. To show how symmetry lowering distortions affect the physical content of the E g doublets and their interactions, we consider JT coupling of Γ 5 symmetry quadrupoles O αβ to the corresponding local distortions e αβ : The local quadrupolar fields (34) modify the E g doublet functions at site i as follows: where s = √ 3/2, c = 1/2 [the normalization factors p 2 ↑ = 1 + s 2 (δ 2 x + δ 2 y ) and p 2 ↓ = 1 + δ 2 z + c 2 (δ 2 x + δ 2 y ) are not shown]. The parameters δ z = g e xy /∆ c , etc quantify the degree of admixture of virtual triplet states T x , T y , and T z into the ground state due to strain e αβ field. This admixture "magnetizes" the E g doublet, by inducing a dipolar component into the s y operator. By calculating matrix elements of total angular momentum J within the modified E g doublet (35,36), we find an induced moment J iα = 4δ iα s y i , illustrating a partial recovery of the dipolar moments due to local distortions. The corresponding magnetic moment, which is carried by the s y operator, is M iα = 2δ iα s y i (using g factor g = 1/2 of J = 2 state). In principle, a direct link between lattice distortions and magnetism is generic to all spin-orbit Mott insulators. In 5d 2 ion systems, where the non-magnetic nature of the E g doublet is protected by cubic symmetry (i.e. independent of covalency, etc), lattice distortions have an especially strong impact on magnetic properties. To illustrate this point further, we may consider uniform strain applied along the [111] axis of a crystal: e xy = e yz = e zx = e/ √ 3. This induces a magnetic moment with g factor g = 2 √ 3 ∆tr ∆c , where ∆ tr = g e is the strain induced field, while the g factors within the [111] plane remain zero. This results in an extreme anisotropy of the magnetic response to lattice distortions in 5d 2 systems.
Non-cubic crystal fields also induce new couplings between pseudospin moments. Projection of the exchange interaction (9) onto the E g doublet with "distorted" wavefunctions (35,36) modifies the J τ term in Eqs. (10) and (16) as follows: where a new term, the s y i s y j coupling which operates in the magnetic channel, appears. Its relative strength is given by a (z) ij = 9(δ iz δ jz + 1 4 δ ix δ jx + 1 4 δ iy δ jy ) for z type bonds (results for γ = x, y follow from symmetry). We see that even rather small strain fields are sufficient to support s y order locally: the new term becomes comparable with quadrupolar τ -coupling already at δ ∼ 1/3. The sign of a ij depends on the relative orientation of local distortions; for antiferro-type distortions (δ i δ j < 0), the s y moments are coupled ferromagnetically, and vice versa (i.e. following Goodenough-Kanamori rules).
The lattice effects on the physical content of the E g doublets should be essential for understanding the magnetic properties of 5d 2 osmates. A qualitative picture is that while pseudospin one-half ordering in these compounds is predominantly of a quadrupolar type, there should also be a weak and spatially inhomogeneous dipolar component of the pseudospin order parameter, induced by the random distortions inevitable in real crystals.

VI. CONCLUSIONS
We have developed a microscopic theory for multipole orders in spin-orbit Mott insulators of non-Kramers d 2 ions, which have a non-magnetic doublet ground state of E g symmetry. The exchange Hamiltonians for various hopping processes are derived and expressed in terms of pseudospin-1/2 operators. Reflecting the spin-orbital mixed nature of the E g wavefunctions, the pseudospin interactions are in general strongly anisotropic and depend on the bond directions. The phase behavior of these models on different lattices are considered by means of analytical and numerical methods.
On a honeycomb lattice, we find that interplay between direct overlap of t 2g orbitals and their hopping via ligand ions gives rise to a competition between two distinct multipole orders: vortex-type quadrupole order and collinear AF octupole order. These two states become degenerate at the parameter point where the model has a hidden SU(2) symmetry and is isomorphic to the FM Heisenberg model. The model also can be mapped to the extended Kitaev model for d 5 systems, which is useful to understand its global phase behavior.
On a triangular lattice, a combination of geometrical and spin-orbital frustrations result in a novel type of ordering which can be viewed as a coherent superposition of vortex-type quadrupole and ferri-type octupole orders. This complex state appears as an intermediate phase between collinear AF and FM quadrupole states, as a compromise to their competition.
Double perovskite compounds of d 2 ions with strong SOC such as osmium Os 6+ or rhenium Re 5+ are discussed in more detail, including also Jahn-Teller coupling of electron quadrupole moments to lattice degrees of freedom. We find that the exchange and JT effects do work cooperatively to support quadrupole order in DP lattices. Static lattice distortions associated with this order are expected to be small, due to a reduction of the order parameters by the dynamical JT effect and pseudospin frustrations caused by bond-dependent interactions on the fcc lattice. Nevertheless, quadrupole order should lead to well detectable changes in phonon spectra as well as in nuclear magnetic resonance lineshapes. Signatures of JT dynamics in Raman and optical data should also be interesting to look for. A possible scenario for TR symmetry breaking due to noncubic crystal fields near defects is suggested. We show that the local distortions modify the ground state wavefunctions, induce local magnetic moments and enhance their exchange interactions, illustrating the importance of the lattice effects for interpretation of the experimental data.
Overall, the models introduced and discussed in this paper are of interest in their own right. From a materials perspective, our findings suggest rich physics yet to be explored in 5d 2 spin-orbit Mott insulators. The present work forms a theoretical basis for the future research of these compounds.