[lectures/latex.git] / posic / publications / sic_prec_reply02.txt

Re: BC11912
    Combined ab initio and classical potential simulation study on the
    silicon carbide precipitation in silicon
    by F. Zirkelbach, B. Stritzker, K. Nordlund, et al.

and

Re: BA11443
    First-principles study of defects in carbon-implanted silicon
    by F. Zirkelbach, B. Stritzker, J. K. N. Lindner, et al.



Dear Dr. Dahal,

thank you for the feedback to our submission.

> We look forward to receiving such a comprehensive manuscript. When you
> resubmit, please include a summary of the changes made, and a detailed
> 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.

Please find below the summary of changes and a detailed response to
the recommendations of the referee.

Some arguments here were already put forward in our previous reply and
are repeated for the sake of clarity. We would be glad to comment at
length on further upcoming, more detailed questions.

Sincerely,

Frank Zirkelbach


--------------- Response to recommendations ----------------

> I am not happy with these two papers for a multitude of reasons,
> and I recommend that the authors rewrite them as a single longer
> paper, to eliminate the criticism of serial publication. I do not
> accept the authors argument that they should be two papers ­ they
> address the same issues, using the same methods. If they were to
> be split into two papers, it would be one for the VASP
> calculations, and one for the MD ­ this is not how I suggest you
> do it, though.

We now combined the two manuscripts to a single comprehensive one.

> do it, though. First, though, the following issues should be
> addressed (some are simply pasted from my previous reviews, where
> I feel that the authors have ignored them, or not responded
> adequately).
>
> 1. I feel that the authors are a bit too convinced by their own
> calculations.  They do not state the error bars that would be
> expected for calculations like this +/- 0.2 eV would be a very
> optimistic estimate, I suggest. That being so, many of their
> conclusions on which structure or migration routes are most
> likely start to look rather less certain.

In literature, very often, differences less than 0.2 eV are obtained
in DFT studies and respective conclusions are derived. For instance,
differences in the energy of formation ranging from 0.05 - 0.12 eV are
considered significant enough to conclude on the energetically most
favorable intrinsic defect configurations in Si (PRB 68, 235205
(2003); PRL 83, 2351 (1999)).  This is due to the fact that existing
errors are most probably of the systematic rather than the random
type. The error in the estimate of the cohesive energy is canceled out
since it is likewise wrong in the defect as in the bulk configuration,
which are substracted in the expression of the defect formation
energy. Even if the defect formation energy is overestimated due to a
too small size of the supercell resulting in a non-zero interaction of
the defect with its images, this is likewise true for other defects.
Although the actual value might be wrong, observed differences in
energy, thus, allow to draw conlcusions on the stability of defect
configurations. This is also valid for diffusion barriers, which are
given by differences in energy of different structures.
In fact, differences of 0.2 eV in DFT calculations are considered
insignificant when being compared to experimental results or data of
other ab initio studies. However, the observed differences in energy
within our systematic DFT study are considered reliable.

> 2. Why is 216 atoms a large enough supercell - many defect
> properties are known to converge very slowly with supercell size.

Of course, choosing a supercell containing 216 atoms constitutes a
tradeoff. It is considered the optimal choice with respect to
computational efficiency and accuracy.

We would like to point out that, both, single defects as well as
combinations of two defects were investigated in such supercells in
successive calculations.

For single defects, the size of the supercell should be sufficient.
This is shown in PRB 58, 1318 (1998) predicting convergence of the
vacancy in silicon - the defect assumed to be most critical due to
the flatness of the total energy surface as a function of the ionic
coordinates - for supercells containing more than 128 atomic sites,
where the defect formation energy is already well estimated using
smaller supercells of 64 atomic sites. Thus, convergence of the
formation energies of single defects with respect to the size of the
supercell is assumed.

A repsective statement was added (Change 3).

> They appear to be separating defects by as large a distance as
> can be accommodated in the supercell to approximate the isolated
> defects, but then they are only separated by a few lattice
> spacings from a whole array of real and image defects ­ how does
> that compare with taking the energies of each defect in a
> supercell.

We would like to remind the referee that the properties of isolated,
non-intertacting defects were modeled in separate simulation runs.  It
is not our purpose to separate defects by a large distance in order to
approximate the situation of isolated defects.  We are rather
interested in interacting defects. However, we did find that for
increasing defect distances, configurations appear, which converge to
the energetics of two isolated defects. This is indicated by the
(absolute value of the) binding energy, which is approaching zero with
increasing distance.  From this, we conclude a decrease in
interaction, which is already observable for defect separation
distances accessible in our simulations. Combinations of defects with
similar distances were already successfully modeled in a supercell
containing 216 atoms as described in PRB 66, 195214 (2002).

