6.9 Excitation energies

The quasiparticle excitation energies are the energies for adding an electron to the system or subtracting one from it. A quasiparticle energy has both real and imaginary parts, the latter giving the quasiparticle lifetime. For the minimum energy electron and hole quasiparticle excitations the imaginary parts of the quasiparticle energies are zero and the quasiparticles have an infinite lifetime. In this case the exact quasiparticle energy gap can be written as

$$E_{\rm g} = (E_{N+1} - E_{N}) + (E_{N-1} - E_{N})\;,$$ (6.16)

where $E_{N+1}$, $E_{N-1}$, and $E_{N}$ are the ground state total energies of the $N+1$, $N-1$ and $N$ electron systems. [3] A general quasiparticle energy gap cannot be written in terms of differences between energies of exact eigenstates of the system, but such an approximation is often accurate for low energy gaps. The quasiparticle energies are measured in photoemission and inverse photoemission experiments. In an optical absorption experiment a different process occurs in which an electron is excited from the valence to the conduction band. This introduces two quasiparticles into the system, the electron and hole, which interact and can form an exciton, in which case the lowest excitation energy is smaller than $E_{\rm g}$ by the exciton binding energy, $E_{\rm b}$.

The HF method gives approximations to the energies of quasiparticles and the interactions between them. According to Koopmans' theorem, if orbital relaxation is neglected, the quasiparticle energies are equal to the HF eigenvalues. The extended Koopmans' theorem [98,99] has been used in conjunction with VMC methods to calculate quasiparticle energies in silicon (see chapter 7 and reference [100]). In both of these methods the ``quasiparticle energies'' are real as they are obtained as approximations to the energy differences between exact eigenstates of the system.

Recently there has been significant progress in applying QMC techniques to calculate approximate excitation energies from eigenstates using ``direct methods''. In these approaches an excitation energy is obtained by performing separate QMC calculations for the ground and excited states. A Slater-Jastrow wave function is used for the ground state, and for the optical gap the excited state is formed by replacing a valence band single particle orbital by a conduction band one. This calculation is referred to as a ``promotion'' calculation, and such calculations have been reported for a nitrogen solid [42], diamond [101,102], and silicon [103]. Photoemission/inverse photoemission gaps may be obtained by using QMC to calculate the ground state energies of the ${N+1}$, ${N-1}$, and ${N}$ electron systems. Wave functions for the ${N+1}$ and ${N-1}$ electron systems may be formed by adding or subtracting an orbital from the up- or down-spin determinants of the $N$-electron wave function - an ``addition/subtraction'' calculation. For calculations with periodic boundary conditions the simulation cell is made charge neutral by adding a compensating uniform background charge density. Calculations of this type have been reported for one- [104] and two-dimensional [92] model systems, while the results results for a three-dimensional system (silicon) are reported in this thesis.

