Interplay between multipolar spin interactions, Jahn-Teller effect and electronic correlation in a J_{eff}=\frac{3}{2} insulator
Dario Fiore Mosca, Leonid V. Pourovskii, Beom Hyun Kim, Peitao Liu, Samuele Sanna, Federico Boscherini, Sergii Khmelevskyi, Cesare Franchini
IInterplay between multipolar spin interactions, Jahn-Teller effect and electroniccorrelation in a J eff = insulator Dario Fiore Mosca, Leonid V. Pourovskii,
2, 3
Beom Hyun Kim, Peitao Liu, Samuele Sanna, Federico Boscherini, Sergii Khmelevskyi, and Cesare Franchini
1, 5, ∗ University of Vienna, Faculty of Physics and Center for Computational Materials Science, Vienna, Austria Centre de Physique Th´eorique, Ecole Polytechnique, CNRS,Institut Polytechnique de Paris, 91128 Palaiseau Cedex, France Coll`ege de France, 11 place Marcelin Berthelot, 75005 Paris, France Korea Institute for Advanced Study, Seoul 02455, South Korea Department of Physics and Astronomy, Alma Mater Studiorum - Universit`a di Bologna, Bologna, 40127 Italy Center for Computational Materials Science, Institute for Applied Physics,Vienna University of Technology, Wiedner Hauptstrasse - , Vienna, Austria (Dated: February 23, 2021)In this work we study the complex entanglement between spin interactions, electron correlationand Janh-Teller structural instabilities in the 5d J eff = spin-orbit coupled double perovskiteBa NaOsO using first principles approaches. By combining non-collinear magnetic calculationswith multipolar pseudospin Hamiltonian analysis and many-body techniques we elucidate the ori-gin of the observed quadrupolar canted antifferomagnetic. We show that the non-collinear mag-netic order originates from Jahn-Teller distortions due to the cooperation of Heisenberg exchange,quadrupolar spin-spin terms and both dipolar and multipolar Dzyaloshinskii-Moriya interactions.We find a strong competition between ferromagnetic and antiferromagnetic canted and collinearquadrupolar magnetic phases: the transition from one magnetic order to another can be controlledby the strength of the electronic correlation ( U ) and by the degree of Jahn-Teller distortions. Double perovskites with strong spin-orbit coupling(SOC) represent a unique playground for the emer-gence of exotic spin-orbital-lattice entangled phases [1–5].Structural and orbital frustrations [6], SOC-induced for-mation of effective
J manifolds [7] and relatively strongelectronic correlation give rise to a wide range of un-usual electronic and magnetic phases including so-calledDirac-Mott insulating states [8, 9], multipolar spin in-teractions [1, 2, 8] and entangled SOC-Jahn-Teller ef-fects [3, 10].A particularly interesting case is that of 5 d dou-ble perovskites characterized by an effective total an-gular momentum J eff = 3 /
2, an example of whichis Ba NaOsO (BNOO). BNOO is Dirac-Mott insula-tor [8, 11] with strong SOC and electronic correlation ef-fects [12–14]. Magnetization measurements indicate thatbelow 6.8 K BNOO exhibits weak ferromagnetic interac-tions and a small ordered spin moment of 0.2 µ B [15–18]. Refined NMR analysis confirmed the theoreticallypredicted canted spin state with two magnetic sublat-tices [2, 19]. In this spin model, depicted in Fig.1(a), thespins align ferromagnetically within the xy plane withalternating canting angle φ in adjacent planes along the[001] axis [see Fig.1(c)]. As the optimal canting angle ex-tracted from the NMR data, φ expt =67 ◦ , is closer to theAFM limit (90 ◦ ) than to the FM one (0 ◦ ) it is more ap-propriate to use the jargon canted-antiferromagnetic (c-AFM), rather than c-FM as initially proposed [2, 19, 20]. ∗ [email protected] The most stunning aspect of BNOO is that this c-AFMphase is associated with a local point symmetry break-ing manifested by local deformations of the BO octa-hedra [19, 21]. However, x-ray diffraction shows thatthe cubic symmetry is maintained even at low temper-ature [16] and no direct evidence of Jahn-Teller (JT)instabilities have been detected by NMR [19, 21], eventhough JT instabilities would be expected in a d elec-tronic configuration [6, 22] (one of the models proposedin Ref. [19] involves a contraction of the OsO octahedra,assimilable to a local Q -JT distortion, see Fig. 1). Thisapparent violation of the JT theorem [23, 24] has beenexplained by the entanglement of vibronic dynamics withSOC, suggesting the presence of a dynamical Jahn-Tellereffect with small static component [3, 10] and orbital se-lective quadrupolar charge ordering [25]. This explainsthe breaking of local point symmetry and suggests thatthe canted magnetic ground state of BNOO should bedescribed in terms of the spin-orbital-lattice entangledstates involving high-rank spin interactions.With this study we aim to clarify the link be-tween JT instabilities and spin couplings by comput-ing the variation of the multipolar interactions due tolocal structural distortions and explain how this JT-multipolar coupling affects the spin deflection and dic-tates the formation of the c-AFM ground state. We doso by employing first principles non-collinear relativis-tic density functional theory (DFT) [29–31] within thePerdew-Burke-Ernzerhof approximation [32] using theVASP code [33–35] and many-body techniques. Theinter-site interactions have been extracted from mag- a r X i v : . [ c ond - m a t . s t r- e l ] F e b FIG. 1. (Color online) Graphical representation of (a) thecanted-AFM magnetic ground state with layer-dependent (1& 2) spin ordering and (b) associated JT vibrational mode.Following Van Vleck notation [26, 27] we obtain Q = 0 .
031 ˚ A and θ = tan − ( Q /Q ) = 79 . ◦ (plane 1) and -79.35 ◦ (plane2) [28]. This structural motif is compatible with model Aproposed by Liu et al. [21]. (c) Visualization of the staggered-type JT-distortion in the two layers. (d) Schematic definitionof the canting angle φ : φ =0 corresponds to the FM limit;by convention 0 < φ < π/ π/ < φ < π/ netically constrained DFT+U+SOC calculations, fromDFT+dynamical mean-field theory (DMFT) [36, 37] inthe Hubbard-I approximation (HI) [38, 39] in conjunc-tion with the magnetic force theorem (FT) to effec-tive intersite interactions (abbreviated, FT-HI) [40], aswell as from exact diagonalization (ED) of a microscopiceffective Hamiltonians [41]. The charge self-consistentDFT+HI approach[42, 43] is implemented using theWien-2k package [44] and TRIQS library [45, 46] usingthe projective technique of Refs. [42, 47]. A more de-tailed description of the computational protocol is givenin the method section and in the SM [27].We find that in the fully undistorted (JT-quenched)limit the spins prefer to align ferromagnetically butthe onset of staggered JT distortions activates bothfirst-order (dipole-dipole) and second order (quadrupole-quadrupole + dipole-octupole) Dzyaloshinskii-Moriya(DM) interactions [48] which drive the development of the observed c-AFM phase at the optimum JT distor-tion. Moreover, our data show that the multipolar c-AFM ground state is highly sensitive to the degree of JTdistortions and to the strength of the Coulomb interac-tion U : slightly lower U or larger JT distortions shiftthe system to a competing canted-ferromagnetic (c-FM)state.We start by comparing the dependence of the energy onthe canting angle for the JT-distorted and JT-quenched(undistorted) phases in Fig.2, obtained by changing step-wise the canting angle φ from the FM ( φ = 0 ◦ ) to theAFM ( φ = 90 ◦ ) limits. The JT structure exhibits twominima, a lower c-AFM ground state found at φ = 66 ◦ ,in excellent agreement with experiment (67 ◦ [19]) and ac-FM phase just 0.05 meV/u.c. higher in energy. Con-versely, the undistorted structure shows a clear minimumfor a parallel spin alignment followed by an AFM state.The size of the spin ( m s ) and orbital ( m o ) moments areessentially the same in both structures ( m s ≈ . µ B and m o ≈ . µ B ). In the ground state JT-distorted structurethe resulting net ordered moment is 0 . µ B , in very goodagreement with measurements (0.2 µ B ) [16, 19]. Theseresults clearly indicate that the JT effect is essential tostabilize the observed canted ground state.To decipher the role of JT distortions on the spin de-flection we have extracted the intersite spin interactionsby fitting the canting energies curves with a pseudospinHamiltonian in the general multipolar form [49, 50]: H ij = (cid:88) K,K (cid:48) (cid:88)
Q,Q (cid:48) I QQ (cid:48) KK (cid:48) O KQ ( J i ) O K (cid:48) Q (cid:48) ( J j ) , (1)where i , j are the site indexes, I QQ (cid:48) KK (cid:48) the coupling con-stants and O KQ ( J i ) and O K (cid:48) Q (cid:48) ( J j ) are the multipolar ten-sor operators of rank K and K (cid:48) , with K ≤ J i , Q = − K, ..., K , K (cid:48) ≤ J j , Q (cid:48) = − K (cid:48) , ..., K (cid:48) . Themultipolar operators O are estimated within the meanfield approximation, i.e. O ( J i ) O ( J j ) ≈ (cid:104) O ( J i ) (cid:105) O ( J j ) + O ( J i ) (cid:104) O ( J j ) (cid:105) − (cid:104) O ( J i ) (cid:105)(cid:104) O ( J j ) (cid:105) . After applying symme-try properties of pseudovectors and multipolar operatorsfor the J eff = (or Γ ) state [1, 51] as well as crys-tal symmetry we obtained the following fitting formulawhich expresses the magnetic energy as a function of thecanting angle φ (the derivation is given in the SM [27]): E ( θ ) = A cos(2 φ ) + A cos(4 φ ) + A cos(6 φ )+ B sin(2 φ ) + B sin(4 φ ) (2)To validate the DFT-mapping and decode the specifictype of spin couplings at play we have performed a directcalculation of the intersite interactions using the FT-HImethod proposed by Pourovskii [40] and ED [41], tak-ing the undistorted phase as a reference (FT-HI and EDcannot be applied to the JT-distorted phase). All possi-ble on-site fluctuations within the paramagnetic ground FIG. 2. (Color online) Dependence of the energy on the canting angle and spin-spin interaction matrix. (a) CalculatedDFT+U+SOC φ dependent energies for the optimized JT distorted phase (circles) and corresponding fitting curve obtainedusing Eq. 2. The fitting curves obtained by removing the A (o-o), the A (q-q) and A -A -B (q-q- and o-o) coefficients arealso displayed; (b) Comparison between DFT, FT-HI and ED canting energies for the undistorted (Q =0) phase. (c) FT-HIinteraction matrix showing the non-zero and non-constant mean values of the tensor operators (cid:104) O (cid:105) in the pseudo spin eigenbasis(see Eq. 1) summed over the out-of plane interacting bonds (see SM for details [27]). To guarantee a quantitatively consistentcomparison between DFT and FT-HI the DFT FM-AFM energy difference is aligned to that obtained by FT-HI (see SM fordetails [27]). state of a shell, encoded by the corresponding multipolemoments O KQ ( J i ), are included; all relevant inter-site cou-plings can thus be extracted. In the present case of the J eff = shell of Os, those are all inter-site interac-tions between dipole, quadrupole, and octupole momentsthat are permitted by the symmetry. With the effectiveHamiltonian thus obtained we evaluated the dependenceof zero-temperature total energy versus φ ; the spin-spincoupling parameters A i are subsequently extracted byfitting it to the form (2).The fitting curves follow very well the first principlesdata for both the distorted and undistorted phases asshown in Fig. 2(a,b). The optimal fitting parameters andthe corresponding FT-HI and ED results are collected inTab. I where each coefficient is associated with the spe-cific type of dipolar or multipolar couplings. Both phases Spin-spin InteractionsUndistorted JT-distortedcFM cAFMCoefficients FT-HI ED DFT DFTA -0.72 -1.06 -0.66 0.13d-d 0.84 0.39q-q -0.03 0.00d-o -1.32 -1.25o-o -0.21 -0.20A q-q -0.40 -0.48 -0.69 -0.38A o-o 0.05 0.08 -0.03 -0.08B DM d-d 0 0 0 -0.99B DM q-q, d-o 0 0 0 -0.08TABLE I. Collection of spin-spin interactions (in meV) asderived from DFT, FT-HI and ED approaches for the undis-torted Q =0 FM phase and corresponding DFT values for thefully distorted c-AFM phase. The coefficients are decomposedover their distinct d-d, d-o, o-o, and q-q character. exhibit relatively large quadrupole-quadurpole (q-q) cou-plings and non negligible octupole-octupole (o-o) inter-actions. The inclusion of o-o terms in the JT-phase isessential to shift the minimum from 70 ◦ to 66 ◦ , whereasthe formation of the two competing canted minima ismostly governed by the q-q term A . We also note thatDM interactions, embedded in coefficients B and B , areinoperative in the undistorted phase (not allowed by sym-metry), but the breaking of the local symmetry causedby the JT effect activates the DM channel which is themain responsible for stabilizing of the canted phases. Ourdata indicate that, besides the standard first order d-dDM interaction, an additional second order DM couplingemerges involving both q-q and d-o terms. This typeof interaction can be classified as multipolar DM inter-action, namely an antisymmetric DM exchange amongmultipolar tensors.After clarifying the essential features of the c-AFMground state we move now to the analysis of the compe-tition between the c-AFM and c-FM phases and how thiscompetition is controlled by the strength of U and the de-gree of JT mode Q . Being almost degenerate in energy,it is expected that small perturbations could change therelative stability of these two phases, providing insightson the mechanism favoring one or the other spin ordering.To explore this possibility we have computed the cantingenergies for different values of the relative JT distortion (cid:102) Q (from the fully distorted structure, here identified as (cid:102) Q = 1 to the undistorted (cid:102) Q =0 limit) and U (from 3.4to 2.4 eV). The results, displayed in Fig. 3, indicate thatboth degrees of freedom act on the spin deflection andcan trigger a magnetic transition from the c-AFM to c-FM.By decreasing U the c-AFM phase is progressivelydestabilized in favor of the c-FM phase which becomes FIG. 3. (Color online) Magnetic canting energy as a function of (a) U (relative global minimum set to zero energy) and (b) (cid:102) Q ( (cid:102) Q = 1 refers to the fully optimized JT phase). (c) Global magnetic phase diagram showing the dependence of the groundstate magnetic ordering (canting angle) for different values of U and (cid:102) Q . Decreasing (cid:102) Q drives the system to a collinear FMstate whereas a reduction of U washes-out the c-AFM system and stabilizes a c-FM configuration. the global minimum. A further decrease of U washes-out completely the AFM canting. Similarly, reducing (cid:102) Q acts on the canting angle by inverting the stability ofthe two canted minima and drives the systems towardsthe collinear limits as (cid:102) Q approaches 0. In order to ex-plore the combined action of JT effect and electronic cor-relation we have repeated similar calculations and con-structed a complete magnetic phase diagram as a func-tion of (cid:102) Q and U , displayed in Fig. 3(c). The phase dia-gram shows three different phases characterized by spe-cific value of the canting angle: (i) a FM phase which isfavored in the limit of vanishing JT-distortion and small U ; (ii) a c-FM phase emerges beyond a small critical JTdistortion in the whole considered U -regime and finally(iii) the c-AFM phase, corresponding to the experimen-tally observed ground state, which covers the top portionof the phase diagram and is favored by increasing the JTmode and U . For (cid:102) Q =1 the c-AFM phase represents theground state for U larger than 3.1 eV.The interplay between (cid:102) Q and U affects the exchangecouplings. The transition from one phase to the other canbe rationalized by inspecting the variation of the spin-spin interactions upon U and (cid:102) Q shown in Fig. 4. Firstwe note that the tendency to move toward a FM-typecanting is mostly triggered by the steep decrease of thed-d interaction embedded in the exchange coefficient A ,which changes sign rather quickly when U or (cid:102) Q are re-duced, thereby driving the AFM-to-FM transition. Wenote that coefficients A and B besides being relativelysmall remain essentially unaffected by tuning U at theoptimized JT distorted structure, implying that the o-ocharacter is mostly linked to the structural instabilities.This is confirmed by the observation that A is rapidly di-minished with decreasing (cid:102) Q thus revealing a direct linkbetween o-o interaction and structural distortion. TheJT distortion is therefore a determinant factor for theonset of the dipolar and multipolar DM interactions ex-pressed by the coefficient B and B , respectively. Sim- ilarly to the multipolar DM component also the dipolarone ( B ), which is zero by symmetry in the undistortedlimit, rapidly increases vs. (cid:102) Q and reaches a large valuein the fully JT-distorted phase ( B = − .
99 meV). Fi-nally, from the trend of A we infer that q-q interactionsbecome stronger with decreasing (cid:102) Q and are responsi-ble for the preservation of the two minima in the undis-torted phase. Conversely, reducing U leads to a linearreduction of A , which contributes to the disruption ofthe large- φ minimum. The above analysis can give in-sights on the role of dynamical JT effects in BNOO, thatshould be manifested by two types of magnetic fluctua-tions: FM fluctuations, controlled by the changes in A ,and canted-FM to canted-AFM fluctuations governed bythe variation of B . FIG. 4. (Color online) Evolution of the DFT+U spin-spinexchange coefficients as a function of (left) U and (right) (cid:102) Q . In conclusion, by combining many-body effectiveHamiltonians with material specific DFT-based descrip-tions we have decoded the complex entanglement be-tween structural frustration and electronic correlation inthe spin-orbit coupled and spin-canted double perovskiteBNOO. We show that this exotic magnetic order orig-inates from the anomalous Dzyaloshinskii-Morya inter-actions generated by JT -distortions, and elucidate thedirect effect of U and JT on the spin-spin interactions.This work contributes to the raising interest on atypicalmultipolar orderings in relativistic oxides [52, 53] with atransparent analysis that helps to light up the intricatephysics at play in this class of materials.
METHODS
To study the cross coupling between spin exchange in-teraction, Jahn-Teller distortions and electronic correla-tion we have integrated three different schemes: mag-netically constrained non-collinear DFT+U, DMFT inthe Hubbard-I approximation (FT-HI) combined withthe magnetic force theorem and exact diagonalization.Since the direct computation of spin exchanges at FT-HIand ED for JT-distorted structures is at present not pos-sible, we have first verified the validity of the DFT+Udata for the undistorted phase by comparing the spin ex-changes extracted from DFT+U total energies with thosecalculated explicitly at FT-HI and ED level (Fig.2(b)).In DFT+U the spin interactions are estimated by fittingthe total energies curves obtained at different canting an-gle (via magnetically constrained calculations) with thespin Hamiltonian given in Eq. 2 (derived in the SM [27].After normalizing the DFT+U data to the FT-HI/EDvalues (see SM [27] for further details) we have then ad-dressed the JT-induced modifications at DFT+U level,by performing a series of calculation stepwise from theideal undistorted phase to the fully distorted real struc-ture (Fig.3(a-b)). In order to account for the role ofelectron-electron correlation and construct the phase di-agram displayed in Fig. 3(c), these calculations have beenconducted for different values of U , from 2.4 eV to 3.4 eV.In studying the JT distortion we have followed the VanVleck notation for the definition of the main JT distor-tions Q , Q and Q . Since the value of Q is signif-icantly larger than the corresponding Q and Q value(see SM [27]) we have considered the variation of thethe spin interactions as a function of the JT parameter˜ Q , defined as the ratio between Q for the interpolatedstructure and the corresponding value for the optimalone. Additional details can be found in the SM [27],which includes Refs. [54–61]. ACKNOWLEDGEMENTS
We gratefully acknowledge the helpful discussions withV. Mitrovic. C.F. thanks S. Streltsov for useful insightson the Jahn-Teller effect. DFM acknowledges supportfrom the Vienna Doctoral School in Physics and from theUniversity of Bologna (Thesis Abroad Scholarship). LVP acknowledges the support by the European ResearchCouncil grants ERC-319286-”QMAC” and is grateful tothe computer team at CPHT for support. B.H.K wassupported by a KIAS Individual Grant (CG068701) atKorea Institute for Advanced Study. The computationalresults were achieved by using the Vienna Scientific Clus-ter (VSC). [1] G. Chen, R. Pereira, and L. Balents, Phys. Rev. B ,174440 (2010).[2] H. Ishizuka and L. Balents, Phys. Rev. B , 184422(2014).[3] N. Iwahara, V. Vieru, and L. F. Chibotaru, Phys. Rev.B , 075138 (2018).[4] J. Romh´anyi, L. Balents, and G. Jackeli, Phys. Rev.Lett. , 217202 (2017).[5] C. Martins, M. Aichhorn, and S. Biermann, Journal ofPhysics: Condensed Matter , 263001 (2017).[6] D. I. Khomskii and M. V. Mostovoy, Journal of PhysicsA: Mathematical and General , 9197 (2003).[7] J. B. Goodenough, Phys. Rev. , 466 (1968).[8] K.-W. Lee and W. E. Pickett, Europhysics Letters (EPL) , 37008 (2007).[9] H. J. Xiang and M.-H. Whangbo, Phys. Rev. B ,052407 (2007).[10] L. Xu, N. A. Bogdanov, A. Princep, P. Fulde, J. van denBrink, and L. Hozoi, npj Quantum Materials , 16029(2016).[11] H. L. Feng, S. Calder, M. P. Ghimire, Y.-H. Yuan, Y. Shi-rako, Y. Tsujimoto, Y. Matsushita, Z. Hu, C.-Y. Kuo,L. H. Tjeng, T.-W. Pi, Y.-L. Soo, J. He, M. Tanaka,Y. Katsuya, M. Richter, and K. Yamaura, Phys. Rev. B , 235158 (2016).[12] S. Gangopadhyay and W. E. Pickett, Phys. Rev. B ,045133 (2015).[13] S. Gangopadhyay and W. E. Pickett, Phys. Rev. B ,155126 (2016).[14] C. Wang, Y. Xu, W. Hao, R. Wang, X. Zhang, K. Sun,and X. Hao, Phys. Rev. B , 035126 (2019).[15] K. E. Stitzer, M. D. Smith, and H.-C. zur Loye, SolidState Sciences , 311 (2002).[16] A. S. Erickson, S. Misra, G. J. Miller, R. R. Gupta,Z. Schlesinger, W. A. Harrison, J. M. Kim, and I. R.Fisher, Phys. Rev. Lett. , 016404 (2007).[17] A. J. Steele, P. J. Baker, T. Lancaster, F. L. Pratt,I. Franke, S. Ghannadzadeh, P. A. Goddard, W. Hayes,D. Prabhakaran, and S. J. Blundell, Phys. Rev. B ,144416 (2011).[18] K.-H. Ahn, K. Pajskr, K.-W. Lee, and J. Kuneˇs, Phys.Rev. B , 064416 (2017).[19] L. Lu, M. Song, W. Liu, A. P. Reyes, P. Kuhns, H. O. Lee,I. R. Fisher, and V. F. Mitrovic, Nature Communications , 14407 (2017).[20] K. Willa, R. Willa, U. Welp, I. R. Fisher, A. Rydh, W.-K.Kwok, and Z. Islam, Phys. Rev. B , 041108 (2019).[21] W. Liu, R. Cong, A. P. Reyes, I. R. Fisher, and V. F.Mitrovi´c, Phys. Rev. B , 224103 (2018).[22] D. I. Khomskii, Transition Metal Compounds (Cam-bridge University Press, 2014). [23] H. A. Jahn and W. H. Bragg, Proceedings of the RoyalSociety of London. Series A - Mathematical and PhysicalSciences , 117 (1938).[24] G. A. Gehring and K. A. Gehring, Reports on Progressin Physics , 1 (1975).[25] R. Cong, R. Nanguneri, B. Rubenstein, and V. F.Mitrovi´c, Phys. Rev. B , 245141 (2019).[26] J. H. Van Vleck, The Journal of Chemical Physics , 72(1939), https://doi.org/10.1063/1.1750327.[27] Supplementary Material.[28] J. Kanamori, Journal of Applied Physics , S14 (1960),https://doi.org/10.1063/1.1984590.[29] D. Hobbs, G. Kresse, and J. Hafner, Phys. Rev. B ,11556 (2000).[30] P. Liu, S. Khmelevskyi, B. Kim, M. Marsman, D. Li,X.-Q. Chen, D. D. Sarma, G. Kresse, and C. Franchini,Phys. Rev. B , 054428 (2015).[31] S. L. Dudarev, P. Liu, D. A. Andersson, C. R. Stanek,T. Ozaki, and C. Franchini, Phys. Rev. Materials ,083802 (2019).[32] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[33] G. Kresse and J. Hafner, Phys. Rev. B , 558 (1993).[34] G. Kresse and J. Furthm¨uller, Phys. Rev. B , 11169(1996).[35] P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994).[36] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozen-berg, Rev. Mod. Phys. , 13 (1996).[37] V. I. Anisimov, A. I. Poteryaev, M. A. Korotin, A. O.Anokhin, and G. Kotliar, Journal of Physics: CondensedMatter , 7359 (1997).[38] J. Hubbard, Proc. Roy. Soc. (London) A 276 , 238 (1963).[39] A. I. Lichtenstein and M. I. Katsnelson, Phys. Rev. B ,6884 (1998).[40] L. V. Pourovskii, Phys. Rev. B , 115117 (2016).[41] B. H. Kim, G. Khaliullin, and B. I. Min, Phys. Rev.Lett. , 167205 (2012).[42] M. Aichhorn, L. Pourovskii, V. Vildosola, M. Ferrero,O. Parcollet, T. Miyake, A. Georges, and S. Biermann,Phys. Rev. B , 085101 (2009).[43] M. Aichhorn, L. Pourovskii, and A. Georges, Phys. Rev.B , 054529 (2011).[44] P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka, J. Luitz,R. Laskowski, F. Tran, and L. D. Marks, WIEN2k,An augmented Plane Wave + Local Orbitals Programfor Calculating Crystal Properties (Karlheinz Schwarz,Techn. Universit¨at Wien, Austria,ISBN 3-9501031-1-2,2018).[45] O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann,I. Krivenko, L. Messio, and P. Seth, Computer PhysicsCommunications , 398 (2015).[46] M. Aichhorn, L. Pourovskii, P. Seth, V. Vildosola,M. Zingl, O. E. Peil, X. Deng, J. Mravlje, G. J.Kraberger, C. Martins, et al. , Computer Physics Com-munications , 200 (2016).[47] B. Amadon, F. Lechermann, A. Georges, F. Jollet, T. O.Wehling, and A. I. Lichtenstein, Phys. Rev. B , 205112(2008).[48] M. Hosoi, T. Mizoguchi, T. Hinokihara, H. Hiroyasu Mat-suura, and M. Ogata, arXiv:1804.04874 [cond-mat.str-el].[49] A. Arima, M. Harvey, and K. Shimizu, Physics LettersB , 517 (1969). [50] P. Santini, S. Carretta, G. Amoretti, R. Caciuffo,N. Magnani, and G. H. Lander, Rev. Mod. Phys. ,807 (2009).[51] R. Shiina, O. Sakai, H. Shiba, and P. Thalmeier, Jour-nal of the Physical Society of Japan , 941 (1998),https://doi.org/10.1143/JPSJ.67.941.[52] D. D. Maharaj, G. Sala, M. B. Stone, E. Kermarrec,C. Ritter, F. Fauth, C. A. Marjerrison, J. E. Greedan,A. Paramekanti, and B. D. Gaulin, Phys. Rev. Lett. , 087206 (2020).[53] A. Paramekanti, D. D. Maharaj, and B. D. Gaulin, Phys.Rev. B , 054439 (2020).[54] S. M. Winter, Y. Li, H. O. Jeschke, and R. Valent´ı, Phys.Rev. B , 214431 (2016).[55] B. H. Kim, D. V. Efremov, and J. van den Brink, Phys.Rev. Materials , 014414 (2019).[56] Y. T. T. Inui and Y. Onodera, Group theory and its ap-plications in physics (Springer-Verlag Berlin Heidelberg,1990).[57] A. Horvat, L. Pourovskii, M. Aichhorn, and J. Mravlje,Phys. Rev. B , 205115 (2017).[58] C. Franchini, T. Archer, J. He, X.-Q. Chen, A. Filippetti,and S. Sanvito, Phys. Rev. B , 220402 (2011).[59] T. Archer, C. D. Pemmaraju, S. Sanvito, C. Franchini,J. He, A. Filippetti, P. Delugas, D. Puggioni, V. Fioren-tini, R. Tiwari, and P. Majumdar, Phys. Rev. B ,115114 (2011).[60] L. V. Pourovskii, private communication.[61] L. V. Pourovskii and S. Khmelevskyi, Phys. Rev. B ,094439 (2019). SUPPLEMENTARY MATERIALMETHODS
In this section the computational approaches used toinvestigate the multipolar interactions are described.
Computational details
Density Functional Theory (DFT) calculations weredone within the projector augmented wave method [35] asimplemented in the Vienna ab initio Simulation Package(VASP) [33, 34]. The Perdew-Burke-Ernzerhof approx-imation [32] for the exchange-correlation functional wasused, and, in order to take into account the relativisticeffects, we considered the fully relativistic scheme that in-cludes Spin Orbit Coupling (SOC), with a SOC energy of482 meV (computed by mapping the DFT energy differ-ence calculated with and without SOC onto a relativisticatomic Hamiltonian for d orbitals [31]). Furthermore, weadded an on-site Hubbard U correction as implementedin the Dudarev’s approach for noncollinear spin configu-rations [31]: E DF T + U = E DF T + U (cid:88) σ (cid:104)(cid:16) (cid:88) m n σm ,m (cid:17) − (cid:16) (cid:88) m ,m ˆ n σm ,m ˆ n σm ,m (cid:17)(cid:105) , (3)where ˆ n σm ,m is density matrix with m and m runningover the orbital angular momentum and σ over the spinangular momentum and U is the effective Hubbard Uinteraction. We used a U of 3.4 eV, in line with previousexperimental estimations [16]. As a reference we showin Fig. 8 the optical spectra obtained with this value of U at DFT+U+SOC and ED level.The supercell used is shown in Figure 5, where the lat-tice parameters are a/ √ × a/ √ × a , with a = 8.287 ˚Abeing the experimental room-temperature value [15]. Itcontains two formula units, or, equivalently, two osmiumatoms on parallel planes along the z direction. The recip-rocal space was sampled with a k -point mesh of 6 × × ∼ − eV.The local structure was studied by means of selectivedynamics routine, fixing the positions of the Na, Os andBa sites and allowing the oxygens to be relaxed to theirground state through the direct inversion in the iterativesubspace (RMM-DIIS) algorithm. The magnetic config-uration in the canted phase was set with noncollinearmagnetic moments constrained as in Figure 1 (b) (maintext), with φ = 67 o . Both DFT, DFT + U and DFT+ U + SOC were used in order to probe the role of thedifferent degrees of freedom. We verified that the moststable local structure is the one defined as Model A inRef. [21]. The canted procedure used for studying the exchangeinteractions – keeping the structure fixed and changingthe direction of the local magnetic moments – allows tomake use of the magnetic force theorem. To do so, wetilted the angles between magnetic moments in parallelplanes in steps of 2.5 degree starting from a AFM [¯110]and ending in the FM [110] configuration (see Figure 6)This can be achieved in VASP by means of the con-strained magnetic moment approach which constrainsthe direction of the local spin magnetic moment byadding an extra penalty energy, which has the followingexpression [30] E ( λ ) = (cid:88) I λ (cid:104) ˜M I − ˆM I (cid:0) ˆM I · M I (cid:1)(cid:105) . (4)Here the sum runs over the ionic sites I , ˆM I is the unitvector that gives the constrained direction of the mag-netic moment at site I , and (cid:126)M I is the integrated magneticmoment inside the Muffin-tin sphere Ω I . The constant λ gives the strength of the penalty energy and was chosenmanually in order to get a value of E ( λ ) (cid:28) E ( DF T ),a necessary requirement for having meaningful DFT to-tal energies. Our values of E ( λ ) were always lower than10 − eV.The result is a energy variation as a function of thecanting angle, that, as we will show in the next section,can be mapped onto a Pseudo Spin Hamiltonian (PSH)with multipolar interactions. Calculations of exchange interactions withinDFT+HI
Our charge self-consistent DFT+DMFT calculationsusing the Hubbard-I (HI) approximation for Os 5 d , ab-briviated as DFT+HI, were carried out for the paramag-netic phase and the undistored cubic structure of BNOOwith a = 8 .
287 ˚A and the Os-O distance of 1.869 ˚A.We employed the Wien-2k electronic structure package inconjunction with ”TRIQS” library implementations forthe DMFT cycle and Hubbard-I method. The spin-orbitcoupling was included in Wien2k within the standardsecond-variation treatment. The Brillouin zone (BZ) in-tegration was carried out using 400 k -points in the fullBZ.The Wannier orbitals representing t g states of Os5 d were constructed by the projective technique ofRefs. [42, 47] using the bands enclosed by the energywindow [ − .
09 : 2 .
04] eV; this window thus encloses all t g -like bands. The Coulomb interaction was specifiedby the Slater parameter F =3.2 eV and the Hund’s rulecoupling J H =0.5 eV. We employed the rotationally in-variant Kanamori Coulomb repulsion between Os t g ,which was correspondingly defined by the intra-orbitalrepulsion U K = F + 4 J H / FIG. 5. Magnetic unit cell of BNOO.FIG. 6. Canting angles for BNOO rule J K = 3 J H /
5. The double counting was evaluatedusing the fully-localized-limit expression and the atomicoccupancy N = 1 of Os 5 d shell. The DFT+HI calcula-tions were converged to 0.01 mRy in the total energy.The DFT+HI calculations for the paramagneticBNOO phase predict the ground state quadruplet J eff =3 / d shell, as expected; the excited doublet J eff = 1 / λ for the t g subshell equal to 0.3 eV. We subsequentlyemployed the FT-HI method of Ref. [40] to evaluateall nearest-neighbour exchange interactions between the J eff = 3 / V RR (cid:48) between sites R and R (cid:48) read: (cid:104) M M | V RR (cid:48) | M M (cid:105) = Tr (cid:34) G RR (cid:48) δ Σ at R (cid:48) δρ M M R (cid:48) G R (cid:48) R δ Σ at R δρ M M R (cid:35) , (5)where M is the projection quantum number, M = − / ... /
2, for the multiplet J eff = 3 / ρ M i M j R is thecorresponding element of the J eff = 3 / R , δ Σ at R δρ MiMj R is the derivative of atomic (Hubbard-I) self-energy Σ at R over a fluctuation of the ρ M i M j R ele-ment, G RR (cid:48) is the inter-site Green’s function evaluatedwithin the DFT+HI. Once all matrix elements (5) arecalculated, they are transformed to the couplings I QQ (cid:48) KK (cid:48) between on-site moments (eq. 1 of the main text) as fol-lows: I QQ (cid:48) KK (cid:48) ( RR (cid:48) ) = (cid:88) M M M M (cid:104) M M | V RR (cid:48) | M M (cid:105) (cid:2) O KQ ( J ) (cid:3) M M (cid:104) O K (cid:48) Q (cid:48) ( J ) (cid:105) M M , where (cid:2) O KQ ( J ) (cid:3) M M is the M M matrix element of thecorresponding multipolar tensor.The resulting 15 ×
15 interaction matrix I QQ (cid:48) KK (cid:48) forthe [1/2,1/2,0] Os-Os bond (omitting the irrelevantmonopole term) is shown in Table III. Interaction matri-ces for other NN bonds can be obtained from it using thecorresponding symmetry operations; the interaction ma-trix summed over all interlayer bonds and omitting theterms irrelevant for the planar magnetic order of BNOOis graphically depicted in Fig. 2c of the main text. Exact Diagonalization
To calculate the multipolar interactions between neigh-boring Os ions, we considered two-site t g -orbital sys- tems with following Hamiltonian H = (cid:88) iµνσ (cid:15) iµν c † iµσ c iνσ + λ (cid:88) iµνσσ (cid:48) ( l · s ) µσ,νσ (cid:48) c † iµσ c iνσ (cid:48) + 12 (cid:88) iσσ (cid:48) µν U µν c † iµσ c † iνσ (cid:48) c iνσ (cid:48) c iµσ + 12 (cid:88) iσσ (cid:48) µ (cid:54) = ν J µν c † iµσ c † iνσ (cid:48) c iµσ (cid:48) c iνσ + 12 (cid:88) σµ (cid:54) = ν J (cid:48) µν c † iµσ c † iµ, − σ c iν, − σ c iνσ + (cid:88) µνσ (cid:16) t µν c † µσ c νσ + h.c. (cid:17) , (6)where c † iµσ is the creation operator of t g electrons with µ ( µ = xy, yz, zx ) orbital and σ ( σ = + , − ) spin statesat the i -the site ( i = 1 , t g orbitals due to the crystalfield and spin-orbit coupling, respectively. The spin-orbitcoupling strength λ is set to be 0 . U µµ = U , U µ (cid:54) = ν = U − J H , and J µν = J (cid:48) µν = J H . We set U = 3 . J H = 0 . t g orbitals at1 and 2 sites. The hopping integral t µν and local crystalsplitting (cid:15) iµν are estimated with the Wannier projectionof the DFT calculation. We took into account all possi-ble states of two electrons in two-site cluster. Employingthe exact diagonalization method, we solved Eq. 6 andgot eigenvalues and eigenstates.Because U is much larger than strength of hoppingintegrals, lowet 16 eigenstates are dominantly deter-mined by unperturbative eigenstates | M M (cid:105) ’s ( − / ≤ M , M ≤ / E n and | Ψ n (cid:105) be n -th eigenvalueand eigenstate, respectively. The spin-orbital dynamicsof J eff = 3 / H eff = (cid:88) M M M (cid:48) M (cid:48) h M M M (cid:48) M (cid:48) | M M (cid:105)(cid:104) M (cid:48) M (cid:48) | . (7)The overlap matrix S is defined as S nm = (cid:80) M M (cid:104) Ψ n | M M (cid:105)(cid:104) M M | Ψ m (cid:105) . Coefficient h M M M (cid:48) M (cid:48) isapproximately estimated as h M M M (cid:48) M (cid:48) = (cid:88) n =1 16 (cid:88) m =1 16 (cid:88) m (cid:48) =1 E n [ S − / ] mn [ S − / ] nm (cid:48) (cid:104) M M | Ψ m (cid:105)(cid:104) Ψ m (cid:48) | M (cid:48) M (cid:48) (cid:105) , (8)where S − / is the inverse of square-root matrix of S [54, 55].Let T KQ ( J ) be the irreducible tensor operator in a J manifold ( K = 0 , · · · , J and Q = − K, · · · , K ). Accord-ing to the Wigner-Eckart theorem, T KQ ( J ) is given as T KQ ( J ) = (cid:88) MM (cid:48) (cid:104) J || T K ( J ) || J (cid:105) C JMJ (cid:48) M (cid:48) KQ √ J + 1 | J, M (cid:105) (cid:104)
J, M (cid:48) | , (9)where, (cid:104) J || T K ( J ) || J (cid:105) is the reduced matrix elementof rank K tensors and C JMJ (cid:48) M (cid:48) KQ is the Clebsch-Gordan (CG) coefficients defined as C JMJ (cid:48) M (cid:48) KQ =0 TABLE II. The interaction matrix I QQ (cid:48) KK (cid:48) calculated by DFT+HI for the nearest-neighbour Os-Os bond [1/2,1/2,0], in meV.The order of K and Q is indicated in the first two columns.-1 -0.96 -0.03 -0.74 0 0 0 0 0 0.14 -0.01 -1.05 0 0.19 0 -0.541 0 -0.03 1.36 -0.02 0 0 0 0 0 0 -0.83 -0.01 0.44 0 0.01 01 -0.74 -0.02 -0.85 0 0 0 0 0 0.52 0 0.17 0 -1.05 0 -0.09-2 0 0 0 -1.04 0 -1.02 -0.02 -0.01 0 0 0 0 0 0 0-1 0 0 0 0 -1.65 0.01 0.10 -0.01 0 0 0 0 0 0 02 0 0 0 0 -1.02 0.01 1.43 0.01 0.07 0 0 0 0 0 0 01 0 0 0 -0.02 0.10 0.01 -1.76 0.02 0 0 0 0 0 0 02 0 0 0 -0.01 -0.01 0.07 0.02 -1.42 0 0 0 0 0 0 0-3 0.14 0 0.52 0 0 0 0 0 1.69 -0.01 -0.20 0.01 0.49 -0.01 0-2 -0.01 -0.83 0 0 0 0 0 0 -0.01 -1.14 0 0.32 0.02 0 -0.02-1 -1.05 -0.01 0.17 0 0 0 0 0 -0.20 0 0.19 -0.01 0.58 0.01 -0.493 0 0 0.44 0 0 0 0 0 0 0.01 0.32 -0.01 1.34 0 -0.09 01 0.19 0 -1.05 0 0 0 0 0 0.49 0.02 0.58 0 0.32 0 0.082 0 0.01 0 0 0 0 0 0 -0.01 0 0.01 -0.09 0 -1.49 0.013 -0.54 0 -0.09 0 0 0 0 0 0 -0.02 -0.49 0 0.08 0.01 1.68 (cid:104) JM | J (cid:48) M (cid:48) KQ (cid:105) [56]. Note that CG coefficients are sat-isfied with following orthogonality relation (cid:88) KQ K + 12 J + 1 C JM JM KQ C JM (cid:48) JM (cid:48) KQ = δ M M δ M (cid:48) M (cid:48) . (10)We can easily show the following relation | J, M (cid:105) (cid:104)
J, M (cid:48) | = (cid:88) KQ K + 1 √ J + 1 C JMJM (cid:48) KQ (cid:104) J || T K ( J ) || J (cid:105) T KQ ( J ) = (cid:88) KQ A JMM (cid:48) KQ T KQ ( J ) , (11)where A JMM (cid:48) KQ = K +1 √ J +1 C JMJM (cid:48) KQ (cid:104) J || T K || J (cid:105) . The multipolar tensoroperator O KQ ( J ) of real-valued (tesseral) harmonics arebuilt as following O K ( J ) = T K ( J ) (12a) O KQ ( J ) = 1 √ (cid:16) T K − Q ( J ) + ( − Q T KQ ( J ) (cid:17) ( Q > O K − Q ( J ) = i √ (cid:16) T K − Q ( J ) − ( − Q T KQ ( J ) (cid:17) ( Q > . (12c)The projection operator | J, M (cid:105) (cid:104)
J, M (cid:48) | is given as | J, M (cid:105) (cid:104)
J, M (cid:48) | = (cid:88) KQ B JMM (cid:48) KQ O KQ ( J ) , (13)where B JMM (cid:48) K = A JMM (cid:48) K , B JMM (cid:48) KQ = A JMM (cid:48) K, − Q +( − Q A JMM (cid:48) KQ √ , and B JMM (cid:48) K, − Q = A JMM (cid:48) K, − Q − ( − Q A JMM (cid:48) KQ i √ ( Q > O KQ operators as following H eff = I I + (cid:88) (cid:48) KQ I K Q O KQ ( J )+ (cid:88) (cid:48) KQ I K Q O KQ ( J )+ (cid:88) (cid:48) K Q (cid:88) (cid:48) K Q I K K Q Q O K Q ( J ) O K Q ( J ) , (14)where I K K Q Q = (cid:80) M M M (cid:48) M (cid:48) h M M M (cid:48) M (cid:48) B JM M (cid:48) K Q B JM M (cid:48) K Q and (cid:88) (cid:48) KQ refers to the summation over all K and Q values except for K = Q = 0. OPTIMAL STRUCTURAL DATA
The energy minimum is found for an expansion of theOs-O octahedra along the x, y, and z directions withdifferent magnitude along the different crystallographicaxis. Between parallel planes along z, the distortionsare staggered in the xy plane, meaning that the Os(1)-Obond length along x (y) is equal to the Os(2)-O bondlength along y (x). Furthermore, the angles betweenbonds in the xy plane are slightly changed from the per-fect 90 degrees, forming alternated angles of ∼ o and ∼ o .The corresponding Point Group is D h − or Pmnnand allows for anti-symmetric exchange interactions. Thestructural CIF file for this compound is in Appendix.1 TABLE III. The interaction matrix I K K Q Q calculated by ED for the nearest-neighbour Os-Os bond [1/2,1/2,0], in meV. Theorder of K and Q is indicated in the first two columns.-1 -0.68 0 -0.61 0 0 0 0 0 0.06 0 -1.07 0 0.11 0 -0.281 0 0 0.90 0 0 0 0 0 0 0 -0.53 0 0.15 0 0 01 -0.61 0 -0.68 0 0 0 0 0 0.28 0 0.11 0 -1.07 0 -0.06-2 0 0 0 -1.17 0 0.80 0 0 0 0 0 0 0 0 0-1 0 0 0 0 -1.63 0 0.08 0 0 0 0 0 0 0 02 0 0 0 0 0.80 0 1.51 0 0 0 0 0 0 0 0 01 0 0 0 0 0.08 0 -1.63 0 0 0 0 0 0 0 02 0 0 0 0 0 0 0 -1.42 0 0 0 0 0 0 0-3 0.06 0 0.28 0 0 0 0 0 1.09 0 -0.07 0 0.28 0 0-2 0 -0.53 0 0 0 0 0 0 0 -1.08 0 0.31 0 0 0-1 -1.07 0 0.11 0 0 0 0 0 -0.07 0 0.33 0 0.47 0 -0.283 0 0 0.15 0 0 0 0 0 0 0 0.31 0 1.30 0 0 01 0.11 0 -1.07 0 0 0 0 0 0.28 0 0.47 0 0.33 0 0.072 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.26 03 -0.28 0 -0.06 0 0 0 0 0 0 0 -0.28 0 0.07 0 1.09 Jahn-Teller distortions
The JT vibrational mode Q , Q and Q , are definedfollowing the standard Van Vleck notation [26]: Q = dx + dy + dz √ Q = dx − dy Q = dx + dy − dz √ θ = tan − (cid:16) Q Q (cid:17) (18) Q = (cid:113) Q + Q , (19)where dx (dy, dz) is x i - x (y i - y , z i - z ) i.e. thedifference between the distorted Os-O bond length andthe undistorted one.To construct the phase diagram we have introduced therelative JT parameter ˜ Q , defined as the ratio betweenQ for the interpolated structure and the correspondingvalue for the optimal one. The values of Q , Q and Q ,as well as of θ and Q as a function of ˜ Q are in Table V. Phase Diagram
The phase diagram was derived by combining ab initiocalculations and linear interpolation of the PSH param-eters. The first-principles calculations were done by relaxingthe structures with the cAFM configuration for U=2.8,3.0, 3.2 and 3.4 eV. Then, all the local structural param-eters were linearly interpolated from the distorted to theundistorted (experimental-high-temperature) case. Theparameters that are interpolated include the xy anglesand the variation of the bond length which define the ˜ Q parameter in Figure 3(c) of the main text.We have performed a linear interpolation of the fittingparameters, by linearly interpolating the fitting coeffi-cients from the distorted to the undistorted structure.Being a linear model suitable for the transition (See Fig-ure 3(b) of the main text), we extrapolated the full phasediagram at higher values of ˜ Q .2 Atoms x y zNa 0 0 00.5 0.5 0.5Os 0.5 0.5 00 0 0.5Ba 0.5 0 0.750 0.5 0.750.5 0 0.250 0.5 0.25O 0.27126 0.26896 0.00.72874 0.73104 0.00.72667 0.27112 0.00.27333 0.72888 0.00.22667 0.22888 0.50.77333 0.77112 0.50.77126 0.23104 0.50.22874 0.76896 0.50 0 0.728490 0 0.271510.5 0.5 0.771510.5 0.5 0.22849TABLE IV. Optimized structural positions in direct coordi-nates.TABLE V. Values of Q , Q and Q as a function of ˜ Q inplane 1 & 2.˜ Q Q (˚A) Q (˚A) Q (˚A) θ (deg) Q (˚A) Q (˚A) Q (˚A) Q (˚A) θ (deg) Q (˚A)Plane 1 (0 0 0) Plane 2 (0 0 0.5)0 0 0 0 0 0 0 0 00.25 0.00789 0.00219 0.00041 79.35 0.00222 0.00789 -0.00219 0.00041 -79.35 0.002220.5 0.01579 0.00437 0.00082 79.35 0.00445 0.01579 -0.00437 0.00082 -79.35 0.004450.75 0.02369 0.00656 0.00123 79.35 0.00667 0.02369 -0.00656 0.00123 -79.35 0.006671 0.03159 0.00875 0.00164 79.35 0.00890 0.03159 -0.00875 0.00164 -79.35 0.00890 BAND STRUCTURE AND OPTICAL SPECTRA
The Dirac-Mott nature of this compound is confirmedby our ab-initio calculations with DFT, DFT+SOC,DFT+U and DFT+U+SOC methods (see Figure 7).Even with a high U of 5 eV we are not able to open thegap within the DFT+U scheme and only by including therelativistic effects and a U larger than 1.5 eV, is possibleto open the gap.The band structure within DFT+U+SOC does notshow significant changes upon different magnetic config-urations, as well as by changing the local structure. Thecorresponding optical spectra at DFT+U*SOC level isshown in Fig. 8.
PSEUDO SPIN HAMILTONIAN
The exchange interactions were studied by mappingthe DFT total energies onto Pseudo Spin Hamiltonian(PSH). The PSH used is given as a multipolar expansionof tensor products of dipolar, quadrupolar and octupolaroperators that can be written, in its general form, as H ij = (cid:88) K,Q I K i ,K j Q i ,Q j ( i, j ) O K i Q i ( J ) O K j Q j ( J ) , (20)where i,j label ionic sites, K ≤ J + 1, Q = − K, ..., K and O KQ are multipolar tensor operators written in thespherical harmonic superbasis. The rank of such oper-ators is given by K; where for K = 1 we have a dipole ,K=2 identifies a quadrupole and K=3 the octupole . Thecorresponding operators are written in Table VI.The PSH calculation proceeded in two steps. The firstwas to write the Pseudo Spin Matrices for a J = 3/2Pauli pseudo spin and to project them onto the directionof the osmium magnetic moment. Due to the cantingprocedure described above, the coordinate system usedis the global one of Figure 9, and therefore, the unitvectors that describe the two osmium pseudo spins aregiven byˆ n = ( sin ( α ) , cos ( α ) ,
0) ˆ n = ( cos ( α ) , sin ( α ) , . (21)where α is defined in Figure 9 and it is related to φ via α = 45 o − φ . Once this projection has been done, the newpseudo spin matrix was diagonalized and the eigenstatewith highest eigenvalue was chosen.Then we applied the mean field analysis, summarizedin equation O (1) O (2) ≈ (cid:104) O (1) (cid:105) O (2) + O (1) (cid:104) O (2) (cid:105) −(cid:104) O (1) (cid:105)(cid:104) O (2) (cid:105) , which allows to rewrite the total energyas a function of the mean values over pseudo eigenstatesof the tensors operators, leading to E = 12 (cid:88) (cid:104) i,j (cid:105) (cid:88) K,Q I K i ,K j Q i ,Q j ( i, j ) (cid:104) O K i Q i ( J ) (cid:105)(cid:104) O K j Q j ( J ) (cid:105) . (22) From Table VI it is clear that every tensor containspowers of j x , j y and j z , and therefore, powers of thecanting angle. Due to time reversal symmetry, only theeven powers of pseudo spins are allowed, and, throughalgebraic simplifications, Eq. (22) can be rewritten as E ( φ ) = A + N (cid:88) n =1 A n cos (2 nφ ) + M (cid:88) m =1 B m sin (2 mφ ) . (23)Other restrictions come from the symmetry of the mag-netic cell, which allows symmetric exchange-interactionmatrices only in-plane (i.e. when Os(1) and Os(2) belongto the same plane), while anti-symmetric exchanges areallowed between out-of-plane interactions.The energy can be further separated into E = A + E d − d + E q − q + E d − o + E o − o , (24)where E d − d is the dipole-dipole contribution, E q − q is thequadrupole-quadrupole contribution, E d − o is the dipole-octupole contribution, E o − o is the octupole-octupole con-tribution. Converting Eq. (24) into angular dependenciesgives E ( φ ) = A + (cid:88) n =1 A n cos (2 nφ ) + (cid:88) m =1 B m sin (2 mφ ) , (25)where the coefficients A n are symmetric, while the coef-ficients B m come from broken inversion symmetry. Comparison DFT vs HI vs ED
The comparison between DFT and HI/ED wasachieved by reproducing the DFT+U+SOC cantingcurve starting from the calculated exchange coefficientsobtained by HI and ED. This was done by reducing thefull 15 ×
15 HI and ED interaction matrices to the rele-vant coefficients acting within the xy plane canting, basedon a pseudo spin analysis. As previously specified, themean values (cid:104) O K i Q i ( J ) (cid:105)(cid:104) O K j Q j ( J ) (cid:105) are powers of cos ( φ ) and sin ( φ ). These mean values were multiplied by the corre-sponding coefficients, collected and added, so to obtainthe final expression with coefficients A n and B m .As mentioned in the main text, in order to achievea quantitatively consistent comparison between DFT+Uand DFT+HI, the DFT+U+SOC AFM ( φ =90) -FM ( φ = 0) energy difference has been aligned tothe corresponding DFT+HI value. Specifically, theDFT+U+SOC curves (and therefore the correspondingmagnetic interactions) have been divided by a factorof 10. The ED results have not been rescaled, whichguarantees a full consistency between DFT+HI and EDdata, manifested by the very similar interaction matri-ces (see Tab. II and Tab. III) and corresponding curvesand exchange coefficients discussed in the main text. In4 FIG. 7. Band structure of BNOO with DFT (a), DFT+SOC (b), DFT+U (c) and DFT+U+ SOC (d). U is 3.4 eV.FIG. 8. DFT+U+SOC calculated optical spectra for Ba NAaOsO compared with the ED prediction of the main Hubbardpeak (dominated by t g − t g transitions.).FIG. 9. Global reference frame for the canting angles. t g Wannier functions inDFT+HubI and ED; moreover, (ii) while the DFT+Umethod calculates the total energies of competing mag-netic phases, the ED and HubI-based approaches extracteffective Hamiltonians by a perturbation-theory treat-ment of the local-moment paramagnetic phase.Despite this fundamental difference between these twoapproaches, for localized d system (in particular 3 d andto a lesser extent 4 d ) there is generally a good corre-spondence between these two basis and the results ob-tained are typically consistent, see for instance the re-sults for SrTcO [57, 58] and NiO [59, 60]. However, formore extended d orbitals, where hybridization effects aregenerally large the degree of comparability between thetwo models is less robust. In effective Hamiltonian ap-proaches the occupation of the orbital subspace is welldefined and corresponds to the ’formal charge’ (in ourcase one electron in the t g orbitals), whereas in DFTthe ’full occupation’ is difficult to identify objectivelyas it depends on the choice of the integrated spheresand hybridization effects. This ambiguity between ’for-mal charge’ and (DFT) ’full charge’ is particularly evi-dent in Os-double perovskites as discussed by Pickett forBaCaOs O (BCOO) [13]. In BCOO the formal chargeis 2, but DFT leads to a full occupation of about 5 elec-trons, due to the highly dispersive character of the Os5 d orbital and the large degree of hybridization with theO- p states. In our calculations we fall precisely in thissituation as the full charge computed by VASP is 6 elec-trons, significantly larger than the corresponding formalcharge considered in the model calculations (1 electron).Provided these significant differences in the cor-related subspace occupancies between VASP andDFT+HubI/ED, the resulting large difference in the or-dered state energetics is not surprising. We would like tounderline that formal and full d -shell occupation are vir-tually identical in NiO and (to a lesser extent) in SrTcO ;in the latter DFT+U occupancy is 3.43 as compared tothe formal 4 d occupancy of 3 [57]. The similarity of for-mal and full d -shell occupancies in those cases results inthe exchange energies being very similar in the both ap-proaches. Conversely, for the 5 d oxide UO quantitativedifference between DFT+U+SOC and HI were alreadyreported in literature [31, 61].Support on the quantitative validity of the effectiveHamiltonian approaches can be inferred from the com-parison of calculated and experimental N´eel temperature T N . The mean-field value of T N =17.5 K obtained fromthe calculated FT-HI interactions for the undistortedphase is in a reasonable agreement, taking into accountthe usual mean-field overestimation, with its experimen- tal value of 6.7 K. For lattice distortions to occur at theAFM transition, the ordering energy for the distortedstructure must be larger than that for the undistortedone. Hence, one would expect the ordering energy (and,similarly, T N ) for the distorted phase to be at least notsmaller than that of the undistorted one. The orderingenergy and (to a lesser extent) transition temperature ofthe undistorted case thus provide a lower limit for thosequantitites in the distorted case. Since undistorted FT-HI T N already has a correct (and even somewhat overes-timated) magnitude, one would not expect the magneticenergy to be significantly larger than that given by FT-HI(and ED).From these arguments the overall magnitude of FT-HI/ED interactions seems to be in a better correspon-dence to experiment as compared to DFT+U, hence, wenormalized the magnitude of DFT+U interactions ac-cordingly. However, notice that all three approaches pre-dict virtually the same non-trivial shape of magnetic en-ergy vs. angle φ (Fig. 2b). Correspondingly, the relativemagnitudes of coefficients A , A and A in the undis-torted phase also agree well. Basing on these observa-tions one may be confident that the qualitative shape ofmagnetic energy vs. φ and the relative magnitude of vari-ous multipolar contributions into it are correctly resolvedby the DFT+U method.6 TABLE VI.
Multipolar operators for the BNOO Hamiltonian.
The Multipole moments within a cubic Γ quartet. Barsover symbols correspond to permutations of dipole operators, e.g. j x ( j y ) = j x j y j y + j y j x j y + j y j y j x .Moment Spherical Superbasis Moment Spherical SuperbasisDipole O − = j y Octupole O − = (cid:112) / j x j y − j y ) O = j z O − = (cid:112) / j x j y j z O = j x O − = (cid:112) / / j z j y − / j x j y − j y )Quadrupole O − = √ / j x j y O = j z − / j x j z + j y j z ) O − = √ / j y j z O = (cid:112) / / j z j x − / j y j x − j x ) O = [3( j z ) − j ( j + 1)] / O = (cid:112) /