Computational Modeling of the Dielectric Function of Silicon Slabs

Jan 28, 2014 - The dynamical dielectric function of a silicon slab, in the region from near IR to UV light frequencies, is expected to vary with its t...
0 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/JPCC

Computational Modeling of the Dielectric Function of Silicon Slabs with Varying Thickness Tijo Vazhappilly and David A. Micha* Quantum Theory Project, Departments of Chemistry and of Physics, University of Florida, Gainesville, Florida 32611-8435, United States S Supporting Information *

ABSTRACT: The dynamical dielectric function of a silicon slab, in the region from near IR to UV light frequencies, is expected to vary with its thickness, and it is important to know whether its optical properties are similar to those of bulk silicon. Slabs of varying thickness are modeled starting from their atomic structure. Modeled Si(111) surfaces are terminated by hydrogen atoms to compensate the dangling bonds, and optical properties have been obtained for Si slabs with 4, 8, and 12 layers. Real and imaginary parts of the dielectric function are obtained from the polarization of the slab, expressed in terms of a delayed response function constructed from a reduced density matrix (RDM) which includes electronic dissipative effects due to coupling of photoexcited electrons to the substrate lattice and to electronic excitations. The related index of refraction and absorption coefficient have also been calculated from the above treatment. These optical properties are obtained using density functional theory (DFT) and plane wave basis sets to construct the equation of motion of a RDM, solved for steady light absorption. Both GGA (PBE) and hybrid (HSE) DFT exchange-correlation density functionals are employed to calculate the optical properties from the RDM. The imaginary part of the dielectric function is related to the light absorbance, and has been compared to measurements showing that better agreement is obtained with the HSE hybrid functional containing part of the exact short-range electronic exchange energy. We present a procedure by which one can reproduce the HSE results for the dielectric function from the computationally less expensive GGA PBE functional calculations, using a single photon energy shift parameter and results from PBE calculations. Our treatment shows that the onset of light absorption and strong diffraction are similar for thin slabs and bulk silicon, and that they have similar peak structure as functions of photon energy. Both properties increase with slab thickness at most photon energies. This makes silicon slabs reliable structures for photovoltaic applications.

1. INTRODUCTION The dielectric properties of semiconductor slabs are important in many aspects of electronics and optics, and for related devices. Their response to an applied electric field of given frequency contains delay effects and its frequency analysis leads to the introduction of a complex valued dynamical susceptibility function of frequency, with real and imaginary parts related to physically measurable properties. Optical properties are expected to vary with their thickness, and it is important to know whether they are similar to those of bulk silicon, particularly in connection with the search for thin film photovoltaic materials. These properties can be measured for bulk silicon materials but are seldom available for slabs. These on the other hand can be obtained from computer modeling as presented here. The real and imaginary parts of the dielectric function of semiconductors can be measured using inelastic light scattering or spectroscopic ellipsometry.1 The dielectric function varies with the frequency ω of the applied electric field, with features dependent on the electronic band structure of the solid, on © 2014 American Chemical Society

confinement effects due to physical boundary conditions, and on the presence of medium excitations such as plasmons, excitons, and phonons. The dielectric function is a tensor with Cartesian components εξη. In the applications presented here, the crystal supercells used in the modeling of slabs with surfaces perpendicular to a z-axis are orthorhombic, and as a result the dielectric tensor has only diagonal components εxx , εyy , εzz . A local electric vector (or electric displacement vector) with frequency ω for light along the z-direction and propagating along the x-direction in a dissipative medium, is of the form (loc) , (loc) z (x,t) = , z0 exp(ikx − ωt), where k(ω) = ωn(ω)/c, with c being the speed of light and n(ω) = η (ω) + iκ(ω) being a complex refractive index which can be measured from reflection experiments for varying photon energies.1,2 The imaginary part of the index is related to the absorption coefficient, and the dielectric function follows from the field wave equation as Received: October 25, 2013 Revised: January 27, 2014 Published: January 28, 2014 4429

dx.doi.org/10.1021/jp410579k | J. Phys. Chem. C 2014, 118, 4429−4436

The Journal of Physical Chemistry C

Article

