Oneshot calculation of temperaturedependent optical spectra
and
phononinduced bandgap renormalization
Abstract
Recently, Zacharias et al. [Phys. Rev. Lett. 115, 177401 (2015)] developed a new ab initio theory of temperaturedependent optical absorption spectra and band gaps in semiconductors and insulators. In that work the zeropoint renormalization and the temperature dependence were obtained by sampling the nuclear wavefunctions using a stochastic approach. In the present work, we show that the stochastic sampling of Zacharias et al. can be replaced by fully deterministic supercell calculations based on a single optimal configuration of the atomic positions. We demonstrate that a single calculation is able to capture the temperaturedependent band gap renormalization including quantum nuclear effects in directgap and indirectgap semiconductors, as well as phononassisted optical absorption in indirectgap semiconductors. In order to demonstrate this methodology we calculate from first principles the temperaturedependent optical absorption spectra and the renormalization of direct and indirect band gaps in silicon, diamond, and gallium arsenide, and we obtain good agreement with experiment and with previous calculations. In this work we also establish the formal connection between the WilliamsLax theory of optical transitions and the related theories of indirect absorption by Hall, Bardeen, and Blatt, and of temperaturedependent band structures by Allen and Heine. The present methodology enables systematic ab initio calculations of optical absorption spectra at finite temperature, including both direct and indirect transitions. This feature will be useful for highthroughput calculations of optical properties at finite temperature, and for calculating temperaturedependent optical properties using highlevel theories such as GW and BetheSalpeter approaches.
pacs:
78.40.q, 71.15.m, 71.38.kI Introduction
The electronphonon interaction plays a central role in the optical properties of solids. For example, electronphonon couplings lead to the temperature dependence and the quantum zeropoint renormalization of the critical point energies, to temperaturedependent broadening of light absorption and emission lineshapes, and to indirect optical transitions.
Recently, it has become possible to study these effects using ab initio calculations. The phononinduced renormalization of band gaps and band structures was investigated from first principles in Refs. Allen and Cardona, 1981; Marini, 2008; Giustino et al., 2010; Gonze et al., 2011; Cannuccia and Marini, 2011; Poncé et al., 2014a, b; Antonius et al., 2014; Monserrat and Needs, 2014; Poncé et al., 2015; LloydWilliams and Monserrat, 2015; Monserrat, 2016a; Nery and Allen, starting from the theory of Allen and Heine (AH).Allen and Heine (1976) The optical absorption in indirectgap semiconductors was studied from first principles in Ref. Noffsinger et al., 2012 using the classic theory of Hall, Bardeen, and Blatt (HBB),Hall et al. (1954) and in Ref. Zacharias et al., 2015 using the theory of WilliamsWilliams (1951) and LaxLax (1952) (WL). A review of the standard formalism and the of most recent literature can be found in Ref. Giustino, .
In this manuscript we focus on the WL theory, and on how to perform accurate and efficient ab initio calculations of temperaturedependent band gaps and optical spectra in semiconductors and insulators using the WL formalism. In its original formulation the WL theory was employed to study the vibrational broadening of the photoluminescence spectra of defects in solids.Williams (1951); Lax (1952) In a recent work we showed that the same theory can successfully be employed for predicting temperature dependent optical spectra and band gaps in semiconductors, including phononassisted indirect absorption.Zacharias et al. (2015) The reason for this success is that a perturbative treatment of the WL theory naturally leads to the adiabatic approximations of the AH and the HBB theories. In fact, as it was shown in Ref. Zacharias et al., 2015, the AH theory of temperaturedependent band structures can alternatively be derived from the WL theory by neglecting the optical matrix elements. Similarly, the adiabatic limit of the HBB theory of indirect optical absorption can be derived from the WL theory by retaining only onephonon processes. The relations between the WL, the HBB, and the AH theory will be analyzed in detail in Sec. V.
In the WL theory,Williams (1951); Lax (1952); Patrick and Giustino (2013, 2014); Zacharias et al. (2015) the effect of quantum nuclear motion on the optical properties is described by first calculating the optical spectrum evaluated at clamped nuclei, and then taking the expectation value of this quantity over a given nuclear wavefunction. The temperature is introduced by performing a canonical average over all possible nuclear quantum states. Formally, this approach corresponds to a ‘semiclassical’ FranckCondon approximation, whereby the initial quantum states are described by BornOppenheimer products of electronic and nuclear wavefunctions, and the final quantum states are replaced by a classical continuum. This approach is related to but does not coincide with the standard adiabatic approximation. An extensive discussion of the formalism and its limit of applicability can be found in Refs. Lax, 1952; Patrick and Giustino, 2014.
The key advantages of the WL method are (i) the calculations are simple and can be performed as a postprocessing step on top of any electronic structure code. (ii) The formalism is agnostic of the level of theory used to describe optical excitations at clamped nuclei; therefore the same procedure can be used with any level of theory (e.g. independentparticle approximation, randomphase approximation, GW/BetheSalpeter), so long as the optical process can be described by means of Fermi’s Golden Rule. (iii) The method seamlessly combines the AH theory of temperaturedependent band structure and the HBB theory, of phononassisted indirect optical absorption.
The two main disadvantages of the WL method are (i) the calculations require the use of supercells in order to accommodate phonon wavevectors within the first Brillouin zone. (ii) The evaluation of expectation values over the nuclear wavefunctions requires calculations for many nuclear configurations. In Refs. Patrick and Giustino, 2013; Zacharias et al., 2015 the latter issue was addressed by using a stochastic approach based on importancesampling Monte Carlo integration. In this manuscript we further improve the configurational averaging by replacing the stochastic approach of Ref. Zacharias et al., 2015 with a fully deterministic method. In particular, we demonstrate that it is possible to choose a single configuration of the nuclei yielding at once the band structure renormalization and indirect optical absorption at a given temperature. In order to demonstrate this method we report applications to silicon, diamond and gallium arsenide. Our calculated spectra and temperature dependent band gaps compare well with previous calculations and with experiment. For completeness we also provide a detailed analysis of the relation between the WL theory, the AH theory, and the HBB theory.
The organization of the manuscript is as follows: in Sec. II we briefly outline the WL expression for the temperaturedependent dielectric function, and summarize the ‘oneshot’ procedure for evaluating this expression using a single atomic configuration. In this section we also show our main results for the optical absorption spectra of Si, C, and GaAs in order to emphasize the simplicity and effectiveness of the formalism. In Sec. III we develop the formalism which is used to select the optimal atomic configuration in the oneshot calculations of Sec. II. In particular, we prove that our optimal configuration yields exact results in the limit of infinite supercell size. In Sec. IV we extend the concepts of Sec. III by showing that it is possible to deterministically select further atomic configurations in order to control and systematically reduce the error resulting from the configurational averaging. In Sec. V we discuss the link between the WL theory of temperaturedependent optical spectra, the AH theory of temperaturedependent band structures, and the HBB theory of indirect optical absorption. In Sec. VI we present our calculations of temperaturedependent band gaps for silicon, diamond, and gallium arsenide. Section VII reports all computational details of the calculations presented in this manuscript. In Sec. VIII we summarize our key findings and indicate avenues for future work. Lengthy formal derivations and further technical details are left to Appendices AD.
Ii Oneshot method and main results
In this section we outline the procedure for calculating temperaturedependent optical spectra using oneshot frozenphonon calculations. For clarity we also anticipate our main results on silicon, diamond, and gallium arsenide, leaving all computational details to Sec. VII.
In the WL theory the imaginary part of the dielectric function of a solid at the temperature is given by:Zacharias et al. (2015)
(1) 
In this expression, denotes the energy of a nuclear quantum state evaluated in the BornOppenheimer approximation, is the Boltzmann constant, and is the canonical partition function. The function is the imaginary part of the macroscopic, electronic dielectric function, evaluated at clamped nuclei. For notational simplicity we indicate the set of all atomic coordinates by . In the following we denote by the total number of atomic coordinates. In Eq. (1) each expectation value is taken with respect to the quantum nuclear state with energy , and involves a multidimensional integration over all atomic coordinates. A detailed derivation of Eq. (1) can be found in Sec. 9.2 of Ref. Patrick and Giustino, 2014.
In order to focus on quantum nuclear effects and temperature shifts, we here describe the dielectric function at clamped nuclei using the simplest possible approximations, namely the independentparticle approximation and the electric dipole approximation:
(2) 
In this expression is the electron mass, is the number of electrons in the crystal unit cell, is the plasma frequency, and the photon frequency. The factor 2 is for the spin degeneracy. The sum extends to the occupied KohnSham states of energy , as well as the unoccupied states of energy . The superscripts are to keep in mind that these states are evaluated for nuclei clamped in the configuration labelled by . The matrix elements of the momentum operator along the polarization direction of the photon is indicated as . In the present case we use nonlocal pseudopotentials and a scissor operator, therefore the momentum matrix elements are modified following Ref. Starace, 1971, as described in Sec. VII. In all the calculations presented in this manuscript the dielectric functions are obtained by first evaluating Eqs. (1) and (2) for each Cartesian direction, and then performing the isotropic average over the photon polarizations.
In principle Eq. (1) could be evaluated using the nuclear wavefunctions obtained from the solution of the nuclear Schrödinger equation with electrons in their ground state. This choice would lead to the automatic inclusion of anharmonic effects. However, for conciseness, in the present work we restrict the discussion to the harmonic approximation.
In the harmonic approximation, every manybody nuclear quantum state can be expressed as a product of Hermite functions, and the atomic displacements can be written as linear combinations of normal coordinates. Sakurai (1994) By exploiting the property of Hermite polynomials and Mehler’s formula, Watson (1933) the summation in Eq. (1) is exactly rewritten as follows:Patrick and Giustino (2014)
(3) 
Here the product runs over all the normal coordinates . In this and all following expressions it is understood that the three translational modes with zero vibrational frequency are skipped in the sums. We indicate the vibrational frequency of the th normal mode by . The corresponding zeropoint vibrational amplitude is given by , where is a reference mass that we take equal to the proton mass. Using these conventions, the Gaussian widths in Eq. (3) are given by:
(4) 
where is the BoseEinstein occupation of the th mode. In the reminder of this manuscript we will concentrate on the expression for the WL dielectric function given by Eq. (3).
The configurational average appearing in Eq. (3) was evaluated in Ref. Zacharias et al., 2015 using importancesampling Monte Carlo integration.Patrick and Giustino (2013) More specifically, the Monte Carlo estimator of the integralCaflisch (1998) was evaluated by averaging over a set of atomic configurations in a Bornvon Kármán supercell. Each configuration in the set was generated according to the importance function . In Ref. Zacharias et al., 2015 it was remarked that, in the case of the optical spectrum of silicon, random samples were sufficient in order to converge the integral in Eq. (3). Furthermore, calculations performed using a single sample were found to be of comparable accuracy to fullyconverged calculations. Motivated by these observations, we decided to investigate in detail why the stochastic evaluation of Eq. (3) requires only very few samples.
In Sec. III we provide a formal proof of the fact that, in the limit of large supercell, only one atomic configuration is enough for evaluating Eq. (3). In the reminder of this section we only give the optimal configuration and outline the calculation procedure, so as to place the emphasis on our main results.
In order to calculate the optical absorption spectrum (including band gap renormalization) at finite temperature using a oneshot frozenphonon calculation, we proceed as follows:

