Point defect energetics in silicon using the LDA+U method
% insert abstract here We present the first principles results of point defect energetics in silicon calculated using the LDA+U method: a Hubbard type on-site interaction added to the local density approximation (LDA). The on-site Coulomb and exchange parameters were tuned to match the experimental band gap in Si. The relaxed configuration was obtained using the LDA; LDA+U was used only to calculate the energies. Our values of point defect energetics are in very good agreement with both recent experimental results and quantum Monte Carlo (QMC) calculations.
ICAM2003lanl.tex
%% ****** Start of file template.aps ****** %
%%
%%
%%   This file is part of the APS files in the REVTeX 4 distribution.
%%   Version 4.0 of REVTeX, August 2001
%%
%%
%%   Copyright (c) 2001 The American Physical Society.
%%
%%   See the REVTeX 4 README file for restrictions and more information.
%%
%
% This is a template for producing manuscripts for use with REVTEX 4.0
% Copy this file to another name and then work on that file.
% That way, you always have this original template file to use.
%
% Group addresses by affiliation; use superscriptaddress for long
% author lists, or if there are many overlapping affiliations.
% For Phys. Rev. appearance, change preprint to twocolumn.
% Choose pra, prb, prc, prd, pre, prl, prstab, or rmp for journal
%  Add 'draft' option to mark overfull boxes with black boxes
%  Add 'showpacs' option to make PACS codes appear
%  Add 'showkeys' option to make keywords appear
\documentclass[aps,prb,twocolumn,groupedaddress,showkeys,amsmath,amssymb]{revtex4}
%\documentclass[aps,prl,preprint,superscriptaddress]{revtex4}
%\documentclass[aps,prl,twocolumn,groupedaddress]{revtex4}

% You should use BibTeX and apsrev.bst for references
% Choosing a journal automatically selects the correct APS
% BibTeX style file (bst file), so only uncomment the line
% below if necessary.
%\bibliographystyle{apsrev}

\usepackage{longtable}
\usepackage{graphicx}
\usepackage{dcolumn}
\usepackage{bm}

\begin{document}

% Use the \preprint command to place your local institutional report
% number in the upper righthand corner of the title page in preprint mode.
% Multiple \preprint commands are allowed.
% Use the 'preprintnumbers' class option to override journal defaults
% to display numbers if necessary
\preprint{Submitted to the Journal of Applied Physics}

%Title of paper
\title{Point defect energetics in silicon using the LDA+U method}

% repeat the \author .. \affiliation  etc. as needed
% \email, \thanks, \homepage, \altaffiliation all apply to the current
% author. Explanatory text should go in the []'s, actual e-mail
% address or url should go in the {}'s for \email and \homepage.
% Please use the appropriate macro foreach each type of information

% \affiliation command applies to all authors since the last
% \affiliation command. The \affiliation command should follow the
% other information
% \affiliation can be followed by \email, \homepage, \thanks as well.
\author{Panchapakesan Ramanarayanan}
\email[]{panchram@stanford.edu}
\affiliation{Department of Mechanical Engineering, Stanford University,
  California 94305, USA}
\author{Renat F. Sabirianov}
\email[]{rsabirianov@mail.unomaha.edu}
\affiliation{Department of Physics, University of Nebraska, Omaha,
Nebraska 68182, USA}
\author{Kyeongjae Cho}
\email[]{kjcho@stanford.edu}
%\homepage[]{Your web page}
%\thanks{}
%\altaffiliation{}
\affiliation{Department of Mechanical Engineering, Stanford University,
  Stanford, California 94305, USA}
%Collaboration name if desired (requires use of superscriptaddress
%option in \documentclass). \noaffiliation is required (may also be
%used with the \author command).
%\collaboration can be followed by \email, \homepage, \thanks as well.
%\collaboration{}
%\noaffiliation

%\date{\today}

\nopagebreak

\begin{abstract}
% insert abstract here
We present the first principles results of point defect
energetics in silicon calculated using the LDA+U method: a
Hubbard type on-site interaction added to the local density approximation
(LDA).  The on-site Coulomb and exchange parameters were tuned to match
the experimental band gap in Si.  The relaxed configuration was
obtained using the LDA; LDA+U was used only to calculate the energies.
Our values of point defect energetics are in very good agreement with both
recent experimental results and quantum Monte Carlo (QMC) calculations.
\end{abstract}

% insert suggested PACS numbers in braces on next line
% \pacs{61.72.Ji, 61.72.Ss, 66.30.Hs}
% insert suggested keywords - APS authors don't need to do this
\keywords{Silicon, self-diffusion, point defects, local density approximation, 
Hubbard, LDA+U}

%\maketitle must follow title, authors, abstract, \pacs, and \keywords
\maketitle

% body of paper here - Use proper section commands
% References should be done using the \cite, \ref, and \label commands
%\section{}
% Put \label in argument of \section for cross-referencing
%\section{\label{}}
%\subsection{}
%\subsubsection{}

\section{\label{sec:intro}Introduction}
\input{introlanl}

\section{\label{sec:method}Method}
\input{methodlanl}

\section{\label{sec:results}Results and Discussion}
\input{resultslanl}

\section{\label{sec:summary}Summary}
\input{summarylanl}

% If in two-column mode, this environment will change to single-column
% format so that long equations can be displayed. Use
% sparingly.
%\begin{widetext}
% put long equation here
%\end{widetext}

