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

\section{Introduction}

In this chapter we will extend the discussion of chapter
\ref{chap:real-transport} to Hamiltonians defined in nonorthogonal bases.  The
reason why this is important is due to the fact that quite frequently the
basis of wavefunctions used in the solution of the electronic structure
problems is composed of localized atomic-like wavefunctions.  These
wavefunctions are by choice not orthogonal and although one could in principle
construct an orthogonal basis by using the Gram-Schmidt method one would quite
likely end up with a basis that is delocalized and the strong correlations
that were local in the orthogonal basis become non-local and much harder to
deal with.

The best way to start familiarizing oneself with the formalism of band
structure calculations is to consider the problem in real space, where
we construct localized wavefunctions $\chi_{l}(r)$ for each atom (here
l is a multi-index representing angular momentum, spin and individual
atoms within the unit cell) and consider the matrix elements of the
Hamiltonian in the basis of these wavefunctions and their translations
throughout the crystal lattice.  Thus the wavefunctions basis consists
of the $\{\chi_{lR}\}$ where $\chi_{lR}(r) = \chi_{l}(r-R)$ where $R$
is a lattice translation vector. If we denote the states of these
wavefunctions by $|\chi_{lR}\rangle$ then we look at the matrix
elements of the Hamiltonian and self energy:
\begin{equation}
H_{ll'}(R-R') = \langle\chi_{lR}|H|\chi_{l'R'}\rangle
\qquad\mbox{and}\qquad
\Sigma_{ll'}(R-R') = \langle\chi_{lR}|\Sigma|\chi_{l'R'}\rangle.
\end{equation}
In addition we have a matrix of overlaps due to the nonorthogonality
of the basis:
\begin{equation}
O_{ll'}(R-R') = \langle\chi_{lR}|\chi_{l'R'}\rangle.
\end{equation}
As we have indicated the matrix elements depend only on the relative
distance $t = R-R'$ and thus we can write our matrix elements in terms
of $t$.  These matrices elements are given to us in the Fourier
transformed form, 
\begin{equation}
H_{ll'}(k) = \sum_{t}e^{ikt}H_{ll'}(t)
\end{equation}
etc.  Similarly we can go back:
\begin{equation}
H_{ll'}(t) = \sum_{k}e^{-ikt}H_{ll'}(k).
\end{equation}
In our calculations we assume that the correlations are local and lead
to a local self energy within the dynamical mean field theory and thus
$\Sigma_{ll'}(t) = \Sigma_{ll'}\delta_{t,0}$.  So let us assume that
the self energy is given by a local operator which is diagonal in the
space of orbitals, and simply consist of multiplying the state by a 
given complex number, $\Gamma_{l}$.  Then we have 
\begin{equation}
\Sigma_{ll'}(t) = \Gamma_{l}\delta_{ll'}\delta_{t,0}
                  \langle\chi_{l0}|\chi_{l0}\rangle.
                = \Gamma_{l}\delta_{ll'}\delta_{t,0}O_{ll}(0).
\end{equation}
When we Fourier transform we end up with a momentum independent self energy
\begin{equation}
\Sigma_{ll'}(k) = \Gamma_{l}\delta_{ll'}O_{ll}(0) = 
                  \Gamma_{l}\delta_{ll'}\sum_{k}O_{ll}(k).
\end{equation}
This equation is important because it tells us how to scale the self
energy when it is added to the Hamiltonian.  We can in fact also scale
both the overlap matrix and the Hamiltonian in such a way that the
diagonal on-site matrix elements of the overlap matrix are equal to 1.0.
This we can obtain by defining:
\begin{equation}
O_{l} = \sum_{k}O_{ll}(k)
\end{equation}
and then scaling the matrix elements such that
\begin{equation}
H_{ll'}(k) \longrightarrow {H_{ll'}(k) \over \sqrt{O_{l}O_{l'}}}
\qquad\mbox{and}\qquad
O_{ll'}(k) \longrightarrow {O_{ll'}(k) \over \sqrt{O_{l}O_{l'}}}.
\end{equation}
This scaling has of course no effect on the actual calculations but
makes the connection between the matrix elements of the Hamiltonian
and the actual energy eigenvalues clearer.

\section{Many Body Theory in Nonorthogonal Bases}
\label{sec:no_mbt}
Our starting point here is a representation of the kinetic term of the
Hamiltonian in an orthogonal basis, $\{|i\rangle\}$ and we assume that
this basis is related to the non-orthogonal basis,
$\{|\alpha\rangle\}$ by the transformation matrix:
\begin{equation}
|i\rangle = \sum_{\alpha}|\alpha\rangle S_{\alpha i}
\qquad\mbox{and}\qquad
\langle i|      = \sum_{\alpha}\langle\alpha|S^{*}_{\alpha i}
                = \sum_{\alpha}S^{+}_{i\alpha}\langle\alpha|
\end{equation}
The Hamiltonian is now given by:
\begin{equation}
H       = \sum_{ij}\langle i|H|j\rangle c^{+}_{i}c_{j}
        = \sum_{ij\alpha\beta}S^{+}_{i\alpha}
           \langle\alpha|H|\beta\rangle S_{\beta j} c^{+}_{i}c_{j}
        = \sum_{\alpha\beta}H_{\alpha\beta}c^{+}_{\alpha}c_{\beta}.
\end{equation}
The last term in the equation above is a requirement that we place on
the creation and destruction operators in the non-orthogonal basis and
thus we find that:
\begin{equation}
c^{+}_{\alpha} = \sum_{i}c^{+}_{i}S^{+}_{i\alpha}
\qquad\mbox{and}\qquad
c_{\alpha} = \sum_{j}S_{\alpha j}c_{j}.
\end{equation}
The non-orthogonality of the basis is encoded in the overlap matrix,
$O_{\alpha\beta} = \langle\alpha|\beta\rangle$ and this matrix can be
related to the transformation matrix, $S$ in the following manner
\begin{equation}
\delta_{ij}     =  \langle i|j\rangle 
                = \sum_{\alpha\beta}
        S^{+}_{i\alpha}\langle\alpha|\beta\rangle S_{\beta j}
                = \sum_{\alpha\beta}S^{+}_{i\alpha}O_{\alpha\beta}S_{\beta j}
\end{equation}
Therefore we see that the overlap matrix is given by:
\begin{equation}
O = (SS^{+})^{-1}.
\end{equation}
We should note here that the creation operator $c^{+}_{\alpha}$ does
not create a particle in the state $|\alpha\rangle$ when acting on 
the vacuum, since as we see
\begin{equation}
c^{+}_{\alpha}|0\rangle 
= \sum_{i}c^{+}_{i}S^{+}_{i\alpha}|0\rangle
= \sum_{i}|i\rangle S^{+}_{i\alpha}
= \sum_{i\beta}|\beta\rangle S_{\beta i} S^{+}_{i\alpha}
= \sum_{\beta}|\beta\rangle O^{-1}_{\beta\alpha}.
\end{equation}
It is however worth noting that this state has unit overlap with the
state $|\alpha\rangle$ and zero overlap with all of the other
basis-states!
The commutation relationship of these operators are the same as for
regular Fermi operators except that we get:
\begin{equation}
\{c^{+}_{\alpha},c_{\beta}\} = 
\sum_{ij}S_{\beta j}\{c^{+}_{i},c_{j}\}S^{+}_{i\alpha} 
= S_{\beta i}S^{+}_{i\alpha} = O^{-1}_{\beta\alpha}
\end{equation}
Let us finally obtain the expression for the Green's function in the non
orthogonal basis:
\begin{equation}
G_{\alpha\beta}(\tau) = -\langle c_{\alpha}(\tau)c^{+}_{\beta}(0)\rangle.  
\end{equation}
The easiest way to calculate this Green's function is by looking at
the Lagrangian for the system in the orthogonal basis and then simply
transform it into the non-orthogonal one. We have (summation over
repeated indices implied):
\begin{eqnarray}
{\cal L} &=& c^{+}_{i}{\partial \over \partial\tau}c_{i}
         - c^{+}_{i}H_{ij}c_{j}
        \nonumber \\
         &=& S^{-1}_{i\beta }
           c^{+}_{\alpha}{\partial \over \partial\tau}c_{\beta}
           (S^{+})^{-1}_{\alpha i}
         - c^{+}_{\alpha}H_{\alpha\beta}c_{\beta}
        \nonumber \\
         &=& c^{+}_{\alpha}O_{\alpha\beta}
         {\partial \over \partial\tau}c_{\beta}
         - c^{+}_{\alpha}H_{\alpha\beta}c_{\beta}.
\end{eqnarray}
The free Matsubara Green's function can now be obtained by simply
Fourier transforming the operators in the Lagrangian and then the
inverse of the Green's function, $G^{0}_{\alpha\beta}(i\omega)$ is
simply the term multiplying $c^{+}_{\alpha}c_{\beta}$.  Thus we obtain
\begin{equation}
G^{0}_{\alpha\beta}(\omega) = \left[i\omega O-H\right]^{-1}_{\alpha\beta}.
\end{equation}
The renormalized Green's function is now obtained just as in the
orthogonal case by adding a self energy to the Hamiltonian and thus
\begin{equation}
G(i\omega) = \left[i\omega O - H - \Sigma\right].
\end{equation}
We should remark here that these Green's function do not share the
same properties as their cousins in the orthogonal bases do and in
particular the total density is not given by trace of $G(\tau =
0^{-})$.  To see that we go back to the orthogonal basis where we know
how things work and write the density operator as:
\begin{equation}
\rho    = \sum_{i}c^{+}_{i}c_{i} 
        = \sum_{i\alpha\beta}S^{-1}_{i\beta}c^{+}_{\alpha}
                             c_{\beta}(S^{+})^{-1}_{\alpha i}
        =  \sum_{\alpha\beta}O_{\alpha\beta}c^{+}_{\alpha}c_{\beta}.
\end{equation}
Thus the total density of electrons in the system is:
\begin{equation}
n_{tot} = \langle\rho\rangle
        =  \sum_{\alpha\beta}O_{\alpha\beta}
           \langle c^{+}_{\alpha}c_{\beta}\rangle       
        =  \sum_{\alpha\beta}G_{\beta\alpha}(\tau=0^-)O_{\alpha\beta}
\end{equation}
We should note in particular that this means that there seems to be no
good way of assigning a density to a particular orbital in the non
orthogonal case.


\section{Transport Calculation}
\label{sec:n.o.trans}
\marginlabel{sec:n.o.trans}
Now that we know a little bit about many body theory in non orthogonal
bases we can turn to the calculation of the transport function.  At
first one might be tempted to simply use the old expression:
\begin{equation}
\phi^{\mu\nu}(\epsilon) = 
\trace\{v^{\mu}\rho(\epsilon)v^{\nu}\rho(\epsilon)\}
\end{equation}
where the velocity is as before given by the relation:
\begin{equation}
\vec{v}_{\alpha\beta} = \vec{\nabla}H_{\alpha\beta}.
\end{equation}
The following exercise will however convince us that this is not quite
the correct expression for the velocity and will also point the way to
the correct formula.  Let us therefore consider the effect of changing
the zero of energy in the Hamiltonian and moving it up by
$\epsilon_{0}$.  The way this should affect the transport function is
simply by shifting it rigidly up the energy axis.  More concretely let
us denote the original Hamiltonian by $H$ and the shifted Hamiltonian
by $H'$.  Due to the overlap we have $H' = H+\epsilon_{0}O$.  Now
omitting the axis labels and defining $u = \vec{\nabla}_{k}O_{k}$ we
have
\begin{eqnarray}
\phi'(\epsilon) &=& \trace\{v'\rho'(\epsilon)v'\rho'(\epsilon)\}
\nonumber \\ &=&
\trace\{(v+\epsilon_{0}u)\rho(\epsilon-\epsilon_{0})
(v+\epsilon_{0}u)\rho(\epsilon-\epsilon_{0})\}
\nonumber \\ &=&
\phi(\epsilon-\epsilon_{0}) + 
\trace\{\epsilon_{0}u\rho(\epsilon-\epsilon_{0})
(2v+\epsilon_{0}u)\rho(\epsilon-\epsilon_{0})\}.
\end{eqnarray}
Thus we see that the effect of this change was not just a simple
translation of the transport function but in addition we have now an
extra term.  We see from this however that if we define the velocity
in the non orthogonal basis as:
\begin{equation}
\vec{v}^{no}_{k}(\epsilon) = \vec{\nabla}_{k}(H-\epsilon O) 
= \vec{v}-\epsilon\vec{u}
\end{equation}
then the transport function will behave as expected.  Let us now derive
this expression more rigorously.  We follow the exposition of 
\cite{Moreno:1996} which gives gauge-theoretic derivation of the 
heat current for strongly correlated systems.  The starting point for
this derivation is the action for the system which in the non orthogonal
basis is given by:
\begin{equation}
S = \int d\tau {\cal L}(\tau) = 
\int d\tau \sum_{k} c_{k\alpha}^{+}
\left(
O_{k\alpha\beta}\lrpart{\tau}+
H_{k\alpha\beta}
\right)c_{k\beta}.
\end{equation}
Here $\lrpart{\tau} = {1 \over 2}(\rpart{\tau}-\lpart{\tau})$
denotes the anti-symmetrized time derivative.  The particle and heat
currents can now be obtained by considering the invariance of the
action under local phase transformation and local translations in time
respectively.  In the orthogonal case one is lead to the following
expression for the currents:
\begin{equation}
\vec{j} = -{\partial H[\vec{A}_{p}]
\over
\partial \vec{A}_{p}}|_{\vec{A}_{p} = 0} 
\qquad\mbox{and}\qquad
\vec{Q} = -{\partial H[\vec{A}_{h}]
\over
\partial \vec{A}_{p}}|_{\vec{A}_{h} = 0} 
\end{equation}
where $\vec{A}_{p}$ and $\vec{A}_{h}$ are gauge fields conjugate to
the currents and $H[\vec{A}_{p}]$ and $H[\vec{A}_{h}]$ denote the
gauged Hamiltonian, i.e. the Hamiltonian with the replacements
$\vec{k} \rightarrow \vec{k}-\vec{A}_{p}$ and $\vec{k} \rightarrow
\vec{k}+\vec{A}_{h}\lrpart{\tau}$ respectively\footnote {Note that
there is a plus for the second replacement instead of a minus as in
\cite{Moreno:1996}.  This is due to the fact that here we are using
complex time $\tau = it$.}  This replacement is performed in both the
kinetic and the interaction terms but not in the field operators.  In
our case however the overlap matrix appearing in the action depends
also on momentum and therefore the proper generalization of the
currents to non orthogonal basis will also take the overlap matrix
into account.  Thus we obtain:
\begin{equation}
\vec{j} = 
-\frac{\partial (O[\vec{A}_{p}]\lrpart{\tau}+H[\vec{A}_{p}])}
{\partial \vec{A}_{p}}|_{\vec{A}_{p} = 0} 
\qquad\mbox{and}\qquad
\vec{Q} = -\frac{\partial (O[\vec{A}_{h}]\lrpart{\tau}+H[\vec{A}_{h}])}
{\partial \vec{A}_{h}}|_{\vec{A}_{h} = 0} 
\end{equation}
Performing these operations lead to the following expressions:
\begin{eqnarray}
\vec{j} &=& \sum_{k\alpha\beta}
c_{k\alpha}^{+}\left(
\vec{\nabla}_{k}O_{\alpha\beta}\lrpart{\tau}+
\vec{\nabla}_{k}H_{\alpha\beta}
\right)c_{k\beta}
\nonumber \\ &=&
\sum_{k\alpha\beta}
c_{k\alpha}^{+}\left(
\vec{u}_{\alpha\beta}\lrpart{\tau}+
\vec{v}_{\alpha\beta}
\right)c_{k\beta}
\\
\vec{Q} &=& \sum_{k\alpha\beta}
c_{k\alpha}^{+}\left(
-\vec{\nabla}_{k}O_{\alpha\beta}\left(\lrpart{\tau}\right)^{2}
-\vec{\nabla}_{k}H_{\alpha\beta}\lrpart{\tau}
\right)c_{k\beta}
\nonumber \\ &=&
\sum_{k\alpha\beta}
c_{k\alpha}^{+}\left(
-\vec{u}_{\alpha\beta}\left(\lrpart{\tau}\right)^{2} 
-\vec{v}_{\alpha\beta}\lrpart{\tau}
\right)c_{k\beta}
\end{eqnarray}
Let us now define the following operators:
\begin{equation}
B_{\alpha\beta}^{(n)} = 
(-1)^{n}c_{\alpha}^{+}\left(\lrpart{\tau}\right)^{n}c_{\beta}
={1 \over 2^{n}}\sum_{p = 0}^{n}\binom{n}{p}(-1)^{p}
\left(c_{\alpha}^{+}\right)^{(n-p)}c_{\beta}^{(p)} 
\end{equation}
The currents can then be written as:
\begin{eqnarray}
\vec{j} &=& \sum_{k\alpha\beta}
\left(
\vec{v}_{\alpha\beta}B^{(0)}_{\alpha\beta}-
\vec{u}_{\alpha\beta}B^{(1)}_{\alpha\beta}
\right)
\\
\vec{Q} &=& \sum_{k\alpha\beta}
\left(
\vec{v}_{\alpha\beta}B^{(1)}_{\alpha\beta}-
\vec{u}_{\alpha\beta}B^{(2)}_{\alpha\beta}
\right)
\end{eqnarray}
To obtain the transport coefficients we must now calculate the following 
current-current correlation functions, sum over momentum and repeated
orbital indices implied:
\begin{eqnarray}
\langle j^{\mu}(\tau)j^{\nu}(0)\rangle &=&
v^{\mu}_{\alpha\beta}v^{\nu}_{\gamma\delta}
\langle 
B^{(0)}_{k\alpha\beta}(\tau)B^{(0)}_{k\gamma\delta}(0)
\rangle
-v^{\mu}_{\alpha\beta}u^{\nu}_{\gamma\delta}
\langle 
B^{(0)}_{k\alpha\beta}(\tau)B^{(1)}_{k\gamma\delta}(0)
\rangle
\nonumber \\ && -
u^{\mu}_{\alpha\beta}v^{\nu}_{\gamma\delta}
\langle 
B^{(1)}_{k\alpha\beta}(\tau)B^{(0)}_{k\gamma\delta}(0)
\rangle+
u^{\mu}_{\alpha\beta}u^{\nu}_{\gamma\delta}
\langle 
B^{(1)}_{k\alpha\beta}(\tau)B^{(1)}_{k\gamma\delta}(0)
\rangle 
\\ 
\langle j^{\mu}(\tau)Q^{\nu}(0)\rangle &=&
v^{\mu}_{\alpha\beta}v^{\nu}_{\gamma\delta}
\langle 
B^{(0)}_{k\alpha\beta}(\tau)B^{(1)}_{k\gamma\delta}(0)
\rangle
-v^{\mu}_{\alpha\beta}u^{\nu}_{\gamma\delta}
\langle 
B^{(0)}_{k\alpha\beta}(\tau)B^{(2)}_{k\gamma\delta}(0)
\rangle
\nonumber \\ && -
u^{\mu}_{\alpha\beta}v^{\nu}_{\gamma\delta}
\langle 
B^{(1)}_{k\alpha\beta}(\tau)B^{(1)}_{k\gamma\delta}(0)
\rangle+
u^{\mu}_{\alpha\beta}u^{\nu}_{\gamma\delta}
\langle 
B^{(1)}_{k\alpha\beta}(\tau)B^{(2)}_{k\gamma\delta}(0)
\rangle 
\\ 
\langle Q^{\mu}(\tau)Q^{\nu}(0)\rangle &=&
v^{\mu}_{\alpha\beta}v^{\nu}_{\gamma\delta}
\langle 
B^{(1)}_{k\alpha\beta}(\tau)B^{(1)}_{k\gamma\delta}(0)
\rangle
-v^{\mu}_{\alpha\beta}u^{\nu}_{\gamma\delta}
\langle 
B^{(1)}_{k\alpha\beta}(\tau)B^{(2)}_{k\gamma\delta}(0)
\rangle
\nonumber \\ && -
u^{\mu}_{\alpha\beta}v^{\nu}_{\gamma\delta}
\langle 
B^{(2)}_{k\alpha\beta}(\tau)B^{(1)}_{k\gamma\delta}(0)
\rangle+
u^{\mu}_{\alpha\beta}u^{\nu}_{\gamma\delta}
\langle 
B^{(2)}_{k\alpha\beta}(\tau)B^{(2)}_{k\gamma\delta}(0)
\rangle 
\end{eqnarray}
Note here that we have already taken into account that the correlation
functions are diagonal in momentum and thus one of the dummy momentum 
indices is gone.
We must now evaluate the correlation functions:
\begin{equation}
\Gamma^{mn}_{k\alpha\beta\gamma\delta}(\tau) = 
\langle B^{(m)}_{k\alpha\beta}(\tau)B^{(n)}_{k\gamma\delta}\rangle
\end{equation}
and that we do to lowest order by evaluating a single Feynman
diagram:\footnote{See Appendix \ref{app:spectral} for a discussion on the
  Green's function and its derivatives.}
\begin{eqnarray}
\Gamma^{mn}_{k\alpha\beta\gamma\delta}(\tau) &=& 
{1 \over 2^{m+n}}\sum_{p=0}^{m}\sum_{q=0}^{n}
\binom{m}{p}\binom{n}{q}(-1)^{p+q}
\nonumber \\ && \times
\left\langle 
\left(c^{+}_{\alpha}\right)^{(m-p)}(\tau)
c_{\beta}^{(p)}(\tau)
\left(c^{+}_{\gamma}\right)^{(n-q)}
c_{\delta}^{(q)}
\right\rangle
\nonumber \\ &=&
{1 \over 2^{m+n}}\sum_{p=0}^{m}\sum_{q=0}^{n}
\binom{m}{p}\binom{n}{q}(-1)^{p+q}
\nonumber \\ && \times
\left\langle 
\left(c^{+}_{\alpha}\right)^{(m-p)}(\tau)
c_{\delta}^{(q)}
\right\rangle
\left\langle 
c_{\beta}^{(p)}(\tau)
\left(c^{+}_{\gamma}\right)^{(n-q)}
\right\rangle
\nonumber \\ &=& 
{1 \over 2^{m+n}}\sum_{p=0}^{m}\sum_{q=0}^{n}
\binom{m}{p}\binom{n}{q}(-1)^{p+q}
\nonumber \\ && \times
(-1)^{q}G^{(m-p+q)}_{\delta\alpha}(-\tau)
(-1)^{n-q+1}G^{(n-q+p)}_{\beta\gamma}(\tau)
\nonumber \\ &=&
\int d\epsilon d\epsilon' 
\rho_{k\beta\gamma}(\epsilon)\rho_{k\delta\alpha}(\epsilon')
{1 \over 2^{m+n}}\sum_{p=0}^{m}\sum_{q=0}^{n}
\binom{m}{p}\binom{n}{q}\epsilon^{n-q+p}(\epsilon')^{m-p+q}
\nonumber \\ &=&
\int d\epsilon d\epsilon' 
\rho_{k\beta\gamma}(\epsilon)\rho_{k\delta\alpha}(\epsilon')
e^{(\epsilon'-\epsilon)\tau}f(-\epsilon)f(\epsilon')
\left(\frac{\epsilon+\epsilon'}{2}\right)^{m+n}
\end{eqnarray}
Let us now define the transport function as:
\begin{equation}
\phi^{\mu\nu}(\epsilon,\epsilon') = \sum_{k}\trace\left\{
\left[
v^{\mu}_{k}-\left(\frac{\epsilon+\epsilon'}{2}\right)u^{\mu}_{k}
\right]
\rho_{k}(\epsilon)
\left[
v^{\nu}_{k}-\left(\frac{\epsilon+\epsilon'}{2}\right)u^{\nu}_{k}
\right]
\rho_{k}(\epsilon')
\right\}
\end{equation}
Then we can write the current-current correlation functions as
\begin{eqnarray}
\langle j^{\mu}(\tau)j^{\nu}(0) \rangle &=&
\int d\epsilon d\epsilon' 
\phi^{\mu\nu}(\epsilon,\epsilon')
e^{(\epsilon'-\epsilon)\tau}f(-\epsilon)f(\epsilon')
\\  
\langle j^{\mu}(\tau)j^{\nu}(0) \rangle &=&
\int d\epsilon d\epsilon' 
\phi^{\mu\nu}(\epsilon,\epsilon')
e^{(\epsilon'-\epsilon)\tau}f(-\epsilon)f(\epsilon')
\left(\frac{\epsilon+\epsilon'}{2}\right)
\\  
\langle j^{\mu}(\tau)j^{\nu}(0) \rangle &=&
\int d\epsilon d\epsilon' 
\phi^{\mu\nu}(\epsilon,\epsilon')
e^{(\epsilon'-\epsilon)\tau}f(-\epsilon)f(\epsilon')
\left(\frac{\epsilon+\epsilon'}{2}\right)^{2}
\end{eqnarray}
For simplicity let us now symmetrize these functions with respect to
the indices $\mu$ and $\nu$.  The symmetrized transport function will
then be real, due to the cyclicity of the trace and the hermiticity of
the matrices appearing in the trace.  It is now a simple matter to
obtain the retarded version of the symmetrized correlation functions
and that calculation proceeds exactly as before in \article\ref{sec:trans-eq}.  We should note here that the what we usually
refer to as the transport function is the value of $\phi$ on the
diagonal, i.e.
\begin{equation}
\phi^{\mu\nu}(\epsilon) = \phi^{\mu\nu}(\epsilon,\epsilon)
\end{equation}
and this function is of course real and we don't need any symmetrization.
It is only when the two arguments are not equal that we need to symmetrize
to be able to simplify the retarded version of the correlation function.

\section{Transport function calculation}
\label{sec:n.o.transfunc}
\marginlabel{sec:n.o.transfunc} We see from the previous section that the main
task is the evaluation of the transport function, i.e. the computation of the
Brillouin zone integral.  This calculation is almost exactly the same as
discussed in \article\ref{sec:transfunc} and \article\ref{sec:comp-phi} but
since we are dealing with a non orthogonal basis it is worth commenting on the
numerical diagonalization of the Hamiltonian.  The eigenvalue problems that we
have to solve are the following:
\begin{equation}
H_{\alpha\beta}\Psi^{R,p}_{\beta} = 
E^{p}O_{\alpha\beta}\Psi^{R,p}_{\beta}
\qquad\mbox{and}\qquad
\Psi^{L,p}_{\alpha}H_{\alpha\beta} = 
E^{p}\Psi^{L,p}_{\alpha}O_{\alpha\beta}.
\end{equation}
We should note here that we are now solving the Hamiltonian which
includes  the self energy so the Hamiltonian is no longer Hermitian
but the overlap  matrix is of course Hermitian.  The strategy is to
reduce this problem to a regular eigenvalue problem and in principle
one could simply find the inverse of the overlap matrix and calculate
the matrix representation of the Hamiltonian, $\hat{H}$ and then solve
the regular eigenvalue problem.  A better solution is to use the fact
that the overlap matrix is both Hermitian and positive definite and
thus we can use Cholesky decomposition and  write:
\begin{equation}
O = U^{+}U
\end{equation}
where $U$ is an upper triangular matrix.  Now $U$ can be inverted quite
efficiently, so let
\begin{equation}
B = U^{-1}.
\end{equation}
We now define the transformed Hamiltonian as
\begin{equation}
H' = B^{+}HB
\end{equation}
and then solve the regular left and right eigenvalue problems for this
matrix.  If we denote those eigenvectors by $\bar{v}^{L}$ and
$\bar{v}^{R}$ then our back-transformed eigenvectors are given by
\begin{equation}
\bar{v}^{L}B^{+} 
\qquad\mbox{and}\qquad
B\bar{v}^{R}
\end{equation}
We should note that the left and the right eigenvectors corresponding
to different eigenvalues are orthogonal in the following sense:
\begin{equation}
\Psi^{L,p}_{\alpha}O_{\alpha\beta}\Psi^{R,p'}_{\beta} = 0
\end{equation}
when $p \neq p'$.  We can choose to normalize them in such a way that
the equation above becomes:
\begin{equation}
\Psi^{L,p}_{\alpha}O_{\alpha\beta}\Psi^{R,p'}_{\beta} = \delta_{pp'}.
\end{equation}
This is easily seen considering the following:
\begin{equation}
\Psi^{L,p}_{\alpha}H_{\alpha\beta}\Psi^{R,p'}_{\beta} = 
E^{p}\Psi^{L,p}_{\alpha}O_{\alpha\beta}\Psi^{R,p'}_{\beta} =
E^{p'}\Psi^{L,p}_{\alpha}O_{\alpha\beta}\Psi^{R,p'}_{\beta}
\end{equation}
and thus if $p \neq p'$ we arrive at the conclusion above.  Let us now
define the eigenvector matrices $A^{L}$ and $A^{R}$ such that:
\begin{equation}
A^{L}_{p,\alpha} = \Psi^{L,p}_{\alpha}
\qquad\mbox{and}\qquad
A^{L}_{\alpha,p} = \Psi^{P,p}_{\alpha}
\end{equation}
i.e. the rows of $A^{L}$ are the left eigenvectors and the columns of
$A^{R}$ are the right eigenvectors.  Then we see that $A^{L}OA^{R}$
equals the identity matrix.  It is now easy to see by inspection that
the Hamiltonian can be written on the form:
\begin{equation}
H = OA^{R}EA^{L}O
\end{equation}
where $E$ is a diagonal matrix of eigenvalues.  The Green's function
can therefore be written as:
\begin{eqnarray}
G(z) &=& [zO-H]^{-1} = [zO-OA^{R}EA^{L}O]^{-1}
\nonumber \\ &=&
[OA^{R}(zI-E)A^{L}O]^{-1} = [A^{L}O]^{-1}[zI-E]^{-1}[OA^{R}]^{-1} 
\nonumber \\ &=&
A^{R}[zI-E]^{-1}A^{L}.
\end{eqnarray}
Here $z$ can be any complex value outside the real axis.  We note of
course that the self energy and therefore the Hamiltonian depends on
$z$ and thus both the eigenvectors and eigenvalues depend on $z$.  In
particular we can consider $z = \omega+i0$ in which case we obtain the
retarded Green's function and from there we get the spectral density:
\begin{eqnarray}
\rho(\omega) &=& {-1 \over 2\pi i }
\left\{
G^{R}(\omega)-\left(G^{R}\right)^{+}
\right\}
\nonumber \\ &=& 
{-1 \over 2\pi i }
\left\{
A^{R}[\omega I-E]^{-1}A^{L} - 
\left(A^{L}\right)^{+}[\omega I-E^{+}]^{-1}\left(A^{R}\right)^{+}  
\right\}.
\end{eqnarray}
Putting this into the expression for the transport function we find again 
that it can be expressed as:
\begin{equation}
\phi^{\mu\nu}(\omega) = {-1 \over 2\pi^2}\real\trace\left\{
r^{\mu}_{qp}r^{\mu}_{pq}D_{p}D_{q}-
{1 \over 2}\left(s^{\mu}_{qp}t^{\nu}_{pq}+s^{\nu}_{pq}t^{\mu}_{qp}\right)
D_{p}D^{+}_{q}
\right\}
\end{equation}
where the matrices $r$, $s$ and $t$ are given by 
\begin{equation}
r^{\mu} = A^{L}v^{\mu}A^{R}, 
\qquad
s^{\mu} = A^{L}v^{\mu}(A^{L})^{+},
\qquad\mbox{and}\qquad
t^{\mu} = (A^{R})^{+}v^{\mu}A^{R}
\end{equation}
and the diagonal elements $D$ are given in terms of the eigenvalues of
the Hamiltonian as:
\begin{equation}
D_{p} = {1 \over \omega+i0-E_{p}}.
\end{equation}

\section{An alternative look at the velocity}
\label{sec:alt.vel}
\marginlabel{sec:alt.vel}
In this section we will take another look at the definition of the
current operator and instead of arriving at it from a gauge invariance
argument let us simply look at the standard quantum mechanical
representation of the particle current:
\begin{equation}
\vec{j}(r) = \frac{1}{2mi}
\left\{
\psi^{+}(r)\vec{\nabla}_{r}\psi(r)-[\vec{\nabla}_{r}\psi^{+}(r)]\psi(r)
\right\}.
\end{equation}
In our case we are interested in the zero momentum component of the
current which is simply the integral over $r$ of this quantity and
then we see that the current can simply be regarded as the real part 
of the expectation value of the operator:
\begin{equation}
{\cal V} = \frac{1}{im}\vec{\nabla}_{r}.
\end{equation}
Let us now see how the matrix elements of this quantity looks like in
the crystal within the LMTO single electron formalism.  The basis
functions we take are of the form \cite{Skriver:LMTO}:
\begin{eqnarray}
\chi^{\alpha}_{k}(r) &=& 
\sum_{lm}a_{lm}^{\alpha k}\sum_{R}e^{ikR}\chi_{lm}(E,r-R)
\nonumber \\ &=&
e^{ikr}\sum_{lm}a_{lm}^{\alpha k}\sum_{R}e^{-ik(r-R)}\chi_{lm}(E,r-R)
\nonumber \\ &=&
e^{ikr}u^{\alpha}_{k}(r).
\end{eqnarray}
Here we have, in keeping with Bloch's theorem, written the basis
function in terms of a traveling wave multiplied by a function that
has the periodicity of the lattice.
The Hamiltonian we are trying to solve is nothing but a single particle
Hamiltonian:
\begin{equation}
H = -\frac{\nabla_{r}^{2}}{2m}+V(r)
\end{equation}
and here $V(r)$ is the periodic potential due to the atoms in the
lattice and the rest of the electrons.  Acting with this Hamiltonian
on the basis functions gives:
\begin{eqnarray}
H\chi^{\alpha}_{k}(r) &=& 
\left(-\frac{\nabla_{r}^{2}}{2m}+V(r)\right)e^{ikr}u^{\alpha}_{k}(r) 
\nonumber \\ &=&
e^{ikr}\left(
-\frac{(\vec{\nabla}_{r}+ik)^{2}}{2m}+
V(r)\right)u^{\alpha}_{k}(r)
\equiv e^{ikr}H_{k}u^{\alpha}_{k}(r).
\end{eqnarray}
The matrix elements can now be written as:
\begin{equation}
H_{k\alpha\beta}  
= \langle \chi_{k}^{\alpha}|H|\chi_{k}^{\beta}\rangle
= \langle u_{k}^{\alpha}|H_{k}|u_{k}^{\beta}\rangle.
\end{equation}
Let us now take the gradient of this operator with respect to momentum:
\begin{eqnarray}
\grad{k}H_{k\alpha\beta}  
&=& \langle \grad{k}u_{k}^{\alpha}|H_{k}|u_{k}^{\beta}\rangle
+\langle u_{k}^{\alpha}|H_{k}|\grad{k}u_{k}^{\beta}\rangle
+\langle u_{k}^{\alpha}|\frac{\grad{r}+ik}{im}|u_{k}^{\beta}\rangle.
\end{eqnarray}
We see now that:
\begin{equation}
\langle u_{k}^{\alpha}|\frac{\grad{r}+ik}{im}|u_{k}^{\beta}\rangle
= 
\langle \chi_{k}^{\alpha}|\frac{\grad{r}}{im}|\chi_{k}^{\beta}\rangle
= {\cal V}_{k\alpha\beta} = \vec{j}_{k\alpha\beta}
\end{equation}
and thus
\begin{equation}
\vec{j}_{k\alpha\beta} = \grad{k}H_{k\alpha\beta}  
-\langle \grad{k}u_{k}^{\alpha}|H_{k}|u_{k}^{\beta}\rangle
-\langle u_{k}^{\alpha}|H_{k}|\grad{k}u_{k}^{\beta}\rangle.
\end{equation}
Now we do not necessarily assume that the basis is orthogonal and the 
overlaps of the basis functions are given by the overlap matrix:
\begin{equation}
O_{k\alpha\beta} = \langle\chi_{k}^{\alpha}|\chi_{k}^{\beta}\rangle = 
\langle u_{k}^{\alpha}|u_{k}^{\beta}\rangle 
\end{equation}
and the resolution of the identity in terms of the $u_{k}^{\alpha}$
basis is given by
\begin{equation}
|u^{\alpha}_{k}\rangle (O^{-1})_{k\alpha\beta}\langle u^{\beta}_{k}|
\end{equation}
as can easily be seen by simply acting with this operator on one of
the $u$ functions.  We now insert this into the expression for the 
current operator and get:
\begin{equation}
\vec{j}_{k\alpha\beta} = \grad{k}H_{k\alpha\beta}  
-\langle \grad{k}u_{k}^{\alpha}|u_{k}^{\gamma}\rangle 
(O^{-1})_{k\gamma\delta}H_{k\delta\beta}
-H_{k\alpha\gamma}(O^{-1})_{k\gamma\delta}
\langle u_{k}^{\delta}|\grad{k}u_{k}^{\beta}\rangle.
\end{equation}
Now we note that:
\begin{equation}
\grad{k}O_{k\alpha\beta} = 
\langle \grad{k}u_{k}^{\alpha}|u_{k}^{\beta}\rangle +
\langle u_{k}^{\alpha}|\grad{k}u_{k}^{\beta}\rangle
\end{equation}
and thus we get
\begin{eqnarray}
\vec{j}_{k\alpha\beta} &=& 
\grad{k}H_{k\alpha\beta} - 
(\grad{k}O_{k\alpha\gamma})(O^{-1})_{k\gamma\delta}H_{k\delta\beta}
\nonumber \\ &&
+\langle u_{k}^{\alpha}|\grad{k}u_{k}^{\gamma}\rangle 
(O^{-1})_{k\gamma\delta}H_{k\delta\beta}
-H_{k\alpha\gamma}(O^{-1})_{k\gamma\delta}
\langle u_{k}^{\delta}|\grad{k}u_{k}^{\beta}\rangle.
\end{eqnarray}
We have now finally succeeded in writing the current in terms of the
gradient of the Hamiltonian along with some other terms.  We note that
the third and fourth terms in this expression are invariant with
respect to a shift in the energy level and the second term ensures the
same invariance for the first two terms.

The second term is quite similar to the extra term we obtained in the
current in \article\ref{sec:n.o.trans} and plays the same role.  The
other additional terms are similar to what is obtained if one carries
out a momentum dependent rotation in the space of orbitals.

Let us for simplicity define the following quantities:
\begin{equation}
\tilde{H}_{k\alpha\beta} = O^{-1}_{k\alpha\gamma}H_{k\gamma\beta}
\qquad\mbox{and}\qquad
\vec{A}_{k\alpha\beta} = 
\langle u_{k}^{\alpha}|\grad{k}u_{k}^{\beta}\rangle.
\end{equation}
Then the current matrix element can be written as:
\begin{equation}
\label{eq:ClassicalCurrent}
\vec{j}_{k\alpha\beta} = 
\grad{k}H_{k\alpha\beta} - 
(\grad{k}O_{k\alpha\gamma})\tilde{H}_{k\gamma\beta}
+\vec{A}_{k\alpha\gamma}\tilde{H}_{k\gamma\beta}
-\tilde{H}^{+}_{k\alpha\gamma}\vec{A}_{k\gamma\beta}
\end{equation}
Since this expression for the current matrix elements looks surprisingly
similar to the one that we derived previously using the local gauge invariance
method it is worthwhile to try to see if we can show their equivalence.  So
let us consider the complex time, noninteracting, Schr\"odinger equation:
\begin{equation}
H|\Psi\rangle = \rpart{\tau}|\Psi\rangle
\qquad\mbox{and}\qquad
\langle\Psi |H = -\langle\Psi |\lpart{\tau}.
\end{equation}
Now expanding the wave function in terms of the basis functions that
we are using and looking at the matrix elements we get (sum over
repeated indices implied)
\begin{equation}
H_{\alpha\beta}c^{+}_{\beta} = O_{\alpha\beta}\dot{c}^{+}_{\beta}
\qquad\mbox{and}\qquad 
c_{\alpha}H_{\alpha\beta} = -\dot{c}_{\alpha}O_{\alpha\beta}.
\end{equation}
Hence we see that 
\begin{equation}
  O^{-1}_{\gamma\alpha}H_{\alpha\beta}c^{+}_{\beta} 
  = \tilde{H}_{k\gamma\beta}c^{+}_{\beta}  
  = \dot{c}^{+}_{\gamma}
  \qquad\mbox{and}\qquad 
  c_{\alpha}H_{\alpha\beta}O^{-1}_{\beta\gamma} 
  = c_{\alpha}\tilde{H}^{+}_{\alpha\gamma}
  = -\dot{c}_{\gamma}.
\end{equation}
Now plugging these expression into equation (\ref{eq:ClassicalCurrent}) we
get:
\begin{equation}
  \vec{j}_{k\alpha\beta} = 
  \grad{k}H_{k\alpha\beta} - 
  (\grad{k}O_{k\alpha\beta})\rpart{\tau}
  +\vec{A}_{k\alpha\beta}\rpart{\tau}
  +\lpart{\tau}\vec{A}_{k\alpha\beta}.
\end{equation}
If we now decide to anti-symmetrize the time derivative, we see that the last
two terms do indeed vanish and we end up with exactly the same expression as 
we had from the previous derivation of the current.

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

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