# Fileset

[PhysRevB.110.205116.pdf](https://mdr.nims.go.jp/filesets/998981a6-879e-485c-87bf-dce3f54c554f/download)

## Creator

[I. V. Solovyev](https://orcid.org/0000-0002-2010-9877), R. Ono, S. A. Nikolaev

## Rights

Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article’s title, journal citation, and DOI.[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Ferromagnetic ferroelectricity due to the Kugel-Khomskii mechanism of orbital ordering assisted by atomic Hund's second rule effects](https://mdr.nims.go.jp/datasets/b28dbc54-790c-41f9-a038-b432068cc940)

## Fulltext

Ferromagnetic ferroelectricity due to the Kugel-Khomskii mechanism of orbital ordering assisted by atomic Hund&apos;s second rule effectsPHYSICAL REVIEW B 110, 205116 (2024)Editors’ Suggestion Featured in PhysicsFerromagnetic ferroelectricity due to the Kugel-Khomskii mechanism of orbitalordering assisted by atomic Hund’s second rule effectsI. V. Solovyev ,1,* R. Ono ,1 and S. A. Nikolaev 21Research Center for Materials Nanoarchitectonics (MANA), National Institute for Materials Science (NIMS),1-1 Namiki, Tsukuba, Ibaraki 305-0044, Japan2Department of Materials Engineering Science, The University of Osaka, Toyonaka 560-8531, Japan(Received 30 May 2024; accepted 26 August 2024; published 7 November 2024)The exchange interactions in insulators depend on the orbital state of magnetic ions, obeying certain phe-nomenological principles, known as Goodenough-Kanamori-Anderson rules. Particularly, the ferro order ofalike orbitals tends to stabilize antiferromagnetic interactions, while the antiferro order of unlike orbitals favorsferromagnetic interactions. The Kugel-Khomskii theory provides a universal view on such coupling betweenspin and orbital degrees of freedom, based on the superexchange processes: namely, for a given magnetic order,the occupied orbitals tend to arrange in a way to further minimize the exchange energy. Then, if two magneticsites are connected by the spatial inversion, the antiferro orbital order should lead to the ferromagnetic couplingand break the inversion symmetry. This constitutes the basic idea of our work, which provides a pathwayfor designing ferromagnetic ferroelectrics: the rare but fundamentally and practically important multiferroicmaterials. After illustrating the basic idea on toy-model examples, we propose that such behavior can be indeedrealized in the van der Waals ferromagnet VI3, employing for this analysis the realistic model derived fromfirst-principles calculations for magnetic 3d bands. We argue that the intra-atomic interactions responsiblefor Hund’s second rule, acting against the crystal field, tend to restore the orbital degeneracy of the ionic d2state in VI3 and, thus, provide a necessary flexibility for activating the Kugel-Khomskii mechanism of theorbital ordering. In the honeycomb lattice, this orbital ordering breaks the inversion symmetry, stabilizing theferromagnetic-ferroelectric ground state. The symmetry breaking leads to the canting of magnetization, whichcan be further controlled by the magnetic field, producing a huge change of electric polarization.DOI: 10.1103/PhysRevB.110.205116I. INTRODUCTIONIn a broad sense, multiferroics are materials, where theferroelectric (FE) order can coexist with a magnetic one [1,2].These are the key material systems for achieving the cross-control of magnetic and electric properties by applying anelectric or magnetic field [3]. Nevertheless, literally, the mul-tiferroicity implies a somewhat narrower requirement: bothorders should be of the ferro type, so that the material is notsimply magnetic but ferromagnetic [4–6]. This is particularlyimportant for the cross-control applications: if the ferromag-netic (FM) moment, M, is finite and preferably large, it can bemanipulated by a relatively weak magnetic field. The sameholds for the ferroelectric polarization, P, and the electricfield. Thus, from the practical point of view, it is desirableto have materials with large M and P, and strong couplingbetween them. However, the ferroelectricity and ferromag-netism obeys very different principles and very rarely coexist*Contact author: SOLOVYEV.Igor@nims.go.jpPublished by the American Physical Society under the terms of theCreative Commons Attribution 4.0 International license. Furtherdistribution of this work must maintain attribution to the author(s)and the published article’s title, journal citation, and DOI.in nature. Namely, the ferroelectricity implies breaking of theinversion symmetry. However, it cannot be achieved by a sim-ple FM arrangement of spins, which has the same symmetryas the crystallographic one. On the other hand, if the inversionsymmetry breaking results from the intrinsic instability of thecrystal structure, there is no guarantee that the correspondingto it magnetic structure will be ferromagnetic. In fact, most ofinsulating transition-metal oxides are antiferromagnetic.Therefore, the main attention is paid to creation of ar-tificial materials, which would combine the FE and FMcharacteristics within one sample or device [5]. One possibledirection is the synthesis of heterostructures, consisting ofFE and FM layers of two different materials [7]. Anotherpromising direction is the strain engineering. Particularly,some transition-metal oxides can turn into the FE-FM stateby epitaxial strain [8–10]. The main driving force is the in-trinsic FE instability of the so-called d0 materials, related tothe coupling between the occupied bonding and unoccupiedantibonding states of opposite parity [11,12]. For instance,the coupling between the occupied O 2p and unoccupied Ti3d bands in EuTiO3, caused by the FE displacements, canlower the energy [8]. Moreover, the magnetic Eu2+ ions alterthis coupling, making it dependent on magnetic structure ofthe Eu sublattice. Thus, although cubic EuTiO3 is the para-electric antiferromagnet, the epitaxial strain can turn it intothe FE-FM state [8,9]. The partial occupation of antibonding2469-9950/2024/110(20)/205116(14) 205116-1 Published by the American Physical Societyhttps://orcid.org/0000-0002-2010-9877https://orcid.org/0000-0002-1529-3761https://orcid.org/0000-0001-7481-1007https://ror.org/026v1ze26https://ror.org/035t8zc32https://crossmark.crossref.org/dialog/?doi=10.1103/PhysRevB.110.205116&domain=pdf&date_stamp=2024-11-07https://doi.org/10.1103/PhysRevB.110.205116https://creativecommons.org/licenses/by/4.0/I. V. SOLOVYEV, R. ONO, AND S. A. NIKOLAEV PHYSICAL REVIEW B 110, 205116 (2024)transition-metal 3d states weakens the FE instability. Nev-ertheless, the effect can still persist for certain 3d configu-rations, such as d3 [12], as was theoretically proposed forSrMnO3 [10,13], where the same Mn3+ ions are responsiblefor magnetism and participate in the FE displacements, thus,resulting in stronger spin-lattice coupling and larger magnetictransition temperature in comparison with EuTiO3.In this article we propose a completely new and so farunexplored route for designing ferroelectric (or polar) fer-romagnets, which is based on the Kugel-Khomskii (KK)mechanism of the orbital ordering [14].The interatomic exchange interactions between spins de-pend on the orbital state of atoms participating in theseexchange processes: which orbitals are occupied, which areempty, and how they are oriented relatively to each other inthe magnetic bonds, i.e., what is commonly called the orbitalordering [14]. The basic rules describing the character of theseinteractions in insulators are widely know as Goodenough-Kanamori-Anderson (GKA) rules [15–18]. Particularly, theferro orbital order, where electrons occupy the same orbitals,typically leads to the antiferromagnetic (AFM) coupling be-tween the spins. On the other hand, the antiferro orbital order,where occupied orbitals alternate on the lattice, usually fa-vors the FM interactions. These fundamental principles werefurther elaborated by KK [14,19,20] on the basis of superex-change (SE) theory [21], resulting in what is now called theKK mechanism of the orbital ordering, which states that fora given spin order, the orbital degrees of freedom will tend torelax in the direction to further minimize the exchange energy.The KK mechanism was proposed long before the currentera of multiferroic materials and so far has not been consid-ered as a possible source of breaking the inversion symmetry.Typical applications of the KK mechanism are focused on theanalysis of spin and orbital phenomena in compounds, wheremagnetic sites are located in the inversion centers and thematerials remain centrosymmetric irrespectively of the spinor orbital order [22], as in colossal magnetoresistive man-ganites [23] or other perovskite transition-metal oxides [24].In fact, many of these materials do exhibit the antiferroorbital order, which is responsible for the FM character of ex-change interactions, as it happens, for instance, in YTiO3 [24],LaMnO3 [14,16], or BiMnO3 [25,26]. However, the existenceof inversion centers makes most of them antipolar [27].What if the inversion center is located between two mag-netic sites? Then, one can expect that the antiferro orbitalorder across the inversion center will lead to the FM interac-tions between the spins, as required by GKA rules, and breakthe inversion symmetry, giving us a unique possibility for real-izing simultaneously the ferromagnetism and ferroelectricitywithin one phase. This constitutes the main idea of our work,which will be elaborated as follows.First, in Sec. II, we will explore this basic idea by consid-ering toy-model examples of degenerate yz and zx orbitals inthe zigzag chain and honeycomb lattice, where the problemcan be solved analytically providing a transparent expressionfor the exchange energy, which explains the emergence ofthe antiferro orbital order and electric polarization. The keyaspect of the zigzag chain and honeycomb lattice is that bothof them are centrosymmetric. However, the inversion centersare located in the mid-points connecting two magnetic sites.FIG. 1. (a) Fragment of the crystal structure of VI3: each V atomsis surrounded by six I atoms, forming the hexagon of edge-sharingVI6 octahedra. The vanadium sublattices, which are transformed toeach other by the inversion operation are shown by different colorsand denoted as V1 and V2. The inversion centers are denoted by+. (b) Stacking of the honeycomb planes. (c) Densities of states(DOS) in the local-density approximation. Shaded areas show partialcontributions of the V 3d states. The Fermi level is at zero energy.Therefore, these are the structures where the antiferro orbitalorder will simultaneously break the inversion symmetry andstabilize the FM ground state. The role of charge and orbitaldegrees of freedom in assisting the multiferroic behavior wasknown before [1,28]. The new aspect of our proposal is thatthe orbital ordering alone can be the source of both ferroelec-tricity and ferromagnetism.Then, in Sec. III, we will turn to a realistic example ofVI3, which has attracted a considerable attention as a newlayered FM semiconductor with relatively high Curie tem-perature TC � 50 K [29]. The main structural motif of thisquasi-two-dimensional van der Waals ferromagnet is againthe honeycomb planes (Fig. 1). According to formal valencearguments, each V site has two 3d electrons. In the octahe-dral environment they populate two out of three t2g orbitals,indicating the importance of orbital degrees of freedom inthe physics of VI3. The main question is, however, how wellthese orbital degrees of freedom are quenched by the localdistortions of the VI6 octahedra. Indeed, the distortions willtend to split the t2g levels. The fundamental Jahn-Teller the-orem states in this respect that the splitting should lift theorbital degeneracy in the direction to form a nondegenerateground state [30]. Nevertheless, if the splitting is small, otheringredients can come into play. Particularly, two electrons inthe 3d shell are subjected to Hund’s rule effects, which actin the opposite direction and tend to reenforce the groundstate with maximal multiplicity. The corresponding energygain is controlled by the Racah parameter B [31,32]. Usingelectronic structure calculations based on density functionaltheory (DFT), we will evaluate relevant parameters and showthat B in VI3 is sufficiently large to overcome the crystal-fieldsplitting and activate the KK mechanism of the orbital order-ing, as it will follow from the analysis of atomic multipletstructure and dynamical mean-field theory (DMFT) calcula-tions on the honeycomb lattice [33,34]. Then, we will showthat for the realistic parameter range, the antiferro orbitalorder can be indeed established in VI3, resulting in the FM-FE ground state. The relativistic spin-orbit (SO) interactioninterplays with the symmetry breaking caused by the orbitalordering, resulting in a canted magnetic structure, which canbe further controlled by the magnetic field, leading to a hugechange of electric polarization.205116-2FERROMAGNETIC FERROELECTRICITY DUE TO THE … PHYSICAL REVIEW B 110, 205116 (2024)Finally, in Sec. IV, we will summarize our results, dis-cussing their implications to the properties of VI3 as well asmore general aspects of the Hund’s rule physics in solids.II. TOY-MODEL CONSIDERATIONSThe goal of this section is to illustrate the basic idea ofinversion symmetry breaking by the orbital ordering, result-ing in coexistence of ferroelectricity and ferromagnetism. Forthese purposes we consider toy-model examples of degenerateyz and zx orbitals in the zigzag chain and honeycomb lattice.We do not aim to find the correct ground state of the Hubbardmodel in a certain parameter range, which is an interestingproblem on its own right. There may be other possible candi-dates for the ground state, including spin-dimerized or orderedcomplex harmonic states [22]. The analysis of these states isbeyond the scopes of our work. Nevertheless, we would like toemphasize that the scenario of ferromagnetic ferroelectricity,which we propose, should be seriously considered amongother possible solutions of the Hubbard model.A. Ordering of the yz and zx orbitals in the zigzag chainThe simplest model, which explains the basic physics ofhow the KK mechanism can break the inversion symmetry andinduce the electric polarization is the one-dimensional zigzagchain (see Fig. 2). In this case, there are two sites in the unitcell (1 and 2), which can be transformed to each other by thespatial inversion about the midpoint of the bond, connectingthese two sites. Let us assume that there is only one electronper site, which is shared by two atomic states, yz and zx. Thus,in the atomic limit, the ground state is degenerate. Then, theelectron hoppings t̂i j are such that in the neighboring bondsthey will connect zx with zx in the direction x and yz withyz in the direction y [35]: t zx,zxi j||x = t yz,yzi j||y = t . This hopping liftsthe degeneracy, ordering the orbitals in the alternating way, asexplained in Fig. 2, which minimizes the energy of SE interac-tions for the FM state [14]. The same orbital ordering makesthe atomic sites inequivalent and, thus, breaks the inversionFIG. 2. Ordering of the yz and zx orbitals in the zigzag chain,breaking the inversion symmetry and stabilizing the ferromagneticcoupling: side view (upper panel) and top view (lower panel). Theelectron densities across the inversion centers (denoted by +) areplotted by different colors: larger objects are the densities in theatomic limit and smaller objects are the densities transferred fromthe neighboring sites due to the superexchange processes in thedirections, which are shown by arrows. a0 is the lattice parameter.symmetry. The corresponding electric polarization can beevaluated along the same line as in the theory of SE interac-tions [36,37] but starting for these purposes with the generalexpression for P in periodic systems, formulated in terms ofthe Wannier functions [38–40]. Namely, if |αoi 〉 is the occupiedWannier function at site i in the atomic limit, t̂ will induce thetail of this orbital, |αoi→ j〉, spreading to the neighboring site j.It can be evaluated by treating t̂ as a perturbation, the same asin the SE theory, which yields |αoi→ j〉 = − 1�|αuj 〉〈αuj |t̂ ji|αoi 〉,where � is the proper combination of on-site Coulomb re-pulsion U and intra-atomic exchange interaction J describingthe splitting of occupied and unoccupied states with the samespin and |αuj 〉 is the unoccupied Wannier function. It is as-sumed that J is sufficiently large so that the hopping processesresulting in the AFM exchange coupling can be neglected.This yields P||a = ∓e( t�)2 [41], where two signs stand forthe orbital ordering depicted in Figs. 2(a) (+) and 2(b) (−),and e is the minus electron charge. In the one-dimensionalcase, P||a is nothing but the edge charge [39]. Similar modelwas considered in Refs. [36,43] to explain the emergence ofelectric polarization in the E-phase of manganites. The maindifference is that, in manganites, the orbital order is driven bythe Jahn-Teller distortion, which is an external factor in theconsidered electronic model, while here it originates solelyfrom the SE interactions and, formally, no distortion is neededto break the inversion symmetry.B. Ordering of the yz and zx orbitals in the honeycomb planeNow we turn to a more realistic model of the honeycombplane, which may have some relevance to realistic materi-als, such as TiCl3 [44]. In the honeycomb lattice, there arealso two sites in the unit cell, which can be transformed toeach other by the spatial inversion [see Fig. 3(a)]. Again, weassume that there are only two orbitals, yz and zx, and oneelectron per site. The transfer integral operates only betweenorbitals, which are parallel to the bond. For instance, for thebond 1-2 in Fig. 3, these are zx orbitals. The transfer integralsin other bonds can be obtained by threefold rotations, asexplained in the Supplemental Material [41].FIG. 3. Ordering of yz and zx orbitals in the honeycomb plane,breaking the inversion symmetry and stabilizing ferromagnetic inter-actions. The electron densities across the inversion centers (denotedby +) are plotted by different colors. The unit cell is shown by dashedline. a0 is the lattice parameter. (a) Ideal lattice. (b) Distorted lattice,where β �= β ′ (as explained in the inset).205116-3I. V. SOLOVYEV, R. ONO, AND S. A. NIKOLAEV PHYSICAL REVIEW B 110, 205116 (2024)Defining at each site (ν) of the unit cell the occupied (o)and unoccupied (u) orbitals as∣∣αoν〉 = cos φν |yz〉 + sin φν |xz〉,∣∣αuν〉 = − sin φν |yz〉 + cos φν |xz〉,and treating the transfer integrals as a perturbation, it isstraightforward to find the following expression for the totalenergy change [41]:E = −3t24�(2 − cos 2(φ1 − φ2)),which takes the minimum if φ2 = φ1 + π2 mod π . Thus, oneof the angles, φ1, remain unspecified. Then, the polarizationis given byP ≡(PxPy)= ea0(t�)2(− cos 2φ1sin 2φ1),where a0 is the lattice parameter. Thus, the orbital order-ing breaks the spatial inversion, yielding finite |P| = ea0( t�)2.However, the ground state remains degenerate and P can haveany direction in the xy plane, depending on the angle φ1.There can be several scenarios of lifting the degeneracy.For instance, in a more general case, where there are severalunoccupied states, such a degeneracy does not occur (seeSec. III D). In the simplest two-orbital model, consideredhere, the value of φ1 can be decided by the exchange stric-tion effects. Particularly, the orbital ordering in Fig. 3(a) willmake the bond 1-2′ and 1-2′′ different from the bond 1-2 andthe structure will tend to relax in order to further minimizethe energy change. Here, we assume that such deformationof the honeycomb plane can be described by the angle β =2π3 + δβ, formed by the bonds 1-2′ and 1-2′′ with the bonds1-2, which is different from the angle β ′ = 2π3 − 2δβ, formedby the bonds 1-2′ and 1-2′′ with each other, while the bondlengths are assumed to be the same. δβ = 0 corresponds to theundistorted structure. The situation is explained in Fig. 3(b).The transfer integrals in the bond 1-2 operate only betweenzx orbitals, while the ones in the bonds 1-2′ and 1-2′′ areobtained by considering rotations of 1-2 by the angle ∓β [41].For small δβ, the corresponding energy change is given byE = −3t24�(2 − cos 2(φ1 − φ2) + 4√3δβ cos 2(φ1 + φ2)).Then, δβ > 0 strengthens the antiferro orbital order in thebonds 1-2′ and 1-2′′. In this case, the second and third termsin (. . . ) can be minimized independently, yielding φ1 = 0 andφ2 = π2 . Thus, the degeneracy is lifted and the polarization isparallel to the x axis.The considered models dealing with the d1 systems impliesthat the degenerate yz and zx states are split off by the crystalfield to become the lowest energy atomic states, which accom-modate a single electron. Although these models are easy tosolve, and in this sense can be very insightful, they are hardlypractical. The main obstacle for the practical realization ofthe considered scenarios is the direction of the crystal field,which typically acts to form a nondegenerate ground state, asit is required by the Jahn-Teller theorem [30]. For instance,one possible d1 candidate to form the antiferro orbital orderin the honeycomb lattice is TiCl3 [44]. However, the crystalfield in TiCl3 tends to stabilize the nondegenerate z2 orbital,in the direction perpendicular to the honeycomb plane, whilethe degenerate states lie higher in energy. The basic limitationof the d1 system is that the crystal field is the only parameter,which can control the order of the atomic states. Due to theone-electron character of the problem, the atomic Hund’s rulessimply do not apply here. Thus, there is no way to reverse theorder of the crystal-field orbitals in the favor of the degenerateground state. In this sense, a more promising direction is toexplore d2 materials. Then, the crystal field will still tend toselect a nondegenerate ground state. Nevertheless, the newaspect of the d2 systems is that the interaction between two3d electrons is subjected to the Hund’s rule effects, whichact in the opposite direction and tend to stabilize the groundstate with maximal orbital multiplicity. In the next section, wewill argue that VI3 is indeed a good candidate for practicalrealization of such scenario.III. IMPLICATIONS TO THE PROPERTIES OF VI3VI3 exhibits the structural phase transition at Ts � 78 K,which is followed by the FM transition at TC � 50 K [29].Another structural phase transition, at around 32 K, was alsosuggested [45,46]. However, currently there is no clear con-sensus even about the crystallographic symmetry of VI3. Forinstance, the high-temperature (T > Ts) phase was proposedto have trigonal R3 (space group No. 148) [29,45,46], trigonalP31c (163) [47], and monoclinic C2/m (12) [48] symmetry.The low-temperature (T < Ts) phase was proposed to be tri-clinic P1 (2) [46,49], monoclinic C2/c (15) [47], and trigonalR3 (148) [48]. Thus, even the direction of the distortion withlowering temperatures appears to be the subject of contro-versy: some reports suggest symmetry lowering [46,47], whileanother report suggests that the symmetry becomes higher,similar to what is observed in CrI3 [48]. To a certain extent,the value of Ts can be controlled by the magnetic field [45].Moreover, the 32 K transition was also suggested to be mag-netic, due to either disappearance of the magnetic order inone of the V sublattices [50] or reorientation of the magneticmoments [49].In this section, we systematically study the symmetrybreaking in VI3 caused by the orbital ordering, starting forthese purposes with the structure with highest R3 symme-try [46]. The R3 space group can be generated by consideringthe threefold rotations about z in combination with spatialinversion. We will show that the experimentally observedbreaking of the threefold rotation symmetry can be rational-ized by considering the KK mechanism of the orbital orderingin combination with the Hund’s second rule effects. Further-more, we predict the new symmetry pattern, also driven bythe KK mechanism, where both threefold rotation and in-version symmetries are broken by the orbital ordering. Theprediction remains largely intact even for the low-symmetryexperimental P1 structure. This structure has inversion sym-metry. However, it can be again broken by the orbital order.We evaluate the electric polarization, induced by the inversionsymmetry breaking, and propose how it can be controlled bymagnetic field.205116-4FERROMAGNETIC FERROELECTRICITY DUE TO THE … PHYSICAL REVIEW B 110, 205116 (2024)A. MethodThe electronic structure calculations have been performedfor the experimental crystal structure of VI3, reported inRef. [46], and using for these purposes either linear muffin-tinorbital (LMTO) method [51,52] or the pseudopotential Quan-tum ESPRESSO (QE) method [53], which were supplementedwith, respectively, local density approximation (LDA) [54]and generalized gradient approximation (GGA) [55] for theexchange-correlation potential. After that, we construct thefive-orbital model for the magnetic V 3d bands located nearthe Fermi level [Vt2g + Veg bands in Fig. 1(c)]:Ĥ =∑i j∑σσ ′∑abtaσ,bσ ′i j ĉ†iaσ ĉ jbσ ′+ 12∑i∑σσ ′∑abcdU abcd ĉ†iaσ ĉ†icσ ′ ĉibσ ĉidσ ′ , (1)where ĉ†iaσ (ĉiaσ ) stands for the creation (annihilation) of anelectron with spin σ at the Wannier orbital a of site i. Forthe construction of Wannier functions in QE, we use themaximally localized Wannier functions method [56], as im-plemented in the WANNIER90 package [57], while in LMTOwe employ the projector-operator technique [56,58]. Then,the one-electron part in Eq. (1), t̂ = [t aσ,bσ ′i j ], is given bythe matrix elements of the LDA (GGA) Hamiltonian in theWannier basis. The dependence of t aσ,bσ ′i j on the spin indicesis due to the SO interaction, where the main contributionscome from the heavy I sites. In order to include these con-tributions, it is essential that the electronic structure for the V3d bands should be calculated with the SO coupling beforethe construction of the model Hamiltonian. Without SO cou-pling, the matrix elements t aσ,bσ ′i j become t aσ,bσ ′i j = t a,bi j δσσ ′ .Then, the crystal-field splitting is given by the site-diagonalparameters t a,bii .The screened Coulomb interactions are evaluated withinconstrained random-phase approximation [59]. In LMTO, thisis done in an approximate way, basically by considering theself-screening effect of on-site Coulomb interactions in the V3d bands by the same V 3d states, which are admixed into theI 5p and other bands [see Fig. 1(c)], as explained in Ref. [58].After that, the 5 × 5 × 5 × 5 matrix Û = [U abcd ] of screenedintra-atomic Coulomb interactions was fitted in terms of threeparameters, which would describe these interactions in aspherical atomic environment: the Coulomb repulsion U =F 0, responsible for the overall stability of atomic shell withthe given number of electrons; the intra-atomic exchange in-teraction J = (F 2 + F 4)/14, responsible for Hund’s first rule;and the Racah parameter B = (9F 2 − F 4)/441, responsiblefor Hund’s second rule (where F 0, F 2, and F 4 are Slaterintegrals) [31,32]. The obtained parameters U , J , and B arelisted in Table I. The value of U is somewhat smaller inLMTO, while the values of J and B, obtained LMTO and QEmethods, are comparable.B. Hund’s second rule and orbital degeneracyIn this section we consider the interplay between crys-tal field and Coulomb interactions in the atomic limit. Thecrystal-field orbitals, obtained after the diagonalization of theTABLE I. Parameters U , J , and B of intra-atomic Coulomb in-teractions, and trigonal splitting �tr between eπg and a1g levels (allare in eV) in the R3 phase of VI3, as obtained in the model based onthe LMTO and QE methods.Method U J B �trLMTO 1.21 0.74 0.07 0.01QE 1.97 0.64 0.06 0.03site-diagonal part of t̂ are shown in Fig. 4(a). The octahedralfield, 10Dq, splits the 3d levels into the three-dimensionalrepresentation t2g and two-dimensional representation eσg ,∣∣eσ,1g〉 = − sin α|xy〉 + cos α|zx〉,∣∣eσ,2g〉 = − sin α|x2 − y2〉 + cos α|yz〉.The t2g levels are further split by the trigonal field �tr into thedoubly-degenerate states eπg ,∣∣eπ,1g〉 = cos α|xy〉 + sin α|zx〉,∣∣eπ,2g〉 = cos α|x2 − y2〉 + sin α|yz〉,which belong to the same representation as eσg , and the non-degenerate state |a1g 〉 = |z2〉. The numerical value of α isabout 34◦. Since �tr > 0, the eπg states lie lower than a1g [seeFig. 4(d)] and accommodate both 3d electrons. Thus, fromthe viewpoint of the one-electron crystal-field splitting, theground state is expected to be nondegenerate (the so-calledeπg eπg state [60]), in agreement with the Jahn-Teller theo-rem [30], which should be satisfied at the level of LDA/GGAcalculations. Nevertheless, LDA (GGA) is an approximation,which does not properly take into account the interactionsresponsible for Hund’s second rule [61]. The latter are pro-portional to the Racah parameter B and can easily overcomethe trigonal splitting if B > �tr . Realistic estimate of the pa-rameters suggests that this situation is indeed realized in VI3:B is certainly small. However, �tr appears to be even smaller(see Table I). This trend appears to be generic as very similarbehavior was obtained for the triclinic modifications of VI3 atT = 9 and 60 K [46]: the triclinic distortion further splits thet2g levels. However, the splitting (∼10 - 20 meV [41]) remainsto be smaller than B.Thus, the intra-atomic Coulomb interactions will tend toreverse the order of eπg and a1g levels. This is due to the funda-mental property of the Hund’s rule coupling to form a groundstate with the greatest value of multiplicity. The intuitive rea-son can be understood as follow: The main contribution tothe eπg states are associated with the xy and x2 − y2 orbitals(the corresponding cos α ∼ 0.83), which are both located inthe xy plane and, therefore, experience a strong Coulombrepulsion. In order to reduce this repulsion, it is energeticallymore favorable to replace one of the occupied eπg states by thea1g state. Simple mean-field considerations can be found inSupplemental Material [41], which clearly show that the effectis indeed driven by the Racah parameter B.Results of diagonalization of the two-electron Hamilto-nian, combining the on-site Coulomb interactions with thecrystal field are summarized in Fig. 4. First, we enforce205116-5I. V. SOLOVYEV, R. ONO, AND S. A. NIKOLAEV PHYSICAL REVIEW B 110, 205116 (2024)FIG. 4. (a) Crystal-field orbitals for VI3. (b) The effect of Racah parameter B on the two-electron energies without trigonal crystal field.(c) Amplified low-energy part of (b). (d) Parameters of the crystal field splitting and their values obtained in the model based on the LMTOmethod. (e) Two-electron energies, relative to the lowest state, depending on the trigonal splitting δtr = 100(�tr/�0tr ), where �0tr = 12 meV isthe actual value obtained in the LMTO based model. The value of B is set to 70 meV.�tr = 0. Then, without B, the energy diagram is pretty muchsimilar to the one obtained by Tanabe and Sugano for thed2 configuration in octahedral environment [62] [Fig. 4(b)].Namely, the ground state is 3T1, which is ninefold degenerate(being both spin and orbital triplet). Some states, like T1 andT2 remain accidentally degenerate when B = 0. Finite B liftsthis degeneracy not only between T1 and T2, but also withineach of these states. Particularly, the 3T1 ground state is splitinto the sixfold degenerate state: the spin triplet with orbitalmagnetic quantum numbers m = ±1 and threefold degeneratespin triple state with m = 0 [Fig. 4(b)]. Then, we consider theeffect of �tr [Fig. 4(e)]: As expected, for small �tr , the stateswith m = ±1 remain lower in energy, while further increaseof �tr will eventually change the order of the m = ±1 andm = 0 states.Thus, within realistic parameters range, the ground state ofthe V3+ ion in VI3 remains to be orbitally degenerate, openinga possibility for the KK mechanism of orbital ordering, oncewe consider the transfer integrals in the honeycomb lattice.It is also instructive to compare the behavior of VI3 withV2O3. The latter compound, formally hosting the same V3+ions in the trigonal environment, was regarded as the canoni-cal S = 12 Mott insulators, where a1g electrons are dimerizedin the V-V bonds and form a singlet, while the remainingeπg electrons form an orbitally ordered state, resulting in apeculiar AFM structure [63]. However, this picture was laterrevisited on the basis of first-principles calculations, suggest-ing the predominantly eπg eπg configuration of V3+ in V2O3,corresponding to the nondegenerate m = 0 state [60]. Then,why is the behavior of V3+ in VI3 so different from V2O3? Wehave constructed the model for the V 3d bands in V2O3 usingthe experimental R3c structure at T = 175 K [64]. The cor-responding parameters, controlling the order of the m = ±1and m = 0 states, are B = 0.09 eV and �tr = 0.15 eV. Thus,although B > �tr in VI3, it appears that B < �tr in V2O3,reversing the order of states and forcing the m = 0 groundstate. Furthermore, 10Dq is another important parameter inthe problem. The Hund’s rule interactions tend to mix theeπg and eσg states, belonging to the same representation, tofurther minimize the energy of intra-atomic interactions. Thismixing acts against the octahedral field. Thus, the perspec-tives of realization of the orbitally degenerate ground statedepend not only on the values of B and �tr , but also on10Dq. In this respect, it is also important that 10Dq is rel-atively small in VI3 (10Dq = 1.49 eV) in comparison withV2O3 (10Dq = 2.31 eV, see Supplemental Material [41] fordetails).It is needless to say that the use of the five-orbital model,constructed for both the t2g and eσg bands, is essential to incor-porate the effects of Hund’s second rule. Formally, the V t2gbands in VI3 are well separated from other bands (see Fig. 1).Then, it would be straightforward to construct a more compactthree-orbital model for the Vt2g bands, which is certainly eas-ier to solve. However, there can be no separate Hund’s secondrule in the isolated t2g shell, where two-electron configurationwith S = 1 is equivalent to a noninteracting single-hole con-figuration. Therefore, such a model would be meaningless forthe purposes of our work, as it does not take into account theessential piece of physics.C. Basic electronic structureIn most of the applications, the model (1) was solved inmean-field Hartree-Fock (HF) approximation [58], replacingthe interaction part by∑i∑σσ ′∑abVσσ ′i,ab ĉσ†ia ĉσib , (2)and solving the obtained one-electron equations with themean-field potential V̂i = [Vσσ ′i,ab ] self-consistently. In orderto check the validity of the HF method, we also employthe DMFT. The details can be found in the Supplemental205116-6FERROMAGNETIC FERROELECTRICITY DUE TO THE … PHYSICAL REVIEW B 110, 205116 (2024)FIG. 5. Total and partial densities of states in the crystal-fieldrepresentation as obtained in the dynamical mean-field theory forT = 580 and 120 K in the paramagnetic state (top) and Hartree-Fock approximation for the ferromagnetic state (bottom). The Racahparameter B was set to either 0 (left) or 0.06 eV (right). Otherparameters were taken from the QE set (see Table I). The Fermi levelis at zero energy (the middle of the band gap).Material [41]. In order to reproduce the semiconducting be-havior of VI3 in the framework of DMFT, it is important to usethe QE model parameters. The smaller value of U , obtained inthe LMTO method, was insufficient to open the gap.As we will see below, the HF approximation gives us a sta-ble mean-field solution, where the orbitally ordered state has amuch lower energy compared to other possible configurations.Therefore, we believe that the mean-field approach is a reli-able starting point for the analysis of multiferroic behavior inVI3. Furthermore, the HF approximation is a more convenienttool for the analysis of electric polarization in terms of theBerry-phase theory as it can be adapted for the use of theKing-Smith-Vanderbilt formula [38–40]. On the other hand,the quantum orbital fluctuations are known to play an impor-tant role in the physics of d2 materials [65–67]. With this inmind, we employed a more sophisticated DMFT technique,which allows us to incorporate local quantum effects and be-comes exact in the limit of infinite dimensions [33,34]. WhileDMFT calculations in our work have been carried out for theparamagnetic state, the results are very reasonable and supportour main conclusions derived from the HF approximation.The corresponding densities of states are shown in Fig. 5.First, the HF method captures the main tendencies of two-electron calculations in the atomic limit considered in theprevious section. Namely, without B, the doubly degeneratemajority-spin eπg states are occupied, while the nondegeneratea1g states are located in the unoccupied part of the spec-trum, which is totally consistent with the fact that for B = 0the occupation of atomic states is controlled solely by thecrystal field. Nevertheless, the situation changes dramaticallywhen we switch on B. Then, the degeneracy of the eπg statesis lifted, splitting them into the occupied and unoccupiedbands. Furthermore, the a1g states become occupied. Thistendency is supported by DMFT calculations for the paramag-netic state. However, these calculations are performed at finitetemperatures, which additionally affect the distribution ofthe atomic states. Particularly, T = 580, 230, and 120 Kused in the calculations (all are far above TC) correspond tokBT = 0.05, 0.02, and 0.01 eV. The results for T = 580 and120 K are shown in Fig. 5 and the ones for T = 230 K arediscussed in the Supplemental Material [41]. The first value ofkBT is larger than �tr = 0.03 eV, while two other values aresmaller than �tr (but all the values are comparable to �tr andB = 0.06 eV). This readily explains the fact that even for B =0 there is a finite weight of the a1g states in the occupied part.Nevertheless, this weight clearly decreases with the decreaseof T and practically vanishes for T = 120 K. Therefore, itcan be attributed to the finite temperature effects. On the otherhand, when B is finite, the same weight of the occupied a1gstates practically does not depend on T , meaning that in thiscase it is the feature of Hund’s second rule.D. Orbital ordering and magnetic ground stateIn order to study the orbital ordering in VI3, we turn tothe HF calculations. The use of the QE set of model pa-rameters was important at the level of DMFT calculationsto reproduce the semiconducting character of VI3. However,the HF approximation tends to overestimate the band gapin comparison with DMFT (see Fig. 5). To be specific, theenergy gap in DMFT is well consistent with the experimentalvalue of 0.6 eV [29], while the HF band gap, 0.9 eV, is clearlyoverestimated. Therefore, we believe that at the level of HFcalculations, it would be more reasonable to use a smallervalue of U in order to mimic the band gap obtained in DMFT.This is the main reason why we switch to the LMTO setof model parameters (Table I). The corresponding electronicstructure can be found in the Supplemental Material [41].The orbital ordering (the distribution of electron densitiesaround V sites) obtained in the HF approach for the FM stateis displayed in Fig. 6. Particularly, we have performed threetypes of calculations (by properly averaging the HF potentialin the process of calculations): (i) by enforcing the original R3symmetry of the lattice, including threefold rotation about zand inversion symmetries; (ii) by enforcing only the inversionsymmetry and treating two V sites in the unit cell as equivalent(the space group P1); (iii) by fully relaxing the symmetryand allowing a different shape of the electron density aroundtwo V sites, which are crystallographically connected by theFIG. 6. Orbital ordering in the ferromagnetic state of rhombohe-dral R3 phase of VI3 as obtained in the Hartree-Fock calculations byenforcing the original trigonal R3 symmetry, the triclinic P1 sym-metry, and fully relaxing the symmetry (P1). The crystallographicinversion centers are denoted by +. �E is the corresponding energychange relative to the ferromagnetic R3 state.205116-7I. V. SOLOVYEV, R. ONO, AND S. A. NIKOLAEV PHYSICAL REVIEW B 110, 205116 (2024)FIG. 7. Spin-wave stiffness tensors (in meVÅ2) for the orbitalstates of the R3, P1, and P1 symmetry, and corresponding spin-wavedispersions near the  point of the Brillouin zone.spatial inversion (the space group P1). In the R3 case, bothelectrons are forced to occupy two eπg states, resulting in theeπg eπg configuration with small admixture of the eσg states dueto the hybridization effects. In the P1 case, one of the occupiedstates is eπg while another one is a1g (the configuration eπg a1g).Nevertheless, the occupied eπg states are the same at the sitesV1 and V2 (see Fig. 1). Thus, the threefold rotational symme-try is broken, but the spatial inversion is preserved. In the P1case, the occupied eπg states are different at the sites V1 andV2, thus breaking both threefold rotational and inversion sym-metries. The total energy steadily decreases in the directionR3 → P1 → P1. Thus, the R3 phase experiences the internalinstability due to the orbital ordering, which tends to lower thesymmetry.By enforcing B = 0 in the HF calculations, we were ableto obtain only one solution corresponding to the R3 symmetry.Other solutions with the P1 and P1 symmetries, obtained forB �= 0, steadily converge to the R3 one after setting B = 0.Thus, the Hund’s second rule effects play a crucial role inbreaking the symmetry and establishing the antiferro orbitalorder in the P1 phase. Yet, the precise meaning of “antiferroorbital order” in this context needs to be clarified, because theoccupied eπg orbitals in the vanadium sublattices V1 and V2are not orthogonal and, strictly speaking, there are simultane-ously ferro and antiferro components of the orbital ordering.Nevertheless, contrary to the P1 phase, where there is only theferro component, the P1 phase also contains the antiferro one.Employing linear response theory [68], we study the localstability of the FM state. The spin-wave stiffness tensors,D̂, calculated for each of the orbital states of the R3, P1,and P1 symmetry, and the corresponding spin-wave disper-sions near the  point of the Brillouin zone are displayedin Fig. 7. For the R3 symmetry, the tensor D̂ is negativedefinite in the xy plane, meaning that the FM state is unstable.With lowering the symmetry R3 → P1 → P1, the tensor D̂becomes positive-definite. Thus, the FM order is stabilizedby the orbital order in the P1 and P1 phases. Furthermore,lowering the symmetry results in the anisotropy of D̂ in thexy plane, which increases in the direction P1 → P1. TheFIG. 8. Orbital ordering of the P1 symmetry around two V sitesin the rhombohedral cell.AFM alignment of two V sublattices was also considered.However, the obtained AFM phases were substantially higherin energy than the FM ones. The details can be found in theSupplemental Material [41].We also tried to choose various starting conditions in theHF calculations, assuming different populations of the atomicstates in the initial guess for the potential (2) and then solvingthe problem self-consistently. Nevertheless, such calculationssteadily converged to one of the orbital ordering displayed inFig. 6. Thus, the degeneracy of the orbitally ordered states,which was encountered, for instance, in the toy-model analy-sis in Sec. II B, is lifted, even without lattice distortions. Thereason is that the degeneracy of unoccupied eσg states is liftedby electron-electron interactions with the occupied configura-tion eπg a1g. Therefore, the virtual hoppings into the subspaceof unoccupied eσg orbitals, relevant to the SE process, are nolonger equivalent, leading to a specific type of the orbitalordering. In the R3 crystal structure, the only equivalent typesof the orbital order can be obtained by rotating the ones inFig. 6 by ±120◦ about z. Then, the orbital ordering pattern ofthe R3 symmetry will transform to itself, while for each of theP1 and P1 symmetries, such rotations will generate two moreequivalent domains. Then, applying the inversion operationto such three orbital domains of the P1 symmetry, one cangenerate three more domains.In the P1 state, the orbital ordering breaks the inver-sion symmetry and induces the electric polarization P. Thelatter can be evaluated using Berry-phase theory [38–40],which can be adapted for the model Hamiltonian (1) in theHF approximation [69]. Without SO interaction, it yieldsP = (−0.06, 0.18,−0.05) µC/cm2. Thus, P appears to befinite in the honeycomb plane as well as in the directionperpendicular to the plane. In the SE approximation, whichcan be derived starting from the general Berry-phase the-ory [37,69], this polarization takes the pairwise form P =∑〈i j〉 Pi j and can be expressed via weights of the Wannierfunctions, transferred to the neighboring sites, |αi→ j |2:Pi j = eτ i jV(|α j→i|2 − |αi→ j |2), (3)where τ i j = Ri − R j is the vector connecting the atomic sitej with the site i. The crystal structure of VI3 is such thataround each site V1 there are three neighboring sites V2in the honeycomb plane (2‖1 - 2‖3 in Fig. 8, separated by3.95 Å) and one next-nearest-neighbor site V2 in the205116-8FERROMAGNETIC FERROELECTRICITY DUE TO THE … PHYSICAL REVIEW B 110, 205116 (2024)perpendicular direction (2⊥, separated by 6.55 Å). In theR3 structure, the sublattices V1 and V2 are transformed toeach other by the spatial inversion. Therefore, such bondsare centrosymmetric and |αi→ j |2 = |α j→i|2. Nevertheless, theorbital order breaks the inversion symmetry, so that besidesthe symmetric part, each of |αi→ j |2 will also acquire the anti-symmetric one, resulting in finite Pi j . Similar situation occursaround site 2. If for each of the bond 1 j around site 1, 2 j′ is theequivalent to it bond around site 2 (say, 12‖1 and 21‖1, or 12⊥and 21⊥ in Fig. 8), we have |α1→ j |2 = |α j′→2|2, τ1 j = −τ2 j′ ,and therefore P1 j = P2 j′ , resulting in final total P. Thus, thebonds 12⊥ and 21⊥ are responsible for finite Pz, while otherbonds are responsible for finite polarization in the honeycombplane.E. Triclinic distortionIn this section we briefly consider the effect of triclinicdistortion on the orbital ordering, using for these purposesexperimental parameters of the P1 structure at T = 9 K [46].The triclinic distortion additionally splits the eπg levels. How-ever, the magnitude of this splitting (∼10 – 20 meV [41]) issmaller than B, so that the main tendencies obtained for thetrigonal structure remain largely intact.In the HF approximation, we were able to obtain two dis-tinct solutions for the FM state (Fig. 9). The first one hasthe P1 symmetry, the same as the crystal structure, where twoV sublattices are transformed to each other by the spatial in-version. The second solution has the P1 symmetry, where theinversion symmetry is broken by the antiferro orbital order.The solutions are nearly degenerate. The spin-wave stiffnessFIG. 9. Top: Two types of the orbital ordering of the P1 andP1 symmetry, obtained in the triclinic structure of VI3 for the fer-romagnetic state. The crystallographic inversion center is denotedby +. �E is the energy change relative to the orbital state of theP1 symmetry. Bottom: Spin-wave stiffness tensors (in meVÅ2) forthe orbital states of the P1 and P1 symmetry, and correspondingspin-wave dispersions near the  point of the Brillouin zone.tensor D̂ is positive-definite for both symmetries, meaning thatthe FM state is stable. Enforcing B = 0 in the HF calculations,the P1 solution steadily relaxes to the P1 one. Thus, the Racahparameter B is solely responsible for the antiferro orbital orderand breaking the inversion symmetry.F. Spin-orbit interaction and magnetic-field controlof electric polarizationSimilar to more studied CrI3 [70], the main source ofthe SO interaction in VI3 are the heavy iodine atoms. Thecorresponding parameter of the SO coupling, ξ ∼ 0.1 eV, iscomparable to B and larger than �tr . Moreover, unlike inCrI3, the majority-spin t2g shell in VI3 is only partially filled.Therefore, one can generally expect the SO coupling to play avery important role in VI3.The SO interaction practically does not change the dis-tribution of electron density around V sites [see Fig. 10(a)].However, it has a profound effect on other elements of the den-sity matrix, responsible for magnetic properties. Furthermore,besides spin, there is an appreciable orbital magnetization.Breaking the threefold rotation symmetry in the orbitally or-dered states of the P1 and P1 symmetry will lead to a cantedmagnetic structure. On the microscopic level, this canting isrelated to the population of the a1g orbital and one of the eπgorbitals, so that the SO interaction operating between the oc-cupied a1g and unoccupied eπg orbitals will obey the selectionrules: �σ = ±1 and �m = ∓1, leading to rotation of spin(MνS) and orbital (MνL) magnetic moments away from the zaxis. In the P1 state, where two V sites remain equivalent, thisrotation occurs in the same direction for ν = V1 and V2. Thus,the canting is ferromagnetic. The additional breaking of inver-sion symmetry in the P1 state will lead to the AFM canting ofspin and orbital magnetic moments, which can be regardedas an effect of Dzyaloshinskii-Moriya interaction [71,72],induced by the antiferro orbital ordering in the otherwisecentrosymmetric crystal structure. Microscopically, the AFMcanting occurs because, in the P1 state, the unoccupied eπgFIG. 10. Results of Hartree-Fock calculations with the spin-orbitinteraction: (a) Orbital ordering of the P1 symmetry, (b) Spin mag-netic moments (MνS), (c) Orbital magnetic moments (MνL). Thenumerical values of MνS and MνL in two magnetic sublattices are givenin parentheses. The crystallographic inversion center is denotedby +.205116-9I. V. SOLOVYEV, R. ONO, AND S. A. NIKOLAEV PHYSICAL REVIEW B 110, 205116 (2024)FIG. 11. Magnetic-field dependence of magnetization (top) andelectric polarization (bottom). The field is applied parallel to the zaxis. MzS and MzS + MzL are the z components of, respectively, spinand total magnetic moments. Pz is the z component of the electricpolarization.orbital at the sites 1 and 2 are different, resulting in the differ-ent coupling with the occupied a1g orbitals. The FM cantingtakes place mainly in the yz plane [see Figs. 10(b) and 10(c)].First, we note that besides the spin, there is a large orbitalmagnetic moment MzL = 0.5 µB, which is consistent with theexperimental value of 0.6 µB derived from the x-ray magneticcircular dichroism spectra [73]. The polar angles formed bythe spin, orbital, and total (MνS + MνL) magnetic momentswith the z axis can be estimated as ϑS = 24◦, ϑL = 160◦,and ϑS+L = 25◦, respectively. The latter is consistent with theexperimental estimate of ϑexpS+L ∼ 36◦ [49]. Then, the AFMcanting takes place mainly in the xy plane, where the magneticmoments of the sublattices V1 and V2 are additionally rotatedrelative to each other by �ϕS = 5◦ and �ϕL = 36◦, for thespin and orbital counterparts, respectively.Thus, the SO interaction largely modifies the magneticstructure of VI3. This changes the electric polarization dra-matically. Using Berry-phase theory, one can readily evaluatethe electronic part of P, associated with the change of theelectronic structure after taking into account the SO coupling.It yields P = (−1.86, 1.59,−3.85) µC/cm2, which exceedsthe same value obtained without SO coupling by more thanone order of magnitude. This change is associated with the ad-ditional contributions to the magnetoelectric coupling, whichare activated by the SO interaction. First, the orbital magne-tization can additionally contribute to P [74,75], which is aquite plausible scenario in the present case because ML islarge. Another contribution to P is due to the noncollinearmagnetic alignment [76–78].Then, we consider how the electric polarization can becontrolled by external magnetic field H . The basic idea is thatby applying the magnetic field one can control the cantingof magnetic moments and the degree of mixing of the a1gand eπg characters in the ground state, which plays a crucialrole in establishing the antiferro orbital ordering and devel-oping the electric polarization. Thus, P can be eventuallycontrolled by H . Then, we apply the magnetic field parallelto z, H = (0, 0, H ), and monitor the behavior of spin andorbital magnetic moments as well as the electric polarization,derived from the HF calculations. The results are summarizedin Fig. 11: the magnetization reveals the specific hysteresisloop, while the polarization has a butterflylike shape. Themagnetic field H ∼ 10 T applied in the direction of MzS issufficient to saturate both magnetization and the electric polar-ization (Pz ∼ 3 µC/cm2). The magnetic field in the oppositedirection gradually decreases MzS and MzS + MzL and increasesthe xy components of these moments. Then, H ∼ −7 T causesthe reorientation of MzS along the field. The correspondingpolarization undergoes the jump �Pz ∼ 2.4 µC/cm2, whichis comparable to �Pz ∼ 1.7 µC/cm2 in CaBaCo4O7 and sofar regarded as the largest experimentally observed change ofelectric polarization induced by the magnetic field [79].IV. SUMMARY AND OUTLOOKWe have proposed a route for designing ferroelectric ferro-magnets: a fundamentally and practically important subclassof multiferroic materials, which are not simply magnetic,but ferromagnetic. Our basic idea is that the antiferro orbitalordering across the inversion center should not only producethe FM interactions between the spins, as it follows from theGKA rules, but can also break the inversion symmetry. Thevitality of this idea was illustrated on the toy models of theorbital ordering in the zigzag chain and honeycomb plane.Then, we have proposed that such a scenario can be indeedrealized in the van der Waals ferromagnet VI3, where theHund’s second rule effects tend to form the atomic groundstate with the greatest possible multiplicity, thus unquenchingthe orbital degrees of freedom and activating the KK mech-anism of the orbital ordering. This mechanism is responsiblefor the antiferro orbital order in VI3, which breaks not only thethreefold rotation, but also inversion symmetry, resulting inthe FE-FM ground state. Thus, the orbital degeneracy is lifted,as it is required by the Jahn-Teller theorem [30]. However,this degeneracy lifting occurs via the SE processes, whereasthe crystal distortions probably play a secondary role. Thisis in line with general symmetry considerations suggestingthat the pseudo-Jahn-Teller mechanism, which results in thenoncentrosymmetric ionic displacements, is not operative inthe d2 systems, such as VI3 [12].The relativistic SO interaction, collaborating with the sym-metry breaking, results in the canting of magnetization, whichcan be further manipulated by the magnetic field. This opensa possibility for controlling the electric polarization, whichundergoes a huge change in the magnetic field.The available experimental information about the crystalstructure of VI3, especially regarding the stacking of the hon-eycomb planes as well as the symmetry of these planes, isvery controversial, as several different structures have beenproposed for the room-temperature as well as low-temperaturephases [45–50]. We consider such fragility of the crystal205116-10FERROMAGNETIC FERROELECTRICITY DUE TO THE … PHYSICAL REVIEW B 110, 205116 (2024)structure to be a manifestation of the orbital phenomena inVI3: it appears that the orbital degrees of freedom in VI3remain flexible and there may be several scenarios of liftingthe orbital degeneracy depending on the experimental condi-tions. We hope that our scenario, where the orbital orderingnot only lifts the degeneracy but also breaks the inversionsymmetry, leading to the ferroelectric ferromagnetism, can beeventually realized in VI3. From this perspective, particularlyinteresting are the results of Ref. [47], where the P31c andC2/c structures were proposed for, respectively, the high-and low-temperature phases of VI3. These structures are stillcentrosymmetric. However, the vanadium sites V1 and V2,forming the honeycomb planes, become inequivalent. Thus,within each honeycomb plane, the inversion symmetry ap-pears to be broken and this is consistent with our scenarioof the antiferro orbital order. Another, again indirect, indica-tion of the inversion symmetry breaking in the honeycombplane is the different behavior of the magnetic sublattices V1and V2 reported in certain temperature range (36 K < T <51 K) [50].There is a number of theoretical studies reporting thethreefold rotation breaking in VI3 at the level of DFT + Ucalculations with the SO coupling, due to partial popula-tion of the a1g state [73,80–83]. Nevertheless, none of thesestudies reported the inversion symmetry breaking. The sym-metry breaking in Refs. [73,80–83] is solely related to the SOinteraction: the eπg eπg configuration, respecting the threefoldrotation symmetry, would yield only small ML (being about−0.1 µB, along the z axis, emerging due to the mixing of eπgwith the unoccupied eσg states by the SO coupling). Therefore,in order to increase ML and thus maximize the energy gaincaused by the SO interaction, it is essential to rotate ML awayfrom the z axis by breaking the threefold rotation symmetryand populating the a1g state. Another theoretical scenario ofthreefold inversion breaking is based on the electron-latticeinteractions, resulting in three distinct V-V bond lengths inthe honeycomb plane, as was theoretically proposed for theFM VCl3 monolayer without SO coupling [84]. Then, theelectric polarization can be induced by considering a substrateeffect, i.e., pretty much similar to the standard procedureemployed in heterostructures [7]. What we propose here isfundamentally different: according to our scenario, the sym-metry can be broken by ordering the orbitals, which wouldremain degenerate in the atomic limit. Formally, neither SOinteraction nor lattice distortion are needed in our case, thoughthey can play an important role by further facilitating thesymmetry breaking and for establishing the magnetic-fieldcontrol of P. Whether such behavior can be indeed achieved atthe level of DFT + U calculations depends on the implemen-tation, which must include all necessary terms proportionalto the Racah parameter B. Although it was considered onearlier stages, where the DFT + U functional was formulatedin terms of all Slater integrals and, therefore, explicitly in-cluded the dependence on the Racah parameter B [85,86],the latter, commonly used but simplified, versions are for-mulated in terms of only one parameter Ueff = U − J and,thus, disregard the contributions responsible for Hund’s sec-ond rule [87,88].In the most general formulation, DFT should be able toincorporate the exchange-correlation interactions responsiblefor atomic Hund’s second rule. However, such effects areomitted in many popular approximations supplementing DFT,such as LDA or GGA, which take the functional form ofthese interactions from the limit of homogeneous electrongas, where, strictly speaking, the atomic Hund’s rules are nolonger applicable, resulting in a number of fundamental issuesfor atomic systems [61]. Therefore, a very popular directionaround 1990s was to simulate the Hund’s rule physics on thetop of LDA by introducing a phenomenological correction,proportional to some appropriate Racah parameter (B for 3delectrons), with the aim to reproduce the orbital magnetizationin solids, which was severely underestimated in the localspin-density approximation [89–91]. For the Mott insulators,the logic behind was that the Coulomb repulsion U shouldsplit the occupied and unoccupied 3d states [92]. However, theRacah parameter B is another important ingredient to decidethe correct symmetry of states, which will be further split byU [90]. Nevertheless, in most of the cases, such symmetry isalready decided by the crystal field and SO interaction. More-over, the orbital magnetization strongly depends on the valueof Coulomb repulsion U , as it controls the strength of the hy-bridization between occupied and unoccupied states [93,94].In many cases, the value of orbital magnetization can bereproduced by the Coulomb U alone, especially when it istreated as an adjustable parameter. From this perspective, theunique aspect of VI3 is that both parameters, U and B, appearto be important for finding the correct ground state.While the importance of intra-atomic exchange couplingJ in the physics of strongly correlated materials is well rec-ognized today [95], the more delicate effects, driven by theRacah parameter B, remain largely unexplored. It is true thatB is much smaller than J (typically, B ∼ 0.1J), reflectingthe well-known hierarchy of atomic Hund’s rules, when thesecond rule always follows the first one. Nevertheless, if Bis larger or comparable to the characteristic crystal field, theHund’s second rule effects can lead to a number of interestingand so far unexplored effects. The ferromagnetic ferroelec-tricity in VI3, which we propose in this work, is certainly oneof them.ACKNOWLEDGMENTSMANA is supported by World Premier International Re-search Center Initiative (WPI), MEXT, Japan. R.O. wassupported by JSPS KAKENHI Grant No. JP23KJ2165. Thecomputations in this study were partly performed on the Nu-merical Materials Simulator at the NIMS.We are deeply indebted to the late Prof. Daniel Khomskiifor his profound contributions to this field.I.V.S. conceptualized the work, performed most of the cal-culations (except specified below), and wrote the manuscript.R.O. and S.A.N. performed the electronic structure calcula-tions and constructed the model Hamiltonian using the QEmethod. S.A.N. performed the DMFT calculations. All au-thors discussed the results and commented on the manuscript.The authors declare no competing interests.205116-11I. V. SOLOVYEV, R. ONO, AND S. A. NIKOLAEV PHYSICAL REVIEW B 110, 205116 (2024)[1] S.-W. Cheong and M. Mostovoy, Multiferroics: A magnetictwist for ferroelectricity, Nat. Mater. 6, 13 (2007).[2] D. Khomskii, Classifying multiferroics: Mechanisms and ef-fects, Physics 2, 20 (2009).[3] Y. Tokura, S. Seki, and N. Nagaosa, Multiferroics of spin origin,Rep. Prog. Phys. 77, 076501 (2014).[4] N. Hill, Why are there so few magnetic ferroelectrics? J. Phys.Chem. B 104, 6694 (2000).[5] W. Eerenstein, N. D. Mathur, and J. F. Scott, Multiferroic andmagnetoelectric materials, Nature (London) 442, 759 (2006).[6] Y. Tokura, Multiferroics as quantum electromagnets, Science312, 1481 (2006).[7] C. A. Fernandes Vaz and U. Staub, Artificial multiferroic het-erostructures, J. Mater. Chem. C 1, 6731 (2013).[8] C. J. Fennie and K. M. Rabe, Magnetic and electric phasecontrol in epitaxial EuTiO3 from first principles, Phys. Rev.Lett. 97, 267602 (2006).[9] J. H. Lee et al., A strong ferroelectric ferromagnet createdby means of spin-lattice coupling, Nature (London) 466, 954(2010).[10] J. H. Lee and K. M. Rabe, Epitaxial-strain-induced multifer-roicity in SrMnO3 from first principles, Phys. Rev. Lett. 104,207204 (2010).[11] I. B. Bersuker, On the origin of ferroelectricity in perovskite-type crystals, Phys. Lett. 20, 589 (1966).[12] I. B. Bersuker, Pseudo Jahn-Teller origin of perovskite multi-ferroics, magnetic-ferroelectric crossover, and magnetoelectriceffects: The d0-d10 problem, Phys. Rev. Lett. 108, 137202(2012).[13] A. Edström and C. Ederer, First-principles-based strain andtemperature-dependent ferroic phase diagram of SrMnO3, Phys.Rev. Mater. 2, 104409 (2018).[14] K. I. Kugel’ and D. I. Khomskii, The Jahn-Teller effect andmagnetism: Transition metal compounds, Sov. Phys. Usp. 25,231 (1982).[15] P. W. Anderson, Antiferromagnetism. Theory of superexchangeinteraction, Phys. Rev. 79, 350 (1950).[16] J. B. Goodenough, Theory of the role of covalence in theperovskite-type manganites [La, M(II)]MnO3, Phys. Rev. 100,564 (1955).[17] J. B. Goodenough, An interpretation of the magnetic propertiesof the perovskite-type mixed crystals La1−xSrxCoO3−λ, J. Phys.Chem. Solids 6, 287 (1958).[18] J. Kanamori, Superexchange interaction and symmetry prop-erties of electron orbitals, J. Phys. Chem. Solids 10, 87(1959).[19] K. I. Kugel and D. I. Khomskii, Superechange ordering ofdegenerate orbitals and magnetic structure of dielectrics withJahn-Teller ions, Zh. Eksp. Teor. Fiz. 15, 629 (1972) [JETP Lett.15, 446 (1972)].[20] D. I. Khomskii and K. I. Kugel, Orbital and magnetic structureof two-dimensional ferromagnets with Jahn-Teller ions, SolidState Commun. 13, 763 (1973).[21] P. W. Anderson, New approach to the theory of superexchangeinteractions, Phys. Rev. 115, 2 (1959).[22] D. I. Khomskii and S. V. Streltsov, Orbital effects in solids:Basics, recent progress, and opportunities, Chem. Rev. 121,2992 (2021).[23] R. Maezono, S. Ishihara, and N. Nagaosa, Phase diagram ofmanganese oxides, Phys. Rev. B 58, 11583 (1998).[24] M. Mochizuki and M. Imada, Orbital-spin structure and latticecoupling in RTiO3 where R=La, Pr, Nd, and Sm, Phys. Rev.Lett. 91, 167203 (2003).[25] A. Moreira dos Santos, A. K. Cheetham, T. Atou, Y. Syono, Y.Yamaguchi, K. Ohoyama, H. Chiba, and C. N. R. Rao, Orbitalordering as the determinant for ferromagnetism in biferroicBiMnO3, Phys. Rev. B 66, 064425 (2002).[26] I. V. Solovyev and Z. V. Pchelkina, Orbital ordering and mag-netic interactions in BiMnO3, New J. Phys. 10, 073021 (2008).[27] P. Baettig, R. Seshadri, and N. A. Spaldin, Anti-polarity in idealBiMnO3, J. Am. Chem. Soc. 129, 9854 (2007).[28] K. Yamauchi and P. Barone, Electronic ferroelectricity inducedby charge and orbital orderings, J. Phys.: Condens. Matter 26,103201 (2014).[29] T. Kong, K. Stolze, E. I. Timmons, J. Tao, D. Ni, S. Guo,Z. Yang, R. Prozorov, and R. J. Cava, VI3 – a new lay-ered ferromagnetic semiconductor, Adv. Mater. 31, 1808074(2019).[30] H. A. Jahn and E. Teller, Stability of polyatomic moleculesin degenerate electronic states - I–Orbital degeneracy, Proc. R.Soc. A 161, 220 (1937).[31] G. Racah, Theory of complex spectra. II, Phys. Rev. 62, 438(1942).[32] J. C. Slater, Quantum Theory of Atomic Structure (McGraw-Hill, New York, 1960).[33] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Dynamical mean-field theory of strongly correlated fermionsystems and the limit of infinite dimensions, Rev. Mod. Phys.68, 13 (1996).[34] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O.Parcollet, and C. A. Marianetti, Electronic structure calculationswith dynamical mean-field theory, Rev. Mod. Phys. 78, 865(2006).[35] J. C. Slater and G. F. Koster, Simplified LCAO method for theperiodic potential problem, Phys. Rev. 94, 1498 (1954).[36] I. V. Solovyev and S. A. Nikolaev, Double-exchange theoryof ferroelectric polarization in orthorhombic manganites withtwofold periodic magnetic texture, Phys. Rev. B 87, 144424(2013).[37] R. Ono, S. Nikolaev, and I. Solovyev, Fingerprints of spin-current physics on magnetoelectric response in the spin- 12magnet Ba2CuGe2O7, Phys. Rev. B 102, 064422 (2020).[38] R. D. King-Smith and D. Vanderbilt, Theory of polarization ofcrystalline solids, Phys. Rev. B 47, 1651(R) (1993).[39] D. Vanderbilt and R. D. King-Smith, Electric polarization as abulk quantity and its relation to surface charge, Phys. Rev. B 48,4442 (1993).[40] R. Resta, Electrical polarization and orbital magnetization:the modern theories, J. Phys.: Condens. Matter 22, 123201(2010).[41] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevB.110.205116 for details of the toy-modelanalysis, model parameters for VI3, and further details ofDMFT and HF calculations, which includes Ref. [42].[42] O. Arcelus, S. Nikolaev, J. Carrasco, and I. Solovyev, Mag-netism of NaFePO4 and related polyanionic compounds, Phys.Chem. Chem. Phys. 20, 13497 (2018); M. Wallerberger et al.,w2dynamics: Local one- and two-particle quantities from dy-namical mean field theory, Comput. Phys. Commun. 235, 388(2019); R. Levy, J. P. F. LeBlanc, and E. Gull, Implementation205116-12https://doi.org/10.1038/nmat1804https://doi.org/10.1103/Physics.2.20https://doi.org/10.1088/0034-4885/77/7/076501https://doi.org/10.1021/jp000114xhttps://doi.org/10.1038/nature05023https://doi.org/10.1126/science.1125227https://doi.org/10.1039/c3tc31428fhttps://doi.org/10.1103/PhysRevLett.97.267602https://doi.org/10.1038/nature09331https://doi.org/10.1103/PhysRevLett.104.207204https://doi.org/10.1016/0031-9163(66)91127-9https://doi.org/10.1103/PhysRevLett.108.137202https://doi.org/10.1103/PhysRevMaterials.2.104409https://doi.org/10.1070/PU1982v025n04ABEH004537https://doi.org/10.1103/PhysRev.79.350https://doi.org/10.1103/PhysRev.100.564https://doi.org/10.1016/0022-3697(58)90107-0https://doi.org/10.1016/0022-3697(59)90061-7http://jetpletters.ru/ps/1754/article_26676.pdfhttps://doi.org/10.1016/0038-1098(73)90362-1https://doi.org/10.1103/PhysRev.115.2https://doi.org/10.1021/acs.chemrev.0c00579https://doi.org/10.1103/PhysRevB.58.11583https://doi.org/10.1103/PhysRevLett.91.167203https://doi.org/10.1103/PhysRevB.66.064425https://doi.org/10.1088/1367-2630/10/7/073021https://doi.org/10.1021/ja073415uhttps://doi.org/10.1088/0953-8984/26/10/103201https://doi.org/10.1002/adma.201808074https://doi.org/10.1098/rspa.1937.0142https://doi.org/10.1103/PhysRev.62.438https://doi.org/10.1103/RevModPhys.68.13https://doi.org/10.1103/RevModPhys.78.865https://doi.org/10.1103/PhysRev.94.1498https://doi.org/10.1103/PhysRevB.87.144424https://doi.org/10.1103/PhysRevB.102.064422https://doi.org/10.1103/PhysRevB.47.1651https://doi.org/10.1103/PhysRevB.48.4442https://doi.org/10.1088/0953-8984/22/12/123201http://link.aps.org/supplemental/10.1103/PhysRevB.110.205116https://doi.org/10.1039/C8CP01961Dhttps://doi.org/10.1016/j.cpc.2018.09.007FERROMAGNETIC FERROELECTRICITY DUE TO THE … PHYSICAL REVIEW B 110, 205116 (2024)of the maximum entropy method for analytic continuation, ibid.215, 149 (2017); J. Kanamori, Theory of the magnetic proper-ties of ferrous and cobaltous oxides, I, Prog. Theor. Phys. 17,177 (1957).[43] P. Barone, K. Yamauchi, and S. Picozzi, Ferroelectricity due toorbital ordering in E -type undoped rare-earth manganites, Phys.Rev. Lett. 106, 077201 (2011).[44] G. Natta, P. Corradini, and G. Allegra, The different crystallinemodifications of TiCl3, a catalyst component for the polymer-ization of α-olefins. I: α-, β-, γ -TiCl3. II: δ-TiCl3, J. Polym.Sci. 51, 399 (1961).[45] P. Doležal, M. Kratochvílová, V. Holý, P. Čermák, V.Sechovský, M. Dušek, M. Míšek, T. Chakraborty, Y. Noda, S.Son, and J.-G. Park, Crystal structures and phase transitionsof the van der Waals ferromagnet VI3, Phys. Rev. Mater. 3,121401(R) (2019).[46] Th. Marchandier, N. Dubouis, F. Fauth, M. Avdeev, A.Grimaud, J.-M. Tarascon, and G. Rousse, Crystallographic andmagnetic structures of the VI3 and LiVI3 van der Waals com-pounds, Phys. Rev. B 104, 014105 (2021).[47] S. Son, M. J. Coak, N. Lee, J. Kim, T. Y. Kim, H. Hamidov,Hwanbeom Cho, C. Liu, D. M. Jarvis, P. A. C. Brown, J. H.Kim, C.-H. Park, D. I. Khomskii, S. S. Saxena, and J.-G. Park,Bulk properties of the van der Waals hard ferromagnet VI3,Phys. Rev. B 99, 041402(R) (2019).[48] S. Tian, J.-F. Zhang, C. Li, T. Ying, S. Li, X. Zhang, K. Liu, andH. Lei, Ferromagnetic van der Waals Crystal VI3, J. Am. Chem.Soc. 141, 5326 (2019).[49] Y. Hao, Y. Gu, Y. Gu, E. Feng, H. Cao, S. Chi, H. Wu, andJ. Zhao, Magnetic order and its interplay with structure phasetransition in van der Waals ferromagnet VI3, Chin. Phys. Lett.38, 096101 (2021).[50] E. Gati, Y. Inagaki, T. Kong, R. J. Cava, Y. Furukawa, P. C.Canfield, and S. L. Bud’ko, Multiple ferromagnetic transitionsand structural distortion in the van der Waals ferromagnet VI3at ambient and finite pressures, Phys. Rev. B 100, 094408(2019).[51] O. K. Andersen, Linear methods in band theory, Phys. Rev. B12, 3060 (1975).[52] O. Gunnarsson, O. Jepsen, and O. K. Andersen, Self-consistentimpurity calculations in the atomic-spheres approximation,Phys. Rev. B 27, 7144 (1983).[53] P. Giannozzi, S. Baroni, N. Bonini et al., QUANTUMESPRESSO: A modular and open-source software project forquantum simulations of materials, J. Phys.: Condens. Matter 21,395502 (2009).[54] S. H. Vosko, L. Wilk, and M. Nusair, Accurate spin-dependentelectron liquid correlation energies for local spin densitycalculations: A critical analysis, Can. J. Phys. 58, 1200(1980).[55] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradientapproximation made simple, Phys. Rev. Lett. 77, 3865 (1996);78, 1396 (1997).[56] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D.Vanderbilt, Maximally localized Wannier functions: Theory andapplications, Rev. Mod. Phys. 84, 1419 (2012).[57] A. A. Mostofi, J. R. Yates, G. Pizzi, Y. S. Lee, I. Souza,D. Vanderbilt, and N. Marzari, wannier90: A tool for ob-taining maximally-localised Wannier functions, Comput. Phys.Commun. 185, 2309 (2014).[58] I. V. Solovyev, Combining DFT and many-body methods tounderstand correlated materials, J. Phys.: Condens. Matter 20,293201 (2008).[59] F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S.Biermann, and A. I. Lichtenstein, Frequency-dependent localinteractions and low-energy effective models from electronicstructure calculations, Phys. Rev. B 70, 195104 (2004).[60] S. Yu. Ezhov, V. I. Anisimov, D. I. Khomskii, and G. A.Sawatzky, Orbital occupation, local spin, and exchange inter-actions in V2O3, Phys. Rev. Lett. 83, 4136 (1999).[61] M. Weinert, R. E. Watson, and G. W. Fernando, Density-functional theory and atomic multiplet levels, Phys. Rev. A 66,032508 (2002).[62] Y. Tanabe and S. Sugano, On the absorption spectra of complexions II, J. Phys. Soc. Jpn. 9, 766 (1954).[63] C. Castellani, C. R. Natoli, and J. Ranninger, Magnetic struc-ture of V2O3 in the insulating phase, Phys. Rev. B 18, 4945(1978).[64] P. Rozier, A. Ratusznab, and J. Galy, Comparative Structuraland Electrical Studies of V2O3 and V2−xNixO3 (0 < x < 0.75)Solid Solution, Z. Anorg. Allg. Chem. 628, 1236 (2002).[65] G. Khaliullin, P. Horsch, and A. M. Oleś, Spin order due toorbital fluctuations: Cubic vanadates, Phys. Rev. Lett. 86, 3879(2001).[66] P. Horsch, G. Khaliullin, and A. M. Oleś, Dimerization versusorbital-moment ordering in a mott insulator YVO3, Phys. Rev.Lett. 91, 257203 (2003).[67] C. Ulrich, G. Khaliullin, J. Sirker, M. Reehuis, M. Ohl, S.Miyasaka, Y. Tokura, and B. Keimer, Magnetic neutron scatter-ing study of YVO3: Evidence for an orbital peierls state, Phys.Rev. Lett. 91, 257202 (2003).[68] I. V. Solovyev, Linear response theories for interatomic ex-change interactions, J. Phys.: Condens. Matter 36, 223001(2024).[69] S. A. Nikolaev and I. V. Solovyev, Microscopic theory of elec-tric polarization induced by skyrmionic order in GaV4S8, Phys.Rev. B 99, 100401(R) (2019).[70] J. L. Lado and J. Fernández-Rossier, On the origin of mag-netic anisotropy in two dimensional CrI3, 2D Mater. 4, 035002(2017).[71] I. Dzyaloshinsky, A thermodynamic theory of “weak” ferro-magnetism of antiferromagnetics, J. Phys. Chem. Solids 4, 241(1958).[72] T. Moriya, Anisotropic superexchange interaction and weakferromagnetism, Phys. Rev. 120, 91 (1960).[73] D. Hovančík, J. Pospíšil, K. Carva, V. Sechovský, and C.Piamonteze, Large orbital magnetic moment in VI3, Nano Lett.23, 1175 (2023).[74] A. Malashevich, S. Coh, I. Souza, and D. Vanderbilt, Full mag-netoelectric response of Cr2O3 from first principles, Phys. Rev.B 86, 094430 (2012).[75] A. Scaramucci, E. Bousquet, M. Fechner, M. Mostovoy, andN. A. Spaldin, Linear magnetoelectric effect by orbital mag-netism, Phys. Rev. Lett. 109, 197203 (2012).[76] H. Katsura, N. Nagaosa, and A. V. Balatsky, Spin current andmagnetoelectric effect in noncollinear magnets, Phys. Rev. Lett.95, 057205 (2005).[77] I. V. Solovyev, Superexchange theory of electronic polarizationdriven by relativistic spin-orbit interaction at half filling, Phys.Rev. B 95, 214406 (2017).205116-13https://doi.org/10.1016/j.cpc.2017.01.018https://doi.org/10.1143/PTP.17.177https://doi.org/10.1103/PhysRevLett.106.077201https://doi.org/10.1002/pol.1961.1205115602https://doi.org/10.1103/PhysRevMaterials.3.121401https://doi.org/10.1103/PhysRevB.104.014105https://doi.org/10.1103/PhysRevB.99.041402https://doi.org/10.1021/jacs.8b13584https://doi.org/10.1088/0256-307X/38/9/096101https://doi.org/10.1103/PhysRevB.100.094408https://doi.org/10.1103/PhysRevB.12.3060https://doi.org/10.1103/PhysRevB.27.7144https://doi.org/10.1088/0953-8984/21/39/395502https://doi.org/10.1139/p80-159https://doi.org/10.1103/PhysRevLett.77.3865https://doi.org/10.1103/PhysRevLett.78.1396https://doi.org/10.1103/RevModPhys.84.1419https://doi.org/10.1016/j.cpc.2014.05.003https://doi.org/10.1088/0953-8984/20/29/293201https://doi.org/10.1103/PhysRevB.70.195104https://doi.org/10.1103/PhysRevLett.83.4136https://doi.org/10.1103/PhysRevA.66.032508https://doi.org/10.1143/JPSJ.9.766https://doi.org/10.1103/PhysRevB.18.4945https://doi.org/10.1002/1521-3749(200206)628:5<1236::AID-ZAAC1236>3.0.CO;2-Chttps://doi.org/10.1103/PhysRevLett.86.3879https://doi.org/10.1103/PhysRevLett.91.257203https://doi.org/10.1103/PhysRevLett.91.257202https://doi.org/10.1088/1361-648X/ad215ahttps://doi.org/10.1103/PhysRevB.99.100401https://doi.org/10.1088/2053-1583/aa75edhttps://doi.org/10.1016/0022-3697(58)90076-3https://doi.org/10.1103/PhysRev.120.91https://doi.org/10.1021/acs.nanolett.2c04045https://doi.org/10.1103/PhysRevB.86.094430https://doi.org/10.1103/PhysRevLett.109.197203https://doi.org/10.1103/PhysRevLett.95.057205https://doi.org/10.1103/PhysRevB.95.214406I. V. SOLOVYEV, R. ONO, AND S. A. NIKOLAEV PHYSICAL REVIEW B 110, 205116 (2024)[78] I. Solovyev, R. Ono, and S. Nikolaev, Magnetically inducedpolarization in centrosymmetric bonds, Phys. Rev. Lett. 127,187601 (2021).[79] V. Caignaert, A. Maignan, K. Singh, Ch. Simon, V. Pralong,B. Raveau, J. F. Mitchell, H. Zheng, A. Huq, and L. C.Chapon, Gigantic magnetic-field-induced polarization andmagnetoelectric coupling in a ferrimagnetic oxide CaBaCo4O7,Phys. Rev. B 88, 174403 (2013).[80] K. Yang, F. Fan, H. Wang, D. I. Khomskii, and H. Wu,VI3: A two-dimensional Ising ferromagnet, Phys. Rev. B 101,100402(R) (2020).[81] L. M. Sandratskii and K. Carva, Interplay of spin magnetism,orbital magnetism, and atomic structure in layered van derWaals ferromagnet VI3, Phys. Rev. B 103, 214451 (2021).[82] T. P. T. Nguyen, K. Yamauchi, T. Oguchi, D. Amoroso, andS. Picozzi, Electric-field tuning of the magnetic properties ofbilayer VI3: A first-principles study, Phys. Rev. B 104, 014414(2021).[83] L. Camerano and G. Profeta, Symmetry breaking in vanadiumtrihalides, 2D Mater. 11, 025027 (2024).[84] D. Guo, C. Wang, L. Wang, Y. Lu, H. Wu, Y. Zhang, and W. Ji,Orbital-ordering driven simultaneous tunability of magnetismand electric polarization in strained monolayer VCl3, Chin.Phys. Lett. 41, 047501 (2024).[85] I. V. Solovyev, P. H. Dederichs, and V. I. Anisimov, Correctedatomic limit in the local-density approximation and the elec-tronic structure of d impurities in Rb, Phys. Rev. B 50, 16861(1994).[86] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Density-functional theory and strong interactions: Orbital ordering inMott-Hubbard insulators, Phys. Rev. B 52, R5467 (1995).[87] I. Solovyev, N. Hamada, and K. Terakura, t2g versus all 3dlocalization in LaMO3 perovskites (M=TiCu): First-principlesstudy, Phys. Rev. B 53, 7158 (1996).[88] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys,and A. P. Sutton, Electron-energy-loss spectra and the structuralstability of nickel oxide: An LSDA+U study, Phys. Rev. B 57,1505 (1998).[89] O. Eriksson, B. Johansson, and M. S. S. Brooks, Meta-magnetism in UCoAl, J. Phys.: Condens. Matter 1, 4005(1989).[90] M. R. Norman, Orbital polarization and the insulating gap inthe transition-metal oxides, Phys. Rev. Lett. 64, 1162 (1990).[91] M. R. Norman, Crystal-field polarization and the insulating gapin FeO, CoO, NiO, and La2CuO4, Phys. Rev. B 44, 1364 (1991).[92] K. Terakura, T. Oguchi, A. R. Williams, and J. Kübler, Bandtheory of insulating transition-metal monoxides: Band-structurecalculations, Phys. Rev. B 30, 4734 (1984).[93] I. V. Solovyev, A. I. Liechtenstein, and K. Terakura, Is Hund’ssecond rule responsible for the orbital magnetism in solids?Phys. Rev. Lett. 80, 5758 (1998).[94] I. V. Solovyev, Orbital polarization in itinerant magnets, Phys.Rev. Lett. 95, 267205 (2005).[95] A. Georges, L. de’ Medici, and J. Mravlje, Strong correlationsfrom hund’s coupling, Annu. Rev. Condens. Matter Phys. 4, 137(2013).205116-14https://doi.org/10.1103/PhysRevLett.127.187601https://doi.org/10.1103/PhysRevB.88.174403https://doi.org/10.1103/PhysRevB.101.100402https://doi.org/10.1103/PhysRevB.103.214451https://doi.org/10.1103/PhysRevB.104.014414https://doi.org/10.1088/2053-1583/ad3137https://doi.org/10.1088/0256-307X/41/4/047501https://doi.org/10.1103/PhysRevB.50.16861https://doi.org/10.1103/PhysRevB.52.R5467https://doi.org/10.1103/PhysRevB.53.7158https://doi.org/10.1103/PhysRevB.57.1505https://doi.org/10.1088/0953-8984/1/25/012https://doi.org/10.1103/PhysRevLett.64.1162https://doi.org/10.1103/PhysRevB.44.1364https://doi.org/10.1103/PhysRevB.30.4734https://doi.org/10.1103/PhysRevLett.80.5758https://doi.org/10.1103/PhysRevLett.95.267205https://doi.org/10.1146/annurev-conmatphys-020911-125045