% figures should be put into the text as floats.
% Use the graphics or graphicx packages (distributed with LaTeX2e)
% and the \includegraphics macro defined in those packages.
% See the LaTeX Graphics Companion by Michel Goosens, Sebastian Rahtz,
% and Frank Mittelbach for instance.
%
% Here is an example of the general form of a figure:
% Fill in the caption in the braces of the \caption{} command. Put the label
% that you will use with \ref{} command in the braces of the \label{} command.
% Use the figure* environment if the figure should span across the
% entire page. There is no need to do explicit centering.

% \begin{figure}
% \includegraphics{}%
% \caption{\label{}}
% \end{figure}

% Surround figure environment with turnpage environment for landscape
% figure
% \begin{turnpage}
% \begin{figure}
% \includegraphics{}%
% \caption{\label{}}
% \end{figure}
% \end{turnpage}

% tables should appear as floats within the text
%
% Here is an example of the general form of a table:
% Fill in the caption in the braces of the \caption{} command. Put the label
% that you will use with \ref{} command in the braces of the \label{} command.
% Insert the column specifiers (l, r, c, d, etc.) in the empty braces of the
% \begin{tabular}{} command.
% The ruledtabular enviroment adds doubled rules to table and sets a
% reasonable default table settings.
% Use the table* environment to get a full-width table in two-column
% Add \usepackage{longtable} and the longtable (or longtable*}
% environment for nicely formatted long tables. Or use the the [H]
% placement option to break a long table (with less control than 
% in longtable).
% \begin{table}%[H] add [H] placement to break table across pages
% \caption{\label{}}
% \begin{ruledtabular}
% \begin{tabular}{}
% Lines of table here ending with \\
% \end{tabular}
% \end{ruledtabular}
% \end{table}

% Surround table environment with turnpage environment for landscape
% table
% \begin{turnpage}
% \begin{table}
% \caption{\label{}}
% \begin{ruledtabular}
% \begin{tabular}{}
% \end{tabular}
% \end{ruledtabular}
% \end{table}
% \end{turnpage}

% Specify following sections are appendices. Use \appendix* if there
% only one appendix.
%\appendix
%\section{}

% If you have acknowledgments, this puts in the proper section head.
\begin{acknowledgments}
This work was funded by the United States Department of Energy, Basic
Energy Sciences Grant No. DE-FG03-99ER45788.  The computations were
performed on the NPACI IBM p690, NERSC IBM SP RS/6000 (as part CMSN),
and on the PC
clusters of the Multiscale Simulation Laboratory, Stanford: {\sc multipod}
and {\sc xepod}. Grants from NPACI and NERSC are sincerely appreciated.
\end{acknowledgments}

% Create the reference section using BibTeX:
\bibliography{ICAM2003lanl}

\end{document}
%
% ****** End of file template.aps ******

introlanl.tex
Modeling techniques empower the investigator with the capability to
explore into the details of natural phenomena with great dexterity, aptly
complementing their experimental counterparts.  Hierarchical multiscale
modeling (HMM) is one such technique.  We have recently investigated self
diffusion in silicon-germanium alloys in considerable detail by developing
a database of first principles energetics and using the database to
perform kinetic Monte Carlo simulations: a typical HMM scheme
\cite{rama03a}.  The exponential growth in computational 
processing power makes such schemes more viable; the explosive growth in
nanotechnology provides technologically relevant platforms that are
amenable to such schemes. 

Understandably, the accuracy of the results such as those presented in
Ref.~\onlinecite{rama03a} hinges to a great extent on the correctness of 
the energetics database.  In that study, we used the popular local density
approximation (LDA) to develop the energetics database.  However, the
(approximately 1 eV) discrepancy \cite{souz98a} between the theoretical activation energy for Si
self diffusion computed using the LDA (or the generalized gradient
approximation (GGA)) and the experimental values was an important reason
that precluded us from making direct comparisons of our results with
experimental studies of self diffusion in SiGe alloys such as those in
Refs.~\onlinecite{zang01a, stro02a}.

Our present work is directed essentially at finding a means of developing
a reliable point defect energetics database (like in
Ref.~\onlinecite{rama03a}) using an as parameter free a technique as
possible.  (While considerable progress has been made on techniques like
the quantum Monte Carlo (QMC) method, it is beyond the current computational
capacity to use it to develop a detailed energetics database as required
for our HMM scheme.)  We focus on the activation energy for self-diffusion
in pure Si. For self diffusion in Si, we consider the vacancy mechanism,
the interstitialcy mechanism, and the concerted exchange mechanism. These
three have been shown to be the most likely diffusion mechanisms in
silicon \cite{ural98a,ural99a}. For the interstitialcy mechanism, we
consider the hexagonal to split $\langle 110 \rangle$ mechanism
as this has shown to be the most likely mechanism \cite{need99a}. 
 
This article is organized as follows: In Sec. \ref{sec:method}, we provide the details
of the LDA and the LDA+U calculations.  In Sec. \ref{sec:results}, we present the results
of our calculations.  We provide an intuitively appealing explanation for
the underperformance of the LDA.  We then provide a brief overview of the
LDA+U technique and discuss the results obtained using this method for the
different diffusion mechanisms. We summarize our article in Sec. \ref{sec:summary}.

