Use of energy bandgap adjustment when simulating optical properties to maintain correct energy dissipation: Application to Ruddlesden–Popper oxysulfides

Calculations of the optoelectronic properties of materials used in photo-electrochemical conversion systems and solar cells are frequently necessary, especially in the case that these properties are not accessible by experimental means. One common example is that of particulate semiconductors for photo-electrochemical water splitting, for which particle sizes are typically too small to allow absorption coefficients to be determined experimentally.

Optical properties can be computed from frequency dependent dielectric coefficients using density functional theory (DFT).¹ In fact, DFT has been widely adopted as a means of evaluating ground state electron densities and accessing ground state properties such as binding energies, elastic constants, and forces required for molecular dynamics simulations.² When employing a single-particle model, DFT can also be employed to calculate one-body electron densities and the corresponding eigenenergies, which can then be used to obtain electronic energy levels including those of excited states.³ However, one specific limitation associated with DFT predictions of excited states is the underestimated energy bandgap. This limitation can be overcome by using many-body perturbation theory, otherwise termed the GW approximation (GWA) based on Green’s function. Even so, the GWA is costly and so simpler approximations, such as the local density approximation (LDA), are more widely employed. The use of a scissor operation to shift the conduction bands has been introduced based on the observation that the dispersion of the conduction and valence bands as computed with the LDA is very similar to that obtained with the GW, although the energy bandgap is considerably shifted.^4–8 The scissor operators were initially introduced to calculate band structures and static dielectric constant.^4–8 Optical absorption spectra have been calculated using the scissor operator approach to correct Kohn–Sham energies or local-density approximation (LDA).^9–14 One practical approach to calculating the dispersion of the conduction and valence bands while obtaining accurate energy bandgap values is to use the hybrid exchange–correlation functional (HSE06), although the computational cost is high. Note that, in the case of such DFT/HSE06 calculations, the Hartree–Fock (HF) exchange correlation energy is mixed with a parameter.

Recently, absorption coefficient spectra have been accurately calculated using the semilocal generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE) method together with a high-density k-point mesh, albeit with an underestimated energy bandgap. In this prior work, the imaginary part of the dielectric coefficient was blueshifted to compensate for the underestimated energy bandgap and the amplitude was determined so as to satisfy the moment sum rule. We refer to this method as the PBE + HSE06 + sum rule (PHS).^14–16 These previous studies showed that absorption spectra computed by a cost-effective DFT technique accurately reproduce the experimental absorption spectra of various photovoltaic materials.

Although the PHS method computed the absorption spectra accurately for seven solar cell materials (GaAs, InP, CdTe, CuInSe₂, Cu₂ZnGeSe₄, chalcogenide perovskites, and MgSiN₂), the physical basis for this approach was not clear.^15–19 The PHS method was originally developed to satisfy the first moment sum rule, $∫−∞∞dEEϵi(E)=constant$⁠, where ϵ_i(E) is the imaginary part of the dielectric coefficient. It should be noted that, in such cases, the energy scale is shifted to compensate for the difference between the true energy bandgap and the predicted energy bandgap. Strictly speaking, other forms of ϵ_i(E) might satisfy the first moment sum rule because integration is performed with respect to E. In addition, the first moment sum rule is not sufficient to justify the scaling of ϵ_i(E).

The present work provides an alternative justification for scaling the amplitude of the imaginary part of the dielectric coefficient after introducing a blueshift. In particular, instead of imposing the sum rule, we impose the condition that the power loss and/or the optical conductivity is independent of the energy shift. Both the power loss and the optical conductivity result from scattering of carriers, which in turn is affected by the dispersion of the conduction or valence band but could be insensitive to the energy bandgap. By taking into account the energy bandgap shift associated with inter-band transitions, the power loss and optical conductivity can be evaluated given that the dispersion of the conduction and valence bands can be accurately calculated using DFT. Here, we show that scaling the imaginary part of the dielectric coefficient is based on the fact that the power loss and optical conductivity are dictated by the energy dispersion, and the energy bandgap shift associated with inter-band transitions is considered independently. Other optical properties, such as the real part of the dielectric coefficient, can be obtained by using the Kramers–Kronig transformation together with additional well-known relationships among the optical coefficients.^15,16 In our work, the PHS method is introduced to maintain correct energy dissipation or the electrical conductivity while the energy bandgap is shifted. Although the PHS method is originally proposed for the PBE method against the HSE06 method, the method is generally applicable for adjusting the energy bandgap of DFT approach as long as the dispersion of the conduction and valence bands is accurately predicted.

