text:
\section{Introduction}

During the last decades, the size of electronic circuits has continuously
been reduced. Today, systems like quantum wires and quantum dots are routinely
produced on the nanometer scale. Recently, the seemingly ultimate limit of
minituarization has been achieved by several experimental groups who were
able to place single molecules between two macroscopic electrodes.
\cite{ReedZhouMullerBurginTour:97,
SansDevoretDaiThessSmalleyGeerligsDekker:97}
To describe transport properties on such a small scale, a quantum theory
of transport is required.\cite{jortner:97,book}

A cornerstone of such a theory is the Landauer-B\"uttiker formalism
\cite{Landauer:57,Buettiker:86} which provides a method to compute the
steady-state current of non-interacting electrons for meso- or nanoscopic
systems connecting two (or more) macroscopic electrodes.

Alternatively, the technique of nonequilibrium Green functions (NEG)
\cite{KadanoffBaym:62,Keldysh:65}
has been used to tackle quantum transport. Studies using the NEG
approach typically use tight-binding-like model Hamiltonians to describe
the combined system electrodes plus ``device''.
A well known scheme is the one introduced by Caroli {\em et al.}
\cite{caroli1,caroli2} In the remote past the left and
right electrodes are disconnected and in equilibrium at two different
chemical potentials; the conducting part of the Hamiltonian is switched
on adiabatically and eventually a steady-state develops.
Within this contacting approach
also time-dependent transport phenomena have been investigated.
\cite{JauhoWingreenMeir:94}
Caroli {\em et al.} discussed non-interacting systems only. Their
approach has later been extended to account for short-range
electron-electron interaction and for interaction with vibrations in
the device region.\cite{MeirWingreen:92} An excellent overview
of the field can be
found in the book by Haug and Jauho \cite{HaugJauho:98} and in
Ref.~\onlinecite{book}. Despite its appeal,
the above scheme has limitations since the
time-dependent perturbation is the tunneling Hamiltonian, connecting
the electrodes to the device, rather than the external electric
field.

Cini proposed another scheme \cite{cini} also based on NEG. Here,
the system electrodes plus ``device'' is
connected and in equilibrium in the remote past.
The time-dependent perturbation is the external scalar potential.
It has been shown \cite{saprb04} that for non-interacting systems the
contacting approach and the Cini approach
yield the same current in the long-time limit and that in the dc case the
steady state current does not depend on the history of the applied potential.
Moreover, the Cini scheme is well suited for a density functional
extension since the electrons are driven out of equilibrium by a local
potential rather than by a non-local one (see below).

With recent experimental progress to place single molecules as devices between
macroscopic electrodes there also has been considerable activity to describe
transport through these systems on an {\em ab initio} level.
Most approaches are based on a self-consistency procedure first
proposed by Lang.\cite{lang} In this steady-state approach based on density
functional theory (DFT), exchange and correlation is approximated by the
static Kohn-Sham (KS) potential and the charge density is obtained
self-consistently in the presence of the steady current.
However, the original justification
involved subtle points such as different
Fermi levels deep inside the left and right electrodes and the
implicit reference of non-local perturbations such as tunneling
Hamiltonians within a DFT framework. (For a
detailed discussion we refer to Ref. \onlinecite{saprb04}.)
The steady-state DFT approach has been further developed
\cite{DerosaSeminario:01,BrandbygeMozosOrdejonTaylorStokbro:02,XueDattaRatner:02,Calzolari:04}
and the results have been most useful for understanding the
qualitative behavior of measured current-voltage characteristics. Quantitatively,
however, the theoretical I-V curves often differ from the experimental
ones by several orders of magnitude.\cite{DiVentraPantelidesLang:00}
Several explanations are possible for such a mismatch:
models are not sufficiently refined,
parasitic effects in measurements have been underestimated,
the characteristics of the molecule-contact interface
are not well understood and difficult to address given their atomistic complexity.
Adding to the theoretical reason for this discrepancy is the fact that the
transmission functions computed from static DFT have resonances at the
non-interacting Kohn-Sham excitation energies which in general do not
coincide with the true excitation energies. Furthermore, different
exchange-correlation functionals lead to DFT-currents that vary by
more than an order of magnitude.\cite{Krstic:03}

Excitation energies of interacting systems are accessible via time-dependent
(TD) DFT.\cite{RungeGross:84,PetersilkaGossmannGross:96} In
this theory, the time-dependent density of an interacting system moving in
an external, time-dependent local potential can be calculated via a fictitious
system of non-interacting electrons moving in a local, effective
time-dependent potential. Therefore this theory is in principle well suited
for the treatment of nonequilibrium transport problems.\cite{StefanucciAlmbladh:04}
A basic issue is that most exchange-correlation functionals
have been derived under equilibrium conditions and their application to
non-equilibrium problems should be analyzed in more detail. However, this is beyond
the scope of the present work.

Before a TDDFT calculation of transport can be tackled, a number of
technical problems have to be addressed. In particular, one needs a practical
scheme for the propagation of the time-dependent Schr\"odinger equation for
an infinitely large system. Of course, since one can in practice only deal
with finite systems this can only be achieved by applying the correct
boundary conditions. The problem of so-called ``transparent boundary
conditions'' for the time-dependent Schr\"odinger equation has been
attacked by many authors. For a recent overview, the reader is referred to
Ref.~\onlinecite{Moyer:04}.

In this paper we present a propagation scheme which is particularly designed
to be used for the calculation of time-dependent transport problems. In
Section \ref{secII}, we combine the Cini scheme with TDDFT and we
develop a general formalism based on the
propagation of Kohn-Sham orbitals in open systems.
In Section \ref{secIII} we will address the question of how
to obtain the correct initial states for the propagation.
An algorithm for the time-evolution of open systems
is proposed in Section \ref{secIV}. It is based on a
modified version of the Crank-Nicholson algorithm.
Section~\ref{details} describes some details of our numerical
implementation and Section~\ref{examples} gives results
for several one-dimensional model systems.
We draw our conclusions in Section \ref{conc}.

\section{General Formulation}
\label{secII}

We consider an electrode-junction-electrode system which is initially
in equilibrium ($t<0$). The system is contacted and no current flows through
the junction, the charge density of the electrodes being perfectly
balanced, see Fig. \ref{system}.
Therefore, the system initially is in its ground state which, due to the
Hohenberg-Kohn theorem,\cite{HohenbergKohn:64} is a functional of the
density. This density can then be computed in the usual way by
$n({\bf r},0)=\sum_{occ} |\psi_{s}({\bf r},0)|^{2}$ where the sum is over
the occupied Kohn-Sham orbitals $\psi_{s}({\bf r},0)$, i.e., the eigenfunctions
of the Kohn-Sham Hamiltonian $\mbox{\boldmath $H$}(0)$ with eigenergy below the Fermi energy.
Here and in the following we use boldface notation to denote operators
in one-electron Hilbert space.
\begin{figure}[t]
\includegraphics[scale=0.9]{fig1.eps}
\caption{\label{system} Schematic diagram of the system described in the
main text for the calculation of charge transport through the central
constriction.}
\label{curr1}
\end{figure}

To observe a current we drive the system out of equilibrium by
exposing the electrons to an external electric potential (bias).
Without loss of generality we will assume that this
external bias vanishes for times $t\leq0$. According to the Runge-Gross
theorem \cite{RungeGross:84}, the time-dependent density can be
calculated by evolving the KS orbitals according to
the KS equation of TDDFT $i\dot{\psi}_{s}(t)=\mbox{\boldmath $H$}(t)\psi_{s}(t)$ where $\mbox{\boldmath $H$}(t)$ is
the time-dependent KS Hamiltonian. Thus,
$n({\bf r},t)=\sum_{occ}|\psi_{s}({\bf r},t)|^{2}$ and the
continuity equation is
$\frac{\partial}{\partial t}n({\bf r},t)=
-\nabla\cdot{\bf j}_{\rm KS}({\bf r},t)$,
where ${\bf j}_{\rm KS}({\bf r},t)=-\sum_{occ}{\rm Im}[
\psi_{s}^{*}({\bf r},t)\nabla\psi_{s}({\bf r},t)]$ is the KS current
density. Due to the Runge-Gross theorem and the continuity equation one can
deduce that the longitudinal part of the KS current density equals the
longitudinal part of the true current density. This need not be true
for the transverse part of the current density. However, the transverse
part of the current density does not contribute to the total current which can
then be calculated by a surface integral
\begin{equation}
I_{S}(t)=e\sum_{s=1}^{N}\int_{S}{\rm d}\sigma\; \hat{{\bf n}}\cdot{\rm Im}[
\psi_{s}^{*}({\bf r},t)\nabla\psi_{s}({\bf r},t)],
\label{curror}
\end{equation}
where $\hat{{\bf n}}$ is the unit vector perpendicular to the surface
element ${\rm d}\sigma$ and the surface $S$ is perpendicular to the longitudinal
geometry of our system. In order to propagate the KS
orbitals we need to solve the Schr\"odinger equation for a macroscopic
and non-periodic system. This goal is hopeless unless we know
the dynamics of the remote parts of the system.
We restrict ourselves to metallic electrodes. Then, the external
potential and the disturbance introduced by the device region are
screened deep inside the electrodes.
As the system size increases, the remote parts are less disturbed by the
junction and the density inside the electrodes approaches
the equilibrium bulk-density. Thus, the macroscopic size of the
electrodes leads to an enormous simplification since the initial-state
self-consistency is not disturbed far away from the constriction.

