text:
# Unconventional superconducting states induced in a ferromagnet by a 𝑑-wave superconductor

## Abstract

We develop a quasi-classical theory for the superconducting proximity effect in a ballistic ferromagnetic layer in contact with a $`d`$-wave superconductor. In agreement with recent experiments we find that the density of states oscillates around the normal state value with varying the thickness of the ferromagnetic layer. We show that the phase, the amplitude, and the period of these oscillations depend on the orientation of the superconductor. This effect reveals spatial oscillations and anisotropy of the induced superconducting correlations in the ferromagnet.

Proximity effects in hybrid structures of superconductors and ferromagnets provide the possibility for the controlled studies of the coexistence of ferromagnetism and superconductivity. Superconducting correlations penetrate into normal metals by means of a special scattering process at the normal metal-superconductor (NS) interface known as Andreev reflection Andreev : an electron propagating in the normal metal can cross the NS boundary by reflecting as a hole and transferring a Copper pair into the superconductor. The Andreev reflected hole is correlated with the incident electron and carries information about the phase of the order parameter of the superconductor. The existence of correlated electron-hole pairs corresponds to a non-vanishing superconducting pair amplitude (order parameter), which gives to the normal metal the typical characteristics of the superconducting states Tink .

The peculiarity of this proximity effect in ferromagnetic metals comes from the fact that the Andreev reflection from a singlet pairing superconductor is accompanied by an inversion of the spin direction, and consequently the correlated electron-hole occupy opposite spin bands. The spin splitting exchange field of the ferromagnet causes a phase shift of the correlated Andreev electron-hole pairs, which results in a spatially oscillatory behavior of the induced pair amplitude paos . Such an inhomogeneous superconducting state coexisting with ferromagnetism is strikingly similar to an FFLO state in a bulk superconductor with spin splitting originally predicted by Fulde and Ferrell and by Larkin and Ovchinnikov FFLO .

One of the most interesting manifestations of the oscillations in the proximity pair amplitude is the oscillatory behavior of the density of states (DOS) in ferromagnet-superconductor bilayers as function of the ferromagnetic layer thickness. Kontos et al aprili01  have observed this effect experimentally in thin ferromagnetic layers connected to a conventional $`s`$-wave superconductor. Theoretically, the superconducting proximity effect in ferromagnetic layers has been studied extensively paos ; zbn02 . It has been shown zbn01  that the theory based on the quasi-classical formalism Eil  is in a quantitative agreement with the experiments.

In parallel, there also has been much attention to the proximity of a ferromagnet to high-$`T_c`$ superconductors, which are widely believed to have a dominant $`d`$-wave pairing symmetry tsuei00 . Most of the studies have concentrated on the influence of the exchange field on the zero bias conductance peak (ZBCP) as the most remarkable feature in the tunneling spectroscopy of a junction between a high-$`T_c`$ superconductor and a normal metal hu94 ; tanaka95 . ZBCP is the result of a zero-energy Andreev peak (ZEAP) in the DOS, which originates from change in the sign of the $`d`$-wave order parameter under a $`\pi /2`$ rotation. Many effects including a splitting of the ZBCP and the developing of a zero bias conductance dip (ZBCD) have been found and attributed to the exchange splitting larkin98 ; zhu00 ; stef02 .

Most recently, Freamat and Ng Ng03  have performed tunneling spectroscopy measurement on multilayered junctions of $`d`$-wave superconductor and ferromagnets. They observed an oscillatory behavior in the tunneling conductance spectra with the thickness of the ferromagnetic domain. This effect has been considered as the $`d`$-wave analog of the effect observed in aprili01  and explained as the signature for an FFLO state induced in the ferromagnetic layers by the $`d`$-wave superconducting layer. However one would expect that the anisotropy and the sign change of the $`d`$-wave superconducting order parameter leads to an anisotropy and phase effects in the oscillation of the DOS, which distinguish it from the case of a conventional $`s`$-wave superconductor. The goal of the present work is to study this effect theoretically.

In this letter we study the proximity effect at the interface between a ferromagnetic metal and a $`d`$-wave superconductor. We model the ferromagnet as a thin ballistic film with rough boundaries which is connected to the superconductor through a disordered interface of finite transparency. We make use of the quasi-classical formalism Eil  to calculate the density of states in the ferromagnetic layer. In correspondence with the recent experiments, we find that at sufficiently strong exchange field the DOS oscillates as a function of the thickness of the ferromagnetic layer. We show that the phase, the amplitude and the period of the DOS oscillations are affected by changing the orientation of the $`d`$-wave order parameter. This effect, which reveals the highly anisotropic nature of the induced superconducting correlations in the ferromagnet, can also be used as a further test of the $`d`$-wave scenario in high-$`T_c`$ superconductors.

We consider the proximity system shown in Fig. 1. A thin ferromagnetic layer (F) of thickness $`d`$ is contacted by a $`d`$-wave superconductor (S) on one side and covered on the other side by an insulator. We consider a clean structure at both of the S and F sides of the FS-interface. At the same time the FS-interface itself may contain disorder and band mismatch which enhance backscattering of quasi-particles from this boundary. F is characterized by a mean exchange field $`h`$ which is assumed to be constant in the F-layer and vanishes in S. S is characterized by a $`d_{x^2y^2}`$-wave order parameter of the form $`\mathrm{\Delta }_d(\theta )=\mathrm{\Delta }\mathrm{cos}[2(\theta \chi )]`$, where $`\theta `$ is the angle that the Fermi velocity $`𝐯_\text{F}`$ makes with the normal to the interface and $`\chi `$ is the angle between the crystallographic $`a`$ axis of S and the interface normal. For calculation we adopt the Eilenberger equations Eil  in the clean limit for this system. In the absence of spin-flip scattering in F, transport of quasi-particles with spin $`\sigma `$ ($`=\pm 1`$) is described by an independent equation for the corresponding matrix Green’s function $`\widehat{g}_\sigma (E,𝐯_\text{F},𝐫)`$, which reads

$`i𝐯_\text{F}\widehat{g}_\sigma (E,𝐯_\text{F},𝐫)=[(E+\sigma h(𝐫))\widehat{\tau }_3`$
$`i\widehat{\tau }_2\mathrm{\Delta }(𝐫,𝐯_\text{F}),\widehat{g}_\sigma (E,𝐯_\text{F},𝐫)].`$ (1)
Here $`\widehat{\tau }_i`$ denote the Pauli matrices and $`\mathrm{\Delta }(𝐫,𝐯_\text{F})`$ is the superconducting order parameter, which is taken to be real. We neglect the spatial variation of the order parameter close to the FS-interface and take $`\mathrm{\Delta }(𝐫,𝐯_\text{F})=\mathrm{\Delta }_d(𝐯_\text{F})`$ in all points of S.

Equation (1) can be conveniently solved along a classical trajectory. A typical trajectory is shown in Fig. 1. It starts at S and extends into the F layer making an angle $`\theta _1`$ with the normal to the FS-interface. After several reflections from the insulator and the FS interface, which makes a path of total length $`l`$ in F, the trajectory returns into S at an angle $`\theta _2`$. The superconducting order parameter that an electron feels in the beginning and at the end of the trajectory are respectively $`\mathrm{\Delta }_d(\theta _1)`$ and $`\mathrm{\Delta }_d(\theta _2)`$. We have solved Eq. (1) and found that the DOS on a given trajectory only depends on $`\theta _1`$, $`\theta _2`$, and $`l`$ and given by

$`N(E,\theta _1,\theta _2,l)={\displaystyle \frac{N_0}{4}}{\displaystyle \underset{\sigma =\pm 1}{}}\text{Re}[\text{tr}\widehat{\tau }_3\widehat{g}_\sigma (E,𝐯_\text{F},𝐫)]={\displaystyle \frac{N_0}{2}}{\displaystyle \underset{\sigma =\pm 1}{}}\text{Re}\{i\mathrm{tan}[k_\sigma l/v_F+\alpha (\theta _1)+\alpha (\theta _2)+\eta (\theta _1,\theta _2)/2]\},`$ (2)
where $`N_0`$ is DOS at the Fermi level in the normal state, $`\alpha (\theta )=\mathrm{arcsin}[E/|\mathrm{\Delta }_d(\theta )|]/2`$, $`k_\sigma =(E+\sigma h)/v_F`$ and $`\eta (\theta _1,\theta _2)=\pi [1\text{sign}(\mathrm{\Delta }_d(\theta _1)\mathrm{\Delta }_d(\theta _2))]/2`$ is the $`\pi `$ phase shift resulting from the possible change in the sign of the order parameter at the beginning and the end of the trajectory.

