[312], [313]).

Above, we have considered the $U\to\infty$ limit. This means that the
spin–spin, spin–orbital, and orbital–orbital correlations between neighboring
sites, which take place in two-band Hubbard Hamiltonians (204) and (211), are
not taken into account. Such correlations can affect the actual types of
magnetic and orbital orderings. However, the corresponding terms are of the
order of $t^{2}/U\ll t,\,\Delta,\,J_{H}$. Hence, they are rather small and
should not significantly change the phase diagrams. Using the mean-field
approach, we also neglect the fluctuation effects, but it is known that such
approach captures the main qualitative features of the phase diagrams.

The long-range Coulomb interaction related to the charge inhomogeneity in the
phase-separated state can reduce the doping range corresponding to the phase
separation and affect the size of inhomogeneities. The minimization of both
the long-range Coulomb repulsion and the surface energy provides an
opportunity ro reveal the geometry of the phase-separated state (e.g., the
inhomogeneities can take the shapes of droplets, alternating layers, or rods)
and to determine the characteristic sizes of inhomogeneities and the
percolation threshold for the doping, at which the insulator–metal transition
occurs [94, 218, 191].

###  6.6 Conclusions

In this section, we have formulated a simple model capturing the main physics
of strongly correlated electron systems with spin-state transitions. Based on
this model, we predict the possibility of electronic phase separation for
perovskite cobaltites under doping. We determined the corresponding phase
diagrams exhibiting broad ranges of the phase separation. The specific shape
of phase diagrams exhibits significant changes with the variation of the ratio
of the electron hopping integral $t$ and and the Hund’s rule coupling constant
$J_{H}$.

We did not analyze in detail the possible structure of the phase-separated
state. However, as we have discussed in the previous sections in the case
manganites, the calculations [94] and numerical simulations [257] taking into
account the surface and long-range Coulomb contributions to the total energy
lead to the characteristic size of nanoscale inhomogeneities of the order of
several lattice constants. One can expect that the inhomogeneities caused by
the phase separation in cobaltites with spin-state transitions would be of a
similar scale.

Note that the existence of the electronic phase separation in the systems with
the spin-state transitions is rather robust phenomenon appearing in different
models. For example, the numerical studies reported in [314] revealed a
possibility of the phase separation in the two-band model with the bands
differing in the spin state. In such a model, the phase separation appears
upon doping if the widths of these band are comparable to each other.

It is also worth noting rather active line research related to the
photoinduced spin-state transitions and the formation of the corresponding
inhomogeneous states, which can be revealed,for example, by the pump-probe
techniques. In particular, a bound state of optically exited hole and the HS
state can be created in cobaltites within the LS phase. This bound state
brings about a characteristic peak structure in the optical pump-probe spectra
which is distinct from the spectra in thermally excited states [315, 316].
This brings us a clue for the new optical spin control technique.

The nanoscale phase separation in materials with spin state transitions can
involve some additional degrees of freedom. The detailed studies of
La2-xSrxCoO4 layered cobaltites near $x=1/3$ in the experiments including the
resonant micro X-ray diffraction, neutron scattering, and muon spin rotation
spectroscopy revealed nanosize undoped and hole doped islands involving the
charge order correlations in addition to the magnetic ones [317, 318].

##  7 Effects of orbital degrees of freedom, different shapes of polarons

###  7.1 Introduction

In this section, we study the effects related to one more degree of freedom,
namely with orbitals and their ordering. The orbital ordering (OO) plays an
important role in the physics of magnetic oxides, especially of those, which
contain the ions with the orbitally degenerate ground state, so called
Jahn–Teller (JT) ions. The compounds with JT ions usually exhibit different
kinds of superstructures (OO) [38, 319]. In fact, the OO is a specific feature
of magnetic insulators. The additional charge carriers introduced by doping
are able to transfer an OO insulator to the metallic state. Hence, here we are
dealing again with the competition between the charge carrier localization and
their itinerancy, which can give rise to the phase separation, especially at
low doping levels [112, 111, 28]. In the analysis of phase separation in
magnetic materials such as manganites, the OO is usually not addressed (except
maybe the treatment of orbital and magnetic polarons [294, 28, 320, 321]). In
this section, we study the role of OO for the electronic phase separation
employing minimal models including itinerant charge carriers interacting with
the OO superstructure. We demonstrate that at a low doping level, the
orbitally ordered systems can exhibit the phase separation, which not directly
related to the possible magnetic structure.

Such phase separation is treated using two approaches. In subsection 7.2, we
base our analysis on a symmetrical model analogous to the considered above
double exchange model where the orbital variables play a role of local spins.
We assume that owing to the orbital degeneracy, localized electrons induce
lattice distortions of JT type, thus giving rise to OO. Itinerant charge
carriers appearing due to the doping move on OO background. Within the other
approach (subsection 7.3), we employ the $e_{g}$ symmetry of the wave
functions for the doped charge carriers. Using both approaches, we reveal the
tendency to the phase separation related to the existence of the orbital
degrees of freedom. In other words, itinerant charge carriers introduced by
doping stimulate the formation of inhomogeneous states with nanosize domains
having the orbital order differing from that taking place in the absence of
doping. In subsection 7.4, we analyze the possible shapes and sizes of such
domains and demonstrate that at different relations between the electron
hopping integral $t$ and the parameter $J$ characterizing the interaction
between orbitals, such shapes can be quite different. For example, in the
layered compounds, the two-dimensional case, below some critical value of of
$t/J$, nearly circular domains can undergo the transformation to the needle-
like ones. This is a consequence of anisotropic charge distribution
characteristic of specific orbitals, which brings forward additional
possibilities to the physics of phase separation.

In subsection 7.5, we are dealing with a picture more relevant to the actual
JT compounds, in which the same charge carriers are involved both in OO and,
after doping, in conduction. We assume that at zero doping, the orbitals are
ordered due to the electron–lattice coupling. The doped charge carriers are
moving on the OO background. In subsection 7.6, we consider another model,
where the same $e_{g}$ electrons are responsible for the band structure and
for the orbital ordering.

###  7.2 Symmetrical model

Here, we analyze the situation typical of compounds JT ions having a double-
degenerate ground state. Different kinds of lattice distortions (e.g.
compression or stretching) lift such degeneracy making favorable two possible
ground states for each ion, say, either $a$ or $b$. The $a$ and $b$ states at
lattice site $\mathbf{n}$ are, describe, in fact, the occupation of electron
orbitals at this site. Usually we deal with the linear combination of $a$ and
$b$ orbitals, which can be expressed in terms of angle $\theta$

$|\theta\rangle=\cos\frac{\theta}{2}|a\rangle+\sin\frac{\theta}{2}|b\rangle\,.$
(214)

The orbitals and the related local distortions at different sites can form
regular arrays (orbital ordering). In the case of two basis states, the
orbitals can be put into correspondence to some spin-1/2 variables, and the
interaction between orbitals in its simplest version can be written in the
form of Heisenberg Hamiltonian

$H_{\text{OO}}=J\sum_{\langle{\bf n}{\bf m}\rangle}\bm{\tau}_{\bf
n}\bm{\tau}_{\bf m}\,,$ (215)

where $\bm{\tau_{n}}=\\{\tau^{x}_{\mathbf{n}}$, $\tau^{z}_{\mathbf{n}}\\}$ are
the Pauli matrices, and $a$ and $b$ states of the ion ${\bf n}$ correspond to
eigenvectors of operators $\tau^{z}_{\bf n}$, with eigenvalues $1$ and $-1$,
respectively. Hamiltonian (215), provides a possibility of two simplest types
of orbital ordering: ferro-OO (similar orbitals at each site) and antiferro-OO
(alternating orbitals at neighboring sites). At zero doping (no itinerant
charge carriers), the system can have either the antiferro-OO ground state (at
$J>0$) or the ferro-OO one (at $J<0$). In actual materials containing
Jahn–Teller ions, orbital Hamiltonians are not so simple. Nevertheless, model
(215) appears to be sufficient for capturing main physical effects related to
the orbital ordering.

Let us now pass to the case of nonzero doping and hence of nonzero density of
itinerant charge carriers in the system, $n\neq 0$. These itinerant charge
carriers are moving on the background of orbitally ordered state formed by
localized electrons. (Below, we demonstrate that the main qualitative results
remain the same, if itinerant charge carries such as $e_{g}$ electrons, are
responsible both for the orbital ordering and for the electrical
conductivity). The electron hopping integrals are determined by the specific
electron orbitals located at the neighboring lattice sites. The Hamiltonian
describing only the electron states have the form

$H_{\text{el}}=-\sum_{\langle{\bf n}{\bf
m}\rangle,\alpha,\beta,\sigma}t^{\alpha\beta}\left(P_{{\bf
n}\alpha\sigma}^{{\dagger}}a^{{\dagger}}_{{\bf n}\alpha\sigma}a_{{\bf
m}\beta\sigma}P_{{\bf m}\beta\sigma}+h.c.\right),$ (216)

where, $a^{{\dagger}}_{{\bf n}\alpha\sigma}$, $a_{{\bf n}\alpha\sigma}$ are
creation and annihilation operators cottespondingg to charge carriers at site
${\bf n}$ with spin projection $\sigma$ at orbital $\alpha$. Projection
operators $P$ in (216) take into account that the double occupancy of the
lattice sites by electrons is highly unfavorable in energy due to the strong
on-site Coulomb repulsion. Below, we are dealing mostly with the case $n<1$.
In this section, our analysis of the electron contribution to the total energy
is based on the tight-binding approximation. To avoid unnecessary
complications, we illustrate the essential physics of orbitally ordered
systems considering the square lattice in the two-dimensional (2D) case and
simple cubic lattice in the three-dimensional (3D) case. In the tight-binding
approximation, the spectrum of charge carriers has the form

$E^{\alpha\beta}({\bf k})=-t^{\alpha\beta}\frac{1}{D}\sum_{i=1}^{D}\cos
k_{i}=t^{\alpha\beta}\xi({\bf k})\,,$ (217)

where $D$ (2 or 3) is the number of spatial dimensions. Components $k_{i}$ of
the wave vector ${\bf k}$ are expressed in the units of inverse lattice
constant $1/d$.

Within the approach adopted in this section, the charge carriers introduced by
doping are built in to the JT distorted background and occupy orbital states
$|a\rangle$ or $|b\rangle$. Then, their intersite charge transfer is described
by three hopping integrals: $t^{aa}$, $t^{bb}$, and
$t^{ab}=t^{ba}=t^{\prime}$. For simplicity, we can assume that
$t^{aa}=t^{bb}=t$. Such system is characterized by an interplay two
contributions to the total energy: the possible gain in the kinetic energy can
come from the increase in the electron bandwidth, whereas the potential energy
can be reduced by the optimization of the orbital ordering.

In a more general case, when there is a superposition of $|a\rangle$ or
$|b\rangle$ states at site $\mathbf{n}$ characterized by angle $\theta$, (see,
by Eq. (214)), we can write the expression for the hopping integral between
the sites characterized by orbital states $|\theta_{1}\rangle$ and
$|\theta_{2}\rangle$

$t^{\theta_{1}\theta_{2}}=t\cos{\frac{\theta_{1}-\theta_{2}}{2}}+t^{\prime}\sin{\frac{\theta_{1}+\theta_{2}}{2}}$
(218)

This is in fact a generalization to the OO case of the semiclassical
expression for the effective hopping integral in the double exchange model
[120].

We shall start from the homogeneous phase, in which orbitals
$|\theta_{1}\rangle$ and $|\theta_{2}\rangle$ alternate at different lattice
sites ($J>0$). Within the mean-field approach, the total energy per site is
given by the expression

$E_{tot}(\theta_{1},\theta_{2})=zt^{\theta_{1}\theta_{2}}\varepsilon_{0}(n)+\frac{zJ}{2}\cos(\theta_{1}-\theta_{2})\,,\;\;\varepsilon_{0}(n)<0\,.$
(219)

Here, $z$ is the number of nearest-neighbor sites and $\varepsilon_{0}(n)$ is
the dimensionless kinetic energy, which depends on the type of crystal
lattice. In this subsection, we consider an isotropic case, for which $t$,
$t^{\prime}$, and, hence, $t^{\theta_{1}\theta_{2}}$, are independent on the
direction in space. For such isotropic system, $\varepsilon_{0}(n)$ does not
depend on $\theta_{1,2}$, and we can easily calculate the orbital structure by
minimization of total energy (219) with respect to angles $\theta_{1}$ and
$\theta_{2}$. At high doping levels, namely for $t|\varepsilon_{0}(n)|>2J$,
the ferro-OO state with $\theta_{1}=\theta_{2}=\pi/2$ becomes favorable. In
the opposite limiting case, $t|\varepsilon_{0}(n)|<2J$, the orbital structure
corresponds to $\theta_{2}=\pi-\theta_{1}$, and

$\theta_{1}=\arcsin\left(\frac{t|\varepsilon_{0}(n)|}{2J}\right)\,,\;\;\frac{t|\varepsilon_{0}(n)|}{2J}<1\,.$
(220)

Such alternation of orbitals can be interpreted as a kind of canting for the
directions of pseudospins $\tau$ in the orbital space. For the energy of such
state, we have

$E_{tot}=zt^{\prime}\varepsilon_{0}(n)-\frac{zt^{2}}{4J}\varepsilon_{0}^{2}(n)-\frac{zJ}{2}\,.$
(221)

If it is possible to represent the kinetic energy as
$\varepsilon_{0}(n)=nf(n)$, where $f(n)$ a slowly changing function of $n$,
then $E_{tot}$ can exhibit a negative curvature, at least, at small $n$, which
signals on an instability of a homogeneous orbitally ordered state (negative
compressibility).

We can find a specific form of $\varepsilon_{0}(n)$ and determine the shape of
the curve depicting the behavior of the total energy with doping. For the
electron energy (217), the density of states $\rho_{0}(E)$, can be written as

$\rho_{0}(E)=\\!\int\\!\frac{d\mathbf{k}}{(2\pi)^{D}}\delta(E-\xi(\mathbf{k}))=\\!\int\limits_{0}^{\infty}\\!\frac{ds}{\pi}\cos(Es)J_{0}^{D}\\!\left(\frac{s}{D}\right)\\!,$
(222)

where $J_{0}$ is the Bessel function. The expression for $\varepsilon_{0}(n)$
n has the form

$\varepsilon_{0}(n)=\int_{-1}^{\mu(n)}\\!\\!\\!\\!\\!dE\,E\rho_{0}(E)\,,$
(223)

where the chemical potential $\mu$ can be found based on the usual relation
$n=\int_{-1}^{\mu}dE\rho_{0}(E)$.

At small density of itinerant charge carriers, $n\ll 1$, we can write an
explicit expression for $\varepsilon_{0}(n)$. For example, in a 2D case, we
have $\varepsilon_{0}(n)\approx-n+\pi n^{2}/2$, and the corresponding
expression for the total energy is

$E_{tot}\approx-zt^{\prime}n-z\left(\frac{t^{2}}{4J}-\frac{\pi
t^{\prime}}{2}\right)n^{2}-\frac{zJ}{2}\,.$ (224)

Using Eq. (224), we can find the range of parameters corresponding to the
negative curvature of the energy plot, i.e. $d^{2}E_{tot}/dn^{2}<0$. Such
situation takes place if

$\frac{t}{J}>\frac{2\pi t^{\prime}}{t},$ (225)

This means that the homogeneous state with canted orbitals is unstable with
respect to the phase separation into a inhomogeneous states with ferro- and
antiferro-orbital ordered domains. Here, we have a close analogy to the
situation typical of the double exchange (in the simplest single-band case)
[135], which corresponds to $t^{\prime}=0$. In the opposite limiting case,
that is, for $2\pi t^{\prime}/t>t/J$, a homogeneous state is favorable within
the whole doping range – in contrast to the usual double exchange.

A rough estimate for the phase separation range can be found if we neglect the
quadratic term in Eq. (220), $\varepsilon_{0}(n)\approx-n$. Thus, we find

$0<n\lesssim\frac{2J}{t}.$ (226)

Actually, the homogeneous orbitally canted state appears to be unfavorable in
energy within almost the whole doping range, where the difference
$\theta_{2}-\theta_{1}=\pi-2\theta_{1}$ with $\theta_{1}$ from Eq. (220) does
not vanish. The situation does not qualitative change, if we use more accurate
calculations for $E_{tot}$ based on $\varepsilon_{0}(n)$ determined employing
density of states found from (222). The corresponding $E_{tot}(n)$ plot in the
2D case is shown in Fig. 47.

Figure 47:  (Color online) Total energy of homogeneous orbitally ordered state
state (221) in 2D versus doping $n$. Two types of behavior are shown: with the
range corresponding to the negative curvature (red solid line, $J/t=0.05$),
and without it (blue dashed line, $J/t=0.35$); $t^{\prime}=0.5t$ for both
curves. The range of instability for the homogeneous state (with the negative
curvature) is located on the left-hand side of vertical line. The energy of
the phase-separated state is drawn by the dot-and-dash line. The blue dashed
curve represents the energy of the homogeneous state, which at $J/t=0.35$ is
stable within the whole doping range [39].

In the 3D case, the situation is not so simple. At low doping levels, the
corresponding expression for the kinetic energy has the form
$\varepsilon_{0}(n)\approx-n+an^{5/3}$, where

$a=\frac{3}{5}\left(\frac{\pi^{2}}{\sqrt{6}}\right)^{2/3}\\!\\!\\!,$

and the total energy is

