Spin-Orbital Order Modified by Orbital Dilution in Transition Metal Oxides: From Spin Defects to Frustrated Spins Polarizing Host Orbitals

We study the $3d$ substitution in $4d$ transition metal oxides in the cases of $3d^3$ doping at either $3d^2$ or $4d^4$ sites which realize orbital dilution. We derive the effective $3d-4d$ (or $3d-3d$) superexchange in a Mott insulator with different ionic valencies, underlining the emerging structure of the spin-orbital coupling between the impurity and the host sites and demonstrate that it is qualitatively different from that encountered in the host itself. This derivation shows that the interaction between the host and the impurity depends in a crucial way on the type of doubly occupied $t_{2g}$ orbital. One finds that in some cases, due to the quench of the orbital degree of freedom at the $3d$ impurity, the spin and orbital order within the host is drastically modified by doping. The impurity acts either as a spin defect accompanied by an orbital vacancy in the spin-orbital structure when the host-impurity coupling is weak, or it favors doubly occupied active orbitals (orbital polarons) along the $3d-4d$ bond leading to antiferromagnetic or ferromagnetic spin coupling. This competition between different magnetic couplings leads to quite different ground states. We find that magnetic frustration and spin degeneracy can be lifted by the quantum orbital flips of the host but they are robust in special regions of the incommensurate phase diagram. The spin-orbit coupling can lead to anisotropic spin and orbital patterns along the symmetry directions and cause a radical modification of the order imposed by the spin-orbital superexchange. Our findings are expected to be of importance for future theoretical understanding of experimental results for doped $4d$ transition metal oxides doped with $3d^3$ ions. We suggest how the local or global changes of the spin-orbital order induced by such impurities could be detected experimentally.