The total density of states $`N(E)`$ is obtained by averaging the DOS per trajectory Eq. (2) over all the possible values of $`\theta _1`$,$`\theta _2`$, and $`l`$:

$$N(E)=𝑑l𝑑\theta _1𝑑\theta _2\overline{p}(l,\theta _1,\theta _2)N(E,l,\theta _1,\theta _2),$$
(3)
where $`\overline{p}(l,\theta _1,\theta _2)`$ is the distribution function of the trajectories.

To determine $`\overline{p}(l,\theta _1,\theta _2)`$, we model F as a weakly disordered thin layer bounded by a rough surface and a rough FS-interface with a constant transparency $`t`$. We introduce $`g(\theta _1,\theta _2)`$ as the correlation function between the incoming and outgoing directions $`\theta _1`$ and $`\theta _2`$. This function determines the specularity of the scattering from the boundaries. For perfectly flat boundaries with purely specular scattering, $`\theta _1`$ and $`\theta _2`$ are completely correlated and $`g\delta (\theta _1+\theta _2)`$. For the boundaries with strong roughness the dominantly diffusive scattering destroys any correlation between $`\theta _1`$ and $`\theta _2`$, and $`g(\theta _1,\theta _2)`$ is constant. In the general case of an arbitrary strength of the roughness we assume $`g=C^1\text{exp}((\theta _1+\theta _2)^2/(\pi z)^2)`$, where $`C=_{\pi /2}^{\pi /2}𝑑\theta _1𝑑\theta _2g(\theta _1,\theta _2)`$ is a normalization factor and $`z`$ measures the degree of the roughness. The limits of $`z1`$ and $`z1`$ correspond to weak and strong roughness respectively. As we will see below the parameter $`z`$ play a crucial role to reveal the anisotropic nature of the induced superconducting correlations in the F-layer.

The full distribution function $`\overline{p}(l,\theta _1,\theta _2)`$ can be written in terms of $`g(\theta _1,\theta _2)`$ as

$`\overline{p}(l,\theta _1,\theta _2)=g(\theta _1,\theta _2)[t\delta (l{\displaystyle \frac{d}{\mathrm{cos}\theta _1}}{\displaystyle \frac{d}{\mathrm{cos}\theta _2}})`$
$`+rp(l{\displaystyle \frac{d}{\mathrm{cos}\theta _1}}{\displaystyle \frac{d}{\mathrm{cos}\theta _2}})],`$ (4)
where $`p(l)`$ is the length distribution of the trajectories regardless of the values of $`\theta _1`$ and $`\theta _2`$, and $`r=1t`$. For $`t=1`$ the delta function in Eq. (4) assure that the total length of the trajectory is given by $`l=d/\mathrm{cos}\theta _1+d/\mathrm{cos}\theta _2`$. In the limit of small transparency of the FS-interface $`t1`$, $`p(l)`$ can be approximated by the exponentially decaying function zbn02

$$p(l)=\frac{1}{\overline{l}}e^{l/\overline{l}}\mathrm{\Theta }(l2d),$$
(5)
where $`\overline{l}2d\mathrm{ln}(\mathrm{}_{\text{i}mp}/d)/t=2d_t`$ is the mean trajectory length in the limit of small $`d/\mathrm{}_{\text{i}mp}`$ ($`\mathrm{}_{\text{i}mp}`$ is the elastic mean free path).

Combining Eqs. (2), (3), (4), and (5) we find the following result for the total DOS

$$N(E)=\frac{N_0}{2}\underset{\sigma =\pm 1}{}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}[t+rp(2nk_\sigma d)]J_n,$$
(6)
where $`p(k)=e^{i2kd}/(ik\overline{l}+1)`$ is the Fourier transform of $`p(l)`$ and

$`J_n(k)={\displaystyle _{\pi /2}^{\pi /2}}𝑑\theta _1𝑑\theta _2g(\theta _1,\theta _2)f_n(\theta _1)f_n(\theta _2)e^{in\eta (\theta _1,\theta _2)},`$
$`f_n(\theta )=\{\begin{array}{cc}e^{i2nk_\sigma d/\mathrm{cos}\theta +i2n\varphi (\theta )}\text{for}\hfill & \hfill |E||\mathrm{\Delta }_d(\theta )|\\ e^{i2nk_\sigma d/\mathrm{cos}\theta 2|n|\beta (\theta )}\text{for}\hfill & \hfill |E|>|\mathrm{\Delta }_d(\theta )|\end{array}.`$ (9)
Here $`\varphi (\theta )=\mathrm{arccos}[E/|\mathrm{\Delta }_d(\theta )|]/2`$ is the Andreev phase and $`\beta (\theta )=\text{acosh}[E/|\mathrm{\Delta }_d(\theta )|]/2`$.

Equation (6) expresses the total DOS in terms of known functions. In real samples the F-layer is expected to have a nonuniform thickness due to the large scale roughness of the boundaries. To take this into account we average (6) over a Gaussian distribution of the thicknesses around a mean value $`d`$. This leads to a smearing of the sharp features in the DOS which arise from the lower cutoff in the length distribution. The qualitative behavior, however, will not be changed. In practice, we have taken the width of the distribution to be of order 10%, which corresponds to the conditions of the experiments aprili01 ; Ng03 .

Let us start by analyzing the effect of a weak exchange field on the proximity DOS of a thin F-layer with $`dv_\text{F}/\mathrm{\Delta }`$ and $`t1`$. For the orientation $`\chi =\pi /4`$ the main feature in the DOS of a normal layer ($`h=0`$) is a sharp ZEAP for any strength of the roughness $`z`$. In Fig. 2a the normalized DOS $`N(E)/N_0`$ versus $`d_th/v_\text{F}`$ is plotted for different energies and values $`\chi =\pi /4`$, $`z=0.5`$. The effect of the exchange field is to split the ZEAP into distinct peaks at finite energies. Increasing $`h`$ further leads to the shifting of the splitted peaks toward higher energies, which is associated with decreasing the height of these peaks. The DOS at $`E=0`$ decreases with $`h`$ and goes through a minimum as $`d_th/v_\text{F}`$ is increasing. Thus the exchange field can suppress the ZEAP in the DOS and develop a dip at $`E=0`$. This will lead to a ZBCD in the conductance spectra at low temperatures, which also has been observed in the experiments larkin98 .

For the orientation $`\chi =0`$ the ZEAP is absent for perfectly specular boundaries ($`z=0`$). The ZEAP, however, appears when $`z`$ is finite. This is the result of diffusive scattering caused by the roughness at the boundaries, which allows for the sign change of the superconducting order parameter at two sides of a sizable fraction of trajectories even when $`\chi =0`$ fogel97 . In Fig. 2b we plot the normalized DOS $`N(0)/N_0`$ at $`E=0`$ versus $`d_th/v_\text{F}`$ when $`\chi =0`$ and for different values of $`z`$. For the normal layer ($`h=0`$) the height of the ZEAP increases with increasing $`z`$. The effect of the exchange field for a given value of $`z`$ is to suppress the superconducting features (zero DOS for $`z=0`$ and ZEAP for finite $`z`$ ) in the DOS. At higher exchange fields when $`d_th/v_\text{F}\genfrac{}{}{0pt}{}{>}{}1`$, the zero energy DOS for all values of $`z`$ approach each other and becomes close to the normal state values $`N_0`$.

At strong exchange fields when $`dh/v_\text{F}1`$, the DOS at all energies oscillates around the normal state values $`N_0`$ as a function of $`dh/v_\text{F}`$. This oscillatory behavior of induced superconducting correlations has been observed in the experiments Ng03 . Fig. 3 represents these oscillations resulting from our calculations. The energy dependence of the normalized DOS $`N(E)/N_0`$ is plotted for different values of the F-layer thickness when the FS-interface is highly transparent ($`t=1`$), $`h=3\mathrm{\Delta }`$ and $`\chi =\pi /4`$. We have taken $`z=1`$, modeling rough boundaries. In this case the period of the DOS oscillations is roughly $`\pi /4`$.