$E_{tot}\approx-
zt^{\prime}n-z\left(\frac{t^{2}}{4J}-\frac{a}{n^{1/3}}\right)n^{2}-\frac{zJ}{2}\,.$
(227)

The second derivative of $E_{tot}$ is positive at $n\to 0$ and becomes
negative at

$n_{c}\approx\left(\frac{5aJt^{\prime}}{9t^{2}}\right)^{3}.$ (228)

The same line of reasoning as that in the 2D case leads to the following
estimate for the phase-separation range

$n_{c}\lesssim n\lesssim\frac{2J}{t}.$ (229)

According to Eq. (228), a nonzero value of the nondiagonal hopping integral
$t^{\prime}$ manifests itself in a nonzero critical concentration $n_{c}$ for
phase separation. The phase separation range obtained by the Maxwell
construction is in fact broader and starts from some $n_{0}<n_{c}$. Here, the
orbital state is described by a classical vector, but nevertheless we have the
lower critical concentration for the canted state and hence for the phase s
eparation, which in the usual double exchange model appears only, if we take
into account the quantum nature of core spins [135]. Inequalities (226) and
(229) are approximate ones and are valid only at quite small values of
parameter $J/t$.

###  7.3 Anisotropic model

Let us now analyze the effects related to the specific form of $e_{g}$
orbitals. To make the underlying physics clearer, we discuss in this
subsection the two-dimensional case (square lattice). This case is relevant,
e.g., to layered cuprates, like K2CuF4, or manganites (La2MnO4 or La2Mn2O7).
Another simplification is that we choose the Heisenberg-like form (215) for
the orbital exchange. Such simplification do not lead to a loss in generality.
At least, in the analysis based on the superexchange mechanism of the orbital
ordering [38] in the case of of $e_{g}$ orbitals demonstrates the the main
physics of the phase separation discussed below remains qualitatively the
same.

As we have mentioned above, any wave function can be written as a linear
combination of basis orbitals. For $e_{g}$ orbitals, we have
$|a>=|x^{2}-y^{2}>$ and $|b>=|2z^{2}-x^{2}-y^{2}>$ and
$|\theta>=\cos(\theta/2)|x^{2}-y^{2}>+\sin(\theta/2)|2z^{2}-x^{2}-y^{2}>$.
Hopping integrals $t^{\alpha\beta}$ in Eq. (216) now become anisotropic and
can be presented in the matrix form

$\displaystyle(t_{x,y})^{\alpha\beta}=\frac{t_{0}}{4}\left(\begin{array}[]{cc}3&\mp\sqrt{3}\\\
\mp\sqrt{3}&1\\\ \end{array}\right)\,,$ (232)

where minus and plus signs correspond to the electron hopping along the axes
$x$ and $y$, respectively.

Again, we choose the the orbital ground state in the form of alternating
$|\theta_{1}\rangle$ and $|\theta_{2}\rangle$ orbitals. In this case, the
spectrum of charge carriers takes the form

$E(\mathbf{k})=-t_{0}\left(A_{x}(\theta_{1},\theta_{2})\cos
k_{x}+A_{y}(\theta_{1},\theta_{2})\cos k_{y}\right)\,,$ (233)

where

$A_{x,y}(\theta_{1},\theta_{2})=\left|\cos\left(\frac{\theta_{1}-\theta_{2}}{2}\right)+\cos\left(\frac{\theta_{1}+\theta_{2}}{2}\pm\frac{\pi}{3}\right)\right|.$
(234)

The corresponding expression for the total energy is

$\displaystyle E_{tot}(\theta_{1},\theta_{2})$ $\displaystyle=$ $\displaystyle
t_{0}\left(A_{x}(\theta_{1},\theta_{2})+A_{y}(\theta_{1},\theta_{2})\right)\times$
(235)
$\displaystyle\varepsilon(n;\theta_{1},\theta_{2})+2J\cos(\theta_{1}-\theta_{2})\,,$

where
$\varepsilon(n;\theta_{1},\theta_{2})=\int_{-1}^{\mu}dE\,E\rho(n;\theta_{1},\theta_{2})$,
and the electronic density of states has the form

$\rho(n;\theta_{1},\theta_{2})=\\!\int\limits_{0}^{\infty}\\!\frac{ds}{\pi}\cos(Es)J_{0}\\!\left(\frac{sA_{x}}{A_{x}+A_{y}}\right)J_{0}\\!\left(\frac{sA_{y}}{A_{x}+A_{y}}\right)\\!.$
(236)

Density of states (236) is the function of orbital angle $\theta_{1}$ and
$\theta_{2}$ . Therefore, to determine the optimum orbital configuration, it
is necessary to minimize $E_{tot}$, Eq. (235), with respect to angles
$\theta_{1}$ and $\theta_{2}$. As a result of such minimization, we find that
at low doping levels ($n$ less than some threshold value, $n_{1}$, which
depends on parameter $J/t_{0}$), the ground state configuration corresponds to
$\theta_{1}=0$, $\theta_{2}=\pi$. This is the homogeneous phase with
alternating $|x^{2}-y^{2}>$ and $|2z^{2}-x^{2}-y^{2}>$ orbitals, i.e. the
antiferro-OO. We do not take in to account here anharmonic terms and higher-
order interactions, which favor locally elongated octahedra corresponding in
our notation to the angles equal to $\theta=\pi$, $\pm 2\pi/3$, see Refs.
[322, 323]). After the minimization, the total energy takes the form

$E_{tot}=t_{0}\sqrt{3}\varepsilon_{0}(n)-2J\,.$ (237)

The plot of this energy has the positive curvature,
$\partial^{2}E_{tot}/\partial n^{2}>0$, corresponding to the stability of such
state.

At $nn_{1}$, the system passes in a jump-like manner to another state
characterized by alternating orbital angles $\theta_{2}=-\theta_{1}$, where

$\theta_{1}=\arccos\left(\frac{t_{0}|\varepsilon_{0}(n)|}{4J}\right)\,,$ (238)

and $E_{tot}(n)$ exhibits a kink at $n=n_{1}$. This canted state arising at
$n>n_{1}$ has the energy

$E_{tot}=t_{0}\varepsilon_{0}(n)-\frac{t_{0}^{2}}{4J}\varepsilon_{0}^{2}(n)-2J\,.$
(239)

The angle $\theta_{1}$ decreases with the increase in $n$ and eventually
vanishes at $n=n_{2}$. The value of $n=n_{2}$ can be found from the equation
$t_{0}|\varepsilon_{0}(n_{2})|/4J=1$. The resulting state,
$\theta_{1}=\theta_{2}=0$, is the ferro-orbital one formed by
$|x^{2}-y^{2}\rangle$ orbitals. The plot total energy versus doping is
presented in in Fig. 48. In this figure, we see that the plot of energy (239)
can exhibit both positive and negative curvature (see the inset of Fig. 48).
For the positive curvature ensuring the stability, the phase separation can
nevertheless arise due to the kink within the doping range close to $n=n_{1}$.
At the negative curvature, the phase separation naturally appears. Thus, both
the negative curvature of $E_{tot}$ plot and the kink can give rise to
different kinds of inhomogeneous states [219].

Figure 48:  (Color online) Energy Eq. (235) versus doping at $J/t_{0}=0.025$
(red solid curve). Within the doping range close $n_{1}\approx 0.08$, the
homogeneous state is unstable, and a phase separation can arise. In the inset,
we plot the function $E_{tot}(n)+t_{0}\sqrt{3}n$ (for the better visuality, we
subtract the linear term in $E_{tot}(n)$) in the vicinity of $n<n_{1}$ at
different values of parameter $J/t$. The red solid curve (blue dashed curve)
corresponds to $J/t_{0}=0.01$ ($J/t_{0}=0.035$), and have a negative
(positive) value of $\partial^{2}E_{tot}/\partial n^{2}$ in the region
$n>n_{1}$ close to $n_{1}$. The phase separation appears in both cases. The
dot-and-dash line is the result of the Maxwell construction [39].

###  7.4 Inhomogeneities in the orbitally ordered structures

In the previous subsections, we have shown that itinerant charge carriers
introduced by doping can favor the formation of charge inhomogeneities in the
case of orbital ordering. Now, let us discuss possible shapes of such
inhomogeneities. We base our discussion on the case of $e_{g}$ electrons
forming a regular array at the square lattice. We assume that the host lattice
is characterized by the antiferro-OO structure with alternating
$|x^{2}-y^{2}>$ and $|2z^{2}-x^{2}-y^{2}>$ orbitals considered in the previous
subsection. An additional charge carrier favors the formation of another type
of orbital ordering within a limited region close to it. The optimum structure
can be found by minimization of the additional energy related to the doping.
Without a loss in generality, we can consider alternating $|\theta_{1}\rangle$
and $|\theta_{2}\rangle$ orbitals.

In such a situation, the energy spectrum of charge carriers is specified by
Eq. (233). Since we are dealing with a limited region, we can write an
effective Hamiltonian by expanding Eq. (233) up to up to the second order in
$\mathbf{k}$. As a result, we have

$\hat{H}_{eff}=-t_{0}\left(A_{x}+A_{y}\right)+\frac{t_{0}}{2}\left(A_{x}\frac{\partial^{2}}{\partial
x^{2}}+A_{y}\frac{\partial^{2}}{\partial y^{2}}\right),$ (240)

where $A_{x}$, $A_{y}$ are defined in Eq. (234). Then, we assume that the
region under study is an ellipse with semiaxes $\sqrt{A_{x}}\rho_{0}$ and
$\sqrt{A_{y}}\rho_{0}$ and search for the solution of the Schrödinger equation
within this ellipse employing Hamiltonian (240). In such elliptical droplet,
the kinetic energy of a charge carrier can be written as

$E_{kin}=-t_{0}\left(A_{x}+A_{y}\right)+\frac{t_{0}j_{0,1}^{2}}{2\rho_{0}^{2}}\,.$
(241)

Here, $j_{0,1}\cong 2.405$ is the first root of the Bessel function $J_{0}$.

For the orbitally ordered state, we can write the potential energy $E_{pot}$
related to the orbital ordering, which is proportional to the volume $v$ of
the elliptical droplet ($v=\pi\sqrt{A_{x}A_{y}}\rho_{0}^{2}$). It includes the
energy of the new canted OO state $zJ\cos(\theta_{1}-\theta_{2})/2$ and the
loss in energy of the antiferro-OO host structure occurring due to arising
droplets with the modified OO structure ($zJ/2$, $z=4$). For a set of
identical droplets, the sum of these two contributions is

$E_{pot}=4\pi\rho_{0}^{2}J\sqrt{A_{x}A_{y}}\cos^{2}\left(\frac{\theta_{1}-\theta_{2}}{2}\right)\,.$
(242)

Minimizing the droplet energy $E_{kin}+E_{pot}$ with respect to $\rho_{0}$, we
find

$\rho_{0}=\left(\frac{t_{0}j_{0,1}^{2}}{8\pi
J\sqrt{A_{x}A_{y}}\cos^{2}\left(\frac{\theta_{1}-\theta_{2}}{2}\right)}\right)^{1/4}\\!\\!.$
(243)

The total energy (per lattice site) then reads

$\displaystyle E_{tot}$ $\displaystyle=$
$\displaystyle-2J+E(\theta_{1},\,\theta_{2})\,n\,,$ (244) $\displaystyle
E(\theta_{1},\,\theta_{2})$ $\displaystyle=$
$\displaystyle-2t_{0}(A_{x}+A_{y})\phantom{\frac{1}{2}}+j_{0,1}\left(8\pi
t_{0}J\sqrt{A_{x}A_{y}}\right)^{1/2}\left|\cos\left(\frac{\theta_{1}-\theta_{2}}{2}\right)\right|\,.$
(245)