The proposed approach is validated herein using the perovskite semiconductor Y₂Ti₂O₅S₂ as a model. This compound has recently received attention with regard to possible applications in visible-light-driven overall water splitting, although the relatively low absorption coefficient for this material in the vicinity of the energy bandgap could limit its performance.^20–26 Y₂Ti₂O₅S₂ has also attracted interest as an anode material for Li-ion batteries because it exhibits rapid lithium ion diffusion in the Ruddlesden–Popper phase.^27,28 The optical properties following the scaling method described above, calculated using DFT/PBE, are close to those determined by DFT/HSE06, although the latter has a high computational cost. The optical properties following scaling have already been shown to be equivalent to those of certain solar cell materials such as GaAs, InP, CdTe, CuInSe₂, Cu₂ZnGeSe₄, and chalcogenide perovskites.^15–19 However, this method has not yet been applied to Ruddlesden–Popper oxysulfides (Y₂Ti₂O₅S₂). The optical properties of the other Ruddlesden–Popper oxysulfide (Gd₂Ti₂O₅S₂) have been calculated using the above-mentioned method based on DFT/PBE.²⁹ 

For simplicity, tensor notation is not used to specify the Cartesian coordinates herein unless otherwise stated. Note that this tensor notation can be easily recovered. As an example, the electrical conductivity, σ, should be written as σ_α,β using Cartesian coordinates denoted by α = x, y, z and β = x, y, z, where x, y, and z indicate the three axes in the Cartesian coordinate system.

When using DFT, the energy bandgap is often underestimated. DFT can be used to calculate the imaginary part of the dielectric coefficient at a given energy, E, defined as $ϵi(bg)E$⁠, where the superscript (bg) indicates that the energy bandgap associated with E_bg is obtained using DFT and thus may differ from the true value E_BG. In such cases, the energy difference is denoted by ΔE = E_BG − E_bg. Herein, we denote the correct imaginary part of the dielectric coefficient at energy E as $ϵi(BG)E$⁠, where the energy bandgap is E_BG. In some instances, $ϵi(BG)E$ can be obtained either experimentally or using a high-cost DFT method while employing, for example, the hybrid HSE06 functional. However, accurate determination of energy bandgap values for compound semiconductors composed of more than three elements based on DFT techniques remains challenging. The following discussion demonstrates a recently developed empirical method that scales the magnitude of the imaginary part of the dielectric coefficient obtained using DFT, for which the energy bandgap is shifted. This approach can be justified by considering that the energy loss and the electrical conductivity both result from scattering, which in turn is affected by energy dispersion. In addition, the energy associated with inter-band transitions is modified according to the energy bandgap difference.

The relationship given by Eq. (5) can also be inferred from the photoconductivity. For metals, the imaginary part of the dielectric coefficient is related to the electrical conductivity, σ, according to the equation $ϵ(BG)E=1+4πiσ(ω)/ω$⁠.

Based on the fluctuation–dissipation (FD) theorem, the electrical conductivity can be expressed by the correlation function for the random part of the current that results from carrier scattering. The electrical conductivity is affected by the dispersion of the conduction band or valence band but, again, can be insensitive to the energy bandgap.

At 0 K, $fFDEn,k⃗$ in the numerator of Eq. (16) will be non-zero only for the Fermi–Dirac distribution at the valence band (occupied states); this is not affected by the energy bandgap shift in the DFT results. Here, we should also note that $fFDℏω+En,k⃗$ is zero for ℏω > E_BG, where we have used $En′,k⃗=ℏω+En,k⃗$ according to the delta function in Eq. (15). Even though the energy bandgap is shifted in DFT, the numerator of Eq. (16) is not influenced by the energy bandgap shift. When the energy bandgap is shifted, $En′,k⃗(bg)−En,k⃗(bg)$ in Eqs. (17) and (18) can be regarded as $En′,k⃗−En,k⃗$⁠, where we have $En′,k⃗(bg)−En,k⃗(bg)=En′,k⃗−En,k⃗−ΔE$⁠. By introducing such a replacement with ω′ and ω satisfying ℏω′ = E − ΔE and ℏω = E, we find $σr(inter,bg)(ω′)=σr(inter)(ω)$⁠, which is nothing but Eq. (9). On this basis, Eqs. (9) and (5) can be reproduced using $ϵ(BG)E=1+4πiσr(ω)/ω$⁠. In the Bethe–Salpeter equation (BSE) approach to many-body perturbation theory, Eq. (13) represents a so-called diagonal term for which off-diagonal electron–hole exchange interactions are neglected,^35,36 which might represent a limitation of the correction method considered herein. As long as the e–h interaction is taken into account within quasiparticle approximations, the energy bandgap correction by the PHS method can be applicable.