methodlanl.tex
Self-consistent electronic structure calculations were performed using
the plane-wave pseudopotential code VASP \cite{VASP1,VASP2,VASP3,VASP4}
with the projector augmented-wave (PAW) potentials \cite{VASPAW1,VASPAW2} at a kinetic energy cutoff
of 20 Rydberg.  A 64 atom supercell 
was used.  Electronic minimization was carried out to a tolerance of 
$2.7 \times 10^{-5}$ eV
and structures were relaxed until the maximum force on any atom
was less than 0.015 eV/{\r{A}}.  Saddle point configuration for the hexagonal
to split $\langle 110 \rangle$ interstitialcy mechanism was determined using the nudged
elastic band method \cite{NEB}.  Structural relaxation was performed in two stages:
initially with a $2^3$ Monkhorst-Pack \cite{monk76a} {\bf k}-point sampling followed
by a $6^3$ Monkhorst-Pack \cite{monk76a} {\bf k}-point sampling.  Energy calculations were done
using tetrahedron method with a $6 \times 6 \times 6$ division of Brillouin zone.
The lattice constant of systems
containing point defects (vacancy or interstitial) was determined by
fitting the total energy versus the supercell volume to Murnaghan's
equation of state\footnote{The Murnaghan's equation of state \cite{murn44a} is a
relation between the energy $E$ and the volume $V$
of the supercell and has the following form:
$E = E_0 + [(V-V_0)/b] - \{[(V^{1-a} - V_0^{1-a})V_0^a]/[(1-a)b]\}$.
$E_0$, $V_0$, $a$, and $b$ are the fitting parameters.}.

The plane-wave psedopotential code VASP \cite{VASP1,VASP2,VASP3,VASP4} also has
the option to perform calculations by the LDA+U method. (A brief overview of the
method itself will be presented in Sec. \ref{sec:results}. Here we merely give the details of the
parameters used in our calculations.) We have used the rotationally invariant LDA+U scheme 
according to Liechtenstein et al. \cite{liec95a}. The on-site interactions were added
for the {\em p}-orbitals. The effective on site Coulomb interaction parameter ($U$)
was set to 0 eV and the effective on site exchange interaction parameter ($J$) was 
set to 4 eV. (We provide justification for this choice in Sec. \ref{sec:results}.) The LDA
was used to perform structural relaxation; the LDA+U was used only to
perform energetics calculation on the structure obtained using the LDA. We
adopted this approach\footnote{We would like to point out that 
Leung {\em et. al.} \cite{leun99a} take a 
similar approach: they calculate their structure using LDA and use QMC
to compute the energetics.}
because the LDA gives structutal properties
(eg. lattice constant, see Table \ref{tab:resultstab}) that are closer to experimental values
than the LDA+U.