To determine the optimum shapes of OO droplets, we perform the minimization of
energy(244) with respect to orbital angles $\theta_{1}$ and $\theta_{2}$ using
expressions (234)[for $A_{x,y}$. The behavior of$E(\theta_{1},\,\theta_{2})$
at two different values of $J/t_{0}$ is illustrated in Fig 49. In this figure,
we can clearly see that the function $E(\theta_{1},\,\theta_{2})$ exhibits
several minima, and the specific angles $\theta_{1}$ and $\theta_{2}$
corresponding to the lowest energy strongly depend on the value of $J/t_{0}$.
At small $J/t_{0}$ [Fig 49(a)], the lowest energy is achieved at to
$\theta_{1}=\theta_{2}=0$ (ferro-OO state with the filled $|x^{2}-y^{2}>$
orbitals). For such state, the optimum shape of droplets is a circle (see the
left panel of Fig. 50). At $J/t_{0}$ exceeding a certain critical value
($J_{cr}/t_{0}\simeq 0.0075$), the minimum $\theta_{1}=\theta_{2}=0$
corresponds to a metastable state, and the energy $E(\theta_{1},\,\theta_{2})$
attains the lowest value at four degenerate structures: with the alternation
of either $\theta_{1}=\pi/3$ and $\theta_{2}=2\pi/3$, or $\theta_{1}=-\pi/3$
and $\theta_{2}=-2\pi/3$, and the other two minima correspond to the similar
structures with the replacement $\theta_{1}\leftrightarrow\theta_{2}$ [see Fig
49(b)]. Such four orbital arrays are the chains with alternating
$|x^{2}-z^{2}>$ and $|2y^{2}-x^{2}-z^{2}>$ (or $|y^{2}-z^{2}>$ and
$|2x^{2}-y^{2}-z^{2}>$) orbitals. This, in fact, leads to the formation of
nearly one-dimensional (cigar-shaped) droplets directed along $y$ or $x$ axes
(right panel of Fig. 50). At the further increase in $J/t_{0}$, the metastable
$\theta_{1}=\theta_{2}=0$ array splits into two more colicated structures with
$\theta^{*}_{1}=-\theta^{*}_{2}$, in which $\theta^{*}_{1}$ can be either
positive and (see Fig 49(b)). The inhomogeneities related to the latter arrays
should have the circular shape, but the OO structure inside them is canted.

(a) $J/t_{0}=0.005$ (b) $J/t_{0}=0.02$

Figure 49: Energy Eq. (244) of an orbital versus angles $\theta_{1}$ and
$\theta_{2}$ at $J/t_{0}$ smaller (a) and larger (b) than the critical value
$\simeq 0.0075$ [39].

Figure 50:  (Color online) Possible types of orbital inhomogeneities arising
due to doping: circular (left panel) and needle-like (right panel) droplets.
The motion of an itinerant charge carrier results in the formation of a ferro-
or canted-OO structure within a limited region inside the host material with
the antiferro-OO. [39].

The physical mechanism underlying the formation of two types of droplets with
different shapes is quite clear. Indeed, the ferro-OO state with
$|x^{2}-y^{2}>$ orbitals provides the maximum gain in the kinetic energy. This
gain prevails at small $J/t_{0}$ favoring circular droplets with such the
ferro-OO. In the opposite case of dominant potential energy (large $J/t_{0}$),
nearly one-dimensional structures turn out to be favorable. Therefore, the
optimum droplets have now cigar- or needle-like shape. The orbital
configuration in a needle-like droplet provides the maximum gain in the
kinetic energy in the absence of hopping between neighboring chains.

It can be shown that in a certain doping range, $0<n<n^{*}_{c}$, the
inhomogeneous state including circular or needle-like OO droplets within the
antiferro-OO host material corresponds to the energy lower than that of the
homogeneous state. With the increase in the charge carrier density, the
droplets can overlap; therefore, at some threshold value $n=n^{*}_{c}$, the
the considered droplet type disappears. However, the phase separation can
exist in a broader doping range (see the previous subsection). For the
threshold charge carrier density in the case of circular droplets, we obtain
an estimate $n^{*}_{c}\sim 1/\pi\rho_{0}^{2}$. Thus, at $J/t_{0}=0.005$, we
have $\rho_{0}\approx 2$ (in units of lattice constant) and $n^{*}_{c}\sim
0.08$.

For the needlelike droplets, $A_{x}=0$ (or $A_{y}=0$), and Eq. (243) gives us
the chains of infinite length (but with zero volume $v$), $\rho_{0}=\infty$.
This rather unphysical result stems from the assumption that the potential is
proportional only to the droplet volume only. To estimate the characteristic
length $L$ of the needle, we should consider also the surface contribution to
the energy (proportional to the length of droplet) to the potential energy
$E_{pot}$. For such an estimate, we can consider, for example, an array of
$|y^{2}-z^{2}>$ and $|2x^{2}-y^{2}-z^{2}>$ orbitals stretched along the $x$
axis. For this array, we have $A_{x}=\sqrt{3}$, $A_{y}=0$. The corresponding
effective Hamiltonian (240) takes the form

$\hat{H}_{eff}=-t_{0}\sqrt{3}+\frac{t_{0}\sqrt{3}}{2}\frac{\partial^{2}}{\partial
x^{2}}\,,$

and the kinetic energy of a charge carrier in the chain of length $L$ can be
written as $E_{kin}=-t_{0}\sqrt{3}(1-\pi^{2}/2L^{2})$. The surface energy
results from the interaction between the orbitals in the chain and in the host
material. It The minimum surface energy corresponds to the configuration
illustrated in Fig. 50, for which each $|y^{2}-z^{2}>$
($|2x^{2}-y^{2}-z^{2}>$) orbital in the chain has its nearest neighbor
$|x^{2}-y^{2}>$ ($|2z^{2}-x^{2}-y^{2}>$) orbital in the host. Using the
continuum approximation, we can write the potential energy as $E_{pot}=9JL/4$.
Then, we minimize the total energy $E_{kin}+E_{pot}$ with respect to $L$, and
obtain the expression for the characteristic length of a needle-like droplet

$L_{0}=\left(\frac{4\pi^{2}t_{0}\sqrt{3}}{9J}\right)^{1/3}\\!\\!.$ (246)

For$J/t_{0}=0.05>J_{cr}/t_{0}$, we obtain $L_{0}\approx 5.5$. If we assume a
random distribution of the a needle-like droplet in host (with the needles
directed both along $x$ and $y$ axes), the critical density is about
$n^{*}_{c}\sim 1/L_{0}^{2}$. Note, however, that this density can be higher in
a more regular array such as that corresponding to stripes.

Let us recall that for spin polarons discussed in the previous sections, the
existence of the gapless (Goldstone) excitations leads to slowly decaying spin
structure distortions around magnetic defect (magnetic polaron) [120].
Therefore, it is natural to expect that orbital polarons should give rise to a
long-range strain fields will be created around them. In analogy to the case
of spin polarons, these long range fields should not affect the existence of
orbital polarons themselves. They can only modify the parameters of
inhomogeneous state and lead to an additional long-range interaction between
inhomogeneities.

In this and in the previous subsections, we were dealing with a two-
dimensional case of the anisotropic model. However, in contrast to the
isotropic case treated in subsection 7.2, the space dimensionality can
drastically affect the obtained results. Indeed, the hopping integrals for
$e_{g}$ electrons strongly anisotropic. This can produce a pronounced effect
on the OO type in the homogeneous state as well as on the shape and orbital
order of the inhomogeneities. That is why, the detailed study of the three-
dimensional case is a separate and quite complicated task.

###  7.5 Simplified two-band model with the orbital ordering

####  7.5.1 Formulation of the problem

In the previous subsections, we studied the effect of doping on already
existing orbital ordering. Now, we shall treat a more natural situation when
the same charge carriers are involved both in the orbital ordering and in
charge transfer. We start with the analysis of the system with JT ions having
double-degenerate ground state. As usual, such degeneracy is lifted by local
lattice distortions. Then, each lattice site, can correspond to to two
different ground states, $a$ or $b$ (in perovskites, for example $a$ ($b$)
state can be related to stretching (compression) of the surrounding anion
octahedron). The $a$ and $b$ states at site $\mathbf{n}$ are, in fact,
possible orbitals filled by a charge carrier. Of course, the local distortions
are not independent, and their interaction can lead to certain regular
structure. The simplest Hamiltonian describing the orbital interaction can be
described in an Ising-like form [324]

$H_{\text{OO}}=J\sum_{\langle{\bf n}{\bf m}\rangle}\tau^{z}_{\bf
n}\tau^{z}_{\bf m}\,.$ (247)

Here, the $\tau^{z}_{\bf n}$ operators have the eigenvalues $1$ or $-1$
depending on the state of site ${\bf n}$ ($a$ or $b$, respectively). In such
system, we can consider two simplest types of the orbital ordering: either
ferro-OO (with the same orbitals occupying each site) and antiferro-OO
(corresponding to the alternation of orbitals).

The doping gives rise to nonzero density, $n\neq 0$, of itinerant charge
carriers appear in the system. As a result, the system includes localized
charge carriers responsible for OO and itinerant ones. The electron motion is
determined by hopping integrals depending on the states of the neighboring
lattice sites. In the model under study, there are three possible hopping
integrals, namely, $t^{aa}$, $t^{bb}$, and $t^{ab}=t^{ba}$. For simplicity,
let us assume that $t^{aa}=t^{bb}>t^{ab}$. At both types of orbital ordering
(ferro- and antiferro-OO) the charge carriers are described by the Hamiltonian

$H_{\text{el}}=-t\sum_{\langle{\bf n}{\bf
m}\rangle,\sigma}\left(a^{{\dagger}}_{{\bf n}\sigma}a_{{\bf
m}\sigma}+h.c.\right)+\frac{U}{2}\sum_{{\bf n}\sigma}n_{{\bf n}\sigma}n_{{\bf
n}\bar{\sigma}}\,.$ (248)

The first term in $H_{\text{el}}$ describes the kinetic energy of itinerant
charge carriers. Symbol $t$ denotes the hopping integral $t^{aa}$ in for the
ferro-OO and $t^{ab}$ for antiferro-OO. The second (Hubbard) term in Eq. (248)
takes into account the onsite Coulomb repulsion between electrons and the bar
above $\sigma$ means the spin component of the opposite sign. Angular brackets
$\langle...\rangle$ denote the summation over nearest neighbor sites. We also
assume that the parameters of Hamiltonian (248) are independent of the doping
level $n$.

####  7.5.2 Homogeneous states and phase separation

Hamiltonian (248), is the usual Hubbard Hamiltonian for charge carriers moving
on the ferro- or antiferro-OO background. We begin our analysis from the
possible homogeneous states corresponding to such Hamiltonian. At $J>0$ and in
the absence of doping($n=0$), we have the uniform antiferro-OO ground state,
which is the ground state if $J>0$ (the case of ferro-OO state can be treated
in a similar way). In the antiferro-OO state we can consider only the
interorbital hopping $t=t^{ab}$ in Eq. (248). We determine the energy spectrum
of charge carriers using the Hubbard I approximation [23] in equation of
motion for one-electron Green’s function $G_{\sigma}({\bf
n,n}_{0};\,t-t_{0})=-i\langle Ta_{{\bf n}\sigma}(t)a^{{\dagger}}_{{\bf
n}_{0}\sigma}(t_{0})\rangle$. In the frequency–momentum representation, this
Green’s function can be written as [23]

$G_{\sigma}({\bf k},\omega)=\frac{\omega+\mu-U(1-n/2)}{\left(\omega+\mu-
E_{1}({\bf k})\right)\left(\omega+\mu-E_{2}({\bf
k})\right)}=\\!\frac{\omega\\!+\\!\mu\\!-\\!U(1\\!-\\!n/2)}{E_{1}({\bf
k})-E_{2}({\bf k})}\\!\left\\{\\!\\!\frac{1}{\omega\\!+\\!\mu\\!-\\!E_{1}({\bf
k})}\\!-\\!\frac{1}{\omega\\!+\\!\mu\\!-\\!E_{2}({\bf
k})}\\!\\!\right\\}\\!\\!,$ (249)

where

$E_{1,2}({\bf k})=\frac{U+\varepsilon({\bf
k})}{2}\mp\sqrt{\left(\frac{U-\varepsilon({\bf
k})}{2}\right)^{2}+U\varepsilon({\bf k})\frac{n}{2}}\,,$ (250)

and $\varepsilon({\bf k})$ is the ${\bf k}$-dependence of electron energy in
the absence of correlations ($U=0$). We do not consider any spin ordering and
assume that the numbers of electrons with each spin projections are equal,
$n_{\sigma}=n_{-\sigma}=n/2$. In the case of simple cubic lattice, we have
$\varepsilon({\bf k})=-2t^{ab}(\cos k_{x}+\cos k_{y}+\cos k_{z})$, where for
$k_{i}$, we use the the units of inverse lattice constant $d$.

The charge carrier density can be determined using Green’s functions (249)

$Cn=-i\int\limits\frac{d\omega}{\pi}\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}G_{\sigma}({\bf k},\omega+i0\;\textrm{sign}\omega)e^{i\omega
0}.$ (251)

Changing in (249) $\omega$ by $\omega+\mu$ and using spectrum (250), we find

$n=2\int_{\tilde{w}}^{\mu}d\omega\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}\frac{U\left(1-n/2\right)-\omega}{\sqrt{(U-\varepsilon({\bf
k}))^{2}+2\varepsilon({\bf k})nU}}\left\\{\delta\left(\omega-E_{1}({\bf
k})\right)-\delta\left(\omega-E_{2}({\bf k})\right)\right\\}\,,$ (252)

where $\tilde{w}=\tilde{w}(n)$ is the bottom of the lower band. When the lower
band is not fully occupied filled, the second $\delta$-function in (252) does
not work. contribute to the result. In fact, Eq. (252) allows us to find the
the dependence of chemical potential on the doping level $n$. The, the kinetic
energy per lattice sire has the form

$E_{\text{kin}}=2\int_{\tilde{w}}^{\mu}\omega d\omega\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}\frac{U\left(1-n/2\right)-\omega}{\sqrt{(U-\varepsilon({\bf
k}))^{2}+2\varepsilon({\bf k})nU}}\left\\{\delta\left(\omega-E_{1}({\bf
k})\right)-\delta\left(\omega-E_{2}({\bf k})\right)\right\\}\,,$ (253)

After that, we can calculate the total energy of the homogeneous antiferro-OO
state

$E_{AFOO}=E_{kin}(t^{ab})-Jz/2.$ (254)

The plot of such energy as function of doping $n$ shown in Fig. 51, where it
is compared to the corresponding plot for the ferro-OO state (see below)

Figure 51:  (Color online) Dependence of the energies of homogeneous ferro-OO
(blue solid line) and antiferro-OO (red dashed line) states onthe doping
level$n$. Green dashed lines describe the energy of an inhomogeneous state.
The plots correspond to the following choice of parameters:
$t^{ab}/t^{aa}=0.3$, $U/t^{aa}=3$, and $J/t^{aa}=0.1$ [324].

In the case of a homogeneous ferro-OO state (the charge carriers occupy
similar orbitals), the energy loss per site related to the unfavorable orbital
ordering is $zJ/2$. However, the gain in kinetic energy is larger because the
hopping integral between similar orbitals exceeds that between different
orbitals ($t^{aa}>t^{ab}$). Here, the kinetic energy is described by the same
expression as (253) but with the different hopping integral, $t^{aa}$. Thus,
instead of Eq. (254), we have

$E_{FMOO}=E_{kin}(t^{aa})+Jz/2.$ (255)

In Fig. 51, we illustrate the behavior of energies for the system in ferro-OO
and antiferro-OO states within the $0<n<1$ doping range. We can clearly see
that at some dopings, antiferro-OO is unfavorable and the ferro-OO becomes the
ground state owing to a larger kinetic energy gain. In the vicinity of the
crossover range, we can see that the homogeneous state turns out to be
unstable toward the phase separation into a mixture of ferro-OO and antiferro-
OO domains characterized by different values of charge carrier density (see
green dashed lines in Fig. 51).

###  7.6 A model with charge carriers occupying $e_{g}$ orbitals: 2D case

####  7.6.1 A model

Let us now generalize the previous analysis to the case of some specific
orbitals such as $e_{g}$ orbitals on the square 2D lattice [324]. This
situation is relevant, e.g. to such compounds as layered cuprates,( K2CuF4,)
or manganites (La2MnO4 and La2Mn2O7). As it was discussed above, any $e_{g}$
orbital can be represented as a linear combination of two basis states
$|a>=|x^{2}-y^{2}>$ and $|b>=|2z^{2}-x^{2}-y^{2}>$

$|\theta>=\cos\frac{\theta}{2}\,|x^{2}-y^{2}>+\sin\frac{\theta}{2}\,|2z^{2}-x^{2}-y^{2}>.$
(256)

In the case under study, the same charge carries are responsible both for the
formation of orbital order and for the charge transfer. This means that the
Hamiltonian of the model should involve a, all possible types of of the
intersite hopping integrals. Of course, such hopping integrals
$t^{\alpha\beta}$ ($\alpha,\beta=a,b$) should depend on the direction of of
hopping. We can write the corresponding Hamiltonian in the form

$H_{el}=-\sum_{\langle{\bf
nm}\rangle\alpha\beta\sigma}t_{\mathbf{n}-\mathbf{m}}^{\alpha,\beta}d^{{\dagger}}_{{\bf
n}\alpha\sigma}d_{{\bf m}\beta\sigma}+\frac{U}{2}\sum_{{\bf
n}\alpha\sigma}n_{{\bf n}\alpha\sigma}n_{{\bf
n}\alpha\bar{\sigma}}+\frac{U}{2}\sum_{{\bf
n}\alpha\sigma\sigma^{\prime}}n_{{\bf n}\alpha\sigma}n_{{\bf
n}\bar{\alpha}\sigma^{\prime}}\,,$ (257)

where for the the hopping integrals, we have
$t_{\mathbf{n}-\mathbf{m}}^{\alpha\beta}=t_{x}^{\alpha\beta}$
($t_{\mathbf{n}-\mathbf{m}}^{\alpha\beta}=t_{y}^{\alpha\beta}$) for
$\mathbf{n}-\mathbf{m}$ directed along the $x$ ($y$) axis. As usual, in Eq.
(257), $\bar{\sigma}$ and $\bar{\alpha}$ denote not $\sigma$ and not $\alpha$,
respectively. The hopping integrals $t_{x,y}^{\alpha\beta}$ can be represented
in the matrix form

$\displaystyle(t_{x,y})^{\alpha\beta}=\frac{t_{0}}{4}\left(\begin{array}[]{cc}3&\mp\sqrt{3}\\\
\mp\sqrt{3}&1\\\ \end{array}\right)\,.$ (260)

Here, signs minus and plus correspond to $x$ and $y$ hopping directions,
respectively. Further on, we limit ourselves to the case of large $U/t_{0}$ in
Hamiltonian (257).

####  7.6.2 Antiferro-orbital ordering near half filling

The case of one electron per site (half filling) is of special interest. Here,
we study the situation near half filling, $x=1-n\ll 1$. It can be demonstrated
that at $n=1$ and large $U$, Hamiltonian (257) has the insulating the ground
state with the antiferro-orbital ordering. Indeed, the second order of
perturbation theory in terms of $t_{0}/U$ gives rise to the effective
Hamiltonian in the form [38]

$\displaystyle H_{OO}$ $\displaystyle=$ $\displaystyle
J\sum_{\langle\mathbf{nm}\rangle}\left\\{\mathbf{s_{n}s_{m}}\left[\tau^{z}_{\mathbf{n}}\tau^{z}_{\mathbf{m}}-2\tau^{z}_{\mathbf{n}}+1\pm
2\sqrt{3}\tau^{z}_{\mathbf{n}}\tau^{x}_{\mathbf{m}}\mp
2\sqrt{3}\tau^{x}_{\mathbf{m}}+3\tau^{x}_{\mathbf{n}}\tau^{x}_{\mathbf{m}}\right]\right.$
(262)
$\displaystyle\left.+\frac{1}{2}\left[\tau^{z}_{\mathbf{n}}\tau^{z}_{\mathbf{m}}+2\tau^{z}_{\mathbf{n}}\pm{\sqrt{3}\tau^{x}_{\mathbf{m}}}2\sqrt{3}\tau^{z}_{\mathbf{n}}\tau^{x}_{\mathbf{m}}\pm
2\sqrt{3}\tau^{x}_{\mathbf{m}}+3\tau^{x}_{\mathbf{n}}\tau^{x}_{\mathbf{m}}\right]\right\\}\,,$

where $J=t_{0}^{2}/U$,
$\mathbf{s_{n}}=\sum_{\alpha\sigma\sigma^{\prime}}d^{{\dagger}}_{{\bf
n}\alpha\sigma}(\bm{\sigma})_{\sigma\sigma^{\prime}}d_{{\bf
n}\alpha\sigma^{\prime}}/2$ is the operator determining electron spin at site
$\mathbf{n}$, and
$\tau^{z,x}_{\mathbf{n}}=\sum_{\alpha\beta\sigma}d^{{\dagger}}_{{\bf
n}\alpha\sigma}(\sigma^{z,x})_{\alpha\beta}d_{{\bf n}\beta\sigma}/2$ describe
the occupation of orbitals, and $\bm{\sigma}$ are the Pauli matrices.
Similarly to Eq. (260), signs minus and plus correspond to $x$ and $y$
directions of the $\mathbf{n}-\mathbf{m}$ bond, respectively. In the absence
of spin ordering, $\mathbf{s_{n}s_{m}}=0$, and the above Hamiltonian takes a
simpler form

$H_{OO}=\frac{J}{2}\sum_{\langle\mathbf{nm}\rangle}\left[\tau^{z}_{\mathbf{n}}\tau^{z}_{\mathbf{m}}+2\tau^{z}_{\mathbf{n}}\pm
2\sqrt{3}\tau^{z}_{\mathbf{n}}\tau^{x}_{\mathbf{m}}\pm
2\sqrt{3}\tau^{x}_{\mathbf{m}}+3\tau^{x}_{\mathbf{n}}\tau^{x}_{\mathbf{m}}\right]\,.$
(263)

Hamiltonian Eq. (263) has the antiferro-OO ground state corresopnding to the
alternation of $|\theta_{0}>$ and $|-\theta_{0}>$ orbitals, where

$\theta_{0}=\arccos\left(-\frac{1}{4}\right)\,.$ (264)

The orbital ordering leads to the energy gain per lattice site
$E_{OO}=-13J/4$.

At the hole doping, there appear some empty lattice sites. It is natural to
assume that the long-range antiferro-orbital order is not destroyed at low
doping levels, $x=1-n\ll 1$, but the energy gain related to the ordering
becomes smaller. In the framework of the mean-field approximation, we can
write $\langle\tau^{z,x}_{\mathbf{n}}\rangle\propto 1-x$ and

$E_{OO}=-\frac{13+11x^{2}}{4}J+6J\,x\,.$ (265)

To determine the kinetic energy of mobile holes in antiferro-OO structure, we
can define the creation $b^{{\dagger}}_{\mathbf{n}\sigma}$ and annihilation
$b_{\mathbf{n}\sigma}$ operators for such holes

