\chapter{Transport in Real Multiband Systems: Orthogonal Case}
\label{chap:real-transport}
\marginlabel{chap:real-transport}

\section{Introduction}
In this chapter we will continue to focus on the computation of transport
properties of correlated materials but instead of restricting ourselves to
simple models like we did in previous chapters we will extend our discussion
to realistic models.  Thus we will be working with Hamiltonians obtained from
first principle bandstructure calculations (LMTO) and we will derive the
equations for the lowest order contributions to the transport parameters.  As
such the Hamiltonian will be given at each point in the Brillouin zone by a
matrix in the space of valence orbitals that describes both the energies and
the mixing of the orbitals in question.

As an example of such a Hamiltonian we can take for instance the Hamiltonian
for LaTiO.  In this material the bands that are closest to the Fermi surface
and participate in the transport are transition metal d-orbitals which
hybridize with oxygen p-orbitals that are slightly below the Fermi-surface.
The rest of the orbitals are inert and can be either folded down or to first
approximation, ignored.  We thus obtain an effective Hamiltonian of the form:
\begin{equation}
H_{0} =
\sum_{k\alpha\beta}\epsilon^{p}_{k\alpha\beta}p^{+}_{k\alpha}p_{k\beta}
+ 
\sum_{k\alpha\beta}\epsilon^{d}_{k\alpha\beta}d^{+}_{k\alpha}d_{k\beta}
+
\sum_{k\alpha\beta}\left(
V_{k\alpha\beta}p^{+}_{k\alpha}d_{k\beta}
+
V^{*}_{k\alpha\beta}d^{+}_{k\beta}p_{k\alpha}
\right).
\end{equation}
Here we have suppressed the spin index.  To this Hamiltonian we want to add an
on-site Hubbard type interaction for the d-electron.  Since this on-site
interaction has been taken into account in some average fashion in the LMTO
calculation some changes should be made also to the original Hamiltonian,
$H_{0}$.  We will however not worry about that for the moment.  Thus our final
Hamiltonian is:
\begin{equation}
H = H_{0}+{U \over 2}\sum_{i}\sum_{(\alpha,\sigma)\neq(\alpha',\sigma')}
d^{+}_{i\alpha\sigma}d_{i\alpha\sigma}
d^{+}_{i\alpha'\sigma'}d_{i\alpha'\sigma'}.
\end{equation}
To simplify the notation we will in the following write the LMTO Hamiltonian
as follows:
\begin{equation}
H_{0} = \sum_{kmn}\varepsilon_{kmn}c^{+}_{km}c_{kn}
\end{equation}
where $\varepsilon_{kmn}$ would equal $\epsilon^{p}_{kmn}$ and
$\epsilon^{d}_{kmn}$ in the diagonal p-orbital sector and d-orbital sector
respectively and $V_{kmn}$ in the off diagonal.

Now in order to compute the transport properties we also need the self energy
and that is where the DMFT comes in.  In recent years some efforts have been 
put into extending the DMFT to deal with real materials and the general DMFT
approach to real materials and clusters has recently been formalized and been
shown to be causally consistent \cite{Kotliar:2000}.

In what follows here we will simply assume that the self energy is available
as an output of a realistic DMFT calculation.

%---------------------------------------------------------------------------

