text:
\section{Introduction}

Graphene, a two-dimensional hexagonal lattice of carbon atoms isolated
in 2004 by Geim and Novoselov~\cite{geim2004}, represents a very interesting material
in several fields of science and technology. Since its low-energy
properties can be described by a Dirac equation, it is considered to be an
ideal test bed for the investigation of relativistic effects at non-relativistic
velocities~\cite{katsnelson2007}.
Moreover, its unique properties make it suitable for applications in many different
fields of technology~\cite{geim2009}. Its atomic thickness and high room-temperature mobility,
for example, make it a candidate material for postsilicon electronics.

For applications in digital electronics the presence of a sufficiently large energy gap is
fundamental~\cite{iannaccone2009}. A gap can be induced
by the lateral confinement
in narrow ribbons with a transverse
size of a few nanometers or of tens of nanometers, which can be efficiently
modeled with atomistic techniques, such as tight-binding approaches.
Fabrication of such nanowires is, however, very challenging, and therefore also
alternative and/or complementary approaches to open up a gap are being
pursued, such as the usage of bilayer graphene~\cite{mccann2006,ohta2006,
castro2007,zhang2009}, of chemical
functionalization~\cite{ouyang2008,echtermeyer2009,elias2009,boukhvalov2009},
and of doping~\cite{acs2012,iwceroche,tang2012,wei2009}.
Large graphene devices (with a size of several hundreds
of nanometers or of microns) can, however, be convenient
(or mandatory) in radio-frequency or sensor applications,
which do not necessarily require
an energy gap~\cite{novoselov2012,schwierz2010,mina2013,hill2011,bonaccorso2010,xia2009}.
In these cases an atomistic analysis is not feasible, since it is numerically
too expensive, and more approximate techniques (such as envelope-function descriptions)
have to be adopted. The availability of techniques suitable for the simulation of
large-area graphene structures is also essential for a direct comparison with
transport and noise measurements
performed on micron-sized flakes (see for example Refs.~\cite{cooper2012,peres2010,danneau2008,balandin2013}).

Numerical efficiency is particularly important for transport calculations,
which are usually self-consistently coupled with the solution of the Poisson
equation, and  therefore
typically need hundreds or thousands of iterations to reach
global convergence.
As a consequence, numerical performance can often make
the difference between a feasible
calculation and a computationally  impossible task. For this reason,
it is important to develop reliable algorithms that are more
approximate but more efficient than the atomistic ones.

A continuum approach has been often used for the numerical simulation
of transport in ribbons made up of transverse regions
with constant potential~\cite{wurm2009,li2009,zhou2010}, for which
analytical expressions for the wave function in each region are available
\cite{castroneto2009,brey2006,kp2011}.
However, for a generic potential, the envelope function equation (which
is a Dirac equation) has to be solved numerically, and
the relativistic dispersion relation introduces some complications
in the standard discretization schemes.

In order to avoid these difficulties, the first numerical studies~\cite{Bardarson,Nomura}
adopted a momentum space regularization of the Dirac equation.
For ribbons with a large aspect ratio, such as
those studied in Refs.~\cite{Bardarson,Nomura},
the boundary conditions are expected to be largely uninfluential and thus periodic
boundary conditions were used. Tworzid{\l}o {\it et al.}~\cite{tworzydlo2008}
later performed a real-space transport analysis of large aspect ratio graphene
ribbons by adopting the Stacey discretization scheme~\cite{stacey1982,
bender1983}.
More recently, Hern\'andez {\it et al.}~\cite{hernandez2012} studied
the transport properties of zigzag and armchair ribbons in the direct space,
with physical boundary conditions.  For the zigzag ribbon, in the transverse direction
they adopted the discretization by Susskind~\cite{susskind1977}
and in the longitudinal direction the discretization by Stacey;
for the armchair ribbon they used the Stacey discretization in both directions.
Snyman {\it et al.}~\cite{snyman2008} previously adopted the alternative
method of mapping the Dirac equation onto a Chalker-Coddington network
model~\cite{chalker1988},
which in the past was used to describe percolation in disordered samples
in the quantum Hall effect~\cite{ho1996,connolly2012,logoteta2012}.

After the seminal contributions by Bardarson \textit{et al.}~\cite{Bardarson} and by
Nomura \textit{et al.}~\cite{Nomura} the attention of the graphene community
mainly focused on the finite difference methods, for which
the physical boundary conditions can be easily
implemented~\cite{CCnote}. In this paper we will show how to extend
the reciprocal space algorithm to the case of armchair boundaries, and we will
present compelling evidence that this method is numerically much more efficient than the finite difference schemes.

For the solution of the transport problem we will adopt a
scattering-matrix approach. The ribbon is partitioned into
a series of thin slices in the direction along the current flow
(which we will refer to as the longitudinal direction).
In each slice the potential is approximated with a longitudinally constant function,
and the Dirac eigenvalue problem is solved
in the Fourier transformed space. The conductance of the whole structure is then
evaluated applying a mode-matching procedure at the interfaces between adjacent slices.

The paper is organized as follows.
In Sec.~\ref{s:problem} we introduce the $\vec k\cdot \vec p$ approximation
for graphene ribbons and present the equations that describe the
transverse slices of the device. In order to perform a comparison with
the approach we have decided to adopt, in
Sec.~\ref{s:direct} we outline a few finite-difference techniques
that could be employed to solve the envelope function equation
in the direct space. In Sec.~\ref{s:reformulation} we
introduce a mapping of the armchair problem into one with periodic boundary conditions
and in Sec.~\ref{s:fourier} we provide a detailed discussion of the reciprocal space technique
used in Ref.~\cite{PT}, where numerical precision was crucial.
In Sec.~\ref{s:comparison} we compare the numerical efficiency of the
methods in the direct and reciprocal space.
Section \ref{s:transport} is devoted to the solution of the two-dimensional
transport problem for an armchair graphene ribbon in the realistic situation
of a potential that varies also in the longitudinal direction.

\section{The numerical problem}\label{s:problem}

The wave function of graphene can be approximated by means of a linear
combination of the $2 p^z$ orbitals of the carbon atoms of its two inequivalent
sublattices (see, e.g., Refs.~\cite{castroneto2009,ando2005,kp2011}).
Using the subscripts $A$ and $B$
to distinguish the quantities associated with different sublattices, the wave
function is written as
\begin{equation} \label{total}
\psi (\vec{r}\,)=\sum_{\vec{R}_A}\psi_A (\vec{R}_A)\varphi(\vec{r} -\vec{R}_A)+
\sum_{\vec{R}_B} i\,\psi_B (\vec{R}_B)\varphi(\vec{r} -\vec{R}_B) \, ,
\end{equation}
where the sums are over the atom positions, and $\varphi(\vec r\,)$ is the $2 p^z$ orbital.