We consider a supercell of the primitive unit cell. We determine the interatomic force constants Giustino () by means of densityfunctional perturbation theory calculations in the primitive unit cell, using a Brillouinzone grid.Baroni et al. (2001a); Giannozzi et al. (2009)

By diagonalizing the dynamical matrix obtained from the matrix of force constants, we determine the vibrational eigenmodes and eigenfrequencies ( and indicate the atom and the Cartesian direction, respectively).

For a given temperature , we generate one distorted atomic configuration by displacing the atoms from equilibrium by an amount , with:
(5) In this expression is the mass of the th nucleus, and the sum runs over all normal modes. The vibrational modes are assumed to be sorted in ascending order with respect to their frequencies. In order to enforce the same choice of gauge for each vibrational mode, the sign of each eigenvector is chosen so as to have the first nonzero element positive. The prescription given by Eq. (5) will be motivated in Sec. III.

We calculate the dielectric function using the atomic configuration specified by Eq. (5). The result will be the temperaturedependent dielectric function at the temperature .

We check for convergence by repeating all previous steps using increasingly larger supercells.
In Fig. 1 we present the roomtemperature optical absorption coefficients of Si, C, and GaAs calculated using the procedure just outlined (red solid lines), and we compare our results with experiment (grey discs). For completeness we also show the absorption coefficients evaluated with the atoms clamped at their equilibrium positions (blue solid line). The calculations were performed on supercells, using densityfunctional theory (DFT) in the local density approximation (LDA) and a scissor correction; all computational details are provided in Sec. VII. Our present results for Si in Fig. 1(a), obtained using a single atomic configuration, are in perfect agreement with those reported in Ref. Zacharias et al., 2015 using importancesampling Monte Carlo integration. In Fig. 1(a) we see that the spectrum calculated with the nuclei clamped in their equilibrium positions exhibit absorption only above the direct gap, as expected. At variance with these calculations, our oneshot calculation based on the WL theory correctly captures indirect transitions. In particular, this calculation is in good agreement with experiment throughout a wide range of photon energies.Green and Keevers (1995); Phillip and Taft (1964) This agreement surprisingly extends over seven orders of magnitude of the absorption coefficient. However, we should point out that our adiabatic theory does not capture the fine structure close to the absorption onset: there the nonadiabatic HBB theory gives two different slopes for phonon absorption and emission processes.Macfarlane et al. (1958); Noffsinger et al. (2012) Near the onset for direct transitions, our calculation underestimates the experimental data. This behavior is expected since we are not including electronhole interactions, which are known to increase the oscillator strength of the peak near 3.4 eV. Benedict et al. (1998)
In Fig. 1(b) we show our WL calculation for diamond. In this case our method correctly captures the absorption in the indirect range, however we observe more pronounced deviations between theory and experiment than in the case of Si. We assign the residual discrepancy to the inability of DFT/LDA to accurately describe the joint density of states of diamond. In fact, in contrast to the case of silicon, in diamond a simple scissor correction is not enough to mimic quasiparticle corrections.Lambert and Giustino (2013) For example, the GW corrections to the , , and states of Si are all in the narrow range between 0.66 and 0.75 eV; instead the GW corrections to the same states of diamond span a broader range, between 1.47 and 2.04 eV.Lambert and Giustino (2013) This is expected to lead to a significant redistribution of spectral weight precisely in the range of photon energies considered in Fig. 1(b). Also in this case the strength of the peak near 7.3 eV is underestimated due to our neglecting electronhole interactions.Benedict et al. (1998) In Fig. 1(b) we are reporting two sets of experimental data. Phillip and Taft (1964); Clark et al. (1964) These data exhibit different intensities near the absorption edge (at energies eV). According to Ref Phillip and Taft, 1964, the intensity of the absorption coefficient for energies below 6.5 eV is not fully reliable. Our calculated spectrum is in closer agreement with the data from Ref Clark et al., 1964, which exhibit a sharper absorption edge. This comparison suggests that our present method might prove useful for the validation of challenging experiments, especially near the weak absorption edge. As in the case of Fig. 1(a), the fine structure features close to the absorption onset are absent in our calculation.
In Fig. 1(c) we compare the absorption coefficient calculated for GaAs with experiment. Aspnes and Studna (1983); Sturge (1962) This example clearly demonstrates that our WL calculation correctly describes the absorption spectrum of a directgap semiconductor. In this case the shape of the absorption coefficient is not altered, as expected, but the spectrum is redshifted as a result of the zeropoint renormalization of the band structure. Also in the case of GaAs the calculated absorption coefficient underestimates the measured values. This is partially a consequence of our neglecting of excitonic effects,Rohlfing and Louie (1998) but most importantly it is a consequence of the inability of DFT/LDA to accurately describe the effective masses of GaAs. In fact, according to the standard theory of absorption in directgap semiconductors,Bassani and Parravicini (1975) the absorption coefficient scales as , where is the average isotropic effective mass. Since DFT calculations for GaAs yield masses which are up to a factor of 2 smaller than in experiment, we expect a corresponding underestimation of the absorption coefficient by up to a factor of . This estimate is in line with our results in Fig. 1(c). We note that this issue is not found in the case of Si [Fig. 1(a)], since the DFT/LDA effective masses of silicon are in surprising agreement with experiment.Giustino (2014) We also note that for GaAs our calculation correctly predicts phononassisted absorption below the direct gap. This phenomenon is related to the Urbach tail.Urbach (1953)
After discussing the comparison between our calculations and experiment, we briefly comment on the computational effort and the numerical convergence. Figure 2 shows the imaginary part of the WL dielectric function of Si, C, and GaAs evaluated at zero temperature using the procedure outlined in the previous page. In order to achieve convergence we performed calculations for supercells of increasing size, from to . It is clear that large supercells are required in order to obtain converged results. Increasing the supercell size has the twofold effect of refining the sampling of the electronphonon coupling in the equivalent Brillouin zone, and of approaching the limit where Eq. (5) becomes exact. In Fig. 2 we also report the band gaps of Si, C, and GaAs as extracted from using the standard Tauc plots.Tauc et al. (1966) These results are in agreement with previous work and will be discussed in Sec. VI. From this figure we see that our methodology correctly describes the zeropoint renormalization of the band gap of both direct and indirectgap semiconductors.
From Figs. 1 and 2 it should be apparent that, apart from the inherent deficiencies of the DFT/LDA approximation, with our new method it is possible to compute optical spectra and band gaps in both direct and indirect semiconductors including electronphonon interactions, at the cost of a single supercell calculation with clamped nuclei.
In the following sections we develop the theory underlying our computational approach, and we provide an extensive set of benchmarks.
Iii Formal justification of Eq. (5)
In this section we provide the rationale for the choice of the optimal configuration given by Eq. (5). We start from a heuristic argument, and then we provide a formal justification.
iii.1 Heuristic approach
According to Eq. (3), the WL dielectric function can be interpreted as the average of over the standard normal random variables . Let us consider the sum of the squares of these variables, . The random variable follows by construction the distribution.Koch (1999) Owing to the central limit theorem, the distribution tends to a normal distribution for . More specifically, in the limit of large the variable follows a standard normal distribution. As a result, as increases, the variable becomes strongly peaked at , with a standard deviation .Abramowitz and Stegun (1972) As a sanity check, we verified these limits numerically by generating one million random atomic configurations for a supercell of diamond. Based on these considerations, we infer that for large the integration in Eq. (3) is dominated by atomic configurations such that . This same conclusion can alternatively be reached be rewriting the integral in Eq. (3) as the product of an integral over the ‘radial’ variable , and an integral over a generalized ‘angular’ variable which runs over the dimensional sphere of radius .
In the absence of information about the electronphonon coupling constants of each vibrational mode, we must assume that all modes are equally important in the evaluation of the integral in Eq. (3). Therefore the most representative sets of coordinates for evaluating the integral in Eq. (3) are those satisfying the condition with . This is precisely what we observed in the importancesampling Monte Carlo calculations reported in Ref. Zacharias et al., 2015. A similar conclusion was reached in Ref. Monserrat, 2016a, where the concept of ‘thermal lines’ was introduced. Here we do not follow up on the idea of thermal lines since, as we prove below, there exists one atomic configuration which yields the exact temperaturedependent dielectric function in the limit of large . If we were to choose a single configuration to evaluate the integral in Eq. (3), in absence of information about the electronphonon couplings the leastbiased choice would correspond to taking random signs for each normal coordinate. This reasoning formed the heuristic basis for the choice made in Eq. (5).
iii.2 Formal proof
We now proceed to demonstrate that Eq. (5) is not just a sound approximation, but it is indeed the optimal configuration for evaluating Eq. (3) using a oneshot calculation. To this aim we perform a Taylor expansion of in the variables , and then evaluate each integral in Eq. (3) analytically. The result is:
(6) 
In this expression denotes the dielectric function evaluated for nuclei clamped in their equilibrium positions, and the term is a shorthand notation to indicate all terms of the kind and higher powers.
We now consider the dielectric function calculated with the nuclei clamped in the positions specified by Eq. (5). We denote this function by , with ‘1C’ standing for ‘one configuration’. Another Taylor expansion in the normal mode coordinates yields:
(7) 
By comparing Eqs. (6) and (III.2) we see that and do coincide up to if the following conditions hold: (i) the summations on the first and third lines of Eq. (III.2) vanish; (ii) all terms in the second line of the same equation vanish.
In general the conditions (i) and (ii) do not hold, therefore calculations of and will yield very different results. However, these conditions are fulfilled in the limit , as we show in the following. To this aim let us consider the summation on the first line of Eq. (III.2). We focus on two successive terms in the sum, and . Assuming a vibrational density of states (vDOS) which is nonvanishing up to the highest vibrational frequency, when , then , therefore these modes are effectively degenerate, and hence must exhibit the same electronphonon coupling coefficients. Under these conditions we can write:
(8) 
This reasoning can be repeated for every pair of vibrational modes appearing in the first line of Eq. (III.2). If is even, then this proves that the sum on the first line vanishes. If is odd, then there is one mode left out, but the contribution of this one mode is negligibly small for . If there are gaps in the vDOS, then the above reasoning remains valid by considering separately the frequency ranges where the vDOS is nonzero. This completes the proof that the sum in the first line of Eq. (III.2) vanishes in the limit of large .
The summation in the third line of Eq. (III.2) can be analyzed along the same lines, after noticing that one can rearrange the sum as follows:
(9) 
Also in this case we can consider any pair of successive eigenmodes and and repeat the reasoning made above for the first line of Eq. (III.2). The result is that for the entire sum must vanish.
If we now consider the second line of Eq. (III.2), the sum of the terms with must vanish for large . In fact, the second derivatives and enter the sum with opposite signs, therefore their contribution vanishes in the limit . On the other hand, when , two successive terms in the sum contribute with the same sign, yielding . These contributions lead precisely to the second term on the r.h.s. of Eq. (6).
Taken together, Eqs. (6)(9) and the above discussion demonstrate that our singleconfiguration dielectric function, , and the exact WL dielectric function, , do coincide to for . It is not difficult to see that this result can be generalized to all orders in , therefore the following general statement holds true: in the limit of large supercell the dielectric function evaluated with the nuclei clamped in the configuration specified by Eq. (5) approaches the WL dielectric function, Eq. (3). In symbols:
(10) 
This result forms the basis for the methodology presented in this manuscript. The importance of the equivalence expressed by Eq. (10) resides in that it allows us to calculate dielectric functions at finite temperature using a single atomic configuration. This represents a significant advance over alternative techniques such as for example pathintegrals molecular dynamics, importancesampling Monte Carlo, or the direct evaluation of each term in Eq. (6) using frozenphonon calculations for each vibrational mode.
We emphasize that, while we have proven the limit in Eq. (10), we have no information about the convergence rate of towards the exact result . In principle this rate could be estimated a priori by inspecting the convergence of the Eliashberg function Allen and Cardona (1981) with the sampling of the Brillouin zone in a calculation within the primitive unit cell. In practice we found it easier to directly calculate for supercells of increasing size. Convergence tests for Si, C, and GaAs were reported in Fig. 2. It is seen that, for these tetrahedral semiconductors, converged results are obtained for supercells. We emphasize that in Fig. 2, is given in logarithmic scale: in a linear scale the differences between a and an calculation would be barely discernible.
Calculations using the largest supercells in Fig. 2 are obviously timeconsuming. However, one should keep in mind that each line in this figure corresponds to a single calculation of the dielectric function at clamped nuclei, therefore this method enables the incorporation of temperature at low computational cost, and removes the need of configurational sampling.
In the insets of Fig. 2 we also presented the zeropoint renormalization of the fundamental gaps of Si, C, and GaAs. These gaps were directly obtained from the calculated using the standard Tauc plots.Tauc et al. (1966) Details will be discussed in Sec. VI; here we only point out that also the band gaps are calculated by using a single atomic configuration in each case. To the best of our knowledge this is the first time that calculations of temperaturedependent band gaps using a single atomic configuration have been reported.
Iv Improvements using configurational averaging
There might be systems for which calculations using large supercells are prohibitively timeconsuming, and the limit in Eq. (10) is practically beyond reach. For these cases it is advantageous to extend the arguments presented in Sec. III to calculations using more than one atomic configuration. In this section we show how it is indeed possible to construct a hierarchy of atomic configurations in order to systematically improve the numerical evaluation of Eq. (3) at fixed supercell size.
We restart by considering atomic configurations specified by the normal coordinates , where the signs or are yet to be specified. A Taylor expansion as in Eq. (III.2) yields:
(11) 
Here the superscript denotes the entire set of signs, . With this notation, the choice expressed by Eq. (5) corresponds to setting .
In the language of stochastic sampling, the configurations specified by and are called an ‘antithetic pair’.Caflisch (1998) It is immediate to see that the dielectric function calculated using both and does not contain any odd powers of . In fact from Eq. (IV) we have:
(12) 
This result is already very close to the exact expansion in Eq. (6), independent of the size of the supercell.
In order to obtain Eq. (6) exactly, one would need to further eliminate all terms in the sum on the second line of Eq. (IV). This elimination can be achieved systematically by considering an additional configuration defined as follows:

