\chapter{The Hubbard model}
\label{chap:HubbardModel}

\section{Introduction}
\label{sec:hubbard-model}
The Hubbard model was first proposed in 1963 independently by Hubbard
\cite{Hubbard:1963} and Gutzwiller \cite{Gutzwiller:1963} in an attempt to
describe the effects of correlations for d electrons in metals.  Transition
metals where the d-band is partially full display characteristics which in
some circumstances can be best understood in terms of localized atomic-like
electrons but in other cases the experimental features are better described in
terms of band theory.  The Hubbard model is essentially the simplest model
that captures this duality and it can be derived from the full electronic
Hamiltonian for the metal by keeping only interaction terms between electrons
that are situated on the same site, or more specifically between electrons of
different spins which are in the same Wannier state.  For an s-band the
Hubbard Hamiltonian is given by:
\begin{equation}
H = -\sum_{\langle ij\rangle\sigma}t_{ij}c_{i\sigma}^{+}c_{j\sigma}
+ U\sum_{j}n_{j\uparrow}n_{j\downarrow}. 
\end{equation}
Here $c_{i\sigma}^{+}$ creates an electron in a Wannier state localized at
site $i$ with spin $\sigma$, $U$ is the local on site interaction between
different spins and $t_{ij}$ is the hopping integral between sites $i$ and
$j$.  For our purposes we will want to generalize this Hamiltonian by
increasing the spin degeneracy from 2 to $N$ and since we will be working in
the grand canonical ensemble we introduce a chemical potential that we add to
the Hamiltonian.  Thus we have:
\begin{equation}
\label{eq:Hubbard-Hamiltonian}
H = -\sum_{\langle ij\rangle\sigma}t_{ij}c_{i\sigma}^{+}c_{j\sigma}
+ {U \over 2}\sum_{j\sigma\neq\sigma'}
n_{j\sigma}n_{j\sigma'}
-\mu\sum_{j\sigma}n_{j\sigma}.
\end{equation}
Here $\sigma$ runs from $-j$ to $j$ where $2j+1= N$.  This model has a higher
spin degeneracy than the original Hamiltonian of Hubbard and can thus be
expected to be a better representation of the physics in the materials of
interest, i.e. compounds of transition and rare earth metals which have partly
filled d and f bands. 
% It should be noted though that here it is assumed that the interaction is the
% same between electrons of different total angular momenta irrespective of what
% they might be.  However in real materials, electrons of opposite spin in the
% same orbital usually interact more strongly than electron in different
% orbitals, but this complication will be neglected.
In spite of being one of the simplest model of correlated electrons on a
lattice the Hubbard model has proven to be quite difficult to tackle
theoretically.  It is only during the last few years that some real progress
has been made, especially with the development of the dynamical mean field
theory \cite{Georges:1996}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{An overview of dynamical mean field theory}
\label{sec:dmft}
The dynamical mean field theory can be regarded as the quantum analog of the
classical Curie-Weiss mean field theory for magnets.  As its classical
counterpart it can be shown to be the exact formulation of the model in
question in the limit when the number of dimensions or the coordination number
of the lattice goes to infinity, with an appropriate scaling of parameters.
When the number of neighbors of a given lattice site increases it is
reasonable to assume that its interactions with the neighbors can be replaced
by an effective interaction with some average field that represents the
dynamic state of the neighbors.  These ideas can be concretely realized in
the following way: A given lattice site is selected, $o$ say, and an effective
action is derived for the degrees of freedom living on this site by
integrating out the degrees of freedom of the other sites in the lattice.  The
effective action can be represented in terms of the Green's functions,
$G_{ij}^{0}$, of the so called cavity Hamiltonian, i.e. the lattice
Hamiltonian obtained by deleting the site $o$ and all the bonds associated
with it.  By scaling the hopping integral such that $z^{|i-j|}t_{ij}^{2}$ is a
constant, where $z$ is the coordination number of the lattice and $|i-j|$ is
the Taxi-cab distance between the sites $i$ and $j$, and then taking the limit
when $z \rightarrow \infty$ the action simplifies tremendously and we obtain:
\begin{equation}
\label{eq:effective-action}
S_{eff} = -\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau'
\sum_{\sigma}c_{o\sigma}^{+}\mathcal{G}^{-1}(\tau-\tau')c_{o\sigma}
+{U \over 2}\int_{0}^{\beta}d\tau\sum_{\sigma\neq\sigma'}
n_{o\sigma}n_{o\sigma'}
\end{equation}
where the function $\mathcal{G}$ is the so called Weiss function and plays the
role of the Weiss field of the Curie-Weiss theory and is given by,
\begin{equation}
\label{eq:bare-prop}
\mathcal{G}^{-1}(i\omega_{n}) = i\omega_{n}+\mu-\sum_{\langle ij \rangle}
t_{io}t_{jo}G_{ij}^{0}(i\omega_{n}).
\end{equation}
Given an explicit form for the action, the on-site interacting Green's
function $G$ can in principle be calculated (in practice this is a nontrivial
problem) by evaluating the following expectation value:
\begin{equation}
G(\tau-\tau') = -\langle
Tc_{o}(\tau)c_{o}^{+}(\tau')\rangle_{S_{eff}}.
\end{equation}
The self consistency condition that is now imposed is that this Green's
function be the same as the diagonal, on site lattice Green's function, $G_{oo}
\equiv G_{loc}$.  To compute the full, lattice Green's function it is useful
to switch over to momentum representation and introduce a self energy $\Sigma$
which is momentum independent in the infinite dimensional limit and given by,
\begin{equation}
\label{eq:self-energy}
\Sigma(i\omega_{n}) = \mathcal{G}^{-1}(i\omega_{n})-G^{-1}(i\omega_{n}).
\end{equation}
Then the momentum dependent Green's functions are given by,
\begin{equation}
G(k,i\omega_{n}) = 
{1 \over i\omega_{n}+\mu-\varepsilon_{k}-\Sigma(i\omega_{n})}
\end{equation}
where $\varepsilon_{k}$ is the dispersion relation of the lattice. The general
inter-site lattice Green's function can be computed from this expression by
Fourier transform and the diagonal term is of course simply the momentum
integral of the expression above.  At this point it is useful to stop and look
over the equations that we have so far.  Let us for a moment assume that we
have a guess for the self energy, $\Sigma$.  From there we can compute the
local Green's function by momentum integration over $G(k,i\omega)$ and then
the Weiss field can be computed from equation (\ref{eq:self-energy}).  We now
have all the ingredients to compute the local Green's function from the
action, and then we can use equation (\ref{eq:self-energy}) again to get a new
solution for the self energy.  This process can now be started over and
iterated until self consistency is attained.  Thus it seems that we have all
that is necessary to solve our problem but we still haven't used the cavity
Green's function, $G^0$.  In the limit $z\rightarrow\infty$ the cavity Green's
function, $G_{ij}^{0}$, is very simply related to the full lattice Green's
function , $G_{ij}$.  This relation is,
\begin{equation}
\label{eq:cavity-relation}
G_{ij}^{0} = G_{ij}-{G_{io}G_{jo} \over G_{oo}}
\end{equation}
and it is quite simple to understand: using the interpretation of the Green's
function as the likelihood of finding an electron at site $j$ if it is
created at $i$ we see that this likelihood is the same in the full lattice as
in the lattice with a cavity at $o$, except that in the cavity lattice we
cannot go through the site $o$ and thus we must subtract those contributions
off after normalization for processes that start at $o$ and come back to $o$.
Putting this relation into the equation (\ref{eq:bare-prop}) and using the
form for the momentum dependent Green's functions we end again up with
equation (\ref{eq:self-energy}).  So, somehow this relation between the cavity
Green's function and the full Green's function seems to be equivalent to the
momentum independence of the self energy.  Here it is useful to use the
representation of the local Green's function in terms of the density of
states:
\begin{equation}
G(i\omega_{n}) = \int {D(\epsilon)d\epsilon \over \zeta-\epsilon}
\end{equation}
where $D(\epsilon)$ is the bare density of states of the lattice and
$\zeta = i\omega_{n}+\mu-\Sigma(i\omega_{n})$.
We see that $G$ is essentially given by a Hilbert transform over the 
density of states and if we can invert that we can express the self 
energy in term of the Green's function and then the Weiss function can
be expressed solely in terms of the Green's function through equation
(\ref{eq:self-energy}).

In the calculations that follow we will use the Bethe lattice of infinite
connectivity with only nearest neighbor hopping given by $t_{o}$, $zt_{o}^{2}
= t^{2}$.  In that case the self consistency condition becomes especially
simple and we get,
\begin{equation}
\mathcal{G}^{-1}(i\omega) = i\omega+\mu-
t^{2}G_{loc}(i\omega)-G^{-1}_{loc}(i\omega).
\end{equation}
Given the physical picture embodied in the DMFT, i.e. of an electron on a
single site interacting with a bath of electrons on the rest of the lattice it
is not surprising that a connection can be made between the dynamical mean
field theory described above and impurity models like the Anderson impurity
model,
\begin{equation}
H^{AIM} = \sum_{k\sigma}\tilde{\varepsilon}_{k}
c_{k\sigma}^{+}c_{k\sigma}-\mu\sum_{\sigma}
f_{\sigma}^{+}f_{\sigma}+\sum_{k\sigma}
(V_{k}c_{k\sigma}^{+}f_{\sigma}+h.c.) + 
{U \over 2}\sum_{\sigma\neq\sigma'}n_{\sigma}n_{\sigma'}
\end{equation}
by noting that if the conduction electrons $c_{k\sigma}$ are integrated out an
action of the form (\ref{eq:effective-action}) is obtained with
\begin{equation}
\mathcal{G}^{-1}(i\omega) = i\omega+\mu-
\int d\varepsilon {\Delta(\varepsilon) \over i\omega-\varepsilon}
\quad\mbox{where}\quad
\Delta(\varepsilon) = \sum_{k}V^{2}_{k}
\delta(\varepsilon-\tilde{\varepsilon}_{k})
\end{equation}
is the hybridization function of the Anderson model.  The full
impurity Green's function will now be equal to $G_{loc}$ provided
$\Delta$ is chosen such that the self-consistency condition is 
fulfilled.  On the Bethe-lattice that requirement gives us
\begin{equation}
\Delta(\varepsilon) = t^{2}\rho_{loc}(\varepsilon) = 
-{t^{2} \over \pi}\imag G_{loc}(\varepsilon+i0).
\end{equation}
The impurity model connection is very fortunate since impurity models 
have been studied intensely in the past, their physics is well known 
and many tools and methods exist to deal with such models such as quantum
Monte-Carlo methods, exact diagonalization methods, iterated perturbation
methods and the so called non-crossing approximation (NCA) which we will
mainly focus on here.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{The non-crossing approximation}
\label{sec:nca}
The NCA solution to the Anderson model applies specifically to the case where
the interaction strength $U$ goes to infinity which prevents more than one
electron from occupying each site in the lattice.  Then the interaction term in
the Anderson Hamiltonian goes away but is replaced by the auxiliary condition
that the occupancy of the impurity be less than or equal to one.  In this case
the Hamiltonian can be rewritten in the form \cite{Coleman:1984},
\begin{equation}
H_{AIM} = \sum_{k\sigma}\tilde{\varepsilon}_{k}c_{k\sigma}^{+}c_{k\sigma}
-\mu\sum_{\sigma}f_{\sigma}^{+}f_{\sigma}
+\sum_{k,\sigma}[V_{k}c_{k\sigma}^{+}b^{+}f_{\sigma} + h.c.]
\end{equation}
where $b$ is a bosonic operator that represents the empty impurity and
$f_{\sigma}$ is a pseudo-fermion with the original fermion now being
represented by the combination $f_{\sigma}^{+}b$.  We can think of it in such
a way that in order to create an electron at an empty impurity site we first
need to destroy the empty impurity and then create the electron from some new
unphysical vacuum state.  The physical subspace is the one in which the number
operator, $Q = b^{+}b+\sum_{\sigma}f^{+}f$ takes the value 1.  This condition
can be imposed by introducing a chemical potential associated with $Q$ and
then let it tend to infinity in order to project out the physical subspace.
The non-crossing approximation to the solution of the impurity problem
consist, as the name implies, in calculating the self energies for the boson
$b$ and the fermion $f$ by summing up all the Feynman diagrams that do not
have any crossing lines.  This procedure leads to the equations\footnote{All
  Green's functions and self energies are assumed to be retarded unless
  explicitly indicated otherwise.}:
\begin{alignat}{2}
G_{b}(\omega) &= {1 \over \omega - \Sigma_{b}(\omega)},
& \quad 
\Sigma_{b}(\omega) &= N\int_{-\infty}^{\infty}d\varepsilon
\Delta(\varepsilon)f(\varepsilon)G_{f}(\omega+\varepsilon), 
\\
G_{f}(\omega) &= {1 \over \omega + \mu - \Sigma_{f}(\omega)},
& \quad 
\Sigma_{f}(\omega) &= \int_{-\infty}^{\infty}d\varepsilon
\Delta(\varepsilon)(1-f(\varepsilon))G_{b}(\omega-\varepsilon).
\end{alignat}
where $f(\varepsilon)$ is the Fermi function.
The impurity Green's functions is now given in terms of the auxiliary 
Green's functions by:
\begin{equation}
G_{imp}(\omega) = \int_{-\infty}^{\infty} d\varepsilon
\left(
\eta_{b}(\varepsilon)G_{f}(\varepsilon+\omega) 
- \eta_{f}(\varepsilon)G_{b}^{\star}(\varepsilon-\omega)
\right)
\end{equation}
Here the $\eta$'s are given by:
\begin{equation}
\eta_{f}(\omega) = {1 \over Z}e^{-\beta\omega}\rho_{f}(\omega),
\quad\quad
\eta_{b}(\omega) = {1 \over Z}e^{-\beta\omega}\rho_{b}(\omega).
\end{equation}
with
\begin{equation}
\rho_{f}(\omega) = -{1 \over \pi}\imag G_{f}(\omega), \quad\quad
\rho_{b}(\omega) = -{1 \over \pi}\imag G_{b}(\omega).
\end{equation}
and $Z$ is the impurity partition function.
To avoid numerical problems at negative energies when $\eta$
is calculated from $\rho$ it is beneficial to have explicit NCA 
equations for the $\eta$'s.  These are easily found from the
original NCA equations and turn out to be
\begin{eqnarray}
\phi_{b}(\omega) \equiv \eta_{b}(\omega)|G_{b}(\omega)|^{-2} &=& 
N\int_{-\infty}^{\infty}d\varepsilon 
\Delta(\varepsilon)(1-f(\varepsilon))\eta_{f}(\omega+\varepsilon)
\\
\phi_{f}(\omega) \equiv \eta_{f}(\omega)|G_{f}(\omega)|^{-2} &=& 
\int_{-\infty}^{\infty}d\varepsilon 
\Delta(\varepsilon)f(\varepsilon)
\eta_{b}(\omega-\varepsilon)
\end{eqnarray}
Note that the equations for $\eta_{b}$ and $\eta_{f}$
determine these functions only up to a constant and thus we need
some additional relation to determine them uniquely.  This is 
furnished by the following equation for the partition function:
\begin{equation}
Z = \int_{-\infty}^{\infty}d\varepsilon
\left(
\rho_{b}(\varepsilon) + N\rho_{f}(\varepsilon)
\right)e^{-\beta\varepsilon}
\end{equation}
from which we obtain
\begin{equation}
\int_{-\infty}^{\infty}d\varepsilon
\left(
\eta_{b}(\varepsilon) + N\eta_{f}(\varepsilon)
\right) = 1.
\end{equation}
For numerical calculations the following sum rules are quite useful as a check
on the validity of the solution \cite{Bickers:1987}:
\begin{equation}
\label{eq:sumrule-1}
\int_{-\infty}^{\infty}\rho_{b}(\varepsilon)d\varepsilon = 1,
\qquad
\int_{-\infty}^{\infty}\rho_{f}(\varepsilon)d\varepsilon = 1,
\end{equation}
\begin{equation}
\label{eq:sumrule-2}
\int_{-\infty}^{\infty}\rho_{imp}(\varepsilon)f(\varepsilon)
d\varepsilon = {n \over N},
\qquad
\int_{-\infty}^{\infty}\rho_{imp}(\varepsilon)(1-f(\varepsilon))
d\varepsilon = 1-n,
\end{equation}
where $n$ is the occupancy of the impurity.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Numerical implementation of the NCA+DMFT}
\label{sec:num-calc}
Before starting the numerical calculations it is useful to 
define new functions as follows: 
\begin{alignat*}{2}
\mytilde{G}_{f}(\omega) &= G_{f}(\omega-\mu),
& \quad
\mytilde{G}_{b}(\omega) &= G_{b}(\omega),
\\
\mytilde{\Delta}(\omega) &= \Delta(\omega-\mu),
& \quad
\mytilde{G}_{loc}(\omega) &= G_{loc}(\omega-\mu),
\end{alignat*}
with the same convention for the other boson and pseudo-fermion functions
($\phi,\rho,\eta$).  In terms of these functions the set of equations becomes:
\begin{alignat}{2} 
\mytilde{G}_{b}(\omega) &= 
{1 \over \omega - \mytilde{\Sigma}_{b}(\omega)}, 
& \quad 
\mytilde{\Sigma}_{b}(\omega) &= N\int_{-\infty}^{\infty}d\varepsilon
\mytilde{\Delta}_{p}(\varepsilon)
\mytilde{G}_{f}(\omega+\varepsilon), 
\\
\mytilde{G}_{f}(\omega) &= 
{1 \over \omega - \mytilde{\Sigma}_{f}(\omega)},
& \quad
\mytilde{\Sigma}_{f}(\omega) &= \int_{-\infty}^{\infty}d\varepsilon
\mytilde{\Delta}_{h}(\varepsilon)
\mytilde{G}_{b}(\omega-\varepsilon),
\end{alignat}
\begin{eqnarray}
\mytilde{\phi}_{b}(\omega) \equiv
\mytilde{\eta}_{b}(\omega)
|\mytilde{G}_{b}(\omega)|^{-2} &=& 
N\int_{-\infty}^{\infty}d\varepsilon 
\mytilde{\Delta}_{h}(\varepsilon)
\mytilde{\eta}_{f}(\omega+\varepsilon)
\\
\mytilde{\phi}_{f}(\omega) \equiv
\mytilde{\eta}_{f}(\omega)
|\mytilde{G}_{f}(\omega)|^{-2} &=& 
\int_{-\infty}^{\infty}d\varepsilon 
\mytilde{\Delta}_{p}(\varepsilon)
\mytilde{\eta}_{b}(\omega-\varepsilon)
\end{eqnarray}
\begin{equation}
\mytilde{G}_{imp}(\omega) = \int_{-\infty}^{\infty} d\varepsilon
\left(
\mytilde{\eta}_{b}(\varepsilon)\mytilde{G}_{f}(\varepsilon+\omega) 
- \mytilde{\eta}_{f}(\varepsilon)
\mytilde{G}_{b}^{\star}(\varepsilon-\omega)
\right)
\end{equation}
and 
\begin{equation}
\mytilde{\Delta}(\omega) = t^{2}\mytilde{\rho}_{imp}(\omega)
\end{equation}
Here we have defined
\begin{equation}
\mytilde{\Delta}_{p}(\varepsilon) = 
\mytilde{\Delta}(\varepsilon)f(\varepsilon-\mu)
\quad\mbox{and}\quad
\mytilde{\Delta}_{h}(\varepsilon) = 
\mytilde{\Delta}(\varepsilon)[1 - f(\varepsilon-\mu)].
\end{equation}
The equations above constitute a closed set of self consistent equations for
the impurity Green's function $G_{imp} = G_{loc}$.  The energy scale of the
model is set by the hopping integral, $t$, which we will in all our
calculations set equal to ${1 \over 2}$ which corresponds to a bare DOS of
half-width $D = 1.0$ on the Bethe-lattice.  We solve these equations
numerically by cutting off the integrals and then introduce a uniform mesh of
points between the cutoffs.  Traditionally the NCA equations are solved on
logarithmic meshes \cite{Cox:thesis,Pruschke:1989,Anders:1995} due to
singularities in the spectrum but here these singularities are less pronounced
due to the smoothing effects of the mean field consistency condition.  For
this reason and also to be able to use fast Fourier transform techniques we
use uniform meshes.  Because of this we could not, however, go to very low
temperatures in our calculations but as we will see later this is not the
limiting factor in the usefulness of the method.

On a uniform mesh the convolution integrals can either be done by
straightforward integration or by the Fast Fourier Transform.  We actually
implemented both methods.  The FFT method is of course faster but it turned
out that the direct method gave better results as determined by the
convergence of the chemical potential (or thermopower) to the value given by
Heikes formula (see \article\ref{sec:res-transport}).  The upshot was that the
FFT method is used first to find an approximate solution for a given
temperature and density and then the direct method is used to improve upon
this solution.

For the FFT solution we define our Fourier transforms as follows.
\begin{equation}
g(\omega) = \int_{-\infty}^{\infty}dte^{i\omega t}g(t), \quad\quad
g(t) = \int_{-\infty}^{\infty}{d\omega \over 2\pi}
e^{-i\omega t}g(\omega).
\end{equation}
After having imposed an energy cutoff, $\omega_{max}$, the Fourier 
transforms are approximated by
\begin{equation}
g(t) \simeq \int_{-\omega_{max}}^{\omega_{max}}{d\omega \over 2\pi}
e^{-i\omega t}g(\omega).
\end{equation}
and the cutoff Fourier integral is then calculated by using cubic 
splines and the Fast Fourier Transform as described in \cite{Press:NR} and 
Appendix \ref{app:transforms}.
After the Fourier transformation the equations for the self energies
are
\begin{equation}
\mytilde{\Sigma}_{b}(t) = 
2N\pi\mytilde{\Delta}_{p}^{*}(t)\mytilde{G}_{f}(t), \quad\quad
\mytilde{\Sigma}_{f}(t) = 
2\pi\mytilde{\Delta}_{h}(t)\mytilde{G}_{b}(t).
\end{equation}
To take care of the long ${1 \over \omega}$ tails of the Green's
functions we also subtract off the real part of an almost free 
Green's functions\footnote{We only subtract off the real part since
the imaginary part of the Green's function already falls off faster
than ${1 \over \omega^{2}}$.} 
\begin{equation}
G_{free}^{0}(\omega) = {1 \over \omega+i\Gamma}
\end{equation}
before doing the transform and then add by hand its Fourier transforms
which is equal to half of the Fourier transforms of the Green's
function above,
\begin{equation}
G_{free}^{0}(t) = -i\theta(t)e^{-\Gamma t}.
\end{equation}
Here $\Gamma > 0$ is a quantity of the order of the mesh spacing.  Note that
all the functions we are dealing with are retarded so they are zero at
negative times.  For the auxiliary Green's functions we have the condition
that $G_{b}(t = 0^{+}) = G_{f}(t = 0^{+}) = -i$, and in the numerical
calculation these values are put in by hand.

The direct integration is much simpler to implement.  We approximate 
the functions by linear interpolation between mesh points and assume 
that they are identically zero outside of the cut-off which is chosen 
the same as for the FFT solution. The convolutions of the
interpolating functions are then done analytically.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Results: Green's functions and self energies}
\label{sec:res-green-self}
In figures \ref{fig:greenfig1} and \ref{fig:greenfig2} we have displayed the
results from the numerical calculations of the local Green's functions as
described in the previous section.  The calculations were done for several
densities between zero and one, for degeneracies $N = 2,4,6$ and temperatures
in the range of 0.025-5.0 in most cases\footnote{Remember that all energies
  are measured in terms of the half-bandwidth $D$ unless otherwise
  specified.}.
%----------------------------
\begin{figure}[htb]
\centering
\begin{tabular}{c@{\quad}c}
\mbox{\epsfig{file=Figures-2/imgreent0.05d2.eps,width=0.45\linewidth}}&
\mbox{\epsfig{file=Figures-2/imgreent0.05d4.eps,width=0.45\linewidth}}\\
\mbox{\epsfig{file=Figures-2/imgreent0.05d6.eps,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-2/chempotst0.05.eps,width=0.45\linewidth}}
\end{tabular}
\caption{Results from NCA+DMFT calculations.  I: imaginary part
  of the Green's functions for different densities, $n$, for degeneracy $N =
  2$ and temperature $T = 0.05$. II: Same as I for $N = 4$.  III: Same as I
  for $N = 6$.  IV: The chemical potential corresponding to the Green's
  functions in I-III.  }
\label{fig:greenfig1}
\end{figure}
%----------------------------
\begin{figure}[htb]
\centering
\begin{tabular}{c@{\quad}c}
\mbox{\epsfig{file=Figures-2/imgreen-n0.90d2.eps,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-2/imgreen-n0.90d4.eps,width=0.45\linewidth}} \\
\mbox{\epsfig{file=Figures-2/imgreen-n0.90d6.eps,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-2/chempots-n0.90.eps,width=0.45\linewidth}}
\end{tabular}
\caption{Results from NCA+DMFT calculations.  I: imaginary part
  of the Green's functions for different temperatures, $T$, for degeneracy $N
  = 2$ and density $n = 0.90$. II: Same as I for $N = 4$.  III: Same as I for
  $N = 6$.  IV: The chemical potential corresponding to the Green's functions
  in I-III.  }
\label{fig:greenfig2}
\end{figure}
%----------------------------
The two main features of the spectra are the broad band corresponding to the
incoherent lower Hubbard band and a quasiparticle peak or Kondo peak just
above the Fermi-surface.  There is no sign of an upper Hubbard band since the
interaction strength, $U$, is infinite and thus no charge excitations exist in
the system -- the upper Hubbard band has essentially moved off to infinity.
But it takes a lot of spectral weight with it as we can see by computing the
total integrated spectral weight of the system from the sum rules in equation
(\ref{eq:sumrule-2}):
\begin{equation}
\int_{-\infty}^{\infty}\rho_{loc}(\varepsilon)d\varepsilon = 
1-n+{n \over N}.
\end{equation}
As the density increases more and more spectral weight is lost and at half
filling, i.e. $n = 1$ it is down to $1/N$.  We note also that as the
degeneracy of the system, $N$, increases the total spectral weight decreases.
From the data we see however that this spectral weight is mostly taken from
the incoherent band -- both the width and the height of the quasiparticle peak
show much weaker dependence on the degeneracy.

If we look at figure \ref{fig:greenfig1} we note that for the lowest values of
density the spectral density is quite close to being equal to the free
electron spectral density on the Bethe lattice.  As the density increases we
start seeing the formation of the quasiparticle peak.  The quasiparticle peak
corresponds to carriers that have a well defined momentum and are thus in a
delocalized low entropy state.  In order to form this state the carriers must
get rid of their entropy and thus we can expect this state to appear earlier
in the lower degeneracy systems that have a smaller entropy per carrier.  This
is exactly what can be observed in the data (see figure \ref{fig:greenfig1})
-- the quasiparticle peak is already evident in the lowest density curve for
$N=2$ whereas only a slight shoulder can be seen for the corresponding density
in the $N=6$ system.

From figure \ref{fig:greenfig2} we see that as the temperature increases the
Kondo peak melts away, the weight is transferred into the Hubbard band which
becomes increasingly more symmetric around $\omega = -\mu$ as the temperature
increases.

From the NCA equations it is easy to see that if we start the DMFT iterations
with a hybridization function $\mytilde{\Delta}$, and initial guesses for
$\mytilde{G}_{f}$ and $\mytilde{G}_{b}$ which are symmetric around zero then
the self consistent solutions for $\mytilde{G}_{f}$ and $\mytilde{G}_{b}$ are
also symmetric around zero to lowest order in $\beta$ and therefore the self
consistent solution, $\mytilde{G}_{loc}$, to the mean field equations will be
symmetric to lowest order in $\beta$.  To find how the chemical potential
depends on temperature we use the sum rules (\ref{eq:sumrule-2}) and write:
\begin{eqnarray}
{n \over N} &=& \int_{-\infty}^{\infty}\rho_{imp}(\varepsilon)f(\varepsilon)
d\varepsilon \simeq A f(-\mu) 
\nonumber \\
1 - n  &=& \int_{-\infty}^{\infty}\rho_{imp}(\varepsilon)f(-\varepsilon)
d\varepsilon \simeq A f(\mu)
\nonumber 
\end{eqnarray}
from which we get,
\begin{equation}
{\mu \over T} = \ln({n \over N(1-n)}).
\end{equation}
Thus the chemical potential is a linear function of $T$ at large $T$.  This
formula can clearly be regarded as the result of the competition between the
empty site with density $1-n$ and the occupied site with density per spin
$n/N$. If the empty site density is higher than the density per spin the
chemical potential is forced to become negative and vice versa.  This equation
can also derived very simply by considering the high temperature entropy of a
system of $N_{e}$ electrons with degeneracy $N$ on a lattice with $N_{s}$
sites, only one electron per site \cite{Chaikin:1976}.  The partition function
is obviously
\begin{equation}
Z = {N^{N_{e}}N_{s}! \over N_{e}!(N_{s}-N_{e})!}
\end{equation}
and using the Stirling approximation we find that the entropy density is
\begin{equation}
\label{eq:entropy}
S = {\ln(Z) \over N_{s}} = -n\ln(n/N)-(1-n)\ln(1-n)
\end{equation}
where $n = {N_{e} \over N_{s}}$ is the density.  Now using the fact that ${\mu
  \over T} = {\partial S \over \partial n}$ we obtain the formula for $\mu$
which is given above.  If we look at the self energy for the same parameters
as in figure \ref{fig:greenfig2} we obtain the graphs in figure
\ref{fig:sigmafig1}.  The main feature to observe is the sharp dip at low
temperatures around the Fermi-level.  This correspond to the quasiparticle
peak in the local spectral density and as we can see the self energy spectral
density actually becomes negative at that point if the temperature is low
enough.
%----------------------------
\begin{figure}[htb]
\centering
\begin{tabular}{c@{\quad}c}
\mbox{\epsfig{file=Figures-2/imsigma-n0.90d2.eps,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-2/imsigma-n0.90d4.eps,width=0.45\linewidth}} \\
\mbox{\epsfig{file=Figures-2/imsigma-n0.90d6.eps,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-2/tbreak.eps,width=0.45\linewidth}}
\end{tabular}
\caption{Results from NCA+DMFT calculations.  I: imaginary part
  of the self energy for different temperatures, $T$, for degeneracy $N = 2$
  and density $n = 0.90$. II: Same as I for $N = 4$.  III: Same as I for $N =
  6$.  IV: The breakdown temperature for different densities (see main text).
  }
\label{fig:sigmafig1}
\end{figure}
%----------------------------
This pathology is of course quite serious since it prevents the calculation of
certain quantities like the transport parameters at temperatures below certain
breakdown temperature at which the imaginary part of the self energy becomes
zero.  The breakdown temperature depends both on the density and degeneracy as
we can see in figure \ref{fig:sigmafig1}(IV).  The question that comes
naturally up is; how does this pathology occur, and how is it that one can
solve the model and obtain the spectral density if the self energy, through
which the lattice self consistency proceeds, is so pathological.  To see this
we need to look at how the self consistency and the calculation of the self
energy is carried out on the Bethe lattice:
\begin{equation}
\Delta(\omega) = t^{2}\rho_{loc}(\omega), \qquad
\imag\Sigma(\omega) = \imag G_{loc}(\omega)\left(
{1 \over |G_{loc}(\omega)|^{2}} - t^{2}
\right).
\end{equation}
First of all we note that, since the solution of the associated impurity model
furnishes the spectral density, $\rho_{loc}$ from which the new hybridization
function is calculated, we never have to use the self energy in the solution
to the model.  It is not until one wants to calculate lattice dependent
response functions that the self energy has to be used.  We also see from the
equations above why the self energy becomes pathological.  Namely as long as
$|G_{loc}|$ is less than $1/t$ everything is fine but if this doesn't hold,
which happens at the quasiparticle peak at low temperatures, things go awry.
It is also easy to understand why the breakdown temperature decreases with
density and degeneracy.  As we mentioned earlier the total weight of the
spectral density is $1-n+n/N$ so as $n$ and $N$ increase this weight becomes
less and therefore the likelihood of crossing the bound on $|G_{loc}|$ is
diminished.

We should note that faulty numerics can be ruled out as a possible source for
this breakdown -- it is known from analytic results
\cite{Muller-Hartmann:1984,Parcollet:1998} that the solution of the NCA
equations does indeed exceed this bound for low temperatures.

We should note that other lattice place similar condition on the local Green's
function.  To see this let us compute the local Green's function as prescribed
by the DMFT:
\begin{eqnarray}
  \label{eq:densityBound}
  -\imag G_{loc}(\omega) &=& \imag\sum_k
   \frac{-1}{\omega+\mu-\epsilon_k-\Sigma(\omega)}
   \nonumber \\ &=& 
   \int_{-\infty}^{\infty}d\epsilon D(\epsilon)
   \imag\left(\frac{-1}{\omega+\mu-\epsilon-\Sigma(\omega)}\right)
   \nonumber \\ &\leq&
   D(\epsilon)_{max}\int_{-\infty}^{\infty}d\epsilon
   \imag\left(\frac{-1}{\omega+\mu-\epsilon-\Sigma(\omega)}\right)
   \nonumber \\ &=&
   \pi D(\epsilon)_{max}.
\end{eqnarray}
Here we have assumed that the self energy is properly retarded.  We see now
that the local spectral density $\rho(\omega) = -\imag G_{loc}(\omega)/\pi$
can never be greater than maximum value of the noninteracting density of
states, $D(\omega)$.  On the Bethe lattice this gives us the condition,
$\rho_{loc}(\omega) \leq 1/t$.  It would seem strange that the DMFT should
place such restrictions on the spectral density; we are after all trying to
capture the effects of strong correlations and thus one of the features that
we can expect to see narrow resonances of large magnitude in the density of
states.  

We will not dwell on this pathology here, it is quite sufficient for our
purposes to restrict the use of the NCA to higher temperatures and we will
employ other methods in the low temperature regime in the computation of
transport properties in the following chapters.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 