εxx(ω) = n(ω)2. This identifies η with the index of refraction of the medium, and ωκ/c with a path absorbance of light along the x-direction. The optical absorption in a semiconductor is significant for photon energies larger than the band gap between its valence band (VB) and conduction band (CB), and occurs within an optical penetration length for the material. For silicon slabs limited by two surfaces, the surface sensitivity is higher for incident photon energies near the band gap, and surface properties show up in the dependence of the absorption with photon energies. The optical penetration depth for Si is much larger than the thickness of the slabs considered here, and we can ignore the effect of intensity reduction of the applied electric field inside the slab. Light absorption and photovoltages do not follow simple scaling laws with the thickness of the slab.3 The dependence of photoinduced surface electric dipole on the slab thickness has been discussed for Si slabs in our previous publication on light absorption4 and additional detailed information about scaling of band gap, absorption, surface dipole and surface photovoltage with slab thickness can be obtained from the present atomic structure models, also for amorphous silicon slabs.5 In most experimental studies, applied electromagnetic fields are weak and the induced polarization is linearly dependent on the field. This leads to a linear response,6 which can conveniently be obtained in terms of the statistical density operator for extended systems.7 In what follows, the standard linear response treatment is further extended to include electronic dissipation phenomena due to medium effects, by means of steady state solutions of equations of motion (EOM) for a reduced density operator.8,9 The polarization dependence on light frequency can be obtained from transition electric dipoles between electronic states which enter the calculation of the dielectric function, starting from the atomic structure of the slab. The present study deals with both real and imaginary parts of the dielectric function, their dependence on the light frequency and on their thickness of model slabs, to obtain its general dependence on absorbed photon energies, and slab thickness. Here we are combining ab initio electronic structure calculations and the statistical density operator, using density functional treatments, to generate basis sets of electronic orbitals, to construct the EOM of motion of a reduced density matrix (RDM), and to describe the linear response to steady state excitation by light absorption. Electronic structures are obtained from a density functional theory (DFT) suitable for extended systems with many atoms, and capable of generating excited orbitals with energies well above the excitation photon energies.10 Optical excitation in semiconductors can also be treated within a linear response theory using time-dependent DFT (TD-DFT) or Green’s function methods.10 However, in our present treatment of slabs with a periodic structure of large atomic supercells, it has been computationally advantageous to generate orbitals at the DFT stage with periodic boundary conditions, and then use a RDM to calculate several physical properties (absorption and surface photovoltage among them) in a unified way. Optical properties are obtained here within the Kohn−Sham (KS) DFT formalism and provides valuable information on quantum confinement effects.11,12 We use a large number of KS orbitals (KSOs) as an expansion basis set to obtain equations of motion and for our calculations of response functions and

dynamical susceptibilities. Details are given in section 2, where the equations are described in detail and expressions are given for stationary electronic states. The dissipative dynamics of photoexcited electrons, due to electron−boson (phonon or exciton) coupling with the medium, is described by rates obtained from perturbation theory or from measurements. There are more accurate theoretical treatments for the optical properties of semiconductors such as the Bethe-Salpeter and GW approaches including electronic correlation, and TD-DFT approaches corrected for electron self-interaction13,14 but they are computationally expensive for the large supercells of present interest.13,14 We have used in our DFT calculations the GGA functional PBE15 and the hybrid functional ωPBEh also known as HSE.16−18 It is well-known that in DFT, when using LDA/ GGA functionals, there is an electron self-interaction problem that affects the accuracy of results. To improve on this, the HSE functional has been introduced as a modified PBE functional containing part of the exact Hartree−Fock electronic exchange at short distances. The HSE functional has two parameters: a substituted 1/4 fraction of H−F exchange and an electron short separation range, called μ here, which are however independent of atomic structure. The HSE functional has been found to provide accurate results for properties of solids.19 In the GW approximation,20,21 a better description of electronic excitations can be obtained from the perturbative inclusion of many body effects through the self-energy. This GW self-energy can be used to correct either PBE or HF energies, but it involves substantially larger computational efforts. A simpler alternative to HSE and GW is a correction based on the partial subtraction of the electron−electron Coulomb interaction.20 This can rectify a problem with derivative discontinuities present in LDA/GGA functionals, but is not accurate for our applications. The HSE band gaps are close to the band gaps obtained from GW method for bulk silicon.22 From ref 22, the HSE band gap for bulk silicon is 1.14 eV, GW gives 1.23 eV, and from GW-BSE it is 1.23 eV (optical gap). The experimental quasiparticle and optical band gaps for bulk silicon are 1.17 and 1.16 eV, respectively. From ref 22, the HSE band gap for bulk silicon is 1.14 eV, GW gives 1.23 eV, and from GW-BSE it is 1.23 eV (optical gap). The experimental quasiparticle and optical band gaps for bulk silicon are 1.17 and 1.16 eV, respectively. In recent studies, it has been shown that hybrid functionals produce excellent results for static and dynamic screening properties of bulk Si.23 Therefore, one of the main aims of this paper is to compare the PBE and HSE results and to show how the PBE results can be corrected to provide better agreement with HSE and with experimental values from spectroscopic ellipsometry measurements for the dielectric function versus photon energies. In particular, the threshold of the absorbance function and locations of its peaks can be used as a criterion for accuracy of the DFT calculations.24,25 In what follows, section 2 presents the derivation of the polarization in terms of a time-dependent delayed response function, and its Fourier transform, obtained from a RDM formulation. This provides expressions for the dielectric function of frequency, absorbance, and refraction coefficient, including dissipative medium effects. Section 3 presents the real and imaginary parts of the dielectric function, the path absorbance, and index of refraction for varying thicknesses, using the PBE and HSE functionals for their comparison. Energy shifted results for optical properties are presented in 4430

dx.doi.org/10.1021/jp410579k | J. Phys. Chem. C 2014, 118, 4429−4436

The Journal of Physical Chemistry C

Article