(13) 
where for the sake of clarity we considered the case . In practice is simply obtained by swapping all the signs of the secondhalf of the vector . It is immediate to verify that the dielectric function calculated by averaging the 4 atomic configurations specified by , , , and contains only half of the terms that appear in Eq. (IV). Since each individual configuration satisfies the asymptotic relation in Eq. (10), it is clear that the above calculation using 4 atomic configurations will approach the exact WL dielectric function faster as the supercell size increases.
This strategy can also be iterated by partitioning each subset in Eq. (13) in two halves, and applying a sign swap on two out of four of the resulting subsets. At each level of iteration the number of atomic configurations doubles, and the number of remaining terms in Eq. (IV) halves. For example, it is easy to verify that, by using 4, 16, and 64 configurations generated in this way, it is possible to eliminate 50%, 87.5%, and 96.875% of the terms in Eq. (IV), respectively.
In Fig. 3(a) we show the effect of using Eqs. (5) and (IV) for the calculation of the dielectric function of silicon at 300 K via 2 distinct atomic configurations. In Fig. 3(b) we repeat the calculations, this time using 4 distinct configurations, according to Eqs. (5), (IV), and (13). Here we see that increasing the number of atomic configurations suppresses spurious fluctuations in the spectra; the effect is most pronounced for the smallest supercell, which corresponds to Si unit cells. We note that, in Fig. 3, the curves corresponding to the supercell have been rigidly shifted for clarity: without such a shift the curves are almost indistinguishable from the those calculated using a supercell.
Figure 3 clearly shows that, irrespective of the configurational averaging, too small supercells may not be enough to accurately evaluate the WL dielectric function. This is easily explained by considering that, in order to correctly describe an indirect absorption onset, we need a supercell which can accommodate phonons connecting the band extrema.
V Relation between the WilliamsLax theory, the AllenHeine theory, and the HallBardeenBlatt theory
Having described the conceptual basis of our methodology, we now establish the link between the WL dielectric function given by Eq. (3) and the standard theory of Hall, Bardeen, and Blatt of indirect optical absorption,Hall et al. (1954) as well as the theory of temperaturedependent band structures of Allen and Heine.Allen and Heine (1976) In this section we only present the main results, leaving the mathematical details to Appendices A and B.
The HBB theory describes indirect optical absorption by means of timedependent perturbation theory, and the final expression [see Eq. (24) below] involves the momentum matrix elements evaluated for the nuclei clamped in their equilibrium positions, , and the linear electronphonon matrix elements:
(14) 
Here denotes a KohnSham state and is the variation of the KohnSham potential with respect to the normal mode coordinate . In order to make these quantities explicit in Eq. (2), we expand the momentum matrix elements to first order in the atomic displacements, and the energies to second order in the displacements. Using RaleighSchrödinger perturbation theory we find:
(15) 
where the primed summation indicates that we skip terms such that or . Similarly the expansion of the energies yields, for example:
(16)  
with being the ‘DebyeWaller’ electronphonon matrix element:Allen and Heine (1976); Allen and Cardona (1981); Patrick and Giustino (2014); Giustino ()
(17) 
To make contact with the AH theory of temperaturedependent band structures and with the HBB theory of indirect absorption, we analyze separately the cases of (i) direct gaps and (ii) indirect gaps.
v.1 Direct gaps
When the gap is direct, the optical matrix elements in Eq. (15) are dominated by the term . In fact, if we denote by the characteristic electronphonon matrix element and by the minimum gap, and we note that , then the terms in the square brackets are , and hence can be neglected next to . This approximation is implicitly used in all calculations of optical absorption spectra which do not include phononassisted processes. By using this approximation in Eqs. (2) and (3) we obtain:
(18) 
where we have defined , , and . Equation (V.1) can be rewritten by exploiting the Taylor expansion of the Dirac delta distribution in powers of . The derivation is laborious and is reported in Appendix A. The final result is:
(19) 
In this expression, denotes the temperaturedependent electron energy in the AllenHeine theory:Allen and Heine (1976)
(20) 
In particular, the first term in the square brackets is the Fan selfenergy correction, while the second term is the DebyeWaller correction.Antončík (1955); Walter et al. (1970); Allen and Heine (1976); Allen and Cardona (1981); Cardona and Thewalt (2005) The quantity in Eq. (V.1) is the width of the optical transition, and is defined as follows:
(21) 
Equation (V.1) shows that, in the case of direct absorption processes, the WL theory yields a dielectric function which exhibits normalized peaks at the temperaturedependent excitation energies . We note that the width of the optical transitions is similar to but does not coincide with the electronphonon linewidth obtained in timedependent perturbation theory, compare for example with Eq. (169) of Ref. Giustino, . This subtle difference is a direct consequence of the semiclassical approximation upon which the WilliamsLax theory is based.
To the best of our knowledge the connection derived here between the WilliamsLax theory and the AllenHeine theory, as described by Eqs. (V.1) and (20), is a novel finding. In particular, the present analysis demonstrates that, to lowest order in perturbation theory, the WL theory of optical spectra yields the adiabatic version of the AH theory of temperaturedependent band structures.
v.2 Indirect gaps
In the case of indirect semiconductors, owing to the momentum selection rule, optical transitions near the fundamental gap are forbidden in the absence of phonons.Yu and Cardona (2010) By consequence, the optical matrix elements at equilibrium must vanish, and we can set in Eq. (15). This observation represents the starting point of the classic HBB theory of indirect optical processes.Hall et al. (1954) In this case Eq. (3) becomes:
(22) 
We stress that in this expression we are only retaining those terms in Eq. (15) which are linear in the normal coordinates. This choice corresponds to considering only onephonon processes, as it is done in the HBB theory. While multiphonon processes are automatically included in the WL theory, in the following we do not analyze them explicitly.
By performing a Taylor expansion of the Dirac delta appearing in Eq. (V.2) in powers of , as we did for the case of direct gaps, we arrive at the following expression:
(23) 
The derivation of this result is lengthy, and is reported for completeness in Appendix B.
Equation (V.2) demonstrates that the WL theory correctly describes indirect optical absorption. In fact, if we neglect the broadening and the temperature dependence of the band structure, we obtain the theory of Hall, Bardeen, and Blatt:Hall et al. (1954); Yu and Cardona (2010); Noffsinger et al. (2012)
(24)  
The only difference between this last expression and the original HBB theory is that here the Dirac delta function does not contain the phonon energy. Physically this corresponds to stating that the WL theory provides the adiabatic limit of the HBB theory of indirect absorption.
To the best of our knowledge, this is the first derivation of the precise formal connection between the WL theory of indirect absorption, as expressed by Eq. (V.2), and the theories of Hall, Bardeen, and Blatt and of Allen and Heine, as given by Eqs. (24) and (20).
We emphasize that the adiabatic HBB theory does not incorporate the temperature dependence of the band structure, as it can be seen from Eq. (24). Instead, the WL theory includes band structure renormalization by default, see Eq. (V.2). This point is very important in view of performing predictive calculations at finite temperature.
From Eq. (V.2) we can also tell that, in order to generalize the HBB theory to include temperaturedependent band structures, one needs to incorporate the temperature dependence only in the energies corresponding to real transitions, i.e. in the Dirac deltas in Eq. (24), and not in the energies corresponding to virtual transitions, i.e. the energy denominators in the same equation.
In order to illustrate the points discussed in Secs. V.1 and V.2, we show in Fig. 4 a comparison between the absorption spectra of silicon and diamond calculated using either the HBB theory or the WL theory. For the HBB calculations we employed the adiabatic approximation and we recast Eq. (24) in the following equivalent form:
(25) 
The equivalence between this expression and Eq. (24) is readily proven by using Eqs. (15) and (B). For the evaluation of the second derivatives for each vibrational mode we used finitedifference formulas; this operation requires frozenphonon calculations. In the examples shown in Fig. 4 we employed supercells, corresponding to 768 calculations on supercells containing 128 atoms. We emphasize that the WL spectrum requires instead only 1 calculation.
Figure 4(a) compares the measured absorption coefficient of silicon (grey dots) with those calculated using the HBB theory (red line) and the WL theory (blue line). All data are for the temperature K. Here we see that the WL and HBB calculations yield similar spectra. However, while the optical absorption onset in the HBB theory coincides with the indirect band gap of silicon at equilibrium, the WL spectrum is redshifted by an amount 0.1 eV, corresponding to the zeropoint renormalization and the temperature shift of the gap. As a result, the WL calculation is in better agreement with experiment. We emphasize again that the agreement between theory and experiment remarkably extends over a range spanning 6 orders of magnitude, and the curves shown in Fig. 4 do not carry any empirical scaling factors.
The effect of band gap renormalization becomes more spectacular in the case of diamond, as shown in Fig. 4(b). In this case the HBB calculation misses the absorption onset by as much as 0.5 eV, while the WL theory yields an onset in agreement with experiment. The discrepancy between the prediction of the adiabatic HBB theory and experiment is understood as the result of the very large electronphonon renormalization of the band gap of diamondGiustino et al. (2010); Cannuccia and Marini (2011); Antonius et al. (2014). In addition, the adiabatic version of the HBB theory also misses the small redshift associated with phononemission processes.Noffsinger et al. (2012)
We stress that the WL theory does not come without faults. The main shortcoming is that, unlike the HBB theory, it does not capture the fine structure corresponding to phonon absorption and emission processes near the absorption onset. This is a direct consequence of the semiclassical approximation underpinning the WL theory, whereby the quantization of the final vibrational states is replaced by a classical continuum.Lax (1952)
Vi Band gap renormalization
vi.1 Comparison between calculations using the WilliamsLax theory and HallBardeenBlatt theory
In Sec. V we demonstrated that the imaginary part of the temperaturedependent dielectric function in the WL theory exhibits an absorption onset at the temperaturedependent band gap given by the AH theory, see Eqs. (V.1) and (V.2). This result suggests that it should be possible to extract temperaturedependent band gaps directly from calculations of or , as it is done in experiments.
Using Eqs. (V.1) and (V.2) and the standard parabolic approximation for the band edges of threedimensional solids,Bassani and Parravicini (1975) the following relations can be obtained after a few simple manipulations:
direct:  (26)  
indirect:  (27) 
Here is the temperaturedependent band gap, and these relations are valid near the absorption onset. Equations (26) and (27) simply reflect the jointdensity of states of semiconductors, and form the basis for the standard Tauc plots which are commonly used in experiments in order to determine band gaps.Tauc et al. (1966) These relations are employed as follows: after having determined , one plots for directgap materials, or for indirect gaps. The plot should be linear near the absorption onset, and the intercept with the horizontal axis gives the band gap .
In the case of indirectgap semiconductors it is not possible to use Eq. (26) in order to determine the direct band gap. Nevertheless, in these cases one can still determine the gap by analyzing secondderivative spectra of the real part of the dielectric function, . These spectra exhibit characterisic dips, which can be used to identify the direct gap. This is precisely the procedure employed in experiments in order to measure direct band gaps.Lautenschlager et al. (1987a); Logothetidis et al. (1992); Lautenschlager et al. (1987b) In the following we use this lineshape analysis to calculate the band gaps of Si, C, and GaAs.
vi.2 Indirect semiconductors: silicon
Figure 5(a) shows the Tauc plots obtained for silicon using the WL theory. As expected from Eq. (26), we obtain straight lines over an energy range of almost 1 eV from the absorption onset. From the intercept of linear fits taken in the range 01.8 eV we obtain the indirect band gaps at several temperatures; the results are shown in Fig. 5(b) as red discs, and compared to experiment (grey discs). The agreement with experiment is good, with the exception of a constant offset which relates to our choice of scissor correction for the band structure of silicon (cf. Sec. VII).
In order to minimize numerical noise, we determine the zeropoint renormalization of the band gap, , as the offset between the square root of the joint density of states calculated at equilibrium and that at K. This refinement is discussed in Appendix C. In this case we find meV. Small changes of this value are expected for larger supercells and denser Brillouinzone sampling. For completeness we show in Fig. 5(c) the convergence of the zeropoint renormalization of the indirect gap with the number of points in the supercell.
Our calculated zeropoint renormalization of the indirect gap of silicon is in good agreement with previous calculations and with experiment. In fact, Ref. Monserrat and Needs, 2014 reported a renormalization of 60 meV using finitedifferences supercell calculations; Ref. Monserrat, 2016a reported 58 meV using Monte Carlo calculations in a supercell. Ref. Poncé et al., 2015 obtained a zeropoint renormalization of 56 meV using the perturbative AllenHeine approach. Measured values of the renormalization range between 62 and 64 meV.Cardona (2001, 2005); Cardona and Thewalt (2005)