I. INTRODUCTION
The investigation of strongly correlated electrons in transition metal oxides (TMOs) focuses traditionally on 3d materials [1], mainly because of high-temperature superconductivity discovered in cuprates and more recently in iron-pnictides, and colossal magnetoresistance manganites, but also because the more extended 4d orbitals would a priori suggest a weaker correlation due to a reduced ratio between the intraatomic Coulomb interaction and the electron bandwidth. Nevertheless, the extension of the 4d-shells couples the 4d orbitals and the neighboring oxygen orbitals, implying that these TMOs have the tendency to form distorted structures with respect to the ideal ones. As a consequence, the change in the 3d-2p-3d bond angle often leads to narrowing of the (effective) d-bandwidth, bringing the system on the verge of a metal-insulator transition, or even into the Mott insulating state. Hence, 4d materials share common features with 3d systems having additionally a significant sensitivity of the electronic states to the lattice structure and to relativistic effects due to larger spin-orbit or other magneto-crystalline couplings.
Based on the large variety of physical phenomena exhibited by 4d and 3d-based oxides, the synthesis of hybrid compounds having both transition metal ions may open a novel route for unconventional effects in complex materials. There are several reasons for expecting new scenarios in mixed 3d − 4d spin-orbital-lattice materials. The difference between the 3d and 4d electronic shells makes the chemical substitution a physical parameter for tuning the amplitude of the intraatomic Coulomb interaction when compared to the effective d − d electron transfer. Then, the atomic numbers for the 4d and 3d atoms would also naturally point towards a different role of the spin-orbit coupling. Finally, the possibility of having nonequivalent magnetic moment and orbital degree of freedom at the 3d and 4d sites allows one to realize a unique doping by spin and/or orbital dilution with a significant impact on the coupling between the spin-orbital-lattice degrees of freedom.
As a demonstration of the highly nontrivial physics emerging in 3d − 4d oxides, remarkable effects have been observed when, for instance, Ru ions are replaced by Mn, Ti, Cr or other 3d elements. The role of Mn doping in the arXiv:1408.1838v1 [cond-mat.str-el] 8 Aug 2014 SrRuO Ruddlesden-Popper series is strongly linked to the dimensionality through the number n of RuO 2 layers in the unit cell. The Mn doping of the SrRuO 3 cubic member drives the system from the itinerant ferromagnetic state to an insulating antiferromagnetic (AF) configuration in a continuous way via a possible unconventional quantum phase transition [23]. Doping by Mn ions in Sr 3 Ru 2 O 7 leads to a metal-to-insulator transition and AF long-range order for more than 5% Mn concentration [24]. Subtle orbital rearrangement can occur at the Mn site, as for instance the inversion of the crystal field in the e g sector observed via x-ray absorption spectroscopy [25]. As a function of doping neutron scattering studies indicate the occurrence of an unusual E-type antiferromagnetism (planar order with zigzag FM chains coupled antiferromagnetically) with moments aligned along the c axis within a single bilayer [26].
Similarly, for the isovalent Ca-based ruthenate bilayer, small amount of Ti doping in Ca 3 Ru 2 O 7 leads to an insulating state with a G-type nearest neighbor(NN) AF order [27]. In contrast to Mn doping, Cr ions can substitute Ru either in the trivalent (3d 3 ) or the tetravalent (3d 2 ) states providing alternative way of spin and orbital modification in the 4d host. Cr substitution for Ru generally results in no structural change but does induce interesting magnetic behavior with a tendency to itinerant ferromagnetism as demonstrated by SrRu 1−x Cr x O 3 [28][29][30], CaRu 1−x Cr x O 3 [31], and Ca 3 (Ru 1−x Cr x ) 2 O 7 [32]. Surprisingly, Cr doping for Ru in Ca 2 RuO 4 not only leads to spin canting and ferromagnetism but it also strongly reduces the structural distortions. Anomalous negative volume expansion is observed in Ca 3 (Ru 1−x Cr x ) 2 O 7 with a surprising Invar effect in the region where the system is supposed to be orbitally and magnetically ordered [33]. Hence, the 3d chemical substitution in 4d oxides naturally provides paths for unusual metal-to-insulator transitions accompanied by magnetic and orbital orders with a large potential for exploring the complexity of the electronic ground state in strongly correlated materials.
The competition of different and complex types of order is ubiquitous in strongly correlated TMOs mainly due to frustrated exchange. Doping competes with this frustrated superexchange interactions and generates complex types of spin or charge order [1,34]. This happens, for instance, in doped cuprates, where stripe order emerges as a compromise between the magnetic and kinetic energy [35,36], and remarkable evolution of the stripes and their spectral properties is found under doping [37][38][39][40]. Also in systems with orbital order stripe phases occur by a similar mechanism which involves orbital polarons [41,42]. A prerequisite to these phenomena are the undoped systems with active orbital degrees of freedom and unpolarized spins, where the low-energy physics and spin-orbital order are dictated by effective spin-orbital superexchange [43][44][45][46][47][48][49]. In such cases an interesting interplay between spin and orbital degrees of freedom takes place [50] and certain exotic types of spin order may arise even without doping as a consequence of frustrated interactions [51][52][53][54][55][56][57].
In doped La 1−x (Sr,Ca)Mn x O 3 manganites several different types of order are observed at increasing doping [58][59][60] which require careful modeling that takes into account both 3d 3 and 3d 4 electronic configurations of Mn ions [61][62][63][64][65][66]. The orbital order changes radically at the phase transitions found between different types of magnetic order at increasing doping [67]. Situations such as in the FM metallic phase where the spin and orbital degrees of freedom separate and spin excitations are driven by the orbital liquid [68,69] are exceptional, and in general spin-orbital entangled states play a role [50], as in the case of cubic vanadates where orbital fluctuations lead to several exotic phenomena [70][71][72][73][74][75][76][77]. Doping by Ca (Sr) ions replacing Y (La) ones in YVO 3 (LaVO 3 ) changes the magnetic order [78,79] and generates spectral weight within the Mott-Hubbard gap [78,80]. It is quite remarkable that due to rather different magnetic and kinetic energy scales, even low doping may suffice for changing drastically the magnetic order, as for instance in electron-doped manganites [81,82].
In the present work we focus on the insulating phases of 3d − 4d hybrid structures, where doping happens directly at transition metal sites, and the value of the spin is locally changed. In fact, this doping not only modifies the number of electrons in the ionic states but also influences directly the orbital degrees of freedom in the undoped compound. To simplify the analysis we assume that onsite Coulomb interactions are strong and charge degrees of freedom are thus projected out, and only virtual charge transfer can occur between 3d and 4d ions via the oxygen ligands. As an example of the physical situation in such systems, we deal below with 3d 3 doping for 4d 4 atomic configurations which can be of interest for the case of Mn or Cr doped ruthenates. The undoped system is a hole analogue of cubic vanadates [70,73], with t 2g orbital degree of freedom and S = 1 spins per site in both cases. For convenience, we shall define the orbital degree of freedom as a doublon (double occupancy) in the d 4 configuration which corresponds to an empty t 2g orbital in the d 2 configuration in a vanadate. This doping leads thus effectively to the removal of a t 4 2g ion with the corresponding double occupancy (doublon) in one of t 2g orbitals which we label as {a, b, c} (this notation is explained below) and to replacing it by a t 3 2g ion. Unlike in electron doped manganites where orbital degree of freedom is added [81,82], doping by d 3 ions of a d 4 system removes the orbital degree of freedom at some sites and results in orbital dilution and creating vacancies in the doublon order within the host. Such defects, on one hand, can weaken the spin-orbital coupling in the host and, on the other hand, may open a new channel of interaction between the magnetic and orbital degree of freedom through the host-impurity exchange, see Fig. 1.
Before presenting the details of the quantitative analysis, let us concentrate of the main idea of superexchange modified by doping in a spin-orbital system. to Hund's exchange. While a pair of d 3 ions, e.g. in SrMnO 3 , is similar to the t-J model and has AF superexchange [69], the superexchange for the d 3 − d 4 has a rather rich structure and may also be FM. The spin exchange depends then on whether the orbital degree of freedom is active and participates in charge excitations along a considered bond or electrons of the doublon cannot move along this bond due to the symmetry of t 2g orbital, as explained in Fig. 1. This qualitative difference to systems without active orbital degrees of freedom is investigated in detail in Sec. II. The main outcomes of our analysis are: (i) the determination of the effective spin-orbital exchange Hamiltonian describing the low-energy sector for the 3d − 4d hybrid structure, (ii) establishing that a 3d 3 impurity without an orbital degree of freedom modifies the orbital order in the 4d 4 host, (iii) providing the detailed way how the microscopic spin-orbital order within the 4d 4 host order is modified around the 3d 3 impurity, and (iv) suggesting possible spin-orbital patterns that arise due to a periodic and finite substitution (doping) of 3d atoms in the 4d host. The emerging physical scenario is that the 3d impurity acts as an orbital vacancy when the host-impurity coupling is weak and as an orbital polarizer of the bonds active t 2g doublon configurations when it is strongly coupled to the host. The tendency to orbital polarize the host around the impurity turns out to be independent of the spin configuration. Otherwise, it is the resulting or-bital arrangement around the impurity and the strength of Hund's coupling at the impurity that set the character of the host-impurity magnetic exchange.
The strategy we adopt is to firstly analyze the ground state properties of a single 3d impurity surrounded by 4d atoms by investigating how that spin-orbital pattern can be modified in the host at the NN sites to the 3d atom. This study is performed for different spin-orbital patterns of the 4d host with special emphasis on the FM chains alternating antiferromagnetically (C-type order) with alternating orbital (AO) order. Starting from the single-impurity solution we address the problem of a periodic arrangement of 3d doping at different concentrations. We demonstrate that the spin-orbital order in the host can be completely modified by the presence of impurities leading to striped patterns with alternating FM/AF domains and islands of fully FM states.
The remaining of the paper is organized as follows. In the Sec. II we introduce the effective model describing the spin-orbital superexchange at the 3d−4d bonds which serves to investigate the changes of spin and orbital order around individual impurities and at finite doping. We arrive at a rather general formulation which emphasizes the impurity orbital degree of freedom, being a doublon, and present some technical details of the derivation in the Appendix. Section III is devoted to the analysis of a single impurity case. We focus on the local effects in the spin-orbital order of the host induced around the impurity and present the respective cluster approximation in Sec. III A. As explained in Sec. III B, there are two nonequivalent cases which depend on the precise modification of the orbital order by the 3d impurity, doped either to replace a doublon in a orbital (Sec. III C) or in c orbital (Sec. III D). In Sec. IV A we consider the changes of spin-orbital order which arise at periodic doping with macroscopic concentration. Here we limit ourselves to two representative cases: (i) commensurate x = 1/8 doping in Sec. IV B, and (ii) incommensurate x = 1/5 and x = 1/9 doping in Secs. IV C and IV D. A general discussion and main conclusions are presented in Sec. V.