Let us now analyze the effect of the $`d`$-wave symmetry of the superconducting order parameter on the oscillations of the DOS in F. For this we examine the oscillatory behavior of the DOS at the Fermi level ($`E=0`$) for different orientations of S. We consider the case of $`t1`$ where for $`hv_\text{F}/d`$ the deviations of the DOS from the normal state value $`N_0`$ are of order $`t`$. In Fig. 4a we have plotted the reduced DOS $`[N(0)/N_01]/t`$ versus $`dh/v_\text{F}`$ for the case of specular boundaries ($`z=0`$) at different $`\chi `$. The variation of $`\chi `$ affects the phase and the amplitude of the oscillations. Varying $`\chi `$ from $`0`$ to $`\pi /4`$, the phase is shifted by $`\pi `$ and the amplitude is decreased. The $`\pi `$ phase shift originates from the sign change of the order parameter, which occurs for all the trajectories in specular boundaries if $`\chi =\pi /4`$. We note that the period of the oscillations is almost the same for different $`\chi `$ and equals to $`\pi /2`$, the period for the case where S has a $`s`$-wave order parameter zbn01 .

The effect of an anisotropic order parameter is more pronounced for the case of diffusive boundaries. This is shown in Fig. 4b where the DOS oscillations at $`E=0`$ are presented for the case of $`z=1`$. Now, in addition to the variation of the phase and the amplitude, by changing the S orientation from $`\chi =0`$ to $`\chi =\pi /4`$ the period of the oscillations is also changed. Fig. 4c shows how the change in the period occurs when the strength of the roughness $`z`$ varies for $`\chi =\pi /4`$. We can understand these effects by noting that the spatially oscillating order parameter in the F-layer has a direction-dependent amplitude and phase resulting from the anisotropy and the sign change of the $`d`$-wave order parameter of S. The change in the phase, amplitude and the period of the DOS oscillations is the result of unconventional induced superconducting correlations

In conclusion, we have investigated theoretically the superconducting proximity effect in a thin ferromagnetic layer in contact with a $`d`$-wave superconductor. In correspondence with recent experiments Ng03  we have found that at sufficiently strong exchange fields the density of states in the ferromagnet oscillates around the normal state value as a function of its thickness. The phase, the amplitude, and the period of the oscillations depend on the orientation of the superconductor. This direction-dependence is the signature of an unconventional oscillatory superconducting state induced in the ferromagnet by the proximity to the $`d`$-wave superconductor.

We thank M. Freamat, P. Fulde, Yu. V. Nazarov, and H. Schomerus for useful discussions.


format: text

text:
# Energetics of clusters in the two-dimensional Ising spin glass

## I Introduction

According to one of the principal scenarios for the spin glass state, the “droplet picture”Fisher and Huse (1986, 1987, 1988); Bray and Moore (1986), the minimum energy excitation of linear dimension $`L`$ containing a given spin, a droplet, has an energy proportional to $`L^\theta `$. The “stiffness exponent” $`\theta `$ is positive if one assumes that the transition temperature $`T_c`$ is finite. It follows that, in the thermodynamic limit, excitations that flip a finite fraction of the spins cost an infinite energy. However, the results of several numerical calculationsKrzakala and Martin (2000); Palassini and Young (2000); Marinari and Parisi (2000); Katzgraber et al. (2001); Katzgraber and Young (2001, 2002); Houdayer et al. (2000); Houdayer and Martin (2000), on small system sizes, imply that the amount of energy needed to generate system-size droplet excitations is independent of size, so that $`\theta `$ for droplet-like excitations is zero. This is consistent with the alternative replica symmetry breaking (RSB)Parisi (1979, 1980, 1983); Mézard et al. (1987) picture. Note that numerics simultaneously finds that the stiffness exponent for domain-wall excitations is positive and therefore apparently different from the exponent for droplets, an unexpected feature.

To explain the discrepancy between those two different estimates of the stiffness exponent, it has been proposedMoore (2002); Hartmann and Moore (2002) that the droplet theory is correct on large scales, i.e. the stiffness exponent $`\theta `$ is the same for droplets and domain walls, and that the apparently contradictory numerical data can be explained by corrections to scaling. More precisely Refs. Moore, 2002; Hartmann and Moore, 2002 propose that the droplet energy $`\mathrm{\Delta }E`$ has the form

$$\mathrm{\Delta }E=AL^\theta +BL^\omega ,$$
(1)
where $`\theta >0`$ (if $`T_c>0`$ which is the case for $`d3`$), $`\omega `$ is a correction to scaling exponent, and $`A`$ and $`B`$ are positive constants. If $`\theta >0`$, $`\mathrm{\Delta }E`$ has a minimum at some value of $`L`$ and if this occurs in the range of sizes where numerical data is taken, the numerical values could be roughly independent of size in the range studied. In general, only a modest range of sizes, $`L10`$, can be studied in dimension $`d`$ equal to 3.

To quantitatively check this proposal, simulations have recently been carried outKatzgraber and Young (2003) on a one-dimensional model with long range interactions which fall off as a power of the distance, in a region of parameters for which $`T_c>0`$. A much larger range of sizes can be studied than in three dimensions but even up to $`L=512`$ the energy of system-size excitations was found to be independent of size. Hence, for Eq. (1) to be correct with $`\theta >0`$ the values of $`\theta `$, $`\omega `$, $`A`$ and $`B`$ must conspire to give an almost constant $`\mathrm{\Delta }E`$ over a very wide range of sizes.

Another model where large sizes can be studied is the two-dimensional Ising spin glass at $`T=0`$, where efficient algorithmsBieche et al. (1980); Hartmann and Rieger (2001); Rieger et al. (1996) can be applied to determine exact ground states. Now $`\theta <0`$ in $`d=2`$, with $`T=0`$ “domain-wall” calculations givingMcMillan (1984); Bray and Moore (1984); Rieger et al. (1996); Palassini and Young (1999); Hartmann and Young (2001); Carter et al. (2002) $`\theta 0.28`$. However, finite-$`T`$ simulationsLiang (1992); Kawashima et al. (1992) and $`T=0`$ studies in a magnetic fieldKawashima and Suzuki (1992), both of which excite droplets rather than domain walls, find $`\theta 0.47`$. A good review of this situation is given in Ref. Kawashima and Aoki, 2000. Since $`\theta <0`$ in $`d=2`$, both the $`L^\theta `$ and $`L^\omega `$ contributions to the droplet energy in Eq. (1) decrease with increasing $`L`$, but Refs. Moore, 2002; Hartmann and Moore, 2002 argue that they could combine to give an effective exponent of about $`0.47`$ for small $`L`$ crossing over to the asymptotic value of about $`0.28`$ for larger sizes. Hartmann and MooreHartmann and Moore (2002) find evidence for this crossover in their numerical studies in which a certain prescription was used to generate droplets, though earlier Hartmann and one of usHartmann and Young (2002), who generated droplets in a different way, did not find evidence for an effective exponent close to $`0.47`$. This difference may indicate that the amplitude of the correction term $`B`$ in Eq. (1) depends significantly on the precise way in which the droplets are generated. In yet another way of generating droplets, Picco et al.Picco et al. (2002a, b) find $`\theta `$ close to the $`0.47`$ value but with no crossover to the domain wall value of around $`0.28`$ for the larger sizes.

Recently, in an interesting paper, Lamarcq et al.Lamarcq et al. (2002) have calculated the energy of droplet (cluster)clu  excitations in three dimensions, in a manner very similar to the spirit of the droplet modelFisher and Huse (1986, 1987, 1988). Interestingly they find that the energy actually decreases, although very slightly, with increasing size. They also find that the volume of the droplets has a non trivial fractal dimension less than the space dimension. In Ref. Moore, 2002, this slow decrease of the droplet energy was taken as a possible evidence of the relevance of the correction to scaling term in Eq. (1). Note however that no sign of a crossover towards the supposedly correct value of the stiffness exponent was reported, even when larger clusters could later be included in the analysisLamarcq et al. (2003).

