% Example LaTeX file for preparing camera-ready manuscripts for
% Nanotech conference proceedings.
% Nanotech.cls (LaTeX2e) style file required.
% Ver. Nov. 2003. % Please suggest changes to:
% bfr@nsti.org

\documentclass{iwce_13} % For LaTeX2e
\usepackage{epsfig}
\usepackage{amsmath}
\usepackage{amsfonts}


\title{\ \\[-36pt]Nanoelectronic Device Simulation Using\\[5pt]Extended H\"uckel Theory
(EHT) and NEGF}
\author{Zhiping Yu, Ximeng Guan, Ming Zhang, and Qiushi Ran
\\Institute of Microelectronics, Tsinghua University, Beijing 100084, China
\\yuzhip@tsinghua.edu.cn}

\begin{document}
\maketitle

\section*{\sc Abstract}

Nanoelectronic devices can be, in one way, characterized by the
large surface/volume ratio in addition to the central role of quantum
effects. This paper describes a computationally efficient way of obtaining
the bandstructure of the intrinsic device including the interface with
metal contacts using the extended H\"uckel theory (EHT). Carrier quantum
transport is then computed by NEGF (non-equilibrium Green's function).
GNR-FETs (graphene-nanoribbon FET) are simulated using this approach as
an application example.

\section{\sc Introduction}

Recent experimental work on GNR-FETs \cite{hdai_prl_may_08}
reconfirms the critical role of Schottky contact between the
semiconductive channel and metal source/drain (S/D) regions in
the transistor operation for carbon crystal, e.g. CNT (carbon nanotubes)
and GNR, based field-effect devices. It is actually a common fact
that all nano-electronic devices are characterized by a large
surface-volume ratio and the correct modeling of the material
interface becomes an essential part of device simulation.

%Due to the size confinement, the bandstructure in the intrinsic
%region of nano-devices departs from its bulk counterpart
%appreciably. Various semi-empirical methods have been applied to the
%{\em in situ} calculation of bandstructure, e.g., $sp^3d^5s^*$
%tight-binding (TB) \cite{klimeck_prb_04} approach is known to be
%computationally efficient and with good accuracy. The drawback is,
%however, TB relies on empirical parameters which blur the underlying
%atomic composition.

%Another semi-empirical approach, EHT for extended H\"uckel theory
%originally developed as a quantum mechanical method in calculating
%electronic structure for molecules, uses explicit atomic-like
%orbital (AO) basis functions to construct all matrix elements of the
%Hamiltonian. It is thus more physics transparent and in particular
%is suitable for hetero-atomic systems such as material interface and
%surface bond termination.

%There are recently extensive efforts in extending the EHT for the
%bandstructure calculation in 1D (e.g., CNT), 2D (e.g., surface of a
%solid), and 3D crystals, see for example \cite{Kienle_06}.

%This paper describes the basic principles for EHT with the emphasis on the
%calculation of electronic structure for realistic device structures
%at non-equilibrium state, i.e., under biasing condition. To obtain
%$I$-$V$ characteristics, NEGF method is applied to simulate the
%carrier quantum transport in nano-devices. Simulation of a GNR-FET
%is finally given as an application example.

\section{Principles of EHT Method}

EHT can be characterized as one of LCAO (linear combination of
atomic orbitals) methods in solving\\ Schr\"odinger equation for a system to
obtain its electronic structure: either energy levels for molecules
or bandstructure of $E(\mathbf{k})$ for crystals. Different from orthogonal
TB method, the basis functions in EHT are directly tied to the atomic orbitals
with a more localized treatment, the Slater-type orbitals (STO) \cite{Slater_30}
as follows,
\begin{equation}
 \phi_{nlm}(\mathbf{r})=\sigma r^{n-1}e^{-\zeta r}Y_{lm}(\theta,\varphi)
(=\mathrm{STO})
\end{equation}
where spherical coordinates $r,\theta,\varphi$ have been used for
$\mathbf{r}$, and $n\ (=1,2,\ldots),l\ (=0,1,\ldots,n-1), m\ (=-l,
-l+1, \ldots, 0, 1, \ldots, l-1, l)$ are principal, angular
momentum, and magnetic quantum numbers, respectively, and
$Y_{lm}(\theta,\varphi)$ is the real-valued spherical harmonics.
$\zeta$ is called the orbital exponent and is related to the
effective nuclear charge (the nuclear charge partially shielded by
electrons in the core of an atom). $\sigma$ is the normalization
constant and is obtained from
\begin{equation}
 \sigma^2\int_0^\infty \left(r^{n-1}e^{-\zeta
 r}\right)^2r^2dr=1\Longrightarrow
 \sigma=\frac{(2\zeta)^{n+1/2}}{\sqrt{(2n)!}}
\end{equation}
With the above AOs, we consider first molecules or a primitive cell
in crystals. The complete set of basis functions are
\begin{eqnarray}
 &&[\phi_1(\mathbf{r})\ \phi_2(\mathbf{r}) \cdots
 \phi_N(\mathbf{r})]^T= \\
&&[\phi_{11}\ \phi_{21} \cdots
 \phi_{N_\mathrm{ao}1}\ \phi_{12} \cdots \phi_{N_\mathrm{ao}2}
 \cdots \phi_{1N_\mathrm{at}} \cdots
 \phi_{N_\mathrm{ao}N_\mathrm{at}}]^T \nonumber
\end{eqnarray}
where we have omitted the explict function argument $\mathbf{r}$ on
the right-hand side, $N=N_{\rm ao}\times N_{\rm at}$ with $N_{\rm
ao}$ the number of AOs associated with each atom (assumed the same
for all atoms) and $N_{\rm at}$ that of atoms (in the molecule or
cell).
\begin{equation}
 \phi_{ij}(\mathbf{r})=\mathrm{STO}_i(\mathbf{r}-\mathbf{d}_j),\ i=1,\ldots,
N_\mathrm{ao}, j=1,\ldots,N_\mathrm{at}
\end{equation}
where $\mathbf{d}_j$ is the displacement of atom $j$ relative to
atom 1 ($\mathbf{d}_1=\mathbf{0}$) in the cell.

%For crystals, because of the translational symmetry we need only consider atoms
%in a primitive cell and take into consideration Born-von Karman
%periodic boundary condition (BC) for the crystal wavefunction. The
%AOs are replaced by Bloch sums as follows

%\begin{equation}
% \phi_{ij}(\mathbf{r};\mathbf{k})=\frac1{N_\mathrm{cell}}\sum_{\mathbf{R}}
%e^{i\mathbf{k}\cdot
%\mathbf{R}}\phi_{ij}(\mathbf{r}-\mathbf{d}_j-\mathbf{R})
%\end{equation}
%where $N_\mathrm{cell}=N_1\times N_2\times N_3$ with $N_1, N_2, N_3$
%numbers of cells in each dimension for 3D block as the repeated
%region for the periodic BC of crystal wavefunctions.  $\mathbf{R}$
%is the vector for the Bravais lattice and sum runs over all
%primitive cells in the repeated region. Note that the Block sum
%$\phi_{ij}(\mathbf{r})$ have wavevector $\mathbf{k}$ dependence.

%Now we can expand wavefunctions as the linear combination of these basis
%functions:
%\begin{equation}
% \Psi(\mathbf{r})
% =\mathbf{c}\cdot\left(\begin{array}{c}
%  \phi_1(\mathbf{r}) \\
%  \phi_2(\mathbf{r}) \\
%  \cdots \\
%  \phi_N(\mathbf{r})
% \end{array}
% \right)
%\end{equation}
%where $\mathbf{c}$ is the vector of coefficients,
%\begin{equation}
% \mathbf{c}=(c_{1}\ c_{2}\ \cdots\ c_{N})^T\in \mathbb{C}^{N\times 1}
%\end{equation}
%Now, we need to solve Schr\"odinger equation (and in the following
%we leave off the explicit $\mathbf{k}$
%dependence)\\$H\Psi_n=\varepsilon_n\Psi_n$.