\section{The currents and the transport equations}
\label{sec:trans-eq}
\marginlabel{sec:trans-eq} 
The classical approach to the computation of the transport properties in real
materials is based on the Boltzmann equation \cite{Ziman:Principles,Kim:1997}.
Here we however use the same approach as in the previous chapters, i.e. we use
the Kubo formalism and relate the transport parameters to retarded
current-current correlation functions.  Hence the first thing we must do is 
to obtain expressions for the currents.  In the multiorbital case the best
framework for deriving the expressions for the currents is by using gauge
invariance principles as explained by Moreno and Coleman \cite{Moreno:1996}.
Since the Hubbard type interaction is a density-density interaction it will
not contribute to the expression for currents which can be derived from a
gauge-principle as the particle current and the thermal current.  Thus in the
following we don't really have to worry about the interaction term except for
the fact that all correlation functions that we encounter have to be
calculated with the full Hamiltonian.  Hence the currents are given by: 
\begin{equation}
\label{eq:real-heatcur}
\vec{Q} = {1 \over 2}\sum_{kmn}\vec{v}_{kmn}
\left(\dot{c}^{+}_{km}c_{kn}-c^{+}_{km}\dot{c}_{kn}\right).
\end{equation}
where $\vec{v}_{kmn} = \vec{\nabla}_{k}\varepsilon_{kmn}$.
The particle current has a similar expression and is given by:
\begin{equation}
\label{eq:real-parcur}
\vec{j} = \sum_{kmn}\vec{v}_{kmn}
c^{+}_{km}c_{kn}.
\end{equation}
To find the transport coefficients we must calculate the following correlation
functions:
\begin{eqnarray}
Z_{jj}^{\nu\lambda}(\tau) &=& 
\langle j^{\nu}(\tau)j^{\lambda}(0)\rangle 
= \sum_{k}\sum_{mn}\sum_{pq}v^{\nu}_{kmn}v^{\lambda}_{kpq}
\Xi^{(1)}_{mnpq}(\tau,k) \\
Z_{jQ}^{\nu\lambda}(\tau) &=& 
\langle j^{\nu}(\tau)Q^{\lambda}(0)\rangle 
= \sum_{k}\sum_{mn}\sum_{pq}v^{\nu}_{kmn}v^{\lambda}_{kpq}
\Xi^{(2)}_{mnpq}(\tau,k) \\
Z_{QQ}^{\nu\lambda}(\tau) &=& 
\langle Q^{\nu}(\tau)Q^{\lambda}(0)\rangle 
= \sum_{k}\sum_{mn}\sum_{pq}v^{\nu}_{kmn}v^{\lambda}_{kpq}
\Xi^{(3)}_{mnpq}(\tau,k)
\end{eqnarray}
where we have defined the following expectation values:
\begin{gather}
\Xi^{(1)}_{mnpq}(\tau,k) = 
\langle c^{+}_{km}(\tau)c_{kn}(\tau)c^{+}_{k'p}c_{k'q} \rangle \\
\Xi^{(2)}_{mnpq}(\tau,k) = {1 \over 2}\langle c^{+}_{km}(\tau)c_{kn}(\tau)
\left(
\dot{c}^{+}_{k'p}c_{k'q} - c^{+}_{k'p}\dot{c}_{k'q}
\right) \rangle \\
\Xi^{(3)}_{mnpq}(\tau,k) = 
{1 \over 4}\langle \left(
\dot{c}^{+}_{km}(\tau)c_{kn}(\tau) - c^{+}_{km}(\tau)\dot{c}_{kn}(\tau)
\right)
\left(
\dot{c}^{+}_{k'p}c_{k'q} - c^{+}_{k'p}\dot{c}_{k'q}
\right) \rangle 
\end{gather}
These expectation values will be evaluated shortly but first we need
to establish some notation.
We define our Green's functions as follows $(\tau > 0)$:
\begin{equation}
G_{mn}(\tau,k) = -\langle c_{km}(\tau)c^{+}_{kn}\rangle 
= {-1 \over Z}Tr\left\{ 
e^{-\beta H}e^{\tau H}c_{km}e^{-\tau H}c^{+}_{kn}
\right\}
\end{equation}
The Fourier transform of this can be calculated and is given by:
\begin{equation}
G_{mn}(i\omega,k) = \left[ 
{1 \over i\omega+\mu-\hat{\varepsilon}_{k}-\hat{\Sigma}(i\omega,k)}
\right]_{mn}
\end{equation}
Note that here both $\hat{\varepsilon}_{k}$ and
$\hat{\Sigma}(i\omega,k)$ are matrices.  It is useful to have a
spectral representation of the Green's functions (see appendix
\ref{app:spectral}),
\begin{equation}
G_{mn}(i\omega,k) = \int_{-\infty}^{\infty}
{\rho_{kmn}(\epsilon) \over i\omega-\epsilon}d\epsilon
\end{equation}
where the spectral density $\hat{\rho}_{k}(\epsilon)$ is the discontinuity
of the Green's function at the branch cut on the real axis, i.e.
\begin{equation}
\hat{\rho}_{k}(\epsilon) = {-1 \over 2\pi i}\left(
\hat{G}^{ret}(\epsilon,k)-\hat{G}^{adv}(\epsilon,k)
\right)
= {-1 \over 2\pi i}\left(
\hat{G}^{ret}(\epsilon,k)-\left(\hat{G}^{ret}(\epsilon,k)\right)^{+}
\right).
\end{equation}
In terms of the spectral density we can write down an explicit formula
for the time dependent Green's function:
\begin{equation}
G_{mn}(\tau,k) = -\int_{-\infty}^{\infty}d\epsilon\rho_{kmn}(\epsilon)
e^{-\epsilon\tau}\left(
\theta(\tau)f(-\epsilon)-\theta(-\tau)f(\epsilon)
\right).
\end{equation}
We now look at the expectation values above and evaluate them using a diagram
expansion keeping only connected diagrams with no vertex corrections.  The
expectation value of time derivatives of the creation and annihilation
operators can be related to the time derivative of the Green's function as
discussed in Appendix \ref{app:spectral} and using those relations we get:
\begin{eqnarray}
\Xi^{(1)}_{mnpq}(\tau,k) &=& -G_{np}(\tau,k)G_{qm}(-\tau,k) 
\nonumber\\ &=&
\int_{-\infty}^{\infty}
d\epsilon d\epsilon'\rho_{knp}(\epsilon)\rho_{kqm}(\epsilon')
e^{(\epsilon'-\epsilon)\tau}
f(-\epsilon)f(\epsilon')
\\
\Xi^{(2)}_{mnpq}(\tau,k) &=& 
\dot{G}_{np}(\tau,k)G_{qm}(-\tau,k) 
- G_{np}(\tau,k)\dot{G}_{qm}(-\tau,k) 
\nonumber\\ &=&
\int_{-\infty}^{\infty}
d\epsilon d\epsilon'\rho_{knp}(\epsilon)\rho_{kqm}(\epsilon')
e^{(\epsilon'-\epsilon)\tau}
f(-\epsilon)f(\epsilon')\left({\epsilon+\epsilon' \over 2}\right)
\\
\Xi^{(3)}_{mnpq}(\tau,k) &=& 
\dot{G}_{np}(\tau,k)\dot{G}_{qm}(-\tau,k)
- G_{np}(\tau,k)\ddot{G}_{qm}(-\tau,k) 
\nonumber\\&-&
\ddot{G}_{np}(\tau,k)G_{qm}(-\tau,k)
+ \dot{G}_{np}(\tau,k)\dot{G}_{qm}(-\tau,k) 
\nonumber\\ &=&
\int_{-\infty}^{\infty}
d\epsilon d\epsilon'\rho_{knp}(\epsilon)\rho_{kqm}(\epsilon')
e^{(\epsilon'-\epsilon)\tau}
f(-\epsilon)f(\epsilon')\left({\epsilon+\epsilon' \over 2}\right)^{2}
\end{eqnarray}
The current-current correlation functions are then:
\begin{eqnarray}
Z^{\nu\lambda}_{jj}(\tau) &=& \sum_{k}\int_{-\infty}^{\infty}
d\epsilon d\epsilon'\phi^{\nu\lambda}(\epsilon,\epsilon')
e^{(\epsilon'-\epsilon)\tau} f(-\epsilon)f(\epsilon')
\\
Z^{\nu\lambda}_{jQ}(\tau) &=& \sum_{k}\int_{-\infty}^{\infty}
d\epsilon d\epsilon'\phi^{\nu\lambda}(\epsilon,\epsilon')
e^{(\epsilon'-\epsilon)\tau}
f(-\epsilon)f(\epsilon')\left({\epsilon+\epsilon' \over 2}\right)
\\
Z^{\nu\lambda}_{QQ}(\tau) &=& \sum_{k}\int_{-\infty}^{\infty}
d\epsilon d\epsilon'\phi^{\nu\lambda}(\epsilon,\epsilon')
e^{(\epsilon'-\epsilon)\tau}
f(-\epsilon)f(\epsilon')\left({\epsilon+\epsilon' \over 2}\right)^{2}.
\end{eqnarray}
Where we have defined the transport function
$\phi^{\nu\lambda}(\epsilon,\epsilon')$ as
\begin{equation}
\phi^{\nu\lambda}(\epsilon,\epsilon') 
\equiv \sum_{k}\trace\left\{
\hat{v}^{\nu}_{k}\hat{\rho}_{k}(\epsilon)
\hat{v}^{\lambda}_{k}\hat{\rho}_{k}(\epsilon')
\right\}. 
\end{equation}
Here Tr denotes trace in the space of orbitals at a given point in the
Brillouin zone.  We note that all the operators on the trace are
hermitian so that the trace is a real number.  It is now easy to take
the Fourier transform of the $Z$ functions and we obtain:
\begin{eqnarray}
Z^{\nu\lambda}_{jj}(i\nu_{n}) &=& \int_{-\infty}^{\infty}
d\epsilon d\epsilon'\phi^{\nu\lambda}(\epsilon,\epsilon')
{f(\epsilon)-f(\epsilon') \over i\nu_{n}+\epsilon'-\epsilon}
\\
Z^{\nu\lambda}_{jQ}(i\nu_{n}) &=& \int_{-\infty}^{\infty}
d\epsilon d\epsilon'\phi^{\nu\lambda}(\epsilon,\epsilon')
{f(\epsilon)-f(\epsilon') \over i\nu_{n}+\epsilon'-\epsilon}
\left({\epsilon+\epsilon' \over 2}\right)
\\
Z^{\nu\lambda}_{QQ}(i\nu_{n}) &=& \int_{-\infty}^{\infty}
d\epsilon d\epsilon'\phi^{\nu\lambda}(\epsilon,\epsilon')
{f(\epsilon)-f(\epsilon') \over i\nu_{n}+\epsilon'-\epsilon}
\left({\epsilon+\epsilon' \over 2}\right)^{2}
\end{eqnarray}
Finally we obtain the kinetic coefficients by looking at the imaginary
part of the retarded version of the $Z$ functions, dividing by
$\beta\omega$ and taking the limit $\omega\rightarrow 0$.  Thus we get:
\begin{eqnarray}
\label{eq:omegaZeroLimit}
\zeta^{\nu\lambda}_{jj} \equiv 
{1 \over \beta}\lim_{\omega\rightarrow 0}{1 \over \omega}\imag 
Z^{\nu\lambda}_{jj}(\omega+i0) &=& \pi\int_{-\infty}^{\infty}
d\epsilon\phi^{\nu\lambda}(\epsilon)
\left(-k_{b}T{df(\epsilon)\over d\epsilon}\right)
\nonumber\\
\zeta^{\nu\lambda}_{jQ} \equiv 
{1 \over \beta}\lim_{\omega\rightarrow 0}{1 \over \omega}\imag 
Z^{\nu\lambda}_{jQ}(\omega+i0) &=& \pi\int_{-\infty}^{\infty}
d\epsilon\phi^{\nu\lambda}(\epsilon)
\left(-k_{B}T{df(\epsilon)\over d\epsilon}\right)
\epsilon
\nonumber\\
\zeta^{\nu\lambda}_{QQ} \equiv 
{1 \over \beta}\lim_{\omega\rightarrow 0}{1 \over \omega}\imag 
Z^{\nu\lambda}_{QQ}(\omega+i0) &=& \pi\int_{-\infty}^{\infty}
d\epsilon\phi^{\nu\lambda}(\epsilon)
\left(-k_{B}T{df(\epsilon)\over d\epsilon}\right)
\epsilon^{2}
\end{eqnarray}
where for simplicity we denote $\phi^{\nu\lambda}(\epsilon,\epsilon)$
by $\phi^{\nu\lambda}(\epsilon)$, i.e.  
\begin{equation}
\phi^{\nu\lambda}(\epsilon) 
\equiv \sum_{k}\trace\left\{
\hat{v}^{\nu}_{k}\hat{\rho}_{k}(\epsilon)
\hat{v}^{\lambda}_{k}\hat{\rho}_{k}(\epsilon)
\right\}. 
\end{equation}
\ We should note here that if we make a momentum independent unitary
transformation in the space of orbitals, the transport function is left
invariant.  Hence if the self energy is momentum independent which is the case
in dynamical mean field theory we can always make \label{page:diagonalSigma} a
unitary transformation to a basis where at least the antihermitian damping
part of the self energy is diagonal.

%---------------------------------------------------------------------------

\section{Transport function}
\label{sec:transfunc}
\marginlabel{sec:transfunc}
It is useful for further numerical calculations to investigate the
transport function a little bit closer.  For this purpose let us look
at the energy and write
\begin{equation}
\hat{\varepsilon}_{k}+\hat{\Sigma}(z,k)
= \hat{U}(z,k)\hat{E}(z,k)\hat{U}^{-1}(z,k)
\end{equation}
where the matrix $\hat{E}_{mn} = \delta_{mn}E_{m}$ is a diagonal
matrix of eigenvalues and the columns of $\hat{U}$ contain the
corresponding eigenvectors.  Note that the eigenvalues are complex and
thus the eigenvectors are not necessarily orthogonal.  Then we can
write:
\begin{equation}
\hat{G}(z,k) = \hat{U}(z,k)\hat{D}(z,k)\hat{U}^{-1}(z,k)
\qquad\mbox{where}\qquad
\hat{D}_{mn} = {\delta_{mn} \over z+\mu-\hat{E}_{m}}.
\end{equation}
Suppressing the frequency and momentum dependence, assuming that
$\hat{G}$ is the retarded Green's function, using the cyclic property
of the trace and the fact that the velocity is hermitian we get:
\begin{eqnarray}
\trace\left\{
\hat{v}^{\nu}\hat{\rho}\hat{v}^{\lambda}\hat{\rho}
\right\}
&=&
{1 \over (2\pi i)^{2}}\trace\left\{
\hat{v}^{\nu}(\hat{G}-\hat{G}^{+})
\hat{v}^{\lambda}(\hat{G}-\hat{G}^{+})
\right\} \nonumber \\ &=& 
{-1 \over (2\pi)^{2}}\left(
\trace\left\{\hat{v}^{\nu}\hat{G}\hat{v}^{\lambda}\hat{G}\right\}-
\trace\left\{\hat{v}^{\nu}\hat{G}^{+}\hat{v}^{\lambda}\hat{G}\right\}
\right. \nonumber \\
&-&
\left.
\trace\left\{\hat{v}^{\nu}\hat{G}\hat{v}^{\lambda}\hat{G}^{+}\right\}+
\trace\left\{\hat{v}^{\nu}\hat{G}^{+}\hat{v}^{\lambda}\hat{G}^{+}\right\}
\right)
\nonumber \\ &=& 
{-1 \over (2\pi)^{2}}\left(
\trace\left\{\hat{v}^{\nu}\hat{G}\hat{v}^{\lambda}\hat{G}\right\}-
\trace\left\{\hat{v}^{\nu}\hat{G}^{+}\hat{v}^{\lambda}\hat{G}\right\}
\right. \nonumber \\
&-&
\left.
\trace\left\{\hat{G}\hat{v}^{\lambda}\hat{G}^{+}\hat{v}^{\nu}\right\}+
\trace\left\{\hat{G}^{+}\hat{v}^{\lambda}\hat{G}^{+}\hat{v}^{\nu}\right\}
\right)
\nonumber \\ &=& 
{-1 \over (2\pi)^{2}}\left(
\trace\left\{\hat{v}^{\nu}\hat{G}\hat{v}^{\lambda}\hat{G}+
\left[\hat{v}^{\nu}\hat{G}\hat{v}^{\lambda}\hat{G}\right]^{+}
\right\}
\right. \nonumber \\
&-&
\left.
\trace\left\{\hat{v}^{\nu}\hat{G}^{+}\hat{v}^{\lambda}\hat{G}+
\hat{v}^{\lambda}\hat{G}^{+}\hat{v}^{\nu}\hat{G}
\right\}
\right)
\nonumber \\ &=& 
{-2 \over (2\pi)^{2}}\real\left(
\trace\left\{L^{\nu\lambda}\right\} - 
\trace\left\{M^{\nu\lambda}\right\}
\right)
\end{eqnarray}
Where we have defined:
\begin{eqnarray}
L^{\nu\lambda} &=& \hat{v}^{\nu}\hat{G}\hat{v}^{\lambda}\hat{G}
= \hat{v}^{\nu}\hat{U}\hat{D}\hat{U}^{-1}
  \hat{v}^{\lambda}\hat{U}\hat{D}\hat{U}^{-1}
\nonumber \\
M^{\nu\lambda} &=& {1 \over 2}\left[
\hat{v}^{\nu}\hat{G}^{+}\hat{v}^{\lambda}\hat{G}+
\hat{v}^{\lambda}\hat{G}^{+}\hat{v}^{\nu}\hat{G}
\right]
\nonumber \\
&=& {1 \over 2}\left[
\hat{v}^{\nu}(\hat{U}^{-1})^{+}\hat{D}^{+}\hat{U}^{+}
\hat{v}^{\lambda}\hat{U}\hat{D}\hat{U}^{-1}+
\hat{v}^{\lambda}(\hat{U}^{-1})^{+}\hat{D}^{+}\hat{U}^{+}
\hat{v}^{\nu}\hat{U}\hat{D}\hat{U}^{-1}
\right]
\end{eqnarray}
Note that taking the real part of the $M$ matrix is really redundant
since it is already real.
Let us now define the following three operators:
\begin{equation}
\hat{r}^{\nu} = \hat{U}^{-1}\hat{v}^{\nu}\hat{U},
\qquad
\hat{s}^{\nu} = \hat{U}^{-1}\hat{v}^{\nu}\left(\hat{U}^{-1}\right)^{+}
\qquad\mbox{and}\qquad
\hat{t}^{\nu}_{3} = \hat{U}^{+}\hat{v}^{\nu}\hat{U}.
\end{equation}
Then our trace can be written as
\begin{equation}
\trace\left\{
\hat{v}^{\nu}\hat{\rho}\hat{v}^{\lambda}\hat{\rho}
\right\} = {-1 \over 2\pi^{2}}\real\trace\left\{
\hat{r}^{\nu}\hat{D}\hat{r}^{\lambda}\hat{D}-
{1 \over 2}\left(
\hat{s}^{\nu}\left(\hat{D}\right)^{+}\hat{t}^{\lambda}\hat{D}+
\hat{s}^{\lambda}\left(\hat{D}\right)^{+}\hat{t}^{\nu}\hat{D}
\right)\right\}
\end{equation}
We can thus write the transport function on the form
\begin{eqnarray}
\label{eq:transportfunctiondiagonalform}
\phi^{\nu\lambda}(\epsilon) &=& {-1 \over 2\pi^{2}}\sum_{k}
\real\left(\trace\left\{
\hat{r}^{\nu}_{k}(\epsilon)\hat{D}_{k}(\epsilon)
\hat{r}^{\lambda}_{k}(\epsilon)\hat{D}_{k}(\epsilon)\right\}\right)
\nonumber\\
&-&
{1 \over 4\pi^{2}}\sum_{k}
\left(\trace\left\{
\hat{s}^{\nu}_{k}(\epsilon)\left(\hat{D}_{k}(\epsilon)\right)^{+}
\hat{t}^{\lambda}_{k}(\epsilon)\hat{D}_{k}(\epsilon)+
\hat{s}^{\lambda}_{k}(\epsilon)\left(\hat{D}_{k}(\epsilon)\right)^{+}
\hat{t}^{\nu}_{k}(\epsilon)\hat{D}_{k}(\epsilon)\right\}
\right)
\nonumber \\ &=&
{-1 \over 2\pi^{2}}\sum_{k}\real\sum_{pq}\left(
\hat{r}^{\nu}_{kpq}(\epsilon)D_{kq}(\epsilon)
\hat{r}^{\lambda}_{kqp}(\epsilon)D_{kp}(\epsilon)\right)
\nonumber \\
&-&
{1 \over 4\pi^{2}}\sum_{k}\sum_{pq}\left(
\hat{s}^{\nu}_{kpq}(\epsilon)D_{kq}^{*}(\epsilon)
\hat{t}^{\lambda}_{kqp}(\epsilon)D_{kp}(\epsilon)
+\hat{s}^{\lambda}_{kpq}(\epsilon)D_{kq}^{*}(\epsilon)
\hat{t}^{\nu}_{kqp}(\epsilon)D_{kp}(\epsilon)
\right)
\nonumber \\ &=&
{-1 \over 2\pi^{2}}\sum_{k}\real\sum_{pq}\left(
{\alpha^{\nu\lambda}_{kpq}(\epsilon) \over 
(\epsilon+\mu-E_{kp}(\epsilon))(\epsilon+\mu-E_{kq}(\epsilon))}
\right)
\nonumber \\ &+&
{1 \over 2\pi^{2}}\sum_{k}\sum_{pq}\left(
{\beta^{\nu\lambda}_{kpq}(\epsilon) \over 
(\epsilon+\mu-E_{kp}(\epsilon))(\epsilon+\mu-E^{*}_{kq}(\epsilon))}
\right)
\end{eqnarray}
where we have defined:
\begin{equation}
\alpha^{\nu\lambda}_{kpq}(\epsilon) = 
\hat{r}^{\nu}_{kpq}(\epsilon)\hat{r}^{\lambda}_{kqp}(\epsilon)
\qquad\mbox{and}\qquad
\beta^{\nu\lambda}_{kpq}(\epsilon) = {1 \over 2}\left(
\hat{s}^{\nu}_{kpq}(\epsilon)\hat{t}^{\lambda}_{kqp}(\epsilon)+
\hat{s}^{\lambda}_{kpq}(\epsilon)\hat{t}^{\nu}_{kqp}(\epsilon)
\right)
\end{equation}

%---------------------------------------------------------------------------

\section{Small scattering limit}
\label{sec:smallscattering}
It is interesting to consider here the limit of the transport function when
the scattering tends to zero.  It is of course quite clear that the transport
function diverges in this limit but it is nevertheless worthwhile to pursue
this limit and see where the strongest contribution comes from.  So we will
assume that the self energy is given by 
\begin{equation}
  \label{eq:smallscatteringSigma}
  \Sigma(\omega) = \gamma\Sigma_0(\omega)
\end{equation}
where $\gamma$ is a small parameter which approaches zero and $\Sigma_0$ is
diagonal\footnote{Note that there is no loss of generality in assuming that
  the self energy is diagonal, see page \pageref{page:diagonalSigma}.}.  In
the limit that $\gamma$ goes to zero the eigenvalues of the full Hamiltonian,
$H_k+\Sigma$ will approach the eigenvalues of the bare Hamiltonian, $H_k$, and
in particular the imaginary part of the eigenvalues will tend to zero linearly
to first order in $\gamma$.  The eigenvectors will also approach the
eigenvectors of the bare Hamiltonian.  Looking at the diagonalized expression
for the transport function (\ref{eq:transportfunctiondiagonalform}) we see that
the expected divergence will come from the energy denominators and therefore
is is enough to evaluate the numerators to zeroth order in $\gamma$.  To
zeroth order the eigenvector matrix $U$ is simply equal to the bare
eigenvector matrix $U_B$ which is unitary and thus we find that in this
approximation the matrices $r^{x}$, $s^x$ and $t^x$ are all equal and are
simply:
\begin{equation}
  r^x = s^x = t^x = U_B^+ v^x U_B.
\end{equation}
Before we take this any further let us look again at the expression for the
velocity.  Using the diagonal representation of the bare Hamiltonian, $H_k =
U_B E_B U_B^+$ where $E_B$ is the diagonal matrix of eigenvalues of the bare
Hamiltonian we can write the velocity as
\begin{equation}
  \vec{v} = \grad{k}H_k = 
U_B\left\{\vec{v}_B-[E_B,\vec{K}]\right\}U_B^+
\end{equation}
where we define the bare velocity and the matrix $\vec{K}$ as 
\begin{equation}
  \vec{v}_B = \grad{k}E_b
  \qquad\mbox{and}\qquad
  \vec{K} = U_B^+\grad{k}U_B.
\end{equation}
Hence we see that 
\begin{equation}
  r^x = v_B^x-[E_B,K^x].
\end{equation}
We note that the matric $E_B$ is diagonal so we can write the matrix elements
of $r^x$ as
\begin{equation}
   r^x_{pq} = \delta_{pq}v^x_p-K^x_{pq}(E_p-E_q).
\end{equation}
The matrix elements $\alpha$ and $\beta$ from equation
(\ref{eq:transportfunctiondiagonalform}) are now given by: 
\begin{eqnarray}
  \alpha^{xx}_{pq} &=& \beta^{xx}_{pq} = r^x_{pq}r^x_{qp}
  = \delta_{pq}(v_p)^2+K^x_{pq}(E_p-E_q)K^x_{qp}(E_q-E_p)
  \nonumber \\
  &=& \delta_{pq}(v_p)^2+|K^x_{pq}|^2(E_p-E_q)^2.
\end{eqnarray}
where we have used the fact that $K^x$ is antihermitian.  The equation for the
transport function simplifies now tremendously and becomes:
\begin{equation}
  \phi^{xx} = \frac{1}{\pi^2}\sum_{kpq}
  \left[\delta_{pq}(v_p)^2+|K^x_{pq}|^2(E_p-E_q)^2\right]
  \imag\left(\frac{1}{\epsilon+\mu-E_{kp}}\right)
  \imag\left(\frac{1}{\epsilon+\mu-E_{kq}}\right).
\end{equation}
Here we have still kept the complex energies in the denominator of the
spectral function, but as the parameter $\gamma$ tends to zero these spectral
functions become delta functions.  The off diagonal terms, $p\neq q$, in this
expression will vanish on account of the delta function if the energies in the
denominators are not the same, and if they are degenerate then the energy
difference term in the matrix element will see to that they vanish.  Hence we
are left with only the diagonal term term involving the bare velocity and this
term of course is formally divergent.  This result is however quite often used
in practice and then one takes one of the spectral functions as a delta
function and the other gives the lifetime in the particular orbital.  We can
think of the lifetime of an orbital as given by:
\begin{equation}
  \tau_{kp} = \frac{1}{2\pi|\imag E_{kp}|}
\end{equation}
and even though we have assumed that this quantity is infinite in our
derivation the final result is quite often written as:
\begin{equation}
  \phi^{xx} = \sum_{kp}
  (v_{kp})^2\tau_{kp}\delta(\epsilon+\mu-E_{B,kp}).
\end{equation}
This expression is nothing but the Boltzmann equation expression for the
transport function \cite{Ziman:Principles,Kim:1997} with the exception
that our definition of the lifetime does not include vertex corrections
\cite{Mahan}.  This was of course clear from the start since the calculation 
here only includes the lowest order diagram contributing to the transport
coefficients.  This is consistent with the dynamical mean field theory
\cite{Georges:1996} where vertex corrections drop out in the infinite
dimensional limit.  Undoubtedly vertex correction will play a role for real
three dimensional materials but the numerical calculations would quickly
become quite unwieldy.

Since the computation of the small scattering limit is much simpler than
computing the full transport function it is useful to have this computation as
a reference for the full transport function computation.  By scaling it by the
parameter $\gamma$ we obtain a finite result which can be compared to a
similarly scaled full computation.  Thus we define 
\begin{equation}
  \label{eq:BareTransportFunction}
  \phi^{xx}_{Bare} = \sum_{kp}
  (v_{kp})^2\tau_{kp}\delta(\epsilon+\mu-E_{B,kp}).
\end{equation}
where the lifetime is computed by adding a very small scattering to the
Hamiltonian and diagonalizing it and then we used the equation
\begin{equation}
  \label{eq:BareTransportLifetime}
  \tau_{kp} = \frac{\gamma}{2\pi|\imag E_{kp}|}.
\end{equation}.

%---------------------------------------------------------------------------

\section{The diagonal case}
\label{sec:DiagonalCase}
Before we go into the details on the computation of the transport
function in the general case it is worth looking at the case where
both the LMTO Hamiltonian and the self energy are diagonal.  In this
case as in the small scattering limit the computation of the transport
function simplifies tremendously as we will see.  In this case the
velocity will be real and of course diagonal and the spectral density
will also be diagonal and is simply be given by $\frac{-1}{\pi}\imag
\hat{G}$.  Thus we get:
\begin{eqnarray}
\phi^{\nu\mu}(\omega) &=& 
        \sum{k}\trace{\hat{v}^{\nu}\hat{\rho}_{k}(\omega)
        \hat{v}^{\mu}\hat{\rho}_{k}(\omega)}
\nonumber \\
&=& \sum_{kn}
        \left(\frac{\partial \epsilon_{kn}}{\partial k_{\nu}}\right)
        \left(\frac{\partial \epsilon_{kn}}{\partial k_{\mu}}\right)    
        \rho^{2}_{kn}(\omega)
\\
&=& \frac{1}{\pi^2}\sum_{kn}
        \left(\frac{\partial v_{kn}}{\partial k_{\nu}}\right)
        \left(\frac{\partial v_{kn}}{\partial k_{\mu}}\right)   
        \left(\imag\left[
        \frac{1}{\omega+\mu-\epsilon_{kn}-\Sigma_{n}(\omega)}
        \right]\right)^2
\nonumber\\
&=& \frac{1}{\pi^2}\sum_{n}\int_{-\infty}^{\infty}d\epsilon
        \Phi^{\nu\mu}_{n}(\epsilon)
        \left(\imag\left[
        \frac{1}{\omega+\mu-\epsilon-\Sigma_{n}(\omega)}
        \right]\right)^2
\nonumber
\end{eqnarray}
where $\Phi^{\nu\mu}_{n}(\epsilon)$ is the so called transport density
(c.f. equation \ref{eq:phifunction}) defined by:
\begin{eqnarray}
\Phi^{\nu\mu}_{n}(\epsilon) &=& \sum_{k}
        \left(\frac{\partial v_{kn}}{\partial k_{\nu}}\right)
        \left(\frac{\partial v_{kn}}{\partial k_{\mu}}\right)   
        \delta(\epsilon-\epsilon_{kn})
        \nonumber \\
        &=& \frac{-1}{\pi}\imag\sum_{k}
        \left(\frac{\partial v_{kn}}{\partial k_{\nu}}\right)
        \left(\frac{\partial v_{kn}}{\partial k_{\mu}}\right)   
        \frac{1}{\epsilon-\epsilon_{kn}+i0}.
\end{eqnarray}
This function is much simpler than the original transport function and
can easily be computed by the tetrahedron method (c.f. appendix
\ref{app:atm} and ref.~\cite{Lambin:1984}) and in some cases can be done
almost exactly using fast Fourier transforms.

%---------------------------------------------------------------------------

\section{A detailed look at the integral}
After this brief aside on the small scattering limit and the simple
case of a diagonal Hamiltonian we now focus our attention on the task
at hand, i.e. computing the integral:
\begin{equation}
  \label{eq:bzintegral}
  \phi^{xx}
  = \sum_{k \in BZ} \trace\left\{v^{x}_{k}\rho_{k}v^{x}_{k}\rho_{k}\right\}
\end{equation}
where BZ denotes the Brillouin zone.  In the numerical computation the
Brillouin zone is the volume inside the cube spanned by the three basis
vectors, $\vec{k}_1$, $\vec{k}_2$ and $\vec{k}_3$ of the reciprocal lattice.
We construct from these vectors the so called minicell vectors, $\vec{p}_i$ by:
\begin{equation}
  \label{eq:minicell}
  \vec{p}_i = \frac{1}{N_i}\vec{k}_i
\end{equation}
where $N_1$ is an integer.  A grid is then constructed in the Brillouin zone
from the minicell vectors as follows:
\begin{equation}
  \label{eq:kmesh}
  \vec{k}_{i_1,i_2,i_3} = i_1 \vec{p}_1 + i_2 \vec{p}_2 + i_3 \vec{p}_3
  \qquad i_m = 0,\ldots,N_m-1.
\end{equation}
In this manner the Brillouin zone is split up into a collection of minicells
that can also be indexed by $i_1$, $i_2$ and $i_3$ and whose corners are at
the points $\vec{k}_{i_1\pm 1,i_2\pm 1,i_3\pm 1}$.  Note that here we use the
periodicity of the crystal to identify $\vec{k}_{N_1,i_2,i_3}$ with
$\vec{k}_{0,i_2,i_3}$ and similarly for the other directions.  Each of these
minicells can now be split up into six tetrahedra (see figure
\ref{fig:CubeSplit}) and those we will use for the computation of the
transport function.  There are four different ways to do this tetrahedral
splitting and usually it is done such that the common diagonal is the shortest
diagonal in the minicell.  Once the splitting is decided the $m-th$ corner of
a given tetrahedron, $t$, $t = 1,\ldots,6$ in the minicell indexed by $i_1$,
$i_2$ and $i_3$ is given by
\begin{equation}
  \label{eq:tetcorners}
  (j_1,j_2,j_3) = (i_1,i_2,i_3) + (kcut(m,t,1),kcut(m,t,2),kcut(m,t,3))
\end{equation}
where $kcut(m,t,n)$ takes the values zero or one.  If we denote this
tetrahedron by $\Delta(i_1,i_2,i_3,t)$ the transport function integral can now
be written as follow:
\begin{equation}
  \label{eq:tetintegral}
  \phi^{xx}
  = \sum_{i_1=0}^{N_1}\sum_{i_2=0}^{N_2}\sum_{i_3=0}^{N_3}\sum_{t=1}^{6}
  \sum_{k \in \Delta(i_1,i_2,i_3,t)} 
  \trace\left\{v^{x}_{k}\rho_{k}v^{x}_{k}\rho_{k}\right\}
\end{equation}

In the LMTO computations of the band structure only the irreducible Brillouin
zone (IBZ) is used and we get data for the Hamiltonian only at the so called
irreducible mesh points.

The points in the irreducible zone are found in the following way.  We
enumerate all the mesh points in some fashion and then take the first mesh
point to be irreducible.  The rest of the irreducible mesh points are then
found recursively by applying all the group operations of the system to each
of the mesh points, if a given operation takes it into an irreducible mesh
point (up to a translation by a reciprocal lattice vector) then it is
reducible, otherwise we label it as irreducible.  

In this fashion we end up with a total of $N_{kp}$ irreducible k-points and
for each mesh point in the whole Brillouin zone we get a pointer
$kirrptr(i_1,i2,i3)$ whose value is the index of the corresponding irreducible
k-point and another pointer $gropptr(i_1,i_2,i_3)$ whose value is the number
of the group operation, $\gamma$, that takes the irreducible point to the mesh
point, possibly modulo a reciprocal lattice vector.

Let us now consider a given tetrahedron whose corners are $\vec{k}_i$, $i =
1,\ldots,4$.  To these four points correspond four or less irreducible mesh
points $\vec{k}^{irr}_{i}$ and we have $\vec{k} =
R^{\gamma_i}\vec{k}^{irr}_{i}$ where $R^{\gamma_i}$ is a rotation matrix
corresponding to the group operation $\gamma_i$.   The four irreducible k
points are said to be an irreducible tetrahedron.  

The problem is now how to relate the spectral densities and velocities at the
corners in the irreducible tetrahedron to the corresponding quantities at the
corners of the real tetrahedron.  The Hamiltonian transforms in the following
manner \cite{Savrasov:notes}:
\begin{eqnarray}
  \label{eq:rotatedhamiltonian}
  H_{l' m'\tau',lm\tau}(R^{\gamma}\vec{k})
  &=& e^{i\vec{k}\cdot(\vec{R}^{g}_{\tau'}-\vec{R}^{g}_{\tau})}
  \sum_{m_{1}m_{2}}
  U^{l'}_{m'm_{1}}(\gamma)
  H_{l'm_{1}\hat{g}^{-1}(\tau),lm_{2}\hat{g}^{-1}(\tau')}(\vec{k})
  U^{l*}_{mm_{2}}(\gamma)
  \nonumber \\ &=&
   e^{i\vec{k}\cdot(\vec{R}^{g}_{\tau'}-\vec{R}^{g}_{\tau})}
  \left(
    U(\gamma)
    H_{\hat{g}^{-1}(\tau),\hat{g}^{-1}(\tau')}(\vec{k})
    U^{+}(\gamma)
  \right)_{l'm',lm}
\end{eqnarray}       
Here $\hat{g}^{-1}$ denotes the operation reciprocal (i.e. in real space) to
$\gamma$ plus a possible shift (by a direct lattice vector) for non-symmorphic
group.  Hence $\hat{g}^{-1}(\tau)$ is the atom corresponding to $\tau$ after
the rotation and if the operation does not directly rotate the atom to $\tau$
into $\hat{g}^{-1}(\tau)$ then it must be shifted back to the original cluster
by a direct lattice vector given by $R^{g}_{\tau}$.  This shift is responsible
for the extra exponential factor in the transformation formula.

We note that the spectral density, $\rho_{k}$, being essentially the inverse
of the Hamiltonian, transforms in the same way as the Hamiltonian but the
velocity does not due to the extra exponential factor in the transformation.
The velocity is given by the gradient of the Hamiltonian and therefore the
velocity matrix in the $\tau',\tau$ sector is given by
\begin{eqnarray}
  \label{eq:rotatedvelocity}
  v^{x}_{\tau',\tau}(R^{\gamma}\vec{k}) &=&
  \frac{\partial}{\partial(R^{\gamma}\vec{k})_x}
  H_{\tau',\tau}(R^{\gamma}\vec{k})
  =
  \frac{\partial}{\partial k_j}
  H_{\tau',\tau}(R^{\gamma}\vec{k})
  \frac{\partial k_j}{\partial(R^{\gamma}\vec{k})_x}
  \nonumber \\ &=&
  e^{i\vec{k}\cdot(\vec{R}^{g}_{\tau'}-\vec{R}^{g}_{\tau})}U(\gamma)
  \left(
    i[(\vec{R}^{g}_{\tau'})^j-(\vec{R}^{g}_{\tau})^j]
    H_{\hat{g}^{-1}(\tau),\hat{g}^{-1}(\tau')}(\vec{k})
  \right)
  U^{+}(\gamma)(R^{\gamma})^{-1}_{jx}
  \nonumber \\ &+&
  e^{i\vec{k}\cdot(\vec{R}^{g}_{\tau'}-\vec{R}^{g}_{\tau})}U(\gamma)
  \left(
    v^{j}_{\hat{g}^{-1}(\tau),\hat{g}^{-1}(\tau')}(\vec{k})
  \right)
  U^{+}(\gamma)(R^{\gamma})^{-1}_{jx}.
\end{eqnarray}
Let us now define the semirotated velocity as:
\begin{equation}
  \label{eq:semirotatedvelocity}
  v^{\gamma,x}_{\hat{g}^{-1}(\tau),\hat{g}^{-1}(\tau')}(\vec{k})
  = \left(
    i[(\vec{R}^{g}_{\tau'})^j-(\vec{R}^{g}_{\tau})^j]
    H_{\hat{g}^{-1}(\tau),\hat{g}^{-1}(\tau')}(\vec{k})+
    v^{j}_{\hat{g}^{-1}(\tau),\hat{g}^{-1}(\tau')}(\vec{k})
    \right)(R^{\gamma})^{-1}_{jx}.
\end{equation}
Then we can write the rotated velocity in terms of the semirotated velocity 
as:
\begin{equation}
  \label{eq:rotatedvelocity2}
  v^{x}_{\tau',\tau}(R^{\gamma}\vec{k}) =
  e^{i\vec{k}\cdot(\vec{R}^{g}_{\tau'}-\vec{R}^{g}_{\tau})}U(\gamma)
  v^{\gamma,x}_{\hat{g}^{-1}(\tau),\hat{g}^{-1}(\tau')}(\vec{k})
  U^{+}(\gamma)
\end{equation}
We can now proceed to evaluate the transport function trace at each of the
corners in the tetrahedron and to make clear the atomic dependence we write
explicitly the atomic labels (repeated indices are summed over):
\begin{eqnarray}
  \label{eq:cornertrace}
  \trace\left\{
    v^{x}(\vec{k}_i)\rho(\vec{k}_i)v^{x}(\vec{k}_i)\rho(\vec{k}_i)
  \right\} &=&
  \trace\left\{
    v^{x}_{\tau_1\tau_2}(\vec{k}_i)\rho_{\tau_2\tau_3}(\vec{k}_i)
    v^{x}_{\tau_3\tau_4}(\vec{k}_i)\rho_{\tau_4\tau_1}(\vec{k}_i)
  \right\}
  \nonumber \\ &=&
  \trace\left\{
    v^{\gamma,x}_{\hat{g}^{-1}(\tau_1)\hat{g}^{-1}(\tau_2)}(\vec{k}^{irr}_i)
    \rho_{\hat{g}^{-1}(\tau_2)\hat{g}^{-1}(\tau_3)}(\vec{k}^{irr}_i)
    \right.   \nonumber \\ && \hspace{2ex}\times\left.
    v^{\gamma,x}_{\hat{g}^{-1}(\tau_3)\hat{g}^{-1}(\tau_4)}(\vec{k}^{irr}_i)
    \rho_{\hat{g}^{-1}(\tau_4)\hat{g}^{-1}(\tau_1)}(\vec{k}^{irr}_i)
  \right\}
  \nonumber \\ &=&
  \trace\left\{
    v^{\gamma,x}(\vec{k}^{irr}_i)\rho(\vec{k}^{irr}_i)
    v^{\gamma,x}(\vec{k}^{irr}_i)\rho(\vec{k}^{irr}_i)
  \right\}. 
\end{eqnarray}
Here we used the cyclicity of the trace to get rid of the Wigner matrices, the
exponential factors cancel automatically and in the end we remove the group
operation dependence on the dummy atomic labels.  Thus we have managed to
express the trace at each corner in terms of quantities at the corresponding
irreducible k point but the velocity, as we see from equation
(\ref{eq:semirotatedvelocity}), acquires an extra term from the Hamiltonian.
We should remark here that since the Wigner matrices and the exponential
factors have been canceled independently at each corner of the tetrahedron
the individual factors of the trace at each corner do not match with the
corresponding factors from the other corners. Thus it would make no sense for
instance to linearize each of these factors inside the tetrahedron and try to
do the resulting integral analytically.  It is the whole expression that is
continuous from one corner to the next, not the individual factors.  The
eigenvalues of the Hamiltonian on the other hand are the same at the real
tetrahedron corner and the irreducible corners.  Hence if we express the trace
in terms of the eigenvalues of the Hamiltonian as we did in the last section
the $\alpha$ and $\beta$ terms, evaluated at the irreducible corners but with
the semirotated velocities, are continuous from one corner to the next and we
are justified in using any approximation we want for those terms individually
from the terms involving the eigenvalues.  This is quite useful numerically
since this means that we do not have to find the Hamiltonian all throughout
the Brillouin zone and diagonalize it separately at each mesh point -
diagonalization on the irreducible mesh is enough.  We also do not need to
compute the velocity at each mesh point, we just have to find it on the
irreducible mesh and then the semirotated velocity can be generated as needed
from the velocity and the Hamiltonian.

Hence we must compute $6\times N_1\times N_2\times N_3$ tetrahedral
contributions to the transport function.  The computation of the velocity and
the diagonalization of the Hamiltonian is just done at the irreducible mesh
points which are much fewer.  Since each mesh point occurs in a total of six
tetrahedra the $\alpha$ and $\beta$ terms are best precomputed before we start 
summing up all the individual contributions.
  
%---------------------------------------------------------------------------

\section{Numerical evaluation of $\phi^{\nu\lambda}$}
\label{sec:comp-phi}
\marginlabel{sec:comp-phi} 

In this section we discuss the numerical computation of the
tetrahedral contributions to the transport function.  The calculation
of the contribution from each tetrahedron involves the numerical
evaluation of a triple integral of the form
\begin{equation}
T_{SS} = \sum_{k\in \Delta}{F(k) \over (z-E_{k\alpha})(z-E_{k\gamma})}
\qquad\mbox{or}\qquad
T_{OS} = \sum_{k\in \Delta}{F(k) \over (z-E_{k\alpha})(z-E_{k\gamma}^*)}
\end{equation}
where $\Delta$ denotes the tetrahedron.  The subscripts $SS$ and $OS$ refer to
the sign of the imaginary part of the energies in the denominator.  In
$T_{SS}$ they both have the same sign and in $T_{OS}$ they have the opposite
sign.  It is useful when doing integrals like these to symmetrize the function
$F$ using the same symmetry operations as we use to collapse the whole
Brillouin zone down to the irreducible part.  Even though we are doing the
integral over the whole Brillouin zone and therefore this symmetrization is
not really necessary the results will probably be better and less dependent on
the choice of tetrahedron splitting of the Brillouin zone minicells.  This
symmetrization is easily done and is formally given by the equation,
\begin{equation}
F_S(k) = {1 \over N_S}\sum_{\sigma}{\cal O}_{\sigma}F(k)
\end{equation}
where $N_S$ is the number of symmetry operations and $\sigma$ denotes
a symmetry operation with the corresponding operator being
${\cal O}_{\sigma}$.  In our case here this symmetrization is done right from
the start by symmetrizing the combination of the velocity factors in the
transport function.  So even though we are really computing $\sigma^{xx}$ we
would in a cubic system compute $(\sigma^{xx}+\sigma^{yy}+\sigma^{zz})/3$.

The functions $E_{k\alpha}$, $E_{k\gamma}$ and $F(k)$ are of course
only known at the corners of the tetrahedron.  Let us denote their
values at those corners by $\alpha_i$, $\gamma_i$ and $f_i$, $i =
1,2,3,4$ and let $k_i$ denote the coordinate of corner $i$.  Then a
vector, $k$ inside the tetrahedron can be written on the form,
\begin{equation}
k = k_1+a(k_2-k_1)+b(k_3-k_1)+c(k_4-k_1)
\qquad
a,b,c \geq 0,
\quad\mbox{and}\quad
a+b+c \leq 1.
\end{equation}
Using linear interpolation on the functions $F$ and $E$ inside the
tetrahedron we can write them in terms of $a,b,c$ as follows,
\begin{gather}
F(a,b,c) = f_1+a(f_2-f_1)+b(f_3-f_1)+c(f_4-f_1)
\\
E_{\alpha}(a,b,c) = 
\alpha_1+a(\alpha_2-\alpha_1)+b(\alpha_3-\alpha_1)+c(\alpha_4-\alpha_1)
\end{gather}
It is useful to redefine the energies so as to get rid of the
external frequency in the integral.
\begin{equation}
\alpha_i-z \rightarrow \alpha_i, \qquad
\end{equation}
and we note that since $z$ is assumed real this does not change the imaginary
part of the energies.  
%In the numerical calculations we do in fact add a small
%imaginary part to $z$ to increase stability but its sign is chosen to be
%opposite to the sign of the imaginary part of $\alpha$ and $\gamma$.  This
%means of course that $z$ can contain imaginary parts of different signs 
%in the same expression.  
Now since we use a linear interpolant for the function $F$
the integral will be linear in the values $f_i$ so the answer will be on the
form:
\begin{equation}
w_1f_1+w_2f_2+w_3f_3+w_4f_4
\end{equation}
where the symbol $w_i$ is called a weight factor for corner $i$.  

There are four terms that contribute to the transport function:
\begin{enumerate}
\item SSDD: The diagonal part ($\alpha = \gamma$) in the SS sector
\item SSOD: The off-diagonal part ($\alpha \neq \gamma$) in the SS sector
\item OSDD: The diagonal part ($\alpha = \gamma$) in the OS sector
\item OSOD: The off-diagonal part ($\alpha \neq \gamma$) in the OS sector
\end{enumerate}
and we will use different methods depending on which contribution we
are evaluating.  The SSDD integral can be done analytically and that
calculation is carried out in appendix \ref{app:ssdd}.  The other
contributions can to some extent be carried out analytically,
depending on the degeneracy of the corners on the tetrahedron, but the
lion's share of the tetrahedra are in the interior of the irreducible
Brillouin zone and have four different irreducible corners and that
case has not yielded an analytical solution.We will therefore explore
alternative approximations.
\begin{enumerate}
\item OSDD: This case can be written in the following form 
\begin{equation}
T 
= \sum_{k\in \Delta}{F(k) \over |z-E_{k\alpha}|^2}
= \sum_{k\in \Delta}{F(k) \over \varepsilon_k^2+\gamma_k^2}
= \sum_{k\in \Delta}{F(k) \over \gamma_k}
\imag\left({1 \over \varepsilon_k-i\gamma_k}\right)
\end{equation}          
Here we have split the denominator $z-E_{k\alpha}$ into its real and imaginary
parts.  The imaginary part, $\gamma$ is due to the self energy and since the
self energy in the DMFT is momentum independent it is reasonable to expect
that $\gamma$ will not depend very strongly on momentum.  Hence we propose to
approximate the factor $F(k)/\gamma$ by a linear function and then we can
simply calculate the integral using the analytical tetrahedron method of
Lambin and Vigneron \cite{Lambin:1984}.  This way of splitting up the
integrand also makes the dependence on the scattering rate evident by
factoring it out and in the limit of small scattering rates it is clear that
the resulting transport function will scale as the inverse of the scattering
rate.
\item SSOD and OSOD: The method used for these two cases are the same.  We
  want to try to reduce this case down to the analytical tetrahedron method
  and we will do that by identifying the part of the denominator that has
  greater modulus and then write 
\begin{equation}
T 
= \sum_{k\in \Delta}{F(k) \over (z-E_{k\alpha_1})(z-E_{k\alpha_2})}
= \sum_{k\in \Delta}\left({F(k) \over (z-E_{k\alpha})_{large}}\right)
{1 \over (z-E_{k\alpha})_{small}}.
\end{equation}
The denominator in the parenthesis has a large modulus and thus the
expression in the parentheses will change relatively little in the
tetrahedron.  We will therefore approximate the expression in the
parentheses by a linear function and in that way we have managed to
reduce the integral down to a form suitable for the analytical
tetrahedron method as we had hoped for.
\end{enumerate}


These methods seem to give quite good results as we can see in
\article\ref{sec:analytical_transport_test} here below, except for the fact
that the OSDD integral sometimes needs to be treated a little differently
depending on the behavior of the scattering.  The reason for this is that in
uncorrelated orbitals the scattering is only due to the hybridization with the
correlated orbitals and at high symmetry points in the Brillouin zone some of
the uncorrelated orbitals will not aquire any scattering in this manner..
This leads to a pole on the real axis in the spectral density and the
transport function will diverge at that point.  When this occurs far away from
the Fermi-surface the philosophy is that we can simply ignore it knowing that
when calculating the transport parameters the transport function will be
weighted by a function which is exponentially decreasing away from the
Fermi-surface.  If, on the other hand this occurs close to the Fermi-surface
we must go back and rethink our calculations because the transport function
will diverge and there is nothing we can do better numerically to prevent
that.

Going back to the calculation of the OSDD case we see by inspecting
the formula and our strategy for evaluating it that we are liable to
run into problems with divergence even though we are far from any real
poles.  The reason is that when we write:
\begin{equation}
{1 \over \epsilon_{k}^{2}+\gamma_{k}^{2}} = 
{1 \over \gamma_{k}}\imag\left({1 \over \epsilon_{k}-i\gamma_{k}}\right)
\end{equation}
we see that there is a delicate cancellation between the two factors
in the limit that $\gamma$ goes to zero.  If we however treat these
two factors differently in the numerical calculation as we have
proposed to do, this cancellation will not occur and we end up with
diverging results where there was no divergence in the first place.

This problem was studied in detail with several different strategies for 
overcoming it such as writing the integral on the following form:
\begin{equation}
T 
= \sum_{k\in \Delta}{F(k) \over |z-E_{k\alpha}|^2}
= \sum_{k\in \Delta}\left({F(k) \over (z-E_{k\alpha})^{*}}\right)
{1 \over (z-E_{k\alpha})}
\end{equation}  
and linearize completely the factor in parenthesis and then use the analytical
tetrahedron method.  Using this method proved problematic, since the switching
between methods tended to introduce discontinuities in the results.  It seems
that the most consistent way to deal with this problem is to add a small
artificial broadening to the uncorrelated orbitals but one should still keep
in mind that any orbital that is close to the Fermi level will affect the
transport properties and must be regarded as correlated in this context and
should be assigned some reasonable self energy.

%---------------------------------------------------------------------------

\section{Numerical calculation of the velocity: Method I}
\label{sec:cmp1.vel}
\marginlabel{sec:cmp1.vel}
The velocity matrix for the k-points in the irreducible Brillouin zone
has to be calculated in order to be able to calculate the transport
function.  Thus we need to calculate the derivatives,
\begin{equation}
v_{kpq}^i = {\partial \over \partial k^i}\epsilon_{kpq}.
\end{equation}
This is done numerically by using the formula,
\begin{equation}
f'(0) \simeq {1 \over 12}(f(-2)-8f(-1)+8f(1)-f(2)).
\end{equation}
Using this formula means of course that we need to know what the
neighbors of a given k-point in the irreducible zone are.  These
neighbors are possibly not inside the irreducible zone so we need to
know to what points in the irreducible zone they correspond.  This
information can be obtained from the German LMTO program although a
little modification is necessary in order to do so.  The lattice
points in the Brillouin zone are given by the equation:
\begin{equation}
\label{eq:BZmeshpoints}
\vec{k} = \sum_{j=1}^{3}p_j{\vec{b}_j \over N_j},
\qquad 
p_j = 0,\ldots,N_j-1
\end{equation}
where the vectors $\vec{b}_j$ are the generators of the reciprocal
lattice.  In terms of these we get,
\begin{equation}
v_{kpq}^i = {\partial \over \partial k^i}\epsilon_{kpq}
= \sum_{j=1}^{3}{\partial p^j \over \partial k^i}
{\partial \over \partial p^j}\epsilon_{kpq}.
\end{equation}
Using the generators $\vec{a}_i$ of the Bravais lattice and the fact
that $\vec{a}_i\cdot\vec{b}_j = 2\pi\delta_{ij}$ we have,
\begin{equation}
p_j = {N_j \over 2\pi}\vec{a}_j\cdot\vec{k}
\end{equation}
and thus we get
\begin{equation}
v_{kpq}^i = \sum_{j=1}^{3}{N_j \over 2\pi}a_j^i
{\partial \over \partial p^j}\epsilon_{kpq}.
\end{equation}
We now define the matrix
\begin{equation}
\label{eq:Amatrix}
A_{ij} = N_j a_j^i
\end{equation}
and thus we can write
\begin{equation}
v_{kpq}^i = {1 \over 2\pi}\sum_{j=1}^{3}A_{ij} 
{\partial \over \partial p^j}\epsilon_{kpq}.
\end{equation}
The energy is now defined on a grid of integers in terms of the
$p$-variables and the derivative is found by using the formula quoted
above.

\section{Numerical calculation of the velocity: Method II}
\label{sec:cmp2.vel}
\marginlabel{sec:cmp2.vel}
There is an alternative way to calculate the velocity which we should
mention also.  As discussed previously the Hamiltonian is approximated
linearly within each tetrahedron and thus if we wanted to be faithful
to that strategy the velocity should be taken as constant within each
tetrahedron.  Thus let us consider a tetrahedron with corners whose
coordinates are given by $k_1$, $k_2$, $k_3$ and $k_4$ with $H_1$,
$H_2$, $H_3$ and $H_4$ being the corresponding values of the
Hamiltonian.  Thus inside the tetrahedron the we can represent the
momentum and approximate the Hamiltonian as follows:
\begin{gather}
k = k_1+a\Delta k_{21}+b\Delta k_{31}+c\Delta k_{41}
\nonumber \\
H = H_1+a\Delta H_{21}+b\Delta H_{31}+c\Delta H_{41}
\nonumber \\
a,b,c \geq 0 \quad\mbox{and}\quad a+b+c \leq 1
\end{gather}
where $\Delta k_{ji} = k_j-k_i$ and $\Delta H_{ji} = H_j-H_i$.
Let us now define the following vectors:
\begin{eqnarray}
\Delta p_{21} &=& \frac{\Delta k_{31}\times\Delta k_{41})}
{\Delta k_{21}\cdot(\Delta k_{31}\times\Delta k_{41})},
\\
\Delta p_{31} &=& \frac{\Delta k_{41}\times\Delta k_{21})}
{\Delta k_{21}\cdot(\Delta k_{31}\times\Delta k_{41})}
\\
\Delta p_{41} &=& \frac{\Delta k_{21}\times\Delta k_{31})}
{\Delta k_{21}\cdot(\Delta k_{31}\times\Delta k_{41})}.
\end{eqnarray}
Now we see that we can write
\begin{equation}
a = (k-k_1)\cdot\Delta p_{21} 
\qquad\mbox{and}\qquad
\vec{\nabla}\!_{k}a = \Delta p_{21} 
\qquad\mbox{etc.}
\end{equation}
The gradient of the Hamiltonian inside the tetrahedron is therefore given 
by:
\begin{eqnarray}
\vec{\nabla}\!_{k}H 
&=& 
\left(\frac{\partial}{\partial a}H\right)\vec{\nabla}\!_{k}a+
\left(\frac{\partial}{\partial b}H\right)\vec{\nabla}\!_{k}b+
\left(\frac{\partial}{\partial c}H\right)\vec{\nabla}\!_{k}c
\nonumber \\ &=& 
\Delta H_{21}\Delta p_{21}+
\Delta H_{31}\Delta p_{31}+
\Delta H_{41}\Delta p_{41}
\end{eqnarray}
We should note here that this method, although closer in spirit to the
approximations that we are making is not used in our calculations for the
simple reason that one would then have to compute the velocity for every
single tetrahedron in the Brillouin zone.  With the other method we only need
to compute the full velocity in the irreducible part of the Brillouin zone and
then the semirotated velocity in the rest of the Brillouin zone. This would
also mean that the cancellation that is described in connection with equation
(\ref{eq:cornertrace}) would not be possible and we would thus have to carry
out a rotation of the spectral density for each different value of the self
energy whereas with the other method this is unnecessary.

%---------------------------------------------------------------------------

\section{Testing the methods}
\label{sec:analytical_transport_test}
Most of the formulas that we have derived in the previous pages have
been tested with some simple data using Mathematica.  It is however
also useful to be able to test the whole formalism on a case that can
be calculated using some other means.  For that purpose we have
created a simple Hamiltonian based on the tight-binding Hamiltonian in
a three dimensional cubic lattice.  It is given by the following
formula:
\begin{eqnarray}
H &=& \sum_{l=1,2}\sum_{k}c_{kl}^{+}c_{kl}(D_{low}\epsilon_{k}+E_{low})
+ \sum_{k}d_{k}^{+}d_{k}(D_{high}\epsilon_{k}+E_{high})
\nonumber \\ &+&
\sum_{l=1,2}\sum_{k}(Vc_{kl}^{+}d_{k}+V^{*}d_{k}^{+}c_{kl})
\end{eqnarray}
with the dispersion relation given by:
\begin{equation} 
\epsilon_k =  {1 \over 3}\sum_{i=1}^{3}\cos(ak_i).
\end{equation}
This Hamiltonian describes a doubly degenerate low lying itinerant band of
half-width $D_{low}$ coupled to a localized band of half-width
$D_{high}<D_{low}$ which we assume lies close to the Fermi-surface.  To this
we add a frequency and momentum independent scattering, $\Gamma$, in the
localized level.  This model is a simplified version of what we expect to be
dealing with for real Hamiltonians, i.e. a low lying itinerant band
hybridizing with a localized correlated orbital close to the Fermi-surface.
%----------------------------
\begin{figure}[htb]
\centering
\mbox{\epsfig{file=Figures-4/test-dos.ps,angle=-90,width=0.90\linewidth}}
\caption{The total density of states of the model Hamiltonian. The 
parameters of the Hamiltonian used are: $E_{low} = -7.0$, $E_{high} =
0$, $\Gamma = 0.5$ and $V = 1.0$.} 
\label{fig:phi-test-dos}
\end{figure}
%----------------------------
In figure (\ref{fig:phi-test-dos}) we have displayed the DOS for this model
corresponding to the parameters, $E_{low} = -7.0$, $E_{high} = 0$, $\Gamma =
0.5$, $V = 1.0$, $D_{low} = 1.0$ and $D_{high} =0.1$.  We can see the low
lying band with fairly well defined band edges and then the higher correlated
band with longer tails Due to the hybridization the bands are pushed slightly
apart and the correlated level transfers a little bit of it's weight to the
lower level and vice versa.  Note however that for this value of scattering
the width of the correlated peak is entirely governed by the magnitude of the
scattering and is not due to the inherent width of the band or the
hybridization with the itinerant band.

\subsection{Analytical calculation of the transport function}
\label{subsec:analytic}
We see that the velocity corresponding to the test Hamiltonian is
diagonal and is given by,
\begin{equation}
\hat{v}_{k}^{x} = \hat{V}
{\partial \epsilon_k \over \partial k^x} 
\end{equation}
where $\hat{V}$ is diagonal and equal to the $D_{low}$ in the uncorrelated
sector and $D_{high}$ in the correlated sector. The transport function is
therefore given by,
\begin{eqnarray}
\phi^{xx}(\omega) &=& 
\sum_{k}\trace\{
\hat{v}^{x}_{k}\hat{\rho}_{k}(\omega)
\hat{v}^{x}_{k}\hat{\rho}_{k}(\omega)
\}
\nonumber\\ &=& 
\int d\epsilon \trace\{
\hat{V}\hat{\rho}(\epsilon,\omega)
\hat{V}\hat{\rho}(\epsilon,\omega)
\}\left[
\sum_{k}\left({\partial \epsilon_k \over \partial k^x}\right)^2
\delta(\epsilon-\epsilon_k)
\right]
\end{eqnarray}
On the three dimensional cubic lattice the function in brackets can be
calculated accurately using Fourier transforms, and then the transport
function can be calculated from there using numerical integration over
$\epsilon$.  In order to be able to test each part of the integration
code is proves useful to expand the trace in terms of the Green's
functions as we did earlier in the chapter and then the contribution
from each term can be compared to the numerical code.  We note that in
order to make exact the correspondence between the analytic
calculation and the numerical calculation we need to do same kind of
diagonalization as we did in the numerical case.
%----------------------------
\begin{figure}[htb]
\centering
\begin{tabular}{c@{\quad}c}
\mbox{\epsfig{file=Figures-4/test-phi.ps,angle=-90,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-4/test-l2log.ps,angle=-90,width=0.45\linewidth}} 
\end{tabular}
\caption{The transport function, $\phi$, calculated using the exact 
  solution and the Brillouin zone integration with 8 and 16 points per axis.
  The parameters of the Hamiltonian used are: $E_{low} = -7.0$, $E_{high} =
  0$, $\Gamma = 0.5$, $V = 1.0$, $D_{low} = 1.0$ and $D_{high} =0.1$.  Also
  shown is the relative L$_2$ deviation,
  $||\phi-\phi_{exact}||_2/||\phi_{exact}||_2$ computed over this interval.}
\label{fig:phi-test1}
\end{figure}
%---------------------------------------------------------------------------
\begin{figure}[htb]
\centering
\begin{tabular}{c@{\quad}c}
\mbox{\epsfig{file=Figures-4/test-ssdd.ps,angle=-90,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-4/test-ssod.ps,angle=-90,width=0.45\linewidth}} \\
\mbox{\epsfig{file=Figures-4/test-osdd.ps,angle=-90,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-4/test-osod.ps,angle=-90,width=0.45\linewidth}}
\end{tabular}
\caption{The four contributions to the transport function, $\phi$, 
calculated using the exact solution and the Brillouin zone integration
with 8 and 16 points per axis.  The parameters of the Hamiltonian
used are the same as in figure (\ref{fig:phi-test1}).}
\label{fig:phi-test2}
\end{figure}
%---------------------------------------------------------------------------
In figure (\ref{fig:phi-test1}) we have displayed the transport function for
the test Hamiltonian in the region around the position of the correlated
level.  We do not compute the transport function around the position of the
uncorrelated level since we did not assign any scattering to that orbital.  As
we can see the Brillouin zone integration method agrees very nicely with the
exact calculation and hardly any difference can be seen when using integration
mesh corresponding to 16 points per axis.  We have also displayed the relative
L$_2$ error as a function of the number of points per axis and as we can see
the error decreases here as the inverse fourth power of the number of points
per axis.

In figure (\ref{fig:phi-test2}) where each of the four contributions to the
transport function are displayed we can see that the convergence is equally
good for each contribution.  We also not that the contributions from the SS
sector are almost mirror images of the contributions from the OS sector and 
therefore the final result depends on delicate cancellation between
contributions from theses sectors.
%---------------------------------------------------------------------------
\begin{figure}[htb]
\centering
\begin{tabular}{c@{\quad}c}
\mbox{\epsfig{file=Figures-4/test-gamma1.ps,angle=-90,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-4/test-gamma2.ps,angle=-90,width=0.45\linewidth}} \\
\end{tabular}
\caption{The scaled transport functions, $\Gamma\phi$, calculated with 16 PPA.
  In the left figure the parameters of the Hamiltonian are the same as in
  figure (\ref{fig:phi-test1}) but in the left figure we have $D_{high} =
  2.0$.}
\label{fig:phi-test3}
\end{figure}
%---------------------------------------------------------------------------

The transport function diverges of course in the limit of zero scattering and
we expect that it will go as $1/\Gamma$ for small values of $\Gamma$. Thus if
we scale the transport function with the scattering then the scaled transport
function should be approximately independent of the scattering, in fact one
would expect that the result should not be very much different from the bare
transport function computed with a constant lifetime equal to $1/2\pi$.  In
figures \ref{fig:phi-test3} we have displayed the scaled transport function
using several values of $\Gamma$ and two different bandwidths for the
correlated orbital.  As we can see the scaling in the center of the broad band
is quite good whereas the scaling for the narrow band is much worse.  The
reason why we see deviations from the scaling behavior is of course due to
the broadening effect of the scattering.  In the center of a broad featureless
band the scattering does not alter the density of states significantly and in
that case one can expect that the scaling behavior would be manifest.
However in narrow bands and close to band edges the scattering dramatically
alters the density of states, creating long tails that extend outside the band
and there the transport function will not exhibit scaling behavior.  In fact
outside the bands we should expect exactly the opposite scaling -- the
transport function will decrease with the value of the scattering.

%---------------------------------------------------------------------------
\begin{figure}[htb]
\centering
\begin{tabular}{c@{\quad}c}
\mbox{\epsfig{file=Figures-4/test-bare1.ps,angle=-90,width=0.45\linewidth}} &
\mbox{\epsfig{file=Figures-4/test-bare2.ps,angle=-90,width=0.45\linewidth}} \\
\end{tabular}
\caption{The scaled transport functions, $\Gamma\phi$, calculated with 16 PPA.
  compared with the bare transport function using a lifetime of $1/2\pi$. 
  The parameters of the Hamiltonian are the same as used in figure
  \ref{fig:phi-test3}.}
\label{fig:phi-test4}
\end{figure}
%---------------------------------------------------------------------------

Finally in figure \ref{fig:phi-test4} we have displayed the scaled transport 
function along with the bare transport function.  As expected in light of 
the previous arguments we see that for the broad band the scaled transport 
function agrees fairly well with the bare one but for the narrow band the
quantitative agreement is much worse.

%---------------------------------------------------------------------------

\section{Afterthoughts: Further refinements and comments}
\label{sec:afterthoughts}
As we noted earlier in \article\ref{sec:comp-phi} the integration weights for
the transport function cannot be computed exactly within the linear
approximation except for a few degenerate cases.  Nevertheless we were able to
obtain fairly good estimates of the integration weights by making further
approximations.  However it is clear that in order to obtain accurate
numerical estimates of the transport function we have to use quite a fine grid
in the Brillouin zone.  It is interesting to investigate whether it would be
possible to generate this grid refinement dynamically and maybe develop an
adaptive approach to computing the transport function.

Here we will explore one such approach which is based on iterative subdivision
of the tetrahedra.  In this approach we will take the original tetrahedron and
subdivide it into eight smaller tetrahedra, see figure \ref{fig:subdivisions}.

\begin{figure}[t]
\centering
\epsfig{file=Figures-4/subdivisions.ps,angle=-90,width=\halfpage}
\caption{The subdivisions of the tetrahedron described by the matrices
  $A_{i}$.  The figure shows three of the smaller subtetrahedra along with one
  of the three choices for the common diagonal of the remaining five
  subtetrahedra.}
\label{fig:subdivisions}
\end{figure}                                                

The integral that we wish to compute can be represented as follows:
\begin{equation}
  F(\vec{\alpha},\vec{\epsilon},\vec{\gamma}) = 
  6\int_{\Delta}dadbdc
  \frac{\alpha_{1}+a\alpha_{21}+b\alpha_{31}+c\alpha_{41}}
  {(\epsilon_{1}+a\epsilon_{21}a+b\epsilon_{31}+c\epsilon_{41})
   (\gamma_{1}+a\gamma_{21}+b\gamma_{31}+c\gamma_{41})}
\end{equation}
where $\vec{\alpha} = (\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4})$ denotes
the value of $\alpha$ at the corners of the tetrahedron and $\alpha_{ij} =
\alpha_{i}-\alpha_{j}$ with the same conventions for $\epsilon$ and $\gamma$.
The first subdivision of this integral is now done in the following manner:
\begin{equation}
    F(\vec{\alpha},\vec{\epsilon},\vec{\gamma}) = \frac{1}{8}\sum_{i=1}^{8}
    F(A_{i}\vec{\alpha},A_{i}\vec{\epsilon},A_{i}\vec{\gamma})
\end{equation}
where the eight matrices $A_{i}$ describe the tetrahedron splitting.
Similarly the second subdivision can be written as:
\begin{equation}
    F(\vec{\alpha},\vec{\epsilon},\vec{\gamma}) = 
    \frac{1}{8^{2}}\sum_{j=1}^{8}\sum_{i=1}^{8}
    F(A_{j}A_{i}\vec{\alpha},A_{j}A_{i}\vec{\epsilon},A_{j}A_{i}\vec{\gamma}).
\end{equation}
The matrices $A_{i}$ can be chosen in three different ways with one choice
being:
\begin{eqnarray}
A_{1} = \frac{1}{2}
\begin{pmatrix}
  1  & 1 & 0 & 0  \\
  0  & 1 & 1 & 0  \\
  0  & 1 & 0 & 1  \\
  0  & 2 & 0 & 0  \\
\end{pmatrix}    
\quad
&&
A_{2} = \frac{1}{2}
\begin{pmatrix}
  1  & 0 & 1 & 0  \\
  0  & 1 & 1 & 0  \\
  0  & 0 & 1 & 1  \\
  0  & 0 & 2 & 0  \\
\end{pmatrix}
\qquad
A_{3} = \frac{1}{2}
\begin{pmatrix}
  1  & 0 & 0 & 1  \\
  0  & 1 & 0 & 1  \\
  0  & 0 & 1 & 1  \\
  0  & 0 & 0 & 1  \\
\end{pmatrix}
\nonumber\\
A_{4} = \frac{1}{2}
\begin{pmatrix}
  0  & 1 & 0 & 1  \\
  1  & 0 & 1 & 0  \\
  1  & 1 & 0 & 0  \\
  2  & 0 & 0 & 0  \\
\end{pmatrix}    
\quad
&&
A_{5} = \frac{1}{2}
\begin{pmatrix}
  0  & 1 & 0 & 1  \\
  1  & 0 & 1 & 0  \\
  1  & 0 & 0 & 1  \\
  2  & 0 & 0 & 0  \\
\end{pmatrix}
\qquad
A_{6} = \frac{1}{2}
\begin{pmatrix}
  0  & 1 & 0 & 1  \\
  1  & 0 & 1 & 0  \\
  1  & 0 & 0 & 1  \\
  0  & 0 & 1 & 1  \\
\end{pmatrix}
\nonumber\\
A_{7} = \frac{1}{2}
\begin{pmatrix}
  0  & 1 & 0 & 1  \\
  1  & 0 & 1 & 0  \\
  1  & 1 & 0 & 0  \\
  0  & 1 & 1 & 0  \\
\end{pmatrix}    
\quad
&&
A_{8} = \frac{1}{2}
\begin{pmatrix}
  0  & 1 & 0 & 1  \\
  1  & 0 & 1 & 0  \\
  0  & 1 & 1 & 0  \\
  0  & 0 & 1 & 1  \\
\end{pmatrix}.
\end{eqnarray}
It is clear that further subdivision are easily generated and moreover one is
free to use different levels of subdivision in different regions in the
original tetrahedron.  This tetrahedral subdivision has not been explored
numerically in any detail so far but it might prove to be useful in some cases
and possible one could develop some simple criteria to determine if
subdivision is beneficial.

Another issue that deserves some further thought is eigenvalue degeneracy and
its effect on the calculation of Brillouin zone integrals.  In our calculation
the eigenvalues are assigned a band index depending on the value of their real
parts.  In cases where there are no band crossings this is the correct way to
go but when bands do cross one can end up with some nasty kinks in the energy
bands.  This can lead to errors in the evaluation of Brillouin zone integrals
and the need to use very fine grids for computations.  In case of complete
band degeneracy at some points in the Brillouin zone we also run the risk of
obtaining band eigenvectors that are not continuous at the degeneracy points,
the reason being that there is no unique way of choosing the eigenvectors in
the multidimensional eigenspace of the degenerate eigenvalue.

We will see an example of this kind of degeneracy later in this thesis in
chapter \ref{chap:latio3_revisited} where we consider a simple three
dimensional t$_{2g}$ model for LaTiO$_{3}$.  In this case there are lines in
the Brillouin zone where we get double degeneracy of eigenvalues.  In this
model if we start out the calculation of the transport function in the cubic
harmonic basis in which the t$_{2g}$ band are diagonal the Hamiltonian matrix
along with the self energy is already diagonal and we are not likely to run
into numerical problems with discontinuities of the eigenvectors.  However,
since our computer program so far only works in the spherical harmonic basis
there is nothing that will a priory secure the continuity of the eigenvectors.
The Hamiltonian along with the self energy will have to be diagonalized at
each k-point; at the nondegenerate points we will obtain the correct cubic
harmonic eigenvectors but at the degenerate points we will in general obtain
eigenvectors that are some linear combination of the cubic harmonics.  We
believe that this kind of discontinuity is likely to lead to errors or noise
in the evaluation of the transport function and since these discontinuities
occur along lines and not just at single points in the Brillouin zone the
errors can be quite severe.  We hypothesize that calculating the transport
function from the start in the cubic harmonic basis will lead to more accuracy
and less noise in the numerical results.

Based on this observation we believe that in general the orbital basis that is
used to compute the transport function should be chosen such that inherent
band degeneracies have already been diagonalized.  This means that if there
exists a momentum independent transformation (such as from spherical to cubic
harmonics) that will diagonalize the LMTO Hamiltonian, one should work in that
basis.


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