resultablelanl.tex
\newcommand{\rb}[1]{\raisebox{1.5ex}[0pt]{#1}}
\begin{table}
\caption{\label{tab:resultstab}
  Theoretical results and experimental values for the following
  properties[symbol](units) of silicon: 
  lattice constant[$L_0$]({\r{A}}),
  vacancy formation volume$^{\dagger}$[$V^f_V$]({\r{A}}$^3$),
  vacancy formation energy[$E^f_V$](eV),
  vacancy migration energy[$E^m_V$](eV),
  interstitial formation volume$^{\dagger\dagger}$[$V^f_I$]({\r{A}}$^3$),
  hexagonal interstitial formation energy[$E^f_{HI}$](eV),
  split $\langle 110 \rangle$ interstitial formation energy[$E^f_{SI}$](eV),
  interstitial migration energy[$E^m_I$](eV),
  concerted exchange migration energy[$E^m_{CX}$](eV). Superscripts indicate
  reference numbers. Quantum Monte Carlo results are from Ref. \onlinecite{leun99a}.}
\begin{ruledtabular}
\begin{tabular}{lcccc}
Prop. & LDA & Expt. & QMC & LDA+U\\
\hline
$L_0$      & 5.4   & 5.43 \cite{kitt96a}                   & ---  & 5.27 \\
$V^f_V$    & 4.41  & ---                                   & ---  & ---  \\
$E^f_V$    & 3.49  &                                       & ---  & 4.59 \\
$E^m_V$    & 0.03  & \rb{4.86 \cite{ural99b}$^{\ddagger}$} & ---  & 0.42 \\
$V^f_I$    & $-$2.77 & ---                                   & ---  & ---  \\
$E^f_{HI}$ & 3.37  &                                       & 4.82 & 4.61 \\
$E^f_{SI}$ & 3.34  & 4.68 \cite{ural99b}$^{\ddagger}$      & 4.96 & 4.70 \\
$E^m_I$    & 0.18  &                                       & ---  & 0.44 \\
$E^m_{CX}$ & 4.56  &                                       & 5.78 & 5.82 \\
\end{tabular}
$^{\dagger}$$\;V^f_V$ = Relaxation volume $+$ Atomic volume\\
$^{\dagger\dagger}$$\;V^f_I$ = Relaxation volume $-$ Atomic volume\\
$^{\ddagger}$ These values are the experimental estimates of the activation energy
i.e., the sum of formation and migration energies.
\end{ruledtabular}
\end{table}

resultslanl.tex
Table \ref{tab:resultstab} summarizes the results of both the LDA and LDA+U calculations.
It also contains experimental values and values from quantum Monte Carlo
simulations where available.

There is consensus on the values obtained within the LDA by different
theoretical groups \cite{prob03a,pusk98a,need99a}. It can be seen that the
LDA gives activation energy that is systematically lower than the experimental
values. 

\subsection{\label{sec:ldaresult} Local density approximation}

One of the main reasons for the underperformance of the LDA is that it
tends to overbind.  Because of this, the dangling bonds that are created
due to the defects, tend to bind with each other and with other orbitals.
Because of this overbinding, the energy of the system with the defect is
lower than it should be in reality.  Thus the difference between the
energy of the pure crystal and the one with the defects (which is in fact
the formation or the migration energy, depending on the situation) is
lower than it should be in reality.  The reason for the LDA's overbinding
can be traced to the spherical nature of the exchange hole. This leads to
a larger negative value of the exchange energy and hence lowers the energy
of the system much more than it should be in reality \cite{koch00a}. 

A related reason for the overbinding nature is the well-known
underestimation of the band gap by LDA.  From our LDA based calculations,
we find that, for example, the Si with a vacancy has a band gap of only
0.88 eV in contrast to experimental 
values of 1.12 eV.  (We note that the band gap in pure Si from LDA
calculations is even less: 0.56 eV (Figs. \ref{fig:si2},
\ref{fig:si64}). The system with 
the impurity has a 
higher band gap because of the interaction of the defect with the other
atoms in the lattice. This can be seen from a simple tight binding
analysis. We consider the band gap of the
system with the defect because that is the environment that the defect
``perceives''.)  Because of the small band gap, the defect states that
are (usually) created in the band gap tend to interact more strongly with the
valence states. This therefore leads to overbinding too. The band
structure and the density of states of the system with a vacancy,
hexagonal interstitial and split $\langle 110 \rangle$ interstitial are
shown respectively in Figs. \ref{fig:si63}, \ref{fig:si65h}, and
\ref{fig:si65s}.
 
\input{resultablelanl}

There are a couple of interrelated reasons for the LDA's underestimation
of the band gap.  The discontinuity of the one-electron potential for
localized states, which is a characterisitc of the exact density functional, is absent
in the LDA as was shown by Perdew {\em et. al.} \cite{perd82a}. Because
the band gap can be expressed as
\begin{eqnarray}
E_{gap} = E[N+1]+E[N-1]-2E[N]
\end{eqnarray}
where $E[x]$ is the energy of the system with $x$ electrons,
the absence of the discontinuity of one-electron potential in the LDA
causes the LDA to give an incorrect band gap.  The second reason for the
underestimation of band 
gap within LDA is the absence of self interaction correction.  The third
reason has to do with the Kohn-Sham approach \cite{kohn65a} itself.  In
the Kohn-Sham theory, there is no direct relationship between the orbital
energies and the ionization energies (unlike that in Hartree-Fock theory)
except for the highest occupied orbital.
              
\begin{figure}
\includegraphics[width=8.5cm]{si2}
\caption{\label{fig:si2} {\small Si band structure computed using the LDA
(solid line) and the LDA+U (dashed line) techniques. The (indirect) band
gap using LDA (LDA+U) is 0.56 (0.95) eV. }}
\end{figure}

\subsection{\label{sec:ldauresult}Local density approximation with Hubbard type correction: LDA+U}

The LDA+U method has been widely applied to metallic systems containing
localized electrons.  We refer the reader to the review article on LDA+U
by Anisimov {\em et. al.} \cite{anis97a}.  Here we provide a brief
explanation of 
the technique along the lines of ``rectifying'' the deficiencies of the
local density approximation that were pointed out in Sec. \ref{sec:ldaresult}
. 
The basic idea of the LDA+U method is to redistribute the
electron-electron interaction.  Because, as we explained in Sec. \ref{sec:ldaresult},
overbinding of the localized electrons associated with the dangling bonds
is a primary reason for the underperformance of the LDA, a direct way of
rectifying this deficiency would be to add an on-site Coulomb repulsion term
while subtracting an average overall electron-electron interaction.
%
\begin{widetext}
\begin{eqnarray}
H_{LDA+U} = H_{LDA} 
          + \sum_{i=1}^{N}\sum_{j=1;j\neq i}^{N}U_{ij}
        \int |\psi_i(\bm{r}_1)|^2\;d\bm{r}_1
	\int |\psi_j(\bm{r}_2)|^2\;d\bm{r}_2
        - U\frac{N(N-1)}{2}
\end{eqnarray}
\end{widetext}
%
where $N=\sum_{i=1}^{N}\int|\psi_i(\bm{r})|^2\;d\bm{r}$ corresponds to the
total number of electrons in the system. The parameter $U$ is adjusted by
trial and error so that the resulting band gap equals the experimentally
observed band gap of the system in question.  The above technique, which
only has the Coulomb term, is
generalized to include the exchange term as well.  We refer the reader to
Refs. \onlinecite{liec95a,solo94a} for more details. The $UN(N-1)/2$ introduces the
discontinuity of the of the one-electron potential. We wish to point out
that the LDA+U does not correct for the self interaction.

As mentioned in Sec. \ref{sec:method}, we have used the scheme as implemented in VASP for
our calculations. 
Ideally, the Coulomb ($U$) and exchange ($J$) interaction parameters would be chosen
so that the band gap for the bulk system would 
be equal to the experimentally observed band gap.  However, because our
calculations are based on a supercell geometry, there is an inevitable but
unphysical interaction between the defect and the matrix.  This alters the
band gap 
%
%
%
from that in the bulk. However, because the system ``sees'' only this band
gap, we have chosen $U$ and $J$ so that the band gap in the system with
the defect is close to the experimentally observed band gap.  This, of
course, would mean that the parameter be tuned for each type of defect.
However, for the sake of demonstration, we have worked with a single set of
parameters. 

We find that the defect formation and migration energies
increase almost linearly with the difference: $U-J$. The best agreement
between the
above energies and the experimental data occurs when the calculated band
gap corresponds to the experimental one. (This occurs when we choose $U=0$
and $J=4eV$ as mentioned in Sec. \ref{sec:method}.) This supports our argument that the
correction of the density functional for the bandgap (discontinuity of the
potential on the number of electrons) should improve the results. Table \ref{tab:resultstab}
summarizes our energetics 
calculations. Figs. \ref{fig:si2}-\ref{fig:si65s} show band structures and
densities of states of the various systems considered.

\begin{figure*}
\includegraphics[width=17cm]{si64}
\caption{\label{fig:si64}{\small Band structure and density of states of
bulk Si computed using the LDA (left) and the LDA+U
(right) methods. This figure is for reference to compare with similar
figures for systems with defects.}}
\end{figure*}

\subsection{\label{sec:defstresult}Defect structure, symmetry, states and energy}
\subsubsection{\label{sec:vacresult}Vacancy} 
There is an inward relaxation of the
atoms around the vacancy site.  The symmetry of the vacancy determines the
impurity states' symmetry. The singlet s-state is located deep inside the
valence band (not shown in Fig. \ref{fig:si63}) while partly occupied triplet
p-states are located in the 
bandgap (Fig.\ref{fig:si63}). There is a splitting of p- states into twofold degenerate lower
band and a single upper band (almost unoccupied). There are a few major differences
between the LDA results and the LDA+U results. First, there is the
expected increase in 
the bandgap in the LDA+U compared to the LDA leading to the separation
between the vacancy states and valence and 
conduction bands. Second, the dispersion of vacancy p-states is
smaller in the LDA+U calculation. Thus, these states are more localized in
the LDA+U calculation. The 
dispersion cannot be removed in the supercell geometry completely, but the
LDA definitely seems to overestimate the itinerancy of electrons in the 
vacancy states. There is no noticeable Jahn-Teller distortion of the atoms
surrounding the vacancy in our 64-atom supercell geometry. This is
consistent with the observation by Zywietz {\em et. al.} \cite{zywi98a}
who suggest that a 128-atom supercell is required to observe the
distortion. The formation and the migration energy for the vacancy defect
from LDA+U calculations are 4.59 eV and 0.42 eV respectively. This is in
good agreement with experimental results \cite{ural99b,watk00a}

\begin{figure*}
\includegraphics[width=17cm]{si63}
\caption{\label{fig:si63}{\small Band structure and density of states of a
system containing a vacancy computed using the LDA (left) and the LDA+U
(right) methods.  The defect states are shown as dashed lines.  The
(indirect) band gap using the LDA (LDA+U) is 0.88 (1.21) eV and the vacancy
dispersion using the LDA (LDA+U) is 0.77 (0.71) eV. The bottom vacancy
state is doubly degenerate.}}
\end{figure*}

\subsubsection{\label{sec:hexintresult}Hexagonal interstitial}
The self interstitial sites provide
four additional electrons to the system and 
create similar ``dangling'' bonds. Like the vacancy, they create an
almost doubly degenerate band near the bottom of the band gap and a single
band near the top of the band gap (Fig. \ref{fig:si65h}). The occupied defect states are very
close to the bottom of the band gap, in agreement with the observation made by
Needs \cite{need99a} about the defect states being quite shallow.
We find that the atoms surrounding the
interstitial move outwards in contrast to the inward movement reported by
Needs \cite{need99a}.  This is probably because, while we have optimized
the lattice constant of the system with the defect, Needs has maintained
the lattice constant at the experimental value. We suspect that Needs' approach might
cause an unphysical restriction on the relaxation. In addition to opening
the band gap, the 
introduction of the on-site repulsion quite dramatically reduces the
mixing of the 
defect states with the valence states. 
The formation energy
for the hexagonal interstitial using the LDA+U method is 4.61 eV which is
in good agreement with 
experimental \cite{ural99a} and quantum Monte Carlo \cite{leun99a} results.

\begin{figure*}
\includegraphics[width=17cm]{si65h}
\caption{\label{fig:si65h}{\small Band structure and density of states of a
system containing a hexagonal interstitial computed using the LDA (left)
and the LDA+U
(right) methods.  The defect states are shown as dashed lines.  The
(indirect) band gap using the LDA (LDA+U) is 0.70 (1.01) eV.
The amount of
mixing of the defect states with the valence states drops in the LDA+U
method, signifying lesser binding.}}
\end{figure*}


\subsubsection{\label{sec:spltintresult}Split $\langle 110 \rangle$ interstitial}
The split
$\langle 110 \rangle$ interstitial creates defect states in the band gap
similar to the hexagonal interstitial excpet that there is a greater
splitting of the two lower bands (Fig. \ref{fig:si65s}). We attribute this to the lower symmetry
of the split $\langle 110 \rangle$ interstitial compared to the hexagonal
interstitial. Another difference is the location of the higher defect
state. It is not as close to the top as in the case of hexagonal
interstitial. We do not have a simple explanation for this. The
introduction of the on-site repulsion has a similar effect with respect to
reducing the mixing of the defect states with the bulk states. The defect
formation energy using the LDA+U method is 4.70 eV which is in good
agreement with the experimental \cite{ural99a} and quantum Monte Carlo
\cite{leun99a} results. The migration energy from the hexagonal to the
split $\langle 110 \rangle$ configurations of the interstitial was
computed using the nudged elastic band menthod \cite{NEB}. We get a value
of 0.44 eV using the LDA+U method. We are not aware of any experimental
result for this specific value.

\begin{figure*}
\includegraphics[width=17cm]{si65s}
\caption{\label{fig:si65s}{\small Band structure and density of states of a
system containing a split $\langle 110 \rangle$ interstitial computed
using the LDA (left) and the LDA+U
(right) methods.  The defect states are shown as dashed lines.  The
(indirect) band gap using the LDA (LDA+U) is 0.99 (1.12) eV.
The amount of mixing of the defect states with the valence
states drops in the LDA+U method signifying lesser binding.}}
\end{figure*}


\subsection{\label{sec:selfdiffresult}Self diffusion in Si}

The contribution of the mechanism to the diffusion process is determined
by the activation energy and the prefactor. Because all the three
mechanisms considered in this article (viz., vacancy, interstitial, 
and concerted exchange) involve similar number of atoms, one can, to a first
approximation, assume that all of them have similar entropic effects and
hence similar prefactors. Thus, based on our calculations of
the activation energies using the LDA+U method, the vacancy mechanism 
[which has an activation energy of 5.01 eV (sum of formation (4.59 eV) 
and migration (0.42 eV) energies)]
and the hexagonal to split $\langle 110 \rangle$ interstitialcy mechanism
[which has an activation energy of 5.14 eV (sum of formation (4.7 eV) and
migration (0.44 eV) energies)]
are quite close in activation energy (within the error of the method), and
should contribute equally to self-diffusion. The concerted exchange
mechanism (which has an activation energy of 5.82 eV) is less significant 
in this respect. Experimental results give a
similar indication \cite{ural99a,fahe89b}.

summarylanl.tex
The present work identified the cause for the poor description of point
defect energetics by the LDA. It corrected for the deficiencies of the LDA
by using the LDA+U method. This gave much better agreement of the
calculated activation energies with experimental observations. This method
can therefore be used for better description of diffusion in similar
semiconductor materials.

ICAM2003lanl.bbl
\begin{thebibliography}{30}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
  \def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
  \def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname citenamefont\endcsname\relax
  \def\citenamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem[{\citenamefont{Ramanarayanan et~al.}(2003)\citenamefont{Ramanarayanan,
  Cho, and Clemens}}]{rama03a}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Ramanarayanan}},
  \bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Cho}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{B.~M.} \bibnamefont{Clemens}},
  \bibinfo{journal}{J.\ Appl.\ Phys.} \textbf{\bibinfo{volume}{94}},
  \bibinfo{pages}{174} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Souza and Narayanan}(1998)}]{souz98a}