This has been presented in detail in ref 4, where the elements have been obtained going from the Schrödinger picture (SP) representation of their time evolution to a representation with a frame rotating at frequency Ω. After a time average of the EOM over a measurement time Tm ≫ 2π/Ω, this leads to a rotating wave approximation (RWA) which is equivalent to using a basis set of electron-field states including one or no photons, and is sufficient for our calculations. The RWA elements, with the average over Tm indicated by an upper bar, are ρ̅(RW) (t) = σjk(t) kj and in the SP lead to ρkj(t) = σjk(t) exp(−isjkΩt), sjk = sgn(εj − εk), with sjj = 0. The RWA elements are obtained in a steady state limit dσjk/dt = 0, where these elements do not change with time. Their values are given for a system initially at thermal equilibrium by

Section 4 to improve the PBE results. Finally, the Conclusions describe how results from the combined electronic structureRDM treatment provide insight on the studied Si slab dielectric properties and how they compare to bulk Si properties.

2. THEORETICAL TREATMENT OF OPTICAL PROPERTIES A. The Dielectric Function of a Slab. The average polarizability P(t), or induced electric dipole D(t) per unit volume, of a sample in an electric field , (t) points along its direction and is obtained from the retarded response function χ(r)(t′) which accounts for delay effects from time t′ to t, as in7 t

P(t ) =

∫−∞ d′χ (r) (t′) ,(t − t′) +∞

=

∫−∞

dt ′χ

(+)

(t ′) ,(t − t ′)

σjj(ss) = Γ−j 1(σkkeq − σjjeq) ∑ Ωjk 2γjk[γjk 2 + Δjk 2 ]−1 + σjj(eq) k

(1)

σjk(ss)

where χ(+)(t′) = χ(r)(t′) θ(t′) is the dynamical susceptibility in terms of the step function θ(t) and the polarizability is along the field direction. Using a symmetric Fourier transform, the frequency transform of the above functions of time, indicated with a tilde, gives P̃ (ω) = (2π)1/2χ̃(+)(ω) Ẽ (ω) where χ (ω) = (2π)1/2χ̃(+)(ω) is the susceptibility function of frequency ω. In our case, we consider light of frequency Ω and constant electric field amplitude given by , (t) = 2, 0 cos(Ωt) = , 0⌊e−iΩt + eiΩt⌋. The relevant susceptibility function is a complex valued function for ω = Ω, or χ(Ω) = χ′(Ω) + iχ″(Ω), with χ′,χ″ being the real and imaginary parts of the susceptibility function. From the literature,7 the relation between the frequency dependent dielectric function ε(ω) and the dielectric susceptibility χ(ω) is, in the cgs Gaussian system of units, ε(ω) = 1 + 4πχ(ω). Here we omit the dielectric tensor indices and equations apply equally to all three diagonal components. This is also a complex quantity ε(ω) = ε′(ω) + iε″(ω) such that (ε(ω))1/2 = η(ω) + iκ(ω) where η(ω) is the refractive index of the medium and κ(ω) relates to the light absorption coefficient. From the above two equations,26 ε′(ω) = η2(ω) − κ2(ω) and ε″(ω) = 2κ(ω) η(ω) with their inverse relations given by ⎡ η(ω) = ⎢ ⎣

1/2 ε ′ (ω)2 + ε ″ (ω)2 + ε ′ (ω) ⎤

⎡ κ(ω) = ⎢ ⎣

1/2 ε ′ (ω)2 + ε ″ (ω)2 − ε ′ (ω) ⎤

⎥⎦

2

where Γj is a depopulation rate and γjk is a quantum decoherence rate, both resulting from medium phonons and excitons, and available from measurements or from molecular dynamics simulations. At the temperature of 0 K and also with good accuracy at room temperature, σeq jj = 2 for j ≤ HOMO and 4 σeq jj = 0 for j ≥ LUMO. Here Ωjk is the Rabi frequency calculated from the applied electric field amplitude , 0 and transition dipole moments as Ωjk = −Djk·, 0/ℏ. Also Δjk(Ω) = Ω − ωjk is the detuning frequency, using ωjk = (εj − εk)/ℏ for the frequency of the transition, and Ω is the frequency of the light, as before. The transition dipole vector Djk for the j→k transition is constructed from plane wave combinations in the KSOs and given by27 Djk = ⟨j|D̂ |k⟩ = ΣG,G′ 0, from which we extract real and imaginary parts of the dynamical susceptibility, as 4431

dx.doi.org/10.1021/jp410579k | J. Phys. Chem. C 2014, 118, 4429−4436

The Journal of Physical Chemistry C

Article

Figure 1. Electronic density of states for different slab thicknesses with PBE and HSE are given. The band gap is seen to be larger for HSE calculations compared to PBE.