Figure 5(d) shows the secondderivative spectra, , calculated for silicon at two temperatures. Following Ref. Lautenschlager et al., 1987a, we determine the energy of the transition using the first dip in the spectra. In Fig. 5(e) we compare the direct band gaps thus extracted with the experimental values. Apart from the vertical offset between our data and experiment, which relates to the choice of the scissor correction (cf. Sec. VII), the agreement with experiment is good. In order to determine the zeropoint renormalization of the direct gap, we took the difference between the dips of the secondderivative spectra calculated at equilibrium and using the WL theory at K. We obtained a zeropoint renormalization of meV, which compares well with the experimental range meV.Lautenschlager et al. (1987a) Values in the same range were reported in previous calculations, namely 28 meV from Ref. Monserrat, 2016b and 42 meV from Ref. Poncé et al., 2015. We note that, in Ref. Monserrat, 2016b, a GW calculation on a 444 supercell yielded a renormalization of 53 meV. For completeness we show in Fig. 5(f) the convergence of the zeropoint renormalization with respect to Brillouinzone sampling.
vi.3 Indirect semiconductors: diamond
Figure 6(a) shows the Tauc plots calculated for diamond using the WL approach. As in the case of silicon, we determined the indirect gap as a function of temperature by means of the intercept of the linear fits with the horizontal axis; the results are shown in Fig. 6(b) together with the experimental data of Ref. Clark et al., 1964. From this comparison we see that the agreement with experiment is reasonable, however the calculations underestimate the measured renormalization. This effect has been ascribed to the fact that DFT/LDA underestimates the electronphonon matrix elements as a consequence of the DFT band gap problem.Antonius et al. (2014) Our calculated zeropoint renormalization of the indirect gap is meV. This result was obtained following the procedure in Appendix C. Our value for the zeropoint renormalization is compatible with previously reported values based on DFT/LDA, namely 330 meV,Poncé et al. (2015) 334 meV,Monserrat and Needs (2014) 343 meV,LloydWilliams and Monserrat (2015) and 344 meV.Monserrat (2016a). Our calculations are in good agreement with the experimental values of 340 meV and 370 meV reported in Refs. Cardona, 2001, 2005. However, the use of a more recent extrapolation of the experimental data yields a renormalization of 410 meV, Monserrat et al. (2014) which is 65 meV larger than our result. Since it is known that GW quasiparticle corrections do increase the zeropoint renormalization as compared to DFT/LDA,Antonius et al. (2014); Monserrat (2016b) we expect that by repeating our WL calculations in a GW framework our results will be in better agreement with the experimental data. For completeness in Fig. 6(c) we show the convergence of the zeropoint renormalization with Brillouinzone sampling (for an 888 supercell).
Figure 6(d) shows the secondderivative spectra of the real part of the dielectric function of diamond, as calculated from the WL theory. The direct gap of diamond was obtained in Ref. Logothetidis et al., 1992 using the deep minimum in the experimental curves; here we follow the same approach, and we report our results in Fig. 6(e). For comparison we also show the experimental data from Ref. Logothetidis et al., 1992. Apart from the vertical offset which reflects our choice of scissor correction, the calculations are in reasonable agreement with experiment. In particular we determined a zeropoint renormalization of 450 meV, to be compared with the experimental ranges of meV (sample II of Ref. Logothetidis et al., 1992) and meV (sample II of Ref. Logothetidis et al., 1992). Our value of the zeropoint renormalization is compatible with previous calculations at the DFT/LDA level, for example Ref. LloydWilliams and Monserrat, 2015 reported 400 meV (when using an 888 Supercell), Ref. Poncé et al., 2014b reported 409 meV, Ref. Monserrat, 2016b reported 410 meV (value extracted from Fig.3) and Ref. Poncé et al., 2015 reported 416 meV.
We also point out that a previous work by one of us on the electronphonon renormalization in diamond reported a correction of 615 meV for the direct gap using DFT/LDA.Giustino et al. (2010) The origin of the overestimation obtained in Ref. Giustino et al., 2010 might be related to numerical inaccuracies when calculating electronphonon matrix elements for unoccupied KohnSham states at very highenergy, although this point is yet to be confirmed. Recent work demonstrated that GW quasiparticle corrections lead to an increase of the zeropoint renormalization of the direct gap.Monserrat (2016b); Antonius et al. (2014) Calculations of optical spectra using the WL theory and the GW method are certainly desirable, but lie beyond the scope of the present work.
For completeness, in Appendix D we also investigate multiphonon effects. In particular, we prove that WL calculations correctly yield the generalization of the adiabatic AH theory to the case of twophonon processes, and we demonstrate that multiphonon effects provide a negligible contribution to the band gap renormalization.
We also note incidentally that, in the case of small supercells, the band gap renormalization includes an additional spurious contribution when the band extrema are degenerate. This effect arises from a linearorder electronphonon coupling which lifts the band degeneracy, and has already been observed in pathintegral molecular dynamics calculations on diamond using a supercell with 64 atoms. Ramírez et al. (2006) This effect can easily be seen in the density of states in Fig. 7 as three separate peaks near the valence band edge. In the limit of large supercells this effect vanishes, since it arises from zonecenter phonons, whose weight becomes negligibly small as .
vi.4 Direct semiconductors: gallium arsenide
Figure 8 shows our WL calculations of the temperaturedependent absorption onset and band gap of GaAs. Since GaAS is polar, we included the nonanalytical part of the dynamical matrix in the calculations of the vibrational frequencies and eigenmodes.Wang et al. (2012) In this case we also took into account the thermal expansion of the lattice, which is not negligible for this semiconductor.Walter et al. (1970) To this aim we performed calculations within the quasiharmonic approximation,Baroni et al. (2001b) that is we repeated the calculations of the normal modes for each temperature, by varying the lattice constant according to the measured thermal lattice expansion coefficient.Clec’h et al. (1989); Pierron et al. (1967)
Indirect gap  Direct gap  