%In the representation of non-orthogonal basis functions, $\phi_1,\ldots,\phi_N$,the Schr\"odinger equation becomes
%\begin{equation}
% H\mathbf{c}_n=\varepsilon_n S\mathbf{c}_n
%\end{equation}
%where $n$ is the energy level/band index and both Hamiltonian matrix
%$H$ and overlap matrix $S$ are of dimension
%\begin{equation}
% H, S\in \mathbb{C}^{N\times N}
%\end{equation}
%The overlap matrix elements are simply
%\begin{equation}
% S_{ij}=\int \phi_i^*(\mathbf{r})\phi_j(\mathbf{r}) d\mathbf{r},\ i,j=1,\ldots,N
%\end{equation}
%And it can be noticed that $S_{ii}=1$ in the case of molecules.

The elements of Hamiltonian $H$ are calculated in principle from
\begin{equation}
 H_{ij}=\int \phi_i^*H\phi_j d\mathbf{r}
\end{equation}
where $H$ is the Hamiltonian for the entire system. In practice, EHT uses
simple formulas to evaluate Hamiltonian matrix elements in the following way:
\begin{eqnarray}
 H_{ii}&=&\varepsilon_i+U_{i,\mathrm{shift}}-qV_i \label{eq:hii}\\
 H_{ij}&=&\frac1{2}S_{ij}[K_\mathrm{EHT}(\varepsilon_i+\varepsilon_j)+
U_{i,\mathrm{shift}}+U_{j,\mathrm{shift}} \nonumber \\
 && \qquad -q(V_i+V_j)],\quad i\ne j \label{eq:hij}