An explanation of the binding energy and the relation to the
interaction of defects was added (Change 8).

> 3. Constant pressure solves some problems, but creates others ­
> is it really a sensible model of implantation? What differences
> 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.

A respective statement was added to the methodology section
(Change 4).

> 4. What method do they use to determine migration paths? How can
> they convince us that the calculations cover all possible
> migrations paths ­ that is, the paths they calculate are really
> the lowest energy ones?  This is a major issue ­ there are a
> number of methods used in the literature to address it ­ are the
> authors aware of them? Have they used one of them?

The constrained relaxation technique is used to determine migration
pathways. The method is named and a reference is given in the
methodology section. The method not necessarily unveils the lowest
energy migration path. The supposed saddle point structure needs to be
attested by investigating the vibrational modes. However, reasonable
results are obtained for the specific system. In fact, so far, the
best quantitative agreement with experimental findings has been
achieved concerning the interstitial carbon mobility (PRB 82, 094110
(2010)) utilizing the constrained relaxation technique. Thus, obtained
migration paths are assumed to be valid without investigating the
vibrational modes of every single supposed saddle point configuration.

For clarity we added a statement that, of course, the true minimum
energy path may still be missed (Change 7).

> 5. I have some serious reservations about the methodology
> employed in the MD calculations. The values given for the basic
> stabilities and migration energies in some cases disagree
> radically with those calculated by VASP, which I would argue
> (despite 4 above) to be the more reliable values. The main

Indeed, discrepancies exist. However, both methods predict the C-Si
100 DB configuration to be the ground-state structure.  The
underestimated energy of formation of substitutional C for the EA
potential does not pose a problem in the present context. Since we
deal with a perfect Si crystal and the number of particles is
conserved, the creation of substitutional C is accompanied by the
creation of a Si interstitial. The formation energies of the different
structures of an additional C atom incorporated into otherwise perfect
Si shows the same ground state, i.e. the C-Si 100 DB structure, for
classical potential as well as ab initio calculations.

This is discussed in full detail in section V in the combined
manuscript.

> problems is the huge over-estimate of the C interstitial
> migration energy (a process which is at the heart of the
> simulations) using the potential used in the paper. I am not
> convinced that the measures they take to circumvent the problems
> in the method do not introduce further uncertainties, and I would
> need a bit more convincing that the results are actually valid.

We hope to be able to convince by responding to the following
statement of the referee.

> The authors' circumvention of this is to do the simulations at
> much heightened temperatures.  However, this only gives a good
> model of the system if all cohesive and migration energies are
> over-estimated by a similar factor, which is demonstratably
> untrue in this case. For this reason, despite the reputation and
> previous work with Tersoff (and similar) potentials, the results
> need a critical scrutiny, which I am not very convinced by in
> this case.

There is not necessarily a correlation of the cohesive and migration
energies. You can always add a constant to the cohesive energies of
respective structures. It is the difference in the cohesive energies
of structures within the migration path, which determines the
migration barrier.

In fact, cohesive energies are most often well described by the
classical potentials since these are most often used to fit the
potential parameters.

The overestimated migration barrier, however, is due to the short
range character of the potential, which drops the interaction to
zero within the first and next neighbor distance using a special
cut-off function as explained in PRB 76, 224103 (2007). The
overestimated barrier and slightly different pathway (however,
starting and final configuration/orientation agree) is indeed 
demonstrated for the carbon interstitial within the present study.
Since the reason of overestimation is inherent to the short range
potential, migration pathways among other configurations are
likewise overestimated.

Since most of the defect structures show atomic distances below the
critical distance, for which the cut-off function is taking effect,
the respective formation energies are quite well described, too (at
least they are not necessarily overestimated in the same way).

Thus, increased temperatures result in an increased probability of
transition. Obviously, this enables the structural transformation
into energetically less stable structures of substitutional carbon and
interstitial silicon that are observed in the high temperature
simulations. Being in nice agreement with experimental findings, these
results suggest the usage of increased temperatures to constitute a
necessary condition to deviate the system out of the ground state as
it is the case in the ion beam synthesis process.

A respective statement and a more detailed comparison with experiment
was added to the combined version of the manuscript (Change 22).

Again, we would like to repeat the arguments that legitimate the usage
of increased temperatures although cohesive and formational energies
are not ovrestimated in the same way than the migration barriers.
While the properties of some structures near the equilibrium position
are well described, the above mentioned effects increase for
non-equilibrium structures and dynamics. Thus, for instance, it is not
surprising that short range potentials show overestimated melting
temperatures. This is not only true for the EA but also (to an even
larger extent) for Tersoff potentials, one of the most widely used
classical potentials for the Si/C system. The fact that the melting
temperature is drastically overestimated although the cohesive
energies are nicely reproduced indicates that there is no reason why
the cohesive and formational energies should be overestimated to the
same extent in order to legitimate the increase in temperature to
appropriately consider the overestimated barrier heights for
diffusion.