\bibinfo{author}{\bibfnamefont{M.~M.~D.} \bibnamefont{Souza}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{E.~M.~S.} \bibnamefont{Narayanan}},
  \bibinfo{journal}{Defect and Diffusion Forum}
  \textbf{\bibinfo{volume}{153-155}}, \bibinfo{pages}{69}
  (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Zangenberg et~al.}(2001)\citenamefont{Zangenberg,
  Hansen, Fage-Pedersen, and Larsen}}]{zang01a}
\bibinfo{author}{\bibfnamefont{N.~R.} \bibnamefont{Zangenberg}},
  \bibinfo{author}{\bibfnamefont{J.~L.} \bibnamefont{Hansen}},
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Fage-Pedersen}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.~N.}
  \bibnamefont{Larsen}}, \bibinfo{journal}{Phys.\ Rev.\ Lett.}
  \textbf{\bibinfo{volume}{87}}, \bibinfo{pages}{125901}
  (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Strohm et~al.}(2002)\citenamefont{Strohm, Voss, Frank,
  Laitinen, and R{\"a}is{\"a}nen}}]{stro02a}
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Strohm}},
  \bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Voss}},
  \bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Frank}},
  \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Laitinen}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{R{\"a}is{\"a}nen}},
  \bibinfo{journal}{Z. Metallkd.} \textbf{\bibinfo{volume}{93}},
  \bibinfo{pages}{737} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Ural et~al.}(1998)\citenamefont{Ural, Griffin, and
  Plummer}}]{ural98a}
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Ural}},
  \bibinfo{author}{\bibfnamefont{P.~B.} \bibnamefont{Griffin}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{J.~D.}
  \bibnamefont{Plummer}}, \bibinfo{journal}{Appl.\ Phys.\ Lett.}
  \textbf{\bibinfo{volume}{73}}, \bibinfo{pages}{1706} (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Ural et~al.}(1999{\natexlab{a}})\citenamefont{Ural,
  Griffin, and Plummer}}]{ural99a}
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Ural}},
  \bibinfo{author}{\bibfnamefont{P.~B.} \bibnamefont{Griffin}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{J.~D.}
  \bibnamefont{Plummer}}, \bibinfo{journal}{J.\ Appl.\ Phys.}
  \textbf{\bibinfo{volume}{85}}, \bibinfo{pages}{6440}
  (\bibinfo{year}{1999}{\natexlab{a}}).