$b_{\mathbf{n}\sigma}=\cos\frac{\theta_{0}}{2}\,d^{{\dagger}}_{na\sigma}\pm\sin\frac{\theta_{0}}{2}\,d^{{\dagger}}_{nb\sigma}\,,\;\;b^{{\dagger}}_{\mathbf{n}\sigma}=(b_{\mathbf{n}\sigma})^{{\dagger}}\,,$
(266)

where signs plus and minus sign correspond to the sites with $|\theta_{0}>$
and $|-\theta_{0}>$ orbitals, respectively. Then, the kinetic energy described
by the first term in Hamiltonian (257)) can be written as

$H_{el}=\frac{t_{0}}{8}\sum_{\langle{\bf n}{\bf m}\rangle\sigma}\left(P_{{\bf
n}\sigma}^{{\dagger}}b^{{\dagger}}_{{\bf n}\sigma}b_{{\bf m}\sigma}P_{{\bf
m}\sigma}+h.c.\right),$ (267)

where we projection operators $P$ are introduced to exclude a double
occupation of each sites. Let us analyse the situation with the 2D square
lattice. Then, the energy spectrum of free holes is

$\varepsilon({\bf k})=-\frac{t_{0}}{4}\left(\cos k_{x}+\cos k_{y}\right)\,.$
(268)

In this case, the energy of of the doped antiferro-OO array can be written as

$E_{AFOO}(x)=-\frac{13+11x^{2}}{4}J+\int_{-t_{0}/2}^{\mu(x)}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!dE\,E\rho_{0}(E)+6Jx\,,$
(269)

where the hole density of states $\rho_{0}(E)$ for holes in given by the
expression

$\rho_{0}(E)=\\!\int\\!\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}\delta(E-\varepsilon(\mathbf{k}))=\\!\int\limits_{0}^{\infty}\\!\frac{ds}{\pi}\cos(Es)J_{0}^{2}\\!\left(\frac{t_{0}s}{4}\right)\\!.$
(270)

The chemical potential $\mu$ is determined by the condition
$x=\int_{-t_{0}/2}^{\mu(x)}dE\rho_{0}(E)$. The calculated energy (269) versus
$n=1-x$ is plotted in Fig. 52, see the blue dashed curve and compared to the
energy of the orbitally disordered state analyzed below.

Figure 52:  (Color online) Dependence of energies on the total number $n$ of
electrons for of orbitally disordered (red solid curve) and antiferro-OO
(dashed blue curve) states. $J/t_{0}=0.08$. The green dot-ans-dash line shows
the energy of a possible phase-separated state [324].

####  7.6.3 Orbitally disordered state

Now, let us discuss the situation characteristic of relatively high doping
levels $x$, the the orbital ordering discussed above disappears. We shall
consider the state corresponding to Hamiltonian (257) in the absence of the
long-range OO at the electron densities $0<n<1$ and compare the corresponding
energy (269). The one-particle Green’s have the form similar to that
considered above in subsubsection 7.5.2 and analyze them using the Hubbard I
approximation [23]. The corresponding equations of motion can be written as

$(\omega+\mu)G_{\alpha\beta}({\bf
n},\omega)=g_{\alpha}\left(1-\sum_{{\bf\Delta}\gamma}t_{\bf\Delta}^{\alpha\gamma}G_{\gamma\beta}({\bf
n+\Delta},\omega)\right)\,,$ (271)

where

$g_{\alpha}=1-\frac{n_{\alpha}}{2}-n_{\bar{\alpha}},$ (272)

and the summation is taken over sites $\mathbf{m}$ , wihich are the neighbors
of the chosen site $\mathbf{n}$, ($\bm{\Delta}=\mathbf{n}-\mathbf{m}$). We
assume the absence of magnetic ordering,
$n_{\alpha\uparrow}=n_{\alpha\downarrow}=n_{\alpha}/2$ and work at in the
limit of strong correlations, $U\rightarrow\infty$.

Passing to the momentum representation, we can write Eq.(271) as

$(\omega+\mu)G_{\alpha\beta}({\bf
k},\omega)=g_{\alpha}\left(1-\sum_{\gamma}\varepsilon^{\alpha\gamma}({\bf
k)}G_{\gamma\beta}(\omega,{\bf k})\right)\,,$ (273)

where $\varepsilon^{\alpha,\beta}({\bf k)}$ is represented in the matrix form

$\hat{\varepsilon}({\bf k})=\frac{-t_{0}}{2}\left(\begin{array}[]{cc}3(\cos
k_{x}+\cos k_{y})&\sqrt{3}(\cos k_{y}-\cos k_{x})\\\ \sqrt{3}(\cos k_{y}-\cos
k_{x})&\cos k_{x}+\cos k_{y}\\\ \end{array}\right).$ (274)

Then, introducing the matrix

$\hat{S}=\sqrt{\hat{g}}=\left(\begin{array}[]{cc}\sqrt{g_{a}}&0\\\
0&\sqrt{g_{b}}\\\ \end{array}\right)\,,$ (275)

it is possible to represent the Green’s function in the compact operator form

$\hat{G}=\hat{S}\frac{1}{\omega+\mu-\hat{\varepsilon}_{g}}\hat{S}\,,$ (276)

where

$\hat{\varepsilon}_{g}=\left(\begin{array}[]{cc}g_{a}\varepsilon_{aa}&\sqrt{g_{a}g_{b}}\varepsilon_{ab}\\\
\sqrt{g_{a}g_{b}}\varepsilon_{ab}&g_{b}\varepsilon_{bb}\\\
\end{array}\right).$ (277)

The energy spectrum is determined by the solution of equation

$\det(\hat{\varepsilon}_{g}-E)=0.$ (278)

Eventually, we find

$E_{\pm 1}({\bf
k})=-\frac{t_{0}(3g_{a}+g_{b})}{4}\left[\phantom{\sqrt{\frac{1}{2}}}\\!\\!\\!\\!\\!\\!(\cos
k_{x}+\cos k_{y})\pm\sqrt{(\cos k_{x}+\cos
k_{y})^{2}-\frac{48g_{a}g_{b}}{(3g_{a}+g_{b})^{2}}\cos k_{x}\cos
k_{y}}\,\right]\,.$ (279)

Using the quasiparticle spectrum (279), we can write the expression for the
Green’s function

$G_{\alpha\beta}(\omega,\mathbf{k})=\sum_{s=\pm
1}\frac{\sqrt{g_{\alpha}g_{\beta}}v^{(s)}_{\alpha}({\bf
k})v^{(s)}_{\beta}({\bf k})}{\omega+\mu-E_{s}({\bf k})}\,.$ (280)

In Eq. (280), $v^{(s)}_{\alpha}$ are the matrix elements of the transformation
determined by the $\hat{v}^{(s)}$ obeying the relations

$\left[\hat{\varepsilon}_{g}({\bf k})\hat{v}^{(s)}\right]_{\alpha}=E_{s}({\bf
k})v^{(s)}_{\alpha}\,.$ (281)

As a result of straightforward but cumbersome calculations, we find an
explicit form of matrix $\hat{v}^{(s)}$ corresponding to the basis orbitals
$|b>=|2z^{2}-x^{2}-y^{2}>$ and $|a>$ = $|x^{2}-y^{2}>$

$\hat{v}^{(s)}_{\alpha}=\frac{1}{\sqrt{(g_{a}\varepsilon_{aa}-E_{s})^{2}-g_{a}g_{b}\varepsilon^{2}_{ab}}}\left(\begin{array}[]{cc}-\sqrt{g_{b}g_{a}}\varepsilon_{ab}\\\
g_{a}\varepsilon_{aa}-E_{s}\\\ \end{array}\right),$ (282)

where $\varepsilon_{aa},\varepsilon_{ab}$, and $\varepsilon_{ab}$ are matrix
elements of matrix $-\hat{t}({\bf k})$, see (274).

Using the above expressions, we can calculate the density of states
corresponding to two new bands with indices $s=1$ and -1

$\rho^{(s)}(E)=\sum_{\alpha}\rho^{(s)}_{\alpha\alpha}(E),$ (283)

where

$\rho^{(s)}_{\alpha\alpha}(E)=\frac{1}{(2\pi)^{2}}\int d{\bf
k}\left[v^{(s)}_{\alpha}(\mathbf{k})\right]^{2}\delta(E-E_{s}({\bf k})).$
(284)

Then, we can determine the partial filling of orbital states

$n_{\alpha}=2g_{\alpha}\sum_{s=\pm
1}\int_{-\infty}^{\mu}\rho^{(s)}_{\alpha\alpha}(E)dE.$ (285)

The chemical potential $\mu$ as can be obtained using the condition

$n_{a}+n_{b}=1-x,$ (286)

where $x$ is the hole density (doping level). Substituting $\mu(n_{a})$ to Eq.
(285), we can determine the dependence of $n_{a}$ and $n_{b}$ on the doping
level $x$. Then, the total energy corresponding to Hamiltonian (248) is given
by the expression

$E_{el}=\sum_{\alpha,s=\pm
1}(1+g_{\alpha})\int_{-\infty}^{\mu}\rho^{(s)}_{\alpha\alpha}(E)EdE.$ (287)

In Fig. 53 we illustrate the behavior of density of states for the new bands
$s=\pm 1$, These densities of states exhibit the logarithmic van Hove
singularity, which is related to the transition from the closed to open Fermi
surface in each bands (the evolution of the Fermi surface is represented in
Fig. 55).

Figure 53:  Densities of states for bands $s=-1$ (red solid line) and $s=1$
(blue solid line). The occupied states at the total number of electrons per
site $n=0.71$ are shown by the shaded area [324].

The behavior of occupation numbers $n_{a}$ and $n_{b}$ for $|a>$ and $|b>$
orbital states with the variation of the total number $n$ of electrons per
site is illustrated in Fig. 54. With the increase in $n$, first only
$|a>=|x^{2}-y^{2}>$ orbitals turn out to be occupied. This means that at high
doping levels $1-n$, we have the occupation of only one sort of orbitals (a
kind of ferro-OO) and these $|x^{2}-y^{2}>$ orbitals provide the largest gain
in the kinetic energy of electrons. The further increase in $n$ gives rise to
electrons at $|b>=|2z^{2}-x^{2}-y^{2}>$ orbitals. Eventually, at $n=1$, we get
$n_{a}=1/4=\sin^{2}(\pi/6)$ and $n_{b}=3/4=\cos^{2}(\pi/6)$. This corresponds
to the average occupation of orbitals at $n=1$ equal to $|\theta=\pm\pi/3>$.
However, the 0btained results on the density of states and on the average
occupation numbers $n_{a}$ and $n_{b}$ do not provide a conclusive evidence
concerning the actual alteration of the $|\pm\theta>$ orbitals.

Figure 54:  (Color online) Average occupation numbers of $|a>$ (red solid
line) and $|b>$ (blue dashed line) orbitals versus total number $n$ of
electrons per site. A kink near $n\simeq 0.6$ corresponds to the doping level
at which the chemical potential crosses the bottom of the upper renormalized
band with $s=1$ [324].

Evolution of the Fermi surface with the increase in the total number $n$ of
electrons is presented in Fig. 55.

(a) $n=0.15$ (b) $n=0.36$ (c) $n=0.5$ (d) $n=0.68$

(e) $n=0.8$ (f) $n=0.98$

Figure 55: Evolution of the Fermi surface with the increase in the total
number of electrons $n$. The red solid curves indicate the Fermi surface sheet
corresponding to the $s=-1$ band and the blue solid curves correspond to the
$s=1$ band [324].

At $n=0.8$ [see Fig. 55(e)] the Fermi surface involves several flat portions.
This quite interesting situation. In Fig 54, we see that $n=0.8$ corresponds
to $n_{a}=n_{b}=0.4$, and hence, $g_{a}=g_{b}=0.4$. Substituting $g_{a}=g_{b}$
into Eqs. (279) and (284), we can find out that in such particular case, the
quasiparticle spectrum and and the corresponding density of states coincide
with those for uncorrelated electrons up to a numerical factor. This means
that for $n_{a}=n_{b}$, we are dealing with the compensation of two main
correlation effects, namely, the changes in bandwidths and in number of states
within the bands. Thus, we have a clear similarity to the case of half filling
($n=1$) for noninteracting electrons [325]. As well as in Ref. [325], at
$n_{a}=n_{b}=0.4$, there appears the nesting of Fermi surface sheets, which
corresponds to the translation by vector $\mathbf{Q}=\\{\pi,\,\pi\\}$ implying
a possibility of arising some orbital superstructure. Owing to the suppression
of correlation effects at this specific doping level, the appearing
superstructure could be the same as that discussed in Ref. [325], but, of
course, this issue demands a special analysis.

For the case of the disordered state under study, dependence of energy on the
total number $n$ of electrons per site is is shown in Fig. 52 ( red solid
line). It is interesting that the shape of this curve is similar to that
obtained in subsubsection 7.5.2 for the simplified model, see Fig. 51. In this
case, the correlations themselves do not give rise to the phase separation.
Nevertheless, the comparison of such orbitally disordered model with the model
that discussed in subsubsection 7.6.2 reveals that within a certain doping
range (at small hole doping), the phase separation could be provide a gain in
energy. The energy of such possible phase-separated state is shown by the
green dot-and-dash line in Fig. 52.

###  7.7 Conclusions

In this section, we have analyzed a possibility of electronic phase separation
in the systems with itinerant charge carriers interacting with an orbitally
ordered array. It was shown that such itinerant charge carriers favor the
formation of an inhomogeneous state with nanosize inhomogeneities having the
orbital structure different from that in the host material. For the 2D lattice
with $e_{g}$ orbitals, we have found the shapes and sizes of such
inhomogeneities. It is shown that the inhomogeneities have he shape strongly
dependent on the ratio of the interorbital interaction constant and
characteristic hopping integral for the charge carriers, $J/t_{0}$: there
exists a threshold value of $J/t_{0}$, corresponding to the transition from
the circular inhomogeneities to needle-like ones.

The orbital model under discussion has certain similarities with that
corresponding to the double exchange. Here, the operators describing the
orbital states are similar to the spin operator. In the previous sections, we
demonstrated that the double exchange model exhibits an instability toward
phase separation into inhomogeneities with different types of magnetic
structure. The inhomogeneous state with circular ferro-OO droplets is, in
fact, similar to that involving magnetic polarons (ferromagnetic droplets
within the antiferromagnetic host material. The latter state is inherent to
the double exchange [111, 28, 135]. However, this similarity is not trivial:
the orbital model has its own specific features owing to the contribution of
non-diagonal hopping and the corresponding anisotropy. Such features give rise
to the behavior specific of orbital systems, such as the kink in the energy of
a homogeneous state and canted-OO needle-like droplets. Note also that the
difference between spins and orbitals leads to many interesting effects, see,
e.g., Refs. [326, 327].

Of course, in actual magnetic materials, we are often deal with the usual
mechanism underlying the formation of magnetic polarons (in simple cases -
ferromagnetic inhomogeneities in AF matrix), which ensures a gain in the
kinetic energy of conduction electrons). However, even in this situation, the
possible existence of orbital degeneracy and orbital ordering present, it is
necessary to consider the orbital structure and its effect on the motion of
charge carriers. For example, the magnetic ordering (parallel spins) could not
impede the electron motion, but if it occurs on the in the presence of
antiferro-OO (such as. alternating $|2x^{2}-y^{2}-z^{2}\rangle$ and
$|2y^{2}-z^{2}-x^{2}\rangle$ orbitals in the basal plane of manganites, which
coexists with the ferromagnetic spin ordering in this plane), such antiferro-
OO itself should affect the motion of itinerant charge carriers. Therefore,
this favor the formation of orbital polarons even at the ferromagnetic spin
background. This implies that in the compounds with JT ions, the
ferromagnetism does not exclude the formation of orbital inhomogeneities. Such
reasoning is in agreement with the Monte Carlo simulations for the two-orbital
model in reported in Ref. [328], which demonstrate the possibility of
nanoscale phase separation even in the case of ferromagnetism. In this
section, we have not considered any magnetic structure and spin ordering of
charge carriers A due account taken for the spin degrees of freedom can give
rise to a rich variety of inhomogeneities with different orbital and spin
structures.

In the first part of this section, we considered the localized electrons
forming an orbital order, whereas conduction electrons or holes were supposed
to be two different groups of electrons. However, we have demonstrated that
that our main results are also valid for a model, where the same electrons
take part both in the hopping and in the formation of orbitally ordered
structure. For example, in magnetic oxides with Jahn–Teller ions, lattice
distortions lift the orbital degeneracy, this leading to the formation of the
ground state with an orbital ordering at $n=1$. In this case, it is natural to
assume the long-range orbital ordering does not disappear at low doping
levels, $x=1-n\ll 1$. Thus, we are dealing with holes at an orbitally ordered
background. Using the mean-field approximation for the analysis of orbitally
ordered arrays, we should introduce the factor $1-x$, thus obtaining the
effective interorbital exchange interaction $J(1-x)^{2}$. In the actual
compounds with Jahn–Teller ions, the orbital Hamiltonians have more
complicated forms than the Heisenberg-like one. Calculations involving the
superexchange mechanism of the orbital ordering [38] show that the obtained
results still remain valid.