QMC calculations of excitation energies in extended systems are computationally very demanding because they are `$\frac{1}{N}$' effects, i.e., the fractional change in energy is inversely proportional to the number of electrons in the system. This means that high statistical accuracy is required to obtain good results. The largest system for which excitation energies have been calculated prior to this work is 16 atoms (64 electrons). [103] The total finite size error in the ground state energy for that system was estimated to be about 16 eV per simulation cell, while the energy scale of interest for the excitations is of order 0.1 eV. Like almost all methods for calculating excitation energies, QMC calculations of this type only work because of a strong cancellation of errors between the ground and excited states. Finite size errors tend to reduce the energy gap, while the errors in the trial wave functions are usually larger for the excited states than for the ground state and so increase the energy gap. Although good agreement with experimental excitation energies has been found using small simulation cells, [101,102,103] one is left with the suspicion that if larger simulation cells were used the agreement with experiment might be significantly worse because the finite size effects would be smaller.

In the next sections, the key questions as to the the size and origin of finite size effects in excitation energy calculations and the differences in finite size effects between promotion and addition/subtraction calculations are addressed. Before these QMC techniques can be relied upon for the calculation of excitation energies it is necessary that finite size effects be properly explored, and answers to these questions found.

6.9.1 HF theory of excitation energies

First excitation energies for solids are considered within HF theory. The HF equations with the MPC interaction are obtained by minimizing the HF energy of equation 6.15 with respect to the single-particle orbitals, giving

$$-\frac{1}{2}\nabla^2 \phi_i + \int \rho({\bf r}') v_{\rm E}({\bf r}-{\bf r}') \, {\rm d}{\bf r}' \, \phi_i$$

$$- \sum_{i,j}^N\delta_{s_is_j} \!\! \int \phi_j^*({\bf r}') \phi_j({... ...-{\bf r}') \, {\rm d}{\bf r}' + V_{\rm ext} \phi_i = \epsilon_i \phi_i \;\;\; .$$ (6.17)

If orbital relaxation is neglected, the energy required to excite an electron from the $j^{th}$ (occupied) orbital into the $i^{th}$ (unoccupied) orbital is

$$\Delta E_{ij} = (\epsilon_{i} - \epsilon_{j})$$

$$- \int \rho_i({\bf r}) \rho_j({\bf r}') v_{\rm E}({\bf r}-{\bf r}') \, {\rm d}{\bf r} \, {\rm d}{\bf r}'$$

$$+ \delta_{s_is_j} \int \phi_i^*({\bf r}) \phi_j^*({\bf r}') \phi_j({\bf r}) \phi_i({\bf r}') f({\bf r}-{\bf r}') \, {\rm d}{\bf r} \, {\rm d}{\bf r}'$$

$$+ \frac{1}{2} \int \rho_i({\bf r}) \rho_i({\bf r}') \left[ v_{\rm E... ...bf r}-{\bf r}')-f({\bf r}-{\bf r}')\right] \, {\rm d}{\bf r} \, {\rm d}{\bf r}'$$

$$+ \frac{1}{2} \int \rho_j({\bf r}) \rho_j({\bf r}') \left[ v_{\rm E... ...r}-{\bf r}')-f({\bf r}-{\bf r}')\right] \, {\rm d}{\bf r} \, {\rm d}{\bf r}'\;,$$ (6.18)

where $\rho_k = \vert\phi_k\vert^2$ is the charge density from the $k^{th}$ orbital. The first term is the eigenvalue difference for the excitation while the second and third terms are the Hartree and exchange interactions between the electron and hole. Within this approximation the electron-hole terms go to zero in the limit of an infinitely large simulation cell. The fourth and fifth terms on the right hand side are absent if one uses the Ewald interaction instead of the MPC interaction, i.e., $f$ is replaced by $v_{\rm E}$. When the relaxation of the orbitals is neglected these terms also go to zero when the size of the simulation cell goes to infinity because $v_{\rm E}$ tends to $1/r$ over most of the simulation cell.

The addition/subtraction gap is given by

$$E_{\rm g} = \left(E_{+i}^{\rm HF}-E_0^{\rm HF}\right) - \left(E_0^{\rm HF}-E_{-j}^{\rm HF}\right)$$

$$= (\epsilon_i - \epsilon_j)$$

$$+ \frac{1}{2} \int \rho_i({\bf r}) \rho_i({\bf r}') \left[ v_{\rm E... ...bf r}-{\bf r}')-f({\bf r}-{\bf r}')\right] \, {\rm d}{\bf r} \, {\rm d}{\bf r}'$$

$$+ \frac{1}{2} \int \rho_j({\bf r}) \rho_j({\bf r}') \left[ v_{\rm E... ...{\bf r}')-f({\bf r}-{\bf r}')\right] \, {\rm d}{\bf r} \, {\rm d}{\bf r}' \;\;,$$ (6.19)

where $E_0^{\rm HF}$ is the HF ground state energy of the $N$-electron system, $E_{+i}^{\rm HF}$ is the energy of the state with an electron added to the $i^{th}$ (previously unoccupied) orbital, along with the uniform background charge, and $E_{-j}^{\rm HF}$ is the energy of the state where an electron is removed from the $j^{th}$ orbital, along with the background charge. The standard Koopmans' theorem has been modified and contains two additional terms, which also occur in the promotion energy, $\Delta E_{ij}$. These additional terms, when evaluated using LDA orbitals, are small even for a small simulation cell ($n$ = 2), being in the range $\pm$0.05 eV, and they decrease rapidly with system size. The use of exact HF orbitals or orbital relaxation is not expected to greatly affect these results.

In figure 6.5 the addition/subtraction energies calculated using equation 6.19 for the $\Gamma_{25'} \! \rightarrow \! \Gamma_{15}$ energy gap and the valence band width, calculated with both the Ewald and MPC interactions, are shown along with the LDA values. The results for other low-lying excitation energies show similar behaviour. The promotion energies are not shown figure 6.5 because they differ from the addition/subtraction energies only by the exciton binding energy, which decreases with increasing system size quite rapidly. Figure 6.5 shows that the HF results for the Ewald and MPC interactions are very similar. The band width converges by about $n$ = 7, but the band gap is still slowly increasing at $n$ = 12, and the Ewald and MPC values are not yet equal, which they must be at convergence. For the largest system size studied ($n$ = 12) the MPC gap and valence band width are 7.4 eV and 17.7 eV respectively, which are a little smaller than the HF values of 8.0 eV and 18.9 eV cited in Ref. [105]. The major reasons for these differences are the use of LDA wave functions and LDA-derived pseudopotentials, although as noted above there is clear evidence that in our calculations the HF energy gap has still not fully converged at $n=12$. The LDA excitation energies converge very rapidly with system size.

In a recent study of silicon clusters, Ögüt et al. [106] found large differences between the band gap in the LDA eigenvalues and the band gap calculated by electron addition/subtraction. As shown by Franceschetti et al. [107], these differences are due to the charging of the cluster when an electron is added or subtracted, which does not occur in our calculations because a uniform background is added to preserve charge neutrality. The slow convergence of the HF excitation energies apparent in figure 6.5 therefore arises from the exchange energy. Moreover, because the results with the Ewald and MPC interactions are almost the same, the source of the error is not the interaction with the exchange hole, but the shape of the exchange hole. This is because the excitation energy depends on the change in the exchange hole due to the excitation, which is not strongly localised.

Figure 6.5: Convergence of the $\Gamma_{25'} \! \rightarrow \! \Gamma_{15}$ excitation energy and valence band width of silicon. The data correspond to addition/subtraction energies calculated within HF theory as a function of simulation cell size, $n$, using both the MPC and Ewald interactions. LDA results are also shown.

In a set of calculations on a model two-dimensional system using LDA, $GW$, VMC and DMC techniques, Engel et al. [92] performed a number of VMC calculations with increasing system size and found that the addition/subtraction gap tended to increase with system size. This is the same behaviour as found in the present HF calculations. Engel et al. went on to give an explanation of this effect. Their explanation was that in the $N+1$ ($N-1$) electron systems there is an additional electrostatic energy due to the interaction of the extra electron (hole) with the additional electron (hole) in the other simulation cells. Taking into account the additional compensating uniform background charge that was added to keep each cell neutral, this additional energy is negative and therefore the energies $E_{N+1}$ and $E_{N-1}$ are lower than they should be. Engel et al. showed that the observed finite size effect is much smaller than the Madelung energy for point charges, and to explain this they argued that the effect would be screened by the response of the other electrons. This argument implies that the finite size effects in addition/subtraction calculations are larger than those in promotion calculations.

A reanalysis of the situation follows. HF calculations are considered for simplicity, where the interaction energy can be divided into Hartree and exchange contributions. The significant underestimation of the HF band gaps of small systems is not due to the Hartree terms, which by construction are the same as for our LDA calculations and give very small finite size effects in the band gaps. The finite size error in the HF gaps therefore arises from the exchange energy. By comparing band gaps calculated with the Ewald and MPC interactions we can see whether the problem lies with the interaction or with the shape of the exchange hole. Because the band gaps calculated with the Ewald and MPC interactions are very similar, the form of the interaction is not the important consideration. Therefore the source of the problem must be the finite size errors present in the shape of the exchange hole. This argument implies that the finite size effects in addition/subtraction calculations are similar to those in promotion calculations. This view is supported by the HF results that have been presented in this section and also by the VMC results of the next section.

In summary, the HF excitation energies calculated with the Ewald and MPC interactions are very similar. Within HF theory the largest finite size error in excitation energies arises from the shape of the exchange hole, which leads to slow convergence with system size. The finite size errors in promotion and addition/subtraction HF calculations are of similar size.

6.9.2 QMC theory of excitation energies

We now apply the theory developed in the previous section to QMC calculations of excitation energies. Although the HF gaps converge rather slowly with system size, it was shown in section 6.8.2 and section 6.8.3 that the finite size effects in VMC and DMC ground state energies are smaller than in HF theory. It is important to see whether this also applies to excitation energies.

VMC is computationally cheaper than DMC and so excitation energies may be computed using VMC over a larger range of system sizes. It is expected that the finite size effects in DMC will follow those in VMC, as the VMC calculations retrieve about 90% of the fixed-node correlation energy. The $\Gamma_{25'} \! \rightarrow \! \Gamma_{15}$ excitation energies were computed within VMC in silicon for the system sizes $n$ = 1, 2, 3, 4 using both promotion and addition/subtraction methods. The calculations were performed with ${\bf k}_s \equiv \Gamma$ and the other computational details were the same as for the ground state calculations. The Jastrow factors were left unchanged for the excited state. In tests on the $n$ = 2 system, separately optimizing the Jastrow factors for both ground and excited states did not significantly change the results. The computational cost of the $n$ = 4 calculations is very large; an error bar in the excitation energies of $\pm 0.3$ eV requires an error bar of $\pm 0.0006$ eV per electron. Although the computational effort is large, such a study is necessary to establish the accuracy of QMC excitation energy calculations.

In figure 6.6 excitation energies obtained with the Ewald interaction via the promotion and addition/subtraction methods are shown. Results for the MPC interaction were found to be very similar. The promotion and addition/subtraction results are nearly the same, but the promotion energies are slightly smaller because they include an exciton binding energy, which decreases as the system size increases. The results are consistent with a slow increase in the excitation energy with system size and indicate that reasonable convergence is already obtained at $n$ = 2. The increase in excitation energies with system size is the same trend as in HF calculations, although the finite size errors are smaller in the correlated calculations. The finite size errors in the promotion and addition/subtraction methods are not significantly different at this level of statistical accuracy. On general grounds one would expect the finite size effect in the promotion calculation to be slightly larger. This follows because the trend for both promotion and addition/subtraction calculations is that the excitation energy is reduced for small system sizes and this effect is enhanced in the promotion calculation by the exciton binding energy which is larger for small systems. The calculated excitation energy is roughly 4 eV, which is larger than the experimental value of 3.40 eV, [108] and also a little larger than the DMC value for an $n$ = 2 simulation cell of 3.7 eV. [103] This study demonstrates that the largest contribution to the error in the VMC band gap for $n \geq 2$ arises from the approximate nature of the trial wave functions, and not from finite size effects.

Figure 6.6: The $\Gamma_{25'} \! \rightarrow \! \Gamma_{15}$ excitation energy of silicon calculated within VMC theory as a function of simulation cell size, $n$. Data for the Ewald interaction are shown obtained via both the promotion and addition/subtraction methods.

The exciton binding energy can be calculated as the difference between the promotion and addition/subtraction gaps. The exciton binding energy for the $\Gamma_{25'} \! \rightarrow \! \Gamma_{15}$ excitation is small and it is only statistically resolvable for the smallest ($n=1$) cell, which gives a value of $0.28 \pm 0.01$ eV. In earlier QMC calculations [42,101,102] the exciton binding energy was estimated using the Mott-Wannier formula

$$E_{\rm b} \sim \frac{1}{2\epsilon r}\;\;,$$ (6.20)

where $\epsilon$ is the relative permittivity and $r$ is the radius of the localisation region. Using $\epsilon$ = 11.7 and the appropriate radius for $n$ = 1 of $r$ = 4.0 a.u. gives an exciton binding energy of 0.29 eV. This is extremely close to the VMC value, but the excellent agreement is probably fortuitous since the $n$ = 1 cell is so small that it is appropriate to use a value of $\epsilon$ at finite wave vector, which would be smaller. The exciton binding energy may also be evaluated within HF theory as the sum of the second and third terms on the r.h.s. of equation 6.18. This gives 0.75 eV for the $n=1$ cell using the Ewald interaction, which is considerably larger than the VMC value because the latter calculation includes screening of the electron-electron interaction.

Note that the promotion and addition/subtraction methods differ significantly in the required computational effort. If the intrinsic variance of the local energy is the same for each of the energies, which is a good approximation for silicon wavefunctions, then the most efficient way to achieve the same error bar in an addition/subtraction gap is to perform $2M$ moves for each of the $N+1$ and $N-1$ systems and $4M$ for the ground state, giving a total cost of $8M$ moves, where an acceptable error bar is obtained in a promotion calculation by performing $M$ Monte Carlo moves for both the ground and excited states. It is therefore four times more expensive to calculate an energy gap to some given accuracy by the addition/subtraction method than by the promotion method.

6.9.3 Modified interaction for excitation energies

In this section a modified electron-electron interaction specifically designed to describe excitation energies within periodic boundary conditions simulations in introduced. Two problems arise when trying to model excitations using finite simulation cells subject to periodic boundary conditions. One is that the excitation is ``squeezed'' into the simulation cell, and the other is that there are spurious interaction between the periodic replicas of the simulation cell. Here the problem of the spurious interactions is addressed using the ideas of the MPC interaction. The charge density on promotion or addition/subtraction of an electron can be written as the sum of the ground state charge density and a small deviation, i.e., $\tilde{\rho}({\bf r}) = \rho({\bf r}) + \Delta({\bf r})$. The Hartree term can be modified so that this charge density interacts with the ground state density in the replicas. This leads to the interaction energy

$$\tilde{E}_{\rm e-e} = \int \vert\tilde{\phi}\vert^2\sum_{i>j} f({\bf r}_i-{\bf r}_j) \prod_k {\rm d}{\bf r}_k$$

$$+ \int \tilde{\rho}({\bf r})\rho({\bf r}') \left[ v_{\rm E}({\bf r}-{\bf r}')- f({\bf r}-{\bf r}') \right] \, {\rm d}{\bf r}\, {\rm d}{\bf r}'$$

$$- \frac{1}{2}\int \rho({\bf r})\rho({\bf r}') \left[ v_{\rm E}({\bf... ...{\bf r}')-f({\bf r}-{\bf r}') \right] \, {\rm d}{\bf r}\, {\rm d}{\bf r}' \;\;,$$ (6.21)

where $\vert\tilde{\phi}\vert^2$ generates the charge density $\tilde{\rho}$ and the ground state charge density, $\rho$, is fixed. A HF analysis of this interaction shows that the HF equations are identical to equation 6.17, so that the orbitals and eigenvalues are unaltered. However, the ground and excited state energy expressions are modified. For the excited states analogues of equations 6.18 and 6.19 are obtained, but without the terms involving $(v_{\rm E}-f)$, i.e., the standard Koopmans' theorem is retrieved. Although these terms are small for silicon, they will be significant in cases when the change in the charge density of the exciton is strongly localized. This analysis also provides further evidence that the electrostatic interactions between the simulation cell and its replicas is not necessarily an important source of finite size error in excited state energy calculations; for these silicon systems, the finite size errors in the excitations are small.