\bibitem[{\citenamefont{Needs}(1999)}]{need99a}
\bibinfo{author}{\bibfnamefont{R.~J.} \bibnamefont{Needs}},
  \bibinfo{journal}{J.\ Phys.:\ Condens.\ Matter}
  \textbf{\bibinfo{volume}{11}}, \bibinfo{pages}{10437} (\bibinfo{year}{1999}).

\bibitem[{\citenamefont{Kresse and Hafner}(1993)}]{VASP1}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Kresse}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Hafner}},
  \bibinfo{journal}{Phys.\ Rev.\ B} \textbf{\bibinfo{volume}{47}},
  \bibinfo{pages}{558} (\bibinfo{year}{1993}).

\bibitem[{\citenamefont{Kresse and Hafner}(1994)}]{VASP2}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Kresse}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Hafner}},
  \bibinfo{journal}{Phys.\ Rev.\ B} \textbf{\bibinfo{volume}{49}},
  \bibinfo{pages}{14251} (\bibinfo{year}{1994}).

\bibitem[{\citenamefont{Kresse and
  Furthm{\"u}ller}(1996{\natexlab{a}})}]{VASP3}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Kresse}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Furthm{\"u}ller}},
  \bibinfo{journal}{Comput.\ Mater.\ Sci.} \textbf{\bibinfo{volume}{6}},
  \bibinfo{pages}{15} (\bibinfo{year}{1996}{\natexlab{a}}).

\bibitem[{\citenamefont{Kresse and
  Furthm{\"u}ller}(1996{\natexlab{b}})}]{VASP4}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Kresse}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Furthm{\"u}ller}},
  \bibinfo{journal}{Phys.\ Rev.\ B} \textbf{\bibinfo{volume}{54}},
  \bibinfo{pages}{11169} (\bibinfo{year}{1996}{\natexlab{b}}).

\bibitem[{\citenamefont{Kresse and Joubert}(1999)}]{VASPAW1}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Kresse}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Joubert}},
  \bibinfo{journal}{Phys.\ Rev.\ B} \textbf{\bibinfo{volume}{59}},
  \bibinfo{pages}{1758} (\bibinfo{year}{1999}).

\bibitem[{\citenamefont{Bl{\"o}chl}(1994)}]{VASPAW2}
\bibinfo{author}{\bibfnamefont{P.~E.} \bibnamefont{Bl{\"o}chl}},
  \bibinfo{journal}{Phys.\ Rev.\ B} \textbf{\bibinfo{volume}{50}},
  \bibinfo{pages}{17953} (\bibinfo{year}{1994}).

\bibitem[{\citenamefont{J{\'o}nsson et~al.}(1998)\citenamefont{J{\'o}nsson,
  Mills, and Jacobsen}}]{NEB}
\bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{J{\'o}nsson}},
  \bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Mills}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{K.~W.} \bibnamefont{Jacobsen}}, in
  \emph{\bibinfo{booktitle}{Classical and Quantum Dynamics in Condensed Phase
  Simulations}}, edited by \bibinfo{editor}{\bibfnamefont{B.~J.}
  \bibnamefont{Berne}},
  \bibinfo{editor}{\bibfnamefont{G.}~\bibnamefont{Ciccotti}}, \bibnamefont{and}
  \bibinfo{editor}{\bibfnamefont{D.~F.} \bibnamefont{Coker}}
  (\bibinfo{publisher}{World Scientific}, \bibinfo{address}{Singapore},
  \bibinfo{year}{1998}), chap.~\bibinfo{chapter}{16}.

\bibitem[{\citenamefont{Monkhorst and Pack}(1976)}]{monk76a}
\bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Monkhorst}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Pack}},
  \bibinfo{journal}{Phys.\ Rev.\ B} \textbf{\bibinfo{volume}{13}},
  \bibinfo{pages}{5188} (\bibinfo{year}{1976}).

\bibitem[{\citenamefont{Liechtenstein et~al.}(1995)\citenamefont{Liechtenstein,
  Anisimov, and Zaanen}}]{liec95a}
\bibinfo{author}{\bibfnamefont{A.~I.} \bibnamefont{Liechtenstein}},
  \bibinfo{author}{\bibfnamefont{V.~I.} \bibnamefont{Anisimov}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Zaanen}},
  \bibinfo{journal}{Phys.\ Rev.\ B} \textbf{\bibinfo{volume}{52}},
  \bibinfo{pages}{R5467} (\bibinfo{year}{1995}).

\bibitem[{\citenamefont{Probert and Payne}(2003)}]{prob03a}
\bibinfo{author}{\bibfnamefont{M.~I.~J.} \bibnamefont{Probert}}
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{M.~C.} \bibnamefont{Payne}},
  \bibinfo{journal}{Phys.\ Rev.\ B} \textbf{\bibinfo{volume}{67}},
  \bibinfo{pages}{075204} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Puska et~al.}(1998)\citenamefont{Puska, P{\"o}ykk{\"o},
  Pesola, and Nieminen}}]{pusk98a}
\bibinfo{author}{\bibfnamefont{M.~J.} \bibnamefont{Puska}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{P{\"o}ykk{\"o}}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Pesola}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{R.~M.} \bibnamefont{Nieminen}},
  \bibinfo{journal}{Phys.\ Rev.\ B} \textbf{\bibinfo{volume}{58}},
  \bibinfo{pages}{1318} (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Koch and Holthausen}(2000)}]{koch00a}
\bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Koch}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.~C.} \bibnamefont{Holthausen}},
  \emph{\bibinfo{title}{A chemist's guide to density functional theory}}
  (\bibinfo{publisher}{Wiley-VCH}, \bibinfo{address}{New York},
  \bibinfo{year}{2000}), \bibinfo{edition}{2nd} ed.

\bibitem[{\citenamefont{Leung et~al.}(1999)\citenamefont{Leung, Needs,
  Rajagopal, Itoh, and Ihara}}]{leun99a}
\bibinfo{author}{\bibfnamefont{W.~K.} \bibnamefont{Leung}},
  \bibinfo{author}{\bibfnamefont{R.~J.} \bibnamefont{Needs}},
  \bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Rajagopal}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Itoh}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Ihara}},
  \bibinfo{journal}{Phys.\ Rev.\ Lett.} \textbf{\bibinfo{volume}{83}},
  \bibinfo{pages}{2351} (\bibinfo{year}{1999}).

\bibitem[{\citenamefont{Kittel}(1996)}]{kitt96a}
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Kittel}},
  \emph{\bibinfo{title}{Introduction to Solid State Physics}}
  (\bibinfo{publisher}{John Wiley \& Sons, Inc.}, \bibinfo{address}{New York},
  \bibinfo{year}{1996}), \bibinfo{edition}{seventh} ed.

