[tex] more numerics

tex/bib.bib

@@ -75,4 +75,54 @@ @Article{PhysRevB.64.024426
   publisher = {American Physical Society}
 }
 
+@article{matrixinversion,
+ author = {Wilkinson, J. H.},
+ title = {Error Analysis of Direct Methods of Matrix Inversion},
+ journal = {J. ACM},
+ volume = {8},
+ number = {3},
+ year = {1961},
+ issn = {0004-5411},
+ pages = {281--330},
+ doi = {http://doi.acm.org/10.1145/321075.321076},
+ publisher = {ACM},
+ address = {New York, NY, USA},
+}
+
+@book{matrixperformance,
+  Author = {Roger A. Horn and Charles R. Johnson},
+  Title = {Matrix Analysis},
+  Publisher = {Cambridge University Press},
+  Year = {1990},
+  ISBN = {0521386322},
+}
+
+@article{rgfschmelcher,
+ author = {Drouvelis, P. S. and Schmelcher, P. and Bastian, P.},
+ title = {Parallel implementation of the recursive Green's function method},
+ journal = {J. Comput. Phys.},
+ volume = {215},
+ number = {2},
+ year = {2006},
+ issn = {0021-9991},
+ pages = {741--756},
+ doi = {http://dx.doi.org/10.1016/j.jcp.2005.11.010},
+ publisher = {Academic Press Professional, Inc.},
+ address = {San Diego, CA, USA},
+}
+
+@Article{fischer-lee,
+  title = {Relation between conductivity and transmission matrix},
+  author = {Fisher, Daniel S. and Lee, Patrick A.},
+  journal = {Phys. Rev. B},
+  volume = {23},
+  number = {12},
+  pages = {6851--6854},
+  numpages = {3},
+  year = {1981},
+  month = {Jun},
+  doi = {10.1103/PhysRevB.23.6851},
+  publisher = {American Physical Society}
+}
+
 

tex/numerics.tex

@@ -23,5 +23,80 @@
 \section*{Performance consideration}
 
 These numeric calculations are computationally expensive. If we simulate a
-lattice with a total of $n$ sites, the matrices $H$, $\Sigma_p$, $G^R$ and
-$G^A$ are of dimensions $2n^2 \times 2n^2$
+lattice with a total of $n \x n$ sites, the matrices $H$, $\Sigma_p$, $G^R$ and
+$G^A$ are of dimensions $2n^2 \times 2n^2$. Inversion of an $m \times m $
+matrix and multiplication of two $ m \times m $ matrices both take $O(m^3)$
+steps\footnote{There are matrix product algorithms with slightly better
+asymptotic scaling, but they are usually very complicated, numerically badly
+conditioned, or only advantageous for huge systems; usually all three apply}
+\cite{matrixperformance}. All other steps are faster, and can be ignored for
+an asymptotic analysis.
+
+So the naive approach takes $O(n^6)$ steps for $n$ lattice sites.
+
+An alternative approach, the \emph{Recursive Green's Function method} works by
+decomposing the sample into slices in such a way that the sites in each slice
+only interact with neighbor slices. For short range interactions such a
+decomposition can be found, and for nearest-neighbor hopping the time
+complexity approaches $O(n^3)$ \cite{rgfschmelcher}. However this is payed by
+a significantly increased implementation complexity, and the algorithm is tied
+to the specific form of the Hamiltonian.
+
+\begin{figure}
+    \includegraphics[angle=270,width=0.7\textwidth]{scaling.pdf}
+    \caption{Scaling of run time with system size for the tight binding
+        simulation, implemented with sparse matrices and the SuperLU sparse
+        direct solver. The data was recorded for square samples including
+        Rashba spin-orbit coupling, and 4 leads of the same width as the
+        sample, on a 2.9GHz AMD64 computer with SSE and 8GB RAM.
+        The run time approximately scales as $t(n) = 2
+        n^{3.64}\mu s$}
+        \label{fig:scaling}
+\end{figure}
+
+We chose a middle ground: the "naive" approach, but implemented with sparse
+matrices, and using an efficient sparse direct solver\cite{superlu99} for
+computing the Green's functions.
+
+The run time thus depends largely on the number of non-zeros in the $H$ and
+$\Sigma_p$, and thus on the nature of the interaction. For nearest neighbor
+hopping and Rashba spin-orbit coupling we observed a run time scaling of
+$O(n^{3.64})$ (See Fig. \ref{fig:scaling}).
+
+\section*{Numerical stability}
+
+The numerical
+calculation uses floating point numbers with limited machine precision,
+therefore numerical errors are inevitable.
+
+The transmission matrix for a sample which is connected to multiple uniform
+leads (with same number of modes per lead) has to fulfill the sum rules
+
+\begin{align}
+    \sum_p T_{pq} = \sum_q T_{pq} = M
+\end{align}
+
+Where $M$ is the number of modes (and thusly and integer). We can calculate
+this left hand side of this equation from the numerical simulation, round it
+to the nearest integer and thus estimate the numerical error in $T_{pq}$.
+
+Since generally matrix inversion is numerically worse conditioned than solving
+a linear equation system\cite{matrixinversion}, we don't actually compute
+$G^A$ and $G^R$, but rather the products $G^A\Sigma_p$ and $G^R\Sigma_q$,
+which can be written as solutions of linear equation systems.
+
+\begin{align}
+    X_p &= \Sigma_p G^R\\
+    X_p^\dagger &= (G^R)^\dagger \Sigma_p^\dagger = G^A \sigma_p^\dagger\\
+    \Rightarrow\quad (G^A)^-1 X_p^\dagger &= \Sigma_p^\dagger
+    \label{eq:lin_gls}
+\end{align}
+
+Eq. \ref{eq:lin_gls} is the form with which linear sparse solvers typically
+work. They typically calculate a LU decomposition (that is they write
+$(G^A)^-1 = L \cdot U$, where $L$ is a lower triangular matrix and $U$ an
+upper triangular matrix), and directly solve the equation by elimination.
+Complex conjugation and transposition of $X_p^\dagger$ finally gives us the
+matrix which we need evaluating the Fischer-Lee relation.
+
+% vim: ts=4 sw=4 expandtab spell spelllang=en_us tw=78

tex/theory.tex

@@ -258,13 +258,15 @@ \subsection*{Transmission and Green's functions}
     \textnormal{with } \Sigma^R_p &= \tau^\dagger g_p^R \tau_p
 \end{align}
 
-It can be shown \cite{baranger} that the transmission $T_{pq}$ from lead
+It can be shown \cite{fischer-lee,baranger} that the transmission $T_{pq}$ from lead
 $q$ through the sample into lead $p$ is
 
 \begin{align}
     T_{pq} = \trace(\Sigma_p G_s^R \Sigma_q G_s^A)
 \end{align}
 
+This is called the Fischer-Lee-Relation.
+
 This allows numeric calculation of transmission coefficients for
 arbitrary Hamiltonian operators, and is thus independent of the individual
 phenomena we're looking at.
@@ -355,7 +357,8 @@ \subsection*{Tight Binding Hamiltonian}
 \end{equation}
 
 To map this to a discrete lattice we substitute the derivation by a discrete
-difference, so $\dell_x f(x)|_{x_0}$ becomes $(f(x_0+a) - f(x_0))/a$. For the
+difference, so $\dell_x f(x)|_{x_0}$ becomes $(f(x_0+a) - f(x_0))/a$. $a$ is
+the lattice constant. For the
 second derivation we chose the symmetric difference $\dell_x^2 f(x)_{x_0} =
 (f(x_0+a) - f(x_0-a))/2a$. In the continuum limit $a \mapsto 0$ both reproduce
 the derivation exactly.
@@ -378,7 +381,7 @@ \subsection*{Tight Binding Hamiltonian}
 $\downarrow$, and $a_x$ and $a_y$ denote the shift to the nearest neighbor in
 $x$ and $y$ direction, respectively. (Note that we assume that the lattice is
 equally spaced in $x$ and $y$ direction, $+a_x$ just means "go to the next
-neighbor in $x$ direction). $t$ is the so-called hopping term, and denotes the
+neighbor in $x$ direction). $t = \frac{\hbar^2}{2ma}$ is the so-called hopping term, and denotes the
 probability of an electron traveling to its nearest neighbor. $r$ is the
 Rashba hopping term, and corresponds to a nearest neighbor hop with spin
 flip.