Here we perform a calculation similar to that of Lamarcq et al. but in two dimensions. The aim is to see if the energies of the droplet energies calculated in this way fit Eq. (1) and give an effective exponent close to $`0.47`$ at small sizes crossing over to about $`0.28`$ at larger sizes as has been proposed by MooreMoore (2002) and Hartmann and MooreHartmann and Moore (2002). Since the way our clusters are generated is directly inspired by the original definition of the droplets, one could expect the prediction (1) to be well suited in our case. Instead, we do not find any crossover in $`\mathrm{\Delta }E`$ by applying the definition of Lamarcq et al. in 2$`d`$ for clusters as large as 64 spins, neither do we find a very negative exponent $`0.47`$. We find however an exponent similar to $`0.47`$ at very small sizes if a slight modification of the definition of the clusters is made, namely if the size of the clusters is let free to evolve between $`n`$ and $`2n`$, as suggested by the original “scale invariant” definition of the droplets.

## II The Model and Details of the Numerics

We take the standard Edwards-Anderson spin glass model

$$=\underset{i,j}{}J_{ij}S_iS_j,$$
(2)
where the $`S_i=\pm 1`$ are Ising spin variables at the sites of a simple cubic lattice, and the $`J_{ij}`$ are nearest neighbor interactions with a Gaussian distribution with zero mean and standard deviation unity. Periodic boundary conditions are applied on lattices with $`N=L^2`$ spins. For most of our work we take $`L=64`$ but we also did some calculations with smaller sizes down to $`L=16`$. We use $`N_s=1000`$ realizations of the disorder.

We are interested in the properties of low-energy, droplet-like excitations above the ground state of the Hamiltonian in Eq. (2). In order to determine these, we first find the ground state for each realization of the disorder. This is done using either the Köln spin glass serverspi , or, for small sizes, by parallel tempering Hukushima and Nemoto (1996); Marinari (1998).

Following Ref. Lamarcq et al., 2002, we then generate a droplet by first choosing randomly a “central” spin in the system and reversing it. We then construct a connected cluster of fixed size $`n`$ around this central spin by flipping all the spins in this cluster. We then find the new ground state with the following three constraints:

1. The central spin is always flipped with respect to the ground state,
2. The number of spins in the cluster is always $`n`$,
3. The cluster is always connected.

The new ground state is found by a combination of parallel tempering and a Kawasaki-type dynamics for the spins in the cluster, in order to conserve its size $`n`$ constant. Our algorithm does not satisfy detailed balance, but this is irrelevant since we only want to find the new, constrained ground state. The lowest temperature in the parallel tempering is taken to be sufficiently low, $`T=0.02`$, that the system is always in a local minimum with no random noise kicking the system to higher energy states. Our results for ground states in the presence of a droplet are obtained from spin configurations at this lowest temperature. We study cluster sizes up to $`n=64`$.

In this procedure, which follows Ref. Lamarcq et al., 2002, the size of the cluster is strictly fixed. It is therefore slightly different from the droplets of Fisher and HuseFisher and Huse (1986, 1987, 1988) which are defined in a “scale-invariant” way, i.e. the size is not strictly fixed but allowed to vary over a certain scale. We have therefore also computed the energies of scale-invariant excitations in which the size of the cluster is allowed to vary by a factor of 2, more precisely to lie between $`n/2`$ and $`n1`$ for $`n=4`$, $`8`$, $`\mathrm{}`$, $`64`$. Examples of some randomly chosen clusters of size $`n=64`$ are shown in Fig. 1.

Our criterion for ensuring that we have found the true ground state in the presence of an overturned cluster is as follows. We first make a few (typically 10) extremely long runs with different bond configurations to estimate the typical time scale to reach the new ground state, $`t_{\mathrm{typ}}(n)`$. We record the energy versus time of the lowest temperature for the two copies, for simulations of at least $`10t_{\mathrm{typ}}`$ sweeps. From these graphs, we can evaluate $`t_{\mathrm{typ}}`$ as the typical time it takes those 10 realizations of the disorder to give a result that can be confidently taken as the ground state. Then, for each of the $`N_s`$ samples, we first run two copies of the system with independently drawn initial clusters for $`t_{\mathrm{typ}}(n)`$ Monte Carlo steps. We then require that the two copies have the same energy for a further continuous period of time equal to $`t_{\mathrm{typ}}(n)`$. In this way, we spend more time on the “hard” samples, which need several times $`t_{\mathrm{typ}}(n)`$ to converge, than on the easy samples for which we do not need to run much more than $`t_{\mathrm{typ}}(n)`$ sweeps. Obviously, one can not be absolutely sure that the ground state has really been found since it might happen that both copies stay for a time larger than $`t_{\mathrm{typ}}(n)`$ in the same state which is not the ground state. However, this seems to be fairly unlikely. To be on the safe side we used $`t_{\mathrm{typ}}(n=64)`$, for all sizes, even though, for $`n<64`$, one has $`t_{\mathrm{typ}}(n)t_{\mathrm{typ}}(n=64)`$.

For the scale-invariant clusters we start one copy with the minimum-size cluster, $`n/2`$, and the other with the maximum-size cluster, $`n1`$, to minimize further the probability that the two copies inadvertently spend a long time in the same state which is not the ground state.

We now present our numerical results.

## III results

### III.1 Clusters of Fixed Size

In this section we present results of the simulations where we fix the droplet size precisely, as was done by Lamarcq et al.Lamarcq et al. (2002) on the $`3d`$ model.

#### III.1.1 Droplet Energies

Fig. 2 shows the droplet energy, $`\mathrm{\Delta }E`$, where the average is over $`N_s=1000`$ samples, as a function of $`n`$, the number of spins in the droplet, on a log-log scale. For each sample, $`\mathrm{\Delta }E`$ is the energy difference between the constrained and unconstrained systems, by definition a positive quantity. If we define the fractal dimension of the volume of the droplets to be $`d_f`$, i.e.

$$nR^{d_f},$$
(3)
where $`R`$ is the mean “radius” of the droplet (to be defined below) then the slope of the data is $`\theta /d_f`$, since by definition $`\mathrm{\Delta }ER^\theta `$. The expectation of the droplet theory is that $`d_f`$ is equal to the space dimension, $`d_f=d=2`$. The slope of $`0.115`$ for $`L=64`$ would lead to $`\theta =0.23`$, not quite as negative as the established value of $`0.28`$ from domain wall calculation, but perhaps the difference is not significant given the rather small range of droplet sizes than we can study. Actually, fits to data presented below indicate that $`d_f`$ is somewhat less than two, which would accentuate the difference between our value of $`\theta `$ and $`0.28`$. However, the apparent difference between $`d_f`$ and 2 may itself be due to corrections to scaling.

We note also in Fig. 2 that the data for $`L=16`$ shows an interesting feature. For small cluster sizes, the relation between $`\mathrm{\Delta }E`$ and $`n`$ follows that of the $`L=64`$ and the two sets of data start to depart from each other at large cluster sizes, so that the apparent stiffness exponent becomes more negative at large sizes. This crossover is however exactly opposite to the one expected from Eq. (1). A visual inspection reveals that even the clusters with $`n=64`$ have a radius smaller than the system size $`L=16`$, so that the effect is not trivially due to the fact that boundaries of the clusters interact with each other. Rather, this crossover is most probably due a finite size effect on the initial ground state energy, a smaller size making it too large, resulting then in too small an energy difference $`\mathrm{\Delta }E`$.

An important conclusion drawn from Fig. 2 is that no crossover from a large to a smaller value of the stiffness exponent is visible, although our data cover a reasonably large size window. We are led to the conclusion that the definition of the droplets used here is relatively free of the corrections to scaling expected from Eq. (1), and moreover that it leads, for $`d=2`$, to a determination of the stiffness exponent in fair agreement with domain wall calculations.

#### III.1.2 Distribution of Droplet Energies