to avoid interactions between periodic images along Z. Each supercell contains N layers each perpendicular to Z and with 4 × 6 Si atoms. The Si(111) surfaces were terminated by 12 hydrogen (H) atoms on each surface to compensate the dangling bonds. In the case of eight layers of c-Si slab, labeled Si192H24, 192 Si atoms were arranged in eight layers with 12 H atoms on lower and upper surfaces of the slab. Our structures have been optimized using DFT with a plane-wave basis set as implemented in the Quantum Espresso28 software package. For the Quantum Espresso calculations, we used norm-conserving pseudopotentials,29 and both PBE (a GGA/DFT functional)15 and the hybrid HSE0617,18 functional with a corrected exchange energy functional at short electronic distances. The c-Si slab structures used here are the ones found in our previous studies [ref 25]. More details about how they are generated can be found in ref 27. The optimized geometries obtained from ab initio calculations were used to compute the Kohn−Sham electronic energy levels and orbitals (KSOs). We have calculated electron density of states (EDOS) with PBE and HSE functionals for different slab thicknesses and they are shown in Figure 1. From the EDOS plot, one can see that the band gap is larger for HSE functionals while density of states are similar. We have used a Lorentzian function with a width of 20 meV for EDOS calculations and details are given in ref 27. We have also compared the charge density of different KSO near the Fermi level generated from PBE and HSE functional and they have been found to be similar. For more information, we have provided in the Supporting Information a graph showing partial charge densities along Z for HOMO and LUMO of c-Si4L using PBE and HSE functional calculated from Quantum Espresso. From the graph, one can see that PBE and HSE partial densities match very well. We have obtained the real and imaginary parts of the dynamic susceptibility and dielectric function from the steady stated RDM, with calculated transition dipoles and our chosen depopulation and decoherence rate coefficients . The transition dipole moment for each transition from VB to CB is calculated from KSOs and its squared magnitude has been averaged over X, Y, and Z directions. Electronic levels in the energy range EHOMO - 4 eV to ELUMO + 4 eV have been included for both

2π Re(χ ̅ (+) ) 1 ∑ |Djk|2 Δjk (Ω)[γjk 2 + Δjk 2(Ω)]−1 = ℏVSC j ≠ k

χ ′(Ω) =

εj > εk

×

(σkkss

− σjjss)

(7)

2π Im(χ ̅ (+) ) 1 ∑ |Djk|2 γjk[γjk 2 + Δjk 2(Ω)]−1 (σkkss − σjjss) = ℏVSC j ≠ k

χ ″(Ω) =

εj > εk

(8)

with the double summations extending over initially occupied orbitals k and unoccupied ones j. In the cgs Gaussian system of units presently used, ε(ω) = 1 + 4πχ(ω), and real and imaginary parts of the dielectric function follow as ε′ = Re(ε) = 1 + 4πχ′ and ε″ = Im(ε) = 4πχ″ . The absorbance is related to the imaginary part of the susceptibility and dielectric function as follows,7,27 α(Ω) = 4πΩVSCχ″(Ω)/[cη(Ω)] = 2VSCΩκ(Ω)/c. In our previous treatment for the absorption of light at semiconductor surfaces, a scaled dynamic absorbance α̅ (Ω) has been calculated from excitation energies and transition dipole moments obtained from ab initio calculations. 27 The absorbance derived from the present linear response and the steady state density matrix treatment has the same form as used in our previous and current treatments, but now with an explicit form for the index of refraction . The real and imaginary parts of the dielectric function are calculated by summing over occupied valence band and over empty conduction band states.

3. DIELECTRIC FUNCTION OF Si SLABS WITH VARYING THICKNESSES We model a c-Si slab with limiting Si(111) surfaces perpendicular to a Z-axis separated by a distance L, with an orthorhombic supercell of lengths a, b, c = L along Z and periodic boundary conditions (PBCs) in X and Y directions. We also use PBCs in the Z direction with added vacuum space 4432

dx.doi.org/10.1021/jp410579k | J. Phys. Chem. C 2014, 118, 4429−4436

The Journal of Physical Chemistry C

Article

The results obtained from the present steady state RDM treatment show good agreement with corresponding ones obtained from the Quantum Espresso software for the complex dielectric function, averaged over diagonal components, with regard to the position and height of the peaks. In this software, the real part is obtained from summations over pair excitations and the imaginary part by using the Kramers−Kronig transformation,33 with an empirical choice of the excitation line width . The comparison of real and imaginary parts of the dielectric function from our steady state RDM for a c-Si12L system using HSE functional, and the Quantum Espresso results with the HSE functional and norm conserving atomic pseudopotentals are given in Figure 3. The shapes of curves are

susceptibility and dielectric function calculations. The used depopulation rate Γ corresponds to a fwhm of 4.13 meV (depopulation time of 1 ps)30,31 and the single decoherence rate γ corresponds to 150 meV14 (decoherence time of 27 fs), taken from the literature for similar systems and employed in our previously published calculations. The real and imaginary parts of the dynamic susceptibility are calculated from eq 7 and eq 8 for 4, 8, and 12 layers of c-Si with PBE and HSE functionals. Similarly, the real and imaginary parts of the dielectric function are calculated for 4, 8, and 12 layers of c-Si with PBE and HSE functionals. The results are shown in Figure 2. These figures can be compared with similar figures for

Figure 2. Comparison of real and imaginary parts of the complex dielectric function computed using PBE and HSE functionals and our steady state RDM treatment for 4, 8, and 12 layers of c-Si, obtained from the average of the three diagonal dielectric tensor components.

Figure 3. Comparison of real and imaginary parts of the dielectric function computed using the HSE functional in Quantum Espresso (dashed lines) and from our steady state RDM treatment (full lines) for a slab with 12 layers of c-Si.