We consider an armchair ribbon with $N_D$ dimer lines
of carbon atoms across its width. The distance between the opposite
edges  is $W=(N_D-1)\,a/2$, where
$a=\sqrt{3}a_{C-C}$ is the graphene lattice constant and $a_{C-C}\simeq 0.142$~nm is the distance
between nearest-neighbor atoms. We denote by $x$ and $y$ the longitudinal
and transverse directions, respectively; the Dirac points $\vec K$ and
$\vec K^\prime$ can be expressed as $\vec{K}=-K \hat{y}$ and $\vec{K'}=K
\hat{y}$, with $K=4 \pi/(3 a)$.

For large enough $N_D$ and a potential that is slowly varying on the lattice
scale, atomistic details can be disregarded and the low-energy physics can be
described by means of a $\vec k\cdot\vec p$ approximation. This
is achieved in practice by promoting the coefficients
$\psi_A(\vec R_A)$ and $\psi_B(\vec R_B)$ in \eqref{total}
to continuous functions of the position ${\vec r}$
and by writing the functions $\psi_A (\vec{r})$
and $\psi_B(\vec{r})$ in terms of four envelope functions $F$
as (see, e.g., Refs.~\cite{castroneto2009,ando2005,kp2011})
\begin{equation}\label{psia}
\begin{aligned}
\psi_{\beta} (\vec r) =
e^{i \vec K \cdot \vec r} F^{\vec K}_{\beta} (\vec r)-i\, e^{i \vec K' \cdot \vec r} F^{\vec K'}_{\beta} (\vec r)\, ,
\end{aligned}
\end{equation}
where $\beta=A,B$.

It can be shown (see, e.g., Refs.~\cite{castroneto2009,ando2005,kp2011})
that the envelope functions satisfy the massless Dirac equation
\begin{equation}\label{dirac0}
\begin{aligned}
\displaystyle
-i\gamma(\partial_x\sigma_x+\partial_y\sigma_y)\vec{F}^{\vec{K}} &= E\vec{F}^{\vec{K}} \\
\displaystyle
-i\gamma(\partial_x\sigma_x-\partial_y\sigma_y)\vec{F}^{\vec{K'}} &= E\vec{F}^{\vec{K'}}\,,
\end{aligned}
\end{equation}
where
\begin{equation}
\vec{F}^{\vec{\alpha}}=\left[\begin{array}{c} F_A^{\vec{\alpha}} (\vec{r})\\ F_B^{\vec{\alpha}} (\vec{r})\end{array}\right]\,
\end{equation}
(with $\vec{\alpha}=\vec{K},\vec{K'}$), $\partial_x=\partial/\partial x$, $\partial_y=\partial/\partial y$,
$\sigma_x$, $\sigma_y$ are Pauli matrices, $E$ is the total
energy of a particle in the ribbon, and $\gamma=(\sqrt{3}/2)\gamma_0 a\equiv v_F \hbar$.
The constant $\gamma_0\simeq 2.7$~eV is the modulus of the transfer integral between nearest-neighbor carbon
atoms, $v_F$ is the Fermi velocity of graphene, and $\hbar$ is the reduced Planck constant.
Within the $\vec k\cdot\vec p$ approximation, the presence of an
external electric field is handled by introducing the potential energy $U(\vec r)$ into Eq.~\eqref{dirac0}:
\begin{equation}\label{dirac}
\begin{aligned}
\displaystyle
[ -i\gamma(\partial_x\sigma_x+\partial_y\sigma_y)+U(\vec{r})I ]\vec{F}^{\vec{K}} &= E\vec{F}^{\vec{K}} \\
\displaystyle
[ -i\gamma(\partial_x\sigma_x-\partial_y\sigma_y)+U(\vec{r})I ]\vec{F}^{\vec{K'}} &= E\vec{F}^{\vec{K'}}\,,
\end{aligned}
\end{equation}
where $I$ is the $2\times 2$ identity matrix.

Dirichlet boundary conditions for the wave function $\psi$ have to be imposed on the two dimer lines just
outside the ribbon, at a distance $a/2$ from the ribbon edges, where passivation
approximately takes place. We
choose the origin of the $y$ axis in such a
way that these dimer lines are identified by the conditions $y=0$ and $y=\tilde W\equiv W+a$.
The vanishing of the wave function on the passivation lines leads to the boundary
conditions~\cite{brey2006,castroneto2009,kp2011}
\begin{equation}\label{boundary0}
\begin{aligned}
&\psi_{\beta} (x,y=0)=\psi_{\beta} (x,y=\tilde W)=0
\end{aligned}
\end{equation}
for both sublattices ($\beta=A,B$).

In a waveguide-like configuration in which the potential energy $U$ depends only on the transverse coordinate $y$,
the longitudinal component of the momentum is constant; we will denote the longitudinal wave vector by $\kappa_x$.
The envelope functions that solve the Dirac equation \eqref{dirac} can be decomposed into a propagating wave
along $x$ and a confined component in the transverse direction:
\begin{equation}\label{Phi}
F_{\beta}^{\vec{\alpha}} (\vec{r})=
e^{i \kappa_x x}
\Phi_{\beta}^{\vec{\alpha}} (y)
\end{equation}
($\vec{\alpha}=\vec{K},\vec{K'}$; $\beta=A,B$).
The functions $\Phi_\beta^{\vec\alpha}$ thus satisfy (\emph{cf.} Eq.~\eqref{dirac})
\begin{equation}\label{system}
\begin{aligned}
\displaystyle
\left[\sigma_x f(y)+\sigma_z\frac{d}{d y}\right]\vec{\varphi}^{\vec K}(y) &=-\kappa_x \vec{\varphi}^{\vec K}(y)\\
\displaystyle
\left[\sigma_x f(y)-\sigma_z\frac{d}{d y}\right]\vec{\varphi}^{\vec K'}(y) &=-\kappa_x \vec{\varphi}^{\vec K'}(y)\,,
\end{aligned}
\end{equation}
where we introduced the shorthands
\begin{equation}\label{varphik}
\vec\varphi^{\vec K}(y)=
\left[\begin{array}{c}
\Phi^{\vec K}_A(y)\\
\Phi^{\vec K}_B(y)
\end{array}\right]\, ,\quad
\vec\varphi^{\vec K'}(y)=
\left[\begin{array}{c}
i\,\Phi^{\vec K'}_A(y)\\
i\,\Phi^{\vec K'}_B(y)
\end{array}\right]\, ,
\end{equation}
and $f(y)=[U(y)-E]/\gamma$.

The boundary conditions \eqref{boundary0} become
\begin{equation} \label{boundary1}
\begin{aligned}
\displaystyle
\vec\varphi^{\vec K}(0) &=\vec\varphi^{\vec K'}(0)\\
\displaystyle
\vec\varphi^{\vec K}(\tilde W) &=e^{2 i K \tilde W} \vec\varphi^{\vec K'}(\tilde W)= e^{-i\eta\frac{2\pi}{3}}\vec\varphi^{\vec K'}(\tilde W)
\end{aligned}
\end{equation}
where
\begin{equation}
N_D+1\equiv \eta\!\!\!\mod 3
\end{equation}
and $\eta\in\{-1,0,1\}$. In particular, if $N_D+1=3M+\eta$
(with $M$ an integer), we have
\begin{equation} \label{n0}
\begin{aligned}
2 K \tilde W=2 \pi (2 M+\eta)-\eta\frac{2\pi}{3}\equiv 2 \pi n_0-\eta\frac{2\pi}{3}\, .
\end{aligned}
\end{equation}
The introduction of the discrete variable $\eta$ is convenient
since the product $K \tilde W$ can be very large in the case of wide ribbons,
where the envelope-function approximation is expected to be more reliable.
In fact, $\eta$ is a geometrical property of the lattice structure that goes
beyond the $\vec k\cdot \vec p$ approximation.

It is important to notice that although the envelope functions associated with different Dirac points decouple
in the differential equations \eqref{system}, they are in fact mixed by the boundary conditions \eqref{boundary1}.
This coupling makes it nontrivial to define a symmetric discretization of \eqref{system}.

As previously noted, in the presence of a generic external electric field the differential eigenproblem \eqref{system}-\eqref{boundary1}
cannot be solved analytically: it is necessary to rely on numerical methods in order to obtain approximate
expressions for the transverse components $\Phi$ of the envelope functions and the
corresponding longitudinal wave vectors $\kappa_x$.

\section{Finite difference methods}\label{s:direct}

In this section we describe a few finite difference techniques
that could be adopted to numerically solve Eqs.~\eqref{system}-\eqref{boundary1}
in the direct space.

In a finite difference approach the unknowns are the values of the functions
$\Phi$ on a grid of $N_y$ points along the effective width $\tilde W$ of the
ribbon; here we assume a uniform grid, by setting
$y_i=(i-1) \Delta_y$, with $\Delta_y=\tilde W/(N_y-1)$ and $i=1,\ldots, N_y$.
The derivatives are expressed as linear combinations of the values of the $\Phi$ on a finite number of grid
points, and the boundary conditions
are constraints that reduce the number of unknowns.
As a consequence, the system of equations \eqref{system} is mapped to an algebraic
eigenvalue problem
$A \,\vec v=-\kappa_x \,\vec v$, where the elements of the vector
$\vec v$ are the values of the $\Phi$ at the grid points and the eigenvalues give the longitudinal wave
vectors $\kappa_x$.

In general, as $\Delta_y$ approaches zero, a subset of  eigenvectors of $A$,
together with their respective eigenvalues, converge to the solutions of Eq.~\eqref{system}.
The remaining eigenvectors and eigenvalues are discretization artifacts that have no meaningful continuum limit.
We will refer to them as spurious solutions.

The adopted discretization scheme affects
both the numerical efficiency and the appearance of spurious solutions.
We now discuss the
implications of the simplest discretization schemes:
\begin{enumerate}[(a)]
\item \label{e:asy}Naive asymmetric discretization
\item \label{e:sym}Naive symmetric discretization
\item \label{e:imp}Improved symmetric discretization.
\end{enumerate}

In scheme \eqref{e:asy} the first differential equation of \eqref{system} is evaluated at the
points $y_i$ with $i=1,\ldots,N_y-1$ and the second differential equation of \eqref{system} at $y_i$
with $i=2,\ldots,N_y$. Different representations for the derivatives are used: in the first equation
the two-point forward discretization formula $(d \Phi/d y)|_{y_i}\simeq [\Phi(y_{i+1})-\Phi(y_i)]/\Delta_y$ is used
and in the second one the two-point backward discretization formula $(d \Phi/d y)|_{y_i}\simeq
[\Phi(y_i)-\Phi(y_{i-1})]/\Delta_y$.
The differential equations \eqref{system} are thus mapped to a
$4\,(N_y-1) \times 4\,(N_y-1)$ eigenvalue problem.
We have also considered the alternative scheme of a symmetric discretization formula inside the
ribbon and an asymmetric one at the edges.
In both cases a very slow convergence is observed. Moreover, a large
number of spurious solutions are obtained, which persist also when higher order discretization schemes are used.

As an example, in Fig.~\ref{figureig} we show the eigenvalues $\kappa_x$
obtained for $U(y)=0$ using a three-point
discretization formula, symmetric inside the ribbon and asymmetric at
the edges.
This problem is analytically solvable, and the exact values
of $\kappa_x$ turn out to be either real or purely imaginary
(see, e.g., Refs.~\cite{castroneto2009, kp2011}).
The discretized problem instead also has a large number of complex
solutions with nonzero real and imaginary parts (see Fig.~\ref{figureig}).
In this simple case we can identify them as spurious solutions. However, for a
generic potential energy function $U(y)$,
complex solutions can be physical~\cite{PT}, so they cannot be rejected
\emph{a priori}. If, instead of a three-point formula, we use, for
example, a five-point formula for the spatial discretization, the nonspurious
eigenvalues converge onto the exact eigenvalues more quickly, as the
discretization step is reduced, but the same number of spurious eigenvalues
are present.
\begin{figure}
\begin{center}
\includegraphics[width=8cm]{figureig.eps}
\caption{(Color online) Eigenvalues $\kappa_x$ for a
graphene nanoribbon with 8131 dimer lines ($\approx 1~\mu$m wide) and null
potential energy, for a Fermi energy of 0.1~eV.
The analytical results are compared with those obtained with a
standard three-point discretization scheme (symmetric inside the ribbon
and asymmetric at the edges), using a large number of discretization points
(5000) along the ribbon width.}
\label{figureig}
\end{center}
\end{figure}

In scheme \eqref{e:sym} the grid is modified to have a symmetric discretization in
every point of the grid. Defining $y_0=-\Delta_y$ and $y_{N_y+1}=\tilde W+\Delta_y$, the original
boundary conditions \eqref{boundary1} can be replaced with the relations
\begin{equation}\label{boundary2}
\begin{aligned}
\displaystyle
\vec{\varphi}^{\vec{K}}(y_0) &=\vec{\varphi}^{\vec{K'}}(y_1) \\
\displaystyle
\vec{\varphi}^{\vec{K'}}(y_0) &=\vec{\varphi}^{\vec{K}}(y_1) \\
\displaystyle
\vec{\varphi}^{\vec{K}}(y_{N_y+1}) &=e^{-i\eta\frac{2\pi}{3}}\vec{\varphi}^{\vec{K'}}(y_{N_y}) \\
\displaystyle
\vec{\varphi}^{\vec{K'}}(y_{N_y+1}) &=e^{i\eta\frac{2\pi}{3}}\vec{\varphi}^{\vec{K}}(y_{N_y})\,,
\end{aligned}
\end{equation}
which reduce to \eqref{boundary1} in the continuum limit.
Derivatives in \eqref{system} are evaluated by using the symmetric three-point discretization
$(d\Phi/d y)|_{y_i}\simeq [\Phi(y_{i+1})-\Phi(y_{i-1})]/(2\,\Delta_y)$ in all the points
$y_i$ of the grid ($i=1,\ldots,N_y$) and the differential equations are thus mapped into a
$(4\,N_y)\times (4\,N_y)$ eigenvalue problem.

The eigenvalues of the discretized problem turn out to be always double degenerate.
In detail, each eigenspace is the span of two vectors, let us say  $\vec v_{(c)}$ and $\vec v_{(l)}$,
such that the components of $\vec v_{(l)}$
exhibit even-odd oscillations, while the others have oscillation frequency almost independent of $N_y$.
The eigenvector $\vec v_{(l)}$ cannot have a continuum counterpart, hence the double degeneracy is in
fact a lattice artefact.

This is a clear manifestation of the so-called fermion doubling
problem: a ``naive'' direct space discretization of the Dirac equation results in the appearance of $2^d$
fermions (instead of one) in $d$ space dimensions (in our case $d=1$); this
is an infrared effect, i.e. it does not disappear in the continuum limit $\Delta_y\to 0$. The fermion doubling is a very
well known problem in the field of lattice quantum chromodynamics (see, e.g., Refs.~\cite{rothe, MM}) and is deeply
connected with the chiral anomaly (see, e.g., Ref.~\cite{smilga}), i.e. with the impossibility of regularizing a
theory with massless fermions in a local, chiral symmetric way.
In our simple case, it is caused by the symmetric three-point
discretization formula for the derivative, which involves an incremental step of $2\,\Delta_y$, and hence
decouples odd and even grid points.

Many methods have been developed to overcome the fermion doubling problem;
in scheme \eqref{e:imp} we employ the method proposed in Refs.~\cite{stacey1982, bender1983},
which has a quite simple implementation and was already applied in Ref.~\cite{tworzydlo2008}.
The idea
is to use a symmetric three-point discretization formula for the
derivative, but with an incremental step equal to $\Delta_y$ instead of $2\,\Delta_y$.
This can be done by evaluating the differential equations on
an auxiliary grid,  with nodes at the
center coordinates $y_{i+(1/2)}$ of the cells of the original
grid: $y_{i+(1/2)}=(y_i+y_{i+1})/2$, for $i=1,\ldots,N_y-1$.
The derivative is then approximated by
$(d\Phi/d y)|_{y_{i+(1/2)}}\simeq [\Phi(y_{i+1})-\Phi(y_i)]/\Delta_y$.
The potential energy $U(y)$ is known for every value of $y$ and can be directly evaluated at $y_{i+(1/2)}$, while the
value of the functions $\Phi$ at $y_{i+(1/2)}$ can be estimated by
the average of the values at $y_i$ and $y_{i+1}$:
$\Phi(y_{i+(1/2)}) \simeq [\Phi(y_{i})+\Phi(y_{i+1})]/2$.

The original differential equations \eqref{system} are thus mapped into the
generalized algebraic eigenproblem
$A \,\vec v=-\kappa_x \,B \,\vec v$,
with $A$ and $B$ that are $4\,(N_y-1) \times 4\,(N_y-1)$ matrices.
Since the matrix $B$ is invertible, this problem is in fact equivalent to the
standard eigenproblem
$(B^{-1} \,A) \,\vec v=-\kappa_x \,\vec v$.
We notice, however, that while $A$ and $B$ are sparse
matrices, $B^{-1}A$ is dense. As a consequence, optimized methods
to solve sparse eigenproblems (like the Arnoldi methods) cannot be directly
applied to the standard form.
The sparsity of
$A$ and $B$ can, however, still be exploited in
the multiplication
$\vec{y}=B^{-1}\, A\vec{x}$, which is the fundamental operation
to be performed.
This can be done by carrying out first
the multiplication by a sparse matrix $\vec{z}=A\vec{x}$ and then solving
a sparse linear system
$B \vec{y}=\vec{z}$.
This discretization scheme solves the problems of schemes (\ref{e:sym}) and (\ref{e:imp}): there are neither
spurious eigenvalues nor unphysical double degeneracies.

In the next sections we will show that
this is not the most
efficient way to solve the system of differential equations \eqref{system}.

\section{Reformulation as a problem with periodic boundary conditions}\label{s:reformulation}

The numerical techniques in the real domain that we described in Sec.~\ref{s:direct} do not enable an efficient
numerical analysis of the problem \eqref{system}-\eqref{boundary1}: in order to obtain
high precision results, very large  matrices have to be diagonalized, and the size soon becomes prohibitive.
In this section we reformulate the problem (\ref{system})-(\ref{boundary1}) on a different domain, but
with periodic boundary conditions. In the next section we will show how to solve the resulting numerical problem in the
reciprocal space.

We define the two-component function $\vec\varphi(y)$ by
\begin{equation}\label{varphi}
\vec\varphi(y)=
\left\{
\begin{array}{l @{\quad} l}
\vec\varphi^{\vec K}(y) &
\hbox{ $y \in [0,\tilde W]$}\\
e^{-i\eta\frac{2\pi}{3}}\vec\varphi^{\vec K'}(2\tilde W-y) &
\hbox{ $y \in [\tilde W,2\tilde W]$}\, .
\end{array}
\right.
\end{equation}
From the second of the boundary conditions (\ref{boundary1}) we see that $\vec{\varphi}$ is continuous in its whole
domain, while the first condition gives
\begin{equation}\label{perbc}
e^{-i\eta\frac{2\pi}{3}}\vec\varphi(0)=\vec\varphi(2\tilde{W})  \ .
\end{equation}
Equation~\eqref{perbc} can be interpreted as the requirement
of $2\tilde W$ periodicity
for the function $\exp[i\eta \pi y/(3\tilde W)]\vec\varphi(y)$.

The differential equation satisfied by $\vec \varphi$ can be easily deduced
from Eq.~(\ref{system}) and can
be written in the compact form
\begin{equation}\label{system2}
\begin{aligned}
\displaystyle
\left[\frac{d}{d\,y} \sigma_z+h(y)\sigma_x\right]\vec\varphi(y) &=-\kappa_x \vec\varphi(y)\\
e^{i\eta\frac{2\pi}{3}}\vec\varphi(2\tilde W) &=\vec\varphi(0)\, ,
\end{aligned}
\end{equation}
where
\begin{equation}\label{h}
h(y)=f(\tilde{W}-|\tilde{W}-y|) \qquad  y \in [0,2\,\tilde W] \, .
\end{equation}

In this way we have halved the number of first-order differential equations by doubling the solution domain.
From Eq.~\eqref{system2} we see that
$\kappa_x$ is an eigenvalue of the system [with corresponding
eigenfunction $\vec\varphi_{\kappa_x}(y)$], if and only if $-\kappa_x$, $\kappa_x^*$
and $-\kappa_x^*$ are eigenvalues, as well.
The corresponding eigenfunctions can
be expressed in terms of the eigenfunction of $\kappa_x$:
\begin{equation}\label{eq:Z2}
\begin{aligned}
\displaystyle
\vec\varphi_{-\kappa_x}(y)
&\propto \sigma_y \,\vec\varphi_{\kappa_x}(y)\\
\displaystyle
\vec\varphi_{\kappa_x^*}(y)
&\propto \sigma_x \left[\vec\varphi_{\kappa_x}(2\tilde W-y)\right]^*\\
\displaystyle
\vec\varphi_{-\kappa_x^*}(y)
&\propto \sigma_z \left[\vec\varphi_{\kappa_x}(2\tilde W-y)\right]^*\, .
\end{aligned}
\end{equation}

It can be shown~\cite{PT} that there is no degeneracy when $\eta=\pm 1$, i.e. when $N_D+1$ is not
divisible by three, while this is not generally true for $\eta=0$.
We also note that the solutions of the problems with $\eta=-1$ and $\eta=+1$ can be mapped into
each other by the relation
\begin{equation}
\vec\varphi_{k_x}^{\,(\eta=-1)}(y)=\sigma_x\vec \varphi_{k_x}^{\,(\eta=1)}(2\tilde W-y)\, .
\end{equation}

Solutions with non-real $\kappa_x^2$ (i.e. $\kappa_x$ with nonzero real and imaginary parts)
can be found in the presence of an external electric field.
Their appearance is related to the existence of exceptional points, i.e. points in which
the operator in Eq.~\eqref{system2} is not diagonalizable (notice that the operator is not self-adjoint).
The existence of nonreal $\kappa_x^2$ values is a manifestation of the
$\mathcal{PT}$ symmetry breaking in the system~\cite{PT}.

The methods described in the previous section [in particular, scheme \eqref{e:imp}] could also  be used to solve
problem (\ref{system2}). However, we
did not observe any significant efficiency gain with respect to
the discretization of the original differential problem.

In the Appendix we show that Eq.~\eqref{system2} can
be recast into the form of a complex second-order differential equation for a scalar unknown function.
In that form the discretization in direct space is free from fermion doubling
effects. We did not
actively investigate its numerical solution, since the method
that we are going to
describe in the next section
turns out to be much more efficient than the direct space ones~\cite{BF:2013}.

\section{Solution in the reciprocal space}
\label{s:fourier}

A key feature of the discretization in the direct space (discussed in Sec.~\ref{s:direct}) was the
representation of the derivative.

Let us consider a uniform grid with node spacing $a$ and the $n$-point discretization on it of the first
derivative $\phi^{\prime}(p)$ of a generic function $\phi$, computed at point $p$.
This discretization is constructed by Taylor expanding $\phi(p+i a)$ for different values of the integer $i$
and finding the linear combination of these expansions that is equal to $a\phi^{\prime}(p)$ up to
$a^{n}\phi^{(n)}(p)$ corrections. If the function $\phi$ is smooth, the discretization error of the
$n-$point derivative is thus $\mathcal{O}(a^{n-1})$ for every
$n$. However, if the $\alpha$th order
derivative of $\phi$ in $p$ is discontinuous, we cannot improve the
precision of the discretization for $n>\alpha$ without introducing coefficients that depend on the specific function $\phi$ itself.
As a consequence, the discretization error of a generic $n-$point discretization of the derivative scales as
\begin{equation}\label{eq:errordis}
|\phi^\prime (p)-\phi^{\prime}_{(N_y)}(p)|\gtrsim \mathcal O(N_y^{1-\min(\alpha,n)})\, ,
\end{equation}
where $\phi^{\prime}_{(N_y)}(p)$ is the $n$-point discretization of the first derivative on a grid
with $N_y$ points.

In the specific case~\eqref{system2}, it is enough that the first derivative of the potential is nonzero
at the boundaries (i.e., the external electric field has a nonzero transverse component at the
edges) for the second derivative of the eigenfunctions to be discontinuous at $y=0,\tilde W$
[see Eq.~\eqref{h}];
if so, the accuracy of the approximation is
independent of $n$ for $n\geq 2$,
and Eq.~\eqref{eq:errordis} represents a very severe limitation
both for the precision and for the efficiency of the numerical solution.

The Fourier methods are better behaved in this respect. While the direct space
methods involve a global distortion of the dispersion relation,
in the Fourier case the derivative is exactly reproduced for
the frequencies lower than the cutoff.
The fermion doubling problem is
absent, as one can argue by tracing back its origin to the periodicity of the wave function across the Brillouin
zone induced by the space discretization (a simple topological argument is given in
Ref.~\cite{smilga}, \S 13.1).

Since both $h(y)$ and $\exp[i\eta \pi y/(3\tilde W)]\vec\varphi(y)$ (separately) assume the same value at 0 and $2 \tilde W$,
they can be extended by periodicity with period $2 \tilde W$ without introducing discontinuities.
We define their Fourier coefficients $h_\ell\equiv h_{-\ell}$ and $\vec a_m$ by
\begin{equation}\label{expansion1}
\begin{aligned}
\displaystyle
h(y) &=\sum_{\ell=-\infty}^{\infty} h_\ell e^{i \pi \ell y /\tilde W} \, , \\
\displaystyle
\vec\varphi(y) &=
\sum_{m=-\infty}^{\infty}\vec a_m e^{i \pi(m-\eta/3) y /\tilde W} \, .
\end{aligned}
\end{equation}
We substitute these expressions in the differential equation of \eqref{system2} and then project
onto the exponential functions $e^{i \pi(n-\eta/3) y /\tilde W}$; for the
generic index $n$ we obtain
\begin{equation}\label{eqfourier}
\sum_{m=-\infty}^{+\infty}
\left[i \frac{\pi}{\tilde W}\left(n-\frac{\eta}{3}\right)\sigma_z \delta_{n,m}+
h_{n-m}\,\sigma_x \right]\vec a_{m}=-\kappa_x \vec a_n \, ,
\end{equation}
where $\delta_{n, m}$ is the Kronecker delta function. These equations are still exact and can be
rewritten in the matrix form
\begin{equation}\label{eq:M}
M \vec{a}=-\kappa_x\vec{a}
\end{equation}
where $M$ is a structured infinite matrix whose $2\times 2$ block is given by
\begin{equation}\label{eq:M2x2}
M_{n,m}=P_n \delta_{n,m}+Q_{n,m}
\end{equation}
with
\begin{equation}
P_n=i \frac{\pi}{\tilde W}\left(n-\frac{\eta}{3}\right)\sigma_z \, , \quad
Q_{n,m}=h_{n-m}\,\sigma_x \, .
\end{equation}

While the weight of the diagonal blocks $P_n$ increases with $|n|$,
Parseval's theorem~\cite{korner1988}
ensures that if $h(y)$ is square
integrable the contribution of the blocks $Q_{n,m}$ vanishes for large values
of $|n-m|$,
i.e., sufficiently far from the principal diagonal of the matrix $M$.
Actually, for the regular potentials for which the envelope
function approach gives reliable results, the hypothesis of square integrability
of $h(y)$ is a very weak assumption.

If we consider a sufficiently large positive integer $D$ such that
\begin{equation}
\frac{\pi D}{\tilde W}\gg \max_j |h_j| \, ,
\end{equation}
the matrix $M_0$, with blocks
$[M_0]_{n,m}=M_{n, m}$ and $|n|,|m|\leq D$, contains the main information about
the slowly varying solutions of \eqref{eqfourier}: the Fourier coefficients that describe the low-frequency components of the spectrum
are well approximated by those of the truncated problem~\cite{BF:2013};
the high-frequency components
are instead negligible for the slow varying solutions.
The finite dimensional problem
$(M_0 +\kappa_x I)\vec v_0=0$ is not affected by doubling, since it is just a
truncation of the original problem \eqref{eqfourier}.

The eigenvalues of $M_0$ are accurate estimates of the longitudinal wave
vectors $\kappa_x$. Each eigenfunction $\vec\varphi(y)$ can be reconstructed,
using the corresponding eigenvector $\vec{a}_0$ of $M_0$, as
\begin{equation}\label{antitransform}
\vec\varphi(y)\approx \vec\varphi_D(y)\equiv \sum_{\mu=-D}^{D}[\vec a_0]_\mu \,
e^{i \pi(\mu-\eta/3) y/\tilde W} \, .
\end{equation}
Using (\ref{varphi}) and (\ref{varphik}), the transverse components $\Phi$ of the envelope functions are given by ($\beta=A,B$)
\begin{equation}\label{Phi-an}
\begin{aligned}
\displaystyle
\Phi^{\vec K}_\beta (y) &\approx
\sum_{\mu=-D}^{D} [a_0^\beta]_\mu \,e^{i \pi (\mu-\eta/3) y /\tilde W}\, ,\\
\displaystyle
\Phi^{\vec K'}_\beta (y) &\approx-i
\sum_{\mu=-D}^{D} [a_0^\beta]_\mu \,e^{-i \pi (\mu-\eta/3) y /\tilde W}\, .
\end{aligned}
\end{equation}
Finally, from \eqref{psia} and \eqref{Phi} it follows that
\begin{equation} \label{combisin}
\psi_\beta (x,y)=2\,i \sum_{\mu=-D}^D
\left\{ [a_0^\beta]_\mu \,\sin[(\mu-n_0) \pi y / \tilde W ]\right\}
e^{i \kappa_x x} \, ,
\end{equation}
where  $n_0$ has been defined in \eqref{n0}.

We explicitly note that, from the numerical point of view, all these computations
strongly benefit from the use of optimized fast Fourier transform routines
for the calculation of the Fourier series.

\section{Numerical efficiency: comparison among methods}\label{s:comparison}

In this section we compare the numerical efficiency of the methods
introduced so far, performing an analysis of the convergence rate
for several test cases.

We are going to compare three main strategies:
\begin{enumerate}
\item[$(S)$] method \eqref{e:imp} of Sec.~\ref{s:direct}
\item[$(S_p)$] method \eqref{e:imp} of Sec.~\ref{s:direct} applied to the periodic problem \eqref{system2}
\item[$(F)$] the Fourier method described in Sec.~\ref{s:fourier}.
\end{enumerate}

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{pot.eps}
\caption{(Color online) Plot of the potentials used to test the
numerical efficiency of the
methods: a step potential (solid curve), a Lorentzian potential
(dashed curve), and a parabolic potential
(dotted curve).}
\label{pot}
\end{center}
\end{figure}

We consider a nanoribbon composed  of \mbox{$N_D=4065$} dimer lines (corresponding to $\eta=1$ and to an effective
width \mbox{$\tilde{W}\approx 500\,\mathrm{nm}$}), with the following potentials (shown in Fig.~\ref{pot}):
\begin{itemize}
\item[(1)] Step potential
\begin{equation}
U(y)=\left\{\begin{array}{ll}
0\,\mathrm{eV} \qquad & y\le 200\,\mathrm{nm} \\
0.2\,\mathrm{eV} & y>200\,\mathrm{nm}
\end{array}\right.
\end{equation}
\item[(2)] Lorentzian potential
\begin{equation}
U(y)=A\frac{\Gamma/2}{(y-y_0)^2+(\Gamma/2)^2}
\end{equation}
with $y_0=200\,\mathrm{nm}$, $\Gamma=100\,\mathrm{nm}$ and $A=10\,\mathrm{eV}\,\mathrm{nm}$
\item[(3)] Parabolic potential
\begin{equation}
U(y)=\bar{A}(y-\bar{y})^2
\end{equation}
with $\bar{y}=250\,\mathrm{nm}$ and $\bar{A}=0.2\,\mathrm{eV}/(250 \,\mathrm{nm})^2$\, .
\end{itemize}
We set the electron injection energy to $E=0.1\,\mathrm{eV}$ and study the scaling of the eigenvalue precision as a
function of the execution time on an Intel Xeon CPU E5420 2.50GHz processor. Diagonalization is performed
by means of standard LAPACK routines. In case $(F)$, the coefficients $h_n$ are computed
on an extremely fine grid, independent of the dimension $D$ of the truncated problem, without introducing any
sizable overhead.

Figures refer to the maximum real eigenvalues. We did not find
significant differences in the behavior of the other eigenvalues.
However, eigenvalues
associated with larger $\kappa_x$ values converge faster to the corresponding
eigenvalues of the original problem~\eqref{system2}.
This can be interpreted as a consequence of the fact that large $\kappa_x$ values correspond to a small kinetic energy
in the transverse direction
(the total energy is constant), i.e., to transverse modes with large wavelength, which are less sensitive to the
discretization or the frequency cutoff.

In Figs.~\ref{ssp} and \ref{fsp} we report the relative error on the largest
real eigenvalue for the two spatial approaches $S$ and $S_p$
(in Fig.~\ref{ssp}), and the methods $S_p$ and $F$ (in Fig.~\ref{fsp}) as
a function of the execution time. The data points for different
values of the execution time have been obtained by varying the discretization
step in the case of the spatial methods, and varying the number of considered
Fourier components in the case of the Fourier methods (the smaller the
discretization step or the greater the number of Fourier components, the
larger the execution time).
Method $(S_p)$ is slightly more efficient than method $(S)$, probably due to the better block structure of the
discretization matrix. However, in all the cases we have studied the Fourier methods largely outperform the direct space
ones, often by several orders of magnitude.

Figure \ref{fsp} also shows that the convergence of the
Fourier method
is strongly dependent on the shape and the analytic properties of the
potential, which
influence the number of Fourier coefficients needed to properly expand
the eigenfunctions (and, in turn, the size of the matrices to be
diagonalized).

The better performance of Fourier methods with respect to direct space ones has been recently noticed also
in Ref.~\cite{wolff2012}, where a Schr\"odinger equation with position dependent mass is considered.
The authors used the following ``Fourier-inspired'' discretization for the derivative ($\mathcal{F}$ is the discrete
Fourier transform and $k$ is the reciprocal space variable)
\begin{equation}
\frac{d}{d x}\longrightarrow \mathcal{F}^{-1}\,k\,\mathcal{F}
\end{equation}
and reported a convergence rate exponentially fast in the number of points of the discrete Fourier transform.

If we compare the numerical errors of the methods as a function of the size of the matrices involved in the analysis, we conclude that
the errors deriving from the finite-difference discretization of the
derivatives (and from the resulting distortion of the
dispersion relation) turn out to be much larger than those related
to the cutoff of the high-frequency Fourier components in the reciprocal
space approach.
This shows that, for the type of potentials we are interested in, the Fourier method is drastically
more efficient than the others.

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{ssp.eps}
\caption{(Color online) Relative error on the largest real
eigenvalue as a function of execution time, when using the spatial methods
$S$ and $S_p$ (see text for the definition of abbreviations). Different
execution times correspond to different discretization steps.}
\label{ssp}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{fsp.eps}
\caption{(Color online) Relative error on the largest real eigenvalue as a
function of the execution time. Comparison between the Fourier method and the
spatial $S_p$ one (see text for the definition of abbreviations). In the case
of the parabolic potential the precision very quickly reaches the machine
precision. Different execution times correspond to different values of the
discretization step for the $S_p$ method and to different numbers of
components for the Fourier method.}
\label{fsp}
\end{center}
\end{figure}

\section{Solution of the transport problem}\label{s:transport}

In the previous sections we have described numerical methods to compute the eigenvalues and the eigenfunctions
of the Dirac equation in a longitudinally invariant ribbon.
The total wave functions on the two sublattices $\psi_\beta (\vec r)$ ($\beta=A,B$) are a linear
combination of modes $\psi_{\beta i} (\vec r)$ of the form:
\begin{equation}\label{chi}
\begin{aligned}
\psi_{\beta i} (\vec r) &=
\left[e^{-iKy} \Phi^{\vec K}_{\beta i} (y)-ie^{iKy}
\Phi^{\vec K'}_{\beta i} (y)\right]
e^{i \kappa_{xi} x} \\
&\equiv\chi_{\beta i} (y)\,e^{i \kappa_{xi} x}\ .
\end{aligned}
\end{equation}

In the case of a general potential $U(\vec{r}\,)$
we divide the ribbon into a series of transverse slices, in such a way that within each slice
the potential is approximately independent of $x$.
For each slice we can then apply the previously discussed methods to estimate the modes $\psi_{\beta i}$ and their
longitudinal momenta $\kappa_x$.
At the interfaces between adjacent slices we have to enforce the continuity of the total wave
function (the $\vec{k}\cdot\vec{p}$ approximation is reliable only if the potential varies slowly
on the lattice scale, thus no $\delta$-type potentials are allowed).
We remark here a difference with respect to the standard
Schr\"odinger case: the Dirac equation is a first-order differential
equation, thus we do not have to impose the continuity of  the normal
derivative of the wave function.

Since the atomic orbitals in Eq.~\eqref{total} are strongly localized,
enforcing the continuity of the total wave function
amounts to imposing the continuity of the wave functions on both inequivalent
sublattices separately.
Moreover, since the functions $\psi_\beta (\vec r\,)$
have Fourier components localized around the two inequivalent Dirac
points (which are significantly separated from each other), the continuity of
the functions $\psi_\beta (\vec r\,)$ implies also the continuity of the
envelope functions $F$.

Integrating the probability current density
in the $x$ direction~\cite{kp2011}
\begin{equation}\label{jx}
J_x(y) =v_F\,[
\vec{F}_{\vec{K}}^{\dag}(y)\sigma_x\vec{F}_{\vec{K}}(y)+
\vec{F}_{\vec{K}'}^{\dag}(y)\sigma_x\vec{F}_{\vec{K}'}(y)]
\end{equation}
over the transverse section, and using Eq.~(\ref{Phi-an}),
we can express the longitudinal probability current as follows:
\begin{equation}
\begin{aligned}
I_x &=\int_0^{\tilde W} J_x(y)d y
=4v_F \tilde W\,\Re\left[\sum_{n=-D}^D {(a_n^A)}^* a_n^B\right] \\
&=v_F\,\int_0^{2\tilde{W}} \vec{\varphi}(y)^{\dag}\sigma_x\vec{\varphi}(y)d y \, .
\end{aligned}
\end{equation}

From Eq.~\eqref{eq:Z2} we deduce that $(I_x)_{-\kappa_x}=-(I_x)_{\kappa_x}$,
$(I_x)_{\kappa_x^*}=(I_x)_{\kappa_x}$ and $(I_x)_{-\kappa_x^*}=-(I_x)_{\kappa_x}$.
In particular, if $\kappa_x$ is purely imaginary, we have
$I_x=0$, i.e., modes with purely imaginary eigenvalues do not carry
current. This is in general not true for eigenvalues that have at the same time
a nonzero real and imaginary part.
We classify the modes as right-moving or left-moving, depending on whether they have a
positive or negative longitudinal probability current $I_x$. We also extend the definition of
right-moving (left-moving) to the modes with $I_x=0$ and $\Im(\kappa_x)>0$ ($\Im(\kappa_x)<0$).

In detail, in our simulation code we order the modes
on the basis of the value of the corresponding $\kappa_x$. We first
consider the real $\kappa_x$ (arranged in order of decreasing modulus),
then the complex ones, and finally the purely imaginary ones (sorted in
order of increasing modulus). Since this ordering reflects the expected weight
of the different modes in a transport simulation, in our computations we
consider only the first
$n_{mod}$ right-moving modes and the first $n_{mod}$ left-moving modes of each slice.
Clearly $n_{mod}$ has to be large enough for the final physical result to be insensitive
to its specific value. Moreover, we select the modes in such a way as to preserve the
$Z_2\times Z_2$ symmetry, which means to pick at
the same time the modes with eigenvalues $\kappa_x$, $-\kappa_x$,
$\kappa_x^\ast$, and $-\kappa_x^\ast$.

Let us now sketch the basic steps to compute the scattering matrix
for a single discontinuity of the potential at the interface between adjacent
slices.
We denote by $l/r$ the modes on the left/right of the discontinuity and
by $+/-$ the right/left-moving modes.
We use the index $i$ to denote the mode impinging on the discontinuity, e.g., from the left.
The wave function $\psi_\beta (\vec r)$ on the left side can be written as
\begin{equation}\label{eq:left}
\chi^{l+}_{\beta i} (y)\,e^{i \kappa^{l+}_{xi} (x_{dis}-x_{in})}
+\sum_n r_{ni} \chi^{l-}_{\beta n} (y)\,e^{i \kappa^{l-}_{xn} (x_{dis}-x_{in})} \ ,
\end{equation}
while on the right side it can  be expressed in the form
\begin{equation}\label{eq:right}
\sum_n t_{ni} \chi^{r+}_{\beta n} (y)\,e^{i \kappa^{r+}_{xn} (x_{dis}-x_{out})} \ .
\end{equation}
Here $x_{in}$ and $x_{out}$ are the longitudinal positions of the
boundaries of the considered scattering region, $x_{dis}$ is the position of the discontinuity, while
$r_{ni}$ and $t_{ni}$ are the reflection and transmission coefficients.
By continuity, functions \eqref{eq:left} and \eqref{eq:right} must be equal.
An analogous relation can be established for a mode injected from
the right.
These continuity relations have to be enforced for both sublattices and for
all the $2\,n_{mod}$ modes impinging from the left and from the right.

In order to evaluate all the $4\,n_{mod}^2$ reflection and transmission coefficients,
we can project the $4\,n_{mod}$ continuity constraints onto a set of functions chosen
in such a way as to obtain the correct number of independent equations.
From Eq.~\eqref{combisin} we have
\begin{equation}
\chi_{\beta i} (y)=2\,i \sum_{n=-D}^D
\left\{ a_n^{\beta i}\,\sin[(n-n_0) \pi y / \tilde W ]\right\} \, ,
\end{equation}
hence it is natural to project each continuity relation
on the set of $n_{mod}$ functions
\begin{equation}\label{sinfunc}
S_j (y)=\sin\left((j-n_0) \pi y / \tilde W \right)\ ,
\end{equation}
for $j=-(n_{mod}-1)/2 , \ldots, (n_{mod}-1)/2$
(for the sake of simplicity $n_{mod}$ is assumed odd).
Since we consider values of $n_{mod}$ such
that $(n_{mod}-1)/2 < n_0$, these functions are linearly independent.
It is simple to show that the matrix elements are
\begin{equation}
\langle S_j (y) | \chi_{\beta i} (y) \rangle =
\int_0^{\tilde W} S_j^* (y) \chi_{\beta i} (y) d y=
i \tilde W a_j^{\beta i} \, ;
\label{samew}
\end{equation}
thus all the computations can be performed in the reciprocal space,
avoiding the evaluation of the sums in Eq.~\eqref{Phi-an}.

Once the scattering matrices corresponding to the various interfaces have been computed,
they can be composed according to the standard procedure (see, e.g., Ref.~\cite{datta})
to obtain the total scattering matrix $S$ of the ribbon:
\begin{equation}
S=\left(\begin{array}{cc}
r & \tilde{t} \\
t & \tilde{r}
\end{array}\right)
\end{equation}
Here $r$ and $t$ are the reflection and transmission matrices for the modes
impinging from the left, $\tilde{r}$ and $\tilde{t}$ the corresponding
matrices for the modes
impinging from the right.

For practical purposes it is convenient to introduce the current form $S'$ of the scattering
matrix:
\begin{equation}
S'=\left(\begin{array}{cc}
r' & \tilde{t}' \\
t' & \tilde{r}'
\end{array}\right) \, ,
\end{equation}
relating the ``current amplitudes'' instead of the ``wave amplitudes'' of
the modes~\cite{datta}.
This matrix involves only the modes with $I_x\neq 0$ and its elements are
given by $s'_{nm}=s_{nm}\sqrt{|I_{xn}|/|I_{xm}|}$, with $s=r,t,\tilde r,\tilde t$.
As a result of current conservation, it can be shown that $S'$ is unitary
(see, e.g., Ref.~\cite{datta}), which is a useful check to be performed at the end of the computations.
In all our simulations we checked that numerical violations of the unitarity relation are less
than $10^{-13}$.

\begin{figure}[t]
\begin{center}
\includegraphics[width=8cm]{fig5.eps}
\caption{(Color online) (a) Map of the potential in the nanoribbon,
given by a superposition of Lorentzian functions.
(b) Normalized conductance as a function of the injection energy,
obtained within the envelope function and the
nearest-neighbor semiempirical tight-binding approximations.}
\label{fig5}
\end{center}
\end{figure}

\begin{figure}[t]
\begin{center}
\includegraphics[width=8cm]{fig6.eps}
\caption{(Color online) (a) Map of the potential in the nanoribbon,
represented by a tilted barrier with a Lorentzian profile.
\hfill\break
(b) Normalized conductance as a function of the injection energy,
obtained within the envelope function and the
nearest-neighbor semiempirical tight-binding approximations.}
\label{fig6}
\end{center}
\end{figure}

From $S'$ we can compute the conductance of the ribbon by means of the
Landauer-B\"uttiker formula
\begin{equation}
G=\frac{2\,e^2}{h}\sum_{n,m} |t'_{nm}|^2\ ,
\label{landauer}
\end{equation}
where the sum runs over the modes with $I_x \ne 0$ in the first and last
transverse regions of the ribbon.

In studies of unconfined graphene or of ribbons with large aspect ratio, it is
usual to assume periodic instead of Dirichlet boundary conditions;
the two Dirac points are then completely decoupled, and it is customary to
solve the Dirac equation for just one valley and use a factor of four instead
of two in Eq.~\eqref{landauer}. The physical Dirichlet boundary
conditions introduce instead a coupling between the two inequivalent Dirac points,
requiring the use of the more general formulation~\eqref{landauer}.

For validation purposes, we have performed a transmission calculation for a structure that is small enough
to allow also a treatment with a standard tight-binding code
(in particular we have used NanoTCAD ViDES)~\cite{vides1,vides2}.
We considered an armchair nanoribbon with $60$ dimer lines
($\approx 7.5$~nm wide), in the presence of two different
realistic potential profiles.

In the first test case the electrostatic potential is a superposition of
Lorentzian functions,
and can schematically represent the effect on the ribbon of charged impurities
located in the substrate on which graphene lies.
In particular, we consider five Lorentzian functions, with a peak amplitude
of $0.5$~eV and a half-width at half-maximum equal to $0.64$~nm
(see Fig.~\ref{fig5}(a)).
In Fig.~\ref{fig5}(b) we show the computed behavior of the
conductance $G$ (in units of $2 e^2/h$)
as a function of the injection energy, together with the corresponding
results obtained with ViDES.
We observe a very good agreement between
the two different approaches in the low injection energy regime $E_{in}\lesssim 0.5$~eV, i.e., the one in which
the envelope function method can be safely applied. For larger energies,
the simple Dirac equation, which represents only a
first-order $\vec k\cdot \vec p$ approximation,
does not appropriately describe the physics of graphene any more, and thus
discrepancies between
the two results appear, even though the qualitative behavior
of the conductance is well
reproduced for all the explored injection energies.

The other potential profile we have considered is
a barrier, tilted with respect to the ribbon edges and with a Lorentzian
profile, with a peak amplitude of $0.625$~eV and a half-width at
half-maximum equal to $2$~nm [see Fig.~\ref{fig6}(a)]. The ribbon has a
width of 60 dimer lines, as in the previous example. This
potential can be the representation of the electrostatic effect at
the graphene level of a biased gate or of a line of charge present at a
certain distance from it. In Fig.~\ref{fig6}(b) we report
the behavior of $G$ (in units of $2 e^2/h$) as a function of the electron
energy,
obtained with our envelope-function-based calculation and with the
tight-binding code. Also in this case we notice a good agreement
between the two approaches, especially in the low-energies regime,
in which the $\vec k\cdot \vec p$ approximation is expected
to be more accurate.

Indeed, our numerical analysis is based on the use, in each section with
longitudinally constant potential, of the continuum Dirac
equation \eqref{dirac}. This approximation is valid if the following
conditions are satisfied:
\begin{enumerate}
\item[(1)] The dispersion relation is approximately linear
\item[(2)] The wave function is slowly varying on the scale of the lattice spacing
\item[(3)] The potential is slowly varying on the scale of the lattice spacing.
\end{enumerate}

In the hypothesis of a slowly varying potential~\cite{ashcroft}, the energy of an
electron with wave vector $\vec{\kappa}$ in position $\vec{r}$ can be written
as $E \simeq T(\vec{\kappa})+U(\vec{r}\,)$, where $T$ is the kinetic energy.
Observing the explicit form of $T(\vec{\kappa})$ (i.e., the actual
dispersion relation in the absence of potential
energy, a relationship which can be derived, for
example, with a tight-binding formulation~\cite{castroneto2009, saito}), we see
that the first condition, i.e., the linear approximation for $T(\vec{\kappa})$
($T(\vec{\kappa})\simeq \pm\hbar v_F|\vec \kappa|$), is
valid for $|T(\vec{\kappa})|\lesssim 1~\mathrm{eV}$, which corresponds to
\begin{equation}\label{cond1}
|E-U(\vec{r}\,)|\lesssim 1~\mathrm{eV}\,.
\end{equation}

The second condition can be expressed in the form
$\lambda=2\pi/|\vec \kappa|\gg a$
(where $\lambda$ is the electron wavelength). Exploiting the relation
(valid under the previous approximations)
$E\simeq \pm\hbar v_F|\vec \kappa|+U(\vec r)$,
this inequality translates into the condition
\begin{equation}\label{cond2}
|E-U(\vec r)| \ll (2\pi\hbar v_F)/a\approx 15~\mathrm{eV}\, ,
\end{equation}
which is clearly weaker than Eq.~\eqref{cond1}.

Instead, the requirement on the smoothness of the potential (the third
condition) introduces a limitation on the derivative of $U$ along $y$
(the only spatial variable along which $U$ varies within each section). If
we require the variation of the potential energy over the lattice constant
$a$ to be negligible with respect to $\gamma_0$ [which represents an
order of magnitude of the energies involved, since, in a first
approximation~\cite{saito}, $T(\vec{\kappa})$ has values between
$-3\gamma_0$ and $3\gamma_0$], this further constraint can be expressed as
\begin{equation}\label{cond3}
|\partial U/\partial y|\ll \gamma_0/a\approx 11~\mathrm{eV}/\mathrm{nm}\ .
\end{equation}
In the previous numerical examples, Eq.~\eqref{cond3} was always well
satisfied; thus the only
limit to the application of the Dirac equation was the condition
\eqref{cond1} on the injection
energy. Indeed, as previously noted, for $E\gtrsim 0.5$~eV the
continuum and the tight-binding results
are not in as good an agreement as in the low-energy region.

\section{Conclusions}\label{s:conclusions}

We have presented a numerically efficient approach, including physical boundary conditions, for the evaluation of transport
properties of graphene devices for which the
application of atomistic techniques is computationally prohibitive.
We have focused on ribbons with armchair edges,
which we modeled within a continuum, envelope function approximation.

For the computation of the transmission we have adopted a recursive scattering
matrix approach, which requires the solution of a collection of Dirac equations
in the presence of longitudinally constant potentials.
We have shown that a reciprocal space approach is largely
preferable with respect to
the more commonly adopted finite difference methods, since it can reduce
the computational cost of the procedure by orders of magnitude.

We have compared our results for structures small enough to allow an atomistic
simulation with those obtained by means of tight-binding techniques, finding good agreement
within the range of validity of the $\vec k\cdot\vec p$ approximation.

\section*{Acknowledgments}

It is a pleasure to thank Paolo Pintus for useful discussions.

\newpage

redpajama_set_name: RedPajamaArXiv


Constitution, and every immoral policy with links to 9/11.
Whenever a policy of President Bush's is questioned--he ties it to 9/11. Invaded the wrong country? 9/11. Tortured prisoners? 9/11. Wiretapped American citizens? 9/11. Sold our ports to the UAE? 9/11. And now: blocked a ceasefire agreement in the UN long enough for Israel and Hezbollah to wage war? 9/11.
For President Bush, it's all about 9/11. Translation? If you question what I do, you are a traitor.
How should Democrats respond? By talking about a "Smart and stronger" foreign policy.
The conflict in Lebanon is part of a broader struggle between freedom and terror that is unfolding across the region. For decades, American policy sought to achieve peace in the Middle East by promoting stability in the Middle East. Yet the lack of freedom in the region meant anger and resentment grew, radicalism thrived and terrorists found willing recruits. We saw the consequences on September the 11th, 2001, when terrorists brought death and destruction to our country, killing nearly 3,000 of our citizens.
Did you catch that? Lebanon is now part of the 'broader struggle between freedom and terror.' Deja vu all over again.
Amazing how Bush can say this stuff with a straight face. For the rest of us, we did not see any 'freedom unfolding' in the past few weeks in Lebanon. We saw buildings collapsing. We saw ambulances exploding. We saw children dying. And we saw a conflict between Israel and Hezbollah escalating.
For everyone in the world except George W. Bush, the situation in the Middle East is much worse after the conflict between Hezbollah and Israel flared up into a full scale shock and awe war.
Israelis are less at ease, not more. Hezbollah is more emboldened, not less. The United States is seen with more suspicion, not less. More people died, not less. Ties between Iran, Syria and rockets being launched at civilians are stronger, not weaker.
Bush's decision to block cease fire--to use Condi Rice and John Bolton to prevent diplomacy from happening--was a stupid and weak policy that has resulted in more death and less stability.
A month ago, Iraq was in full meltdown. Now, Lebanon is a smouldering tinderbox on top of it all.
In the midst of all this, Bush has once again claimed 'victory' where there is none.
Victory? The war in Lebanon has turned Hezbollah into the rebel group of the month. They stood up to Israel and the United States and did not roll over and die. Congratulations, President Bush. Your 'war on terror' in Lebanon did little more than kill thousands and thousands of innocent people, destroy the infrastructure of yet another country, and turn a radical Islamic group hell-bent on destroying Israel into a league of heroes.
And now, on top of all this, Lebanon will be yet another location in the Middle East where people wonder if they will have electricity and water on a regular basis--and blame it all on the United States.
Our policy must be smart and strong. And what is smart? Working with other people? What is strong? Using diplomacy instead of military invasion.
President Bush would like us all to believe that America is strongest when we launch invasions and occupy foreign countries. He is wrong. America is strongest when we work with other nations. American is strongest when we persuade our allies to join us. America is strongest when we stand with our friends in the Middle East, not when we stand in the way of cease fires.
President Bush thinks it is smart to encourage an ally to launch a military invasion. He is wrong. We are already seeing endless trouble as a result of what we did in Iraq. A smart person would learn from those mistakes and try something different. Not George Bush. He sees the problems we have in Iraq and says: I'll do the same thing in Lebanon, but this time it will be better because I say so. Smart policy in the Middle East brings results, not destruction. Smart policy in the Middle East creates new and closer allies, not more people who hate Americans. Smart policy in the Middle East tries every possible avenue to avoid war, not every possible excuse to launch another one. Americans are tired of George Bush and his dumb policies.
Here the logic does not look backwards and blame Bush for what happened in Iraq, but looks forward and says: Let's try something new, now.
One thing is or certain, President Bush will keep in talking about 9/11 as election season heats up. He's tried it several times before and he will try it again.
This time it will not work. Americans are tired of a weak and stupid foreign policy.


text:
eHarmony Balks At Chemistry.com Commercials
By Andrew Belonsky May 14, 2007 at 9:05am · Updated on April 27, 2018
Online dating service eHarmony isn't taking Chemisty.com's advertorial attacks lying down. The company – which got its start under James Dobson's Bible-thumping thumb – has hired famed lawyer Lanny J. Davis – who helped Clinton during that whole Lewisnky blow job brouhaha.
In an effort to save his client's hide, Davis has asked NBC and People to stop running Chemistry.com's adverts, which highlight eHarmony's irrational, allegedly discriminatory practices.
As you may recall, one commercial features a man saying eHarmony rejected his application because he's gay.
eHarmony alleges that Chemistry.com – which is owned by Barry Diller's Match.com – and their commercials trivialize – even scandalize – eHarmony's strict guidelines. Washington Post reports:
[eHarmony's] concerned that its rival's ads suggest that eHarmony is being arbitrary — or worse, racially and religiously discriminatory — in turning people away.
The biggest reason for rejection, it says, is that the applicant is married. Stunningly, nearly one-third of the company's rejects (30 percent) fell into this category. Others are blocked because they're younger than the minimum application age of 21 (27 percent) or because the applicant gives inconsistent answers (9 percent), based on responses to eHarmony's 258-question application.
eHarmony big wig Greg Waldorf also claims that eHarmony does not discriminate against the gays. They just haven't figured out how to work them into its romantic, marriage-inclinded web:
Waldorf says eHarmony's matching system is based on psychological research about heterosexual relationships. Because it doesn't have similar data on gay people, he says, the company isn't confident that it can offer successful matches to same-sex couples. "I'm not saying anything precludes us from going into the same-sex market in the future," he says, "but it's not a service we offer now."
Psychology being used to exclude gays? eHarmony sure is cutting edge.
For the record, People has refused to get involved in this lover's spat.
Bill Clinton Discrimination eHarmony News Relationships
Subscribe to Queerty for a daily dose of #billclinton #discrimination #eharmony stories and more
nystudman
Can't understand (seriously!) why a gay man would bother trying to get on this site when there are 999,999 gay dating sites out there. I mean, yes, they're right-wing dingbats. but it's a private company and it's not like we don't have any other choices. Seems to be a manufactured "controversy." (BTW, how many str8 people have complained about not being able to get onto Manhunt or Gay.com?)
May 14, 2007 at 12:05pm
el polacko
my guess is that str8 guys would be welcome at gay sites, as long as they play nice. i'd be fine with eharmony being a str8s-only site if they advertised themselves that way.
otherwise, i'm fine with seeing them exposed as the dobson-run loony bin that is behind the mask of the kindly old spokesperson. it's main rival, match.com, has always tried to find mates of all persuasions.
May 14, 2007 at 5:05pm
bloomsburyboy
for the record, gay.com has "straight" as an option on its personals profiles.
and this article's focus is not that eharmony is discriminating (they're a private company providing an optional service, so they can do pretty much whatever they please) but rather that they're pissed/scared that the competition has a better marketing campaign than they do. and are trying to unfairly stifle chemistry.com.
May 15, 2007 at 2:05am
Indie Movie
Wildhood
We're All Going To The World's Fair
Please Baby Please
Girl Picture

redpajama_set_name: RedPajamaCommonCrawl


<Support the #11 Bus with a Pub Crawl!
We are planning a carpool please call or email info@forwardchicago.org to let us know you would like to come.
Please join us August 5th for a timely conversation about ageism in America, and what can be done to eradicate it.
ASA's Board President Bob Blancato, and ASA President/CEO Bob Stein, will present results from the 2015 research "Gauging Aging: Mapping the Gaps Between Expert and Public Understandings of Aging in America." This report lays the groundwork for a larger effort to develop an evidence-based narrative about aging in our country and the needs and contributions of older adults. By comparing experts' views to those of average Americans, the report detailed a set of communications challenges to efforts to elevate public support for policies and programs that promote the well-being of older adults.
This session will also detail a partnership of eight of the nation's leading aging-focused organizations to create a better public understanding of older adults' needs and contributions to society—and subsequently to improve the lives of all people as they age.
Finally, they will discuss public policies related to addressing ageism.
Join us in this lively and provocative discussion about ageism—an attitude that permeates our society and inhibits optimal aging and quality of life for America's older adults – and what we can do, both collectively and individually to fight it!


text:
Home/News/Technology/Toyota Offers Everything but a Fix for bZ4X EV Recall
Toyota Offers Everything but a Fix for bZ4X EV Recall
United News Team August 4, 2022
A few months ago, Toyota released its first all-electric vehicle, the bZ4X EV, then quickly announced a recall over the wheels potentially falling off. And while the problem sounds easily fixable, with the hub bolts loosening over time, they've yet to offer a solution.
Toyota announced the recall in June, instructing that all owners stop driving the new bZ4X electric vehicle immediately. Thankfully the car is still relatively new, and the problem only affects around 2,700 vehicles on the road.
The electric motors, battery cells, and everything else electric is fine. This isn't some weird problem with new technology. Instead, this is a wheel issue. We expected a quick fix from one of the biggest auto manufacturers on the planet, but that's not the case.
Now, Toyota is offering owners just about everything but a fix. The automaker recently started sending out letters to bZ4X owners that include a wide array of remedies but not a solution. Toyota is going above and beyond to help, but we're still confused as to why they haven't come out with a fix.
In the letter, Toyota corporate offers EV owners a free loaner vehicle, reimbursement for fuel costs, allowing owners to park and store the bZ4X at dealerships for free, a $5,000 credit for those that paid in full, or even a complete vehicle buyback. Yes, Toyota will buy the EV back from those that can't drive it due to the recall.
In case you didn't notice, the recall isn't going very well. Considering only 2,700 vehicles were affected, we all expected a quick and easy fix. And while yes, this is a physical recall where dealers need to replace the hub bolts, not a software update, it's still amazing that Toyota doesn't have any answers yet.
Toyota produces millions of vehicles yearly with wheels that don't fall off. So, it's a bit of a head-scratcher that its latest EV has this problem. Maybe there's more to the situation than we know, but we'll have to wait and see.
Until further notice, Toyota still instructs owners to stop driving the car indefinitely. Furthermore, the company says to contact a Toyota dealership, which will pick it up, transport your vehicle, and store it for free until a fix is available for the recall. In the meantime, Toyota will reach out to each owner to discuss options.
via InsideEVs
fbq('init', '1137093656460433'); fbq('track', 'PageView'); },3000);
bZ4X fix offers recall Toyota
Steam Now Sells a New Start Menu for Windows 11
Biden administration declares monkeypox a public health emergency
Why HBO's 'The Last of Us' rewrote the love story of Bill and Frank
Perspective | 'The Last of Us' tells a new but familiar queer love story
Review | 'The Last of Us' recap: A kinder, more loving tale for Bill and Frank
Check out these emulated calculators at the Internet Archive
The latest Super Mario Bros. Movie teaser showcases Seth Rogen as Donkey Kong


Akdemir B, Benditt DG. Vagal stimulation: An Evolving Adjunctive Treatment for Cardiac disease. Anatol J Cardiol. 2016 Oct;16(10):804-810. doi: 10.14744/AnatolJCardiol.2016.7129.
Akdemir B, Yarmohammadi H, Alraies MC, Adkisson WO. Approach to Management of Premature Ventricular Contractions: Reassure or refer? Cleveland Clinic Journal of Medicine, 2016 Jul ; Volum 83 , number 7 : doi:10.3949/ccjm.83a.15090.
Dickinson O, Akdemir B, Puppala V, Krishnan B, Detloff B, Sakaguchi S, Adkisson W, Benditt DG . Blunted chronotropic response to hypotension in Cough syncope. JACC Clinical Electrophysiology JACCCEP. 2016;():. doi:10.1016/j.jacep.2016.02.017.
Akdemir B, Vakil K, Benditt DG, Adabag S. Could the AutocaptureTM Function of Permanent Pacemakers Initiate Ventricular Tachyarrhythmias?. Pacing Clin Electrophysiol. 2016 Jan;39(1):92-4. doi: 10.1111/pace.12697. Epub 2015 Aug 24cing Clin Electrophysiol. 2015 Jul 30. doi: 10.1111/pace.12697.
Akdemir B, Krishnan B, Senturk T, Benditt DG. Syncope: Assessment of Risk and An Approach to Evaluation in the Emergency Department and Urgent Care Clinic. Indian Pacing Electrophysiol J. 2015 Jul 29;15(2):103-9. doi: 10.1016/j.ipej.2015.07.005.
Akdemir B, Altekin RE, Kucuk M, Yanıkoglu A, Karakas MS, Aktas A, Demir I, Ermis C. The significance of the left atrial volume index in cardioversion success and its relationship with recurrence in patients with non-valvular atrial fibrillation subjected to electrical cardioversion: study on diagnostic accuracy. The Anatolian Journal of Cardiology. 2013;13(1):18-25.
Lü F, Adkisson WO, Chen T, Akdemir B, Benditt DG. Catheter ablation for long-standing persistent atrial fibrillation in patients who have failed electrical cardioversion.J Cardiovasc Transl Res. 2013 Apr;6(2):278-86.
Karakaş MS, Altekin RE, Ozbek SC, Yanıkoglu A, Akdemir B, Er A, Cilli A, Yalcınkaya S. Determination of Subclinical Atherosclerosis in Obstructive Sleep Apnea Syndrome Patients Without Traditional Risk Factors for Atherosclerosis. Nobel Medicus 2013;9(3):27-34.
Oncel RC, Ucar M, Karakas MS, Akdemir B, Yanikoglu A, Gulcan AR, Altekin RE, Demir I. Relation of Neutrophil/Lymphocyte Ratio With GRACE Risk Score to In-Hospital Cardiac Events in Patients With ST-Elevated Myocardial Infarction. Clin Appl Thromb Hemost. 2013 Sep 27. [Epub ahead of print].
Gulcan AR, Karakas MS, Altekin RE, Ozbek SC, Akdemir B, Yilmaz H. The relationship between subclinical atherosclerosis and Left atrial volume index in patients with metabolic syndrome. Gaziantep Med J 2013;19(1):5-12.
Gulcan AR, Karakas MS, Akdemir B, Ucar M, Altekin RE, Yilmaz H. The relationship between subclinical atherosclerosis and red cell distribution width in patients with metabolic syndrome. Dicle Med J 2013;40(4):579-584.
Karakas MS, Ozbek SC, Akdemir B, Er A, Yanikoglu A, Altekin RE, Baktır AO, Yalcinkaya S, Cilli A. Evaluation of Renal Function In normotensive patients with Obstructive Sleep Apnea Syndrome. Erciyes Med J. 2012;34(3):137-40.
Gulcan AR, Karakas MS, Akdemir B, Ucar M, Altekin RE, Yilmaz H. The relationship between subclinical atherosclerosis and the mean platelet volume in patients with metabolic syndrome. Turkish Society of Cardiology Archive.
Karakas MS, Ozbek SC, Er A, Akdemir B, Altekin RE, Gulcan AR, Yanikoglu A, Kucuk M, Ermis C. Right heart thrombus entrapped in patent foramen ovale with pulmonary embolism in a patient with primary hypercoagulable state. Dicle Med J 2012;39(3):440-444.
Altekin RE, Karakas MS, Yanikoglu A, Ozbek SC, Akdemir B, Demirtas, H., Deger, N. & Cekirge, I. Aortic valve endocarditis and cerebral mycotic aneurysm due to brucellosis. Dec 2011 In : Journal of Cardiology Cases. 4, 3.


Get the latest Manchester United 1968 Special Edition Dream League Soccer – Manchester United 1968 Special Edition Kits and create your own dream team. Manchester United logo is very amazing. You can also check All Manchester United Kits. Manchester United is a very famous Premier League Kits. Size of the kit is 512×512. Manchester United has a very amazing Premier League Kits.


text:
The Mind Renewed
Thinking Christianly in a New World Order
Latest News! Conspiracy Theorists Within US Government Conclude Kennedy Shot as Result of Conspiracy! (HSCA)
HSCA: "Latest News! Conspiracy Theorists Within US Government Conclude Kennedy Shot as Result of Conspiracy!", 1976 -1979.
Hot off the press! Conspiracy theorists working at the highest levels of the US Government say they believe President John F. Kennedy was shot as the result of a conspiracy. In an amazing admission of tin-foil-hat wearing hysteria, they state:
"The Committee believes, on the basis of the evidence available to it, that President John F. Kennedy was probably assassinated as a result of a conspiracy." *
Not content with this level of irrationality, the committee also concludes:
"The committee conducted a three-pronged investigation of conspiracy in the Kennedy assassination. On the basis of extensive scientific analysis and an analysis of the testimony of Dealey Plaza witnesses, the committee found there was a high probability that two gunmen fired at President Kennedy." *
What next? While it's a relief that the Committee didn't go down the route of blaming the Secret Service, FBI or CIA - each one a favourite target for conspiracy theorist paranoia - one cannot help wonder what they'll come out with next! Perhaps their second gunman will turn our to be Buddy Holly! Or the Yeti on a sightseeing trip from the Himalayas! Only time will tell.
Image: "Portrait Photograph, President John F. Kennedy. White House, 07/11/1963" : US National Archives and Records Administration, ARC Identifier 194255, [Public Domain]; resized, cropped.
Copyright © 2021 The Mind Renewed : Thinking Christianly in a New World Order.
Site Map │ Cookie Policy


We've all had dreams of the Inca Trail and walking a path so steeped in history and beauty. But the reality is you have to book months in advance and the experience offers very little in the way of real adventure. If the idea of authentic exploration excites you, consider the Choquequirao trek, the best alternative to Machu Picchu and the Inca Trail.
High on the Quriwayrachina mountain range of Peru lays the Choquequirao archaeological site, and the only way in or out is on foot. Now, this trek is not for the fainthearted, but it's almost unsurpassed in Peru. With sweeping views of the Andes and opportunities to observe local plant and wild life, including the sacred Andean condor, this trek is a must for any serious hiker.
With the country's popularity growing among international tourists, it is becoming difficult to find an authentic, uncrowded adventure in Peru. The Choquequirao trek is just this.
The Choquequirao route takes you through the steep Apurimac Canyon; as you descend down into this plunging valley, and hike back out the steep other side, you will pass through several microclimates that vary with the altitude.
Sunset across the Apurimac Valley in Peru.
Surrounded by lush rainforest and snow-capped peaks, the Choquequirao route has new views at every turn and it should come as no surprise that this mountain-top city is being lauded as the "new Machu Picchu".
Choquequirao, meaning Cradle of Gold, is an Incan city that was built in the same era as Machu Picchu. Thought to be built by Pachacuti, the Incan Emperor who envisioned Machu Picchu, the city was extended by his son, Tupaq Inka Yupanki.
Choquequirao is roughly the same size as Machu Picchu and contains a few different architectural styles. Despite covering six square kilometres, only about a third of it has been reclaimed from the jungle.
Choquequirao was built on a flattened hilltop that was levelled centuries ago and ringed with stones to stop it sliding down the mountain. Though much of the city is inaccessible due to the encroaching jungle, a temple and several other buildings around the main square have been cleared, as well as several terraces.
My favourite part of the site is these terraces. On the west side of the summit, huge stepped terraces have the shapes of llamas built into them with carefully carved white rocks. These are accessible through a clambering jungle path that opens out to breathtaking views. On the east side, below the main site, are the iconic animal shaped terraces, thought to mirror the shape of a guinea pig.
With no people around, and no definitive answers about what everything means or was for, enjoy feeling like a true adventurer exploring a city time forgot.
Currently, the only way to get to the Choquequirao ruins is the four-day trek but the local Peruvian government is planning to build a Choquequirao cable car through the Apurimac Valley. In February 2018 a company was awarded the project, but the start date of construction hasn't been announced yet.
Many people wonder what's the main difference between hiking to Choquequirao and Machu Picchu. In two words: visitor numbers.
The Choquequirao ruins receive less than a dozen visitors a day – that's only 5,800 a year compared to Machu Picchu's 1.2 million. Unfortunately, once the cable car is functioning, the government is expecting over 200,000 visitors a year. This fate is already happening to Kuelap, a huge stone fortress perched high in the mountains near Chachapoyas in northern Peru – it won't take long before Choquequirao is overrun with visitors.
A bromeliad growing off a tree along the Choquequirao trek, Peru.
So, if you want to enjoy this site, go now, before it is too late. This trek is still a quiet and unique experience for those brave enough to hit the trail.
The route is 64 kilometres long with a starting altitude of 2,900m and peaking at 3,050m. While there is only one path, the route can be divided up in a variety of ways, but I have laid out the most common four-day route. It can easily be extended to five days, which will give you a full day to explore the ruins and enjoy your time in the mountains.
Fair warning: the Choquequirao trek difficulty level is high. There is no easy day where you can just stroll along.
While the Choquequirao trail is not at high altitude, there is a 1,500m descent into the canyon, followed by a 1,800m climb the next day. If you are up to this challenge, the trek is totally worth it. You will be able to hike unbothered the entire way and you can explore the ruins practically alone with the most stunning scenic backdrop.
This route begins and ends at Cachora, the town just outside the Apurimac Valley. You can get here from Cusco by bus or private taxi.
The trail begins in Cachora (2,900m) and leads down into the Apurimac Canyon. If you want to save yourself some time, you can catch a taxi to the Capuliyoc Mirador (2,915m), cutting 13 kilometres of road off the beginning of your trek. From this viewpoint, the path plunges down into the Canyon, and so begins the steep descent.
You can either stop overnight at Chiquisca (1,835m) before the bridge or carry on to the quieter Santa Rosa (2,905m) to make the next day easier for yourself.
This is the most difficult day, with 13 kilometres of steep ascent comprising of a series of switchbacks and false summits. Be warned the mosquitoes can be murder on this day and the shade fades quickly as you approach noon. Wear cool, long clothing to protect yourself. There are streams running down the mountain along the route for water refilling but to be safe have some way to purify the water (we love the Steripen Adventurer that works really well with a Nalgene reusable water bottle). Otherwise get water refills at Santa Rosa.
Eventually, you will reach the top and find a flag marking Marampata (2,913m), a small town that overlooks the canyon. Pitch camp here and enjoy playing football with the locals or drink a cold beer while admiring the view.
Trekking to Choquequirao: expect sublime mountain views.
If you still have energy after the climb, you can continue on the easy path to Choquequirao (2,950m) and spend the afternoon exploring the ruins. There is a campsite below the archaeological site, or you can return to Marampata.
Wake up early and wander in the ruins in the beautiful morning light. Then set off back the way you came, unfortunately, there is no alternative route for the return journey unless you want to hike all the way to Machu Picchu (see below). If you would like more time in the ruins, just spend the day there so you can tackle the return journey on fully rested legs.
Be careful not to slip on the steep descent back to Playa Rosalina (lightweight hiking poles such as these are a good investment for this and other South America treks). Continue past the river to the Chiquisca campsite so that you can save your energy for the uphill the next day.
More steep switchbacks will take you back up to the mirador; this is a great area for spotting Condors to keep your mind busy. At the top, have a cold drink before you continue on the easy downhill back to Cachora.
If you are looking for an even greater trek, the Choquequirao trail can be extended all the way to Machu Picchu. This route is 100km long with a much greater range in altitude, starting at 2,900m and reaching 4,668m at its peak. Acclimatisation in Cusco is recommended taking on this long trek.
The Choquequirao to Machu Picchu trek can't be done the other way around, so if you are interested in both sites, consider this option. On this route, you will see parts of the original Inca Road and see views very few people have experienced.
The first two days of the trek are the same as the above route, with the third day spent in the ruins. On day four, you continue forward on the path, instead of turning back. Hiking poles are recommended on this route to help on the daily steep climbs and descents.
After a relaxed day exploring the ruins, wake up early and continue over the ridge (3,300m) and down the steep descent to Rio Blanco (1,910m). Here is a great place to swim before tackling the tough ascent which will take you up to Maizal (3,000m).
Continue to climb up to Victoria Mines (4,200m) where you can see some of the originally paved Inca Road. After the Abra San Juan pass, you will start a long descent to Yanama (3,800m) where you will camp for the night.
The route through the Queliqua Valley is a day of light trekking until you reach the Yanama Pass (4,668m) – the highest point on this trek. A slight descent will leave you in Totora (4,100m) for the night.
A sheltered viewpoint at the Choquequirao ruins.
A long, winding descent will take you to Totora River, 2,000m below the village of the same name. Then a slight ascent leaves you at the village of La Playa (2,400m) which has a campsite for you to overnight in. This village has shops for you to restock in as well as hot springs to soothe your tired legs.
Here it is possible to take local transport to Lucma, otherwise, continue on your trek. At Lucma, the trail ascends for a three-hour climb to the newly discovered ruins of Llacpacta. From here you can catch a glimpse of Machu Picchu before descending to Hydroelectrica where the train departs for Aguas Calientes and takes about one hour.
Get up early and head to the Machu Picchu gates. Tours will take the bus up to the site, otherwise, you can hike up the hour and a half of stairs – if your legs are still strong enough. Spend the day exploring the ruins and take the train back to Cusco or Ollantaytambo (you will need to book the train and Machu Picchu ticket in advance).
Go it alone or take a guide for Choquequirao?
Unlike many hikes in Peru, guides on this trek are just to make logistics easier, not to help you navigate. This trek is doable alone, but it is one of the most difficult in Peru and you should make an honest assessment of your own fitness and experience before deciding whether to do the Choquequirao trek solo.
I did not do a tour as I had a Peruvian friend knew the trail and who helped us organise, so I got the best of both worlds. If this had not been the case, I would probably have joined a tour group.
Sunset along the Choquequirao trail in the Apurimac Valley.
Tours should include transport to and from Cusco as well as many necessities for the trek. They organise your food, camping equipment, mules, water, and emergency horses if you can't make it out on your own. It is obviously a more expensive option but means you only need to focus on the task of hiking the route.
Choquequirao tours are less frequent than the more famous treks but you should be able to find tours from around $500 or S/1600 (for just Choquequirao) and $1,500 or S/ 4900 (for the route to Machu Picchu). If you want to negotiate down the price, find a group before approaching the tour company.
Some companies offering the trek are Choquequirao Trek, Apus Peru and Choquequirao Trail. However, there are endless tour companies in Cusco and any that operate the Salkantay or Inca Trail treks will also offer the Choquequirao trek.
Tours operated out of Cachora are far cheaper because they don't include transport from Cusco and they don't get as much business. However, you also don't know what you are going to find there.
Ask in the hotels Inka Dream and CasaNostra about the trek but be sure to check what is included before you agree to leave. Ensure you will have sufficient equipment, meals, snacks and water for the trip.
Unlike the Inca Trail, it is possible to take on the Choquequirao trek without a guide. However, this is a very challenging route and you need to make an honest assessment of whether you are up to the task. Navigating the path is not difficult alone, and it means you will truly have the opportunity to experience the valley and the ruins in almost complete solitude.
Organising transport, food, mules, equipment and so on without any local assistance is definitely doable, it just requires more effort and preparation. If you decide to go it alone, you really need to ensure you are well prepared. I would not advise it unless you have experience in multi-day treks and are confident about your level of fitness.
To make the trek a little easier for yourself, consider hiring a muleteer so you do not have to carry your own bags. You will be expected to pay them about S/ 30 ($10) a day (including the time it will take them to return home) and share your meals with them.
Carrying your own gear will make this trek substantially more difficult, and if you doubt your strength (or sanity), you need to make the decision to hire a muleteer before you leave. Consider investing in a good, lightweight backpacking tent (read our review of the Big Agnes Copper Spur) and a water filter (read our review of the Steripen, a personal favourite) – both items that will make your life a lot easier on the trek.
Cachora is the only spot with mules for hire and once you hit the trail, you won't be able to change your mind.
More terraces, barely extracted from the encroaching jungle at Choquequirao, Peru.
If planning on doing the Choquequirao trek solo, you will need to pre-plan your meals. Though there are little stores and kitchens along the way, they are often closed or short on supplies as everything must be carried in with mules. To save yourself being caught short, bring most of the food you need.
Thinking of trekking to Machu Picchu? Check out hiking the Lares Trek for a unique alternative.
Also, consider a camping stove because you'll be craving warm food at the end of the long days.
When is the weather best for hiking to Choquequirao?
The Choquequirao trek is best done in the dry season, which is May to September, when the weather is clearest. Though this is the high season, the trail is quiet and uncrowded throughout the year – unlike bustling Machu Picchu.
In the rainy season, there are mudslides and the steep descents can be quite treacherous because of a lack of solid footing. Good hiking boots and a waterproof coat are essentials: check out our packing guide for Patagonia which has useful tips for equipment that will serve you well for hiking in Peru too.
A ceremonial plaza in the Choquequirao ruins.
What about entry tickets for Choquequirao?
The entrance ticket to the Choquequirao Archaeological Park costs S/ 55 ($16). They are bought upon entry to Choquequirao and do not need to be booked in advance.
Though the ticket is only valid for one day, they might let you back in the next morning if you arrived late the previous day.
There is a sign-in sheet at the ticket stand, so you will be able to see how many people have hiked in that day. When I arrived, there were just four names on that sheet.
Are you ready to pack your bags and head to the Apurimac Valley? If you have any further questions, please comment below or look at our Peru travel planning service if you want further inspiration and advice for your trip.
Remember, the proposed Choquequirao cable car will be carrying in busloads of tourists in no time, so no time is better than the present.
Don't miss these other beyond-the-beaten-trail adventures in South America for more inspiration.
Hi, I am looking to trek the Choquequirao alone with a muleteer. In this article, you told that we have to share our food with them. Do we also have to share our tent with them?


text:
Wealthy stock trader convicted in man's death in bunker fire
ROCKVILLE, Md. (AP) — A jury convicted a wealthy stock trader of second-degree murder and involuntary manslaughter late Wednesday in the fiery death of a man who was helping him secretly dig tunnels for an underground nuclear bunker beneath his Maryland home.
Jurors deliberated for about 12 hours before delivering their verdict in the case of 27-year-old stock trader Daniel Beckwitt.
Beckwitt had been charged with both offenses in the September 2017 death of 21-year-old Askia Khafra. He faces up to 30 years in prison at sentencing set for June 17.
Beckwitt slumped over a table, covered his face with his hands and wept after the verdict was read. He shook his head as court officers handcuffed him. Judge Margaret Schweitzer agreed to revoke Beckwitt's $100,000 bond.
The fire erupted as Khafra was digging tunnels under Beckwitt's Bethesda home, which was littered with piles of garbage.
A prosecutor accused Beckwitt of recklessly endangering Khafra's life and sacrificing safety for secrecy.
Defense attorney Robert Bonsib had told jurors the fire was an accident, not a crime. Beckwitt did not testify at his trial, which lasted nearly two weeks.
Bonsib said afterward that he would appeal. He said he believes the jurors may have been unfairly influenced by photographs they were shown of extreme hoarding conditions at the home in the upscale Maryland suburb of the nation's capital.
"The problem is the pictures of the hoarding condition did not reveal what the path out of the home looked like before the fire," Bonsib said. "I just told him there's a lot more fight in this case."
Montgomery County prosecutor Marybeth Ayres had said Beckwitt created the conditions that prevented Khafra from escaping the fire. Ayres said Beckwitt engaged in "extreme risk-taking behavior" and created a "death trap" in his family's home, with mounds of trash blocking Khafra's escape.
"This was a survivable fire, and we know that because the defendant survived," she said Tuesday during the trial's closing arguments.
Beckwitt went to elaborate lengths to keep the project a secret. Jurors heard that he tried to trick Khafra into thinking they were digging the tunnels in Virginia instead of Maryland by having him don "blackout glasses" before taking him on a long drive. They also were told Beckwitt also used internet "spoofing" to make it appear they were digging in Virginia.
Bonsib described his client as a "very strange young man" but urged jurors to look past his idiosyncratic personality.
"Being different, living in a different circumstance, is not a crime," he said.
Hours before the fire broke out in the basement, Khafra texted Beckwitt to warn him it smelled like smoke in the tunnels. Ayres said Beckwitt didn't respond for more than six hours before telling Khafra that there had been a "major electrical failure." Instead of getting Khafra out of the tunnels, Beckwitt told him that he "just switched it all over to another circuit," according to the prosecutor.
Bonsib said Beckwitt screamed for help from neighbors after the fire broke out and risked his own safety in a failed attempt to rescue his friend from the blaze.
The fire erupted in the basement as Khafra was digging tunnels lower down under Beckwitt's home in Bethesda, a Washington suburb. Firefighters found his naked, charred body in the basement, only a few steps from an exit.
Khafra met Beckwitt online. Beckwitt had invested money in a company Khafra was trying to launch as he helped Beckwitt dig the tunnels.
A hole in the concrete basement floor led to a shaft that dropped down 20 feet (6 meters) into tunnels that branched out roughly 200 feet (60 meters) in length. Investigators concluded the blaze was ignited by a defective electrical outlet in the basement.
Khafra worked in the tunnels for days at a time, eating and sleeping in there. They had lights, an air circulation system and a heater.
Bonsib said Khafra was a willing participant in the project. He argued there was no evidence, only speculation, to explain why Khafra died in the fire that day.
"This case is a mystery without an answer," Bonsib told jurors.
Accidents and disasters
Residential fires


text:
Kika Young Named One of the "People to Watch for 2018"
Kika Young, who has been a person to watch at Forest City Gear since she was a toddler, landed on the Rockford Register Star's "People to Watch in 2018" list.
The Register Star started the list in 2013, aiming to highlight Rock River Valley residents who are working to make the area better through "raising awareness, striving for greatness or simply giving back to their community."
This year's list was a typically diverse cross-section of professionals in Winnebago and Boone counties. Among those to watch: an aspiring Olympic boxing hopeful, the new head of Starlight Theatre at Rock Valley College, one of Rockford's few African-American attorneys, the leader of a nonprofit aimed at making men better fathers, and the creator of Stroll on State for the Rockford Area Convention and Visitors Bureau.
Kika, who is FCG's director of human resources and corporate administration, was included because of her desire to continue the work — to make FCG the best gear company in the world — started by her grandfather, FCG founder Stetler Young.
"It would be an honor to continue on this great legacy for my grandparents and parents," Kika told the Register Star. "In my mind, to be able to continue on that path … the third generation would be something I would really enjoy doing."
Kika listed her 2018 goal as wanting to consolidate FCG's vendor base.
The old-timers remember Kika making herself useful on the shop floor — painting, sweeping, doing routine maintenance — while she was a student at Keith Country Day School and in college. She joined the company full time in 2009 after finishing her MBA. It was an honor to see the Register Star recognize the potential FCG saw in her years ago.
If you'd like to see the entire piece, go to http://www.rrstar.com/news/20171229/people-to-watch-2018-kika-young-aspires-to-continue-family-legacy.


text:
Watson Valley () is a valley east of Mount Lewis in the Saint Johns Range of Victoria Land. The valley, which opens southward to Victoria Upper Glacier, is ice free except for a small glacier at the head wall. It was named by the Advisory Committee on Antarctic Names in 2007 after HM2 Donald E. Watson, U.S. Navy (Seabees), who at the time was Medical Assistant and member of the construction crew which built the original Little America V Station and the original Byrd Station in the 1955-57 pre-IGY period. He was the medical person on the oversnow Byrd Traverse to Byrd Station, 1956.

References

Valleys of Victoria Land

redpajama_set_name: RedPajamaWikipedia


Clocks And Culture 1300 To 1700 eBook you can also buy from amazon, please click here to purchase this book at amazon.
eBook Clocks And Culture 1300 To 1700 available at elcentrohistorico.com with Format PdF, ePub, Audiobook & Magazine. Please Create a FREE ACCOUNT to read or download Clocks And Culture 1300 To 1700 FOR FREE.
Finally I can download and read Clocks And Culture 1300 To 1700 Thank you!


Anywhere you have dampness and a lack of ventilation, mold is sure to follow. Mold in your home or business is not just unsightly; it can be harmful to those living or working in the space. Attic mold damage mold in basements commonly goes unnoticed for extended periods of time. If you need attic mold cleanup, contact us for attic mold removal services in the Snohomish area.
Even if you're not sure your ready for basement mold removal or mold clean help, feel free to contact us to discuss your basement mold cleanup or mold suspicions. We're happy to share our experience and answer mold removal questions for our neighbors in the Snohomish area.
PureDry is a mold removal company with IICRC certified technicians prepared to remediate and restore commercial mold damaged Snohomish properties back to their pre-loss condition. We understand the trauma of dealing with property and personal property or contents damage and make it our goal to return our customers' lives to normal as quickly as possible.
The Good News? Most Snohomish mold removal and cleanup jobs are covered 100% by your insurance, and we take care of the lengthy, confusing claim process. Additionally, YOU GET TO CHOOSE YOUR CONTRACTOR!
Snohomish was founded around 1858 by Emory C. Ferguson, E. F. Cady and others. It was originally known as Cadyville, and changed its name to Snohomish City in 1871. The name Snohomish comes from the name of the dominant local Native American tribe "sdoh-doh-hohbsh", whose meaning is widely disputed.