In real substances, At high doping level, real materials can exhibit an
orbitally disordered state or an orbital liquid. For manganites, for example,
such state was discussed in Ref. [329]. This problem is rather controversial,
and other possibilities have been also addressed (e.g., complex orbitals
[330]). In the usual spin spin systems, a disordered state is as a rule more
favorable from the viewpoint of of kinetic energy than the AF one (or in our
case, the antiferro-OO state). However, the detailed analysis of spin systems
undertaken in Refs. [112, 111, 28] clearly demonstrates that even in a
disordered (i.e. paramagnetic) state, an additional energy gain can be
achieved due to the formation a ferromagnetic polaron (“ferron”) even at the
paramagnetic background. It is natural to expect similar possibility for the
systems with orbital degrees of freedom, although this such situation demands
a further analysis. To perform it accurately, it is necessary to take into
account unusual specific features of the electron motion within disordered
media. In particular, one should deal with the disorder, which adjusts itself
in the course of electron motion. This is a quite interesting and complicated
separate task.

##  8 Systems with imperfect nesting and phase separation

###  8.1 Introduction

Among the multiband materials with the tendency to electronic phase separation
there exists an important family of the systems with the nesting of the Fermi
surface. The Fermi surface nesting is a very popular and important concept in
condensed matter physics [331]. The existence of two fragments of the Fermi
surface, which can be matched upon translation by a certain reciprocal lattice
vector (nesting vector $\mathbf{Q}_{0}$), entails an instability of a Fermi-
liquid state. A superstructure or additional order parameter related to the
nesting vector is generated due to the instability. The nesting is widely
invoked for the analysis of charge density wave (CDW) states [332, 333], spin
density waves (SDW) states [334, 335], mechanisms of high-$T_{c}$
superconductivity [336, 337, 338], fluctuating charge/orbital modulation in
magnetic oxides [339], chromium and its alloys [340, 341, 342, 100], organic
metals [343, 344, 345, 346], etc. It is important to emphasize that in a real
material, the nesting may be imperfect, i.e. the Fermi surface fragments can
match only approximately. One of the earliest studies of imperfect nesting was
performed by Rice [100] in the context of chromium and its alloys (see also
review articles [347, 348]). Quite recently, it was demonstrated that the
imperfect-nesting mechanism can be responsible for the nanoscale phase
separation in chromium alloys [97], iron-based superconductors [98], organic
metals (Bechgaard salt) [349], and in doped bilayer graphene [99, 350]. In
this context, the studies of spin and charge inhomogeneities related to the
imperfect nesting are currently especially active in the physics of low-
dimensional compounds [343, 351, 352, 353].

Other types of inhomogeneous states (“stripes”, domain walls) were also
discussed in the literature in the framework of similar models [354].
Moreover, it was shown that the possibility of SDW ordering in the systems
with itinerant charge carriers results in very rich and complicated phase
diagrams involving phase-separation regions [355, 207, 206].

Physical mechanisms underlying the nucleation of the inhomogeneous state in
materials with nesting are, in general, common to other systems considered in
this review. The electronic spectrum instability related to the nesting gives
rise to an order parameter and, consequently, to lowering of the system free
energy. The better is nesting, the larger is this energy gain. Thus, it may be
favorable for the system with imperfect nesting to break up into two phases
having different density of itinerant electrons with the better and worse (or
even absent) nesting.

###  8.2 Phase separation in the Rice model for chromium and its alloys

We start with the model proposed by Rice [100] to describe the incommensurate
antiferromagnetism in chromium alloys (see also review article [348]). This
model is rather simple and allows for a detailed examination of its phase
diagram. Desite its simplicity, it can be applied to describe real important
systems. In our consideration, we mostly follow Ref. [97].

The model band structure corresponds to one spherical electron pocket and one
spherical hole pocket with different radii (imperfect nesting). It includes as
well, another band or bands, which do not participate in the magnetic
ordering. All interactions are ignored except the repulsion between electrons
and holes in the Fermi surface pockets giving rise to the ordering, since even
a small coupling in this channel generates a band instability and opening a
gap in the electron spectrum. The Hamiltonian of such model has the form

$\hat{H}=\sum_{\mathbf{k},\sigma,\alpha}\epsilon^{\alpha}(\mathbf{k})n_{\mathbf{k}\sigma}^{\alpha}+\frac{V}{\cal{V}}\sum_{\mathbf{k,k^{\prime},q},\sigma,\sigma^{\prime}}a^{\dagger}_{\mathbf{k}+\mathbf{q}\sigma}a_{\mathbf{k}\sigma}b^{\dagger}_{\mathbf{k}^{\prime}-\mathbf{q}\sigma^{\prime}}b_{\mathbf{k}^{\prime}\sigma^{\prime}}+\textrm{H.c.},$
(288)

where $\alpha$ is equal to either $a$ (electron pocket), $b$ (hole pocket), or
$c$ (nonmagnetic bands), $a^{\dagger}$ ($b^{\dagger}$) are creation operators
for an electron in the $a$ ($b$) pocket, $n$ is the number operator, $V$ is
the Coulomb interaction, and ${\cal V}$ is the volume. The nonmagnetic $c$
band has a finite density of states $N_{r}$ at the Fermi energy. The energy
spectra for the electron and hole pockets measured relative to the Fermi
energy are taken as

$\displaystyle\epsilon^{a}(\mathbf{k})$ $\displaystyle=$ $\displaystyle\hbar
v_{F}(k-k_{Fa})=\hbar v_{F}(k-k_{F})-\mu,$ (289)
$\displaystyle\epsilon^{b}(\mathbf{k}+\mathbf{Q}_{0})$ $\displaystyle=$
$\displaystyle-\hbar v_{F}(k-k_{Fb})=-\hbar v_{F}(k-k_{F})-\mu,$

where $k_{F}=(k_{Fa}+k_{Fb})/2$, $\mu=\hbar v_{F}(k_{Fa}-k_{Fb})/2$ is the
chemical potential, and the wave vector $\mathbf{Q}_{0}$ connects the centers
of the electron and hole pockets in the reciprocal space. We consider the
weak-coupling regime: $VN_{m}\ll 1$, where $N_{m}=k^{2}_{F}/2\pi^{2}\hbar
v_{F}$. We treat Hamiltonian (288) using a mean-field approach, which gives an
accurate result for weakly interacting electrons. In the case of perfect
nesting, which corresponds to $\mu=0$, the radii of the electron and hole
pockets are identical. If we translate the electron pocket by vector
$\mathbf{Q}_{0}$, its Fermi sphere coincides perfectly with the Fermi sphere
of the hole pocket. If we perform the transformations,
$b_{\mathbf{k}}\rightarrow
b^{\dagger}_{\mathbf{k}},b^{\dagger}_{\mathbf{k}}\rightarrow b_{\mathbf{k}}$,
the interaction constant $V$ and the hole pocket dispersion Eq. (289) change
sign. As a result, the Hamiltonian for $a$ and $b$ bands becomes identical to
two copies of the BCS Hamiltonian. This mapping is useful since it allows us
to use the familiar BCS mean-field approach to study the model.

Performing standard BCS-like calculations, it is easy to show that for $\mu=0$
the system is unstable with respect to the ordering with the AFM order
parameter $\Delta_{0}=\frac{V}{{\cal V}}\sum_{\mathbf{k}}\langle
a^{\dagger}_{\mathbf{k}\sigma}b_{\mathbf{k+Q}_{0}-\sigma}\rangle$. In the
weak-coupling limit under study, we have

$\Delta_{0}=\epsilon_{F}\exp{\left(-1/N_{m}V\right)}\ll\epsilon_{F},$ (290)

where $\epsilon_{F}=\hbar v_{F}k_{F}$. The order parameter $\Delta_{0}$
couples electrons with unequal momenta. Consequently, in real space, the order
parameter $\Delta_{0}$ corresponds to the rotation of the magnetization axis
with wave vector $\mathbf{Q}_{0}$. Usually, the $a$ and $b$ pockets are
located at the high-symmetry points of the Brillouin zone and the vector
$\mathbf{Q}_{0}$ is related to the underlying lattice structure. Thus, this
order is usually called commensurate.

When $\mu\neq 0$, the electron and hole Fermi spheres have different radii,
and do not coincide upon translation. However, the difference between the
translated spheres remains small, if $\mu$ is small and the system energy is
optimized by treating the translation vector
$\mathbf{Q}_{1}=\mathbf{Q}_{0}+\mathbf{Q}$ as a variation parameter. The SDW
order parameter has the form

$\Delta=\frac{V}{{\cal V}}\sum_{\mathbf{k}}\langle
a^{\dagger}_{\mathbf{k},\sigma}b_{\mathbf{k+Q}_{1},-\sigma}\rangle.$ (291)

The vector $\mathbf{Q}$ is small: $|\mathbf{Q}|\ll|\Delta|/\hbar
v_{F}\sim|\mathbf{Q}_{0}|$. Hence, the order parameter $\Delta$ describes an
ordering with a slowly rotating AFM magnetization axis. This rotation is
unrelated to the underlying lattice, and such order is referred to as
incommensurate. If the transformation $b_{\mathbf{k}}\rightarrow
b^{\dagger}_{\mathbf{k}},b^{\dagger}_{\mathbf{k}}\rightarrow b_{\mathbf{k}}$
is applied to a system with nonzero $\mu$, the Hamiltonian of our magnetic
system becomes similar to the Fulde–Ferrel–Larkin–Ovchinnikov (FFLO)
Hamiltonian [356, 357] of a superconductor in the Zeeman magnetic field. Note
that we can introduce the order parameters of slightly different structure
than in Eq. (291), however, the general result does not change significantly
(for details see, e.g., Refs. [97, 357, 358]). In particular, the difference
in the free energy and in the chemical potential calculated for different
forms of the order parameter is small, and in actual systems, the balance may
be shifted by factors unaccounted for by the present model (anisotropy,
disorder, etc.) [97].

In the model under discussion, it is convenient to calculate equilibrium
parameters of the system by minimization of the grand thermodynamic potential

$\Omega=-T\ln{\left[\textrm{Tr}\,e^{-(\hat{H}-\mu\hat{N})}\right]},$ (292)

where $\hat{N}$ is the operator of the total number of particles. To evaluate
$\Omega$ in the mean-field approximation, we need the eigenenergies of
Hamiltonian (288) in the mean-field approach. These are

$E_{1,2}(\mathbf{k})=\frac{\epsilon^{a}(\mathbf{k})+\epsilon^{b}(\mathbf{k+Q_{1}})}{2}\pm\sqrt{\Delta^{2}+\left[\frac{\epsilon^{a}(\mathbf{k})-\epsilon^{b}(\mathbf{k+Q_{1}})}{2}\right]^{2}}.$
(293)

Then, the grand potential equals

$\frac{\Omega}{\cal{V}}=\frac{2\Delta^{2}}{V}\\!-\\!2T\\!\sum_{s}\\!\\!\int{\\!\\!\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}\ln{\\!\left(1\\!+\\!e^{-E_{s}/T}\right)}}-2TN_{r}\\!\\!\int{\\!\\!\ln{\\!\left(1\\!+\\!e^{-(\epsilon-\mu)/T}\right)}d\epsilon}.$
(294)

Here, the first and the second terms are the contributions of the bands
responsible for the SDW ordering, while the third term corresponds to the
nonmagnetic bands. To carry out the integration over $\mathbf{k}$, we expand
the band energies in powers of $|\mathbf{Q}|$ and $\delta
k=|\mathbf{k}|-k_{F}$

$\epsilon^{b}(\mathbf{k+Q_{1}})+\epsilon^{a}(\mathbf{k})\approx
2\mu+2Q\eta,\quad\epsilon^{b}(\mathbf{k+Q_{1}})-\epsilon^{a}(\mathbf{k})\approx
2v_{F}\delta k+2Q\eta,$ (295)

where $Q=2v_{F}|\mathbf{Q}|/2$ and $\eta$ is the cosine of the angle between
$\mathbf{k}$ and $\mathbf{Q}$. Performing the integration, one finds the
following expression [97]

$\displaystyle\\!\\!\\!\\!\\!\\!\delta\Omega=\\!\\!\Omega(\Delta,Q,\mu,T)-\Omega(0,Q,\mu,T)$
(296) $\displaystyle\\!\\!\\!\\!\\!\\!=\\!\\!\frac{k_{F}^{3}{\cal
V}}{\pi^{2}\epsilon_{F}}\\!\left\\{\Delta^{2}\\!\left(\ln{\frac{\Delta}{\Delta_{0}}}\\!-\\!\frac{1}{2}\right)\\!+\\!\frac{Q^{2}}{3}\\!+\\!\mu^{2}\\!+\\!\frac{\pi^{2}T^{2}}{3}\\!+\\!T\\!\int_{0}^{\infty}\\!\\!\\!\\!d\xi\int_{-1}^{1}\\!\\!\\!\\!d\eta\ln{\left[f(Q\eta\\!-\\!\mu\\!-\\!\epsilon)f(\mu\\!-\\!Q\eta\\!-\\!\epsilon)\right]}\right\\},$

where $f(\epsilon)=1/(1+\exp{\epsilon/T})$ is the Fermi function and
$\epsilon=\sqrt{\Delta^{2}+\xi^{2}}$.

The equations for the equilibrium values of the gap $\Delta$ and the magnitude
of the structural vector $Q$ are determined by minimizing $\delta\Omega$ [97]

$\displaystyle\ln{\frac{\Delta}{\Delta_{0}}}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{d\xi}{2\epsilon}\int_{-1}^{1}d\eta\left[f(\epsilon+\mu-Q\eta)+f(\epsilon-\mu+Q\eta)\right],$
$\displaystyle\frac{2Q}{3}$ $\displaystyle=$
$\displaystyle-\int_{0}^{\infty}d\xi\int_{-1}^{1}\eta
d\eta\left[f(\epsilon+\mu-Q\eta)+f(\epsilon-\mu+Q\eta)\right].$ (297)

The total number of electrons per unit volume $n(\mu)$ is the sum of the
magnetic $n_{m}(\mu)$ and nonmagnetic electrons $n_{r}(\mu)$ . The doping $x$
is defined as the difference $x=n(\mu)-n(0)$. Since $T,\mu\ll\epsilon_{F}$, we
can write $n_{r}(\mu)-n_{r}(0)=N_{r}\mu$. For magnetic electrons, we have
$n_{m}(\mu)=(2/{\cal V})\sum_{\mathbf{k},s}{f(E_{s}(\mathbf{k}))}$. After
straightforward calculations, we derive

$\frac{x}{x_{0}}=\frac{n\mu}{\Delta_{0}}+\int_{0}^{\infty}\frac{d\xi}{\Delta_{0}}\int_{-1}^{1}d\eta\left[f(\epsilon-\mu+Q\eta)-f(\epsilon+\mu-Q\eta)\right],$
(298)

where $x_{0}=4\Delta_{0}N_{m}$ and $n=N_{r}/2N_{m}$. Equations (296), (8.2),
and (298) allow constructing a phase diagram of the system [97].

The phase diagram of the model can be constructed by numerical solution of
Eqs. (296), (8.2), and (298). First, one can calculate the Néel temperature
$T_{N}(x)$ for homogeneous state of the system and the transition temperature
$T_{Q}$ between the commensurate and incommensurate AFM phases, which
corresponds to the highest doping, at which $Q=0$. Then, we can determine the
region of the phase-separated state and the phases in this state analyzing the
dependence of the grand potential $\Omega$ or chemical potential $\mu$ versus
doping $x$ following the standard procedure described in the previous sections
(see also for details Ref. [97]).

As an illustration, a typical dependence of the chemical potential $\mu$ on
doping $x$ for the considered model is shown in Fig. 56. It is seen from this
figure that the function $\mu(x)$ has a descending part, which means a
negative compressibility and, hence, instability of the homogeneous AFM phase.
In the phase-separated state, the system segregates into two phases with
different doping values $x_{1}$ and $x_{2}$. The values $x_{1}$ and $x_{2}$
can be found using the Maxwell construction (see, e.g., Ref. [359]). Figure 56
illustrates the latter concept: the horizontal dashed line is drawn in such a
manner that the areas of the shaded regions, $S_{1}$ and $S_{2}$, are equal to
each other.

Figure 56:  (Color online) Chemical potential $\mu$ versus doping $x$ for the
homogeneous phase, $T/\Delta_{0}=0.1$ and $n=0$ [solid (red) line]. The
horizontal dashed (black) line shows the Maxwell construction, with shaded
areas $S_{1}=S_{2}$ [97].

The calculated by the latter method boundary $T_{PS}(x)$ between the
homogeneous and phase-separated states is shown by the dashed (red) curves in
the $(x,T)$ phase diagrams drawn in Fig. 57, for different values of $n$. The
phase with lower doping, $x_{1}$, is the commensurate AFM ($Q=0$), whereas the
phase with higher doping, $x_{2}$, is the incommensurate AFM ($Q\neq 0$).
Thus, the phase separation occurs due to the competition between two AFM
states with different structures.

Figure 57: (Color online) Phase diagram of the Rice model for $n=0$ (no
electrons in nonmagnetic pockets) (a), $n=1$ (b), and $n=3$ (c). The solid
(black) curves represent the Néel temperature $T_{N}(x)$, which separates the
paramagnetic ($\Delta=0$) and AFM ($\Delta\neq 0$) phases. The dotted (blue)
curves denote $T_{Q}(x)$, i.e. the boundary between the commensurate ($Q=0$)
and incommensurate ($Q\neq 0$) homogeneous AFM phases. The dashed (red) curve,
$T_{PS}(x)$, is the boundary between the uniform and phase-separated (shaded
areas) phases [97].

Above, we have mentioned the mapping between the Rice model and the FFLO
superconductor. However, the phase separation is absent from the phase diagram
of the FFLO state in a superconductor. This has a very simple explanation: the
Zeeman field in the FFLO superconductor corresponds to doping in the Rice
model. The field, being an intensive thermodynamic quantity, does not allow
for the phase separation. However, for cold atoms in an optical trap (where an
analogue of the FFLO is predicted), it may be possible to control not the
field, but the polarization, which is an extensive quantity. In this case, the
phase separation occurs [360].