experimental dielectric functions for a cleaved Si(111)(2 × 1) surface32 (Figure 5 from ref 24) and bulk Si24,25 (Figure 3 from ref 25 and Figure 2 from ref 26). Comparison of our results with the above-mentioned experiments and other previous theoretical calculations are given after Table 1.

very similar, but some differences are seen due to our alternative RDM formulation which includes depopulation and decoherence rates and incorporates their physical values in our calculations. The bandgap is decreasing with the number of Si layers in both HSE and PBE calculations as shown in Table 1. This is in good agreement with experimental trends.34 We also calculated the change in the position of the first peak in the imaginary part of the dielectric function ε″(Ω) for the two functionals. The band gap difference and first peak difference are quite similar for the calculated thicknesses, but the difference is larger by about 0.22 eV for the location of the first peak, which depends on the values of transition dipoles. For comparisons, the experimentally measured imaginary part of the dielectric function shows strong peaks around 3.46 and 4.32 eV for both a cleaved Si(111)(2 × 1) surface32 and bulk Si.24,25 In our calculations, we find strong peaks for the imaginary part of the dielectric function at roughly 2.8 and 3.74 eV for the PBE calculation, while our HSE peak positions, at about 3.6 and 4.6 eV, are close to values from measurements. This can be compared to the theoretical calculation of dielectric function for bulk Si using the PBE functional in which E1 and E2 transition are at 2.71 and 3.72 eV.35 The calculations for band gaps of bulk Si using the HSE06 functional16 clearly give values closer to the experimentally measured ones. The real and imaginary parts of the index of refraction have been calculated from our RDM dielectric function using eq 2

Table 1. Band Gap and Difference in the Position of First Peak from Imaginary Part of Dielectric Function for PBE and HSE for Different c-Si slab Thickness Are Given no. of Si layers (N)

band gap difference HSE-PBE [eV]

first peak difference from absorbance HSE-PBE [eV]

band gap using PBE [eV]

band gap using HSE [eV]

4 8 12

0.61 0.59 0.59

0.83 0.81 0.80

1.10 0.91 0.86

1.71 1.50 1.45

The dielectric function calculated using HSE06 and PBE functionals show similar structures of peaks and valleys, while peaks are shifted to higher photon energy in the HSE06 calculation by roughly 0.83 eV in the case of c-Si4L. This difference mainly comes from the difference in the band gap obtained from both functionals, which is 0.61 eV for c-Si4L. There is also a reduction in the peak heights for HSE calculation compared to PBE. The c-Si8L and c-Si12L systems also show similar characteristics using HSE06 functionals compared to PBE. The height of the peaks also increases with increasing slab thickness for both functionals. 4433

dx.doi.org/10.1021/jp410579k | J. Phys. Chem. C 2014, 118, 4429−4436

The Journal of Physical Chemistry C

Article

short-range (SR) and long-range (LR) parts, v̂(PBE) = v̂(SR) x x (μ) + (LR) v̂x (μ), with μ being a parameter specifying the interelectronic short-range. The HSE choice is to modify the SR exchange potential by combining 3/4 of the SR PBE exchange with 1/4 of the SR exchange operator from the Hartree−Fock Hamiltonian operator, given by K̂ (SR) (μ) for each occupied b orbital b. The range parameter μ is 0.207 Å−1, or 0.11 a0−1, in the HSE06 implementation.17,40 Therefore, the HSE energy values can be obtained from the PBE Hamiltonian by adding to ̂ (SR) + it a small perturbation term given by v̂(1) x,corr = −(1/4)∑bKb (SR) (1/4)v̂x . This correction can be used in a perturbation treatment to account for orbital energies ε(HSE) = ε(PBE) + εj′. j j The result is a shift in the excitation energies appearing in the − ε(HSE) = ε(PBE) − transition dipole calculations so that ε(HSE) k j k (PBE) εj + Δkj. Furthermore, the susceptibility functions are obtained as sums over initially occupied orbitals j and unoccupied orbitals k so that the overall result involves many energy shifts in denominators weighted by the magnitude of the (jk) momentum matrix elements. To describe the changes in position and intensity of the peaks in the optical properties while using HSE functional compared to PBE, we shift the transition energy in the transition dipole moment calculation by a factor called Δshift as in

for 4, 8, and 12 layers of c-Si with output from the HSE functional. The results are shown in Figure 4 and have been

Figure 4. Comparison of real and imaginary parts of the refractive index computed with the HSE functional using dielectric functions from Quantum Espresso (dashed lines) and from our steady state RDM treatment (full lines) for 4, 8, and 12 layers of c-Si.

Djk = ⟨j|D̂ |k⟩ = ⟨j|p̂ |k⟩iℏe/[me(εk − εj + Δshift )]