The PHS method can be regarded as a kind of the scissor operator approach; in the scissor operator approach, the energy difference between conduction band states and valence band states is shifted.¹³ In the PHS method considered here, the energy shift is taken into account in the expression of the electrical conductivity and ω is shifted accordingly. We considered the electrical conductivity because the electrical conductivity is directly related to energy dissipation caused by scatterings, as can be understood by comparing the left-hand side of Eq. (7) with the right-hand side of Eq. (1). Here, we assumed that energy dissipation caused by scattering could be accurately calculated as long as the band structures were maintained, although the energy bandgap was shifted. The imaginary part of the dielectric coefficient is obtained using the scaling relation given by Eq. (5) rather than applying the scissor operator approach directly to the imaginary part of the dielectric coefficient. The real part of the dielectric coefficient can be calculated from the imaginary part of the dielectric coefficient by applying the Kramers–Kronig transformation. Note that a part of the energy difference should not be modified when the scissor operator is used directly to the imaginary part of the dielectric coefficient.¹² Here, the scaling relation given by Eq. (5) is introduced to maintain correct energy dissipation or the electrical conductivity while the energy bandgap is shifted. Equation (5) was previously considered to ensure that the same momentum sum rule be fulfilled as that of an unshifted dielectric coefficient.^14–16 

Another method sometimes used to correct the energy bandgaps is the Tran–Blaha modified Becke–Johnson (TB-mBJ) potential.^37–42 The exchange and correlation effects were considered self-consistently by the TB-mBJ potential. The TB-mBJ potential contains parameters fitted to reduce the error of the energy bandgaps compared to the experimental values for a set of materials.^37–42 Although the energy bandgaps are accurately calculated using the TB-mBJ potential, electron mass for the conduction band becomes heavier for some materials.^38–40,42 In the PHS method, the functionals are not modified. The PHS method is applied when calculating the imaginary part of the dielectric coefficient. We have shown that the PHS method is applicable when the power loss (optical conductivity) is independent of the energy bandgap shift associated with inter-band transitions.

Even though the energy value will be shifted in the DFT results, the energy bandgap shift associated with inter-band transitions can be taken into account simply by modifying the corresponding energy in the expressions for the power loss and transport coefficients, such as the electrical conductivity. This will be the case assuming that the band structure has been accurately calculated and the scattering properties of charge carriers are unaltered.

Structural optimizations together with band dispersion and dielectric function calculations were performed using the VASP software.⁴⁶ Y₂Ti₂O₅S₂ and Gd₂Ti₂O₅S₂ were considered, each with tetragonal crystal symmetry and space group I4/mmm. The projected augmented wave pseudopotential with the exchange⁴⁷ and correlation functional under the generalized gradient approximation in the Perdew–Burke–Ernzerhof (PBE) form was applied.⁴⁸ The Brillouin zones were sampled from the Γ − centers of 8 × 8 × 1 k-point grids for Y₂Ti₂O₅S₂ in the conventional cell. For Gd₂Ti₂O₅S₂, calculations were performed using the primitive cell, which contains half the number of atoms compared to the conventional cell, to ensure the convergence of the dielectric function with a large number of k-points, and the Brillouin zone was sampled from 8 × 8 × 12 k-point grids. A plane wave basis set with a kinetic energy cutoff of 520 eV was used, and the electronic and force convergence criteria during the structural optimization were set to 10⁻⁸ eV and 10⁻² eV/Å, respectively. The optimized lattice parameters of Y₂Ti₂O₅S₂ were a = b = 3.79 Å and c = 23.0 Å. For Gd₂Ti₂O₅S₂, the optimized lattice parameters in the conventional cell were a = b = 3.81 Å and c = 23.1 Å. The Heyd–Scuseria–Ernzerhof hybrid functional (HSE06 method)⁴⁹ together with a Hartree–Fock screening parameter of 0.2 was employed when computing the dielectric function in the independent-particle representation.