II. THE MODEL
In this Section we consider a 3d impurity in a strongly correlated 4d system and derive the effective 3d 3 − 4d 4 spin-orbital superexchange. It follows from the coupling between 3d and 4d orbitals due to the p − d hybridization. We implement a strict rule that the hopping within the t 2g sector is allowed in a TMO only between two neighboring orbitals of the same symmetry which are active along the bond direction [83][84][85], and neglect the interorbital processes originating from the octahedral distortions such as rotation or tilting. The interorbital hopping elements are smaller by at least one order of magnitude and may be treated as corrections to the overall scenario established below.
In a strongly correlated system it suffices to concen-trate on a pair of atoms forming a bond 1, 2 , as the the effective interactions depend on the bond [86]. For convenience, we label the impurity 3d atom by m = 2. The effective Hamiltonian for a representative 3d-2p-4d bond, after projecting out the oxygen degrees of freedom, can be expressed in the following form, Here we adopt the notation for a bond connecting the 3d transition metal ion at site m = 1 with the 4d one at site m = 2 along a given crystallographic direction. The total Hamiltonian contains the kinetic energy term H t describing the electron charge transfer via oxygen orbitals, the onsite interaction terms H int (m) for the 3d (4d) ion at site m = 1, 2, respectively, and the local potential of the 4d atom, H ion (2), which takes into account the mismatch of the energy level structure between the two atomic species and prevents valence fluctuations when the host is doped even in the absence of local Coulomb interaction. The kinetic energy is given by, where d † mµσ are the electron creation operators at site m = 1, 2 in the spin-orbital state (µσ). The bond 12 points along one of the two crystallographic directions, γ = a, b, in the square lattice. Without distortions, only two out of three t 2g orbitals are active along each bond 12 and contribute to H t (1, 2) [83][84][85], while the third orbital lies in the plane perpendicular to the γ axis thus the hopping via oxygen is forbidden by symmetry. This motivates a convenient notation as follows, |a ≡ |yz , |b ≡ |xz , |c ≡ |xy , with the orbital inactive along a given direction γ ∈ {a, b, c} labelled by the index γ. We consider a twodimensional (2D) square lattice with transition metal ions connected via oxygen orbitals as in a RuO 2 ab plane. In this case |a (|b ) orbitals are active along the b (a) axis, while |c orbitals are active along both axes, a and b.
The Coulomb interaction on an atom at site m = 1, 2 depends on two parameters [87]: (i) intraorbital Coulomb repulsion U m , and (ii) Hund's exchange J H m . The label m = 1, 2 stands for the ion and distinguishes between these terms at the 3d and 4d ion, respectively. The interaction is expressed in the form, The terms standing in the first line of Eq. (4) contribute to the magnetic instabilities in degenerate Hubbard model [87] and decide about spin order, both in an itinerant system and in a Mott insulator. The remaining terms contribute to the multiplet structure and are of importance for the correct derivation of the superexchange which follows from charge excitations, see below.
Finally, we include a local potential on the 4d atom which encodes the energy mismatch between the host and the impurity orbitals close to the Fermi level and prevents valence fluctuations on the 4d ion due to the 3d doping. This term has the following general structure, with µ = a, b, c. The effective Hamiltonian for the low energy processes is derived from H(1, 2) (1) by a second order expansion for charge excitations generated by H t , and treating the remaining part of H(1, 2) as an unperturbed Hamiltonian. We are basically interested in virtual charge excitations in the manifold of degenerate ground states of a pair of 3d and 4d atoms on a bond, see Fig. 2. These quantum states are enumerated as e k 1 with k = 1, . . . , 4 and {e p 2 } with p = 1, . . . , 9 and their number follows from the solution of the onsite quantum problem for the Hamiltonian H int (m). For the 3d atom the relevant states can be classified according to the four components of the total spin S 1 = 3/2 for the 3d impurity atom at site m = 1 and S 2 = 1 for the 4d hor atom at site m = 2 for the three different positions of the double occupied orbital (doublon). Thus, the effective Hamiltonian will contain spin products ( S 1 · S 2 ) between spin operators defined as, for m = 1, 2 sites and the operator of the doublon position at site m = 2, The doublon operator identifies the orbital γ within the t 2g manifold of the 4d ion occupied by the doublon. It is worth noting that the hopping (2) does not change the orbital flavor thus we expect that the resulting Hamiltonian will be diagonal in the orbital degrees of freedom with only D γ 2 operators. Following the standard second order perturbation expansion for spin-orbital systems [45], we can write the matrix elements of the low energy exchange Hamiltonian, H (γ) J (12), for a bond 12 γ along the γ axis as follows, j one which adds another doublon at the 3d site in the intermediate state. The second type of excitations involves more doubly occupied orbitals and has much larger excitation energy. It is therefore only a small correction to the leading term (i), as we discuss in the Appendix.
Similar as in the case of doped manganites [68], the dominant contribution to the effective low-energy spinorbital Hamiltonian for the 3d − 4d bond stems from the j charge excitations, as they do not involve an extra double occupancy and the U 2 Coulomb energy. The 3d 3 i 4d 4 j 3d 4 i 4d 3 j charge excitations can be analyzed in a similar way as the 3d 3 i 3d 4 j 3d 4 i 3d 3 j ones in doped manganites [68]. In both cases the total number of doubly occupied orbitals does not change, so the main contributions come due to Hund's exchange. In the present case, one more parameter plays a role, which stands for the mismatch potential energy (5) renormalized by the onsite Coulomb interactions {U m } and by Hund's exchange. On a general ground we expect ∆ to be a positive quantity, since the repulsion U m should be larger for smaller 3d shells than for the 4d ones and U m is the largest energy scale in the problem. Let us have a closer view on this dominant contribution of the effective low-energy spin-orbital Hamiltonian for the 3d − 4d bond, given by Eq. (A.2). For the analysis performed below and the clarity of our presentation it is convenient to introduce some scaled parameters related to the interactions within the host and between the host and the impurity. For this purpose we employ the exchange couplings J imp and J host , and use their ratio to investigate the influence of the impurity on the spin-orbital order in the host. Here t 4d is the hopping amplitude between two t 2g orbitals at NN 4d atoms, J H 2 and U 2 refer to the host and ∆ (9) is the renormalized ionization energy of the 3d−4d bonds. The results depend as well on Hund's exchange element for the impurity and on the one at host atoms, Note that the ratio introduced for the impurity, η imp (11), has here a different meaning from Hund's exchange used here for the host, η host (12), which cannot be too large by construction, i.e., η host < 1/3. With the parameterization introduced above, the dominant term in the impurity-host Hamiltonian H where the orbital (doublon) dependent spin couplings J S (D 2,γ ) and the doublon energy E D depend on η imp . The evolution of the exchange couplings are shown in Fig. 3. We note that the dominant energy scale is E γ D , so for a single 3d−4d bond the doublon will avoid occupying the inactive (γ) orbital and the spins will couple with J S (D 2,γ = 0) which can be either AF if η imp 0.43 or FM if η imp > 0.43. Thus the spins at η imp = η c imp 0.43 will decouple according to the H (γ) J (12) exchange. Let us conclude this Section by writing the effective spin-orbital Hamiltonian for the 4d host for the bonds ij along the γ = a, b direction [88], with J (γ) ij operators acting only in the orbital space. They are expressed in terms of the pseudospin defined in the orbital subspace spanned by the two orbital flavors active along a given direction γ, i.e., with standing for the multiplet structure in charge excitations, and the orbital operators { τ i } that for the γ = c axis take the form: For the directions γ = a, b in the considered ab plane one finds equivalent expressions by cyclic permutation of the axis labels {a, b, c} in the above formulas. This problem is isomorphic with the spin-orbital superexchange in the vanadium perovskites [73,74], where a hole in the t 2g shell plays an equivalent role to the doublon in the present case. The operators are the doublon (hard core boson) creation operators in the orbital γ = a, b, c, respectively, and they satisfy the local constraint, meaning that exactly one doublon can occupy one of the three t 2g orbitals at each site i. We use the doublon number operators, below in the classical procedure to determine the ground state of single impurities and doped systems. In this Section we describe the methodology that we applied for the determination of the phase diagrams for a single impurity reported below in Secs. III C, and next at macroscopic doping, as presented in Sec. IV. Let us consider first the case of a single 3d impurity in the 4d host. Since the interactions in the model Hamiltonian are only effective ones between NN sites, it is sufficient to study the modification of the spin-orbital order around the impurity for a given spin-orbital configuration of the host by investigating a cluster of N = 13 sites shown in Fig.  4. We assume the C-AF spin order (FM chains coupled antiferromagnetically) accompanied by AO order within the host which is the spin-orbital order occurring for the realistic parameters of a RuO 2 plane [88]. Such spinorbital pattern turns out to be the most relevant case when considering the competition between the host and the impurity as due to the AO order within the ab plane. Other possible configurations with uniform orbital order and antiferromagnetic spin pattern, e.g. G-type configuration, will be considered in the discussion throughout the manuscript. The sites i = 1 to i = 4 inside the cluster in Fig. 4 have active spin and orbital degrees of freedom while the impurity at site i = 0 has only spin degree of freedom. At the remaining sites the spin-orbital configuration is assigned and it does not change along the computation.
To determine the ground state we assume that the spin-orbital degrees of freedom are treated as classical variables. This implies that for the bonds between atoms in the host we use the Hamiltonian (14) and neglect quantum fluctuations, i.e., in the spin sector we keep only the z (Ising) components and in the orbital one only the terms which are proportional to the doublon occupation numbers (21) and to the identity operators. Similarly, for the impurity-host bonds we use the Hamiltonian Eq. (13) by keeping only the z-projections of spin operators. Since we do neglect the fluctuation of the spin amplitude it is enough to consider only the maximal and minimal values of S z i for spin S = 3/2 at the impurity sites and S = 1 in the host atoms. With these assumptions we can construct all the possible configurations by varying the spin and orbital configurations at the sites from i = 1 to i = 4 in the cluster (see Fig. 4). Note that the outer ions in the cluster belong all to the same sublattice, so two distinct cases have to be considered to probe all the configurations. Since physically it is not likely that a single impurity will change the orbital order of the host globally thus we will not compare the energies from these two cases and analyze two classes of solutions, see Sec. III B. Then, the lowest energy configuration in each class provides the optimal spin-orbital pattern for the NNs around the 3d impurity. In the case of degenerate classical states, the spin-orbital order is established by including quantum fluctuations. In the case of a periodic doping analyzed in the Sec. IV, we use a similar strategy in the computation. To use the most general formulation, we employ larger clusters having both size and shape that depend on the impurity distribution and on the spin-orbital order in the host. For this purpose, the most natural choice is to search for the minimum energy configuration in the elementary unit cell that can reproduce the full lattice by a suitable choice of the translation vectors. This is computationally doable for a periodic distribution of the impurities that is commensurate to the lattice because it yields a unit cell of relatively small size for dopings around x = 0.1. Otherwise, for the incommensurate doping the size of the unit cell can lead to a configuration space of a dimension that impedes finding of the ground state. This problem is computationally more demanding and to avoid the comparison of all the energy configurations, we have employed the Metropolis algorithm at low temperature to achieve the optimal configuration iteratively along the convergence process.