obtained from the averaged dielectric functions. The calculated real and imaginary parts η(Ω) and κ(Ω) of the refractive index show good agreement with values given in the literature for bulk Si.36 For example, the values for η(Ω) at photon energies 0.004 and 3.0 eV are 3.42 and 5.22, respectively, for bulk Si in the tabulation of ref 36. The η(Ω) derived from our RDM dielectric functions for the slab Si12L at the above-mentioned photon energies are 3.32 and 5.98 eV respectively. The complex index of refraction derived from our RDM dielectric function is in better agreement with experiment than the ones obtained from the Quantum Espresso dielectric functions. Both the index of refraction and the absorption increase with slab thickness. The absorbance density α/VSC could also be easily calculated from the generated dielectric functions for the three systems with 4, 8, and 12 layers. It is now possible to obtain the frequency dependence of these functions as they change with slab thickness and in particular to incorporate the changes of the refractive index with frequency, which is usually taken to be nearly constant in the literature. There have been attempts in the literature to correct the band gaps for insulators and semiconductors obtained from DFT within the LDA.37,38 In those treatments, the DFT band gap is corrected with an energy term Δ which comes from the discontinuity of the exchange-correlation energy. As a further (computationally expensive) improvement, one could use the GW method derived from many-body perturbation theory.39 We take a simpler approach next, to correct the PBE results for the dielectric function, using a single energy shift parameter.

(9)

and we use in what follows two different choices of the shift: the first one ΔBG correcting for the difference in band gap from the two DFT calculations, or Δpk as the second one correcting a calculated (or otherwise measured) first peak position in the absorption intensity. We replaced the band gap for PBE with the band gap obtained from HSE by ΔBG = 0.61 eV for c-Si4L as obtained in the present calculations. The above-mentioned differences in the susceptibility, and dielectric functions, are almost reproduced by the shift in the band gap. The results for susceptibility with PBE, shifted PBE with ΔBG and HSE are shown for a c-Si4L system in Figure 5. A more accurate reproduction of HSE optical properties can be obtained from PBE using the difference in the position of the first peak in the imaginary part of the dielectric function.

4. ENERGY SHIFTED RESULTS FOR OPTICAL PROPERTIES The semilocal PBE Hamiltonian has the form h(̂ PBE) = k1̂ + v̂(1) el (PBE) v̂(1) , in terms of the one-electron kinetic energy, the ion + v̂xc Coulomb energy from other electrons and from ion cores, and the exchange correlation final term, which is a function of the electron density and its gradient. In the HSE formulation, the exchange correlation term is given as a sum of exchange and correlation terms and the exchange potential is separated into

Figure 5. Both real and imaginary parts of the susceptibility with PBE, HSE, and PBE with a band gap shift are given for a c-Si4L system. The shift in the PBE band gap almost reproduced the HSE results for the susceptibilities. 4434

dx.doi.org/10.1021/jp410579k | J. Phys. Chem. C 2014, 118, 4429−4436

The Journal of Physical Chemistry C

Article

calculations show that the magnitude of light absorption and refraction index increase with increasing slab thickness. All these observations are relevant to the development of Si based thin film materials for photovoltaic applications. In the numerical implementation, PBE and HSE functionals have been found to give similarly shaped DOS and optical spectrum, but values from HSE suggest better agreement with experiments on bulk Si. The real and imaginary parts of the dielectric function versus photon energies obtained from steady state expressions for the RDM show very good agreement with the shapes of functions obtained from the Quantum Espresso software but values of the refractive index are somewhat different at low photon energies. Furthermore, the real and imaginary parts of the complex index of refraction calculated from our RDM dielectric function gave very good agreements with experimental values for bulk Si in the literature. Calculations of the absorbance density versus photon energies (or light frequency) for the three slabs of interest here could be done in detail from information in this paper, incorporating the presently calculated frequency dependence of the index of refraction, which changes with slab thickness. A straightforward procedure has been introduced to correct the underestimation of the optical absorption by the PBE/GGA functional, compared to the more accurate but more computationally expensive HSE calculation, by using a shift parameter. This Δshift parameter can be obtained from theoretical calculation or from experiments. The shift in the absorption peaks and the decrease in intensity from the HSE calculations compared to PBE is mainly due to the higher band gap obtained in the HSE calculation. The Δshift parameter can be chosen to equal the difference in band gap, or the difference in positions of the lowest energy peaks, accurately calculated or obtained from experimental measurements of absorption peak positions. One could also consider generating the energy shift from other DFT calculations, such as the recently proposed and computationally less expensive M11-L functional instead of the HSE functional, that provides band gaps with accuracy close to HSE.41,42 The results presented in this paper relate to the Γ point of the Brillouin zone for the present orthorhombic supercells, corresponding to an electronic wavevector k = 0, at the origin of the electronic energy bands vs wavevectors. We have also calculated the same optical properties for other wavevectors k and have found that optical properties differ in detail but are qualitatively similar. Averaged optical properties obtained from averages over the electronic wavevectors are expected to give results similar to the ones presented here.

The values given in Table 1 for different Si slab thicknesses are used here for the shift Δpk . The calculated real and imaginary parts of the dynamical susceptibility nearly reproduced the HSE results from PBE. The results are shown in Figure 6. Similar

Figure 6. Both real and imaginary parts of the susceptibility with PBE, HSE, and PBE with peak shift are given for 4, 8, and 12 layers of c-Si. The shift in the PBE transition energy nearly reproduced HSE susceptibility values.