Recently, the optical properties of Y₂Ti₂O₅S₂ have been computed using the DFT/HSE06 method, and the predicted reflectivity data were found to agree with experimental spectra.²³ The reflectivity data were used to show that the wavelength at the maximum of diffuse-reflectance spectrum measured for Y₂Ti₂O₅S₂ particles is correlated with the wavelength where the light penetration depth becomes comparable to the thickness of the particle.²⁵ An energy bandgap of 1.91 eV was obtained from the DFT/HSE06 approach, which was also consistent with the experimental value. Meanwhile, a recent study reported that the energy bandgap of Y₂Ti₂O₅S₂ calculated using the quasiparticle self-consistent GW (QSGW) method is 3.05 eV,⁵⁰ which is a significant overestimation compared to the fundamental energy bandgap obtained using the HSE06 functional. By considering the difference of over 1 eV, we adopted the DFT/HSE06 approach for comparison with the scaling approach (PHS method) using the DFT/PBE method. Although the DFT/HSE06 method is known to accurately predict band structures, the computational cost associated with this technique is high. As an alternative, the present work used the DFT/PBE method and Eq. (5), with $ϵi(bg)E$ calculated via the DFT/PBE technique. This approach gave an estimated energy bandgap of 0.78 eV.

Figure 1 demonstrates that the imaginary part of the average dielectric coefficient of Y₂Ti₂O₅S₂ for isotropic incident light as calculated by the PBE method differed from that obtained using the DFT/HSE06 approach. This discrepancy corresponded to the large energy band difference ΔE = 1.91 − 0.78 = 1.13 eV. Following the energy shift $ϵi(bg)E−ΔE$⁠, the dispersion of the imaginary part of the dielectric coefficient calculated by the DFT/PBE method was similar to that calculated by the DFT/HSE06 method, but the magnitude was inaccurate. By applying the transformation given by Eq. (5), the imaginary parts of the average dielectric coefficients as obtained from the DFT/HSE06 and DFT/PBE approaches were almost identical, especially below 3.5 eV.

Figure 2 presents the additional average optical properties of Y₂Ti₂O₅S₂ for isotropic incident light, including absorption coefficient and reflectivity data, as calculated using both the DFT/HSE06 and DFT/PBE techniques. The results for the imaginary part of the average dielectric coefficient are plotted on a semi-log scale in Fig. 2(a). Note that the DFT/HSE06 values are reproduced by the scaled results for the DFT/PBE method over a wide range. By applying the Kramers–Kronig transformation for the scaled imaginary part of the dielectric coefficient based on Eq. (20), where $ϵi(BG)E$ is obtained from Eq. (5), the real part of the dielectric coefficient could be calculated using the DFT/PBE method. When performing the Kramers–Kronig transformation, the minimum and maximum values of E were introduced. The imaginary part of the dielectric coefficient was evidently close to zero below the energy bandgap and thus would be expected to be insensitive to the minimum cutoff energy value in the case that this minimum is set below the energy bandgap. Note also that the minimum value of E was set to ΔE, which was lower than the energy bandgap. The maximum energy cutoff was set sufficiently high so as to affect the imaginary part of the dielectric coefficient by at most 3% over the range of 20–88 eV. Compared to the value that was simply shifted, $ϵr(bg)E−ΔE$⁠, the scaled value was close to that calculated by the DFT/HSE06 method, and, in fact, a remarkable agreement between the two was obtained below 3.5 eV. The refractive index and the extinction coefficient were calculated based on Eq. (21) using $[(E−ΔE)/E]ϵi(bg)E−ΔE$ and applying the Kramers–Kronig transformation. Subsequently, the absorption coefficient and the reflectivity were obtained using Eqs. (23) and (22), respectively. Because both the imaginary and real parts of the dielectric coefficient were accurately determined, the optical properties obtained from these values also matched the results of the DFT/HSE06 method. The absorption coefficient, α^(bg) (E − ΔE) [where α^(bg) (E) was calculated by the DFT/PBE method], matched that calculated by the DFT/HSE06 method. This outcome indicates that the absorption coefficient can be estimated simply by adjusting the energy to compensate for the difference between the calculated energy bandgap values. The simple energy bandgap shift of absorption coefficient reflects the fact that the real part of the dielectric coefficient shown in Fig. 2(b) can be approximated as a constant and is larger than the imaginary part of the dielectric coefficient for the energy below 6 eV. Under the condition of ϵ_r being larger than ϵ_i, we obtain α(E) = 2Ek_e(E)/(cℏ) with $ke≈ϵi/2ϵr$⁠, which yields $α(E)≈Eϵi/cℏϵr$⁠.¹⁴ When $ϵr$ can be approximated as a constant, we approximately obtain α(E) ∝ Eϵ_i apart from a proportionality constant, which is independent of E. Therefore, we have α(E) ∝ σ_r(E) and α(E) can be shifted without scaling such as the electrical conductivity.