\end{eqnarray}
The above expression for $H_{ij}$ is named after Wolfsberg-Helmholtz, where
$K_\mathrm{EHT}$ is a fitting parameter of typical value of 1.75 for
molecules, 2.3 for solids, and 2.8 for 2D graphene sheets \cite{Kienle_06}.

Special discussion is needed for formula of $H_{ii}$ (\ref{eq:hii}).
$\varepsilon_i$ is the on-site energy (also called the valence state
ionization potential, or VSIP) and is an empirical parameter listed,
e.g., in \cite{Herman_04} for all elements in the Periodic Table.
Since EHT-related parameters may be obtained from different energy
reference setting (e.g., the Fermi-level for transitional metals is
set to -10\,eV and the top of valence band for semiconductors is set
to -13\,eV), it is important to make sure that the vacuum level is
set to the same level for all the atoms in the system, thus arising
the parameter $U_{i,\mathrm{shift}}$ in, where all $U$s
corresponding to a single atom have the same value. In addition, the
potential at each atom is a variable determined from solving
Poisson's equation, so it may need to be accounted for in rigorous
calculations. A slightly different recipe for obtaining $H_{ij}$ can
be found in \cite{Kienle_06}.

Thus completed the description of the EHT method.

\section{NEGF Simulation of Quantum Transport}

NEGF is used in simulation of carrier quantum transport in the
system with the open boundary conditions. The surface Green
functions, which characterize the contacts to the system are
obtained by an iteration scheme. The Green's function in the device
region is calculated by
\begin{equation}
 G^R=[(E+\mathrm{i}\eta)S-H(U_H)-\Sigma]^{-1}\ \mathrm{with}\ \eta\to 0^+
\end{equation}
where $G^R$ is the retarded Green's function, $H$ is the Hamiltonian
of the system with closed BCs, whose eigen-energy can be labeled $\epsilon_n\
(n=1,\ldots)$; $E$ is a given energy which is continuous but not equal
to $\epsilon_n$; $U_H$ is the external potential energy (called
Hartree potential) and will affect $H$ through (\ref{eq:hii}-\ref{eq:hij})
(hence $H$ as function of $U_H$in the above expression),
and $\Sigma$ is the self-energy characterizing
the semi-infinite contact regions. $U_H$ depends on the carrier
concentration profile and is solved from the Poisson's equation. So both
Schr\"odinger and Poisson's equations have to be solved self-consistenly.

\section{Simulation Example}

