Orbital Rotations induced by Charges of Polarons and Defects in Doped Vanadates

We explore the competiton of doped holes and defects that leads to the loss of orbital order in vanadate perovskites. In compounds such as La$_{1-{\sf x}}$Ca$_{\,\sf x}$VO$_3$ spin and orbital order result from super-exchange interactions described by an extended three-orbital degenerate Hubbard-Hund model for the vanadium $t_{2g}$ electrons. Long-range Coulomb potentials of charged Ca$^{2+}$ defects and $e$-$e$ interactions control the emergence of defect states inside the Mott gap. The quadrupolar components of the Coulomb fields of doped holes induce anisotropic orbital rotations of degenerate orbitals. These rotations modify the spin-orbital polaron clouds and compete with orbital rotations induced by defects. Both mechanisms lead to a mixing of orbitals, and cause the suppression of the asymmetry of kinetic energy in the $C$-type magnetic phase. We find that the gradual decline of orbital order with doping, a characteristic feature of the vanadates, however, has its origin not predominantly in the charge carriers, but in the off-diagonal couplings of orbital rotations induced by the charges of the doped ions.


I. INTRODUCTION
The discovery that doping holes (or electrons) into Mott insulators (MIs), formed by CuO 2 layers, not just leads to a metallic state and to the decay of the antiferromagnetic (AF) order but also yields high temperature superconductivity [1] was a great surprise. Very early, it was proposed that the mechanism of high temperature superconductivity originates from the strong electron correlations, intrinsic to the MI [2][3][4], rather than from the exchange of phonons. The discovery triggered a systematic study of transition metal oxides [5,6], many of them MIs, also to achieve deeper insights into the fundamental open questions related to doped MIs. More recently, time dependent phenomena came into focus with the discovery of ultracold fermion systems in optical lattices [7][8][9] that represent alternative platforms to study Mott-Hubbard physics.
Spin and/or orbital ordered phases are typical features of orbitally-degenerate MIs. The motion of doped holes in such compounds leads to strings of misplaced spins [10][11][12][13][14][15][16][17][18] or misoriented orbitals [19][20][21] or both [22][23][24][25]. In cuprates, the perturbation induced by doping holes can be efficiently described in terms of spin-polarons [26,27] moving in a twodimensional (2D) spin-1/2 quantum antiferromagnet. Due to the strong quantum fluctuations, their kinetic energy is only weakly reduced such that their binding to charged defects is small. As a consequence, the metal-insulator transition and the breakdown of AF long-range order occur already at quite low doping. In contrast, the perovskite vanadates, that reveal strong quantum orbital fluctuations in certain regimes [28,29], remain insulating up to moderate or even high doping concentrations. For instance, Y 1−x Ca x VO 3 enters a poor metallic state only at x ≈ 0.5 [30]. In the parent compounds RVO 3 , with R ≡ La, Pr, ..., Tb, Y, and Lu, the d 2 configuration of V 3+ ions has an orbital degeneracy of the t 2g electrons and Hund's exchange stabilizes high spin S = 1 states [31].
Orbitally-degenerate Mott insulators display various kinds of spin and orbital order [32][33][34][35][36][37][38][39][40]. RVO 3 has several phases FIG. 1. Orbital rotations due to charges of defects and doped holes in RVO3: (a) CS/GO reference state with C-type order of S = 1 spins of V 3+ (d 2 ) ions and alternating G-type orbital order of occupied a ≡ yz and b ≡ xz vanadium orbitals. The occupied c ≡ xy orbitals (not shown) are stabilized by a crystal field ∆c, see (b). Orbital rotations induced by a charged defect (red dot, Ca 2+ ion) and by the charge of a doped hole (blue dot, V 4+ ion) on nearest neighbor V 3+ ions are shown in (c) and (e), respectively. For clarity we display only the highest rotated (i.e., unoccupied) orbital per V 3+ ion in the strong coupling limit (i.e., ∆c D or D h ), whereas (d) and (f) display the level splittings for moderate values of coupling constants, D and D h . with the respective transition temperatures dependent on the radius of the R-ions. For instance, all compounds from La to Tb have a ground state with C-type spin and G-type orbital order (CS/GO) [41], i.e., with antiferro order in the ab-planes and ferro order along the c axis, respectively, while the occupation of the topmost occupied t 2g orbital alternates between a = yz and b = xz [42] as shown in Fig. 1(a). The same state was found recently in LaVO 3 thin films [43]. Above the Néel temperature T N , there is a paramagnetic phase with G-type orbital order, which disappears above the orbital order temperature T OO . In YVO 3 , and all systems with smaller R ionic radius, there is a further transition at T CG from the CS/GO to a complementary GS/CO low-T state [44][45][46]. The CS/GO order emerges from the intrinsic super-exchange interactions, i.e., driven by strong pseudospin 1/2 orbital quantum fluctuations along the c axis. These orbital fluctuations lead to much stronger FM spin-couplings in the CS/GO phase [47][48][49] than expected for frozen orbitals, i.e., as assumed in Goodenough-Kanamori rules [50]. YVO 3 instead has a GS/CO ground state triggered by also present Jahn-Teller interactions that increase with decreasing R-ion radius [51][52][53][54]. Interestingly, already traces of Ca-doping switch the GS/CO ground state of YVO 3 to the CS/GO-phase [41,55,56], a feature that could be explained in a model for charged defects that we adopt here [57,58].
Here, we explore the stability of the disordered CS/GOphase and its gradual decay at large doping. It is the coupling to the extra orbital degree of freedom in the vanadates that leads to the quenching of the kinetic energy and to strong localization and binding of polarons by the Coulomb potential of defects. Yet, this strong localization creates a new puzzle: how is then the orbital order destroyed in these compounds? As alternative mechanism, the orbital rotations (ORs) at vanadium ions were identified. They are induced by the Coulomb fields of the charged defects [57,58]. It was shown that ORs are an effective perturbation as each defect is surrounded by eight nearest vanadium neighbors, see Fig. 1(b). It yields a natural explanation for the gradual suppression of G-type orbital order in vanadates as function of doping [59], and the absence of clear signatures of collective phase transitions [60]. In this paper, we explore a complementary, a priori equally important, orbital polarization mechanism triggered by the polaron charge, see Fig. 1(c). A mechanism of this kind was found essential by Kilian and Khaliullin [61] in a study of orbital polarons in the orbital liquid regime of e g orbitals in manganites [62][63][64] . For the t 2g orbitals of vanadates we find an abrupt reduction of orbital order caused by doped holes/polarons beyond a critical coupling strength. Yet in combination with OR's induced by the defect charges the gradual decline of orbital order dominates and is amplified by OR's due to the polarons.
Despite the cubic structure, the undoped CS/GO state is highly anisotropic, due to the FM correlations along the caxis, where strong quantum orbital fluctuations boost the super-exchange [28,[47][48][49]. Tokura and coworkers [65] observed that the anisotropy ratio A opt of optical weights along z and x axis changed from about two to one at large doping in the CS phase. A goal of our work is to shed light on this puzzle by studying the asymmetry in the kinetic energy A = K z /K x that is strictly related to A opt [48]. The case of the vanadates is puzzling as the isotropy of kinetic energy is observed at doping concentrations where the anisotropic magnetic CS order still persists.
The article continues in Section II with a brief description of the minimal model for the doped vanadate Mott insulators. The main focus here is on the orbital rotation terms that control the orientation of vanadium orbitals, and are a consequence of the Coulomb fields of defects and doped holes or electrons. In Section III the effect of orbital rotations on the occupation of orbitals, the magnetic and the orbital order is studied. In Section IV we present our conclusions. An Appendix contains further details of the multi-orbital Hubbard-Hund interaction, the Jahn-Teller and other small terms, as well as the derivation of the orbital polarization terms.

II. THE MULTI-ORBITAL MODEL FOR DOPED VANADATE MOTT INSULATORS
The minimal Hamiltonian that describes the t 2g electrons, the Mott gap, and the defect states in R 1−x Ca x VO 3 is [58]: It includes an extended 3-band Hubbard model H Hub [66,67] that describes the electronic multiplet structure of the V 3+ ions [57] and the different phases of the parent compounds. For a first orientation the details of this term can be ignored. They are described, however, in Appendix A. The 2 nd term in Eq. (14) describes the Coulomb potentials of D − defects, that have an effective negative charge and represent, for instance, Ca 2+ substituting R 3+ ions. Defects attract doped holes and strongly repel electrons of V-ions in the vicinity and shift these states from the lower Hubbard band into the Mott gap [24]. The 3 rd term, the e-e interaction, leads to the screening by t 2g electrons and doped holes. Both terms are determined by the Coulomb field, v(r) ≡ e 2 /ε c r, where ε c 5 is the dielectric constant of the core electrons [57] and r is the distance between charges of: (i) a defect D − at site m and t 2g electrons at a V ion at site i, i.e., r mi = |R m − r i |, and (ii) two V ions at sites i and j with r ij = |r i − r j |.n i = ασn iασ is the t 2g electron charge operator withn iασ =d † iασd iασ and orbital flavors α = {yz, zx, xy} ≡ {a, b, c} [42], see Fig. 1(a).
Central to our discussion are the orbital polarization terms, H pol ≡ H (1) pol + H (2) pol . They describe the OR and the redistribution of t 2g electronic charge at V-ions, induced by the Coulomb fields of defects and doped holes. These terms appear in addition to the monopole terms that are already contained in the minimal model. The orbital polarization term H (1) pol , due to the charged defect [58] reads as: where the coupling constant D is defined by the matrix element iα|v(|r−R m |)|iβ ≡ Dλ d αβ in the basis α = {a, b, c}, see Fig. 1(a), with d = r i − R m . We shall treat D as a free parameter; a typical value is D ≈ 50 meV [59]. The effect of orbital rotation is short-ranged and affects only the 8 V ions of the defect cube C m of the defect m. The matrix elements λ d αβ are traceless, like the 3-flavor SU(3) matrices [68,69]. They depend on the diagonal axis d in the defect cube, i.e.,  (111), respectively. In the large D limit, the unoccupied orbital |a = (|a + |b + |c )/ √ 3 shown in Fig. 1(c) follows from H (1) pol alone. It has the largest overlap with the negative defect and thus the highest energy.
The orbital polarizations H (2) pol induced within the t 2g orbitals of V-ions by the quadrupolar components of the Coulomb fields of doped holes is the central issue of this article. This perturbation of the orbital order, as well as the competition with orbital rotations induced by defects, has not been explored before. The perturbation due to the polaron charge results from e-e interactions, where n 0 − n j measures the hole-density relative to the undoped system: Here, n 0 = 2 and we restrict the sum over r i to a neighborhood N j that includes the 6 V neighbors of the hole at j as shown in Fig. 1(e). A detailed derivation is given in the Appendix. The coupling constant D h is defined by the matrix element of the field of the hole at r j : with respect to the orbital basis a = yz, b = zx, and c = xy: where ζ h αβ depends on the axis h||(100), (010) and (001), respectively. That is, for a (001) V-ion neighbor of a hole at r j , the c orbital energy is raised by 2D h , while a(b) are lowered by D h , see Fig. 1(f). This term frustrates the polarization around charged defects (22).
The level splittings induced by the polaron charge are illustrated in Fig. 2 where the hole was inserted into the b orbital of the ion V 0 in the center of the polaron. In the GO state, all nearest neighbors of V 0 have the c and a orbital occupied. The actions of H (2) pol along the different cubic directions are different, however. To see this anisotropy of the orbital polaron it is useful to consider the in matrix notation of ζ h αβ which depends on the axis h||(100), (010) and (001). For instance, at the (001) vanadium-ion neighbor, labeled V 3 , of the doped hole at the ion V 0 the c orbital energy is raised by 2D h , while a(b) are lowered by D h . This term eliminates the a/b-orbital polarization at V 3 ions. The other cubic directions are different, for instance along the b axis the occupation of V 2 does not change at all. Yet along the a-axis the a orbital of V 1 is raised by 2D h while the others are lowered. This leads to switching from a to b occupation as displayed in V 1 . We note that this switching is particularly harmful for the a/b orbital order, as it favors the inverted order.
An important consequence of the orbital excitation from a to b occupation of ion V 1 when D h crosses the switching value is a blocking of orbital super-exchange processes along the z direction. This is indicated by a red cross on the V 1 -V 4 bond. Since the same occurs on the complementary z-bonds centered at two equivalent V 1 ions the total energy increase, or loss of negative super-exchange energy, corresponds to 4 broken bonds. Each bond corresponding to a virtual kinetic energy of ∆ 0 = E 0 kin = 0.10 eV for t = 0.2 eV. Thus the total increase of kinetic energy is then 0.4 eV, which coincides with the kinetic energy change per polaron , as we shall see further below.
Moreover, the occupation change at V 1 does not involve the change of crystal field energies, although there are some changes of Jahn-Teller energies which we neglect here in the discussion (but not in the calculations). The latter are however small compared to the change of super-exchange energy. Therefore, we can conclude that the crossover scale D c h is not determined by crystal field and/or Jahn-Teller terms, but by the interplay of the OR terms and the super-exchange energy, that is, the scale is determined by the quantum dynamics of the electrons, and is confirmed by our numerical study.

III. ORBITAL ROTATIONS INDUCED BY DEFECTS AND POLARONS: RESULTS
In this study, we explore the spin-orbital order of vanadates in the insulating regime, where doped holes are bound to random defects and typically prefer a V site on a defect cube;which site, depends on the Coulomb interactions with all other defects and doped holes. The latter generate polarons and form defect states inside the Mott gap that persist up to high doping [24]. We calculate the disordered electronic structure using a variant of the unrestricted Hartree-Fock method [70,71] that obeys rotational invariance in both spin and orbital space [58,72], emphasized as well in H Hub [73] and in slave-boson theories [74][75][76]; this formulation preserves the multiplet structure of atoms and ions and thereby avoids the shortcomings of the non-rotational invariant formalism. We first discuss the effect of H pol and H pol on the orbital densities, n ab and n c . In Fig. 3, we monitor three cases, namely the separate effects of defect and polaron induced OR, as well as the combined effect, where we use the geometrical relation, D h /D = ξ ∼ 0.87, defined by the ratio ξ = d/a of nearest neighbor V-D (d) and V-V (a) distances. For D = 0 = D h , doped holes go into a and b orbitals due to the crystal field ∆ c , i.e., n ab = 1 − x and n c = 1. The increase of n ab versus D (at D h = 0) can be qualitatively understood from the rotation of occupied states {|c , |b } into |c = (2|c − |a − |b )/ √ 6 and |b = (|b − |a )/ √ 2 in the large D limit and at t = 0, see Fig. 1(c).
The transfer of holes into the c orbitals due to D h (D = 0) is induced by an upward shift of a c-orbital as shown in Fig. 1(f). From H (2) pol , one recognizes that, along a 2 nd polaron axis, there is no change of occupation and, along the 3 rd axis, there is an interchange of a and b orbital occupation. It is this latter mechanism that is particularly harmful for the a/b orbital order.
The change of G-type orbital order parameter with D and D h is shown in Fig. 4(a). It is determined by the spatial modulation of the local occupation numbers, n ia and n ib , where the disorder average is typically taken over M = 100 defect realizations s and Q G ≡ (π, π, π). Typical system sizes are N = 4 3 or 8 3 . For small D h = ξD, the behavior of m o ab is very similar to the pure D (D h = 0) case. The very different behaviors of the pure polaron D h (D = 0) and defect induced rotation results from the diagonal versus off-diagonal nature of ζ h αβ and λ d αβ , respectively. This explains that, in spite of frustration, the qualitative decay of orbital order of the combined action of D h = ξD is similar to the ORs due to defects, yet the OR due to polarons leads to evident effects for large D.
The drop of orbital order in Fig. 4(a), triggered by the po- The magnetic anisotropy of the CS/GO phase of the parent compounds is a manifestation of the strong orbital fluctuations along the z-axis [47,77,78]. This implies a much larger virtual kinetic energy of a and b orbitals along z as compared to the virtual hopping along x. For the undoped parent state, we find a large anisotropy, A ≡ K z /K x ∼ 1. cubic directions. This anisotropy A is similar to the ratio of measured optical weights A opt ∼ 1.84 for LaVO 3 [65]. Next, we explore the change of the kinetic energy components K z and K x as functions of x and D h for fixed D = 0.05 eV. The upward shift of K z with hole-doping x in Fig. 5(a) reflects the loss of super-exchange or binding energy. A loss is further amplified by the D h term. The total K x values in Fig. 5(b) lie close to −0.06 eV and show only a marginal dependence on x and D h . The small changes of the components K ab x and K c x with D h reflect the transfer of holes from {a, b} to c orbitals. Hence, we find that the change of the total kinetic energy K with doping x is almost completely determined by K z , which results from large {a, b} orbital fluctuations along z that favor G-type orbital order. The data in Fig. 5 implies a decrease of the anisotropy A with increasing polaron parameter D h (at D = 0.05 eV), however most of the reduction of A with doping results from D, i.e., the OR clouds induced by the defects. It is worth noting that, it is much harder to reach self-consistency at moderate D h -due to the motion of holes-than in calculations with defect-induced ORs alone.
The contour plot of the anisotropy A = K z /K x of kinetic energies in Fig. 6(a) for D h = 0 displays a strong reduction towards A ∼ 1 when both the defect induced OR parameter D and doping x become sufficiently large. At x = 0.2 and D = 0.06 eV, the asymmetry A ∼ 1.3, whereas in absence of ORs, for D = 0 = D h , the decay of A with x is much weaker. In the latter case, it is caused exclusively by the motion of doped holes bound to defects as small spin-orbital polarons [24]. Figure 6(b) shows the decay of A versus doping x for different values of D and provides a comparison with the optical anisotropy A opt as determined for La 1−x Sr x VO 3 [41], the only system where such data seems to exist. It is remarkable that the theoretical A and experimental A opt almost coincide for the undoped system. The most pronounced discrepancy is a tendency of A opt towards a cooperative transition.
However, the cooperative nature of the decay of orbital order in the (La,Sr) system appears as an exception; experiments for (Pr,Ca), (Nd,Sr), and (Y,Ca) show a gradual decline of the order parameter with x [41,56]. We saw above that ORs induced by defects act non-cooperatively, consistent with such a gradual decline. Yet, ORs induced by polaron charges may well lead to cooperative transitions as they are driven by ee interactions; accordingly, they can be extremely relevant to fine tune the theory to specific compounds.

IV. CONCLUSIONS
The robustness of the insulating state and of the G-type orbital order in the vanadates, observed in several experiments [41], has two main causes: (i) Doped holes are localized by defects and form small spinorbital-polarons. The main kinetic energy gain of a doped hole is a double exchange process on an active bond (a ferromagnetic bond in c-direction) next to a defect. (ii) Orbital order is predominantly suppressed as function of doping x not by the kinetic energy associated with the small polaron but rather by non-cooperative orbital rotations induced by defect charges. An important feature of both orbital rotation mechanisms, namely the orbital rotations due to the charges of defects and TABLE I. Doping dependence of occupation numbers n ab and nc in the atomic limit (t → 0) for different limiting cases for D h and D. Here, D sat h and D sat denoted saturation values where the orbital rotation is complete. of polarons, is the transfer of holes from the a/b orbitals to c orbitals. In the absence of these terms, in the atomic limit t → 0, one finds n c = 1 and n ab = 1 − x. Where the latter relation shows that doped holes go into a/b orbitals. Table I summarizes the transfer of holes in two further limits; for instance, if D h is in the saturation regime (and D = 0) the number of electrons in the a/b sector becomes even larger than one, that is n ab = 1 + x, whereas there are now even more holes in the c-orbitals, i.e., n c = 1 − 2x, than expected from the doping concentration x. Figure 7 displays the decay of the G-type orbital order parameter m o ab as function of D h . Interestingly, it also shows a slight increase of the complementary C-type spin order in the ab orbital sector, m s ab . This perhaps surprising increase of m s ab is related to the transfer of holes from {a, b} to c orbitals due to ORs. It is reminiscent of the peculiar robustness of the long-range CS order in the large doping regime where orbital order is still present, but short ranged [56].
In summary, we explored the competition between the orbital rotations induced by the polaron charges and those induced by the defect charges, a priori both being of similar importance. We found that these rotations are the key mechanisms that control the decay of orbital order of t 2g electrons in doped vanadates, -much more important than the string mechanism [11,[13][14][15] active in high-T c superconductors. When they act together, the qualitative suppression of orbital order appears to be mainly controlled by the off-diagonal rotations (22) due to defect charges. We found that the energy scale for rotation of orbitals is primarily determined by quan-tum orbital fluctuations rather than by classical crystal fields and Jahn-Teller potentials. It is very surprising that the suppression of the anisotropy of kinetic energy under increasing doping x occurs in the still-anisotropic C-type magnetic state. In this Appendix, Subsection A, we summarize the detailed, minimal multi-orbital Hubbard model for the description of the spin and orbital degrees of freedom which give rise to the rich phase diagram of the RVO 3 vanadates with perovskite structure. We confine ourselves to the space spanned by the t 2g electrons of vanadium. An extension of the minimal model, i.e., by implementation of the Coulomb fields of defects and the electron-electron interactions in Subsection B allows to study the doped systems. By the inclusion of the e-e interactions the extra screening of the defect potentials by the doped electrons and/or holes are taken into account within the UHF approach. In Subsection C and D we discuss the derivation of the extra quadrupolar terms that result from the Coulomb fields of defects and among electrons, respectively. Perhaps it is useful here to remind the reader that the monopole Coulomb fields due to defects and e-e interactions are explicitely contained in the minimal model as described in Appendix B.
Our general approach follows a similar route as the manybody treatment of high-T c superconductors where one also starts from a minimal model, namely the planar one-band Hubbard model [2] -although the model for the vanadates here is more complex. The Hubbard model of cuprates contains the spin-1 2 Heisenberg model, the doped holes and their interactions. All states not directly related to Cu(d 9 ), in particular the O(2p) states, have been "integrated out". Their effect is still present in the form of renormalized effective parameters, like for instance the hopping parameters t ij . As long as one is interested in the low energy and low temperature physics this approach is fully justified. Only if one considers spectroscopies in the energy window of pd-transitions, where the oxygen states come into play, then the model is not sufficient to describe those, as they have been integrated out.
A. The three-orbital flavor Hubbard model The three-orbital Hubbard model for the t 2g electrons was introduced for the triangular lattice [66] and adopted later for the pnictide superconductors [67]. It has very rich physics for electronic orders as shown recently [79]. Here, we use it for doped vanadium (La,Y) 1−x Ca x VO 3 perovskites [24,57,58], Its main part consists of the kinetic energy H kin and of the local interactions between the electrons in the three t 2g orbitals, H int . It describes the situation in the vanadium perovskites after supplementing it by rather weak terms [49,53]: the crystal field (CF) splitting H CF , and the Jahn-Teller (JT) interactions H JT [58]. Below we use the definition of t 2g orbital degrees of freedom which selects uniquely a single cubic direction along which the hopping is inactive to label each orbital flavor [42], The kinetic energy for t 2g electrons preserves the orbital flavor in the hopping along the bond ij γ oriented along one of the cubic axes, γ ∈ {a, b, c}. It reads, Here,d † iασ is the electron creation operator in the t 2g orbital α ∈ {xy, yz, zx} with spin σ =↑, ↓ at site i. The effective hopping (t > 0) occurs in two steps, via the hybridization to an intermediate oxygen 2p π orbital, along idealized 180°V-O-V bonds. Therefore, the hopping is: (i) diagonal and conserves the orbital flavor α when the hybridization with the oxygen 2p π orbitals is finite, and (ii) zero otherwise, i.e., t γα ij = 0 if the hybridization with the oxygen 2p π orbitals vanishes by symmetry.
Local interactions at vanadium ions, H int , are rotationally invariant in the orbital space [73] and depend on two Kanamori parameters: (i) intraorbital Coulomb interaction U and (ii) Hund's exchange J H between each (equivalent) pair of t 2g electrons in different orbitals, Interorbital Coulomb interactions ∝n iαniβ are expressed in terms of orbital electron density operators for a pair α < β,n iα = σniασ = σd B. Coulomb fields and orbital rotations due to charged defects and polarons The complete Hamiltonian for t 2g electrons in doped vanadium (La,Y) 1−x Ca x VO 3 perovskites, Eq. (1), reads [24], It includes the three-band Hubbard model H Hub [49]-see above -, the long-range electron-electron and defect-electron interactions ∝ v(r) -see main text -and orbital polarization terms, The first term was analyzed before and is responsible for the collapse of orbital order under doping by charged defects [59].

C. Defect induced orbital polarization
The term of the Hamiltonian for vanadium ions at r i in the Coulomb potential of defects D − at R m is: Introducing V αβ mi = iα|v(|r−R m |)|iβ (17) andV αα mi = one obtains pol .
where the sums in the first term extend over the whole system. The range of the polarization term, due to the short-range nature of the matrix elements, will be restricted to the defect cube of the respective defect C m . The polarization term (20) consists of a diagonal and an off-diagonal terms. For a nondistorted cubic neighborhood only the off-diagonal terms contribute. where the coupling constant D is defined by the matrix element Dλ d αβ ≡ V αβ mi = iα|v(|r−R m |)|iβ .
Here λ d αβ contains the signs of the matrix elements and is displayed in the main text. The signs do depend on the respective diagonal of the defect cube, i.e., parallel to the vector d = r i −R m connecting the respective V-ion and the defect. The rotation operator can then be summarized as A value D ≈ 50 meV has been estimated by simple defectpotential-mediated superposition integrals, given a strength of the defect potential at the vanadium sites (at a given t 2g orbital) of 2 eV.

D. Orbital polarization due to the polaron charge
The orbital polarization induced on the vanadium ion at r i due to the charge of polaron at r j stems from e-e interactions. The leading term H (2) pol has monopole-quadrupole character. Similar to the defect case we can write (23) Here ( n j − n 0 ) represents the negative density of doped holes, and n 0 is the number of electrons per V-ion in the undoped case, i.e., for the d 2 configuration n 0 = 2. The coupling constant D h follows from the two-center matrix element V αβ ji = iα|v(|r−r j |)|iβ .
For the undistorted cubic system only diagonal terms contribute, that depend on the vector h = r i −r j connecting V-ion and the polaron density at r j , where we include only the nearest V-ions in C j . The coupling constant D h is defined by the matrix element of the field of the hole at r j and the orbitals at r i : We find that for an undistorted cubic lattice ζ h αβ is diagonal with respect to the global t 2g basis and depends on the cubic axis γ(h), i.e., parallel to the vector h = r i − r j that connects the polaron density at r j and a n.n. V-ion at r i . In our study, we shall use either D = 0 or 50 meV with D h as a free parameter. Alternatively, we use the geometric relation D h /D (d/a) γ defined by the nearest neighbor V-V and the V-D distance, labeled as a and d, respectively. This corresponds to 0.87 for γ = 1, which we use here.