Fig. 3 shows the distribution of the energy, $`P(\mathrm{\Delta }E)`$, of fixed size droplets for different sizes with $`L=64`$. In Fig. 4, we show that these distributions can be satisfactorily scaled using the form

$$P(\mathrm{\Delta }E)=\frac{1}{\mathrm{\Delta }E}𝒫\left(\frac{\mathrm{\Delta }E}{\mathrm{\Delta }E}\right),$$
(4)
using the mean energy $`\mathrm{\Delta }E`$ which is plotted in Fig. 2. The data indicates that the scaling function $`𝒫(x)`$ varies linearly for small $`x`$. This result was also obtained in three-dimensionsLamarcq et al. (2002), but differs from the usual assumption in droplet theory that the distribution has a finite weight at the origin.

#### III.1.3 Fractal Dimension of the Droplets

To determine the fractal dimension of the volume of the droplets through Eq. (3), we need to give a definition of their “radius”. We have done this in different ways and checked that they give consistent results for $`d_f`$. We measured $`R_x`$ and $`R_y`$, the maximum extent of the droplets in the $`x`$ and $`y`$ direction as well as max($`R_x,R_y`$). In addition we determined the radius of gyration, $`R_g`$, of the droplet defined by:

$$R_g=\sqrt{\frac{1}{n}\underset{i\mathrm{cluster}}{}|𝐫_i𝐫_{\mathrm{cm}}|^2},$$
(5)
where $`𝐫_{\mathrm{cm}}`$ is the position of the center of mass of the cluster.

Figure 5 shows the result of this analysis for $`L=64`$, together with power law fits to the data. The agreement between $`R_g`$ and $`\mathrm{max}(R_x,R_y)`$ is excellent. If we exclude the $`L=4`$ point, the fits give $`d_f`$ closer to 2. For example the fit to $`R_g`$ gives $`d_f1.88`$. Combining these values of $`d_f`$ with $`\theta /d_f`$ obtained from Fig. 2 we get:

