Temperature dependence of thermodynamic quantities can be calculated from simulations performed with temperature[111, 112]. For example, the free energy of a solid can be determined from the local atomic configuration, hence a minimisation of the free energy with respect to atomic coordinates give both the equilibrium structure and free energy of the solid which could contain defects. This requires thermal statistics to be gathered in either a molecular dynamics or Monte Carlo simulation to calculate the canonical average of the thermal property.
It is not necessary to do this here since the phonon modes for a given defect concentration have been calculated directly for a given defect concentration (which is inversely proportional to the size of the supercell containing the defect). Thermodynamic properties using the vibrational spectra for the crystal can be found simply from the partition function calculated from Bose-Einstein statistics as demonstrated in Chapter 4 where a P-T phase diagram was found. The partition function allows the vibrational free energy and vibrational entropy to be calculated.
By taking the differences between the perfect crystal and the crystals containing defects (scaled to the same number of atoms) the effects that the defects have on the thermal properties can be seen. The plots of these excess free energies, total vibrational energy and vibrational entropies varying with temperature are shown in Figure 7.11.
Figure 7.11: Differences in (a) vibrational free energy, (b) vibrational entropy and (c) total vibrational energy are shown as temperature varies. The zero of each scale is taken to be the value of that quantity of the perfect diamond structure at each temperature. (Full line - hexagonal interstitial, dashed line - tetrahedral interstitial, dotted line - vacancy).
These can be compared to Figure 7.12 which shows the total vibrational energy, free energy and entropy for the perfect crystal.
Figure 7.12: Plots of total vibrational energy, vibrational free energy and vibrational entropy per atom against temperature. Only the plots for the perfect crystal are shown since the differences in the curves for the defects are too small to be seen on this scale (see previous figure).
These quantities for the defected crystals are not shown in Figure 7.12 since they are indistinguishable on that scale. It is found that the tetrahedral interstitial has a free energy higher than that of the hexagonal interstitial at any temperature. Note that this excess free energy is times smaller than that of the defect formation energy and hence the tetrahedral interstitial is stable with respect to the hexagonal site. This leads to the conclusion that the defect formation energies for interstitials and vacancies within a crystal in the harmonic approximation are almost independent of temperature. The free energy of a silicon crystal at ambient pressure was found from the Stillinger-Weber potential by summing the logarithms of the frequencies obtained from the dynamical matrix method used here. Broughton and Li find the sum of the logarithms to be 4.7533 per particle, although a full free energy analysis with temperature was not given. A similar calculation using the dynamical matrix for the perfect crystal, described above, gives the value of 4.6442. Calculations on the hexagonal and tetrahedral interstitials and the vacancy give (in the same reduced units) 4.6449, 4.6452 and 4.6444 respectively. It can be seen that the introduction of defects to the crystal only changes the free energy slightly. It is the small difference in this free energy in addition to the difference in defect formation energies that determines properties such as diffusion in defects through the crystal.
The free energies of the defects can also be found, including the anharmonic effects, by running the molecular dynamics simulations at finite temperature over an extended period of time to collect thermal averages. These simulations comprised of the crystal with defect running for 12000 molecular dynamics time steps of 1.0fs which is longer than 10 periods of the lowest frequency vibrations. This allowed the evaluation of the anharmonic terms to a reasonable tolerance. It was found that the defect formation energy is unchanged within a temperature range of 0-1000K. This confirms the validity of obtaining the above results using the harmonic approximation.
Finally, the configurational entropy of the defects can be calculated to give a comparison to the vibrational entropy. The Boltzmann definition of entropy, , where is the number of configurations, can be used to calculate the configurational entropy per atom, after simplification by Stirling's formula, by
where m is the number of possible different defect sites and c is the defect concentration. The configurational entropies for the three defects considered are shown in Table 7.2.
Table 7.2: Configurational entropies for the large supercell and the smaller supercell. The units are in eV/K per defect.
It is found that under the harmonic approximation, the contributions to the entropy of a crystal are given by the configurational entropy and the vibrational entropy in roughly equal proportions at non-zero temperatures.
In the converged calculations the differences in the interstitial formation energy is 0.51eV (Table 7.1). From Figure 7.11, it can be seen that the vibrational free energy is not enough to transform a tetrahedral interstitial to the hexagonal site within the temperature range shown. The relative stability of the anharmonic case is tested by running a molecular dynamics simulation at increasing temperature, allowing the bonding to change throughout the simulation. A transition from one defect site to another is not found until 1360K. Higher temperatures were required to allow the interstitial atom to move freely throughout the crystal. This also indicates that the harmonic approximation taken here is quite accurate up to relatively high temperatures. Recent work using the Stillinger-Weber potential on self interstitial diffusivity in silicon[113, 105] has found a simple migration path for several self interstitials over a temperature range of 733-1473K. In the work by Maroudas and Brown the relaxation was done in a somewhat different manner from the molecular dynamics simulations here. Energy minimisation was carried out by quenching the system by a steepest descents algorithm after relaxing the system by a lengthy Monte Carlo method at 500K with system size of the host crystal ranging from 60 atoms to 512. The work of Maroudas and Brown and the calculations given here suggest that the tetrahedral interstitial configuration is the most stable (by 1.74eV and 0.51eV, respectively) with respect to the hexagonal interstitial, but the Stillinger Weber formation energies are both significantly larger by about 2eV. However, they find that a self-interstitial (which they name an extended interstitial), which has lower symmetry than either of the hexagonal or tetrahedral sites, is an even more stable configuration. This is similar to the ab initio calculations, above, where a low symmetry configuration was also found to be stable. Their calculations suggest a migration mechanism of the extended interstitial where the extra atom moves from this low energy state to an identical one by passing through the state of a tetrahedral interstitial, thus requiring less energy than the tetrahedral-hexagonal transition that was calculated above. The migration paths have not been examined here, but they agree that diffusion of the tetrahedral interstitial via a path though the hexagonal state does not occur until a much higher temperature than the Maroudas-Brown diffusion mechanism.