\bibitem[{\citenamefont{Ural et~al.}(1999{\natexlab{b}})\citenamefont{Ural,
  Griffin, and Plummer}}]{ural99b}
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Ural}},
  \bibinfo{author}{\bibfnamefont{P.~B.} \bibnamefont{Griffin}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{J.~D.}
  \bibnamefont{Plummer}}, \bibinfo{journal}{Physica B}
  \textbf{\bibinfo{volume}{273-274}}, \bibinfo{pages}{512}
  (\bibinfo{year}{1999}{\natexlab{b}}).

\bibitem[{\citenamefont{Perdew et~al.}(1982)\citenamefont{Perdew, Parr, Levy,
  and Jose L.~Balduz}}]{perd82a}
\bibinfo{author}{\bibfnamefont{J.~P.} \bibnamefont{Perdew}},
  \bibinfo{author}{\bibfnamefont{R.~G.} \bibnamefont{Parr}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Levy}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Jose L.~Balduz}},
  \bibinfo{journal}{Phys.\ Rev.\ Lett.} \textbf{\bibinfo{volume}{49}},
  \bibinfo{pages}{1691} (\bibinfo{year}{1982}).

\bibitem[{\citenamefont{Kohn and Sham}(1965)}]{kohn65a}
\bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Kohn}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{L.~J.} \bibnamefont{Sham}},
  \bibinfo{journal}{Phys. Rev.} \textbf{\bibinfo{volume}{140}},
  \bibinfo{pages}{A1133} (\bibinfo{year}{1965}).