$$\theta \{\begin{array}{c}0.23(\text{assuming }d_f=2),\hfill \\ 0.21\text{(fit to }R_g),\hfill \\ 0.22\text{(fit to }R_g,\text{excluding }n=4).\hfill \end{array}$$
(6)
We emphasize again that no trace of a more negative exponent of the order of $`0.47`$ can be detected in our data, thus in disagreement with the predictions by MooreMoore (2002) that this method should reveal the corrections to scaling included in Eq. (1).

### III.2 Clusters with a Range of Size

In this section we allow the size of the droplets to vary over a factor of two, more precisely between $`n/2`$ and $`n1`$. This is more in the spirit of the droplet picture of Fisher and HuseFisher and Huse (1986, 1987, 1988), where droplets are defined, in $`3d`$, as objects with boundaries lying between a cube of linear size $`\mathrm{}`$ and another one of size $`2\mathrm{}`$. We now perform the same analysis as above for these new clusters.

#### III.2.1 Droplet Energies

The mean droplet energy as a function of size $`n`$ is shown in Fig. 6. For the size specified as $`n`$, the droplets were allowed to vary in size between $`n1`$ and $`n/2`$. For very small sizes, where the data for $`L=16`$ and 64 agree, the slope is about $`0.22`$, considerably more negative than the value of $`0.115`$ found for droplets of fixed size in Sec. III.1. For larger sizes, where the data for $`L=16`$ differs from that for $`L=64`$, the slope of the presumably more reliable $`L=64`$ data becomes less negative and is compatible with $`0.115`$, though there is not enough range of size to allow for a precise fitted value. If we insert $`d_f=2`$, then, according to Eq. (3), this would give a crossover from $`\theta 0.44`$ to $`\theta 0.23`$, not very different from the crossover proposed in Refs. Moore, 2002; Hartmann and Moore, 2002.

However, it should be pointed out that no such crossover was found in Sec. III.1 for fixed size droplets, and, furthermore, the values of the droplet radius $`R`$, discussed in Sec. III.2.3 are extremely small, of order one or two lattice spacings, in the region where the data of Fig. 6 suggests $`\theta 0.44`$, so that the meaning of a power law correction of the type $`L^\omega `$ for such sizes is not clear. Note finally that the apparent fractal dimension for these very small clusters is quite far from 2 (see below), so that, if we accept this value for $`d_f`$ the initial apparent stiffness exponent is around $`0.35`$, rather far from the value $`0.47`$ expected by MooreMoore (2002); Hartmann and Moore (2002).

Finally, we have to admit that we do not understand why the crossover in Fig. 6 is absent when the size is exactly fixed, see Fig. 2.

#### III.2.2 Distribution of Droplet Energies

The distribution of droplet energies is shown in Fig. 7 and the scaled data is plotted in Fig. 8, again making use of Eq. (4). As for droplets of fixed size, the distribution is essentially linear for small energy, though in this case there is a very small, but definitely non-zero, intercept.

#### III.2.3 Fractal Dimension of the Droplets

Figure 9 shows the various definitions of the radius of the droplets, discussed in Sec. III.1.3, for different values of $`n`$. The data does not fit well a straight line if the smallest size is included ($`n=4`$, corresponding to droplets with 2 and 3 spins) so this has been omitted in the fits. The values of $`d_f`$ for $`R_g`$ and max($`R_x,R_y`$) agree well, as for droplets of fixed size, though the value of $`d_f(1.7)`$ is a little lower than in that case. However, the trend is for the slope to increase with increasing $`n`$, this is particularly noticeable for the $`R_g`$ data, and so it is possible that $`d_f`$ is actually equal to 2, as we concluded also in the case of the cluster of fixed size.

## IV Conclusions

The main motivation for our work is to test Moore’s proposalMoore (2002) that contradictory results of the literature, both in 2 and 3 dimensions, for the determination of the stiffness exponent $`\theta `$ in spin glasses were mainly due to “corrections to scaling” and captured by Eq. (1). We have performed simulations similar to those of Lamarcq et al.Lamarcq et al. (2002, 2003) but in the $`2d`$ case, where this theoretical conjecture could be checked. From our results, we can draw the following conclusions.

1. Following exactly the procedure of Ref. Lamarcq et al., 2002, we do not find any crossover in the behavior of the energy of the cluster as a function of their size, up to $`n=64`$, contrary to Moore’s prediction;
2. The value of the extracted stiffness exponent, $`\theta 0.23`$, is in fair agreement with domain wall calculations, $`\theta =0.28`$. This result is clearly different from the $`3d`$ case, where the droplet and domain wall exponents are found to be very different.
3. For scale invariant clusters, a crossover not inconsistent with (1) is obtained, but only for very small cluster sizes.

As a conclusion, the problems we wanted to tackle in this paper unfortunately remain unsolved. We still do not understand the $`3d`$ results by Lamarcq et al.Lamarcq et al. (2002) and we do not understand the origin of the two different stiffness exponents reported in $`2d`$. Apart from very small sizes for one of the two methods of generating clusters, we did not find evidence for the crossover implied by Eq. (1). The nature of corrections to scaling in spin glasses therefore remains poorly understood.

Suppose we take the point of view that the differing claims obtained in Refs. Hartmann and Young, 2002; Hartmann and Moore, 2002; Picco et al., 2002a, b and the present work for the existence, or otherwise, of the proposed crossover in two-dimensions, are evidence for a wide range of values of the amplitude of the correction $`B`$ in Eq. (1). It is then surprising that, to our knowledge, all estimates of $`\theta `$ from Monte Carlo simulations at $`TT_c`$ give $`\theta `$ (for droplets) close to (and consistent with) zero, and different from $`\theta `$ for domain walls which is positive. This is true for Ising spin glasses in three and four dimensionsKatzgraber et al. (2001), vector spin glass modelsKatzgraber and Young (2001, 2002), and a one-dimensional model with long range interactionsKatzgraber and Young (2003). In addition, various $`T=0`$ calculations Krzakala and Martin (2000); Palassini and Young (2000); Marinari and Parisi (2000) also find $`\theta \text{(droplets)}0`$. If the amplitude of the correction term is highly non-universal, as perhaps implied by the $`2d`$ results, it is surprising to us that the parameters in Eq. (1) systematically conspire to give a droplet energy for, say, Monte Carlo simulations which is independent of size over the range studied. Note too, that for the $`1d`$ simulationsKatzgraber and Young (2003), the range of sizes is very large, up to $`L=512`$.

Of course, the usual criticism that the truth lies beyond the reach of numerical simulations may be true. In this case, it is worth discussing the relevance of physics which can only be seen on such astronomically long time scales. From a mathematical point of view, understanding this physics is of fundamental importance. However, if the asymptotic behavior would only be seen on time scales which are inaccessible to experiments as well as to simulations, then the physical relevance is rather limited. In fact, experimental systems below the spin glass transition temperature $`T_c`$ are not truly in equilibrium but are only equilibrated up to some coherence length scale, $`\mathrm{}(t)`$, which increases slowly, presumably logarithmically, with the time of the experiment $`t`$ and also depends on $`T`$. One of us and BouchaudBerthier and Bouchaud (2002) have estimated that $`\mathrm{}(T)10`$ (in units of the spacing between spins) at reasonable experimental time scales at $`T0.5T_c`$. Refs. Jonsson et al., 2002; Nordblad, 2003 estimate a somewhat larger value (rather more than 20) at this temperature and quite a lot largerJonsson et al. (2002) closer to $`T_c`$. Hence a reasonable estimate of $`\mathrm{}(t)`$ well below $`T_c`$ in experiments is around 10-20 times the spacing between the spins, not very different from the sizes that can be simulated numerically, and much less than the sizes which appear to be necessary to see droplet behavior.

It is therefore possible that even if the droplet theory is asymptotically correct, the region in which its predictions can be observed quantitatively is not accessible experimentally. Close to $`T_c`$ the length scales which can be equilibrated in experiments may be much largerBerthier and Bouchaud (2002); Jonsson et al. (2002), but in this region critical fluctuations give significant corrections to droplet behavior so even larger length scales are neededMoore et al. (1998); Bokil et al. (2000) to see droplet behavior than at lower temperatures.

###### Acknowledgements.

APY acknowledges support from the National Science Foundation under grant DMR 0086287 and the EPSRC under grant GR/R37869/01. He also thanks David Sherrington for hospitality during his stay at Oxford when part of this work was performed. LB’s work is supported by a European Marie Curie Fellowship No HPMF-CT-2002-01927, CNRS (France) and Worcester College, Oxford. Numerical results were obtained on Oswell at the Oxford Supercomputing Center, Oxford University.


text:
# Orbital ordering in charge transfer insulators

## Abstract

We discuss a new mechanism of orbital ordering, which in charge transfer insulators is more important than the usual exchange interactions and which can make the very type of the ground state of a charge transfer insulator, i.e. its orbital and magnetic ordering, different from that of a Mott-Hubbard insulator. This purely electronic mechanism allows us to explain why orbitals in Jahn-Teller materials typically order at higher temperatures than spins, and to understand the type of orbital ordering in a number of materials, e.g. K2CuF4, without invoking the electron-lattice interaction.

Since the seminal Zaanen, Sawatzky, and Allen paper ZaanenSawatzkyAllen , it is widely accepted that strongly correlated insulators can be divided into two classes according to the nature of their gap. The Mott-Hubbard insulators (MHI) are characterized by the lowest-energy excitations of the type $`d^n+d^nd^{n1}+d^{n+1}`$, obtained by transferring a $`d`$-electron from one transition metal (TM) ion to another, whereas in charge transfer insulators (CTI) the corresponding excitations are holes on ligand sites (oxygen ions in the case of TM oxides): $`d^n+p^6d^{n+1}+p^5`$. MHI are usually described by an effective Hubbard model, in which only the $`d`$-electron states are considered explicitly. On the other hand, CTI are described by models that include also the relevant oxygen $`p`$-states, such as the $`dp`$ (or three-band Hubbard) model of high-$`T_c`$ superconductors Emery . The value of the gap is determined by the smallest of the two parameters: the energy of electron transfer from ligand to TM sites $`\mathrm{\Delta }`$ and the Hubbard repulsion on TM sites $`U_d`$. For CTI (usually, the oxides of heavier $`3d`$ metals) $`\mathrm{\Delta }<U_d`$, while for MHI $`U_d<\mathrm{\Delta }`$.

The distinction between these two types of insulators becomes apparent upon doping and may play an important role in high-$`T_c`$ superconductors Emery ; Varma , though according to other point of view, there is no qualitative difference between doped CTI and MHI ZhangRice . In the absence of doping, both CTI and MHI with one electron or hole per TM site are $`S=\frac{1}{2}`$ Heisenberg antiferromagnets.

The situation turns out to be quite different in Jahn-Teller systems, in which $`d`$-electrons, in addition to spins, have also orbital degrees of freedom. In this Letter we show that in insulators with orbital degeneracy not only the excitation spectrum, but also an orbital and magnetic ordering in the ground state may strongly depend on the ratio $`\mathrm{\Delta }/U_d`$. The difference between the ground states of CTI and MHI results from a peculiar exchange interaction between orbitals on nearest-neighbor $`d`$-sites, which may be called ‘the orbital Casimir force’, as it originates from the vacuum energy of charge fluctuations due to the hopping of electrons between the ligand and TM ions (whereas the usual exchange involves the electron transfer from one TM site to another). We recently noted the importance of such interaction in geometrically frustrated systems with the $`90^{}`$ metal-oxygen-metal bonds, e.g. NaNiO2, where it determines the type of orbital ordering MostovoyKhomskii . Here we show that this orbital exchange of electronic origin is generally present in all TM oxides and, in particular, in materials with the $`180^{}`$ metal-oxygen-metal bonds, like perovskites ABO3 and their layered analogs A2BO4, and it is not included in the widely used exchange model KugelKhomskiiJETP ; KugelKhomskiiUFN  (sometimes called Kugel-Khomskii (KK) model). That model is derived from an effective Hubbard model describing only the degenerate $`d`$-electron orbitals and which, strictly speaking, only applies to MHI with $`\mathrm{\Delta }U_d`$. The new orbital interaction obtained below helps to resolve a difficulty of the KK model, in which an orbital ordering is strongly coupled to a magnetic one, implying close temperatures of orbital and magnetic orderings KugelKhomskiiUFN , whereas in reality orbitals usually order at higher temperatures than spins.

Exchange interactions: We first calculate an exchange Hamiltonian for two TM ions with one hole on the doubly degenerate $`e_g`$-level, e.g. Cu2+ ($`t_{2g}^6e_g^3`$). The two orthogonal $`e_g`$-orbitals $`|3z^2r^2`$ and $`|x^2y^2`$ may be described by, respectively, the up and down states of the isospin-$`\frac{1}{2}`$ operator $`T^z`$ KugelKhomskiiJETP ; KugelKhomskiiUFN . It is convenient to introduce an overcomplete set of operators $`d_{i\alpha \sigma }`$ on each TM site $`i`$ annihiliating a hole in the state $`|3\alpha ^2r^2`$ ($`\alpha =x,y,z`$) with the spin projection $`\sigma `$. These states are hybridized with the $`p_\alpha `$ orbitals on two neighboring oxygen sites $`i\pm \frac{\alpha }{2}`$ (the corresponding oxygen hole operators are denoted by $`p_{i+\frac{\alpha }{2},\sigma }`$). The hole states on TM sites are described by the standard atomic Hamiltonian, which includes the on-site repulsion $`U_d`$ and the Hund’s rule exchange coupling $`J_d`$ for holes on two orthogonal orbitals. The Hamiltonian of oxygen holes hybridized with the $`d`$-states is

$`H_p`$ $`=`$ $`H_{kin}+U_p{\displaystyle \underset{i,\alpha }{}}p_{i+\frac{\alpha }{2}}^{}p_{i+\frac{\alpha }{2}}p_{i+\frac{\alpha }{2}}^{}p_{i+\frac{\alpha }{2}}`$ (1)
$``$ $`t_{pd}{\displaystyle \underset{i,\alpha ,\sigma }{}}[d_{i\alpha \sigma }^{}(p_{i+\frac{\alpha }{2},\sigma }+p_{i\frac{\alpha }{2},\sigma })+h.c.]`$
where the first term is the kinetic energy of holes in the oxygen $`p`$-band, the second term is the Hubbard repulsion on oxygen sites, and the last term describes the hopping between oxygen and TM sites.

Integrating the oxygen degrees of freedom out, one obtains an effective interaction between $`d`$-electrons in the form of an expansion in powers of $`(t_{pd}/\mathrm{\Delta })^2`$. The lowest order term provides a hopping between TM sites, which together with the Coulomb and Hund’s rule interactions between $`d`$-electrons gives an effective $`dd`$-model. For dispersionless oxygen holes, $`H_{kin}=\mathrm{\Delta }_{i,\alpha ,\sigma }p_{i+\frac{\alpha }{2}\sigma }^{}p_{i+\frac{\alpha }{2}\sigma }`$, the effective $`dd`$-hopping occurs between nearest-neighbor TM sites, and its amplitude is $`t_{dd}=t_{pd}^2/\mathrm{\Delta }`$.

The next term in this expansion ($`t_{pd}^{}{}_{}{}^{4}`$) gives an orbital and spin exchange mediated by two oxygen holes. The corresponding diagrams are shown in Fig.1, where in the diagrams (a) and (b) two holes do not interact with each other in the virtual state, while the black squares on diagrams (c) and (d) correspond to one or more Hubbard interactions on oxygen sites. Each diagram gives rise to a particular form of spin and orbital exchange. Thus, if the holes on the TM sites $`i`$ and $`j`$ exchange places \[Fig.1(a)\], the spin exchange is described by operator $`P_{ex}=+1`$ if the total spin of the two holes $`S=1`$, and $`1`$ for $`S=0`$. This operator can also be written in the form $`P_{ex}=12P_0`$, where $`P_0=1/4\left(𝐒_i𝐒_j\right)`$ is the projector on the $`S=0`$ state. If, on the other hand, the holes return back to their starting points, the corresponding ‘spin factor’ is $`1`$. The interacting holes in the diagrams (c) and (d) must have opposite spin projections, so the spin exchange in that case is described by $`S_i^+S_j^{}+S_i^{}S_j^+`$ \[see Fig.1(c)\], while if the spin projections remain unchanged \[see Fig.1(d)\], the spin operator is $`2\left(1/4S_i^zS_j^z\right)`$. Since the diagrams (c) and (d) only differ by the spin factors, their sum is proportional to $`2P_0`$, corresponding to the fact that the on-site Hubbard interaction between two holes is only nonzero when their spins add to the total spin $`0`$.

While the spin-dependence of the exchange interactions, discussed above, is very general, the orbital structure depends on the dispersion of holes and on relative positions of the sites $`i`$ and $`j`$. In particular, for dispersionless holes the exchange in order $`t_{pd}^4`$ only occurs between two nearest neighbor TM sites, e.g. $`i`$ and $`i+\alpha `$, and the diagrams (a), (c), and (d) are proportional to $`_{i\alpha }=\left(\frac{1}{2}+I_i^\alpha \right)\left(\frac{1}{2}+I_{i+\alpha }^\alpha \right),`$ where

$$I^x=\frac{1}{2}T^z\frac{\sqrt{3}}{2}T^x,I^y=\frac{1}{2}T^z+\frac{\sqrt{3}}{2}T^x,I^z=T^z,$$
(2)
and the operator $`\left(\frac{1}{2}+I_i^\alpha \right)=d_{i\alpha }^{}d_{i\alpha }`$ is the projector on the state $`|3\alpha ^2r^2`$ on the site $`i`$. The contribution of the diagram in Fig.1(b) is proportional to $`4_{\beta ,\gamma }\left(\frac{1}{2}+I_i^\beta \right)\left(\frac{1}{2}+I_j^\gamma \right)=9`$, since $`I^x+I^y+I^z=0`$.

Combining terms containing $`P_0`$ with the exchange interactions appearing in the dd-model in the order $`t_{pd}^4`$, we obtain

$`H_{ST}={\displaystyle \frac{4t_{pd}^4}{\mathrm{\Delta }^2}}{\displaystyle \underset{i,\alpha }{}}[({\displaystyle \frac{1}{U_d}}+{\displaystyle \frac{2}{2\mathrm{\Delta }+U_p}})_{i\alpha }P_0`$
$`+{\displaystyle \frac{J_d}{2U_d^2}}({\displaystyle \frac{1}{4}}I_i^\alpha I_{i+\alpha }^\alpha )P_1],`$ (3)
where the last term describes the Hund’s rule coupling on TM sites and $`P_1=1P_0=3/4+\left(𝐒_i𝐒_j\right)`$ is the projector on the $`S=1`$ state. In addition to the coupled spin and orbital interactions, there exists also a spin-independent orbital exchange

$`H_T^{(h)}`$ $`=`$ $`{\displaystyle \frac{2t_{pd}^4}{\mathrm{\Delta }^2}}{\displaystyle \underset{i,\alpha }{}}[{\displaystyle \frac{1}{\mathrm{\Delta }}}({\displaystyle \frac{1}{2}}+I_i^\alpha )({\displaystyle \frac{1}{2}}+I_{i+\alpha }^\alpha )`$ (4)
$``$ $`{\displaystyle \frac{1}{U_d}}({\displaystyle \frac{1}{4}}I_i^\alpha I_{i+\alpha }^\alpha )]`$
The first term is the spin-independent contribution of the processes shown in Figs.1(a) and (b), while the second one results from the $`dd`$-model. It is important to note that in the order $`t_{pd}^4`$ it is impossible to exchange two holes in such a way that they never occupy the same oxygen site (the corresponding terms in the diagrams (a) and (b) cancel each other). The orbital exchange interaction containing only powers of $`\mathrm{\Delta }`$ is actually due to the virtual hopping of holes to all neighboring oxygen sites and back \[Fig.1(b)\], except for the oxygen bridging the two neighboring TM sites. In the limit $`\mathrm{\Delta }U_d`$ Eqs. (3) and (4) reduce to the KK model.

The term in $`H_T`$ containing only $`\mathrm{\Delta }`$ in denominator is similar to the pure orbital exchange for the $`90^{}`$ metal-oxygen-metal bonds MostovoyKhomskii , and its origin can be better understood by noting that it is the only term that does not vanish in the limit of infinite on-site Hubbard repulsion $`U_p`$ and $`U_d`$. It represents the energy of the delocalization of holes over the octahedra of the ligand ions, surrounding the TM sites, or more correctly, the increase of that energy due to the fact that the neighboring octahedra share one ligand ion and for infinite $`U_p`$ the intermediate two-hole states on the shared ion are forbidden. This energy increase depends on orbital states of the neighboring holes, which gives rise to the orbital exchange with the positive sign. Thus the origin of this orbital interaction is similar to that of the Casimir effect Casimir , resulting from the dependence of the energy of zero-point fluctuations of electromagnetic field on positions of two metallic plates. It is also important that this orbital exchange has positive sign, while all other terms in Eqs. (3) and (4) are negative (since they describe the second-order corrections to the ground state energy in the effective $`dd`$-hoping). The result of this new orbital interaction can be formulated as the following “rule of thumb” (applied also to materials with one electron per TM site): to avoid the energetically unfavorable occupation of the bridging oxygen by two holes, the d-holes should occupy orbitals directed away from the common oxygen either on both neighboring TM sites \[see Fig.2(b)\] or, at least, on one site \[see Fig.2(a)\].

The importance of the orbital exchange Eq.(4) for an orbital ordering may be understood already by considering the minimal-energy states of a pair of holes on the neighboring TM sites $`i`$ and $`i+z`$. In MHI ($`\mathrm{\Delta }>U_d`$) the exchange is dominated by $`H_{ST}`$ (3) and the last term in Eq.(4). The AFM exchange between spins of the holes, described by this term, is the strongest for $`I^z=\frac{1}{2}`$ on both sites, i.e., the lowest energy is achieved if both holes occupy the $`|3z^2r^2`$ state, as shown in Fig.2(c). The ground state of the pair is a spin singlet. On the other hand, in CTI ($`\mathrm{\Delta }<U_d`$) the strongest interaction is the orbital exchange $`H_T`$ (see Eq.(4)), and the state shown in Fig.2(c) has the highest orbital energy, since both hole orbitals are oriented towards the common ligand ion. The more energetically favorable states in CTI will be those, in which at least one hole is in the $`|x^2y^2`$ state (see Fig.2a and b). Furthermore, both the last term in Eq.(4) and the relatively weak FM spin exchange due to the Hund’s rule coupling select the state shown in Fig.2a with the total spin of the pair $`S=1`$. Thus, the ground state of the pair in the CT regime is different from that in the MH regime.

Consider now two neighboring TM sites with one $`e_g`$ electron, as e.g. Mn3+ ($`t_{2g}^3e_g^1`$), Cr2+ ($`t_{2g}^3e_g^1`$), and the low-spin Ni3+ ($`t_{2g}^6e_g^1`$)). In the KK model the systems with 1 electron/site and 1 hole/site order in exactly the same way, since the Hubbard model, from which the KK exchange Hamiltonian is derived, has the electron-hole symmetry KugelKhomskiiUFN ; consequently, the term $`H_{ST}`$ (3) has the same form both for electrons and holes. On the other hand, in the $`dp`$-model this symmetry is lost: the new orbital exchange interaction for one-electron case has the form \[cf. Eq.(4)\]

$`H_T^{(e)}`$ $`=`$ $`+{\displaystyle \frac{2t_{pd}^4}{\mathrm{\Delta }^2}}{\displaystyle \underset{i,\alpha }{}}[{\displaystyle \frac{1}{\mathrm{\Delta }}}({\displaystyle \frac{3}{2}}I_i^\alpha )({\displaystyle \frac{3}{2}}I_{i+\alpha }^\alpha )`$ (5)
$``$ $`{\displaystyle \frac{1}{U_d}}({\displaystyle \frac{1}{4}}I_i^\alpha I_{i+\alpha }^\alpha )]`$
In this case for both MHI and CTI the energy of a pair of TM electrons on the sites $`i`$ and $`i+z`$ is minimal if both electrons occupy the $`d_{3z^2r^2}`$ orbital, Fig.2(c), and their spins form the singlet state.

Note that the electron-phonon (Jahn-Teller) interaction also gives rise to a purely orbital exchange of the form $`H_{JT}_{i,\alpha }I_i^\alpha I_{i+\alpha }^\alpha `$, which for both the electron and hole cases favors the antiferro-orbital (AFO) spin-triplet state shown in Fig.2(a). Thus the JT interaction favors the same state as the novel orbital exchange for one-hole systems, but a different one for one-electron systems. This, in principle, gives a possibility to experimentally discriminate between these two mechanisms and to study their relative importance. The Cu2+ pairs (one hole) in K2ZnF4 were found to have the triplet ground state Eremin , showing that the orbital exchange Eq.(4) and/or JT interaction are in this case stronger than the KK exchange. We are not aware of similar data for pairs of one-electron ions.