improved results are found for dielectric function and absorption coefficient calculations. The better improvement by using the shift Δpk instead of ΔBG comes from the fact that the largest transition dipole moments are in general for orbitals from below the HOMO to above the LUMO levels and their energy difference are larger than the band gap difference. The shapes of orbitals obtained in PBE and HSE calculations (not shown here but included in Supportive Information) are very similar and consequently matrix elements of the electron momentum operator are also very similar. This may explain why the simple change of all excitation energies by a single energy shift gives quite good agreement with detailed HSE results. Corrections such as we have described here would require doing calculations with the more computational expensive HSE functional. However, from Table 1, one can see that the relative change in the band gap difference and in the first peak difference are only within a few percentage going from 4L to 12L for both PBE and HSE functionals. In practice, this means that there is a need to perform only one HSE calculation for a smaller system to estimate the values for Δpk and ΔBG from it for the larger systems.

■

5. CONCLUSIONS In this paper, we have derived a reduced density matrix treatment of the dynamical polarizability of Si slabs to incorporate dissipative phenomena due to excitations in the medium, and to obtain the dynamical susceptibility, dielectric function, and light absorption coefficient. This treatment incorporates depopulation and decoherence rates into the equations for the susceptibilities, which can be physically identified and separately calculated when the atomic structure of the systems is varied, or taken from experiments. From our calculations here, it has been found that Si slabs of varying thickness show the same peak structure in absorption and refraction functions of photon energy as in bulk Si, and have peak positions in agreement with those in bulk Si. In addition,

ASSOCIATED CONTENT

S Supporting Information *

Graph showing partial charge densities from HOMO and LUMO Kohn−Sham orbitals from PBE and HSE functionals for c-Si4L (a Si four layer slab). This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*Mailing address: 2318 NPBldg, POB 118435, University of Florida, Gainesville, Florida 32611. Tel. 352-392-6977. E-mail: [email protected]fl.edu. Notes

The authors declare no competing financial interest. 4435

dx.doi.org/10.1021/jp410579k | J. Phys. Chem. C 2014, 118, 4429−4436

The Journal of Physical Chemistry C

■

Article

(24) Aspnes, D. E.; Studna, A. A. Dielectric Functions and Optical Parameters of Si, Ge, GaP, GaAs, GaSb, InP, InAs, and InSb from 1.5 to 6.0 eV. Phys. Rev. B 1983, 27, 985−1009. (25) Lautenschlager, P.; Garriga, M.; Vina, L.; Cardona, M. Temperature Dependence of the Dielectric Function and Interband Critical Points in Silicon. Phys. Rev. B 1987, 36, 4821−4830. (26) Landau, L. D.; Lifshitz, E. M. Electrodynamics of Continuous Media; Pergamon Press: Oxford, 1960; Chapter IX, Vol. 8. (27) Vazhappilly, T.; Kilin, D. S.; Micha, D. A. Photoabsorbance and Photovoltage of Crystalline and Amorphous Silicon Slabs with Silver Adsorbates. J. Phys. Chem. C 2012, 116, 25525−25536. (28) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; et al. Quantum Espresso: A Modular and Open-Source Software Project for Quantum Simulations of Materials. J. Phys.: Condens. Matter 2009, 21, 395502. (29) Hamann, D. R.; Schlüter, M.; Chiang, C. Norm-Conserving Pseudopotentials. Phys. Rev. Lett. 1979, 43, 1494−1497. (30) Kilin, D. S.; Micha, D. A. Relaxation of Photoexcited Electrons at a Nanostructured Si(111) Surface. J. Phys. Chem. Lett. 2010, 1, 1073−1077. (31) Weinelt, M.; Kutschera, M.; Fauster, T.; Rohlfing, M. Dynamics of Exciton Formation at the Si(100) C(4 × 2) Surface. Phys. Rev. Lett. 2004, 92, 126801. (32) Kelly, M. K.; Zollner, S.; Cardon, M. Modelling the Optical Response of Surfaces Measured by Spectroscopic Ellipsometry: Application to Si and Ge. Surf. Sci. 1993, 285, 282−294. (33) Colle, R.; Parruccini, P.; Benassi, A.; Cavazzoni, C. Optical Properties of Emeraldine Salt Polymers from Ab Initio Calculations: Comparison with Recent Experimental Data. J. Phys. Chem. B 2007, 111, 2800−2805. (34) Bassani, F.; Mihalcescu, I.; Vial, J. C.; Arnaud d’Avitaya, F. Optical Absorption Evidence of Quantum Confinement in Si/Caf2 Multilayers Grown by Molecular Beam Epitaxy. Appl. Surf. Sci. 1997, 117−118, 670−676. (35) Gajdoš, M.; Hummer, K.; Kresse, G.; Furthmüller, J.; Bechstedt, F. Linear Optical Properties in the Projector-Augmented Wave Methodology. Phys. Rev. B 2006, 73, 045112. (36) Handbook of Optical Constants of Solids; Palik, E. D., Ed.; Academic Press: New York, 1998. (37) Sham, L. J.; Schlüter, M. Density-Functional Theory of the Energy Gap. Phys. Rev. Lett. 1983, 51, 1888−1891. (38) Lannoo, M.; Schlüter, M.; Sham, L. J. Calculation of the KohnSham Potential and Its Discontinuity for a Model-Semiconductor. Phys. Rev. B 1985, 32, 3890−3899. (39) Hybertsen, M. S.; Louie, S. G. Electron Correlation in Semiconductors and Insulators: Band Gaps and Quasiparticle Energies. Phys. Rev. B 1986, 34, 5390−5413. (40) Paier, J.; Marsman, M.; Hummer, K.; Kresse, G.; Gerber, I. C.; Á ngyán, J. G. Erratum: “Screened Hybrid Density Functionals Applied to Solids” [J. Chem. Phys.124, 154709 (2006)]. J. Chem. Phys. 2006, 125, 249901. (41) Peverati, R.; Truhlar, D. G. M11-L: A Local Density Functional That Provides Improved Accuracy for Electronic Structure Calculations in Chemistry and Physics. J. Phys. Chem. Lett. 2011, 3, 117− 124. (42) Peverati, R.; Truhlar, D. G. Performance of the M11-L Density Functional for Bandgaps and Lattice Constants of Unary and Binary Semiconductors. J. Chem. Phys. 2012, 136, 134704.