As an example of applying EHT-NEGF method for device simulation, we
have investigated how the width of an armchair GNR (aGNR) affects
the bandgap in the electronic structure. In our simulation, it is
found that the bond length is weakly dependent on the number of
carbons across the ribbon width ($N_a$ in Fig.~2, and the $C$-$C$
bonds in the channel
\begin{figure}[ht]
  \centering
  \includegraphics[width=0.4\textwidth]{gnr_fet_structure.eps}
  \caption{Schematics of GNR-FET. (a) The armchair GNR structure
  with $N_a = 7$ and $M_a = 29$ is chosen for the simulated region;
  (b) The DG structure of the simulated device. The Fermi level is set to
 -12.3\,eV at drain and source while that for GNR channel -13.0\,eV
 [6], and $t_\mathrm{ox}=0.5$\,nm.}
 \label{fig:gnr_fet}
\end{figure}
direction at two edges of the ribbon are shortened by 3.4\% to
reflect the effects of the termination of graphene edge. The results
from our simulation (Fig.~5) are similar to those from the TB and
%\begin{figure}[ht]
%  \centering
%  \includegraphics[width=0.4\textwidth]{icsict_08_invited_figs/bandgap_variation.eps}
%  \caption{Bandgap variation. The bandgap $\Delta_a$ vs. ribbon width
%  $W_a$ with the number of atoms $N_a$ (defined in Fig.~\ref{fig:gnr_fet})
%  as parameter. TB and EHT results share a common trend in the bandgap,
%  but differ slightly in value. See [9] for the TB parameter. $a_{cc}$ equals
%  1.424, 1.422, and 1.423\AA, respectively, for $N_a = 3p, 3p+1$, and $3p+2$
%  with $p$ positive integer.}

%  \label{fig:bandgap_var}
%\end{figure}
first-principles calculations. Unlike TB, EHT does not need
structure-dependent fitting parameters, especially for systems whose
physical properties are strongly affected by the geometry, and the
computational cost is acceptable. Fig.~2 shows how bandgaps at
different $N_a$ vary with the bond length at edges. Since the
electronic structure of armchair GNR has a quantized
width-dependence, GNR with $N_a=3p+1$ (where $p$ is a positive
integer) will provide good performance. Fig.~2 shows the double-gate
(DG) GNR-FET in our simulation with $N_a=7$ and $M_a=29$. A
comparison of transport characteristics simulated using the
EHT-based and TB-based NEGF is shown in Fig.~3.
%\begin{figure}[ht]
%  \centering
%  \includegraphics[width=0.4\textwidth]{icsict_08_invited_figs/eht_tb.eps}
%  \caption{The drain current versus gate voltage under different bias
%  voltages derived by TB and EHT methods. The TB results are shown in
%  solid line and dashed line. The EHT results are plotted under different
%  marks. A large subthreshold current is observed in EHT result when
%  $a_{cc}=1.44$\,\AA\ while the TB curve maintain similar to the one of
%  $a_{cc}=1.42$\,\AA, which exhibits a larger subthreshold slope.}

%  \label{fig:eht_tb}
%\end{figure}

Bond length relaxation of GNR is taken into account in the EHT
calculation, and has a significant influence on the transport.
Fig.~4 shows the output characteristics simulated
%\begin{figure}[ht]
%  \centering
%  \includegraphics[width=0.4\textwidth]{icsict_08_invited_figs/simulated_output_char.eps}
%  \caption{Simulated output characteristics} \label{fig:ids_vds}
%\end{figure}
with independent gate operation. A stronger short-channel effect in EHT
simulation is observed, compared to TB, as well as an exceedingly high
drain current in the saturation region with $a_{cc}=1.44$\AA.

%In Fig.~\ref{fig:transmission} the transmission is shown. The transmission curve
%\begin{figure}[ht]
%  \centering
%  \includegraphics[width=0.4\textwidth]{icsict_08_invited_figs/transmission.eps}
%  \caption{The transmission given by TB and EHT of the two bond lengths
%  under $V_{DS}=0.5$\,V, $V_{g,\mathrm{top}}=0.6$\,V, $V_{g,\mathrm{bot}}=0$\,V
% (translated to the same $E_{FS}$). The
%  curves of three cases are plotted in different type of lines. $E_{FS}$ and
%  $E_{FD}$, marked with upper and lower triangles, are the Fermi levels of
%  source and drain ends.  Transmission of EHT with $a_{cc}=1.44$\,\AA\ (plotted
%  in dotted line) is relatively downward shifted, leaving a larger
%  integration within the interval between $E_{FS}$ and $E_{FD}$ than the one of
%  other two cases.}
%  \label{fig:transmission}
%\end{figure}
%is downward shifted relatively in the case of $a_{cc}=1.44$\AA. The
%increase of bond length reduces the coupling of the orbitals, thus
%narrowing the bandgap, which induces more transmission between the
%Fermi levels of the drain and source ($E_{FS}$ and $E_{FD}$). This
%can only be reflected by the EHT-based simulation, which provides us
%the way to predict the performance of GNR devices under structural
%deformation.

\section{Conclusions}

Even though both TB and EHT are semi-empirical, i.e., dpending on the
fitting parameters to some extent, in nature, the latter has more
direct link to the atomic orbitals of the constitutive atoms. This feature
of EHT can help identify the charge transfer for the chemical bonds at the
material interface and has significant applications such as determining the
dipole at the interface.

\section{Acknowledgement}

The authors are indebted to the grants from NSF (\#90307016) and from
Ministry of Science and Technology (\#2006CB302700), both in China.


\begin{thebibliography}{99}

 \bibitem{hdai_prl_may_08} X. Wang, Y. Ouyang, X. Li,
 H. Wang, J. Guo, and H. Dai, ``Room-Temperature
 All-Semiconducting Sub-10-nm Graphene Nanoribbon Field-Effect
 Transistors,'' {\em Physical Review Letters}, {\bf 100}, 206803
 (May 2008).

% \bibitem{klimeck_prb_04} T. B. Boykin, G. Klimeck, and F.
% Oyafuso, ``Valence band effective-mass expressions in the $sp^3d^5s^*$
% empirical tight-binding model applied to a Si and Ge parametrization,''
% {\em Physical Review B}, {\bf 69}, 115201 (2004).

% \bibitem{hoffmann_jcp_63} Ronald Hoffmann, ``An extended H\"uckel
% theory. I. Hydrocarbons,'' {\em The Journal of Chemical Physics},
% vol. 39, no. 6, p. 1397, Sept. 1963.

 \bibitem{Kienle_06} D. Kienle, K. H. Bevan, G.-C. Liang, and L.
 Siddiqui, J. I. Cerda, and A. W. Ghosh, ``Extended H\"uckel theory for
 band structure, chemistry, and transport.  II. Silicon,'' {\em Journal
 of Applied Physics}, {\bf 100}, 043715 (2006).

 \bibitem{Slater_30} J.C. Slater, {\em Atomic Shielding Constants, Phys. Rev.}
 vol. 36, p. 57 (1930).

 \bibitem{Herman_04} Aleksander Herman, ``Empirically adjusted and
 consistent set of EHT valence orbital parameters for all elements of
 the periodic table,'' {\em Modeling Simul. Mater. Sci. Eng.} {\bf 12}
 21-32 (2004).

% \bibitem{Cerda_00} J. Cerda and F. Soria, ``Accurate and transferable
% extended H\"uckel-type tight-binding parameters," {\em Phys. Rev. B}
% $\mathbf{61}$, 7965, March 2000.

% \bibitem{Sancho_85} M. Sancho, {\em et al., J. Phys.F: Met.Phys.}
% vol. 15, p. 851, 1985.

% \bibitem{NEGF} D. Nikonov, \verb@www.nanohub.org/resources/1983@, 2006.

% \bibitem{Song_06} Y. W. Song, {\em et al., Phys. Rev. Lett.,} vol. 97,
% 216803-1, 2006.

% \bibitem{Porezag_95} D. Porezag, Th. Frauenheim, Th. K\"ohler,
% ``Construction of tight-binding-like potentials on the basis of
% density-functional theory: Application to carbon,'' {\em Phys. Rev. B},
% vol. 51, p. 12947, 1995.

% \bibitem{Guan_07} X. Guan, M. Zhang, Q. Liu, and Z.
% Yu, ``Simulation Investigation of Double-Gate CNR-MOSFETs with a Fully
% Self-Consistent NEGF and TB Method,'' {\em Int'l Electron Devices Meeting}
% ({\em IEDM}), p. 761, Washington DC, Dec. 2007.

\end{thebibliography}

\end{document}