It is convenient to partition the system into three main regions:
a central region $C$ consisting of the junction and a few atomic layers of the left
and right electrodes and two regions $L$, $R$ which describe the left
and right bulk electrodes. Only the central device region $C$ will be
treated explicitly. Our scheme accounts for the full dynamical
screening in the central region. It can be further refined by taking
into account screening effects also deeper in the electrodes at the
level of linear response, with a limited increase in numerical
efforts. (These effects might be of importance in the initial transient
phase where long-range plasma oscillations in the electrodes may
occur.)
According to the above partitioning, the original
KS Hamiltonian can be written as a $3\times 3$ block matrix, and
the Schr\"odinger equation reads
\begin{equation}
i\frac{\partial}{\partial t}\left[\begin{array}{c}
\psi_{L} \\ \psi_{C} \\ \psi_{R}
\end{array}
\right]=
\left[
\begin{array}{ccc}
\mbox{\boldmath $H$}_{LL} & \mbox{\boldmath $H$}_{LC} & 0 \\
\mbox{\boldmath $H$}_{CL} & \mbox{\boldmath $H$}_{CC} & \mbox{\boldmath $H$}_{CR} \\
0 & \mbox{\boldmath $H$}_{RC} & \mbox{\boldmath $H$}_{RR}
\end{array}
\right]
\left[\begin{array}{c}
\psi_{L} \\ \psi_{C} \\ \psi_{R}
\end{array}
\right],
\label{tdse}
\end{equation}
where $\psi_{\alpha}$ is the projected wave-function onto the region $\alpha=L,R,C$.
We can solve the diffential equation for $\psi_{L}$ and $\psi_{R}$ by
introducing the retarded Green function $\mbox{\boldmath $g$}_{\alpha}$ for electrode $\alpha$, which
satisfies $\left\{
i\frac{{\rm d}}{{\rm d} t}-\mbox{\boldmath $H$}_{\alpha\a}(t)
\right\}\mbox{\boldmath $g$}_{\alpha}(t,t')=\delta(t-t')$, $\alpha=L,R$,
with boundary conditions $\mbox{\boldmath $g$}_{\alpha}(t^{+},t)=-i$ and $\mbox{\boldmath $g$}_{\alpha}(t,t^{+})=0$.
Then, we have for $\alpha=L,R$
\begin{equation}
\psi_{\alpha}(t)=i\mbox{\boldmath $g$}_{\alpha}(t,0)\psi_{\alpha}(0)+\int_{0}^{t}{\rm d}t'\mbox{\boldmath $g$}_{\alpha}(t,t')
\mbox{\boldmath $H$}_{\alpha C}\psi_{C}(t').
\label{leadeq}
\end{equation}
Using Eq. (\ref{leadeq}), the equation for $\psi_{C}$ can be written as
\begin{eqnarray}
i\frac{\partial}{\partial t}\psi_{C}(t)&=&
\mbox{\boldmath $H$}_{CC}(t)\psi_{C}(t)+\int_{0}^{t}{\rm d}t'\mbox{\boldmath $\Sigma$}(t,t')\psi_{C}(t')
\nonumber \\
&& +i\sum_{\alpha=L,R}\mbox{\boldmath $H$}_{C\alpha}
\mbox{\boldmath $g$}_{\alpha}(t,0)\psi_{\alpha}(0),
\label{evolc}
\end{eqnarray}
where $\mbox{\boldmath $\Sigma$}=\sum_{\alpha=L,R}\mbox{\boldmath $H$}_{C\alpha}\mbox{\boldmath $g$}_{\alpha}\mbox{\boldmath $H$}_{\alpha C}$ is the
self-energy which accounts for the hopping in and out of region $C$.
Thus, for any given KS orbital we can evolve its projection
onto the central region by solving Eq.~(\ref{evolc})
in region $C$. Eq.~(\ref{evolc}) has also been derived elsewhere
(for static Hamiltonians).\cite{HellumsFrensley:94}
To summarize, all the
complexity of the infinite electrode-junction-electrode system shown in
Fig.~\ref{curr1} has been reduced to the solution of an
open quantum-mechanical system (the central region $C$) with proper
time-dependent boundary conditions.

\section{Computation of extended eigenstates}
\label{secIII}

Eq.~(\ref{evolc}) of the previous Section is the central equation of our
approach to time-dependent transport. It is a reformulation of
the original time-dependent Schr\"odinger equation (\ref{tdse}) of the full
system in terms of an equation for the central (device) region only. The coupling
to the leads is taken into account by the lead Green functions $\mbox{\boldmath $g$}_{\alpha},
\alpha = L,R$. Eq.~(\ref{evolc}) has the structure of a time-dependent
Schr\"odinger equation with two extra terms: the term involving the
self energy $\mbox{\boldmath $\Sigma$}$ which we will also denote as the memory integral
since it involves the wavefunction in the central region at previous
times during the propagation.~\cite{memory}
The second term describes the injection of particles induced by a
non-vanishing projection of the initial wave-function onto the leads.

Eq.~(\ref{evolc}) is first order in time, therefore we need to specify an
initial state which is to be propagated. We want to study the time
evolution of systems perturbed out of their equilibrium ground state.
Of course, the ground state of our noninteracting system is the Slater
determinant of the occupied eigenstates of the full, extended Hamiltonian
in equilibrium, $\mbox{\boldmath $H$}(t<0) = \mbox{\boldmath $H$}^s$.
The practical question then arises how one can obtain these eigenstates
without having to deal explicitly with the extended Hamiltonian. Note
that, unlike in a bulk solid, the translational symmetry is broken in
the present device situation.

In the present Section we propose a solution of this problem which is based
on the partitioning approach used in many steady-state transport calculations
(see, {\em e.g.}, Ref. \onlinecite{Datta:95}). The retarded Green function of
the static Hamiltonian in the energy domain is determined by
\begin{equation}
\left[ (E + i \eta) \mathbf{1} - \mbox{\boldmath $H$}^s \right] \mbox{\boldmath ${\cal G}$}(E) = \mathbf{1} \; .
\label{g-energy}
\end{equation}
The Green function $\mbox{\boldmath ${\cal G}$}(E)$ of the full system can be written
in the same block structure as the Hamiltonian
\begin{equation}
\mbox{\boldmath ${\cal G}$}(E) = \left[
\begin{array}{ccc}
\mbox{\boldmath ${\cal G}$}_{LL}(E) & \mbox{\boldmath ${\cal G}$}_{LC}(E) & \mbox{\boldmath ${\cal G}$}_{LR}(E) \\
\mbox{\boldmath ${\cal G}$}_{CL}(E) & \mbox{\boldmath ${\cal G}$}_{CC}(E) & \mbox{\boldmath ${\cal G}$}_{CR}(E) \\
\mbox{\boldmath ${\cal G}$}_{RL}(E) & \mbox{\boldmath ${\cal G}$}_{RC}(E) & \mbox{\boldmath ${\cal G}$}_{RR}(E)
\end{array}
\right].
\end{equation}
Eq.~(\ref{g-energy}) can be solved for the block of the Green function
referring only to the central region
\begin{eqnarray}
\lefteqn{ \mbox{\boldmath ${\cal G}$}_{CC}(E) = } \nonumber \\
&& \!\!\!\!\! \frac{1}{ (E+ i \eta) \mathbf{1}_C - \mbox{\boldmath $H$}_{CC}^s
- \sum_{\alpha=L,R} \mbox{\boldmath $H$}_{C\alpha}^s \mbox{\boldmath $g$}_{\alpha}(E) \mbox{\boldmath $H$}_{\alpha C}^s }
\label{gcc}
\end{eqnarray}
with the retarded Green function of lead $\alpha$
\begin{equation}
\mbox{\boldmath $g$}_{\alpha}(E) = \frac{1}{(E+i \eta) \mathbf{1}_{\alpha} - \mbox{\boldmath $H$}_{\alpha\a}^s}
\label{glead}
\end{equation}
and the unit matrix $\mathbf{1}_{\alpha}$ in region $\alpha$. This Green function
enters as a central ingredient into the Fisher-Lee relation
\cite{FisherLee:81} for the calculation of the transmission function.
Through the coupling to the leads it provides for level broadening of the
isolated central part, but it also contains information on the eigenstates of
the extended system.

In order to illustrate the central idea of our method
to extract the extended eigenstates from the
Green function we consider $\mbox{\boldmath $H$}^s$ to be the discretized form of a
continuous Hamiltonian $\hat{H}^s({\bf r})$.
The continous Green function and the discretized one for uniform lattice
spacing $\Delta x$ are connected by
\begin{equation}
{\cal G}({\bf r}_i,{\bf r}_j,E) = \mbox{\boldmath ${\cal G}$}(E)_{ij}/(\Delta x)^N
\label{g-contdis}
\end{equation}
where $N=1,2,3$ is the number of spatial dimensions of the problem.
We choose the convention that a single-particle orbital
$\psi_{Ej}$ of the Hamiltonian $\hat{H}^s$ is uniquely specified by
its eigenenergy $E$ and a label $j$ for the $d_E$ degenerate orbitals
of this energy.
Using the Lehmann representation and assuming that  $\hat{H}^s$ is
invariant under time-reversal, the imaginary part of ${\cal G}$ is
\begin{equation}
-\frac{1}{\pi} {\rm Im} \,\left[ {\cal G}({\bf r},{\bf r}',E) \right] =
\sum_{E'}\delta(E-E')
\sum_{j=1}^{d_{E'}} \psi_{E'j}({\bf r})
\psi_{E'j}^*({\bf r}') \; .
\label{img-cont}
\end{equation}
Multiplying Eq.~(\ref{img-cont}) by $\psi_{El}({\bf r}')$,
integrating over ${\bf r}'$, using the orthogonality of the single particle
states and dividing by the density of states
$D(E)=\sum_{E'}\delta(E-E')d_{E'}$ we obtain
\begin{equation}
-\frac{1}{\pi D(E)} \int \! \! {\rm d}^N r \; {\rm Im} \,
\left[{\cal G}({\bf r},{\bf r}',E) \right] \psi_{El}({\bf r}') =\gamma(E) \psi_{El}({\bf r})
\label{img-eigen}
\end{equation}
with
\begin{equation}
\gamma(E)=\frac{1}{D(E)}\sum_{E'}\delta(E-E').
\end{equation}
Eq.~(\ref{img-eigen}) has the structure of an eigenvalue equation where
the energy eigenstate $\psi_{El}$ is also an eigenstate of the
integral operator $-{\rm Im}\,[{\cal G}({\bf r},{\bf r}',E)]/(\pi D(E))$ with eigenvalue
$\gamma(E)$. For a given energy $E$, this integral operator has $d_E$ different
degenerate eigenstates.

We note that Eq.~(\ref{img-cont}) is valid for {\em all}
points in space, in particular also for both ${\bf r}$ and ${\bf r}'$
representing points in the central region. In this case we know the
Green function $\mbox{\boldmath ${\cal G}$}_{CC}$ through Eqs.~(\ref{gcc}) and (\ref{g-contdis}).
Below we show that the
eigenfunctions of ${\rm Im}[\mbox{\boldmath ${\cal G}$}_{CC}]$
can be expressed as linear combination of the $\psi_{El}$.
Let us consider the matrix formed by the elements
\begin{equation}
b_{ml} = \int_C \! \! {\rm d}^N r \; \psi_{Em}^{*}({\bf r}) \psi_{El}({\bf r}) \; .
\end{equation}
where the integration is over the central region only. This matrix
is Hermitian and can be diagonalized, {\em i.e.},
\begin{equation}
\sum_{l=1}^{d_E} b_{ml} a_{jl} = \lambda_j a_{jm}
\end{equation}
with $\lambda_j$ real. Next we compute the matrix elements of the Green
function with respect to the functions
\begin{equation}
a_{Ej}({\bf r}) = \sum_{l=1}^{d_E} a_{jl}(E) \psi_{El}({\bf r}) \; .
\end{equation}
After a straightforward manipulation one finds
\begin{eqnarray}
-\frac{1}{\pi} \int_C \! \! {\rm d}^N r \; \int_C \! \! {\rm d}^N r' \;
a_{El}({\bf r}) {\rm Im} \, [{\cal G}({\bf r},{\bf r}',E)] a_{Ej}({\bf r}') \nonumber \\
= \delta_{jl} \lambda_j^2 \sum_{E'} \delta(E-E')
\end{eqnarray}
which shows explicitly that the functions $a_{Ej}({\bf r})$
diagonalize ${\rm Im} \,[\mbox{\boldmath ${\cal G}$}_{CC}]$
in the central region and that the eigenvalues are positive. Since
any linear combination of degenerate eigenstates is again an eigenstate,
diagonalizing ${\rm Im}\,[\mbox{\boldmath ${\cal G}$}_{CC}(E)]$ gives us one set of linearly independent,
degenerate eigenstates of energy $E$.

In our practical implementation described in more detail in Section
\ref{details}, we diagonalize
\begin{equation}
-\frac{1}{\pi D_{C}(E)} {\rm Im}\,[\mbox{\boldmath ${\cal G}$}_{CC}(E)]
\label{img-pract}
\end{equation}
where
\begin{equation}
D_C(E) = -\frac{1}{\pi} {\rm Tr} \, \left\{{\rm Im}\,[\mbox{\boldmath ${\cal G}$}_{CC}(E)]\right\}
\end{equation}
is the total density of states in the central region. If we use $N_g$ grid
points to describe the central region, the diagonalization in principle
gives $N_g$ eigenvectors but only a few have the physical meaning of
extended eigenstates at this energy. It is, however, very easy to identify
the physical states by looking at the eigenvalues: only few eigenvalues
(for the simple examples we studied either one or two)
are nonvanishing  and they always add up to unity.
The corresponding states are the physical
ones. All the other eigenvalues are zero (or numerically close to zero) and
the corresponding states have no physical meaning.

The procedure described above gives the correct extended eigenstates only up
to a normalization factor. When diagonalizing Eq. (\ref{img-pract})
with typical library routines one obtains eigenvectors which are normalized
to the central region. Physically this
might be incorrect. Therefore, the normalization has to be
fixed separately. In the example of Section \ref{details} we fixed the norm
by matching the wavefunction for the central region to the known form
(and normalization) of the
wavefunction in the macroscopic leads.

It should be emphasized that the procedure described here for the extraction
of eigenstates of the extended system from $\mbox{\boldmath ${\cal G}$}_{CC}(E)$ only works in practice
if $E$ is in the continuous part of the spectrum due to the sharp peak
of the delta function in the discrete part of the spectrum.

Eigenstates in the discrete part of the spectrum can be found by
considering the original Schr\"odinger equation for the full system:
\begin{equation}
\mbox{\boldmath $H$}^s \psi = E \psi.
\end{equation}
Using again the block structure of the Hamiltonian this can be transformed
into an effective Schr\"odinger equation for an {\em energy-dependent} Hamiltonian
for the central region only:
\begin{equation}
\left( \mbox{\boldmath $H$}^s_{CC} - \sum_{\alpha = L,R} \mbox{\boldmath $H$}_{C\alpha} \mbox{\boldmath $g$}_{\alpha}(E)
\mbox{\boldmath $H$}_{\alpha C} \right) \psi_C = E \psi_C.
\end{equation}
This equation has solutions only for certain values of $E$ which are
the discrete eigenenergies of the full Hamiltonian $\hat{H}^{s}$.
Therefore, one can find these states
by iteration. We have succesfully tested this idea for systems where
the analytic solutions are known. However, since the main focus of the
present work is transport where the continuum states are the essential
physical ingredient, we will not deal with the states in the discrete
spectrum for the remainder of this paper. Those states might play a role
in the description of charge-accumulation in molecular transport,
as, {\em e.g.}, in Coulomb blockade phenomena.

\section{Algorithm for time evolution}
\label{secIV}

In order to calculate the longitudinal current in an
electrode-junction-electrode system we need to propagate the
Kohn-Sham orbitals. The main difficulty stems from the macroscopic
size of the electrodes whose remote parts,ultimately, are taken
infinitely far away from the central, explicitly treated, scattering
region $C$.

The problem can be solved by imposing transparent boundary
conditions\cite{Moyer:04} at the electrode-junction interfaces.
Efficient algorithms have been
proposed for wave-packets initially {\em localized} in the scattering
region and for Hamiltonians constant in time. In this Section
we propose an algorithm well suited for
delocalized initial states, as well as for localized ones, evolving
with a time-dependent Hamiltonian.

Let $\mbox{\boldmath $H$}(t)$ be the time-dependent KS Hamiltonian.
We partition $\mbox{\boldmath $H$}(t)$ as
in Section \ref{secII}. The explicitly treated region $C$ includes
the first few atomic layers of the left and right electrodes. The boundaries of
this region are chosen in such a way that the density outside $C$
is accurately described by an equilibrium bulk density.
It is convenient to write $\mbox{\boldmath $H$}_{\alpha\a}(t)$, with $\alpha=L,R$, as the sum of a
term $\mbox{\boldmath $H$}^{s}_{\alpha\a}$ which is constant in time and another term $\mbox{\boldmath $U$}_{\alpha}(t)$
which is explicitly time-dependent, $\mbox{\boldmath $H$}_{\alpha\a}(t)=\mbox{\boldmath $H$}^{s}_{\alpha\a}+\mbox{\boldmath $U$}_{\alpha}(t)$.
In configuration space $\mbox{\boldmath $U$}_{\alpha}(t)$
is diagonal at any time $t$ since the KS potential is local in space.
Furthermore, the diagonal elements $U_{\alpha}({\bf r},t)$ are spatially
constant for metallic electrodes. Thus,
$\mbox{\boldmath $U$}_{\alpha}(t)=U_{\alpha}(t){\bf 1}_{\alpha}$ and $U_{L}(t)-U_{R}(t)$ is the
total potential drop across the central region.
We write $\mbox{\boldmath $H$}(t)=\tilde{\mbox{\boldmath $H$}}(t)+\mbox{\boldmath $U$}(t)$ with
\begin{equation}
\tilde{\mbox{\boldmath $H$}}(t)=\left[
\begin{array}{ccc}
\mbox{\boldmath $H$}^{s}_{LL} & \mbox{\boldmath $H$}_{LC} & 0 \\
\mbox{\boldmath $H$}_{CL} & \mbox{\boldmath $H$}_{CC}(t) & \mbox{\boldmath $H$}_{CR} \\
0 & \mbox{\boldmath $H$}_{RC} & \mbox{\boldmath $H$}^{s}_{RR}
\end{array}
\right],
\end{equation}
and
\begin{equation}
\mbox{\boldmath $U$}(t)=\left[
\begin{array}{ccc}
U_{L}(t){\bf 1}_{L} & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & U_{R}(t){\bf 1}_{R}
\end{array}
\right].
\end{equation}
In this way, the only term
in $\tilde{\mbox{\boldmath $H$}}(t)$ that depends on $t$ is $\mbox{\boldmath $H$}_{CC}(t)$.
For any given initial state $\psi(0)=\psi^{(0)}$ we calculate
$\psi(t_{m}=m\Delta t)=\psi^{(m)}$ by using a generalized form of
the Cayley method
\begin{widetext}
\begin{equation}
(1+i\delta
\tilde{\mbox{\boldmath $H$}}^{(m)})\frac{1+i\frac{\delta}{2}\mbox{\boldmath $U$}^{(m)}}{1-i\frac{\delta}{2}\mbox{\boldmath $U$}^{(m)}}\psi^{(m+1)}=
(1-i\delta
\tilde{\mbox{\boldmath $H$}}^{(m)})\frac{1-i\frac{\delta}{2}\mbox{\boldmath $U$}^{(m)}}{1+i\frac{\delta}{2}\mbox{\boldmath $U$}^{(m)}}\psi^{(m)},
\label{prop}
\end{equation}
\end{widetext}
with
$\tilde{\mbox{\boldmath $H$}}^{(m)}=\frac{1}{2}[\tilde{\mbox{\boldmath $H$}}(t_{m+1})+\tilde{\mbox{\boldmath $H$}}(t_{m})]$,
$\mbox{\boldmath $U$}^{(m)}=\frac{1}{2}[\mbox{\boldmath $U$}(t_{m+1})+\mbox{\boldmath $U$}(t_{m})]$ and $\delta=\Delta t/2$.
It should be noted that our propagator is
norm conserving (unitary) and
accurate to second-order in $\delta$, as is the Cayley propagator.~\cite{castro}
Denoting by $\psi_{\alpha}$ the projected wave function
onto the region $\alpha=R,L,C$, we find from Eq.~(\ref{prop})
\begin{equation}
\psi_{C}^{(m+1)}=
\frac{1-i\delta \mbox{\boldmath $H$}_{\rm eff}^{(m)}}{1+i\delta \mbox{\boldmath $H$}_{\rm eff}^{(m)}}
\psi_{C}^{(m)}
+S^{(m)}-M^{(m)}.
\end{equation}
Here, $\mbox{\boldmath $H$}_{\rm eff}^{(m)}$ is the effective Hamiltonian of the central region:
$\mbox{\boldmath $H$}_{\rm eff}^{(m)}=\mbox{\boldmath $H$}_{CC}^{(m)}-
i\delta\mbox{\boldmath $H$}_{CL}(1+i\delta\mbox{\boldmath $H$}^{s}_{LL})^{-1}\mbox{\boldmath $H$}_{LC}
-i\delta\mbox{\boldmath $H$}_{CR}(1+i\delta\mbox{\boldmath $H$}^{s}_{RR})^{-1}\mbox{\boldmath $H$}_{RC}$ with
$\mbox{\boldmath $H$}_{CC}^{(m)}=\frac{1}{2}[\mbox{\boldmath $H$}_{CC}(t_{m+1})+\mbox{\boldmath $H$}_{CC}(t_{m})]$. The source term
$S^{(m)}$ describes the injection of density into the region $C$,
while the memory term $M^{(m)}$ is responsible for the hopping in and
out of the region $C$. In terms of the propagator for the uncontacted
and undisturbed $\alpha$ electrode
\begin{equation}
\mbox{\boldmath $g$}_{\alpha}=\frac{1-i\delta\mbox{\boldmath $H$}^{s}_{\alpha\a}}{1+i\delta\mbox{\boldmath $H$}^{s}_{\alpha\a}},
\label{discgreen}
\end{equation}
the source term can be written as
\begin{equation}
S^{(m)}=-\frac{2i\delta}{1+i\delta \mbox{\boldmath $H$}_{\rm eff}^{(m)}}
\sum_{\alpha=L,R}\frac{\Lambda_{\alpha}^{(m,0)}}{u_{\alpha}^{(m)}}\mbox{\boldmath $H$}_{C\alpha}
\frac{[\mbox{\boldmath $g$}_{\alpha}]^{m}}{1+i\delta\mbox{\boldmath $H$}^{s}_{\alpha\a}}\psi_{\alpha}^{(0)},
\end{equation}
with
\begin{equation}
u_{\alpha}^{(m)}=\frac{1-i\frac{\delta}{2}U_{\alpha}^{(m)}}{1+i\frac{\delta}{2}U_{\alpha}^{(m)}}
\quad {\rm and}\quad
\Lambda_{\alpha}^{(m,k)}=\prod_{j=k}^{m}[u_{\alpha}^{(j)}]^{2}.
\end{equation}
For a wave packet initially localized in $C$ the projection onto the left
and right electrode $\psi_{\alpha}^{(0)}$ vanishes and $S^{(m)}=0$ for any
$m$, as it should be. The memory term is more
complicated and reads
\begin{widetext}
\begin{equation}
M^{(m)}=-\frac{\delta^{2}}{1+i\delta \mbox{\boldmath $H$}_{\rm eff}^{(m)}}\sum_{\alpha=L,R}\sum_{k=0}^{m-1}
\frac{\Lambda_{\alpha}^{(m,k)}}{u_{\alpha}^{(m)}u_{\alpha}^{(k)}}
[\mbox{\boldmath $Q$}_{\alpha}^{(m-k)}+\mbox{\boldmath $Q$}_{\alpha}^{(m-k-1)}]
(\psi_{C}^{(k+1)}+\psi_{C}^{(k)})
\end{equation}
\end{widetext}
where
\begin{equation}
\mbox{\boldmath $Q$}_{\alpha}^{(m)}=\mbox{\boldmath $H$}_{C\alpha}\frac{[\mbox{\boldmath $g$}_{\alpha}]^{m}}{1+i\delta\mbox{\boldmath $H$}^{s}_{\alpha\a}}\mbox{\boldmath $H$}_{\alpha C}.
\label{Qdef}
\end{equation}
The quantities $\mbox{\boldmath $Q$}_{\alpha}^{(m)}$ depend on the geometry of the
system and are independent of the initial state $\psi^{(0)}$.

\begin{figure}[b]
\includegraphics[scale=0.55]{fig2.eps}
\caption{\label{semi} Schematic sketch of an electrode-junction-electrode
system with semiperiodic electrodes.}
\end{figure}

Below we propose a recursive scheme to calculate the $\mbox{\boldmath $Q$}_{\alpha}^{(m)}$'s for
those system geometries having semiperiodic electrodes along the
longitudinal direction, see Fig. \ref{semi}. In this case $\mbox{\boldmath $H$}^{s}_{\alpha\a}$
has a tridiagonal block form
\begin{equation}
\mbox{\boldmath $H$}^{s}_{\alpha\a}=\left[
\begin{array}{cccc}
    \mbox{\boldmath $h$}_{\alpha} & \mbox{\boldmath $V$}_{\alpha} & 0 & \ldots \\
    \mbox{\boldmath $V$}_{\alpha} & \mbox{\boldmath $h$}_{\alpha} & \mbox{\boldmath $V$}_{\alpha} & \ldots \\
    0 & \mbox{\boldmath $V$}_{\alpha} & \mbox{\boldmath $h$}_{\alpha} & \ldots \\
    \ldots & \ldots & \ldots & \ldots
\end{array}
\right],
\end{equation}
where $\mbox{\boldmath $h$}_{\alpha}$ describes a convenient cell and $\mbox{\boldmath $V$}_{\alpha}$ is the
hopping Hamiltonian between two nearest neighbor cells.
Without loss of generality we assume that both $\mbox{\boldmath $h$}_{\alpha}$ and $\mbox{\boldmath $V$}_{\alpha}$
are square matrices of dimension $N_{\alpha}\times N_{\alpha}$.
Taking into account that the central region
contains the first few cells of the left and right electrodes, the
matrix $\mbox{\boldmath $Q$}_{\alpha}^{(m)}$ has the following structure
\begin{equation}
\mbox{\boldmath $Q$}_{L}^{(m)}=\left[
\begin{array}{ccc}
\mbox{\boldmath $q$}_{L}^{(m)} & 0 & 0 \\
0 & 0  & 0 \\
0 & 0 & 0
\end{array}
\right],\quad
\mbox{\boldmath $Q$}_{R}^{(m)}=\left[
\begin{array}{ccc}
0 & 0 & 0 \\
0 & 0  & 0 \\
0 & 0 & \mbox{\boldmath $q$}_{R}^{(m)}
\end{array}
\right].
\end{equation}
The $\mbox{\boldmath $q$}_{\alpha}^{(m)}$'s are square matrices of dimension $N_{\alpha}\times
N_{\alpha}$ and are given by
\begin{equation}
\mbox{\boldmath $q$}_{\alpha}^{(m)}=\mbox{\boldmath $V$}_{\alpha}\left[
\frac{[\mbox{\boldmath $g$}_{\alpha}]^{m}}{1+i\delta\mbox{\boldmath $H$}_{\alpha\a}}
\right]_{1,1}\mbox{\boldmath $V$}_{\alpha},
\end{equation}
where the subscript $(1,1)$ denotes the first diagonal block of the
matrix in the square brackets. We introduce the generating matrix function
\begin{equation}
\mbox{\boldmath $q$}_{\alpha}(x,y)\equiv\mbox{\boldmath $V$}_{\alpha}\left[
\frac{1}{x+iy\delta\mbox{\boldmath $H$}_{\alpha\a}}
\right]_{1,1}\mbox{\boldmath $V$}_{\alpha},
\label{q1def}
\end{equation}
which can also be expressed in terms of continued matrix fractions (see Appendix
\ref{app})
\begin{equation}
\mbox{\boldmath $q$}_{\alpha}(x,y)=\mbox{\boldmath $V$}_{\alpha}\frac{1}{x+iy\delta\mbox{\boldmath $h$}_{\alpha}+y^{2}\delta^{2}\mbox{\boldmath $q$}_{\alpha}(x,y)}\mbox{\boldmath $V$}_{\alpha}.
\label{qalpha}
\end{equation}
The $\mbox{\boldmath $q$}_{\alpha}^{(m)}$'s can be obtained from
\begin{equation}
\mbox{\boldmath $q$}_{\alpha}^{(m)}=\frac{1}{m!}\left.\left[-\frac{\partial}{\partial x}+\frac{\partial }{\partial y}\right]^{m}
\mbox{\boldmath $q$}_{\alpha}(x,y)\right|_{x=y=1}.
\label{gen}
\end{equation}

From Eqs.~(\ref{gen}) and (\ref{qalpha}) one can build up a recursive
scheme. Let us define
$\mbox{\boldmath $p$}_{\alpha}^{-1}(x,y)=x+iy\delta\mbox{\boldmath $h$}_{\alpha}+y^{2}\delta^{2}\mbox{\boldmath $q$}_{\alpha}(x,y)$ and
$\mbox{\boldmath $p$}_{\alpha}^{(m)}=\frac{1}{m!}[-\frac{\partial}{\partial x}+\frac{\partial}{\partial y}]
\mbox{\boldmath $p$}_{\alpha}(x,y)|_{x=y=1}$. Then, by definition,
\begin{equation}
\mbox{\boldmath $q$}_{\alpha}^{(m)}=\mbox{\boldmath $V$}_{\alpha}\mbox{\boldmath $p$}_{\alpha}^{(m)}\mbox{\boldmath $V$}_{\alpha}.
\end{equation}
Using the identity $\frac{1}{m!}[-\frac{\partial}{\partial x}+\frac{\partial}{\partial y}]^{m}
\mbox{\boldmath $p$}_{\alpha}(x,y)\mbox{\boldmath $p$}_{\alpha}^{-1}(x,y)=0$, one finds
\begin{eqnarray}
(1+i\delta\mbox{\boldmath $h$}_{\alpha})\mbox{\boldmath $p$}_{\alpha}^{(m)}=(1-i\delta\mbox{\boldmath $h$}_{\alpha})\mbox{\boldmath $p$}_{\alpha}^{(m-1)}
\quad\quad\quad\quad\quad
\nonumber \\ -
\delta^{2}\sum_{k=0}^{m}(\mbox{\boldmath $q$}_{\alpha}^{(k)}+2\mbox{\boldmath $q$}_{\alpha}^{(k-1)}+\mbox{\boldmath $q$}_{\alpha}^{(k-2)})
\mbox{\boldmath $p$}_{\alpha}^{(m-k)}
\label{recur}
\end{eqnarray}
with $\mbox{\boldmath $p$}_{\alpha}^{(m)}=\mbox{\boldmath $q$}_{\alpha}^{(m)}=0$ for $m<0$. Once
$\mbox{\boldmath $q$}_{\alpha}^{(0)}$ has been obtained by solving Eq.~(\ref{qalpha})
with $x=y=1$, we can calculate $\mbox{\boldmath $p$}_{\alpha}^{(0)}=
[1+i\delta\mbox{\boldmath $h$}_{\alpha}+\delta^{2}\mbox{\boldmath $q$}_{\alpha}^{(0)}]^{-1}$. Afterwards, we can use
Eq.~(\ref{recur}) with $\mbox{\boldmath $q$}_{\alpha}^{(1)}=\mbox{\boldmath $V$}_{\alpha}\mbox{\boldmath $p$}_{\alpha}^{(1)}\mbox{\boldmath $V$}_{\alpha}$
to calculate $\mbox{\boldmath $p$}_{\alpha}^{(1)}$ and hence $\mbox{\boldmath $q$}_{\alpha}^{(1)}$ and so on and
so forth.

This concludes the description of our algorithm for the propagation of the
time-dependent Schr\"odinger equation for extended systems. It is worth
mentioning an additional complication here which arises for the propagation
of a time-dependent Kohn-Sham equation. This complication stems from the fact
that in order to compute $\psi_{C}^{(m+1)}$ at time step $m+1$ one needs to
know the time-dependent KS potential at the same time step which, via the
Hartree and exchange-correlation potentials, depends on the yet unknown
orbitals $\psi_{C}^{(m+1)}$. Of course, the solution is to use a
predictor-corrector approach: in the first step one approximates
$\mbox{\boldmath $H$}_{CC}^{(m)}$ by $\mbox{\boldmath $H$}_{CC}(t_m)$, computes new orbitals
$\tilde{\psi}_{C}^{(m+1)}$ and from those obtains an improved approximation
for $\mbox{\boldmath $H$}_{CC}^{(m)}$.

\section{Implementation Details}
\label{details}

All the methodological discussion above
 is general and can be applied to general device
configurations as long as they can be mapped into a longitudinal-like geometry
as described in Fig.~\ref{semi}.  In order to demonstrate the feasibility
of the scheme described in the previous
Sections we have implemented it for one-dimensional model systems. The
extension to real molecular-device configurations is presently under
development.~\cite{wirtz}
We have used a simple three-point discretization for
the second derivative
\begin{equation}
\frac{{\rm  d}^2}{{\rm d} x^2} \psi(x) \vert_{x=x_i} \approx
\frac{1}{(\Delta x)^2} \left( \psi(x_{i+1}) - 2 \psi(x_i) +
\psi(x_{i-1}) \right)
\end{equation}
with equidistant grid points $x_i$, $i=1,\dots,N_g$ and spacing $\Delta x$.
Within this approximation matrices of the form $\mbox{\boldmath $H$}_{C\alpha} \mbox{\boldmath $M$} \mbox{\boldmath $H$}_{\alpha C}$
which are $N_g \times N_g$ matrices and appear, {\em e.g.}, in Eq.~(\ref{gcc})
or (\ref{Qdef}), have only one nonvanishing matrix element. For $\alpha=L$ this
is the $(1,1)$ element, for $\alpha=R$ it is the $(N_g,N_g)$ element.

In order to proceed we have to specify the nature of the leads and therefore
the lead Green function. Here we choose the simplest case of
semi-infinite, uniform leads at constant potential $U_{\alpha 0}$. In this case,
the Green function (\ref{glead}) in the energy domain can be given in
closed form:
\begin{eqnarray}
\lefteqn{
[\mbox{\boldmath $g$}_{\alpha}(E)]_{kl} = - \frac{i \Delta x}{\sqrt{2 \tilde{E}_{\alpha}}}
\exp\left\{i \sqrt{2 \tilde{E}_{\alpha}} | x_k - x_l|\right\}} \nonumber \\
&+& \frac{i \Delta x}{\sqrt{2 \tilde{E}_{\alpha}}}
\exp\left\{i \sqrt{2 \tilde{E}_{\alpha}}( | x_k - x_{\a0}| + |x_l -
x_{\a0}|)\right\}
\end{eqnarray}
with $\tilde{E}_{\alpha} = E-U_{\alpha 0}$. The abscissa $x_{\alpha 0}$
is the position of the
interface between the lead and the device region and $x_k = x_{\alpha 0}
\pm k \Delta x$, where the plus sign applies for $\alpha=R$ and the minus sign
for $\alpha=L$.

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{fig3up.eps}\\[1.0cm]
\includegraphics[width=0.9\columnwidth]{fig3low.eps}\\[0.5cm]
\caption{\label{ext-states} Continuum states of square potential barrier
at different energies with leads at zero potential. Upper panel:
eigenstates for $\varepsilon=0.45$ a.u., just below the barrier height of $0.5$ a.u..
Lower panel: eigenstates at $\varepsilon=0.6$ a.u..}
\end{figure}

The results of the procedure for calculating extended eigenstates as described
in Section \ref{secIII} is illustrated in Fig.~\ref{ext-states} for a
square potential barrier with zero potential in both leads.
In the upper panel we
have the square modulus of eigenstates at an energy below the barrier height
while in the lower panel eigenstates with energy higher than the barrier
are shown. The states result from diagonalization of Eq.~(\ref{img-pract}).
In order to obtain the normalization constant we compute the logarithmic
derivative at the boundary of the central region numerically and match it
to the analytic form in the lead to obtain the phase shift $\delta_{\alpha}$:
\begin{equation}
\frac{1}{2} \frac{{\rm d}^2}{{\rm d} x^2} \ln( |\psi(x)|^2 )
\bigg\vert_{x=x_{\alpha 0}} = q \cot(\delta_{\alpha})
\end{equation}
where $q= \sqrt{2 \tilde{E}_{\alpha}}$. Knowing the phase shift we can rescale the
wavefunction such that it matches with the analytic form
$\sin(q (x-x_{\alpha 0}) + \delta_{\alpha})$ at the interface. Of course, this form
of the extended states only applies for $\tilde{E}_{\alpha} >0$ but as long as
$E$ is in the continuous part of the spectrum, it is correct at least for
one of the leads. This is sufficient to determine the normalization. The
states obtained numerically with this procedure coincide with the known
analytical results.

We then implemented the propagation scheme presented in
the previous Section. Within  our three-point approximation,
$\mbox{\boldmath $h$}_{\alpha}$, $\mbox{\boldmath $V$}_{\alpha}$ and $\mbox{\boldmath $q$}_{\alpha}$ are
$1\times 1$ matrices. The equation for $q_{\alpha}^{(0)}$ [see
Eqs.~(\ref{qalpha}) and (\ref{gen})] becomes a simple quadratic equation
which can be solved explicitly
\begin{equation}
q_{\alpha}^{(0)} = \frac{- (1 + i \delta h_{\alpha})
+ \sqrt{(1 + i \delta h_{\alpha}) ^2 + 4 (\delta V_{\alpha})^2}}
{2 \delta^2}  \, .
\label{q0}
\end{equation}
Although the quadratic equation has two solutions, the above choice for
$q_{\alpha}^{(0)}$ is dictated by the fact that the Taylor expansions for small
$\delta$ of Eqs.~(\ref{q0}) and (\ref{qalpha}) have to coincide.
Using this result we then
solved the iterative scheme to obtain the $q_{\alpha}^{(m)}$ for $m \geq 1$.

As a first check on the propagation method we propagated a Gaussian
wavepacket which, at initial time $t=0$, is completely localized in the
central device region. (The source terms $S^{(m)}$ then vanish identically).
As can be seen in Fig.~\ref{wavepacket}, the wavepacket correctly propagates
through the boundaries without any spurious reflections.

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{fig4.eps}
\caption{\label{wavepacket} Time evolution of a Gaussian wavepacket with
initial width 1.0 a.u. and initial momentum 0.5 a.u. for various
propagation times. The transparent boundary conditions allow the wavepacket
to pass the propagation region without spurious reflections at the boundaries.}
\end{figure}

For the propagation of the extended initial states (the eigenstates of the
unperturbed system) we also need to implement the source terms $S^{(m)}$. In
the following we assume that the left and right leads are at the same
potential initially so that the analytic form of the initial states
is in both leads given by $\sin(q(x-x_{\alpha 0}) + \delta_{\alpha}) =
\cos(\delta_{\alpha}) \sin(q(x-x_{\alpha 0})) + \sin(\delta_{\alpha})
\cos(q(x-x_{\alpha 0}))$. The propagation of the term proportional to
$\sin(q(x-x_{\alpha 0}))$ is trivial since this is an eigenstate of the lead
Hamiltonian with energy $\varepsilon_q=q^2/2$.
Therefore, if - in discretized form - $\psi_{R,1}^{(0)} =
(\sin(q \Delta x), \sin(2 q \Delta x), \ldots )^T$ we obtain:
\begin{equation}
\mbox{\boldmath $H$}_{CR}\frac{[\mbox{\boldmath $g$}_R]^m}{1 + i \delta \mbox{\boldmath $H$}_{RR}^s} \psi_{R,1}^{(0)} =
V_{R} \frac{(1 - i \delta \varepsilon_q)^m}
{(1 + i \delta \varepsilon_q)^{m+1}}\, \mbox{\boldmath $e$}_R
\end{equation}
and similarly for the left lead. Here,
$\mbox{\boldmath $H$}_{RR}^s$ is the static part of the right-lead Hamiltonian,
$\mbox{\boldmath $g$}_R$ the corresponding Green function according to definition
(\ref{discgreen}) and $\mbox{\boldmath $e$}_R=(0,\ldots,0,1)^{T}$ is a
unit vector.

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{fig5.eps}
\caption{\label{examp1} Time evolution of the current for a system where
initially the potential is zero in the leads and the propagation
region. At $t=0$, a constant bias with opposite sign in the left and right
leads is switched on, $U=U_L=-U_R$ (values in atomic units).
The propagation region extends from
$x=-6$ to $x=+6$ a.u.. The Fermi energy of the initial state is $\varepsilon_F=0.3$
a.u.. The current in the center of the propagation region is shown.}
\end{figure}

The propagation of the part proportional to $\cos(q(x-x_{\alpha 0}))$ is
more complicated since this is not an eigenstate of the lead Hamiltonian.
We define the function $R(x,y)$ from
\begin{equation}
R(x,y) \mbox{\boldmath $e$}_R = \mbox{\boldmath $H$}_{CR} \frac{1}{x + i y \delta \mbox{\boldmath $H$}_{RR}^s} \psi_{R,2}^{(0)}
\; ,
\end{equation}
where $\psi_{R,2}^{(0)} = (\cos(q \Delta x), \cos(2 q \Delta x), \ldots )^T$.
Introducing the tridiagonal matrix
\begin{equation}
\mbox{\boldmath $O$} = \left(
\begin{array}{ccccc}
0 & 1 & 0 & 0 & \cdots \\
1 & 0 & 1 & 0 & \cdots \\
0 & 1 & 0 & 1 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots
\end{array} \right)
\end{equation}
and using
\begin{equation}
V_R \mbox{\boldmath $O$} \psi_{R,2}^{(0)} = (\varepsilon_q - h_R) \psi_{R,2}^{(0)} - V_R
(1,0,0,\ldots)^T
\end{equation}
one arrives at
\begin{equation}
R(x,y) = \frac{1}{x + i y \delta \varepsilon_q} \left[ V_R \cos(q \Delta x) +
i y \delta  q_R(x,y) \right],
\end{equation}
where $q_R(x,y)$ is given by Eq.~(\ref{q1def}). Defining
\begin{equation}
R^{(m)} = \frac{1}{m!} \left(- \frac{\partial}{\partial x} +
\frac{\partial}{\partial y} \right)^m R(x,y) \big\vert_{x=y=1}
\end{equation}
one finds
\begin{eqnarray}
\lefteqn{
R^{(m)} = \frac{(1- i \delta \varepsilon_q)^m}{(1 + i \delta \varepsilon_q)^{m+1}} V_R
\cos(q \Delta x)} \nonumber \\
&& + i \delta  \sum_{k=0}^{m} \frac{(1- i \delta \varepsilon_q)^{m-k}}
{(1+ i \delta \varepsilon_q)^{m+1-k}} \left(q_R^{(k)} + q_R^{(k-1)} \right)
\end{eqnarray}
and finally
\begin{equation}
\mbox{\boldmath $H$}_{CR}\frac{[\mbox{\boldmath $g$}_R]^m}{1 + i \delta \mbox{\boldmath $H$}_{RR}^s} \psi_{R,2}^{(0)} =
R^{(m)} \mbox{\boldmath $e$}_R.
\end{equation}
One may proceeds along the same lines for extracting the left
component.

To test our implementation we have propagated eigenstates of the extended
system. As expected, these states just pick up an exponential phase factor
$\exp(-i E t)$ during the propagation.

\section{Examples}
\label{examples}

We are now in a position to apply our algorithm to the calculation of
time-dependent currents in one-dimensional model systems. The systems
are initially in thermodynamic equilibrium. At time $t=0$,
a bias is switched on in the electrodes.

As a first example we considered a system where the electrostatic potential
vanishes identically both in the left and right leads as well as in the
central region which is explicitly propagated. Initially, all single particle
levels are occupied up to the Fermi energy $\varepsilon_F$. At $t=0$ a
constant bias is switched on in the leads and the time-evolution of the
system is calculated. We chose the bias in the right lead as the negative
of the bias in the left lead, $U_R = - U_L$. The current is
calculated from Eq. (\ref{curror}):
\begin{eqnarray}
\lefteqn{
I(x,t) = 2 \int_{-k_F}^{k_F} \frac{{\rm d} k}{2 \pi} \;\;
{\rm Im}\, \left( \psi_k^*(x,t) \frac{{\rm d}}{{\rm d} x} \psi_k(x,t) \right)} \nonumber \\
&=& 2 \int_{0}^{k_F} \frac{{\rm d} k}{2 \pi} \;\;
{\rm Im}\, \left( \psi_k^* \frac{{\rm d}}{{\rm d} x} \psi_k +
\psi_{-k}^* \frac{{\rm d}}{{\rm d} x} \psi_{-k} \right)
\label{currdens}
\end{eqnarray}
where the prefactor $2$ comes from spin and $k_F=\sqrt{2 \varepsilon_F}$ is the
Fermi wavevector.

The numerical parameters are as follows: the Fermi energy is
$\varepsilon_F=0.3$ a.u., the bias is $U_L=-U_R=0.05,0.15,0.25$ a.u., the central region
extends from $x=-6$ to $x=+6$ a.u. with equidistant grid points with spacing
$\Delta x=0.03$ a.u.. The $k$-integral in Eq.~(\ref{currdens}) is discretized
with 100 $k$-points which amounts to a propagation
of 200 states. The time step for the propagation was
$\Delta t= 10^{-2}$ a.u.

In Fig.~\ref{examp1} we have plotted the current densities at $x=0$ as a
function of time for different values of the applied bias.
As a first feature we notice that a steady state is achieved,
in agreement with the results of Ref. \onlinecite{saprb04}.
The corresponding steady-state current $I$ can be calculated from the
Landauer formula. For the present geometry this leads to the steady current
\begin{widetext}
\begin{equation}
I=8e\int_{\max(U_{L},U_{R})}\frac{d\omega}{2\pi}[f(\omega-U_{L})-f(\omega-U_{R})]
\frac{\sqrt{\omega-U_{L}}\sqrt{\omega-U_{R}}}
{\left[\sqrt{\omega-U_{L}}+\sqrt{\omega-U_{R}}\right]^{2}+U_{L}U_{R}
\left[\frac{\sin(l\sqrt{2\omega})}{\sqrt{\omega}}\right]^{2}},
\label{land}
\end{equation}
\end{widetext}
where $l$ is the width of the central region.
From Eq. (\ref{land}) with $l=12$ a.u. and $U_L=-U_R$, the
numerical values for the steady-state currents
are $0.0316$ a.u. ($U_{L}=0.05$ a.u.), $0.0883$ a.u. ($U_{L}=0.15$ a.u.) and
$0.0828$ a.u. ($U_{L}=0.25$ a.u.).
We see that our algorithm yields the same answers. Second,
we notice that the onset of the current is delayed in relation to the
switching time $t=0$. This is easily explained by the fact that the
perturbation at $t=0$ happens in the leads only, {\em e.g.}, for
$|x|> 6$ a.u.,
while we plot the current at $x=0$. In other words, we see the delay time
needed for the perturbation to propagate from the leads to the center of
our device region. We also note that the higher the bias the more the
current overshoots its steady-state value for small times after switching
on the bias. Finally it is worth mentioning that increasing the bias
not necessarily leads to a larger steady-state current.

In the second example we studied a double square potential barrier with
electrostatic potential $V(x) = 0.5$ a.u. for $5$ a.u. $\leq |x| \leq 6$ a.u.
and zero otherwise. This time we switch on a constant bias in the left lead
only, {\em i.e.}, $U_R=0$. The Fermi energy for the initial state is
$\varepsilon_F=0.3$ a.u.. The central region extends from $x=-6$ to $x=+6$ a.u. with
a lattice spacing of $\Delta x=0.03$ a.u.. Again, we use 100 different
$k$-values to compute the current and a time step of
$\Delta t = 10^{-2}$ a.u..

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{fig6.eps}
\caption{\label{examp2} Time evolution of the current through a double
square potential barrier in response to an applied constant bias (given in
atomic units) in the left
lead. The potential is given by $V(x) = 0.5$ a.u. for $5 \leq |x|\leq 6$ a.u.
and zero otherwise, the propagation region extends from $x=-6$ to $x=+6$ a.u..
The Fermi energy of the initial state is $\varepsilon_F=0.3$ a.u.. The current in the
center of the structure is shown.}
\end{figure}

In Fig.~\ref{examp2} we plot the current at
$x=0$ as a function of time for several values of the applied bias $U=U_L$.
Again, a steady state is achieved for all values of $U$.
As discussed in Fig.~\ref{examp1} the transient current can exceed the
steady current; the higher the applied voltage the larger is this excess
current and the shorter is the time when it reaches its maximum.
Furthermore, the oscillatory evolution towards the steady current solution
depends on the bias. For higher bias the frequency of the transient
oscillations increases. For small bias the electrons at the bottom of the
band are not disturbed and the transient process is
exponentially short. On the other hand, for a bias close to the Fermi
energy the transient process decays as a power law, due to
the band edge singularity. As pointed out in Ref. \onlinecite{saprb04}, for
non-interacting electrons the steady-state current develops by means
of a pure dephasing mechanism. In our examples the transient process occurs
in a femtosecond time-scale, which is much shorter than the
relaxation time due to electron-phonon interactions.

In Fig.~\ref{examp2norm} we plot the time evolution of the total number of
electrons in the device region for the same values of $U_L$. We see that as
a result of the bias a quite substantial amount of charge is added to the
device region. This result has important implications when simulating
the transport through an interacting system as the effective (dynamical)
electronic screening is modified due not only to the external field but also
to the accumulation of charge state in the molecular device. This illustrates
that linear response might not be an appropriate tool to tackle the dynamical
response and that we will need to resort to a full time-propagation
approach as the one of the present work.
Here we emphasize that all our calculations are done
without taking into account the electron-electron interaction. If we had
done a similar calculation with the interaction incorporated in a
time-dependent Hartree or time-dependent DFT framework we would expect the
amount of excess charge to be reduced significantly as compared to
Fig.~\ref{examp2norm}.

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{fig7.eps}
\caption{\label{examp2norm} Time evolution of the total number of electrons
in the region $|x|\leq 6$ for the double square potential barrier of
Fig.~\ref{examp2}. }
\end{figure}

In Fig.~\ref{examp2a} we show time-dependent currents for the same double
barrier as in Fig.~\ref{examp2} for two different ways of applying the
bias in the left lead: in one case the constant bias $U_0$ is switched
on suddenly at $t=0$ (as in Fig.~\ref{examp2}), in the other case the
constant $U_0$ is achieved with a smooth switching $U(t)=U_0
\sin^2(\omega t)$ for $0<t<\pi/(2 \omega)$. As expected and in agreement
with the results of Ref.~\onlinecite{saprb04}, the same steady state
is achieved after the initial transient time. However, the transient
current clearly depends on how the bias is switched on.

In the final example we address the simulation of ac-transport.
We computed the current for a single square potential
barrier with $V(x) = 0.6$ for $|x|<6$ and zero otherwise. Here we
applied a time-dependent bias of the form $U_L(t) = U_{0} \sin(\omega t)$
to the left lead. The right lead remains on zero bias. The numerical
parameters are: Fermi energy $\varepsilon_F=0.5$ a.u., device region from $x=-6$ to
$x=+6$ a.u. with lattice spacing $\Delta x=0.03$ a.u.. The number of
$k$-points is 100 and the time step is $\Delta t=  10^{-2}$ a.u..
In Fig.~\ref{examp3} we plot the current at $x=0$ as a function of time for
different values of the parameter $U_{0}=0.1,0.2,0.3$ a.u. The frequency was
chosen as $\omega=1.0$ a.u. in both cases. Again, as for the dc-calculation
discussed above, we get a transient that overshoots the average current
flowing through the constriction; again, this excess current is larger the
higher the applied voltage. Also, after the transient we obtain a current
through the system with the same period as the applied bias. Note, however,
that (especially for the large bias), the current is not a simple harmonic
as the applied ac bias.

\section{Summary and outlook}
\label{conc}

In the present work we have presented a formally
rigorous approach towards the description of charge transport using
an open-boundary scheme within TDDFT. We have implemented a
specific time-propagation scheme that incorporates transparent
boundaries at  the device/lead interface in a natural way.
In order to have a clear
definition of a device region in Fig.~\ref{system} we assumed
that an applied bias can be described by adding a spatially constant
potential shift in the macroscopic part of the leads. This implies an
effective ``metallic screening'' of the constriction. The screening limits
the spatial extent of the induced density created by
the bias potential or the external field
to the central region. Our treatment can be further refined to
include dynamical screening deep inside the electrodes on the level of
linear response, which might be of importance for the initial
transient.
Our time-dependent scheme allows to extract both ac and dc I/V
device characteristics and it is ideally suited to describe
external field (photon) assisted processes.

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{fig8.eps}
\caption{\label{examp2a} Time evolution of the current for a double
square potential barrier when the bias is switched on in two different
manners: in one case, the bias $U_0$ is suddenly switched on at $t=0$
while in the other case the same bias is achieved with a smooth switching
$U(t) = U_0 \sin^2(\omega t)$ for $0<t<\pi/(2 \omega)$. The parameters
for the double barrier and the other numerical parameters are the same as
the ones used in Fig.~\ref{examp2}. $U_0$ and $\omega$ given in atomic units.}
\end{figure}

In order to illustrate the performance
of the method, we have implemented it for
one-dimensional models and we have shown:
i) How to extract the proper initial extended states to be propagated.
ii) How to incorporate perfect transparent boundaries for the time propagation.
iii) A steady-state current is always reached upon application
of a dc bias. The transient process occurs on a time-scale much
shorter than the relaxation time due to electron-phonon interaction.
In the case of systems without any source of dissipation it
is known that the steady-current is
independent of the history of the process.~\cite{saprb04}
We have explicitly demonstrated this history independence for two different
switching processes of the external bias. However, if we allow for
dissipation either through electron-electron or electron-phonon
interactions, the current versus voltage characteristics might depend
on the history. For instance, hysteresis loops due to different
transient electronic/geometrical device configurations are possible.
This effect will be more dramatic in the case of ac-driving fields of high
frequencies where the system may not have enough time to respond to
the perturbation.
iv) A periodic  ac current is reached upon perturbation with a
monochromatic field.

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{fig9.eps}
\caption{\label{examp3} Time evolution of the current for a square potential
barrier in response to a time-dependent, harmonic bias in the left lead,
$U_L(t) = U_0 \sin(\omega t)$ for different amplitudes $U_0$ (values in a.u.)
and frequency $\omega=1.0$ a.u.. The potential is given by $V(x)=0.6$ a.u. for
$|x|\leq 6.0$ a.u. and zero otherwise. The propagation region extends from
$x=-6$ to $x=+6$ a.u.. The Fermi energy of the initial state is
$\varepsilon_F=0.5$ a.u.. The current at $x=0$ is shown. }
\end{figure}

Previous work on time-dependent quantum transport mainly uses the idea of
Caroli.\cite{caroli1,caroli2}
This approach is at the core of the Landauer derivation and has the problem
of using different chemical potentials in different parts of the system.
This implies that the initial state is not a ground
state of the entire, contacted system. Furthermore, the
time-dependent perturbation is a tunneling Hamiltonian which is nonlocal
in space. Therefore, it cannot be combined with TDDFT since
there the time-dependent potential is local.

Here, we have used the scheme of Cini \cite{cini} which can be combined
with TDDFT: We start the calculation from a well-defined thermodynamic
equilibrium configuration, therefore the scheme is thermodynamically
consistent. Then we apply an external potential that in general is time
dependent. By virtue
of the Runge-Gross theorem, the time-evolution of this quantum system is
completely determined by the knowledge of the time-dependent density.
In the steady state regime the occupation of left and right moving carriers
is dictated by the time-dependent Schr\"odinger equation.

Another thermodynamically consistent scheme has been put forward by
Kamenev and Kohn. \cite{kamenev:01} They used
the microscopic quantum-kinetic formulation of conductivity
and worked with the Kubo formula~\cite{kubo}
in the linear time-dependent Hartree regime.
They consider a closed system (ring) where the current in the system results
from an external driving vector potential. This approach also overcomes
the problem of having two chemical potentials. They have shown that
it is possible to recover the Landauer result, {\em i.e.}, the universal quantum
of conductance $2e^2/h$ independent of the length and material.
Since the Kamenev-Kohn approach uses a vector potential rather than a
scalar potential, time-dependent current DFT rather than TDDFT would be the
natural density-functional extension.

Most theoretical approaches to transport in molecular electronic devices adopt
open boundary conditions and assume that transport is ballistic,
{\em i.e.}, under steady state conditions inelastic collisions are absent from
a molecular structure and its
contacts.~\cite{Landauer:57,Buettiker:86,book,Datta:95}
Dissipation occurs only in the idealized macroscopic reservoirs
connected by leads to the molecular device. This central role of the reservoirs
in the process of dissipation is a valid approximation, particularly
when the applied bias is small and a device operates in a linear regime.
When inelastic scattering dominates this picture is not applicable.
In particular, experiments \cite{kushmerick,smit} indicate that inelastic
scattering with lattice vibrations is present at sufficiently large bias,
causing local heating of contacts and molecular devices.
Energy transfer to the lattice may also cause atomic migration and
result in dramatic changes in the device characteristics (also
they may give phonon replica structure in the measured conductance).
The modelling of a many-electron system out of equilibrium coupled
to lattice vibrations is a real theoretical challenge.~\cite{car}

Electron correlations are also important in molecular conductors,
for example, Coulomb blockade effects dominate the transport in quantum dots.
Short-range electron correlations seems to be relevant in order to get
quantitative description of I/V characteristics in molecular
constrictions.~\cite{delaney:04,everts:04,ferreti}
In particular it is commonly assumed that the energy scales for
electron-electron and electron-phonon interactions are different and
could be treated separately. However, the metallic screening of the
electrodes considerably reduces the Coulomb-charging energy (from eV to
meV scale). In this regime the energy scale for the two
interactions merge and they need to be treated on the same footing
posing some additional theoretical challenge.

Other approaches put forward in the literature directly look
for a homogeneous current-carrying state either based
on a a maximum entropy principle \cite{godby:03} or by a imposing
the current through Lagrange-multipliers.\cite{kosov:04} In those
approaches it is implicitly assumed
 that the origin of the homogeneous current is
independent whether it is introduced by reservoirs or by
external fields. This is indeed the case for independent electrons but once
dissipation is built in the system might exhibit a dependence on the history
of the applied bias,({\em e.g.}, possible appearance of
hysteresis loop in the current versus voltage characteristics).

It is clear that the quality of the TDDFT functionals is of crucial
importance. In particular, exchange and correlation functionals for
the non-equilibrium situation are required.
Time-dependent linear response theory for dc-steady state
has been implemented in Ref.~\onlinecite{baer:04}
within TDLDA assuming jellium-like electrodes (mimicked by complex
absorbing/emitting potentials). It has been shown that the
dc-conductance changes considerably from the standard Landauer value.
Therefore, a systematic study of the TDDFT
functionals themselves is needed. A step beyond
standard adiabatic approximations and exchange-only potentials
is to resort to many-body schemes as recently
done for the characterization of optical properties of semiconductors and
insulators.~\cite{exc} Another path is to explore in depth
the fact that the true exchange-correlation potential is
current dependent.~\cite{vignale,tokatly}

An appealing feature of the present approach is that electron-electron
and electron-ion correlations and dissipation would in principle be described
correctly in two-component TDDFT.~\cite{kreibig}

At present we are implementing our propagation scheme
for real 3D systems~\cite{wirtz}  within the real-space real-time TDDFT
code, {\em octopus}.~\cite{octopus} We are also
exploring the possibility of a semiclassical description of the
electron-ion coupling in order to avoid the complexity of multicomponent
DFT and the problems related to mixed quantum classical approaches
({\em i.e.}, Ehrenfest dynamics) as they fail to describe the long-term
inelastic electron-phonon scattering correctly.~\cite{alberto,todorov,car}

\section*{Acknowledgments}
This work was supported by the European Community 6th framework Network of
Excellence NANOQUANTA (NMP4-CT-2004-500198) and by the Research and Training
Network EXCITING. AR acknowledges support from
the EC project M-DNA (IST-2001-38051),
Spanish MCyT and the University of the Basque Country. We have benefited
from enlightening discussions with  L. Wirtz, A. Castro, H. Appel, M.
A. L. Marques, C. Verdozzi and U. von Barth.

redpajama_set_name: RedPajamaArXiv


SEOUL: Soldiers from North and South Korea crossed into each other's territory peacefully for the first time on Wednesday, Seoul said, as they began checking the dismantlement of guard posts in the Demilitarized Zone.
The North invaded the South in 1950, triggering the Korean War, and Seoul went on to change hands four times as Pyongyang's Chinese-backed forces and the US-led UN troops supporting the South fought their way up and down the peninsula and back again.
The conflict ended in an armistice in 1953, leaving the two technically still in a state of war, but a rapid reconciliation has taken place this year. The South's President Moon Jae-in and Kim Jong Un, leader of the nuclear-armed North, agreed to remove a limited number of guard posts along their heavily-fortified frontier at a September summit in Pyongyang, among other measures.
North Korea blew up 10 of its facilities in November, while the South tore 10 down on its own side using excavators. South Korean inspectors visited each of the guard posts on the North's side on Wednesday to make sure they had been dismantled and all firearms and troops removed, Seoul's defence ministry said.


text:
180 gram audiophile vinyl. Gatefold jacket. 4-page booklet. Limited Edition of 1500 individually numbered copies on translucent blue coloured vinyl. Recorded in support of 2009's Grammy-winning album Already Free, the double live album Roadsongs contains over 100 minutes of music from across their recording career. The recording is full of dynamics and a warm, present sound while capturing the raw, spontaneous energy of a live show.
Label: MUSIC ON VINYL
Roadsongs (Blue) [Colored Vinyl] (Gate) [Limited Edition] [180 Gram] (Hol)
Artist: Derek Trucks Band
Vinyl - Colored Vinyl] (Gate) [Limited Edition] [180 Gram
$14.99 Buy MP3 Album
1. I'll Find My Way (Live)
2. Down In The Flood (Live)
3. Sailing On (Live)
4. Get What You Deserve (Live)
5. Days Is Almost Gone (Live)
6. Already Free (Live)
7. Afro Blue (Live)
1. I Know (Live)
2. Down Don't Bother Me (Live)
3. Don't Miss Me (Live)
4. Rastaman Chant (Live)
5. Key To The Highway (Live)
6. Get Out My Life Woman/Who Knows (Live)
7. Anyday (Live)

redpajama_set_name: RedPajamaCommonCrawl


Christ was born to fulfill God's plan of salvation for mankind. He fulfilled all that was foretold concerning Him by the prophets and suffered by the hands of wicked men. There is one more prophecy concerning Christ's return and His coming Kingdom that is yet to be fulfilled. The same Christ who died and rose again will return in His glory and those who are His will be united with Him. Will you be part of His family who will meet the Lord in the air? Believe in the Lord Jesus Christ and trust Him for your redemption and He will save you. Praise God!


We are extremely proud to have trio of Specialist Oral Surgeons who bring a wealth of experience to our dental team.
As well as providing care for our patients , we accept referrals in excess of 100 per month from surrounding practices.
Our team can provide oral surgery with or without sedation depending on the complexity of the case or if a patient is anxious. Rest assured that our patients are in very experienced caring hands.
Diana graduated from the University of Birmingham in 1992 and completed her vocational training year in Bromsgrove, before embarking on a hospital career in 1994, subsequently working as a Senior House Officer in hospitals across the North West and the West Midlands before attaining her Fellowship in Dental Surgery in 1997.
Bayan qualified in 1977 with a BDS degree from the College of Dentistry- University of Baghdad. After basic dental qualification, he spent one year in general dental practice in Baghdad followed by two years full time postgraduate course at the University of Baghdad and received the High Diploma in Oral Surgery.
Tom was born in Nantwich. He was educated at Sandbach School and subsequently at The London Hospital gaining a B.D.S. in 1986. After a year in Southend Hospital he returned to Cheshire. Tom has worked in departments of Oral and Maxillofacial Surgery in both Chester and Leighton Hospitals for 25 years.


I'm knitting something that isn't a toy or a CounterPain Square.
Ellie needs some mittens. I'm doubling Rowan Wool Cotton on US 5s (discontinued color – I had to find more on Ebay to match her existing hat) and working from the Handy Book of Patterns by Ann Budd. Great book, by the way. I've mentioned it before.
For mittens you need a couple stitch markers. Well, you don't NEED them need them, but if you want to not have to pay a lot of attention, stitch markers are helpful.
Between the mitten leftovers and unraveling her hat, I hope to have enough yarn for a larger hat. I should have just started over with new yarn on the set – but the hat fit when I committed to the mitten yarn.
This entry was posted in Yarn Miracle and tagged ellieivey, knitting, mittens on December 30, 2010 by emily.
This entry was posted in Yarn Miracle and tagged christmas ornament, felted arorn on December 25, 2010 by emily.
This cough is killing me.
This entry was posted in Yarn Miracle and tagged etsy, treasury on December 19, 2010 by emily.
Candy Cane Ice Cream in my coffee.
I'm a coffee genius. That must be why I am beating the pants off of Christmas this year.
I'm going to drink this while I work on my last Custom Order for the holidays.
Then I'm going to knit my kid a big bear and a pair of mittens.
This entry was posted in Yarn Miracle and tagged coffee, recipe on December 16, 2010 by emily.
This entry was posted in Yarn Miracle and tagged etsy, treasury on December 12, 2010 by emily.
One hundred and ten CounterPain squares made and assembled. Only thirty more to go.
There is a teensy glimmer of light at the end of the tunnel.
This entry was posted in Yarn Miracle and tagged epic counterpain, knitting, my patterns on December 8, 2010 by emily.
Three out of three Iveys agree: Faux Thin Mints are a mighty fine cookie.
The recipe (which I saw on Craft) was posted at A Billion Tastes and Tunes. Our only change was to add a little powdered sugar to the glaze and make sandwich cookies instead. In theory, sandwich cookies will store better (they are stackable) and not as messy to eat. In theory. The cookies have a great tender crispy texture that is hard to beat.
Mix the cornstarch and water and set aside.
Mix margarine and sugar together until well-mixed and fluffy (I used a fork).
Add cornstarch mess and extracts and mix until smooth and consistent.
Add flour, cocoa, baking powder, and salt. Mix until just combined (try not to over-mix, it breaks the baking powder). It's stiff, but it's worth it. So persevere.
Divide the dough in half. Put one half on a sheet of wax paper or plastic wrap. Roll it into about a 10 inch log, wrap it up in the wax paper and put it in the 'fridge for about two hours until firm. If you are in a hurry, the freezer for an hour will work too. Do the same for the second half of the dough.
Preheat your oven to 350. Pull out out the dough, unwrap it and slice it into rounds a little thicker than 1/4″. Spread them out on a cookie sheet – they don't need a lot of room since they don't spread much.
Bake about 12 minutes. The cookies should look a little dry on top but still be soft enough to get a dent when you poke them accidentally. Let them cool a few minutes on the cookie sheet before moving them to a wire rack to cool all the way.
While they are cooling, make your filling. Put the chocolate chips and oil in a microwave safe bowl and heat in 30 second increments until the chocolate chips are melty and the mess can be stirred until smooth. Mix in the peppermint extract (if your bowl is still warm, it shouldn't seize). Stir enough powdered sugar to make the mixture spreadable and a little fluffy.
Using a dinner knife, spread the frosting over the bottom of one cookie. Stick another cookie (bottoms in for prettiness) on top. Repeat until all the cookies are sandwiches or eaten.
Speaking of Cookies. My sister is hosting a Virtual Cookie Swap! With PRIZES! Take a look and share a recipe. I need more cookie recipes.
This entry was posted in Yarn Miracle and tagged chocolate, cookies, recipe on December 6, 2010 by emily.
This entry was posted in Yarn Miracle and tagged etsy, treasury on December 5, 2010 by emily.


TOKYOPOP has licensed the 15-volume Digimon series drawn by a Hong Kong-based artist. They expect to begin releasing it in late 2002. They've also confirmed that Nelvana does indeed have the license to the Groove Adventure Rave Anime, while TOKYOPOP holds the manga license. They are also looking to expand their number of Korean, Chinese and Taiwanese comic titles.
TOKYOPOP will be release two hardcover Art of Rayearth books by October. The books will sell for either $14,99 or $19.99.
Regarding the acquisition of older CLAMP tiltes, such as RG Veda and Tokyo Babylon , TOKYOPOP 'S Joel Burrell stated, "Anything we can get, we will get."
Given the fees involved, the chance of for Kodomo no Omocha TV series being licensed for US release is "unrealistic," although Tokyopop is hoping some arrangement can be made.


What is the scope of the IPACS project ?
IPACS wants to define a new basis for benchmarks measuring system performance of distributed systems. The goal is to develop methods of measurement based on new applications, low level, and grid benchmarks, which also allow performance prediction. A further significant element of the new benchmark environment is represented in the usability of the benchmarks. Mechanisms for the automatic installation and execution of benchmarks are being strived for and grid benchmarks, which also allow performance prediction.
More information about the IPACS BOF Session could be found here.
Autopar is an extension of the GNU autotools with respect to different mpi versions, parallel operating environments and batch systems.
More Information are available at http://sourceforge.net/projects/autopar/.
In the IPACS BOF Session at the ISC 2006 in Dresden (www.supercomp.de) the IPACS experts will present the IPACS benchmark execution environment, the comprehensive IPACS benchmark suite and new methods for performance modeling and performance prediction. More Information are available here.
An article about the IPACS Benchmark Suite, Performance Modeling and Prediction Methods and the Benchmarking Environment has been published in the PIK journal (Praxis der Informationsverarbeitung und Kommunikation), Vol. 03/05, pp. 150-160. For more information please visit K.G. Saur Verlag GmbH.
To fascilate the execution of IPACS benchmarks please use our new IPACS client, which is implemented as a Java GUI application. The client software can be downloaded here.
Visit the project partners of IPACS at ISC2005 in Heidelberg, Germany. Ask our experts about the project and learn about the new benchmark client.
The IPACS project has been presented at SC2004 in Pittsburgh, PN during the TOP500 BoF session. Additionally flyers about the project have been distributed at this event.
Presentation at the 2nd Workshop on Distributed Objects Research, Experiences & Applications (DOREA 2), June 16-18, 2004, Las Vegas, USA in 3rd International Symposium on Information and Communication Technologies.
The IPACS project partner has been invited to present at the ISC2004 tutorial 'Performance in HPC: Evaluation, Modeling, Benchmarking and Prediction' in Heidelberg.


text:
An AssuredPartners Company
About Walsh
Walsh Charities
HR Workplace Services
Quoting Options
Quoting Overview
Small Group Census Forms
Large Group Quoting Templates
Carrier Info & Quoting
Compliance / Education
About Us > Walsh Charities
Walsh Family Charities
Each year, Walsh Family Charities hosts an annual St. Patrick's Day fundraiser. Together with your support, Walsh Family Charities has been able to support the extraordinary work of Dr. Patrick Borgen, MD in his fight against Breast Cancer…
Mammography, the cornerstone of breast cancer diagnosis, is a 50 year old test that has saved countless lives by diagnosing breast cancer at a stage when the disease is highly curable. But this venerable screening test is far from perfect. The false negative rate associated with mammography ranges from 10% to 40%. This means that thousands of women in the United States each year are falsely reassured by a mammography report that is 'normal'. Moreover, mammography involves compression of the breast and exposes the patient to low levels of radiation. Why has more progress NOT been made?
The Breast Specific Gamma Imaging (BSGI) technology represents the first true advance in breast cancer screening in half a century. Rather than relying on anatomic changes in the breast, the BSGI identifies FUNCTIONAL changes that are associated with cancer. The technology captures vital tumor information by viewing the metabolism of cancerous growths by measuring the uptake of a radiopharmaceutical tracer. The mitochondria of cancer cells preferentially uptake this unique dye and the Dilon 6800 gamma camera is able to identify this process in the breast. Cancerous growths just millimeters in size can be visualized with this novel device. The test does not involve compression of the breast and does not expose the breast to radiation. The test is FDA approved and is covered by virtually all insurance carriers.
The Brooklyn Breast Cancer Project at the Maimonides Cancer Center, through a generous grant from the Walsh Family Charities and their supporters has brought this exciting technology to our patients in Brooklyn and the Greater New York area. For more information about scheduling a Breast Specific Gamma Imaging test please call (718) 765-2570.
Dr. Patrick Borgen, MD
Director, Brooklyn Breast Cancer Project
Click here for details on this year's event on June 7, 2018.
Module 'WI Wow Gallery' is not activated.
(c) 2018 AP Benefits Advisors LLC dba AssuredPartners WB
Fair Haven , NJ 07704
© 2019 AssuredPartners Inc. All rights reserved.


There is never enough of it. But once harnessed, money can be a great power.
Hi I am Hari..and if you are at a stage where there are more questions than answers, you have clicked the right link...let us see where you are..well approximately...you have moved out of your parents home.....You may have found a steady partner....and thinking of a new life together. And the first step in most cases is buying a new home.
But along with that, many new concepts come into play - relationships to be managed, time to be managed and yes, money to be managed. The famous question - what flows more easily - money or water? - has a slightly discouraging answer at this stage. Because money tends to flow out just as soon as it comes in. But there are ways to master it and make it work.
Money tends to flow out just as soon as it comes in. But there are ways to master it and make it work.
Many new concepts come into play - relationships to be managed, time to be managed and yes, money to be managed.


at the heart of your home.
J&B are proud to supply Neptune kitchens. We offer a full design and fitting service, all you have to do is provide us the blank canvus for us to create you your dream kitchen.
The Chichester range has been designed to produce a kitchen of the very highest standard at a fantastic price that is truly simple to fit.
The versatile style can be used to create either a contemporary or traditional look, matching perfectly with the Chichester free standing furniture.
Our Suffolk range is a traditional shaker style, hand painted in dove grey with honed slate interior, chrome handles and catches. Neptune's Suffolk kitchen is quite simply stunning.
Purchasing a kitchen can often be a daunting process, not with us! We have fabulous Neptune kitchens with matching free-standing furniture to complete the look. We offer a full design service from the smallest to the biggest of kitchens, we can help you plan, design, accessorise and fit your dream kitchen, pop to our showroom for a browse and a chat and see what we can do for you!


text:
Innovation Ohio
A New Kind of Policy Think Tank
Fair Elections
Ohio Lawmakers Working to Ban Child Marriages
Wednesday, March 21st, 2018 by Terra Goodnight
Today, an Ohio House committee will hear from supporters of new legislation (House Bill 511) that would change Ohio's antiquated laws on child marriage.
The law on the books in Ohio currently requires only that boys must be 18 to marry, while girls can marry as early as age 16. That's sexist and must change.
What's even more disturbing, however, is that Ohio law permits girls under 16, with no established minimum, to marry if they are pregnant and have parental approval. This leads to situations in which many girls, already facing an unwanted pregnancy, are pressured to marry the adult men who impregnate them out of a desire to make the birth 'legitimate.' In fact, "parents' rights" groups were the main opposition to similar legislation that was later amended and passed in the Kentucky legislature.
Girls who marry young are more likely to drop out of high school, and four times less likely to finish college, resulting in a lifetime of diminished earning potential as a result.
House Bill 511 is a bipartisan proposal that would raise the marriage age to 18, with an exception for 17 year olds who have been emancipated and prohibits marriage between a minor and a much-older adult, setting a 4-year age difference for a minor allowed to legally marry.
Today, IO's Erin Ryan, who also serves a Managing Director for the Ohio Women's Public Policy Network, will testify in support of HB511. Included in her testimony are these devastating facts:
Evidence suggests that the age differences in most of these marriages between teen girls and adult men can lead to unequal power and control dynamics, creating an environment rife with domestic violence and abuse, which often continues throughout the marriage. Women who marry as minors are significantly more likely to have mental and physical health issues as a result of the abuse they suffer, and girls marrying under 18 are three times more likely to experience domestic violence. Due to the age of the child at the time of the marriage or the age difference between the child and their spouse, some of these situations constituted statutory rape under Ohio's state law. However, because of the outdated laws dictating marriage in Ohio, these cases lead to marriage licenses, not statutory rape charges.
You can read Erin's full testimony here. To contact your State Representative and ask them to support House Bill 511, use our letter-writing tool.
360 S. 3rd Street, 3rd Floor
info@innovationohio.org
© 2016 Innovation Ohio
For more non-partisan policy
reports, visit the 501(c)(3):


OnlineTips.org respects the privacy of every individual who visits our website. This Privacy Statement describes the information collected and how it will be used.
OnlineTips.org will not collect or retain any personal information such as your name, address, telephone number, or e-mail address, unless you choose to provide it. If you provide this information, OnlineTips.org will only use it for our own information. We will not sell or otherwise divulge the information to 3rd parties.
OnlineTips.org collects certain technical information when you visit http://www.OnlineTips.org. For example, the domain name from which you accessed the Site, the pages visited and the date and time you accessed the Site are analyzed to make the Site more useful to visitors. But, we only use IP addresses and these are not linked to personally identifiable information.
OnlineTips.org contains links to merchant sites. Be aware that OnlineTips.org is not responsible for the privacy practices of these sites. We encourage you to read the privacy statements of each and every web site that collects personally identifiable information. This privacy statement applies solely to information collected by this website.


David Zibell emigrated from Canada to Israel with a vision to establish a distillery near the Golan Heights Winery, aptly naming it the Golan Heights Distillery. I met with David on a warm and muggy night in Tel Aviv, and had a chance to both hear his plans and taste the single malt new make, and a single malt and a malted rye spirit six months into the aging process in new oak barrels.
David has established a good working relationship with the Winery, and hopes to start using wine casks from the winery beginning in the next few still runs. In the mean time, I got a very preliminary look at what he's preparing.
At this point, David is doing everything on his own, including the cooperage. He shaves about a half an inch off the staves, then chars the barrel. His method is less than traditional, and his mash is grain heavy and is more a wort (fermented and distilled with the grain) than a wash (hot water mixed with grist, with the wash drained away for fermentation without the actual grain). It is then distilled in a 500 liter still, averaging four distillations a week.
I started with the new make, cut down to 50%. The nose is very earthy and oily with a lot of grain, wood fire smoke and nuts. This new make is far less sweetness than the Scotch new makes I've had. The palate is overly toasted, almost burnt, with a lot of grain and peppery spice. The linger delivers spice on the throat. This is the basic character carried through the whisky. The new make still needs some development, as I found it a bit short the fruity/floral dimension, and I'm sure David will continue perfecting it as he goes along.
The whisky then goes into a cask, in the case of cask 01 this is a new oak that David charred, and after six months it's presented at 63% ABV, is light bronze with a nose of wet earth, cooked buckwheat (kasha), pepper, vanilla and lavender. It's much softer than the new make. On the palate it has faint honey, harsh green chili pepper, as spicy as any rye whiskey you know, with notes of cotton candy on the back of the dram. The linger is long and drying, with heat high up in the back of the throat and pepper on the tongue with notes of cereal.
I commented to David about how strong the cereal is throughout and he said that he's aimed to get an earthy, cereal-y product, and while I'd probably aim for more fruity esters, he has definitely hit the notes he aimed for, which attests to his skill and his learning from Brian Davis, his mentor.
Next up was the Hermon Single Malt Rye, presented at 66% ABV and, once again, aged six months in a new oak cask. The rye was my favorite of the tasting, and is pale gold with thin legs. The nose is fascinating with a peanut butter and jelly sandwich (seriously!), citrus, light honey, lemon, nutty (hence the peanut butter note). Water teases out some sourness. The palate sees honey before a hit of spice, pepper, lemon and the black herbal licorice candy you get in Germany. The linger leaves honey on the tongue, lots of spice in the mouth but almost none in the back of the throat.
While David had none of the corn mash distillate with him, my fellow blogger Yoav (of the excellent Whisky Gospel) had some of it left, so I tasted that too. It's grain in the extreme, but not sweet. It obviously came out of the same still as the single malt and rye, and is sweet, but not in a corn way.
For interim income, he's also making absinthe, which I tasted at the end of the evening. It was tasty, but I have very little knowledge of the drink to actually assess in comparative terms.
This is a space to watch, especially the rye whiskey, and I look forward to further updating you on the distillery's progress.


text:
Q: concatenating string with tuples in python def print_distance(p1, p2):
    dist = (((p2[0]-p1[0])**2)+((p2[1]-p1[1])**2))**0.5
    d = print("{0:.2f}".format(dist))
    result = ("Distance between ", ','.join(p1) ," and ", ','.join(p2) ," is "+ d)
    return result
print_distance((0,1),(-1,-2))

TypeError: sequence item 0: expected str instance, int found
I couldn't figure it out the output "Distance between (0, 1) and (-1, -2) is 3.16" I am beginner to this thing. Help me.

A: your tuple contains ints so you cant join them with str.join. if you wanted to do that you would need to convert them first to str's. however to get the output you mentioned you can just use an fstring
def print_distance(p1, p2):
    dist = (((p2[0]-p1[0])**2)+((p2[1]-p1[1])**2))**0.5
    print(f"Distance between {p1} and {p2} is {dist:.2f}")

print_distance((0,1),(-1,-2))

OUTPUT
Distance between (0, 1) and (-1, -2) is 3.16

A:     result = ("Distance between ", ','.join(p1) ," and ", ','.join(p2) ," is "+ d)

this is line is just a bit of a mess maybe try using format the way you did elsewhere
    result = "Distance between {p1} and {p2} is {d}".format(p1, p2, d)

I'd like to give you more context but the gist of it is you are using join incorrectly.

redpajama_set_name: RedPajamaStackExchange


The GAMMAFLEX-RS is specially designed for the printing of papers from 16 grammes/m2. Its system of exit into reels and formats makes it ideal for the printing of stamped paper in reams and counter reels.
Robust, compact, closed and protected printing units in order to avoid ink splashing even at high speeds.
Equipped with slow rotation and automatic hydraulic lifting in the stops of the machine. It includes quick advancement by gears system and racks. Dosage of the ink by inking roller and anilox mechanical plot or by means of close chamber doctor blade with pneumatic pressure and ceramic roller.
Drying system composed of high performance screens mounted between the printing units and the upper bridge. Impulsion and expulsion fans for the gas extraction.
Air heating by heating resistors/electric, gas burners or thermal oil. It incorporates an automatic control of temperature by PID regulator.
Simple station unwinder for reels up to 800 mm. of diameter, equipped with electromagnetic break operated by a potentiometer. Regulation of the tension by acquisition of non-contact diameter. Fixation of the reels by pneumatic expansible shaft.
Simple station rewinder for reels up to 800 mm. of diameter, equipped with a direct current motor and operated by a digital variator. Automatic regulation of the tension by pre-selection in a potentiometer. Fixation of the reels by pneumatic expansible shaft.
Serrated sheeter at fix format, specially suitable for the cutting of thin papers from 16 grammes/m2 with accumulation and pneumatic expulsion. It includes a table for collection.
IZARFLEX is the new Giave recent generation gearless stack-type flexographic printing press for high quality of printing and fast-changing jobs.
Preparation of the available printing units for the following job in production, allowing the setup job changes at production speed.
In the MIZARFLEX series, positioning and pressure adjustment of the counter-pressure roller and plate cylinders operated by independent servomotors. Printing unit positioning with automatic adjustment according to printing repeat. Giave uses servo motors to facilitate all required operations fast and efficiently.
The integrated control panel allows the centralized control by an industrial PC and a touch screen, including all the necessary control for the use of the machine.
Drying including two completely independent circuits. Module for the drying between printing units and high performance horizontal tunnel assuring the material drying before rewinding. Giave offers up to 60% recirculation of processed air, which lowers the cost of drying operations for our customers.
MIZARFLEX series offer unwinders and rewinders of totally automatic change or double-station for semi-automatic change. Diameters to 1.200 mm.
In-line with MIZARFLEX, it exists the possibility of incorporating equipments for the solvent lacquered or coat, water-based or UV drying register with flexographic printing.
Double turret slitter-rewinder suitable for plastic films, laminates and papers, for working at a speed up to 600 m/min. Its double turret rewinding system and the automatic unloading reduce drastically the reel change cycle, achieving an important production increase.
Available widths: 1.000 and 1.300 mm.
Maximum production speed: 600 m/min.
Shaftless unwinder for diameters up to 1.200 mm.
Diameter of finished reels: 600 mm.
Precise and uniform tension control during the whole rewinding process.
Minimum time of cycle in reel charge.
Industrial PC with an easy handling touch screen and job memorization.
Cutting unit equipped with circular and razor blades.
Double turret rewinder equipped with differential friction shafts.
Leaders in the manufacture of machinery for the rotogravure and flexographic printing, and slitting and rewinding, at Giave we offer a personalized service with which we adapt our products to the specific needs of your company by optimizing performance and productivity.
At the vanguard of the latest advances in the field of the printing and converting of flexible materials, at Giave we adopt and develop the latest technological innovations in order to maximize the performance of our equipments and the production capacity of your business.
Our commitment to excellence is embodied throughout the production process of our printing and converting presses, and even earlier in the selection on brands and partners whom we collaborate. Proof of this is the maintenance services we continue making in machines installed over 20 years.
The development of our equipment is based on innovation to maximize the profitability of the production process and the competitiveness of our customers. Design innovation, services innovation and above all, technical innovation that integrate with ongoing improvements in our machines.


Movie subtitles related to beyond the walls (2012) at subtitlecube.com! Download this subtitles from our download location, see more details bellow. It's free and easy to download subs for movie beyond the walls (2012).
Welcome To Death Row. The Doctor Will See You Now.