Present  Previous  Experiment  Present  Previous  Experiment  
Si  57  56^{1}^{1}1Ref. Poncé et al., 2015 (757575), 58^{2}^{2}2Ref. Monserrat, 2016a (666), 60^{3}^{3}3Ref. Monserrat and Needs, 2014 (555)  62^{4}^{4}4Ref. Cardona, 2001, 64^{5}^{5}5Ref. Cardona,2005  44  28^{6}^{6}6Ref. Monserrat, 2016b (444), 42^{1}^{1}footnotemark: 1  ^{7}^{7}7Ref. Lautenschlager et al.,1987a 
C  345  330^{1}^{1}footnotemark: 1, 334^{8}^{8}8Ref. Monserrat and Needs, 2014 (666), 343^{9}^{9}9Ref. LloydWilliams and Monserrat, 2015 (484848), 344^{2}^{2}footnotemark: 2  370^{5}^{5}footnotemark: 5, 410^{10}^{10}10Ref. Monserrat et al.,2014  450  409^{11}^{11}11Ref. Poncé et al., 2014b, 410^{6}^{6}footnotemark: 6, 416^{1}^{1}footnotemark: 1, 430^{9}^{9}footnotemark: 9  ^{12}^{12}12Ref. Logothetidis et al., 1992, ^{12}^{12}footnotemark: 12 
GaAs  32  23^{13}^{13}13Ref. Antonius et al., 2014 (444)  ^{14}^{14}14Ref. Lautenschlager et al.,1987b 
In Fig. 8(a) we show the direct absorption onset. In this case we extract the gap from straightline fits of , after Eq. (26). Since the range where the function is linear is rather narrow (due to the presence of lowlying conduction band valleys), we refine the calculation of the zeropoint renormalization using the more accurate procedure discussed in Appendix C.
In Fig. 8(b) we show the calculated band gaps as a function of temperature and we compare our results to experiment. For completeness we report calculations performed without considering the thermal expansion of the lattice. Here we see that our calculations are in reasonable agreement with experiment, and that lattice thermal expansion is definitely not negligible. Using an 888 supercell we obtained a zeropoint renormalization of 32 meV. Our value is in line with previous calculations, yielding 23 meV,Antonius et al. (2014) as well as with the experimental range of meV. We expect that GW quasiparticle corrections will further increase the zeropoint correction by 10 meV.Antonius et al. (2014) For completeness in Fig. 8(c) we show the convergence of the zeropoint renormalization with the sampling of the Brillouin zone (in an 888 supercell).
Recently it was pointed out that, in the case of polar semiconductors, calculations based on the AllenHeine theory exhibit a spurious divergence when the quasiparticle lifetimes are set to zero.Poncé et al. (2015); Nery and Allen () The origin of this artifact relates to the Frölich electronphonon coupling,Verdi and Giustino (2015); Sjakste et al. (2015) and has been discussed in Ref. Giustino, . In the present calculations it is practically impossible to test whether we would have a singularity in the limit of very large spercells. However, we speculate that our calculations should not diverge, since the Bornvon Kármán boundary conditions effectively shortcircuit the longrange electric field associated with longitudinal optical phonons. This aspect will require separate investigation.
In Table 1 we summarize all our calculations of zeropoint renormalization for Si, C, and GaAs, and compare our present results with previous theory and experiment.
Vii Computational setup
All calculations were performed within the local density approximation to density functional theory.Ceperley and Alder (1980); Perdew and Zunger (1981) We used normconserving pseudopotentials,Fuchs and Scheffler (1999) as implemented in the Quantum ESPRESSO distribution.Giannozzi et al. (2009) The KohnSham wavefunctions were expanded in planewave basis sets with kinetic energy cutoffs of 40 Ry, 50 Ry, and 120 Ry for Si, GaAs and C, respectively. The interatomic force constants in the Bornvon Kármán supercell were obtained as the Fourier transforms of the dynamical matrices calculated in the primitive unit cells via density functional perturbation theory.Baroni et al. (2001a)
All calculations of band gap renormalization were performed by using two atomic configurations: our optimal configuration, given by Eq. (5), and its antithetic pair, as obtained by exchanging the signs of all normal coordinates. This choice guarantees high accuracy in the lineshapes near the absorption onset. The absorption coefficients shown in Fig. 1 were calculated using , where is the speed of light and the refractive index calculated as . The optical matrix elements including the commutators with the nonlocal components of the pseudopotential Starace (1971) were calculated using Yambo.Marini et al. (2009) In order to compensate for the DFT band gap problem, we rigidly shifted the conduction bands so as to mimic GW quasiparticle corrections. The scissor corrections were taken from previous GW calculations performed using the same computational setup for Si and C,Lambert and Giustino (2013) and from Ref. Remediakis and Kaxiras, 1999 for GaAs. In particular, we used eV, 1.64 eV, and 0.53 eV for Si, C, and GaAs, respectively. The nonlocality of the scissor operator was correctly taken into account in the oscillator strengths Starace (1971) via the renormalization factors ; this ensures that the sum rule is fulfilled.
The dielectric functions were calculated by replacing the Dirac deltas in Eq.(2) with Gaussians of width 30 meV for Si and C, and 50 meV, for GaAs. All calculations presented in this work were performed using 888 supercells of the primitive unit cell, unless specified otherwise. The sampling of the Brillouin zone of each supercell was performed using random points, with weights determined by Voronoi triangulation.Aurenhammer (1991) In order to sample the 888 supercells of Si, C and GaAs we used 40, 40, and 100 random points, respectively.
The expansion of the crystalline lattice with temperature was only considered in the case of GaAs. This choice is justified by the fact that the calculated band gaps of silicon and diamond change by less than 2 meV when using the lattice parameters at K and K. The temperaturedependence of the lattice parameters of Si and C was taken from Ref. Okada and Tokumaru, 1984; Sato et al., 2002.
Viii Conclusions and outlook
In this manuscript we developed a new ab initio computational method for calculating the temperaturedependent optical absorption spectra and band gaps of semiconductors, including quantum zeropoint effects. The present work significantly expands the scope of our previous investigation in Ref. Zacharias et al., 2015, by completely removing the need for stochastic sampling of the nuclear wavefunctions. In particular we demonstrated, both using a formal proof and by means of explicit firstprinciples calculations, that in order to compute dielectric functions and band gaps at finite temperature it is sufficient to perform a single supercell calculation with the atoms in a welldefined configuration, as given by Eq. (5).
Using this new technique we reported the first calculations of the complete optical absorption spectra of diamond and gallium arsenide including electronphonon interactions, and we confirmed previous results obtained for silicon in Refs. Noffsinger et al., 2012; Zacharias et al., 2015. Our calculations are in good agreement with experiment at the level of lineshape, location of the absorption onset, and magnitude of the absorption coefficient. From these calculations we extracted the temperaturedependence and the zeropoint renormalization of the direct band gaps of Si, C, and GaAs, and of the indirect gaps of Si and C. Our calculations are in good agreement with previous theoretical studies.
Our present work relies on the WilliamsLax theory of optical transitions including nuclear quantum effects.Williams (1951); Lax (1952) For completeness we investigated in detail the formal relation between the WL theory, the AllenHeine theory of temperaturedependent band structures, and the HallBardeenBlatt theory of phononassisted optical absorption. We demonstrated that both the AH theory and the HBB theory can be derived as loworder approximations of the WL theory.
We emphasize that our present approach enables calculations of complete optical spectra at finite temperature, including seamlessly direct and indirect optical absorption. This feature is useful in order to calculate spectra which are directly comparable to experiment in absolute terms, i.e. without arbitrarily rescaling the absorption coefficient or shifting the absorption onset to fit experiment. Our methodology will also be useful for predictive calculations of optical spectra, for example in the context of highthroughput computational screening of materials.
It is natural to think that the present approach could be upgraded with calculations of electronic structure and optical properties based on the GW/BetheSalpeter approach. Indeed, as it should be clear from Eq. (3), our methodology holds unchanged irrespective of the electronic structure technique employed to describe electrons at clamped nuclei. Since the present approach requires only one calculation in a large supercell, it is possible that complete GW/BetheSalpeter calculations of optical spectra including phononassisted processes will soon become feasible. In this regard we note that, recently, Ref. LloydWilliams and Monserrat, 2015 proposed the socalled ‘nondiagonal’ supercells in order to perform accurate supercell calculations at a dramatically reduced computational cost. In the future, it will be interesting to investigate how to take advantage of nondiagonal supercells in order to make the present methodology even more efficient.
Finally, it should be possible to generalize our present work to other important optical and transport properties. In fact, the WL theory is completely general and can be used with any property which can be described by the Fermi Golden Rule. For example, we expect that generalizations to properties such as photoluminescence or Auger recombinationKioupakis et al. (2015) should be within reach.
Acknowledgements.
The research leading to these results has received funding from the the UK Engineering and Physical Sciences Research Council (DTA scholarship of M.Z. and grants No. EP/J009857/1 and EP/M020517/1), the Leverhulme Trust (Grant RL2012001), and the Graphene Flagship (EU FP7 grant no. 604391). The authors acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility (http://dx.doi.org/ 10.5281/zenodo.22558), the ARCHER UK National Supercomputing Service under the ‘AMSEC’ Leadership project and the PRACE for awarding us access to the dutch national supercomputer ’Cartesiusâ.Appendix A WilliamsLax expression for direct absorption
In this Appendix we outline the steps leading from Eq. (V.1) to Eq. (V.1). We start by introducing the compact notation:
(28) 
where the coefficients and are obtained from Eq. (16):
(29)  
(30)  
Equation (V.1) can be simplified by using the Taylor expansion of the Dirac delta with respect to the energy argument. In general we have:Feller (1971)
(31) 
therefore we can set and , and replace inside Eq. (V.1). We find:
(32)  
Now we replace from Eq. (28) and carry out the integrals in the coordinates . The resulting expression does not contain the normal coordinates any more, and the various derivatives of the delta function can be regrouped using Eq. (31) in reverse. The result is:
(33) 
where , and is the temperaturedependent electron energy given by Eq. (20).
In order to make Eq. (A) more compact, it is convenient to rewrite the second term inside the square brackets by performing a Taylor expansion of the Dirac delta around , and group the terms proportional to and higher order inside the term at the end. This step leads to:
(34) 
having defined:
(35) 
In the final step we note that the term in Eq. (A) acts so as to broaden the lineshape. This is seen by using the Fourier representation of the Dirac delta, . We find: