# Fully parameter-free calculation of optical spectra for insulators, semiconductors and metals from a simple polarization functional

###### Abstract

We present a fully parameter-free density-functional approach for the accurate description of optical absorption spectra of insulators, semiconductors and metals. We show that this can be achieved within time-dependent current-density-functional theory using a simple dynamical polarization functional. We derive this functional from physical principles that govern optical spectra. Our method is truly predictive because not a single parameter is used. In particular, we do not use an ad-hoc material-dependent broadening parameter to compare theory to experiment as is usually done. Our approach is numerically efficient; the cost equals that of a calculation within the random-phase approximation.

###### pacs:

71.15.Mb, 71.45Gm, 71.10-w, 71.15QeThe calculation of accurate optical absorption spectra within a density-functional approach has been a major challenge for almost two decades. While an accurate approach to calculate optical spectra is by solving the Bethe Salpeter equation (BSE) Albrecht et al. (1998); Benedict et al. (1998); Rohlfing and Louie (1998); Onida et al. (2002), the large numerical effort required precludes its application to large systems. Therefore an approach based on the numerically more affordable time-dependent density-functional theory (TDDFT) Runge and Gross (1984) or time-dependent current-density-functional theory (TDCDFT) Dhara and Ghosh (1987); Ghosh and Dhara (1988); Vignale (2004) that gives optical spectra of similar quality as the BSE is much sought after. The price to pay for the gained computational efficiency is that the auxiliary Kohn-Sham system of TD(C)DFT is governed by an unknown effective potential for which the derivation of accurate approximations is highly nontrivial. The main failures of standard approximations such as the random-phase approximation (RPA) and the adiabatic local-density approximation (ALDA) Zangwill and Soven (1980) are: 1) the underestimation of the onset of the absorption; 2) the underestimation of the intensity of continuum excitons; 3) the absence of bound excitons; 4) the absence of Drude tails in the spectra of metals. While the first problem can be circumvented by replacing the Kohn-Sham eigenvalues by quasiparticle energies obtained within the approximation Hedin (1965); Aulbur et al. (2000), the other failures have proven to be more complicated to solve. Early attempts to go beyond such standard approximations focused on a correct description of continuum excitons by employing a long-range exchange-correlation (xc) kernel either within TDCDFT de Boeij et al. (2001) or TDDFT Reining et al. (2002); Botti et al. (2004). These approaches produced the desired result but at the cost of introducing a material-dependent parameter. The nanoquanta kernel Sottile et al. (2003a); Marini et al. (2003); Adragna et al. (2003); Stubner et al. (2004); von Barth et al. (2005) does not require such a parameter and leads to a correct description of both bound and continuum excitons. Unfortunately, being derived from the BSE, it is almost as expensive to evaluate as solving the full BSE. An efficient and parameter-free TDCDFT approach can be obtained using the Vignale-Kohn functional Vignale and Kohn (1996). While this method describes well the optical spectra of metals Berger et al. (2006) it does not describe correctly excitons in semiconductors and insulators Berger et al. (2007). Recently several new approaches have been proposed. The bootstrap method advances an expression for the TDDFT xc kernel which has to be obtained from a self-consistent procedure Sharma et al. (2011). In Ref. [Nazarov and Vignale, 2011] an xc functional is derived from a meta-generalized-gradient approximation. Finally, Ref. [Trevisanutto et al., 2013] proposes an xc functional that is based on the jellium-with-gap model. Although these new methods improve the description of excitonic effects, they have several shortcomings: 1) they are not parameter free due to an ad-hoc material-dependent broadening parameter that is used to compare theory and experiment; 2) they are static and therefore do not account for memory effects precluding a parameter-free description of the finite width of Drude tails and bound excitons; 3) they require the calculation of the full Kohn-Sham density-density response function, an calculation (with the number of atoms or electrons in the unit cell).

In this work we present an efficient density functional approach that does not suffer from the shortcomings mentioned above. It is parameter free and dynamic, and it produces accurate optical spectra for metals, semiconductors, and insulators. The method we propose is both simple and elegant and is characterized by the following two steps: 1) We calculate the macroscopic dielectric function using an efficient TDCDFT approach Kootstra et al. (2000); Romaniello and de Boeij (2005) within the RPA (including a scissors shift obtained from ); 2) We apply a simple dynamical polarization functional (see Eq. (14) below) that accounts for continuum and bound excitons as well as Drude tails. We will now give all the details of our approach. We use Hartree atomic units throughout.

The optical response of a solid can be obtained as the response to a homogeneous electric field, i.e., a field with vanishing wave vector . However, for , the periodic density is not enough to uniquely describe the response; one also needs the macroscopic polarization Gonze et al. (1995). This is an immediate consequence of the use of periodic boundary conditions. Macroscopic effects due to the surface density have to be accounted for via the macroscopic polarization. An elegant and efficient approach is to use the current density as the fundamental quantity of a DFT for time-dependent phenomena. The self-consistent perturbing potential is a unique functional of the periodic current density since the latter contains both the periodic density of the bulk (via the continuity relation) and the macroscopic polarization, which is obtained as a bulk property using the definition,

(1) |

Here is the volume of the unit cell. The limit can thus be performed explicitly, by setting , and hence knowledge of the current density suffices and no explicit calculation of Kohn-Sham response functions is needed. The macroscopic polarization is induced by a macroscopic electric field which comprises both the externally applied field and the macroscopic induced electric field. The constant of proportionality is the electric susceptibility tensor which is defined as

(2) |

The macroscopic dielectric tensor can be obtained from the electric susceptibility according to

(3) |

Therefore, for a given , knowledge of the induced current suffices to calculate from Eqs. (1)-(3).

Within Kohn-Sham linear response the current density is given by

(4) |

where . We note that the Kohn-Sham current-current and current-density response functions ( and , respectively) are never needed explicitly since we use a sum-over-states representation. Therefore the integrals can be evaluated independently of . In Eq. (4) we used the microscopic Coulomb gauge in which the microscopic potential is lattice periodic and contains the periodic part of the Hartree and longitudinal xc contributions Kootstra et al. (2000). All remaining xc contributions are included in the the vector potential . Here we neglect microscopic transverse xc contributions and therefore only its macroscopic part, , enters Eq. (4). It is related to the induced current through the TDCDFT tensor xc kernel :

(5) |

Substitution of Eq. (4) into Eq. (1), and using , reveals that is linear in the macroscopic Kohn-Sham electric field, . Following Ref. [de Boeij et al., 2001], we define the auxiliary susceptibility according to

(6) |

where is the susceptibility obtained from a calculation with . We note that the contributions of to are fully accounted for. Since is itself a functional of this is done within a self-consistent field (SCF) calculation. We will now show that the macroscopic xc effects can be accounted for post-SCF. Comparison of Eqs. (2) and (6) leads to the following relation between and ,

(7) |

If one chooses equal to zero as is usually done we obtain . We can now go beyond this approximation and obtain a polarization functional for by neglecting microscopic current components in Eq. (5), i.e., we replace by its unit-cell average. We obtain

(8) |

where we used Eq. (1) and in which we defined

(9) |

Substitution of Eq. (8) into Eq. (7) leads to

(10) |

Therefore, for a given , we can simply calculate from .

Here we will use the RPA () to calculate , i.e., . Since continuum excitons are underestimated and bound excitons and Drude tails are absent in RPA optical spectra, we will include these effects through . We first derive a static from the physical principles that govern strongly bound excitons in solids. To simplify our argumentation we will assume for the moment that , and are isotropic, i.e., , etc.. Following arguments similar to those given in Ref. [Sottile et al., 2003b], we derive in the supplemental material sup () that the exact constraint to have a bound exciton is that for some below the band gap the following relation holds for ,

(11) |

However, the application of the above expression requires the knowledge of which is unknown. Therefore, we would like to rewrite the above expression in terms of known quantities. In the supplemental material sup () we derive such an expression for the case of a wide-gap insulator whose optical spectra is dominated by a bound exciton. The result is a static given by

(12) |

Although the derivation applies to wide-gap insulators, Eq. (12) has the additional advantage that is proportional to . It has previously been shown numerically that such a proportionality leads to good optical spectra for semiconductors Botti et al. (2004, 2005). The physical reason is that is a measure of the amount of screening of the electron-hole interaction. Screening effects are more important in semiconductors than in wide-gap insulators. We note that Eq. (12) has a similar form as the bootstrap kernel Sharma et al. (2011) but that no self-consistent procedure is needed to calculate it. Unfortunately in Eq. (12) is static and will therefore not be able to account for the finite width of bound excitons and Drude tails. For this reason we add to Eq. (12) , the long-range part of the dynamical Vignale-Kohn functional Vignale and Kohn (1996) in which we replace the current density by its unit-cell average de Boeij et al. (2001):

(13) |

Here is the longitudinal (transverse) xc kernel of the homogeneous electron gas, is the xc energy per volume of the homogeneous electron gas, is the ground-state density and is its average in the unit cell. For we use the parametrization of Ref. [Qian and Vignale, 2002] in the QVA approximation Berger et al. (2006). The Vignale-Kohn functional is the exact functional of a slightly inhomogeneous electron gas and describes correctly the optical spectra of metals Berger et al. (2006). For this reason Eq. (13) is complementary to Eq. (12) which tends to zero for metallic systems because screening is complete and therefore . will account for the Drude tails and finite width of bound excitons which are exactly the features that are absent in Eq. (12). Moreover, as we show in the supplemental material sup (), will not much influence the spectra of semiconductors. We finally obtain the following approximation for ,

(14) |

where we generalized Eq. (12) to a tensor form. Since and commute, the order of the multiplication is irrelevant. We note that Eq. (14) satisfies the Kramers-Kronig relations. Equation (14) is the main result of this work.

We will now demonstrate our approach by applying it to the calculation of the optical spectra of several materials. We briefly outline the full procedure. The ground-state calculations are done within the local-density approximation (LDA) Kohn and Sham (1965) and we use LDA lattice parameters in all calculations. We apply a scissors operator to shift the unoccupied bands and we modify the current operator accordingly to guarantee that exact constraints such as the continuity equation remain satisfied. The energy shift is calculated with the method Hedin (1965); Aulbur et al. (2000) and is equal to the correction for the direct band gap at the point. We calculated this correction using the effective-energy technique Berger et al. (2010, 2012a, 2012b), which leads to a speed-up of an order of magnitude. The calculations were done with the Abinit code Gonze et al. (2005). We implemented the polarization functional in a modified version of the Amsterdam Density Functional (ADF) code te Velde et al. (2001); Fonseca Guerra et al. (1998); ADF (). We use the TZ2P (triple- + 2 polarization functions) basis set provided by ADF. The -space integrals are done analytically using a Lehmann-Taut tetrahedron scheme Lehmann and Taut (1972). Since we do not include effects due to electron-phonon coupling the spectra obtained with the above approach are predictions of the optical spectra at low temperature where electron-electron scattering dominates electron-phonon scattering. For this reason we will compare our calculated spectra with spectra measured at low temperature where available.

In Fig. 1 we report the optical absorption spectra of bulk silicon and bulk GaP obtained with our polarization functional and compare it to the RPA spectra and to experimental results obtained at low temperature (15 Kelvin). Silicon and GaP are typical examples of materials for which the RPA strongly underestimates the first peak which appears in the experimental spectra around 3.4 eV (Si) and 3.8 eV (GaP). Our polarization functional solves this problem by including the necessary excitonic effects and the first peak compares well with experiment both in position and magnitude. Overall, the spectra are very close to experiment with the exception of the peak around 5.2 eV in the spectrum of GaP which is overestimated.

In Fig. 2 we show the optical absorption spectra of solid argon and LiF obtained with our polarization functional and compare it to the RPA spectra and to experimental results. Solid argon and LiF are typical materials that exhibit strongly bound excitons. We see that these excitons which appear in the experimental spectra around 12 eV (Ar) and 12.5 eV (LiF) are completely absent in the RPA spectra. Our polarization functional describes these bound excitons and also accurately reproduces their position. The magnitude of the peaks is overestimated with respect to experiment. This discrepancy is probably due to electron-phonon broadening in the experiments and to the fact that density-functional approaches tend to overestimate these peaks Sharma et al. (2011); Trevisanutto et al. (2013). Finally, we verify a posteriori that given in Eq. (12) is indeed a good approximation to Eq. (11) for wide-gap insulators exhibiting a dominant bound exciton. With Eq. (12) we obtain and , while Eq. (11) gives and .

In Fig. 3 we report the optical absorption spectra of diamond and copper obtained with the polarization functional and compare it to the RPA spectra and to experimental results. Diamond is another typical test case since the RPA spectrum is quite different from the experimental spectrum. Due to the absence of excitonic effects the RPA spectrum has too much weight at high energy. With our polarization functional the spectral weight is shifted to lower energy and we obtain a very good agreement with experiment. While the RPA spectrum of copper accurately reproduces the part of the spectrum which is due to interband transitions, the Drude tail at low energy, which is due to intraband transitions, is completely absent. Our polarization functional accurately describes the Drude tail while maintaining the good agreement for the interband part.

In conclusion, we presented the first fully parameter-free density-functional approach that gives accurate optical spectra for insulators, semiconductors and metals alike. Our approach is therefore truly predictive and due to its numerical efficiency opens the way for the prediction of optical spectra of large systems.

We thank Paul L. de Boeij, Lucia Reining, and Pina Romaniello for fruitful discussions.

## References

- Albrecht et al. (1998) S. Albrecht, L. Reining, R. Del Sole, and G. Onida, Phys. Rev. Lett. 80, 4510 (1998).
- Benedict et al. (1998) L. X. Benedict, E. L. Shirley, and R. B. Bohn, Phys. Rev. B 57, R9385 (1998).
- Rohlfing and Louie (1998) M. Rohlfing and S. G. Louie, Phys. Rev. Lett. 81, 2312 (1998).
- Onida et al. (2002) G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 74, 601 (2002).
- Runge and Gross (1984) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
- Dhara and Ghosh (1987) A. K. Dhara and S. K. Ghosh, Phys. Rev. A 35, 442 (1987).
- Ghosh and Dhara (1988) S. K. Ghosh and A. K. Dhara, Phys. Rev. A 38, 1149 (1988).
- Vignale (2004) G. Vignale, Phys. Rev. B 70, 201102 (2004).
- Zangwill and Soven (1980) A. Zangwill and P. Soven, Phys. Rev. A 21, 1561 (1980).
- Hedin (1965) L. Hedin, Phys. Rev. 139, A796 (1965).
- Aulbur et al. (2000) W. G. Aulbur, L. Jönsson, and J. W. Wilkins, in Solid State Physics (Academic, New York, 2000), vol. 54, p. 1.
- de Boeij et al. (2001) P. L. de Boeij, F. Kootstra, J. A. Berger, R. van Leeuwen, and J. G. Snijders, J. Chem. Phys. 115 (2001).
- Reining et al. (2002) L. Reining, V. Olevano, A. Rubio, and G. Onida, Phys. Rev. Lett. 88, 066404 (2002).
- Botti et al. (2004) S. Botti, F. Sottile, N. Vast, V. Olevano, L. Reining, H.-C. Weissker, A. Rubio, G. Onida, R. Del Sole, and R. W. Godby, Phys. Rev. B 69, 155112 (2004).
- Sottile et al. (2003a) F. Sottile, V. Olevano, and L. Reining, Phys. Rev. Lett. 91, 056402 (2003a).
- Marini et al. (2003) A. Marini, R. Del Sole, and A. Rubio, Phys. Rev. Lett. 91, 256402 (2003).
- Adragna et al. (2003) G. Adragna, R. Del Sole, and A. Marini, Phys. Rev. B 68, 165108 (2003).
- Stubner et al. (2004) R. Stubner, I. V. Tokatly, and O. Pankratov, Phys. Rev. B 70, 245119 (2004).
- von Barth et al. (2005) U. von Barth, N. E. Dahlen, R. van Leeuwen, and G. Stefanucci, Phys. Rev. B 72, 235109 (2005).
- Vignale and Kohn (1996) G. Vignale and W. Kohn, Phys. Rev. Lett. 77, 2037 (1996).
- Berger et al. (2006) J. A. Berger, P. Romaniello, R. van Leeuwen, and P. L. de Boeij, Phys. Rev. B 74, 245117 (2006).
- Berger et al. (2007) J. A. Berger, P. L. de Boeij, and R. van Leeuwen, Phys. Rev. B 75, 035116 (2007).
- Sharma et al. (2011) S. Sharma, J. K. Dewhurst, A. Sanna, and E. K. U. Gross, Phys. Rev. Lett. 107, 186401 (2011).
- Nazarov and Vignale (2011) V. U. Nazarov and G. Vignale, Phys. Rev. Lett. 107, 216402 (2011).
- Trevisanutto et al. (2013) P. E. Trevisanutto, A. Terentjevs, L. A. Constantin, V. Olevano, and F. Della Sala, Phys. Rev. B 87, 205143 (2013).
- Kootstra et al. (2000) F. Kootstra, P. L. de Boeij, and J. G. Snijders, J. Chem. Phys. 112, 6517 (2000).
- Romaniello and de Boeij (2005) P. Romaniello and P. L. de Boeij, Phys. Rev. B 71, 155108 (2005).
- Gonze et al. (1995) X. Gonze, P. Ghosez, and R. W. Godby, Phys. Rev. Lett. 74, 4035 (1995).
- Sottile et al. (2003b) F. Sottile, K. Karlsson, L. Reining, and F. Aryasetiawan, Phys. Rev. B 68, 205112 (2003b).
- (30) See Supplemental Material [URL inserted by publisher] which includes Refs. Sottile et al. (2003b); Botti et al. (2005); Sottile et al. (2005).
- Botti et al. (2005) S. Botti, A. Fourreau, F. Nguyen, Y.-O. Renault, F. Sottile, and L. Reining, Phys. Rev. B 72, 125203 (2005).
- Qian and Vignale (2002) Z. Qian and G. Vignale, Phys. Rev. B 65, 235121 (2002).
- Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
- Berger et al. (2010) J. A. Berger, L. Reining, and F. Sottile, Phys. Rev. B 82, 041103(R) (2010).
- Berger et al. (2012a) J. A. Berger, L. Reining, and F. Sottile, Phys. Rev. B 85, 085126 (2012a).
- Berger et al. (2012b) J. A. Berger, L. Reining, and F. Sottile, Eur. Phys. J. B 85, 326 (2012b).
- Gonze et al. (2005) X. Gonze, G. Rignanese, M. Verstraete, J. Beuken, Y. Pouillon, R. Caracas, F. Jollet, M. Torrent, G. Zerah, M. Mikami, et al., Z. Kristallogr. 220, 558 (2005).
- te Velde et al. (2001) G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders, and T. Ziegler, J. Comp. Chem. 22, 931 (2001).
- Fonseca Guerra et al. (1998) C. Fonseca Guerra, J. G. Snijders, G. te Velde, and E. J. Baerends, Theor. Chem. Acc. 99, 391 (1998).
- (40) ADF2013 (modified version), SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, www.scm.com.
- Lehmann and Taut (1972) G. Lehmann and M. Taut, Phys. Status Solidi B 54, 469 (1972).
- Lautenschlager et al. (1987) P. Lautenschlager, M. Garriga, L. Vina, and M. Cardona, Phys. Rev. B 36, 4821 (1987).
- Zollner et al. (1993) S. Zollner, M. Garriga, J. Kircher, J. Humlíček, M. Cardona, and G. Neuhold, Phys. Rev. B 48, 7915 (1993).
- Saile et al. (1976) V. Saile, M. Skibowski, W. Steinmann, P. Gürtler, E. E. Koch, and A. Kozevnikov, Phys. Rev. Lett. 37, 305 (1976).
- Roessler and Walker (1967) D. M. Roessler and W. C. Walker, J. Opt. Soc. Am. 57, 835 (1967).
- Phillip and Taft (1964) H. R. Phillip and E. A. Taft, Phys. Rev. 136, A1445 (1964).
- Stahrenberg et al. (2001) K. Stahrenberg, T. Herrmann, K. Wilmers, N. Esser, W. Richter, and M. J. G. Lee, Phys. Rev. B 64, 115111 (2001).
- Hagemann et al. (1975) H. J. Hagemann, W. Gudat, and C. Kunz, J. Opt. Soc. Am. 65, 742 (1975).
- Sottile et al. (2005) F. Sottile, F. Bruneval, A. G. Marinopoulos, L. K. Dash, S. Botti, V. Olevano, N. Vast, A. Rubio, and L. Reining, Int. J. Quantum Chem. 102, 684 (2005).