B. Two nonequivalent 3d doping cases
The single impurity problem is the key case to start with because it shows how the short-range spin-orbital correlations are modified around the 3d atom according to the host and host-impurity interactions. The analysis is performed by fixing the strength between Hund's exchange and Coulomb interaction within the host at η host = 0.1, and by allowing for a variation of the ratio between the host-impurity interaction and of Coulomb coupling at the impurity site. The choice of η host = 0.1 is made here because this value is within the physically relevant range for the case of the ruthenium materials. Small variations of η host do not affect the obtained results qualitatively. Note that the present approach is fully classical, meaning that the spins of the host and impurity are treated as Ising variables and the orbital fluctuations in the host's Hamiltonian of Eq. (14) are omitted.
As we have already discussed in the model derivation, the sign of the magnetic exchange between the impurity and the host depends on the orbital occupation of the 4d doublon around the 3d impurity. The main aspect that controls the resulting magnetic configuration is then given by the character of the doublon orbitals around the impurity, depending on whether they are active or inactive along the considered 3d − 4d bond. To explore quantitatively such competition we consider an AO order for the host with alternation of a and c doublon configurations accompanied by the C-AF spin pattern. Note that the a orbitals are active only along the b axis, while the c orbitals are active along the both axes: a and b. This state has the lowest energy for the host in a wide range of parameters for Hund's exchange, Coulomb element and crystal-field potential [88].
Due to the specific orbital pattern, the 3d impurity can substitute one of two distinct 4d sites which are considered separately below, either with a or with c orbital configurations for the doublon. Since the two 4d atoms have inequivalent surrounding orbitals, not always active along the 3d − 4d bond, we expect that the resulting ground state will have a modified spin-orbital order. Indeed, if the 3d atom replaces the 4d one with the doublon in the a orbital, then all the 4d neighboring sites have active doublons along the connecting 3d−4d bonds because they are in the c orbital configuration. On the contrary, the substitution at the 4d site with c orbital doublon configuration leads to an impurity state with neighbors having both active and inactive doublons. Since the character of the orbitals in terms of their potential to connect the doublon along the 3d − 4d bond determines the most favorable spin-exchange, we do expect a more intricate competition for the latter case of an impurity occupying the 4d site with c orbital configuration. This may lead to frustrated host-impurity interaction.