Figure 3 presents the average optical properties of Gd₂Ti₂O₅S₂ for isotropic incident light as calculated using both the DFT/HSE06 and DFT/PBE techniques. Like the cases of Y₂Ti₂O₅S₂, the optical properties of the scaled values were overall close to those calculated by the DFT/HSE06 method. By looking closely, however, this figure shows a discrepancy in the imaginary part of the dielectric constant and in the absorption coefficient around 2 eV, although the magnitude is small; the agreement is satisfactory in the energy range of 2.2–6 eV. The value that was simply shifted, $ϵr(bg)E−ΔE$⁠, significantly overestimates the value of the real part of the dielectric constant calculated by the DFT/HSE06 method for the energy below 5 eV. The simple energy bandgap shift of the absorption coefficient matched that calculated by the DFT/HSE06 method as in the case of Y₂Ti₂O₅S₂. As already explained, the energy bandgap shift of the absorption coefficient can be applicable when the real part of the dielectric coefficient can be approximated as a constant and is larger than the imaginary part of the dielectric coefficient. The optical properties other than the absorption coefficient are overestimated if the energy is just shifted from E to E − ΔE. The results shown in Fig. 3 are consistent with those shown in Fig. 2 and thus indicate the validity of the concept of scaling, which is realized by the PHS method, for the Ruddlesden–Popper oxysulfides.

Power loss corresponds to energy dissipation associated with the dispersion of the conduction and valence bands. The imaginary part of the dielectric coefficient is related to both the power loss and the electrical conductivity; the latter expresses dissipation associated with electrical current fluctuations based on the FD theorem. Both the power loss and the electrical current can be determined once the dispersion of the conduction and valence bands is accurately obtained, and the energy can be independently adjusted according to the shift of the energy bandgap. By considering that scattering processes are unaffected by the energy bandgap and that the energy of inter-band transitions is shifted to adjust the energy bandgap, the imaginary part of the dielectric coefficient as calculated by DFT can be scaled. Although the concept of scaling has been introduced previously, justification for this method based on an analysis of the power loss and an explicit expression for conductivity derived from the FD theorem in an independent-particle model has not yet been presented.

We have shown that the PHS method is a kind of the scissor operator approach, where the energy difference between conduction band states and valence band states is shifted in such a way that the energy loss expressed by the imaginary part of the dielectric coefficient multiplied by the energy is maintained after the energy bandgap shift. The optical conductivity is related to the energy loss (dissipation); the PHS method corresponds to the energy bandgap shift operation in the expression of the optical conductivity.

The validity of the present scaling was confirmed by employing this method to calculate the optical properties of Y₂Ti₂O₅S₂. This material has attracted attention because of its activity for the overall water splitting reaction under visible light as well as its potential usage as an anode material in Li-ion batteries.^20–28 Previous work showed excellent agreement between reflectivity values computed using the DFT/HSE06 and experimental spectra.²³ The study reported herein established that the optical properties of Y₂Ti₂O₅S₂ can be accurately predicted using the DFT/PBE approach based on scaling the imaginary part of the dielectric coefficient to conserve electrical conductivity after shifting the energy bandgap.

It should be noted that the DFT/HSE06 method is not universally applicable. As an example, the HSE06 functional has recently been found to overestimate the energy bandgaps of certain vanadium oxides by more than 1 eV.⁵¹ Theoretically, the applicability of the present scaling relationship relies on the accuracy with which the dispersion of the conduction and valence bands is predicted but is not affected by the direction of the energy bandgap shift. Assuming that these dispersions are correctly calculated, the optical properties can be obtained via scaling of the imaginary part of the calculated dielectric coefficient based on the negative shift, with ΔE < 0.

This work was financially supported by the Artificial Photosynthesis Project of the New Energy and Industrial Technology Development Organization (NEDO). K.Y. acknowledges the support by MEXT as “Program for Promoting Research on the Supercomputer Fugaku” (Realization of innovative light energy conversion materials utilizing the supercomputer Fugaku, Grant No. JPMXP1020210317). The DFT computations were performed at the Research Center for Computational Science, Okazaki, Japan (Project No. 22-IMS-C039).

The authors have no conflicts to disclose.

Masanori Kaneko: Data curation (equal); Investigation (equal); Software (lead); Writing – review & editing (equal). Vikas Nandal: Investigation (equal); Writing – review & editing (equal). Koichi Yamashita: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Writing – review & editing (equal). Kazuhiko Seki: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Writing – original draft (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.