\begin{document}
-\title{Combined ab initio and classical potential simulation study on the silicon carbide precipitation in silicon}
+\title{Combined {\em ab initio} and classical potential simulation study on the silicon carbide precipitation in silicon}
\author{F. Zirkelbach}
\author{B. Stritzker}
\affiliation{Experimentalphysik IV, Universit\"at Augsburg, 86135 Augsburg, Germany}
Much progress has been made in 3C-SiC thin film growth by chemical vapor deposition (CVD) and molecular beam epitaxy (MBE) on hexagonal SiC\cite{powell90,fissel95,fissel95_apl} and Si\cite{nishino83,nishino87,kitabatake93,fissel95_apl} substrates.
However, the frequent occurrence of defects such as twins, dislocations and double position boundaries remains a challenging problem.
Apart from these methods, high-dose carbon implantation into crystalline silicon (c-Si) with subsequent or in situ annealing was found to result in SiC microcrystallites in Si\cite{borders71}.
-Utilized and enhanced, ion beam synthesis (IBS) has become a promising method to form thin SiC layers of high quality and exclusively of the 3C polytype embedded in and epitactically aligned to the Si host featuring a sharp interface\cite{lindner99,lindner01,lindner02}.
+Utilized and enhanced, ion beam synthesis (IBS) has become a promising method to form thin SiC layers of high quality and exclusively of the 3C polytype embedded in and epitaxially aligned to the Si host featuring a sharp interface\cite{lindner99,lindner01,lindner02}.
However, the process of the formation of SiC precipitates in Si during C implantation is not yet fully understood.
\begin{figure}
A lot of theoretical work has been done on intrinsic point defects in Si\cite{bar-yam84,bar-yam84_2,car84,batra87,bloechl93,tang97,leung99,colombo02,goedecker02,al-mushadani03,hobler05,sahli05,posselt08,ma10}, threshold displacement energies in Si\cite{mazzarolo01,holmstroem08} important in ion implantation, C defects and defect reactions in Si\cite{tersoff90,dal_pino93,capaz94,burnard93,leary97,capaz98,zhu98,mattoni2002,park02,jones04}, the SiC/Si interface\cite{chirita97,kitabatake93,cicero02,pizzagalli03} and defects in SiC\cite{bockstedte03,rauls03a,gao04,posselt06,gao07}.
However, none of the mentioned studies comprehenisvely investigates all the relevant defect structures and reactions concentrated on the specific problem of 3C-SiC formation in C implanted Si.
% but mattoni2002 actually did a lot. maybe this should be mentioned!^M
-In fact, in a combined analytical potential molecular dynamics and ab initio study\cite{mattoni2002} the interaction of substitutional C with Si self-interstitials and C interstitials is evaluated.
+In fact, in a combined analytical potential molecular dynamics and {\em ab initio} study\cite{mattoni2002} the interaction of substitutional C with Si self-interstitials and C interstitials is evaluated.
However, investigations are, first of all, restricted to interaction chains along the \hkl[1 1 0] and \hkl[-1 1 0] direction, secondly lacking combinations of C interstitials and, finally, not considering migration barriers providing further information on the probability of defect agglomeration.
In particular, molecular dynamics (MD) constitutes a suitable technique to investigate their dynamical and structural properties.
It was shown that the Tersoff potential properly describes binding energies of combinations of C defects in Si\cite{mattoni2002}.
However, investigations of brittleness in covalent materials\cite{mattoni2007} identified the short range character of these potentials to be responsible for overestimated forces necessary to snap the bond of two neighbored atoms.
In a previous study\cite{zirkelbach10}, we determined the influence on the migration barrier for C diffusion in Si.
-Using the Erhart/Albe (EA) potential\cite{albe_sic_pot}, an overestimated barrier height compared to ab initio calculations and experiment is obtained.
+Using the Erhart/Albe (EA) potential\cite{albe_sic_pot}, an overestimated barrier height compared to {\em ab initio} calculations and experiment is obtained.
A proper description of C diffusion, however, is crucial for the problem under study.
-In this work, a combined ab initio and empirical potential simulation study on the initially mentioned SiC precipitation mechanism has been performed.
+In this work, a combined {\em ab initio} and empirical potential simulation study on the initially mentioned SiC precipitation mechanism has been performed.
%
By first-principles atomistic simulations this work aims to shed light on basic processes involved in the precipitation mechanism of SiC in Si.
During implantation defects such as vacancies (V), substitutional C (C$_{\text{s}}$), interstitial C (C$_{\text{i}}$) and Si self-interstitials (Si$_{\text{i}}$) are created, which play a decisive role in the precipitation process.
\section{Methodology}
\label{meth}
% ----- DFT ------
-The first-principles DFT calculations have been performed with the plane-wave based Vienna ab initio Simulation package (VASP)\cite{kresse96}.
+The first-principles DFT calculations have been performed with the plane-wave based Vienna {\em ab initio} Simulation package (VASP)\cite{kresse96}.
The Kohn-Sham equations were solved using the generalized-gradient exchange-correlation functional approximation proposed by Perdew and Wang\cite{perdew86,perdew92}.
The electron-ion interaction is described by norm-conserving ultra-soft pseudopotentials\cite{hamann79} as implemented in VASP\cite{vanderbilt90}.
Throughout this work, an energy cut-off of \unit[300]{eV} was used to expand the wave functions into the plane-wave basis.
\multicolumn{10}{l}{Present study} \\
VASP & 3.39 & 3.42 & 3.77 & 4.41 & 3.63 & 1.95 & 3.72 & 4.16 & 4.66 \\
Erhart/Albe & 4.39 & 4.48$^*$ & 3.40 & 5.42 & 3.13 & 0.75 & 3.88 & 5.18 & 5.59$^*$ \\
- \multicolumn{10}{l}{Other ab initio studies} \\
+ \multicolumn{10}{l}{Other {\em ab initio} studies} \\
Ref.\cite{al-mushadani03} & 3.40 & 3.45 & - & - & 3.53 & - & - & - & - \\
Ref.\cite{leung99} & 3.31 & 3.31 & 3.43 & - & - & - & - & - & - \\
Ref.\cite{dal_pino93,capaz94} & - & - & - & - & - & 1.89\cite{dal_pino93} & x & - & x+2.1\cite{capaz94}
\end{tabular}
\end{ruledtabular}
-\caption{Formation energies of C and Si point defects in c-Si determined by classical potential and ab initio methods. The formation energies are given in electron volts. T denotes the tetrahedral and BC the bond-centered configuration. Subscript i and s indicates the interstitial and substitutional configuration. Dumbbell configurations are abbreviated by DB. Formation energies for unstable configurations obtained by classical potential MD are marked by an asterisk and determined by using the low kinetic energy configuration shortly before the relaxation into the more favorable configuration starts.}
+\caption{Formation energies of C and Si point defects in c-Si determined by classical potential and {\em ab initio} methods. The formation energies are given in electron volts. T denotes the tetrahedral and BC the bond-centered configuration. Subscript i and s indicates the interstitial and substitutional configuration. Dumbbell configurations are abbreviated by DB. Formation energies for unstable configurations obtained by classical potential MD are marked by an asterisk and determined by using the low kinetic energy configuration shortly before the relaxation into the more favorable configuration starts.}
\label{table:sep_eof}
\end{table*}
Although discrepancies exist, classical potential and first-principles methods depict the correct order of the formation energies with regard to C defects in Si.
Substitutional C (C$_{\text{s}}$) constitutes the energetically most favorable defect configuration.
Since the C atom occupies an already vacant Si lattice site, C$_{\text{s}}$ is not an interstitial defect.
-The quantum-mechanical result agrees well with the result of another ab initio study\cite{dal_pino93}.
+The quantum-mechanical result agrees well with the result of another {\em ab initio} study\cite{dal_pino93}.
Clearly, the empirical potential underestimates the C$_{\text{s}}$ formation energy.
The C interstitial defect with the lowest energy of formation has been found to be the C-Si \hkl<1 0 0> interstitial dumbbell (C$_{\text{i}}$ \hkl<1 0 0> DB), which, thus, constitutes the ground state of an additional C impurity in otherwise perfect c-Si.
This finding is in agreement with several theoretical\cite{burnard93,leary97,dal_pino93,capaz94} and experimental\cite{watkins76,song90} investigations.
Due to the high formation energy of the BC defect resulting in a low probability of occurrence of this defect, the wrong description is not posing a serious limitation of the EA potential.
A more detailed discussion of C defects in Si modeled by EA and DFT including further defect configurations can be found in our recently published article\cite{zirkelbach10}.
-Regarding intrinsic defects in Si, classical potential and {\em ab initio} methods predict energies of formation that are within the same order of magnitude.
+Regarding intrinsic defects in Si, classical potential and {\em {\em ab initio}} methods predict energies of formation that are within the same order of magnitude.
%
However discrepancies exist.
Quantum-mechanical results reveal the Si$_{\text{i}}$ \hkl<1 1 0> DB to compose the energetically most favorable configuration closely followed by the hexagonal and tetrahedral configuration, which is the consensus view for Si$_{\textrm{i}}$ and compares well to results from literature\cite{leung99,al-mushadani03}.
The activation energy of approximately \unit[2.24]{eV} is needed to turn the C$_{\text{i}}$ \hkl[0 0 -1] DB into the C$_{\text{i}}$ \hkl[1 1 0] DB located at the neighbored lattice site in \hkl[1 1 -1] direction.
Another barrier of \unit[0.90]{eV} exists for the rotation into the C$_{\text{i}}$ \hkl[0 -1 0] DB configuration for the path obtained with a time constant of \unit[100]{fs} for the Berendsen thermostat.
Roughly the same amount would be necessary to excite the C$_{\text{i}}$ \hkl[1 1 0] DB to the BC configuration (\unit[0.40]{eV}) and a successive migration into the \hkl[0 0 1] DB configuration (\unit[0.50]{eV}) as displayed in our previous study\cite{zirkelbach10}.
-The former diffusion process, however, would more nicely agree with the ab initio path, since the migration is accompanied by a rotation of the DB orientation.
-By considering a two step process and assuming equal preexponential factors for both diffusion steps, the probability of the total diffusion event is given by $\exp\left((\unit[2.24]{eV}+\unit[0.90]{eV})/{k_{\text{B}}T}\right)$, which corresponds to a single diffusion barrier that is 3.5 times higher than the barrier obtained by ab initio calculations.
+The former diffusion process, however, would more nicely agree with the {\em ab initio} path, since the migration is accompanied by a rotation of the DB orientation.
+By considering a two step process and assuming equal preexponential factors for both diffusion steps, the probability of the total diffusion event is given by $\exp\left((\unit[2.24]{eV}+\unit[0.90]{eV})/{k_{\text{B}}T}\right)$, which corresponds to a single diffusion barrier that is 3.5 times higher than the barrier obtained by {\em ab initio} calculations.
Accordingly, the effective barrier of migration of C$_{\text{i}}$ is overestimated by a factor of 2.4 to 3.5 compared to the highly accurate quantum-mechanical methods.
This constitutes a serious limitation that has to be taken into account for modeling the C-Si system using the otherwise quite promising EA potential.
Again a single bond switch, i.e. the breaking of the bond of the Si atom bound to the fourfold coordinated C$_{\text{s}}$ atom and the formation of a double bond between the two C atoms, results in configuration b.
The two C atoms form a \hkl[1 0 0] DB sharing the initial C$_{\text{s}}$ lattice site while the initial Si DB atom occupies its previously regular lattice site.
The transition is accompanied by a large gain in energy as can be seen in Fig.~\ref{fig:026-128}, making it the ground state configuration of a C$_{\text{s}}$ and C$_{\text{i}}$ DB in Si yet \unit[0.33]{eV} lower in energy than configuration B.
-This finding is in good agreement with a combined ab initio and experimental study of Liu et~al.\cite{liu02}, who first proposed this structure as the ground state identifying an energy difference compared to configuration B of \unit[0.2]{eV}.
+This finding is in good agreement with a combined {\em ab initio} and experimental study of Liu et~al.\cite{liu02}, who first proposed this structure as the ground state identifying an energy difference compared to configuration B of \unit[0.2]{eV}.
% mattoni: A favored by 0.2 eV - NO! (again, missing spin polarization?)
A net magnetization of two spin up electrons, which are equally localized as in the Si$_{\text{i}}$ \hkl<1 0 0> DB structure is observed.
In fact, these two configurations are very similar and are qualitatively different from the C$_{\text{i}}$ \hkl<1 0 0> DB that does not show magnetization but a nearly collinear bond of the C DB atom to its two neighbored Si atoms while the Si DB atom approximates \unit[120]{$^{\circ}$} angles in between its bonds.
Similar to what was previously mentioned, configurations of C$_{\text{s}}$ and a Si$_{\text{i}}$ DB might be particularly important at higher temperatures due to the low activation energy necessary for its formation.
At higher temperatures the contribution of entropy to structural formation increases, which might result in a spatial separation even for defects located within the capture radius.
-Indeed, an ab initio molecular dynamics run at \unit[900]{$^{\circ}$C} starting from configuration \RM{1}, which -- based on the above findings -- is assumed to recombine into the ground state configuration, results in a separation of the C$_{\text{s}}$ and Si$_{\text{i}}$ DB by more than 4 neighbor distances realized in a repeated migration mechanism of annihilating and arising Si$_{\text{i}}$ DBs.
+Indeed, an {\em ab initio} molecular dynamics run at \unit[900]{$^{\circ}$C} starting from configuration \RM{1}, which -- based on the above findings -- is assumed to recombine into the ground state configuration, results in a separation of the C$_{\text{s}}$ and Si$_{\text{i}}$ DB by more than 4 neighbor distances realized in a repeated migration mechanism of annihilating and arising Si$_{\text{i}}$ DBs.
The atomic configurations for two different points in time are shown in Fig.~\ref{fig:md}.
Si atoms 1 and 2, which form the initial DB, occupy Si lattice sites in the final configuration while Si atom 3 is transferred from a regular lattice site into the interstitial lattice.
\begin{figure}
$t=\unit[2900]{fs}$
\end{center}
\end{minipage}
-\caption{Atomic configurations of an ab initio molecular dynamics run at \unit[900]{$^{\circ}$C} starting from a configuration of C$_{\text{s}}$ located next to a Si$_{\text{i}}$ \hkl[1 1 0] DB (atoms 1 and 2). Equal atoms are marked by equal numbers. Bonds are drawn for substantial atoms only.}
+\caption{Atomic configurations of an {\em ab initio} molecular dynamics run at \unit[900]{$^{\circ}$C} starting from a configuration of C$_{\text{s}}$ located next to a Si$_{\text{i}}$ \hkl[1 1 0] DB (atoms 1 and 2). Equal atoms are marked by equal numbers. Bonds are drawn for substantial atoms only.}
\label{fig:md}
\end{figure}
Separated configurations of \cs{} and \si{} become even more likely if Si diffusion exhibits a low barrier of migration.
Concerning the mobility of the ground state Si$_{\text{i}}$, an activation energy of \unit[0.67]{eV} for the transition of the Si$_{\text{i}}$ \hkl[0 1 -1] to \hkl[1 1 0] DB located at the neighbored Si lattice site in \hkl[1 1 -1] direction is obtained by first-principles calculations.
Further quantum-mechanical investigations revealed a barrier of \unit[0.94]{eV} for the Si$_{\text{i}}$ \hkl[1 1 0] DB to Si$_{\text{i}}$ H, \unit[0.53]{eV} for the Si$_{\text{i}}$ \hkl[1 1 0] DB to Si$_{\text{i}}$ T and \unit[0.35]{eV} for the Si$_{\text{i}}$ H to Si$_{\text{i}}$ T transition.
-These are of the same order of magnitude than values derived from other ab initio studies\cite{bloechl93,sahli05}.
+These are of the same order of magnitude than values derived from other {\em ab initio} studies\cite{bloechl93,sahli05}.
The low barriers indeed enable configurations of further separated \cs{} and \si{} atoms by the highly mobile \si{} atom departing from the \cs{} defect as observed in the previously discussed MD simulation.
\subsection{Summary}
Erhart/Albe & 3.88 & 4.93 & 5.25$^{\text{a}}$/5.08$^{\text{b}}$/4.43$^{\text{c}}$
\end{tabular}
\end{ruledtabular}
-\caption{Formation energies of defect configurations of a single C impurity in otherwise perfect c-Si determined by classical potential and ab initio methods. The formation energies are given in electron volts. T denotes the tetrahedral and the subscripts i and s indicate the interstitial and substitutional configuration. Superscripts a, b and c denote configurations of C$_{\text{s}}$ located at the first, second and third nearest neighbored lattice site with respect to the Si$_{\text{i}}$ atom.}
+\caption{Formation energies of defect configurations of a single C impurity in otherwise perfect c-Si determined by classical potential and {\em ab initio} methods. The formation energies are given in electron volts. T denotes the tetrahedral and the subscripts i and s indicate the interstitial and substitutional configuration. Superscripts a, b and c denote configurations of C$_{\text{s}}$ located at the first, second and third nearest neighbored lattice site with respect to the Si$_{\text{i}}$ atom.}
\label{tab:defect_combos}
\end{table}
Obviously the EA potential properly describes the relative energies of formation.
However, this configuration is unstable involving a structural transition into the C$_{\text{i}}$ \hkl<1 1 0> interstitial, thus, not maintaining the tetrahedral Si nor the substitutional C defect.
Thus, the underestimated energy of formation of C$_{\text{s}}$ within the EA calculation does not pose a serious limitation in the present context.
-Since C is introduced into a perfect Si crystal and the number of particles is conserved in simulation, the creation of C$_{\text{s}}$ is accompanied by the creation of Si$_{\text{i}}$, which is energetically less favorable than the ground state, i.e. the C$_{\text{i}}$ \hkl<1 0 0> DB configuration, for both, the EA and ab initio treatment.
+Since C is introduced into a perfect Si crystal and the number of particles is conserved in simulation, the creation of C$_{\text{s}}$ is accompanied by the creation of Si$_{\text{i}}$, which is energetically less favorable than the ground state, i.e. the C$_{\text{i}}$ \hkl<1 0 0> DB configuration, for both, the EA and {\em ab initio} treatment.
In either case, no configuration more favorable than the C$_{\text{i}}$ \hkl<1 0 0> DB has been found.
Thus, a proper description with respect to the relative energies of formation is assumed for the EA potential.
\subsection{Summary}
Investigations are targeted at the initially stated controversy of SiC precipitation, i.e. whether precipitation occurs abruptly after enough C$_{\text{i}}$ agglomerated or after a successive agglomeration of C$_{\text{s}}$ on usual Si lattice sites (and Si$_{\text{i}}$) followed by a contraction into incoherent SiC.
-Results of the previous ab initio study on defects and defect combinations in C implanted Si suggest C$_{\text{s}}$ to play a decisive role in the precipitation of SiC in Si.
+Results of the previous {\em ab initio} study on defects and defect combinations in C implanted Si suggest C$_{\text{s}}$ to play a decisive role in the precipitation of SiC in Si.
To support previous assumptions MD simulations, which are capable of modeling the necessary amount of atoms, i.e. the precipitate and the surrounding c-Si structure, have been employed in the current study.
In a previous comparative study\cite{zirkelbach10} we have shown that the utilized empirical potential fails to describe some selected processes.
However, we observed a phase transition of the C$_{\text{i}}$-dominated into a clearly C$_{\text{s}}$-dominated structure.
The amount of substitutionally occupied C atoms increases with increasing temperature.
Entropic contributions are assumed to be responsible for these structures at elevated temperatures that deviate from the ground state at 0 K.
-Indeed, in the ab initio MD simulation performed at \unit[900]{$^{\circ}$C} we observed the departing of a Si$_{\text{i}}$ \hkl<1 1 0> DB located next to a C$_{\text{s}}$ atom instead of a recombination into the ground state configuration, i.e. a C$_{\text{i}}$ \hkl<1 0 0> DB.
+Indeed, in the {\em ab initio} MD simulation performed at \unit[900]{$^{\circ}$C} we observed the departing of a Si$_{\text{i}}$ \hkl<1 1 0> DB located next to a C$_{\text{s}}$ atom instead of a recombination into the ground state configuration, i.e. a C$_{\text{i}}$ \hkl<1 0 0> DB.
\section{Conclusions}
> response to all recommendations and criticisms.
We decided to follow yours and the referee's suggestion to merge the
-two manuscripts in a single comprehensive manuscript.
+two manuscripts in a single comprehensive manuscript. Also, according
+to the referee's suggestions, some points were clarified and explained
+in more detail.
Please find below the summary of changes and a detailed response to
the recommendations of the referee.
> are seen for constant volume calculations (on a few simple
> examples, say)?
-% Differences are supposed to be negligible small since only small
-% changes in volume are detected. However, in experiment, substrate
-% swelling is observed. Thus, to allow for full relaxation, simulations
-% were performed in the NpT ensemble. However, for the above-mentioned
-% reason, no fundamental differences are expected for single defect
-% configurations in the canonical and isothermal-isobaric ensemble with
-% respect to energy.
-
In experiment, substrate swelling is observed for high-dose carbon
implantation into silicon. Indeed, for a single defect, the change in
volume is less than 0.2% in simulation. Due to this, results of single
-defects within an isothermal-isobaric simulation are expected to not
-drastically differ to results of constant volume simulations. Based on the
-experimentally observed change in volume for high-dose carbon
+defects within an isothermal-isobaric simulation are not expected to
+differ drastically to results of constant volume simulations. Based on
+the experimentally observed change in volume for high-dose carbon
implantations, however, the respective relaxation is allowed for in
-simulation for both, single defect calulations as well as the high carbon
-concentration simulations.
+simulation for both, single defect calulations as well as the high
+carbon concentration simulations.
A respective statement was added to the methodology section
(Change 4).
+ = line added
- = line removed
-Change 1: added/merged parts of the Abstract of BA11443
+Change 1: added/merged parts of 'Abstract' of BA11443
from: These aime to clarify ...
until: Finally, results of the ...
-Change 2: added/merged parts of the Introduction of BA11443
+Change 2: added/merged parts of 'Introduction' of BA11443
from: A lot of theoretical work has been done ...
until: However, investigations are, first of all, ...
structures.
Change 8: added definition and explanation of the binding energy to
- the methodology section
+ the 'Methodology' section
from: The binding energy of a defect pair ...
until: The interaction strength, i.e. the ...
-Change 9: removed Results section
+Change 9: removed 'Results' section
Change 10: added 'Comparison of classical potential and
first-principles methods' section