C. Doping at the site with a doublon in a orbital
We start by considering the physical situation where the 3d impurity is placed at the 4d site with the doublon occupying the a orbital. The ground state phase diagram and the sketch of the spin-orbital pattern are reported in Fig. 5 in terms of the ratio J imp /J host and the strength η imp (11) that sets the amplitude of Hund's exchange coupling with respect to the excitation energy In the FSa phase the question mark stands for that the frustrated impurity spin within the classical approach but the frustration is released by the quantum fluctuations of the NN c orbitals in the a direction resulting in small AF couplings along the a axis, and spins obey the C-AF order (small arrow). The labels FMa and AFa refer to the local spin order around the 3d impurity site with respect to the host -these states differ by spin inversion at the 3d atom site.
∆ at the 3d site. There are three different ground states that appear in the phase diagram. Taking into account the structure of the 3d − 4d spin-orbital exchange we expect that, in the regime where the host-impurity interaction is greater than that in the host, the 4d neighbors to the impurity tend to favor the spin-orbital configuration set by the 3d − 4d exchange. In this case, since the orbitals surrounding the 3d site already minimize the 3d − 4d Hamiltonian, we expect that the optimal configuration is to have the 4d spins aligned either antiferromagnetically or ferromagnetically with respect to the impurity 3d magnetic moment.
The case of AF neighbors to the 3d spin impurity occurs in the AFa phase while the FMa region is just obtained by reversing the spin at the impurity and having all the 3d − 4d bonds ferromagnetically coupled. It is interesting to note that due to the host-impurity interaction the C-AF spin-pattern of the host is modified in both the AFa and the FMa ground states. Another intermediate configuration which emerges when the host-impurity exchange is weak in the FSa phase where the impurity spin is undetermined and its configuration in the initial C-AF phase is degenerate with the one obtained after the spin-inversion operation. This is a singular physical situation because the impurity does not select a specific direction even if the surrounding host has a given spinorbital configuration. Such a degeneracy is clearly verified at the critical point η c imp 0.43 where the amplitude of the 3d − 4d coupling vanishes when the doublon occupies the active orbital. Interestingly, such a degenerate configuration is also obtained at J imp /J host < 1 when the host dominates and the spin configuration at the 4d sites around the impurity are basically determined by J host . In this case, due to the C-AF spin order, there always occur two bonds with FM character and two other having AF coupling, independently of the spin orientation at the 3d impurity. This implies that both FM or AF couplings along the 3d − 4d bonds perfectly balance each other which results in the degenerate FSv phase.
It is worth pointing out that there is a considerably large region of the phase diagram where the FSa state is stabilized and the spin-orbital order of the host is not affected by doping with the possibility of having large degeneracy in the spin configuration of the impurities. On the other hand, by inspecting the c orbitals around the impurity from the point of view of the full host's Hamiltonian Eq. (14) with orbital flips included, (τ γ+ i τ γ− j +h.c.), one can easily find out that the frustration of the impurity spin can be released by quantum orbital fluctuations. Note that the c orbitals around the impurity in the a direction have different surroundings than the ones in the b direction. The former ones are connected with two active bonds, see Fig. 6, along the b axis with orbitals a, while the latter ones are connected with only one active bond along the b axis with orbitals b. This means that in the perturbative expansion the orbital flips will add some admixture of the inactive (a) orbital to the orbitals c that neighbor the impurity in the a direction while they will not add any admixture of the inactive (b) orbital to the orbitals c that are neighboring the impurity along the b axis.
This fundamental difference can be easily included in the host-impurity bond in the mean-field manner by setting D i±b,γ = 0 for the bonds along the b axis and 0 < D i±a,γ < 1 for the bonds along the a axis. Then one can easily check that for the impurity spin pointing down we get the energy contribution from the spin-spin bonds which is given by E ↓ = α D i±a,γ , and for the impurity spin pointing upwards we have E ↑ = −α D i±a,γ , with α = α(η host ) > 0. Thus, it is clear that any admixture of the virtual orbital flips in the host's wave function polarize the impurity spin upwards so that the C-AF order of the host will be restored.

D. Doping at the site with a doublon in c orbital
Let us move to the case of the 3d atom at the 4d site where the doublon occupies the c-orbital. As anticipated above, this configuration is more intricate because the orbitals surrounding the impurity, as originated by the C-AF/AO order within the host, lead to inequivalent 3d − 4d bonds. There are two bonds where the doublon is occupying an inactive orbital (and has no hybridization with the t 2g orbitals at 3d atom), and two remaining bonds having doublon active orbitals at the 4d sites.
Since the 3d − 4d spin-orbital exchange depends on the orbital polarization on the 4d sites we do expect a competition which may modify significantly the spin-orbital correlations in the host. Indeed, in Fig. 7 the local correlations in the ground state phase diagram are shown and one observes that three configurations compete, denoted as AF1c, AF2c and FMc. In the regime where the hostimpurity exchange dominates the system tends to minimize the energy from the 3d − 4d spin-orbital coupling and, thus, the orbitals become polarized in the active configurations compatible with the C-AF/AO pattern and the host-impurity spin coupling is AF for η imp ≤ 0.43, while it is FM otherwise. This region resembles orbital polarons in doped manganites [62,65]. Also in this case, the orbital polarons minimize the total energy due to the double exchange mechanism [89].
On the contrary, for weak spin-orbital coupling between the impurity and the host there is an interesting cooperation between the 3d and 4d atoms. Since the strength of the impurity-host coupling is not sufficient to polarize the orbitals at the 4d sites, it is preferable to have an orbital rearrangement to the configuration that is inactive on the 3d − 4d bond with a spin flip. In this way the spin-orbital exchange is optimized in the host and also on the 3d − 4d bonds. The resulting state has an AF coupling between the host and the impurity as it should when all the orbitals surrounding the 3d atoms are inactive with respect to the bond direction. The modification of the orbital configuration induces the change in spin orientation. The double exchange bonds (with inactive doublon orbitals) along the b axis are then blocked and the total energy is lowered, in spite of the frustrated the spin-orbital exchange in the host. For the AF1c state the spins surrounding the impurity are aligned and antiparallel to the spin at the 3d site.
Concerning the C-AF/AO order we note that it is a modified only along the direction where the FM correlations develop and spin defects occur within the chain that is doped by the 3d atom. The FM order is broken there by the insertion of a defect with antiferromagnetically coupled spin moments. Note that this phase is purely orbital-vacancy-driven in the sense that the host develops more favorable orbital bonds to gain the energy in the absence of the orbital degrees of freedom at the impurity. At the same time the impurity-host bonds do not gen-erate too big energy losses as: (i) either η imp is small so that the loss due to E D is compensated by the gain from the superexchange ∝ J S (D 2,γ = 1) (all these bonds are AF), see Fig. 3, or (ii) J imp /J host is small meaning that the overall energy scale of the impurity-host exchange remains small. Interestingly, if we compare the AF1c with the AF2c ground states it is possible to observe that the disruption of the C-AF/AO order is anisotropic and it occurs along the C-AF chains for the AF1c phase and perpendicular to the direction of the FM chains in the AF2c configuration.
Finally, we point out that a very similar phase diagram can be obtained assuming that the host has the FM/AO order with a and b orbitals alternating from site to site. Such configuration can be stabilized by a distortion that favors the out-of-plane orbitals. In this case there is no difference in doping at one or the other sublattice. The main difference is in energy scales -for G-AF/AO order the diagram is similar to the one of Fig. 7 if we rescale J imp by half, which means that the FM/AO order is softer than the C-AF one. Note also that in the peculiar AF1c phase the impurity does not induce any changes in the host for the FM/AO ordered host. Thus we can safely conclude that the observed change in the orbital order for the C-AF host in the AF1c phase is due to the presence of the c orbitals which are not directional in the ab plane.
The opposite limit to the distortion-driven order of the host is the G-AF phase with fully polarized c orbitals, namely G-AF/FO phase which is stable if the distortion favors in-plane orbitals. From the point of view of the 3d impurity this case is trivial -the host's order is never altered and for growing η imp the impurity-host spin-spin correlation changes from AF to FM at η imp = η c imp .

A. General remarks on doping dependence
In this Section we analyze the spin-orbital patterns due to a finite concentration x of 3d impurities within the 4d host with C-AF/AO order, by assuming that the 3d atoms have a periodic arrangement. The study is performed for two representative doping distributions which are commensurate and incommensurate with respect to the lattice. The selected values are x = 1/8, x = 1/5, and x = 1/9. This choice allows us to study different regimes of competition between the spin-orbital coupling within the host and the 3d − 4d coupling. The analysis is performed as for the single impurity case, by assuming the spin and the orbital degrees of freedom are classical variables and by determining the lowest energy configuration for the total Hamiltonian of the host and the impurity. For this analysis we set the spatial distribution of the 3d atoms and we determine the optimal spin and orbital profile that minimizes the energy.
The somewhat trivial but notable problem is the periodic x = 1/2 doping where 3d and 4d atoms form a checkerboard. In this case all the bonds are the impurityhost bonds so the host's order can no longer exist. This resembles the CE phase of doped maganites where the orbital order supports the zigzag spin structure [63]. Here the ground state configuration is then entirely governed by the spin exchange constant shown in Fig. 3 meaning that the doublon at every bond will be pushed into one of the two active orbitals (the active orbitals are b and c for the a bonds and a and c for the b ones. In a case of a 2D system, which we discuss here, it is easy to satisfy all the bonds simply by putting all the doublons in the c orbitals but already adding a second layer will cause frustration, because the c bonds will tend to push the doublon into a and b orbitals, and the frustration will grow as we add more layers and move towards fully 3D system. B. C-AF phase with x = 1/8 doping We begin with the ground state diagram obtained at x = 1/8 3d doping, presented in Fig. 8. In the regime of strong impurity-host coupling the 3d − 4d spin-orbital exchange determines the orbital and spin configuration of the 4d atoms around the impurity. The most favorable state is with the doublon occupying the c-orbital at the sites which are NNs to the impurity. The spin correlations between the impurity and the host are AF (FM), if the amplitude of η imp is below (above) the critical value η c imp leading to the AFa and the FMa states. The AFa state has a striped-like profile with AF chains alternated by FM domains(consisting of three chains) along the diagonal of the square lattice. Even if in the AFa state the coupling between the impurity and the host is always AF, the overall configuration has a residual magnetic moment originating by the uncompensated spins and by the cooperation between the spin-orbital exchange in the 4d host and that for the 3d − 4d. Interestingly, at the point where the dominant 3d − 4d exchange tends to zero (i.e., for η imp η c imp ), there is a region in the phase diagram denoted FSa which is analogous to the FSa phase found for a single impurity problem, see Fig. 5. Again the impurity spin is frustrated for a purely classical approach but this frustration can be easily released by the orbital fluctuations in the host so that the C-AF order of the host is restored. This state is stable for the amplitude of η imp that is close to η c imp and when the ratio between the impurity and the host coupling, η imp (J imp /J host ), is not larger than one.
The regime of smaller J imp /J host ratio is qualitatively different -an orbital rearrangement around the impurity takes place, with a preference to have the doublons in the inactive orbitals along the 3d − 4d bonds. Such orbital configurations favor the AF spin coupling at all the 3d − 4d bonds which is stabilized by the 4d − 4d superexchange [61]. This configuration is peculiar because it generally breaks inversion and does not have any plane of mirror symmetry. It is worth pointing out that the original order in the 4d host is completely modified by the small concentration of 3d doping. The spin-orbital pattern is unaffected only in a limited region of the phase diagram. Otherwise, the AF coupling between the 3d impurity and the 4d host generally leads to patterns where FM chains are alternated with AF ones along the diagonal of the square lattice. Another relevant issue is that the cooperation between the host and impurity can lead to a fully polarized FM state. This implies that the doping can release the orbital frustration which was present in the host for the C-AF/AO order.
C. Phase diagram for periodic x = 1/5 doping Let us consider now doping concentration of x = 1/5 with a given periodic spatial profile. For this choice of 3d distribution the impurities are separated by the translation vectors u = (1, 2) and v = (2, −1) (one can show that for general periodic doping x, | u| 2 = x −1 ) so there is a mismatch between the periodicity of impurities and the host's two-sublattice AO order. The ground state diagram is reported in Fig. 9. One can note that the general structure is similar to the case at x = 1/8 with AF correlations dominating for η imp less than about η c imp and FM ones otherwise. Due to the specific doping distribution there are more phases appearing in the ground state diagram. For η imp < η c imp the most stable configuration is provided by having the impurity coupled antiferromagnetically to the host in the spin channel. This happens both in the AFv and the AFp ground states. The difference between the two AF states is due to the orbital arrangement around the impurity. For weak ratio of the impurity to the host spin-orbital exchange, J imp /J host , the orbitals around the impurity are all polarized to inactive ones, as already seen in AF1c phase of the single impurity problem. On the contrary, in the strong impurity-host coupling regime all the orbitals are polarized to be in active states around the impurity.
More generally, for all phases the approximately hyperbolic boundary given by an approximately constant value of η imp (J imp /J host ) separates the phases where the orbitals around the impurities in the c-orbital sublattice are all inactive [small η imp (J imp /J host )] from the ones where all the orbitals are active [large η imp (J imp /J host )]. The inactive orbital around the impurity always stabilize the AF coupling of the impurity spin to the host spins whereas the active orbitals can give either AF or FM exchange depending on η imp (hence η c imp , see Fig. 3). Since the doping distribution is not matching the elementary unit cell, the resulting ground states do not have specific symmetries in the spin-orbital pattern, they are generally FM due to the uncompensated magnetic moments and the impurity feels a screening by the presence of the surrounding host spins which are aligned antiparallel.
By increasing Hund's exchange coupling at the 3d site the system can stabilize a fully FM state for a large region of the ground state diagram due to the possibility of suitably polarizing the orbitals in the 4d host around the impurity. On the other hand, in the limit where the impurity-host bonds are weak, so either for η imp η c imp and J imp /J host large enough so that all orbitals around the impurity are active, or just for small J imp /J host we get the FS phases where the impurity spin at the a-orbital sublattice is undetermined in the classical approach that we use. This is a similar situation to the one found in the However, the situation here is different because due to the doping the host's order is completely altered and have became isotropic, in contrast to the initial C-AF configuration that was breaking the symmetry between the a and b direction. It was precisely this symmetry breaking that favored one impurity spin polarization over the other one. Here this mechanism is absent -one can easily check that the neighborhood of the c orbitals surrounding degenerate impurities is completely equivalent in both directions (see Fig. 10 for the view of these surroundings) so that the orbital flip argument is no longer applicable. This is a peculiar situation as it means that the spins at the impurities marked with a question mark in Fig. 9 will remain degenerate as any quantum mechanism due to short-range order is missing and the degen- eracy cannot be removed. In Fig. 10 we can see that both in FSv and FSp phases the orbitals are grouped in 3 × 3 and 2 × 2 plaquettes, respectively, that enclose the degenerate impurity spins. For FSv phase we can distingiush between two kind of plaquettes with non-zero spin polarization differing by a global spin inversion. In the case of FSp phases we have four plaquettes with zero spin polarization arranged in two pairs related by a point reflection with respect to the impurity site. It is interesting to remark that in the orbital sector these plaquettes are completely disconnected in the sense that there are no orbitally active bonds that connect them (see Fig. 6 for the pictorial definition of orbitally active bonds). This means that the purely orbital quantum effects can appear only at the short range, i.e., inside the plaquettes. However one can expect that if for some reasons the two degenerate spins in a single elementary cell will polarize then they will polarize the same in all other cells to favor long-range quantum fluctuations in the spin sector related to the translational invariance of the system. D. Phase diagram for periodic x = 1/9 doping Finally we investigate low doping concentration of x = 1/9 with a given periodic spatial profile. Here the impurities are separated by the translation vectors u = (0, 3) and v = (3, 0). Once again there is a mismatch between the periodic distribution of impurities and the host's two-sublattice AO order. The ground state diagram, shown in Fig. 11, presents gradually increasing tendency towards FM 3d−4d bonds with increasing η imp . These bonds polarize as well the 4d − 4d ones and one finds an almost FM order in the FMp state. Altogether, we have found the same phases as at the higher doping of x = 1/5, i.e., AFv and AFp at low values η imp , FMv and FMp in the regime of high η imp , separated by the regime of frustrated impurity spins in the classical approach, visible within the phases: FSv, FS1p, and FS2p.
The difference between the two AF (FM) states is due to the orbital arrangement around the impurity. As for other considered doping levels, x = 1/8 and x = 1/5, we find neutral (inactive) orbitals around 3d impurities in the regime of low J imp /J host in AFv and FMv phases which lead to spin defects within the 1D FM chains in the C-AF spin order. We have reported a similar behavior for single impurities in the low doping regime in Sec. III. This changes radically above the orbital transition withing these two types of local magnetic order, where the orbitals reorient into the active ones. One finds that spin orientations are then the same as those of their neighboring 4d atoms, with some similarities to those found at x = 1/5 and shown in Fig. 9.
In the crossover regime between the AF and FM local order induced by the impurities, frustrated spins occur. This follows from the the local configurations around them, which include two ↑-spins and two ↓-spins accompanied by c orbitals at the NN 4d sites. This frustration may be removed by quantum fluctuations and we suggest that this happens again, in the same way as for x = 1/8 doping, as indicated by small arrows in the respective FS phases shown in Fig. 11.

V. SUMMARY AND CONCLUSIONS
We have derived the spin-orbital superexchange model for 3d 3 impurities in the 4d host in the regime of Mott isulating phase. Although the impurity has no orbital degree of freedom, we have shown that it contributes to the orbital physics and influences strongly the orbital order and tends to project out the inactive orbital at the impurity-host bond to maximize the energy gain from virtual charge fluctuations. In this case the interaction along the superexchange bond can be either antiferromagnetic or ferromagnetic, depending on the ratio of Hund's exchange coupling at impurity (J H 1 ) and host (J H 2 ) ions and on the mismatch ∆ between the 3d and 4d atom energy levels modified by the difference in Hubbard U 's and Hund's exchange J H 's at both atoms. This ratio, denoted here as η imp (11), replaces here the conventional ratio η = J H /U often found in the spin-orbital superexchange models of undoped compounds (e.g., in the Kugel-Khomskii model for KCuF 3 [52]) where it quantifies the proximity to ferromagnetism. On the other hand, if the overall coupling between host and impurity is weak in the sense of the total superexchange, J imp , with respect to the host value, J host , the orbitals being next to the impurity may be forced to stay inactive which modifies the magnetic properties -in such cases the impurity-host bond is always antiferromagnetic.
As we have seen in the case of a single impurity, the above two mechanisms can have a nontrivial effect on the host, especially if the host itself is characterized by frustrated interactions, as it happens in the parameter regime where the C-AF phase is stable. For this reason we have focused mostly on the latter phase of the host and we have presented the phase diagrams of a single impurity configuration in the case when the impurity is doped on the sublattice where the orbitals form a checkerboard pattern with alternating c and a orbitals occupied by doublons. The diagram for the c-sublattice doping shows that in some sense the impurity is never weak, because even for a very small value of J imp /J host it can release the host's frustration around the impurity site acting as an orbital vacancy. On the other hand, for the a-sublattice doping when the impurity-host coupling is weak, i.e., either J imp /J host is weak or η imp is close to η c imp , we have observed an interesting quantum mechanism releasing frustration of the impurity spin that cannot be avoided in the purely classical approach. It turned out that in such situations the orbital flips in the host make the impurity spin polarize in such a way that the C-AF order of the host is completely restored.
The cases of the periodic doping studied in this paper show that the host's order can be completely altered already for rather low doping of x = 1/8, even if the J imp /J host is small. In this case we can stabilize a ferrimagnetic type of phase with a four-site unit cell having magnetization S z i = 3/2. We have established that the only parameter range where the host's order remains unchanged is when η imp is close to η c imp and J imp /J host 1. The latter value is very surprising as it means that the impurity-host coupling must be large enough to keep the host's order unchanged -this is another manifestation of the orbital vacancy mechanism that we have already observed for a single impurity. Also in this case the impurity spins are fixed with help of orbital flips in the host that lift the degeneracy which arises in the classical approach. We would like to point out that the quantum mechanism that lifts the ground state degeneracy mentioned above and the role of quantum fluctuations is of particular interest for the periodically doped checkerboard systems with x = 1/2 doping which is an interesting problem for future research.
From the point of view of generic, i.e., non-periodic doping, the most representative cases are those of an incommensurate doping. To uncover the generic rules such cases, we have studied the periodic x = 1/5 and x = 1/9 dopings. One finds that when the period of the impurity positions does not match the period of 2 for both the spin and orbital order of the host, interesting novel types of order emerge. In such cases the elementary cell must be doubled in both lattice directions which clearly gives a chance of realizing more phases than in the case of commensurate doping. Our results show that indeed, the number of phases increases from 4 to 7 and the host's order is altered in each of them. Quite surprisingly, the overall character of the phase diagram remained unchanged with respect to the one for x = 1/8 doping and if we ignore the differences in configuration it seems that only some of the phases got divided into two versions, differing either by the spin bonds' polarizations around impurities (phases around η c imp ), or by the character of the orbitals around the impurities [phases with inactive orbitals for small enough η imp (J imp /J host ) versus phases with active ones for larger values of η imp (J imp /J host )]. Orbital polarization in this latter region resembles orbital polarons in doped manganites -also here such states are stabilized by the double exchange mechanism [89].
A closer inspection of underlying phases reveals however a very interesting degeneracy of the impurity spins at x = 1/5 that arises again from the classical approch but this time it cannot be released by short-range orbital flips. This happens because the host's order is already so strongly altered that it is no longer anisotropic (as it was the case of the C-AF phase) and there is no way to restore the orbital anisotropy around the impurities that could lead to spin-bonds imbalance and polarize the spin. On the other hand the case of x = 1/9 doping where such effect is absent, and the impurity spins are always polarized as it happens for x = 1/8, shows that this is rather a peculiarity of the x = 1/5 periodic doping. Indeed, one can easily notice that for x = 1/5 every atom of the host is a nearest neighbor of some impurity. In contrast, for x = 1/8 we can find three host's atoms per unit cell which do not neighbor any impurity and for x = 1/9 there are sixteen of them. For this reason the impurity effects are amplified for x = 1/5 which is not unexpected although one may find somewhat surprising that the ground state diagrams for the lowest and the highest doping considered here are very similar. This suggests that the cooperative effects of multiple impurities are indeed not very strong in the low-doping regime, so the diagram obtained for x = 1/9 can be regarded as generic for the dilute doping regime with uniform spatial profile.
In summary, we have shown that 3d 3 impurities change radically the spin-orbital order around them in the 4d 4 host, independently of the parameter regime. As a general feature we have found that doped 3d 3 ions within the 4d 4 correlated insulator have frustrated spins and polarize the orbitals of the host when the spin and Hund's exchange at the impurity are sufficiently large. This trend is independent of doping and is expected to lead to the global changes of spin-orbital order in doped materials. While the latter effect is robust, we argue that the longrange spin fluctuations resulting from the translational invariance of the system will likely prevent the ground state from being macroscopically degenerate, so if the impurity spins in one unit cell happens to choose its polarization than the others will follow. On the contrary, in the regime of weak Hund's exchange 3d 3 ions act not only as spin defects which order antiferromagnetically with respect to their neighbors, but also induce at them the doublon occupancy in inactive orbitals. This behavior with switching between inactive and active orbitals by an orbitally neutral impurity may lead to multiple interesting phenomena at macroscopic doping when global modifications of the spin-orbital order are expected to occur. Finally, we note that this generic treatment and the general questions addressed here, such as the relaese of frustration for competing spin structures due to periodic impurities, are relevant to double perovskites [90]. a given direction γ. Differently from the first term, the