Orbital and magnetic ordering: Equations (35) allow us to consider types of ordering in different situations. Thus, e.g., the ground state of a chain of TM ions in the $`z`$-direction with 1 hole/site is AFO with alternating $`d_{3z^2r^2}`$ and $`d_{x^2y^2}`$ orbitals and FM for CTI, as that is favored by the sum of Eqs. (3) and (4), while for MHI, in which $`H_{ST}`$ dominates, we obtain the ferro-orbital ordering with $`d_{3z^2r^2}`$ orbitals occupied on all sites, and the AFM spin exchange lowd . The latter state also minimizes the exchange energy of one-dimensional CTI with 1 electron/site.

Similarly e.g. for CTI ($`U_d\mathrm{}`$) with a square lattice and 1 hole/site, the orbital exchange interaction Eq.(4) has the form

$$H_T=+\frac{2t^4}{\mathrm{\Delta }^3}\underset{i}{}\left(I_i^xI_{i+x}^x+I_i^yI_{i+y}^yI_i^z\right).$$
(6)
The term bi-quadratic in isospin operators favors an antiparallel orientation of neighboring isospins, while the linear term in Eq.(6) acts as ‘the orbital field’ in the $`z`$ direction, favoring the canted orbital state with the canting angle $`\pm 60^{}`$, which corresponds to the two-sublattice orbital ordering with the $`d_{z^2x^2}`$ and $`d_{z^2y^2}`$ hole states. Since these two neighboring orbitals do not overlap, the Hund’s rule coupling (the last term in Eq.(3)) results in the FM exchange between the spins of the holes. This is precisely the orbital and magnetic structure of the layered compound K2CuCl4. For MHI, however, we would get from Eq.(3) a different ground state: ferro-orbitally ordered with $`I^z=\frac{1}{2}`$ ($`d_{x^2y^2}`$-orbitals) and AFM spin ordering (see, however, Ref. lowd, ). Previously, the KK model was successfully used to predict the orbital structure of that compound KhomskiiKugel , assuming that spins order ferromagnetically; but as we mentioned above, such a state is not the ground state of the KK model.

The case of a cubic lattice of TM ions is somewhat special for several reasons. First, using $`_\alpha I_i^\alpha =0`$, we find that the exchange interactions have the electron-hole symmetry. In particular, both Eq.(4) and Eq.(5) can be written in the form

$$H_T=+\frac{2t^4}{\mathrm{\Delta }^2}\left(\frac{1}{\mathrm{\Delta }}+\frac{1}{U_d}\right)\underset{i}{}\left(I_i^xI_{i+x}^x+I_i^yI_{i+y}^y+I_i^zI_{i+z}^z\right).$$
(7)
Second, the orbital interaction Eq.(7) is frustrated FeinerOlesZaanen ; GKhaliullin : apart from the AFO mean-field ground states with the average isospin vector $`𝐓_i`$ parallel (antiparallel) to an arbitrary vector in the $`(T^x,T^z)`$-plane, there are infinitely many disordered ground states, obtained from the AFO states by performing the transformation $`I^\alpha I^\alpha `$, $`I^\beta I^\gamma `$, where $`\alpha `$, $`\beta `$, and $`\gamma `$ are three different spatial indexes, on an arbitrary set of planes transverse to the $`\alpha `$ axis. This ground state degeneracy is similar to the one of the orbital Hamiltonian of $`e_g`$ electrons on a triangular lattice MostovoyKhomskii . The quantum orbital fluctuations select the AFO states $`I^\alpha =\pm \frac{1}{2}`$, where $`\alpha =x,y`$ or $`z`$, which have the softest (two-dimensional) orbital excitations GKhaliullin . For instance, for $`\alpha =z`$, such a state corresponds to the occupation of $`d_{3z^2r^2}`$ and $`d_{x^2y^2}`$ orbitals on A and B sublattices. The corresponding magnetic ordering would be AFM in the $`xy`$ planes and FM in the $`z`$ direction.

If we now would increase the ratio $`\mathrm{\Delta }/U_d`$, moving from CTI to MHI, the coupling between orbitals and spins would ultimately stabilize the canted AFO state with an AFM coupling along the chains in the $`z`$ direction and a FM exchange between spins in the $`xy`$ layers. The detailed treatment of this case, relevant e.g. for the perovskites KCuF3 and LaMnO3, will be given elsewhere.

In conclusion, we showed that the description and properties of orbital ordering in charge transfer insulators are much different from those in Mott-Hubbard insulators. In particular, even in perovskite CTI with the $`180^{}`$-exchange the orbital interactions are stronger than the spin exchange, so that an orbital ordering, in general, occurs at higher temperatures. In this respect the exchange mechanism of an orbital ordering in CTI resembles the ordering due to the JT (electron-phonon) interaction, although is not identical to it. Another difference between the CTI and MHI is that the electron-hole symmetry existing in the MH case is lost in the CT regime, so that the type of an orbital and spin ordering may be different for one-electron and one-hole materials. The effective exchange Hamiltonian for CTI allows us to describe the ordering in a number of such materials, without invoking any extra factors. Most importantly, in orbitally degenerate insulators not only the excitation spectrum, but also the very type of the ground state may be different in Mott-Hubbard and charge transfer regimes, in contrast to the non-degenerate case, in which it is the same.

This work was supported by the MSCplus program and the Stichting voor Fundamenteel Onderzoek der Materie (FOM).