ACKNOWLEDGMENTS This work has been partly supported by the National Science Foundation grant NSF CHE 1011967. Computing support was provided by the High Performance Computing facility of the University of Florida.

■

REFERENCES

(1) Brillson, L. J. Surfaces and Interfaces of Electronic Materials; WileyVCH: Weinheim, 2010; Chapters 10 and 15. (2) Yu, P.; Cardona, M. Fundamentals of Semiconductors: Physics and Materials Properties, 3rd ed.; Springer-Verlag: Berlin, 2005. (3) Kronik, L.; Shapira, Y. Surface Photovoltage Phenomena: Theory, Experiment, and Applications. Surf. Sci. Rep. 1999, 37, 1−206. (4) Kilin, D. S.; Micha, D. A. Surface Photovoltage at Nanostructures on Si Surfaces: Ab Initio Results. J. Phys. Chem. C 2009, 113, 3530− 3542. (5) Vazhappilly, T.; Micha, D. A. Atomic Modeling of Structural and Optical Properties of Amorphous Silicon. Chem. Phys. Lett. 2013, 570, 95−99. (6) Atkins, P. W.; Friedman, R. S. Molecular Quantum Mechanics; Oxford University Press: New York, 1997; Chapter 6. (7) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976; Chapter 21. (8) May, V.; Kühn, O. Charge and Energy Transfer Dynamics in Molecular Systems: Wiley-VCH: Berlin, 2000; Chapter 3. (9) Nitzan, A. Chemical Dynamics in Condensed Phases; Oxford University Press: New York, 2006; Chapter 18. (10) Martin, R. M. Electronic Structure: Basic Theory and Practical Methods; Cambridge University Press: Cambridge, New York, 2004; Chapters 12 and 20. (11) Nozik, A. J. Spectroscopy and Hot Electron Relaxation Dynamics in Semiconductor Quantum Wells and Quantum Dots. Annu. Rev. Phys. Chem. 2001, 52, 193−231. (12) Patel, B. K.; Rath, S.; Sahu, S. N. Surface Photovoltage Studies of Si Nanocrystallites Prepared by Electrochemical Etching. Phys. E 2006, 33, 268−272. (13) Bylaska, E. J.; Tsemekhman, K.; Gao, F. Phys. Scr. 2006, T124, 86−90. (14) Rohlfing, M.; Louie, S. G. Electron-Hole Excitations and Optical Spectra from First Principles. Phys. Rev. B 2000, 62, 4927−4944. (15) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (16) Heyd, J.; Scuseria, G. E. Efficient Hybrid Density Functional Calculations in Solids: Assessment of the HeydScuseriaErnzerhof Screened Coulomb Hybrid Functional. J. Chem. Phys. 2004, 121, 1187−1192. (17) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Erratum: ″Hybrid Functionals Based on a Screened Coulomb Potential” [J. Chem. Phys. 118, 8207 (2003)]. J. Chem. Phys. 2006, 124, 219906. (18) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid Functionals Based on a Screened Coulomb Potential. J. Chem. Phys. 2003, 118, 8207−8215. (19) Moussa, J. E.; Schultz, P. A.; Chelikowsky, J. R. Analysis of the Heyd-Scuseria-Ernzerhof Density Functional Parameter Space. J. Chem. Phys. 2012, 136, 204117−204110. (20) Hedin, L. New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem. Phys. Rev. 1965, 139, A796−A823. (21) Onida, G.; Reining, L.; Rubio, A. Electronic Excitations: Density-Functional Versus Many-Body Green’s-Function Approaches. Rev. Mod. Phys. 2002, 74, 601−659. (22) Jain, M.; Chelikowsky, J. R.; Louie, S. G. Reliability of Hybrid Functionals in Predicting Band Gaps. Phys. Rev. Lett. 2011, 107, 216806. (23) Paier, J.; Marsman, M.; Kresse, G. Dielectric Properties and Excitons for Extended Systems from Hybrid Functionals. Phys. Rev. B 2008, 78, 121201. 4436

dx.doi.org/10.1021/jp410579k | J. Phys. Chem. C 2014, 118, 4429−4436