The FFLO state is very fragile with respect to the impurity scattering [361].
Thus, the incommensurate AFM, which is the mathematical analog of the FFLO
phase, is expected to be susceptible to microscopic imperfections. Therefore,
the disorder-induced modifications to the phase diagram is an open question.

The study of characteristic scales and geometry of the phase-separated state
could be performed by the method similar to that described in the previous
sections of the present review. It requires, as usual, to take into
consideration an interplay between the long-range Coulomb interaction and the
energy of the boundary between different phases. However, such a study is
absent in literature up to date.

###  8.3 The effect of magnetic field on the phase separation in systems with
imperfect nesting

Thus, the electron–electron coupling gives rise to the nucleation of the
magnetic ordering and the phase separation in the system with nesting.
Naturally, the magnetically ordered system is especially sensitive to the
applied magnetic field. This problem was investigated in Ref. [362] in the
framework of the considered above Rice model. We analyze the effects of the
applied magnetic field following the latter paper.

Considering electrons or holes in the applied magnetic field $\mathbf{B}$, we
should, first, replace the momentum operator $\mathbf{\hat{p}}$ by the gauge-
invariant one $\mathbf{\hat{p}}+(e/c)\mathbf{A}$, where $\mathbf{A}$ is the
vector potential of magnetic field. Second, we have to add the Zeeman term
$g_{\alpha}\sigma\hbar\omega_{\alpha}$ in the single electron Hamiltonian,
where $g_{\alpha}$ is the Landé $g$-factor and
$\omega_{\alpha}=eB/cm_{\alpha}$ is the cyclotron frequency of a charge
carrier with mass $m_{\alpha}$. The energy spectrum of the Rice model is
characterized by two single-particle energy scales. The first scale is the
Fermi energy $\epsilon_{F}$, and the second one is $\hbar\omega_{\alpha}$,
which is the distance between the Landau levels. The energy scale associated
with the interactions will be characterized by the value of a spectral gap
$\Delta_{0}$ and in the case of the weak electron–hole coupling
$\Delta_{0}\ll\epsilon_{F}$. The Landau quantization is of importance in the
range of high magnetic field, $\hbar\omega_{\alpha}>\Delta_{0}$, whereas at
low fields, $\hbar\omega_{\alpha}<\Delta_{0}$, it can be neglected. Here, we
limit ourselves by the regime of low magnetic fields, neglecting any
corrections associated with the small ratio
$\hbar\omega_{\alpha}/\epsilon_{F}\ll 1$. In the range of high magnetic
fields, there exists a set of oscillation phenomena, which highly complicates
the phase diagram even of this simple model.

Adding the Zeeman terms, we write instead of Eq. (289) the single-electron
spectrum in the form

$\displaystyle\epsilon^{a}_{\sigma}(\mathbf{k})$ $\displaystyle=$
$\displaystyle\hbar v_{F}(k-k_{F})-\mu+g_{a}\sigma\hbar\omega_{a},$ (299)
$\displaystyle\epsilon^{b}_{\sigma}(\mathbf{k}+\mathbf{Q}_{0})$
$\displaystyle=$ $\displaystyle-\hbar
v_{F}(k-k_{F})-\mu+g_{b}\sigma\hbar\omega_{b}.$

The applied magnetic field lifts the degeneracy with respect to the spin
projection, and now we have two electron and two hole bands, see Fig. 58. In
such case we have to introduce a two-component order parameter corresponding
to the nesting vectors shown by the arrows in Fig. 58. Neglecting for
simplicity the incommensurate SDW ordering, these components can be written as
[362]

$\Delta_{\uparrow}=\frac{V}{\cal{V}}\sum_{\mathbf{k}}\langle
a^{\dagger}_{\mathbf{k}\uparrow}b_{\mathbf{k}\downarrow}\rangle,\quad\Delta_{\downarrow}=\frac{V}{\cal{V}}\sum_{\mathbf{k}}\langle
a^{\dagger}_{\mathbf{k}\downarrow}b_{\mathbf{k}\uparrow}\rangle.$ (300)

The case when the electron and hole bands are perfectly symmetric is the
simplest. In such a situation, however, the effect of Zeeman field on the
electron spectrum of the model is absent. Thus, some electron–hole asymmetry
is necessary to obtain a non-trivial result. Following Ref. [362], we put
$m_{a}=m_{b}=m$, hence, $\omega_{a}=\omega_{b}=\omega_{H}$, but $g_{a}\neq
g_{b}$. Following the procedure described in the previous subsection, now we
can find eigenenergies and calculate the grand potential $\Omega$.
Minimization of $\Omega$ gives the expression for components of the order
parameter [362]

$\ln{\frac{\Delta_{0}}{\Delta_{\sigma}}}=\int_{0}^{\infty}\\!\\!\\!d\xi\,\frac{f(\sqrt{\Delta^{2}_{\sigma}+\xi^{2}}+\mu_{\sigma})+f(\sqrt{\Delta^{2}_{\sigma}+\xi^{2}}-\mu_{\sigma})}{\sqrt{\Delta^{2}_{\sigma}+\xi^{2}}},$
(301)

where $\mu_{\sigma}=\mu-\sigma(g_{a}-g_{b})\hbar\omega_{H}/2$. Neglecting the
existence of the nonmagnetic part of the Fermi surface, one can write the
expression for doping analogous to Eq. (298)

$\frac{x}{x_{0}}=\sum_{\sigma}{\int_{0}^{\infty}\frac{d\xi}{4\Delta_{0}}\left[f(\sqrt{\Delta^{2}_{\sigma}+\xi^{2}}-\mu_{\sigma})-f(\sqrt{\Delta^{2}_{\sigma}+\xi^{2}}+\mu_{\sigma})\right]}.$
(302)

Figure 58:  (Color online) Band structure of the model in an applied magnetic
field. The applied magnetic field lifts the degeneracy of electron-like ($a$)
and hole-like ($b$) bands with respect to the electron spin. The red arrows
indicate the interband coupling giving rise to the order parameters [362].

In Ref. [362], the phase diagram of the model in the plane (magnetic field,
doping) was constructed in the limit $T=0$. The obtained results are
illustrated in Fig. 59, where the dimensionless magnetic field is defined as
$b=|g_{a}-g_{b}|\hbar\omega_{H}/2\Delta_{0}$. The homogeneous PM state exists
in the system at high doping ($x/4x_{0}>2b+\sqrt{2}$) or at high magnetic
field ($b>1/\sqrt{2}+x/8x_{0}$) (the value $x_{0}$ in Ref. [362] corresponds
to $4x_{0}$ in the present notation). In the range of the low doping and low
magnetic field, the ground state of the system is the phase-separated one,
which we refer to as PS1. The phase PS1 is a mixture of the two AFM phases.
The first phase corresponds to SDW with zero doping, $x=0$, and
$\Delta_{\uparrow}=\Delta_{\downarrow}=\Delta_{0}$. The second phase with
nonzero doping, $\Delta_{\uparrow}=\Delta_{0}$, and $\Delta_{\downarrow}=0$
will be further referred to as AF1. The phase PS1 is the ground state of the
model if $|b|<1/\sqrt{2}$ and $|x/4x_{0}|<1/\sqrt{2}$. The uniform SDW phase
AF1 is stable within the intermediate magnetic fields range,
$x/8x_{0}-\/2\sqrt{2}<b<x/8x_{0}+\/2\sqrt{2}$. The inhomogeneous state PS2 is
stable in the regions between homogeneous PM and AF1 phases. The PS2 state
consists of the mixture of the SDW phase AF1 and PM phase. It is interesting
to note that the phase AF1 is some non-common type of the half-metal, since in
this phase only electrons with a single spin projection exist at the Fermi
surface [363].

Thus, we see that the magnetic field lifts the spin-projection degeneracy of
the SDW order parameter in the system with imperfect nesting. As a result the
number of possible homogeneous states increases, which, in turn, drastically
affects the picture of electronic phase separation.

Figure 59:  (Color online) Phase diagram of the model at zero temperature in
the ($b,x/x_{0}$) plane [362]. The value $x_{0}$ in Ref. [362] corresponds to
$4x_{0}$ in the present notation. In the diagram, we see the homogeneous SDW
phase AF1 (the definition is shown in the figure), two regions with the
uniform PM phase, and the phase-separated states PS1 (mixture of two SDW
phases, see the text) and PS2 (mixture of AF1 and the PM phases).

###  8.4 Phase separation in the organic metals under pressure

In the above consideration we demonstrated that the system with nesting is
very sensitive to the position of the Fermi level due to the doping and/or to
the change of the electron spectrum due to application of the magnetic field.
This conclusion can be formulated in a more general form: any impact that
affects the value of the chemical potential or electron spectrum in a system
with nesting can give rise to the crucial changes in its thermodynamic state.
Here we discuss the effect of the pressure. This case is of a special interest
due to the experiments with an archetypical Bechgaard salt (TMTSF)2PF6 under
pressure [343, 344, 345, 346]. This system is commonly considered as a system
with the nesting of the Fermi surface. In these experiments the coexistence of
the SDW and PM metallic (superconducting at low temperatures) phases was
clearly observed. Moreover, in Ref. [343] the SDW and PM metallic phases
manifest itself in the form of phase separation with macroscopic domains of
the high-pressure PM metallic phase embedded into the insulating SDW host,
which are aligned along certain crystallographic axes. A natural explanation
of this effect is the same as in the case of doping or magnetic field: the
application of the pressure shifts the electron system from the position of a
perfect nesting.

A qualitative explanation of the experiment [343] can be obtained by an
incorporation of the pressure effect in Rice model [349]. In the experiments
[343], an evolution with pressure of the geometry of the SDW-PM phase was also
observed. To trace this evolution theoretically the electrostatic energy of
the inhomogeneous state and the energy of interfaces between SDW and PM phases
should be taken into account as it was described in section 5.5.

In general, the pressure $P$ may affect all parameters of the system, e.g.,
electron coupling $V$, electron and hole effective masses, band positions
$\epsilon^{\alpha}_{\text{min}}$, $\epsilon^{\alpha}_{\text{max}}$, etc.
However, if the pressure is not very high, the main effect on the electron
properties of the material with nesting occurs due to the shift of the Fermi
energy closer or farther from the ideal nesting position. Following Ref.
[349], we assume that the change in the Fermi energy arises due to the shift
of the energy bands with respect to each other. For simplicity, in the cited
work only the shift the reservoir band is taken into consideration while the
shift of the magnetic bands is ignored. Below we discuss only the case of the
commensurate SDW, that is, $\mathbf{Q}=0$. It occurs that this simple minimal
model is sufficient to reproduce main features of the experimental
observations reported in Ref. [343].

If the position of the bottom of the reservoir band is shifted on
$\delta\epsilon$ due to pressure, then, the corresponding variation of the the
chemical potential is $\mu(P)-\mu(0)$. Naturally, the total number of
electrons per unit cell (and, consequently, doping $x$) remains unchanged. The
variation in the number electrons per unit cell in the nonmagnetic band is

$\delta n_{r}=2n_{r}[\mu(P)-\mu(0)-\delta\epsilon]\,.$ (303)

while the variation in the number of magnetic electrons is

$\displaystyle\delta n_{m}$ $\displaystyle=$ $\displaystyle
4n_{m}\int\limits_{0}^{\infty}\\!\\!{d\xi\,\\{f[\eta-\mu(P)]-f[\eta+\mu(P)]-f[\eta\\!-\\!\mu(0)]-f[\eta\\!+\\!\mu(0)]\\}}.$
(304)

We assume that the system is perfectly nested when $P=P_{n}$, that is,
$\mu(P_{n})=0$ and $x=0$. In the range of not too high pressures the value
$\delta\epsilon$ depends linearly on the pressure

$\delta\epsilon(P)=\frac{\partial\epsilon^{c}_{\text{min}}}{\partial
P}\left(P-P_{n}\right)\,.$ (305)

We substitute Eqs. (303) – (305) into Eq. (298) with $Q=0$ and make the
replacement $\mu+\delta\mu\to\mu$. Then, taking into account the condition
$x(P)=const=0$, one derives [349]

$r\bar{P}=\frac{r\mu}{\Delta_{0}}+\frac{1}{\Delta_{0}}\\!\int\limits_{0}^{\infty}\\!\\!{d\xi\,[f(\eta-\mu)-f(\eta+\mu)]}\,,$
(306)

where the dimensionless variable $\bar{P}$ is

$\bar{P}=\frac{\partial\epsilon^{c}_{\text{min}}}{\partial
P}\frac{P-P_{n}}{\Delta_{0}}\,.$ (307)

The sign of the derivative $\partial\epsilon^{c}_{\text{min}}/\partial P$
depends on the particular properties of the system. Below we assume that this
value is positive. In the case of zero doping, $x=0$, Eqs. (8.2) at $Q=0$ and
(306) form a closed system of equations determining the functions
$\Delta(\bar{P},T)$ and $\mu(\bar{P},T)$. If $x\neq 0$, one has to substitute
$r\bar{P}\to r\bar{P}+x/x_{0}$ in the left-hand side of Eq. (306), taking into
account that under doping, pressure $P=P_{n}$ corresponds to a perfect
nesting.

Comparing Eq. (306) with Eq. (298), we see that the pressure acts as some
effective doping. The only difference is that this ‘effective doping’ is
proportional to the density of electron states in the reservoir bands, $r$ and
when the reservoir is absent, the pressure changes nothing. This is an
artefact of the approximation used, where possible changes in the shapes of
the Fermi surfaces of the magnetic bands were neglected. Under such assumption
the only possible way to destroy the ideal nesting is to change the difference
in the Fermi momenta $k_{Fa}-k_{Fb}$ of the magnetic bands, which is
prohibited due to the charge conservation law. When the reservoir is present,
the imbalance in the magnetic electrons appears due to the unequal $k_{Fa}$
and $k_{Fb}$ can be compensated by the reservoir electrons.

The equation for the gap, Eq. (8.2), at $Q=0$ and (306) can be solved
numerically to obtain the functions $\Delta(p,T)$ and $\mu(\bar{P},T)$. Then,
assuming that the ground state of the system is homogeneous, one can find the
grand potential $\Omega$ and the free energy $F=\Omega+\mu n$, and construct
the phase diagram of the model in the $p$–$T$ plane. Such calculations were
performed in Ref. [349] and the phase diagram of the model in the $(P,T)$
plane was obtained. The calculated dependence of the gap $\Delta$ on the
pressure occurs similar to the dependence of the gap versus doping: the gap
decreases with $P-P_{n}$ and vanishes zero at some value of the applied
pressure. However, the true ground state of the model is inhomogeneous similar
to the cases considered in sections 8.2 and 8.3.

In the case of the finite doping $x$ Eq. (306) becomes [349]

$\frac{x_{\text{eff}}(P)}{x_{0}}=\frac{r\mu}{\Delta_{0}}+\frac{1}{\Delta_{0}}\\!\int\limits_{0}^{\infty}\\!\\!{d\xi\,[f(\eta-\mu)-f(\eta+\mu)]}\,,$
(308)

where

$x_{\text{eff}}(P)=x+rx_{0}\bar{P}\,.$ (309)

First, we have to calculate the chemical potential $\mu$ as a function of
$x_{\text{eff}}$. The obtained dependence $\mu(x_{\text{eff}})$ is similar to
$\mu(x)$ shown in Fig. 56. Thus, the value $\mu(x_{\text{eff}})$ behaves
nonmonotonically indicating an instability of the homogeneous state toward the
phase separation. The separated phases are the SDW phase with effective doping
$x^{(1)}_{\text{eff}}$ and the PM phase with effective doping
$x^{(2)}_{\text{eff}}$. As above, the values $x^{(1,2)}_{\text{eff}}$ can be
obtained using the Maxwell construction, see Fig. 56. These values depend on
the model parameters and temperature, but do not depend on pressure, while
“bare” doping in the in the phase-separated phases, $x_{1,2}$, depend on $P$.
From Eq. (309) it follows that $x_{1,2}=x^{(1,2)}_{\text{eff}}-rx_{0}p$. If
the volume fraction of the PM phase is $p$, then, from the charge conservation
law, we have

$(1-p)x_{1}+px_{2}=x\,.$ (310)

In the doping range, $x_{1}<x<x_{2}$, the ground state of the model is a
mixture of SDW phase with the electron density
$x_{1}=x^{(1)}_{\text{eff}}-rx_{0}\bar{P}$ and PM phase with the electron
density $x_{2}=x^{(2)}_{\text{eff}}-rx_{0}\bar{P}$. This doping range shifts
toward the smaller densities when the pressure increases. Suppose that
$x<x^{(1)}_{\text{eff}}$, thus, at $P=P_{n}$ (when $\bar{P}=0$), the ground
state is a homogeneous SDW phase. Increasing $\bar{P}$ above the value
$\bar{P}_{1}=(x^{(1)}_{\text{eff}}-x)/(rx_{0})$ makes the system
inhomogeneous. The fraction of the PM phase, $p$, gradually increases with
$\bar{P}$ and at $\bar{P}=\bar{P}_{2}=(x^{(2)}_{\text{eff}}-x)/(rx_{0})$, the
system becomes a homogeneous PM. As an illustration, the phase diagram of the
model in the $(P,T)$ plane is shown in Fig. 60 for $r\ll 1$ (a) and for $r=2$
(b). The region of the possible inhomogeneous phase occupies a significant
part of the phase diagram.