\bibitem[{\citenamefont{Anisimov et~al.}(1997)\citenamefont{Anisimov,
  Aryasetiawan, and Lichtenstein}}]{anis97a}
\bibinfo{author}{\bibfnamefont{V.~I.} \bibnamefont{Anisimov}},
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Aryasetiawan}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.~I.}
  \bibnamefont{Lichtenstein}}, \bibinfo{journal}{J.\ Phys.:\ Condens.\ Matter}
  \textbf{\bibinfo{volume}{9}}, \bibinfo{pages}{767} (\bibinfo{year}{1997}).

\bibitem[{\citenamefont{Solovyev et~al.}(1994)\citenamefont{Solovyev,
  Dederichs, and Anisimov}}]{solo94a}
\bibinfo{author}{\bibfnamefont{I.~V.} \bibnamefont{Solovyev}},
  \bibinfo{author}{\bibfnamefont{P.~H.} \bibnamefont{Dederichs}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{V.~I.}
  \bibnamefont{Anisimov}}, \bibinfo{journal}{Phys.\ Rev.\ B}
  \textbf{\bibinfo{volume}{50}}, \bibinfo{pages}{16861} (\bibinfo{year}{1994}).

\bibitem[{\citenamefont{Zywietz et~al.}(1998)\citenamefont{Zywietz,
  Furthm{\"{u}}ller, and Bechstedt}}]{zywi98a}
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Zywietz}},
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Furthm{\"{u}}ller}},
  \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Bechstedt}},
  \bibinfo{journal}{Phys. Stat. Sol. (b)} \textbf{\bibinfo{volume}{210}},
  \bibinfo{pages}{13} (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Watkins}(2000)}]{watk00a}
\bibinfo{author}{\bibfnamefont{G.~D.} \bibnamefont{Watkins}},
  \bibinfo{journal}{Materials Science in Semiconductor Processing}
  \textbf{\bibinfo{volume}{3}}, \bibinfo{pages}{227} (\bibinfo{year}{2000}).

\bibitem[{\citenamefont{Fahey et~al.}(1989)\citenamefont{Fahey, Iyer, and
  Scilla}}]{fahe89b}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Fahey}},
  \bibinfo{author}{\bibfnamefont{S.~S.} \bibnamefont{Iyer}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{G.~J.} \bibnamefont{Scilla}},
  \bibinfo{journal}{Appl.\ Phys.\ Lett.} \textbf{\bibinfo{volume}{54}},
  \bibinfo{pages}{843} (\bibinfo{year}{1989}).

\bibitem[{\citenamefont{Murnaghan}(1944)}]{murn44a}
\bibinfo{author}{\bibfnamefont{F.~D.} \bibnamefont{Murnaghan}},
  \bibinfo{journal}{Proc.\ Natl.\ Acad.\ Sci.} \textbf{\bibinfo{volume}{30}},
  \bibinfo{pages}{244} (\bibinfo{year}{1944}).

\end{thebibliography}