Indeed the cut-off effect increases if the system is driven away from
the equilibrium (such as by modeling IBS). Since this is to some
extent cured by increasing the simulation temperature, the work-around
is particularly helpful for short range potentials.


--------------- Summary of changes ----------------

Since the new manuscript is a combination of manuscripts BC11912 and
BA11443, the following summary of changes mainly contains the
construction of the new manuscript by text blocks of previous
manuscripts. Please let me know if a more detailed summary of changes
is required.

The title of the new manuscript is that of BC11912. Thus, stated
changes apply to this manuscript.

Description:
+ = line added
- = line removed

Change 1: added/merged parts of the Abstract of BA11443

 from:  These aime to clarify ...
 until: Finally, results of the ...

Change 2: added/merged parts of the Introduction of BA11443

 from:  A lot of theoretical work has been done ...
 until: However, investigations are, first of all, ... 

 from:  By first-principles atomistic simulations ...
 until: Furthermore, highly accurate quantum-mechanical ...

Change 3: convergence of BZ sampling and size of the supercell

-Sampling of the Brillouin zone was restricted to the $\Gamma$-point.
-The defect structures and the migration paths have been modeled in
 cubic supercells containing 216 Si atoms.
+To reduce the computational effort sampling of the Brillouin zone was
 restricted to the $\Gamma$-point, which has been shown to yield
 reliable results\cite{dal_pino93}.
+The defect structures and the migration paths were modelled in cubic
 supercells with a side length of \unit[1.6]{nm} containing $216$ Si
 atoms.
+Formation energies and structures are reasonably converged with
 respect to the system size.

Change 4: only small changes in volume

+The observed changes in volume were less than \unit[0.2]{\%} of the
 volume indicating a rather low dependence of the results on the
 ensemble choice.

Change 5: name algorithm used for structural relaxation
          in DFT calculations
 
+Ionic relaxation was realized by the conjugate gradient algorithm.

Change 6: name reason for reservoir choice

+This corresponds to the definition utilized in another study on C
 defects in Si\cite{dal_pino93} that we compare our results to. 

Change 7: CRT not necessarily predicts the minimum energy path

+While not guaranteed to find the true minimum energy path, the method
 turns out to identify reasonable pathways for the investigated
 structures.

Change 8: added definition and explanation of the binding energy to
          the methodology section

 from:  The binding energy of a defect pair ...
 until: The interaction strength, i.e. the ...

Change 9: removed Results section

Change 10: added 'Comparison of classical potential and
           first-principles methods' section

+In a first step, quantum-mechanical calculations of defects in Si and
 respective diffusion processes are compared to classical potential
 simulations as well as to results from literature.
+Shortcomings of the analytical potential approach are revealed and
 its applicability is discussed.

Change 11: comprehensive Table including all defects and methods

Change 12: added text on unstable hexagonal Si defect for classical
           potentials - necessary due to combination of manuscripts!

 from:  The hexagonal configuration ...
 until: While not completely rendering impossible ...

Change 13: added configurations that require spin polarized
           calculations

 from:  Instead of giving an explicit value ...
 until: No other configuration, within ...

Change 14: 'Carbon mobility' section of BC11912 mapped to 'Mobility of
           carbon defects' section

Change 15: added 'Quantum-mechanical investigations of defect
           combinations and related diffusion processes' section
           corresponding to 'Results' section of BA11443

Change 16: added 'Mobility of silicon defects" section from III A of
           BA11443

Change 17: added 'Summary' section from 'Discussion' section of
           BA11443

Change 18: relocate 'Excursus: Competition of C_i and C_s-Si_i' section
           of BC11912

Change 19: section 'Classical potential calculations on the SiC
           precipitation in Si' and respective glue text added

 from:  The MD technique is used to gain ...
 until: The approach is follwed and, ...
 
 content corresponds to 'Results' section of BC11912

Change 20: 'Summary' section added containing parts of 'Discussion and
            summary' section of BC11912

Change 21: 'Conclusions' section added containing parts of the
           'Discussion' section of BA11443 and the 'Discussion and
           summary' section of BC11912

Change 22: more detailed comparison to experiment added

 starting from:  Moreover, results of the MD simulations ...

Change 23: 'Summary' section added containing parts of the 'Summary'
           section of BA11443 and the 'Discussion and summary' section
           of BC11912