Figure 60: (Color online) Phase diagram of the model in the
$(x_{\text{eff}},T)$ plane in the cases $r\ll 1$ (a) and $r=2$ (b). Solid
(red) line with up triangles shows the boundary of the uniform stable SDW
phase. The shaded area corresponds to the region with tendency to the phase
separation (PS). (Fig. 3 from Ref. [349])

As it was discussed in section 5.5, the change in the phases volume fractions
is accomplished by the evolution of the inhomogeneous phase geometry, see Fig.
35. This evolution occurs due to change of the ratio of the long-range Coulomb
interaction $E_{c}$ arising due to the violation of the charge neutrality and
the energy $E_{s}$ of the interface between SDW and PM phases. [217, 191]

Figure 61: (Color online) $(P,T)$ phase diagram in the case $x=0$ and $r=2$.
Solid (red) lines show the boundary of the uniform stable phases. Blue dashed
lines show the boundaries of the inhomogeneous phases with different shapes of
the inhomogeneities: PM droplets in the SDW host exist in region 1; PM pillars
in the SDW host exist in region 2; alternating PM and SDW slabs exist in
region 3; SDW pillars in the PM host exist in region 4; and SDW droplets in
the PM host exist in region 5. (Fig. 4 from Ref. [349])

The evolution of the inhomogeneous state geometry with pressure was analyzed
in Ref. [349] The content of the PM phase $p$ can be expressed as (for $x=0$)

$p=\frac{rx_{0}\bar{P}-x^{(1)}_{\text{eff}}}{x^{(2)}_{\text{eff}}-x^{(1)}_{\text{eff}}}\,,$
(311)

where the effective dopings $x^{(1,2)}_{\text{eff}}$ are calculated according
to the Maxwell rule. Then, following the procedure described in section 5.5
the phase diagram of the system can be calculated. The results of such
calculations are illustrated in Fig. 61. The volume fraction of the PM phase,
$p$, monotonically increases with pressure. As a result, the sequence of the
geometries of inclusions changing with the growth of pressure is the
following: uniform SDW phase at low pressure is changed by the phase separated
state with PM islands in the SDW host (area 1 in Fig. 61), then, by PM rods or
pillars (region 2), PM slabs (region 3) and, finally, there arise inclusions
of the SDW phase in the PM host (regions 4 and 5). The latter inhomogeneous
states could be seen as a uniform PM metallic state in the transport
measurements due to metal-insulator percolation transition. [349]

Figure 62: (Color online) (a) Experimental observation of the phase separation
in the Bechgaard salt [343] (Fig. 1 from Ref. [343]). The uniform SDW phase is
observed if pressure $P<P_{0}$. The inhomogeneous mixture of the SDW and PM
metallic phases exists if $P_{0}<P<P_{c}$. The PS state includes metallic
droplets in the SDW host if $P_{0}<P<P_{1}$, the metallic pillars in the SDW
matrix if $P_{1}<P<P_{2}$, and the metallic slabs in the SDW host if
$P_{2}<P<P_{c}$. When $P>P_{c}$ the sample is in the metallic state. The metal
is superconducting if $T<T_{c}$, the dependence $T_{c}(P)$ is shown by (blue)
solid line with circles. (b) A part of the phase diagram shown in Fig. 61
redrawn in different coordinates, the notation corresponds to panel (a). (Fig.
5 from Ref. [349]).

In the above consideration, the pressure changes only the values of $k_{Fa}$
and $k_{Fb}$, but does not change the spherical geometry of the Fermi surface
pockets. Naturally, it is a simplification. Below we discuss the phase
separation observed in the Bechgaard salt [343] following Ref. [349]. The
Bechgaard salts are highly anisotropic and one may expect that its Fermi
surface is also anisotropic, and its nested parts are probably not spherical.
The pressure can affect not only the chemical potential, but also the geometry
of the nested pockets. Thus, we need to include into consideration some
additional geometrical factor, which governs the nesting. In this case, the
de-nesting under pressure can occur even in the absence of the nonmagnetic
reservoir. At the same time, if the electron and hole pockets are nested only
partially, one can suggest that the SDW order parameter arises only in the
nested parts of the Fermi surfaces, while other parts remain ungapped and play
a role of the reservoir. Thus, one can expect that the model considered in
this section could be relevant, at least, qualitatively, even in the case of
nonspherical Fermi surfaces. (A more detailed analysis of the effects of the
Fermi surface shape is given in the next section, where the phase separation
in iron pnictide superconductors is discussed.)

The de-nesting of any nature gives rise to the decrease of the SDW order
parameter and to the phase separation. With the increase of the pressure, the
uniform SDW order is destroyed and the phase separation arises at some
$P=P_{0}$. The volume concentration of the metallic phase increases with
pressure and the described above evolution of the inhomogeneous phase geometry
is observed, see Fig. 61. The metallic phase may be superconducting at low
temperature. In Fig. 62(a), the experimental results on the phase separation
in the Bechgaard salt under pressure taken from Ref. [343] are shown. In Fig.
62(b), for comparison, a part of the phase diagram shown in Fig. 61 redrawn in
the proper coordinates is presented. The inhomogeneous phases at high pressure
where the sample is a metal from the experimental point of view are omitted.
Qualitatively, the pictures illustrating the experimental and theoretical
results are quite similar.

###  8.5 Phase separation in iron pnictide superconductors

Superconducting iron-based pnictides form another system with the imperfect
nesting, which we consider in this section. The phase diagram of iron
pnictides is very rich and contains the regions of superconductivity and spin-
density-wave order, both commensurate and incommensurate. Moreover, these
materials often exhibit spin and charge inhomogeneity with the characteristic
features of the systems with electronic phase separation [364, 365, 366, 367,
368, 369, 98]. Superconducting iron chalcogenides also exhibit a phase
separation [370, 371]. However, the physics related to the Fermi surface
nesting may not be directly applicable to iron chalcogenides since some
chalcogenides have no hole pockets.

Figure 63:  (Color online) Schematic illustration of the Fermi surface of iron
pnictides in the unfolded (a) and folded (b) Brillouin zone (BZ). We have
three hole pockets located near the center of the folded BZ [shown by (blue)
solid curves]. In the unfolded BZ, one of the hole pockets moves to the corner
of the BZ. For the unfolded BZ (a), the electron pockets are elliptic [(red)
dashed curves] and located near the ($0,\pi/a$) and ($\pi/a,0$) points, while
in the folded BZ, they are represented by overlapping ellipses located at the
corners. Arrows show the possible nesting vectors between hole and electron
pockets, which give rise to the SDW order. (c) Illustration of the possible
nesting at different doping levels. (Fig. 1 from Ref. [98]).

Unlike superconducting cuprates, which are believed to be in the strong
electron–electron interaction regime, iron pnictides may be described by a
weak-interaction model (see, e.g., Sec.IIIA of Ref. [372]).

Currently, the Fermi surface of iron pnictides is typically described using
two related representations briefly outlined below. In Fig. 63(a), we plot the
Fermi surface within the so-called unfolded Brillouin zone [372, 373, 374],
which corresponds to the square lattice of iron atoms, with one Fe atom per
unit cell and lattice constant $a$. In this representation, two quasi-two-
dimensional nearly circular hole pockets are centered at the $\Gamma(0,0)$
point, one more circular hole pocket is located near the
$\Gamma^{\prime}(\pi/a,\pi/a)$ point, and two elliptically-shaped electron
pockets are centered at the $M(0,\pi/a)$ and $M(\pi/a,0)$ points [see Fig.
63(a)]. At the same time, the actual unit cell of pnictides contains two Fe
atoms, since in the crystal lattice, the pnictogen atoms are located at
nonequivalent positions. The (folded) Brillouin zone corresponding to this
unit cell is obtained by folding the Brillouin zone shown in Fig. 63(a) by
dashed lines and consequent rotation by $\pi/4$. The Fermi surface in the
folded Brillouin zone is shown in Fig. 63(b). In this figure, all three hole
pockets are located near the $\Gamma(0,0)$ point, whereas the electron pockets
represented by overlapping ellipses are located near the
$M(\pm\pi/\bar{a},\pi/\bar{a})$ points, with $\bar{a}=a/\sqrt{2}$. Formulating
this model, we will make several simplifications. First, we neglect the
effects associated with nonequivalent positions of the Fe atoms. Consequently,
the use of unfolded Brillouin zone is sufficient. Second, we will neglect the
3D structure of the material and study only the 2D model.

Following Ref. [98], the Hamiltonian of the model can be written as
$\hat{H}=\hat{H}_{0}+\hat{H}_{\textrm{int}}$, where the kinetic energy term
has the form

$\hat{H}_{0}=\sum_{\mathbf{k}\lambda\sigma}{\epsilon_{\lambda\mathbf{k}}^{h}a^{\dagger}_{\mathbf{k}\lambda\sigma}a_{\mathbf{k}\lambda\sigma}}+\sum_{\mathbf{k}s\sigma}{\epsilon_{s\mathbf{k}}^{e}b^{\dagger}_{\mathbf{k}s\sigma}b_{\mathbf{k}s\sigma}}.$
(312)

Here $a^{\dagger}_{\mathbf{k}\lambda\sigma},a_{\mathbf{k}\lambda\sigma}$
($b^{\dagger}_{\mathbf{k}s\sigma},b_{\mathbf{k}s\sigma}$) are the creation and
annihilation operators for charge carriers in the holelike (electronlike)
bands $\lambda=1,2,3$ ($s=1,2$) with spectra
$\epsilon_{\lambda\mathbf{k}}^{h}$ ($\epsilon_{s\mathbf{k}}^{e}$). Assuming
that the bands have quadratic dispersion laws near the Fermi level, for
circular holelike bands we have

$\epsilon_{1\mathbf{k}}^{h}=-\frac{\hbar
v_{F}^{h}(k^{2}-k_{F}^{2})}{2k_{F}}-\mu,\,\,\epsilon_{2\mathbf{k}}^{h}=-\frac{\hbar
v_{2F}^{h}(k^{2}-k_{F}^{2})}{2k_{F}}-\Delta\epsilon_{2}-\mu,\,\,\epsilon_{3\mathbf{k+\bar{Q}}}^{h}=-\frac{\hbar
v_{3F}^{h}(k^{2}-k_{F}^{2})}{2k_{F}}-\Delta\epsilon_{3}-\mu.$ (313)

Here, $v_{F}$ and $v_{2,3F}$ are the Fermi velocities for the hole bands,
$\mathbf{\bar{Q}}=(\pi/a,\pi/a)$, and the energy shifts $\Delta\epsilon_{2,3}$
determine the difference in radii of the hole pockets. The electronlike
components of the Fermi surface are elliptic and can be presented as follows
[98]

$\epsilon_{1\mathbf{k+Q_{0}}}^{e}=\epsilon_{1\mathbf{k+Q_{0}^{\prime}}}^{e}=\frac{\hbar
v_{F}^{e}(k^{2}-k_{F}^{2})}{2k_{F}}+\frac{\alpha\hbar
v_{F}^{e}}{k_{F}}(k_{x}^{2}-k_{y}^{2})-\mu,$ (314)

where the centers of the elliptic bands are $\mathbf{Q}_{0}=(\pi/a,0)$ and
$\mathbf{Q_{0}^{\prime}}=(0,\pi p/a)$, and $v_{F}^{e}$ is the Fermi velocity
for the electron bands averaged over the Fermi surface. The parameter $\alpha$
defines the ellipticity of the electron pockets. The major axes of the
ellipses are directed toward the $\Gamma$ point when $\alpha>0$ and
perpendicular to this direction when $\alpha<0$.

Due to the multisheet structure of the Fermi surface, the interaction
Hamiltonian $\hat{H}_{\textrm{int}}$, in general, includes a number of
contributions. However, since we are interested in the SDW order and the phase
separation, one can ignore the electron–electron and hole–hole interactions as
in the previous analysis of the Rice model. It is known [364] that the value
of local magnetic moment in pnictides oscillates along one of the crystal
axes. To reproduce this stripy magnetic structure, it is sufficient to couple
one electron and one hole band [375]. For the definiteness sake, one can
assume that these are the hole band $\epsilon^{h}_{1}$ and the electron band
$\epsilon^{e}_{1}$. Thus, we can write

$\hat{H}_{\textrm{int}}=\frac{V}{{\cal
V}}\sum_{\mathbf{kk^{\prime}q}\sigma\sigma^{\prime}}{a^{\dagger}_{\mathbf{k+q}1\sigma}a_{\mathbf{k}1\sigma}b^{\dagger}_{\mathbf{k^{\prime}-q}1\sigma^{\prime}}b_{\mathbf{k^{\prime}}1\sigma^{\prime}}}.$
(315)

The described model for pnictides is a generalization of the Rice model. As in
the latter case, two bands participate in the magnetic transition, and there
is a “reservoir” (nonmagnetic bands). However, the perfect nesting is
impossible for pnictides since the forms of the hole and electron pockets are
different for any doping level.

As we pointed out above, the Coulomb interaction in iron pnictides is weak and
one can assume that $V/\epsilon_{F}\ll 1$, where the Fermi energy
$\epsilon_{F}=(v_{F}^{e}+v_{F}^{h})k_{F}/2$. Then, Hamiltonian (315) can be
considered in the mean-field approximation [98] and, as in the previous
section, we will search the SDW order parameter in the form (291). In general,
different parts of each pocket of the Fermi surface have different orbital
composition. Such feature can be accounted by introducing a momentum-dependent
coupling $V$ and would lead to a pronounced wave-vector dependence of the SDW
order parameter. However, unless the variations of $V$ are extremely strong,
they do not affect the phase separation qualitatively.

The magnetization corresponding to the chosen SDW order lies in the $xy$
plane. For the commensurate SDW case, the magnetization direction remains
constant when one moves along the direction normal to $\mathbf{Q}_{0}$. When
one moves parallel to $\mathbf{Q}_{0}$, the magnetization reverses its
direction from one iron atom to the next iron atom. For incommensurate SDW,
this “stripy” pattern slowly rotates in the $xy$ plane [98].

Performing mean-field transformations in the interaction Hamiltonian (315), we
obtain the energy spectrum of the model similar to Eq. (293) and grand
potential in the form (294), where energies $\epsilon^{a,b}$ should be
replaced by $\epsilon^{e,h}$. As usual, the minimization of the grand
potential gives the equations for the SDW order parameter $\Delta$ and nesting
vector $\mathbf{Q}$. The expression for $\Delta$ was derived in Ref. [98] at
zero temperature

$\ln{\frac{\Delta_{0}}{\Delta}}=\int{\frac{\epsilon_{F}d^{2}\mathbf{k}}{2\pi
k_{F}^{2}}\frac{\Theta[E_{1}(\mathbf{k})]+\Theta[-E_{2}(\mathbf{k})]}{\sqrt{\Delta^{2}+\frac{1}{4}\left(\epsilon^{e}_{1\mathbf{k}+\mathbf{Q}_{0}}-\epsilon^{h}_{1\mathbf{k}-\mathbf{q}}\right)^{2}}}},$
(316)

where $\Theta(x)$ is the step function, $\epsilon^{e,h}_{1\mathbf{k}}$ are
defined by Eqs. (313) and (314), and the usual limits
$\Delta/\epsilon_{F},|\mathbf{q}|/k_{F}\ll 1$ are assumed.

In the weak-coupling limit under discussion, the SDW order can exist only if
the deviation from the perfect nesting is small, that is, $\alpha\ll 1$.
Taking this into account and using Eqs. (313) and (314), we can rewrite Eq.
(316) in the form [98]

$\ln{\frac{1}{\delta}}=\int_{0}^{2\pi}{\frac{d\varphi}{2\pi}Re\left\\{\cosh^{-1}\left[\frac{\nu_{0}(\mathbf{p},\varphi)-\nu}{\delta}\right]\right\\}},$
(317)

where

$\displaystyle\nu_{0}(\mathbf{p},\varphi)$ $\displaystyle=$ $\displaystyle
p_{x}\cos{\varphi}+p_{y}\sin{\varphi}-\frac{\bar{\alpha}}{2}\cos{2\varphi},\quad\bar{\alpha}=\frac{\alpha\kappa\epsilon_{F}}{\Delta_{0}},\quad\kappa=\frac{2\sqrt{v_{F}^{e}v_{F}^{h}}}{v_{F}^{e}+v_{F}^{h}},$
$\displaystyle\delta$ $\displaystyle=$
$\displaystyle\frac{\Delta}{\Delta_{0}},\quad\nu=\frac{\mu}{\kappa\Delta_{0}},\quad\mathbf{p}=\frac{\kappa
v_{F}}{2\Delta_{0}}\mathbf{q}.$ (318)

Similarly, we obtain the equation for the nesting vector
$\mathbf{Q}=\mathbf{Q}_{0}+2\Delta_{0}\mathbf{p}/v_{F}\kappa$:

$\left(\begin{array}[]{c}p_{x}\\\ p_{y}\\\
\end{array}\right)=\int_{0}^{2\pi}{\frac{d\varphi}{\pi}\left(\begin{array}[]{c}\sin{\varphi}\\\
\cos{\varphi}\\\
\end{array}\right)\textrm{sgn}{[\nu_{0}(\mathbf{p},\varphi)-\nu]}Re\sqrt{[\nu_{0}(\mathbf{p},\varphi)-\nu]^{2}-\delta^{2}}}.$
(319)

Equations (317) and (319) determine the SDW band gap and the nesting vector as
functions of $\mu$. Experiments are performed at fixed doping, and we have to
relate the electron density and $\mu$. Following the previously described
approach [see Eq. (298) and the text above it], we find in the approximation
under discussion [98]

$\frac{x}{x_{0}}=n\nu-\int_{0}^{2\pi}{\frac{d\varphi}{2\pi}\textrm{sgn}{[\nu_{0}(\mathbf{p},\varphi)-\nu]}Re\sqrt{[\nu_{0}(\mathbf{p},\varphi)-\nu]^{2}-\delta^{2}}},$
(320)

where $x_{0}$ is the normalized characteristic doping level, $n$ is the
density of electrons in nonmagnetic bands

$x_{0}=\frac{2v_{0}k_{F}^{2}\Delta_{0}}{\pi\kappa\epsilon_{F}},\quad
n=\frac{(v_{2F}^{h})^{-1}+(v_{3F}^{h})^{-1}+(v_{F}^{e})^{-1}}{(v_{F}^{h})^{-1}+(v_{2F}^{e})^{-1}},$
(321)

and $v_{0}$ is the unit cell volume. Note, that
$x\sim\Delta_{0}/\epsilon_{F}\ll 1$ and considered below phases, existing when
$x\sim x_{0}$, correspond to the low-doping regimes.

Figure 64: (Color online) Dependence of $\Delta$ and $|\mathbf{q}|$ on doping
$x$, $n=\kappa=1$. (a) $\bar{\alpha}=0.1$. The transition from commensurate to
incommensurate SDW state is of the first order: both $\Delta$ and
$|\mathbf{q}|$ rapidly change at $x=x_{1}^{*}$. Dashed lines show the
functions $\Delta(x)$ and $|\mathbf{q}(x)|$ corresponding to the metastable
state. (b) $\bar{\alpha}=0.6$. The transition from commensurate to
incommensurate SDW state is of the second order. Panels (c) and (d) show the
schematics of the incommensurate SDW spin structure. If $\alpha>0$, then, (c)
corresponds to the hole doping ($x<0$) and (d) to the electron doping ($x>0$).
If $\alpha<0$, then panel (c) corresponds to the electron doping, while panel
(d) corresponds to the hole doping [98].

The computed functions $\Delta(x)$ and $q(x)$ are shown in Figs. 64(a) and
64(b) for $n=\kappa=1$ and two different values of elipticity $\alpha$. The
numerical analysis reveals that the SDW order arises if $|\bar{\alpha}|<2$ or
$\alpha<\alpha_{c}=2\Delta_{0}/\kappa\epsilon_{F}\ll 1$. At low doping,
$|x|<x^{*}_{1}$, all the extra charge goes to the nonmagnetic bands. As a
result, the order parameter is independent of $x$, $q(x)$ = 0 (a commensurate
SDW order), and the chemical potential increases linearly with $x$. When
$|x|>x^{*}_{1}$, electrons (holes) appear in the band $E_{2}$ ($E_{1}$). The
SDW order becomes incommensurate, $q\neq 0$. If $n<n_{1}\approx 0.38$, the
commensurate–incommensurate transition is of the second order for any
$\bar{\alpha}$, while for larger $n$ it becomes the first order one if
$\bar{\alpha}<\bar{\alpha}_{1}(n)$. The spin configurations for the
incommensurate SDW are schematically shown in Figs. 64(c) and 64(d) (for
details, see Ref. [98]). For positive (negative) $\alpha$, Fig. 64(c)
corresponds to $x>0$ ($x<0$), and Fig. 64(d) corresponds to $x<0$ ($x>0$).

Figure 65: (Color online) Chemical potential $\mu(x)$ calculated at
$\bar{\alpha}=0.1$ (a), $\bar{\alpha}=0.5$ (b), and $\bar{\alpha}=1.2$ (c).
Here, $n=\kappa=1$. (a), (b) The homogeneous state is unstable toward the
phase separation if $x_{1}<x<x_{2}$. The dashed (red) curve corresponds to the
$\mu_{0}$ found using the Maxwell construction. The shaded areas above and
below $\mu_{0}$ are equal to each other. (c) $\mu(x)$ monotonically increases
with $x$, no phase separation occurs. The homogeneous commensurate ($q=0$) and
incommensurate ($q\neq 0$) SDW, paramagnetic (PM), and inhomogeneous
commensurate–incommensurate SDW (PS) states are separated by vertical dotted
lines [98].

The computed in Ref. [98] plots of the chemical potential $\mu(x)$ are shown
in Fig. 65 for three different values of elipticity $\bar{\alpha}$ and
$n=\kappa=1$. The $\mu(x)$ dependence for the hole doping can be determined
using the particle–hole symmetry. For low $\bar{\alpha}$, the function
$\mu(x)$ is nonmonotonic and multivalued near $x_{1}^{*}$, Fig. 65(a). At
higher $\bar{\alpha}$, the multivaluedness disappears; however, the
nonmonotonicity remains, Fig. 65(b). The latter vanishes at higher values of
$\bar{\alpha}$, Fig. 65(c). Thus, if the elipticity is not too large, there
are finite ranges of doping where the chemical potential decreases with doping
and the homogeneous state is unstable with respect to the separation into two
phases. However, the phase transitions between homogeneous commensurate and
incommensurate SDW phases, which we described above, can be masked, at least
partially, by the phase separation. The range of doping $x$, where the phase
separation exists, is the largest when $\alpha=0$, that is, when the model
under study is identical to the two-dimensional Rice model. The range of the
phase separation $x_{1}<x<x_{2}$ shrinks if the elipticity increases, and
disappears at the critical value $\bar{\alpha}_{c}=1.15$.

Figure 66: (Color online) Phase diagram of the model for pnictides in the
($x,\alpha$) plane, for $n=\kappa=1$ and $\alpha,x>0$. It is symmetric with
respect to the replacement $x\rightarrow-x$ and/or $\alpha\rightarrow-\alpha$.
The boundary between the incommensurate SDW and paramagnetic (PM) states shown
by the dashed (red) curve corresponds to the first order phase transition and
and that shown by the solid (red) curve corresponds to the second order phase
transition. Solid (blue) lines indicate the boundaries of the phase-separated
(PS) state. The solid (black) curve (second order transition) and dashed
(black) curve (first order transition) show the boundaries between
commensurate and incommensurate homogeneous SDW phases. Note that the phase
separation partially masks the transition line between the homogeneous SDW
states [98].

In the phase-separated state, there exist two phases, phase 1 with a lower
charge carrier density is the commensurate one and phase 2 with a higher
charge carrier density corresponds to the incommensurate SDW state. The phase
diagram of the system can be constructed following the approach described in
the previous subsections (see, e.g., subsection 8.2). The obtained results are
summarized in the phase diagram in the $(x,\alpha)$ plane shown in Fig. 66
[98]. This phase diagram is calculated for $n=\kappa=1$. It remains
qualitatively the same if $n\neq 0$. If the nonmagnetic bands are absent, the
homogeneous commensurate SDW phase exists only when doping is zero.
Consequently, the electronic density $x_{1}$ vanishes in phase 1 of the phase-
separated state for any $\alpha$.

The phase separation in iron-based superconductors was observed in several
experiments [365, 366, 367, 368]. For example, the inhomogeneous state with a
commensurate AFM and nonmagnetic domains with characteristic sizes about 65 nm
was observed in the hole-doped Ba1-xKxFe2As2 compound [365]. The described
above model predicts that the second phase is an incommensurate SDW rather
than a nonmagnetic one. However, the proposed mechanism of phase separation
can, in general, be consistent with the observations. The thermodynamic
potentials of the incommensurate SDW and the metastable PM phases are very
close to each other in the doping range $x>x_{2}$. The incommensurate SDW
phase can be destroyed by an additional cause not taken into account in the
model, e.g., by a disorder. In this case, the phase separation might occur
between the commensurate SDW and the PM phases.

Here, we consider the SDW order parameter in the form (291). As it was
mentioned in the previous section, in the literature a different order
parameter is also discussed (see, e.g., Ref. [358]). We must remember that the
relative stability of different order parameters is likely a non-universal
quantity, which depends on a variety of microscopic parameters (e.g., details
of the band structure, interaction, disorder). Therefore, the type of suitable
order parameter cannot be assuredly deduced without an input from experiments.

The calculations above show that there are only two possible equilibrium
directions of the vector $\mathbf{q}$ characterizing the incommensurate SDW
phase: it can be either parallel or perpendicular to the nesting vector
$\mathbf{Q}_{0}$, depending on the type of doping and the sign of the
ellipticity parameter $\alpha$. It is clear that both the magnitude and
direction of $\mathbf{q}$ are sensitive to the shape of the Fermi surface,
which is simplified in the present calculations. In actual pnictides, the
shape of the hole pockets deviates from perfect circles and the spectrum of
the bands depends on the transverse momentum $k_{z}$. However, the observation
of the incommensurate SDW phase with $\mathbf{q}$ perpendicular to
$\mathbf{Q}_{0}$ in the electron-doped Ba(Fe1-xCox)2As2 was reported in Refs.
[376, 377, 378] and double-peaked spin fluctuations corresponding to
$\mathbf{q}$ parallel to $\mathbf{Q}_{0}$ have been observed in the hole-doped
Ba1-xKxFe2As2 compound in a wide doping range [379, 380]. These observations
are in agreement with the prediction of the model.

Thus, a simple model approximating the Fermi surface of pnictides implies the
existence of electronic phase separation even in the weak-coupling regime.
This is an important finding for the interpretation of the experimental data
on phase inhomogeneity of iron pnictides: it proves that a purely electronic
model with moderate interaction is sufficient to explain the observed
inhomogeneities.

##  9 Phase separation in graphene-based systems

The formation of non-uniform electron states is a rather common feature for
many complicated systems. In this section, we show that the reviewed above
ideas of the phase separation can be implemented in the physics of graphene.
We consider two examples. The first one is the so-called AA graphene bilayer
and the second one deals with a graphene sheet with attached hydrogen atoms.
If the hydrogenization of graphene is complete, the latter system is usually
referred to as graphane. It turns out that AA graphene bilayer has the Fermi
surface with evident nesting, whereas the electron system of hydrogenized
graphene can be described in terms of the Falicov–Kimball model on the
honeycomb lattice.

###  9.1 AA graphene bilayer

Graphene is a zero-gap semiconductor exhibiting a plethora of unusual
electronic characteristics [381]. In addition to single-layer graphene,
bilayer graphene is also actively studied [96]. This interest is driven by the
desire to extend the family of graphene-like materials, which could be of
interest for applications. Bilayer graphene exists in three stacking
modifications. The most common is the so-called Bernal, or AB, stacking of
bilayer graphene. In such a stacking, half of the carbon atoms in the top
layer are located above the hexagon centers in the lower layer, and half of
the atoms in the top layer lie above the atoms in the lower layer. A different
layer arrangement, in which carbon atoms in the upper layer are located on top
of the equivalent atoms of the bottom layer, is referred to as AA-stacked
bilayer graphene (Fig. 67). The third modification is the twisted bilayer, in
which the top graphene layer is rotated by some angle with respect to the
bottom layer.

Figure 67: (Color online) Crystal lattice of the AA graphene bilayer consists
of two sublattices $A$ (red circles) and $B$ (blue circles). The unit cell of
the crystal structure includes four atoms $A1$ and $B1$ in the top (1) layer
and $A2$ and $B2$ in the bottom (2) layer. The values $t$ and $t_{0}$ are the
in-plane and intrerplane nearest nearest-neighbor hopping integrals,
respectively.

####  9.1.1 Single-electron spectrum

The unit cell of the single-layer graphene crystal lattice consists of two
carbon atoms. Consequently, two electron bands exist near the Fermi level
[381]. The unit cell of the bilayer graphene crystal includes four carbon
atoms, which forms sublattices $A1$ and $B1$ in the bottom (1) layer and $A2$
and $B2$ in the top (2) layers. The tight-binding analysis shows that both AA
and AB bilayer have four bands (two hole bands and two electron bands).
However, the structures of these bands are different [96]. In the undoped AB
bilayer, if we neglect the so-called trigonal warping, two bands touch each
other at two Fermi points, and the low-energy band dispersion is nearly
parabolic. The AA-stacked bilayer graphene has two bands near the Fermi
surface, one electron-like and another hole-like. The low-energy dispersion in
the AA bilayer is linear, similar to the monolayer graphene. Unlike the
latter, however, AA bilayer have Fermi surfaces instead of Fermi points (see
Fig. 68). An important feature of the AA bilayer is that the hole and electron
Fermi surfaces coincide in the undoped material, that is, there exist a
perfect nesting, which becomes imperfect with doping. According the results
reviewed in Section 8, in such a situation even small electron–electron
coupling gives rise to the SDW instability resulting in formation of the AFM
commensurate and incommensurate ordering and to the phase separation (see also
Refs. [350, 96, 382]).

In the tight-binding approximation, the single-electron part of the
Hamiltonian of AA-stacked graphene bilayer can be written as [96]

$\hat{H}_{0}=-t\\!\\!\\!\sum_{\langle\mathbf{nm}\rangle,i,\sigma}{\\!\\!\left(a^{\dagger}_{\mathbf{n}i\sigma}b_{\mathbf{m}i\sigma}\\!+\\!\textrm{H.c.}\right)}-t_{0}\sum_{\mathbf{n},\sigma}{\\!\\!\left(a^{\dagger}_{\mathbf{n}1\sigma}a_{\mathbf{n}2\sigma}+b^{\dagger}_{\mathbf{n}1\sigma}b_{\mathbf{n}2\sigma}\\!+\\!\textrm{H.c.}\right)}-\mu\\!\sum_{\mathbf{n},i,\sigma}{\\!\left(a^{\dagger}_{\mathbf{n}i\sigma}a_{\mathbf{n}i\sigma}+b^{\dagger}_{\mathbf{n}i\sigma}b_{\mathbf{n}i\sigma}\right)},$
(322)

where $a^{\dagger}_{\mathbf{n}i\sigma}$ ($b^{\dagger}_{\mathbf{n}i\sigma}$) is
the creation operator of an electron in sublattice $a$ ($b$) at site
$\mathbf{n}$ in layer $i=1,2$ with spin projection $\sigma$, and the values
$t$ and $t_{0}$ denote the in-plane and interplane nearest-neighbor hopping
integrals, respectively. In graphene, next-nearest-neighbor hopping amplitudes
are much smaller than $t$ and $t_{0}$ and hence, a longer-range hoping can be
neglected [381, 96]. The energy spectrum of Hamiltonian (322) consists of four
bands, which can be expressed in the form [96]

$\varepsilon^{1,2}_{\mathbf{k}}=-t_{0}\pm
t|f_{\mathbf{k}}|,\quad\varepsilon^{3,4}_{\mathbf{k}}=t_{0}\pm
t|f_{\mathbf{k}}|,$ (323)

where

$f_{\mathbf{k}}=1+\exp{\\!\\!\left(\frac{3ik_{x}a_{0}}{2}\right)}\cos{\\!\\!\left(\frac{\sqrt{3}k_{y}a_{0}}{2}\right)},$
(324)

and $a_{0}$ is the lattice constant for single-layer graphene.

The single-electron bands of the AA graphene bilayer (323) are shown in Fig.
68. The bands $s=2$ and 3 intersect the Fermi surface in the vicinity of the
Dirac points $\mathbf{K}=2\pi(\sqrt{3},1)/3\sqrt{3}a_{0}$ and
$\mathbf{K^{\prime}}=2\pi(\sqrt{3},-1)/3\sqrt{3}a_{0}$ [see Fig. 68(b) and
(c)]. When the doping is zero ($\mu=0$, half filling), the Fermi energy is
determined by the equation $|f_{\mathbf{k}}|=t_{0}/t$. Under condition
$t_{0}/t\ll 1$, the function $|f_{\mathbf{k}}|$ can be expanded near the Dirac
points. Then, we obtain that the Fermi surface consists of two circles with
radius $k_{r}=2t_{0}/(3ta_{0})$. In the case of a doped system, these Fermi
surfaces transform into four circles. It is important that at half filling,
the Fermi surfaces of the electron and hole bands are perfectly nested. This
property is robust against changes in the Hamiltonian. For example, it
survives if we take into account the longer-range hopping or consider a system
with two nonequivalent layers [382].

Figure 68: (Color online) (a) Single-electron bands in the AA graphene
bilayer. (b) The momentum dependence of the energy spectra
$\varepsilon^{s}_{\mathbf{k}}$ near the Dirac point located at $\mathbf{K}$.
Here, $\mathbf{k}=\mathbf{K}+\delta k_{y}\mathbf{e}_{y}$. The bands $s=2$ and
3 intersect each other exactly at zero energy, which corresponds to the Fermi
energy of the undoped graphene. (c) The first Brillouin zone (hexagon) and the
reciprocal-lattice unit cell (rhombus) of the AA-stacked graphene bilayer. The
circles near $\mathbf{K}$ and $\mathbf{K}^{\prime}$ Dirac points show the
Fermi surfaces of the doped sample [350]; $\varepsilon^{(s)}_{0\mathbf{k}}$ in
this figure corresponds to $\varepsilon^{s}_{\mathbf{k}}$ in the present
notation.

####  9.1.2 Symmetry-breaking phases

The spectrum with two bands having identical Fermi surfaces is unstable with
respect to spontaneous symmetry breaking in the presence of a small
electron–electron interaction. Such an instability opens a gap in the
electronic spectrum, which gives rise to a decrease of the free energy of the
system. We outline below only theoretical results since the experimental data
on the electronic properties of the AA graphene bilayer now is nearly absent.
The symmetries of the Hamiltonian can be used to make narrower the choices of
the possible symmetry breakings and order parameters [382]. We can find the