Accurate Prediction of Band Structure of FeS₂: A Hard Quest of Advanced First-Principles Approaches

1 Introduction

Iron disulfide FeS₂ was studied extensively in the last century in the desire of understanding the structural and electronic properties of transition metal dichalcogenides (TMDC) featuring localized or band-like d electrons (Hulliger and Mooser, 1965a; Hulliger and Mooser, 1965b; Goodenough, 1972; Wilson, 1972; Li et al., 1974; Schlegel and Wachter, 1976; Folkerts et al., 1987). Since 1980s (Chatzitheodorou et al., 1986; Ennaoui et al., 1986), increasing practical interest has been drawn to pyrite FeS₂ for its potential as a cheap and competitive candidate material for efficient solar energy conversion (Wadia et al., 2009) because of its natural abundance, non-toxicity, suitable optical gap and extraordinarily large absorption coefficient (Ferrer et al., 1990; Ennaoui et al., 1993). This has led to new solutions under various optoelectronic scenarios, including photovoltaics (Khalid et al., 2018), photo-catalysis (Tian et al., 2015; Barawi et al., 2016), solid-state photo-capacitors (Gong et al., 2013a) and photo-detectors (Wang et al., 2012; Gong et al., 2013b). However, practical application of FeS₂-based optoelectronic devices is seriously hampered by its unexpected low efficiency due to a loss of open-circuit voltage V_OC (Cabán-Acevedo et al., 2014). A number of factors possibly responsible for the low V_OC have been suggested and examined, including the intrinsic and defect surface states (Bronold et al., 1994; Sun et al., 2011; Herbert et al., 2013; Lazić et al., 2013; Cabán-Acevedo et al., 2014; Limpinsel et al., 2014; Walter et al., 2017), bulk sulfur deficiency (Birkholz et al., 1991; Cabán-Acevedo et al., 2014; Shukla et al., 2016) and presence of the metastable marcasite phase as a small-gap impurity (Spagnoli et al., 2010; Sun et al., 2011; Schena et al., 2013).

Despite progress towards understanding the origin of the low V_OC in pyrite FeS₂ (Rahman et al., 2020), consensus is still not reached on the fundamental band gaps of the two FeS₂ phases. Experimentally, values varying from 0.6 to 2.6 eV have been reported for pyrite, primarily due to differences in sample preparation, measuring technique, and analytical model of spectra used in experimental studies (Ferrer et al., 1990; Ennaoui et al., 1993). Measurements of the pyrite band gap are generally carried out through optical absorption spectroscopy (Schlegel and Wachter, 1976; Kou and Seehra, 1978; Ennaoui et al., 1993), which features the neutral excitation (exciton) instead of the charged one as in the photo-electron spectroscopy (PES). Therefore, the measured excitation energies are in fact coupled to the electron-hole binding. Careful investigation by absorption spectroscopy for the marcasite phase is done only recently and gives an optical gap similar to pyrite, which essentially precludes the possibility of marcasite being the culprit for the low V_OC of FeS₂ photovoltaics (Sánchez et al., 2016; Wu et al., 2016). Furthermore, even though PES measurements of pyrite FeS₂ have been conducted (Ohsawa et al., 1974; van der Heide et al., 1980; Folkerts et al., 1987; Mamiya et al., 1997; Ollonqvist et al., 1997; Nesbitt et al., 2003), combined studies of direct and inverse PES (IPS) for regions near the Fermi level are rare. Reported relevant works (Folkerts et al., 1987; Mamiya et al., 1997) were done more than 20 years ago and the spectra were not resolved enough to identify a well-defined fundamental band gap.

Difficulties in characterizing band structures of FeS₂ polymorphs are also encountered from the perspective of first-principles calculations. Within the framework of density functional theory (DFT) (Hohenberg and Kohn, 1964), calculations with Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) (Perdew et al., 1996a) predict pyrite to have a band gap of about 0.3 eV smaller than the experimental value of 0.95 eV as generally accepted (Ennaoui et al., 1993; Schena et al., 2013; Kolb and Kolpak, 2013; Li et al., 2015). Considering the well-known band gap problem of local density approximation (LDA) or GGA (Perdew et al., 1982), orbital-dependent functionals in spirit of generalized Kohn-Sham (GKS) DFT (Seidl et al., 1996; Perdew et al., 2017; Zhang et al., 2020) are also employed to tackle the problem, e.g. PBE plus the Hubbard-U correction (DFT + U) and hybrid functionals (Becke, 1993a; Becke, 1993b; Perdew et al., 1996b). Using an ad hocU of 2 eV, the PBE + U method is able to reproduce the experimental band gap (Sun et al., 2011; Hu et al., 2012; Li et al., 2018) but meanwhile deteriorates the simulated optical spectra compared to PBE (Choi et al., 2012; Schena et al., 2013). Furthermore, despite the good performance in predicting band gaps for typical semiconducting materials (Heyd et al., 2005; Paier et al., 2006b,a; Marsman et al., 2008), hybrid functionals such as Heyd-Scuseria-Ernzerhof (HSE) method (Heyd et al., 2003, 2006) have been shown to give large band gaps for pyrite of over 2 eV (Muscat et al., 2002; Sun et al., 2011; Choi et al., 2012; Hu et al., 2012; Schena et al., 2013; Liu et al., 2019). There are also works using beyond-DFT methods, particularly, the GW method based on many-body perturbation theory (MBPT) (Hedin, 1965). However, the GW results for the pyrite phase are rather scattered, ranging from 0.3 to 1.1 eV (Choi et al., 2012; Lehner et al., 2012; Kolb and Kolpak, 2013; Schena et al., 2013). It is worth noting that Schena and coworkers conducted the state-of-the-art all-electron G₀W₀ calculations with the linearized augmented plane-wave (LAPW) basis for both pyrite and marcasite, and report a pyrite band gap only about 0.3 eV (Schena et al., 2013). The GW gap value is smaller than that from PBE, which is rarely observed in GW practices and hence deserves closer investigation.

For GW implementations involving explicit summation of states, it is established recently by a number of works (Friedrich et al., 2006; Friedrich et al., 2011a; Friedrich et al., 2011b; Klimes et al., 2014; Jiang and Blaha, 2016; Nabok et al., 2016; Jiang, 2018; Zhang and Jiang, 2019; Ren et al., 2021) that an accurate description of high-lying empty states is essential to give accurate correlation self-energy operator and consequent QP band structure. In the pseudo-potential framework, one can improve the accuracy by using a norm-conserving potential with specifically tailored projectors at high energies (Klimes et al., 2014; van Setten et al., 2018). In all-electron calculations with the LAPW basis set, local orbitals with large energy parameters (usually 10^1∼2 Ry higher than the Fermi level) are introduced as additional basis functions to remove the linearization error in unoccupied states up high in the conduction band regime (Friedrich et al., 2006, 2011a,b; Jiang and Blaha, 2016; Nabok et al., 2016). The LAPW basis extended by these high-energy local orbitals (HLOs), termed as LAPW + HLOs, has succeeded in helping produce accurate QP band structures in good agreement with experiment for a variety of semiconductors (Jiang and Blaha, 2016) including the conventionally challenging systems such as ZnO (Friedrich et al., 2011a; Friedrich et al., 2011b; Stankovski et al., 2011; Jiang and Blaha, 2016; Nabok et al., 2016), d/f-electron oxides (Jiang, 2018) and cuprous and silver halides (Zhang and Jiang, 2019). Particularly, the effects of including HLOs on the QP correction have been demonstrated quantitatively to be larger for states with stronger metal-d characters (Zhang and Jiang, 2019). For the FeS₂ polymorphs with states of significant Fe-3d characters in both valence and low-energy conduction band regimes, GW with LAPW + HLOs is likely to give better description of the QP energies and dispersion relation than that with the standard LAPW basis.

A competitive alternative in the DFT framework to GW for band structure prediction is the doubly screened hybrid (DSH) functional method (Cui et al., 2018) in the category of hybrid functionals with system-dependent parameters (Zhang et al., 2020). Derived from a model dielectric function (Cappellini et al., 1993; Shimazaki and Yoshihiro, 2008), the exchange-correlation potential in DSH can be regarded as a further approximation to the Coulomb hole and screened exchange (COHSEX) approximation to the GW self-energy, and is able to capture both dielectric and metallic screening in the exchange interaction (Cui et al., 2018). It is shown that the DSH can evaluate band gaps of typical sp semiconductors with accuracy comparable to GW with the LAPW + HLOs basis while only at modest computational cost (Cui et al., 2018). Furthermore, the one-shot variant DSH0 can outperform fixed-parameter hybrid functionals for band gap predictions in a wide range of materials including narrow-gap semiconductors and transition metal mono-oxides (Cui et al., 2018; Liu et al., 2020). Hence we consider DSH as a hopeful approach to solve the FeS₂ band gap puzzle within the GKS framework of DFT.

In the present work, we investigate the electronic band structures of the pyrite and marcasite polymorphs of FeS₂ by applying the state-of-the-art all-electron GW method with the LAPW + HLOs basis. For comparison, we examine the results from GW with the standard LAPW basis as well. We also investigate the performances of several hybrid functionals, including PBE0 (Perdew et al., 1996b), HSE06 (Heyd et al., 2003, Heyd et al., 2006), screened-exchange-PBE hybrid functional (SX-PBE) (Bylander and Kleinman, 1990; Seidl et al., 1996) and DSH (Cui et al., 2018), in attempt to obtain insights into the failure of the conventional fixed-parameter functionals in predicting the band gap of FeS₂.

2 Theory and Methods

2.1 The GW Method

The central task of the GW method is to solve the quasi-particle (QP) equation with the self-energy operator Σ in the frequency domain expressed as (Hedin, 1965)

where G is the time-ordered Green’s function

With ψ_nk and ε_nk being the wave function and energy of the single-particle state $\left|n\mathbf{k}\rangle \right.⟩$ respectively, μ the chemical potential, and δ and η positive infinitesimals. Atomic units are used throughout the paper. The screened Coulomb interaction W writes

where v (r, r′) = 1/|r − r′| is the bare Coulomb interaction and ε(r, r′; ω) is the microscopic dielectric function calculated at the level of random phase approximation (RPA). In principle, Eqs 1–3 have to be solved self-consistently along with the Dyson equation for the Green’s function (Hedin, 1965). However, due to the computational cost and generally unsatisfactory results of the fully self-consistent GW for solids [e.g. Grumet et al. (2018)], one usually turns to the non-self-consistent variant G₀W₀. Considering the resemblance of KS and QP wave functions in weakly correlated systems (Hybertsen and Louie, 1986), the self-energy or QP energy ${\epsilon }_{n\mathbf{k}}^{\mathrm{Q}\mathrm{P}}$can be computed perturbatively upon the acquisition of Σ from the KS states as

where V_xc is the KS exchange-correlation potential and Z_nk a renormalization factor. One can further perform the so-called energy-only self-consistent GW₀ calculations, where QP energies ${\epsilon }_{n\mathbf{k}}^{\mathrm{Q}\mathrm{P}}$ in place of ${\epsilon }_{n\mathbf{k}}^{\mathrm{K}\mathrm{S}}$ in Eq. 2 are updated iteratively while W is kept the same as in G₀W₀ (Shishkin and Kresse, 2007). The GW method has been implemented in various numerical frameworks (Jiang, 2011; Golze et al., 2019). For a detailed explanation of the basic theory and computational techniques used in the present GW implementation, the readers can refer to Jiang et al. (2013).

2.2 All-Electron Calculations With HLOs-Extended LAPW Basis

In the all-electron framework with LAPW, KS wave functions are expanded by the LAPW basis (Andersen, 1975; Singh and Nordström, 2006; Blaha et al., 2020)

where V_α is the region enclosed by the muffin-tin (MT) sphere of atom α centered at r_α with radius $R_{\mathrm{M}\mathrm{T}}^{\alpha }$, r^α = r − r_α, u_αl (E_αl) is the solution of radial KS equation inside V_α at chosen energy E_αl, ${\dot{u}}_{\alpha l}(E_{\alpha l})\equiv {\left.\partial u_{\alpha l}(E)/\partial E\right|}_{E=E_{\alpha l}}$, and $Y_l^m$ is the spherical harmonic function. The coefficients $A_{\alpha lm}^{\mathbf{k}+\mathbf{G}}$ and $B_{\alpha lm}^{\mathbf{k}+\mathbf{G}}$ are determined by enforcing that ${\phi }_{\mathbf{k}+\mathbf{G}}^{\text{LAPW}}(\mathbf{r})$ be smooth at the boundary of V_α. Local orbitals (LOs) which vanish outside the atomic spheres are proposed to supplement the LAPW basis to better describe the semi-core states (Singh, 1991). Inside the atomic sphere V_α, LOs take the following form

where $E_{\alpha l}^{\mathrm{L}\mathrm{O},i}$ is the energy parameter for the ith LO centered on atom α with angular and azimuthal quantum numbers l and m, respectively.

HLOs fall into the category of LOs with $E_{\alpha l}^{\mathrm{L}\mathrm{O}}$ typically 10 ∼ 100 Ry above the Fermi level. Such extra LOs have been found to facilitate accurate description of unoccupied states by remedying the linearization error therein when using the LAPW basis (Krasovskii et al., 1994; Krasovskii, 1997; Friedrich et al., 2006; Michalicek et al., 2013). In ground state calculations with LDA/GGA or hybrid functionals, the error causes no essential difficulties, since only occupied and low-lying unoccupied states are involved which are usually handled in sufficient accuracy with the usual or standard LAPW basis generated as default in popular DFT implementations with LAPW basis (Blaha et al., 2020). However, the error can be detrimental to the numerical accuracy of methods where the summation over unoccupied states is required, e.g. GW and DFT methods with density approximations belonging to the fifth rung of Jacobi ladder (Perdew and Schmidt, 2001) such as the adiabatic-connection dissipation-fluctuation (ACFD) calculation under RPA for ground-state energy (Ren et al., 2012; Cui et al., 2016; Zhang et al., 2018). In these methods, the completeness of summation and quality of unoccupied states play a crucial role. Previous GW studies (Jiang and Blaha, 2016; Jiang, 2018; Zhang and Jiang, 2019; Shen et al., 2020) have suggested that both can be taken into account by including localized orbitals energetically higher than the Fermi level in addition to the standard LAPW basis. HLOs have been shown to effectively improve the optical properties (Krasovskii et al., 1994; Krasovskii, 1997), NMR chemical shifts (Laskowski and Blaha, 2012; Laskowski and Blaha, 2014), GW QP energies (Friedrich et al., 2006; Friedrich et al., 2011a; Jiang and Blaha, 2016; Nabok et al., 2016; Jiang, 2018; Zhang and Jiang, 2019), optimized effective potential (Betzinger et al., 2011, 2012) and RPA correlation energy (Betzinger et al., 2015).

In the current implementation, HLOs are generated systematically by following the way described by Laskowski and Blaha (2012). The quality of LAPW + HLOs is controlled by two parameters besides those for the LAPW basis, namely, the additional number of nodes in the radial function of highest energy local orbital with respect to that of the LAPW function with the same angular quantum number and the maximum angular quantum number of used HLOs, denoted as n_LO and $l_{\mathrm{m}\mathrm{a}\mathrm{x}}^{(\mathrm{L}\mathrm{O})}$, respectively. Generally speaking, the larger n_LO and $l_{\mathrm{m}\mathrm{a}\mathrm{x}}^{(\mathrm{L}\mathrm{O})}$ are, the higher the HLOs can reach in the energy space. We use n_LO = 0 to denote the usual or standard LAPW basis. Since the convergence rate of the QP energy with respect to the two parameters can be different for states featuring distinct atomic characters, careful convergence check is required to obtain numerically accurate GW results.

2.3 Hybrid Functionals

Hybrid functionals have been widely used in first-principles simulations of condensed matter for their good balance between performance and computational cost, and have been actively developed to further exploit the potential of its particular functional form. Readers interested in detailed description on the current status of hybrid functional development are directed to several recent reviews (Kümmel and Kronik, 2008; Baer et al., 2010; Maier et al., 2019; Zhang et al., 2020). Here we briefly introduce the general formalism of the range-separated hybrid functionals and the variants relevant to the current study.

The essential ingredient in hybrid functional methods is the exchange-correlation energy E_xc or potential V_xc composed of non-local orbital-dependent (screened) Hartree-Fock (HF) exchange terms. In the present work, we focus on hybrid functionals with V_xc in the range-separated form as (Zhang et al., 2020)

where $V_{\text{x}}^{\mathrm{H}\mathrm{F},\mathrm{s}\mathrm{r}}$ and $V_{\text{x}}^{\mathrm{H}\mathrm{F},\mathrm{l}\mathrm{r}}$ are the short- and long-ranged Fock exchange potentials, respectively, which, using the reduced density-matrix defined as $\rho (\mathbf{x},{\mathbf{x}}^{\prime })\equiv {\sum }_{i\in \text{occ}}{\psi }_i(\mathbf{x}){\psi }_i*({\mathbf{x}}^{\prime })$, can be written as

In Eq. 8v^sr (r, r′; μ) denotes the short-ranged Coulomb interaction of a certain form characterized by screening parameter μ (Zhang et al., 2020). x denotes collectively the spatial and spin coordinates of an electron, x ≡ (r, σ). $V_{\text{x}}^{\mathrm{S}\mathrm{L},\mathrm{s}\mathrm{r}}$ and $V_{\text{x}}^{\mathrm{S}\mathrm{L},\mathrm{l}\mathrm{r}}$ are the semi-local (SL) counterparts of the exchange potentials in LDA, GGA or meta-GGA. μ and the mixing ratios α_sr and α_lr are the adjustable parameters of the hybrid functional form.

Conventionally, the parameters are determined by either theoretical analysis or fitting against some dataset of particular properties, and then applied to other systems as fixed. Famous examples of the fixed-parameter hybrid functionals include PBE0 α_sr = α_lr = 1/4 (Perdew et al., 1996b) and the HSE series α_sr = 1/4, α_lr = 0, μ = 0.2–0.3 Å^-1 (Heyd et al., 2003; Heyd et al., 2006). Recently, hybrid functionals with system-dependent parameters are developed by several groups (Shimazaki and Yoshihiro, 2008; Marques et al., 2011; Kronik et al., 2012; Koller et al., 2013; Skone et al., 2014; Chen et al., 2018; Cui et al., 2018). Among different methods, the doubly screened hybrid (DSH) functional has been demonstrated as a competitive candidate for accurate description of band structures of both wide- and narrow-gap semiconductors (Cui et al., 2018). The underlying idea of DSH is to approximate the screening effect in solids by employing the Bechstedt model dielectric function (Bechstedt et al., 1992)

where ɛ_M is the macroscopic dielectric constant, q_TF the Thomas-Fermi wave vector and α an empirical parameter chosen for semiconductors (Cappellini et al., 1993). A screened Coulomb interaction can be derived from this model to take both dielectric and metallic screening into account, leading to parameters in Eq. 7 as

The corresponding short-ranged Coulomb interaction in Eq. 8 is

where erfc is the complementary error function. In practice, an initial ɛ_M is required, which can be obtained from the PBE calculation or experimental measurements, to construct the DSH potential and solve the GKS equation. The resulting single-particle states act as the inputs to compute a new ɛ_M, which is in turn used to update the DSH potential. The self-consistent loop stops when ɛ_M is converged. Alternatively, one can break after solving the GKS equation with the initial ɛ_M, leading to the one-shot scheme denoted as DSH0.

2.4 Computational Details

The unit cells of pyrite and marcasite FeS₂ used in our calculations are shown in the left panel of Figure 1. The crystal structure of pyrite FeS₂ (Figure 1A) can be viewed as a faced-centered cubic cell of Fe atoms with S₂ dumbbells occupying the octahedral interstitials and pointing to different <111> crystallographic axes. The anion coordination octahedra (FeS6) are connected only through sharing vertices. In the orthorhombic marcasite phase (Figure 1C), (FeS6) are connected by sharing edges with the two neighbors along c-axis and linked together through sharing vertices on the aOb plane. In terms of lattice parameters, we use a = 5.418 Å, u = 0.3850 for pyrite (space group $Pa\bar{3}$) and a = 4.443 Å, b = 5.425 Å, c = 3.387 Å, u = 0.2005, v = 0.3783 for marcasite (space group Pnnm). These values follow the results from X-ray diffraction experiments at ambient conditions (Brostigen and Kjekshus, 1969; Brostigen et al., 1973; Chattopadhyay and Von Schnering, 1985; Zuñiga-Puelles et al., 2019). The corresponding S-S bond lengths in the two polymorphs are 2.16 and 2.21 Å, respectively.

FIGURE 1. Lattice structures (A,C) and first Brillouin zones (B,D) of pyrite (upper) and marcasite (lower) phases of FeS₂. Brown and yellow spheres represent Fe and S atoms, respectively.

The present all-electron GW calculations are performed by the GW facilities in the GAP2 program (Jiang et al., 2013; Jiang and Blaha, 2016) interfaced to WIEN2k (Blaha et al., 2001, 2020). Results in both G₀W₀ and GW₀ schemes are presented, where KS orbital energies and wave functions calculated with the PBE (Perdew et al., 1996a) GGA are used as the input to construct one-body Green’s function and screened Coulomb interaction. The KS states are obtained by using charge density pre-converged under self-consistent field (SCF) calculation with PBE and the standard LAPW basis. The energy criterion for convergence of SCF iterations is set to 10^–8 Rydberg (Ry). 64 (4 × 4 × 4) and 120 (5 × 4 × 6) k points are sampled in the first Brillouin zones of pyrite and marcasite FeS₂, respectively. All available unoccupied states are considered in the summation of states for screened Coulomb interaction W and self-energy Σ. Mixed product basis is used to describe the wave function products in the two-point functions, e.g. W and Σ (Aryasetiawan and Gunnarsson, 1994; Kotani and van Schilfgaarde, 2002). We choose Q = 0.75 and $l_{\text{max}}^{\text{MB}}=3$for the interstitial plane wave and MT product basis, respectively [Jiang et al. (2013) for the meanings of these parameters]. LAPW and LOs with E_αl < 20 Ry are used to build the MT product basis. Frequency dependence of W is treated explicitly on a 16-point double Gauss-Legendre grids along the positive imaginary axis. Σ on the same grid is calculated and analytically continued to the real axis (Rojas et al., 1995). A rather coarse k/q-point mesh, 2 × 2 × 2 for pyrite and 4 × 2 × 4 for marcasite, is sufficient to converge the direct band gap at the Γ point ${E_{\mathrm{g}}}^{\mathrm{\Gamma }}$ within 0.01 eV. The QP band structure diagrams along particular k-point paths (see the right panel of Figure 1) are calculated by interpolating the QP energies obtained with the above mesh using the Fourier interpolation technique (Pickett et al., 1988).

In terms of the LAPW basis, the usual or standard LAPW basis set is created automatically in the recent version of WIEN2k (Blaha et al., 2001), which is actually a mixture of the APW + lo basis for the valence states (Madsen et al., 2001), the ordinary LAPW basis for higher l channels up to l_max = 10 and additional local orbitals (LOs) for semi-core Fe-3s and Fe-3p states (Blaha et al., 2001). The convergence with respect to the two HLOs parameters n_LO and $l_{\mathrm{m}\mathrm{a}\mathrm{x}}^{(\mathrm{L}\mathrm{O})}$ is investigated, the latter being represented by $\mathrm{\Delta }{l}_{\mathrm{L}\mathrm{O}}=l_{\mathrm{m}\mathrm{a}\mathrm{x}}^{(\mathrm{L}\mathrm{O})}-l_{\mathrm{m}\mathrm{a}\mathrm{x}}^{(\mathrm{v})}$ where $l_{\mathrm{m}\mathrm{a}\mathrm{x}}^{(\mathrm{v})}$ is the largest l of valence orbitals for each element. In the present study, $l_{\mathrm{m}\mathrm{a}\mathrm{x}}^{(\mathrm{v})}=2$ and $l_{\mathrm{m}\mathrm{a}\mathrm{x}}^{(\mathrm{v})}=1$ for Fe and S, respectively. Since the convergence with respect to HLOs parameters are decoupled from the choice of k-point mesh, we choose marcasite with a coarse 2 × 1 × 2 mesh for HLOs convergence test. RKmax ≡ R_MT,minK_max = 7.0 is chosen for the plane-wave cut-off in the interstitial region, where R_MT,min is the minimal muffin-tin radius R_MT used in the lattice. In the present FeS₂ case, R_MT is set to 2.1 Bohr for Fe and 1.9 Bohr for S. Using RKmax = 9.0 will reduce the band gap from GW (LAPW + HLOs) by less than 0.03 eV, indicating that adequate accuracy can be delivered with the current RKmax = 7.0 setup. Due to limited computational resources, RKmax = 6.0 is used for HLOs convergence test. Following Laskowski and Blaha (2012), the linear independence of HLO basis functions is assured by choosing the energy parameters such that the overlap between the HLO radial functions is smaller than a threshold, which is 0.6 in the present work.

For hybrid functional calculations, we consider PBE0 (Perdew et al., 1996b), HSE06 (Heyd et al., 2003, 2006) and screened exchange SX-PBE (Bylander and Kleinman, 1990) methods as well as DSH. All hybrid functional calculations are performed with the projector augmented waves (PAW) method (Blöchl, 1994) implemented in the Vienna ab-initio Simulation Package (VASP) (Kresse and Furthmüller, 1996). The static dielectric function is calculated from the average of diagonal elements of macroscopic dielectric tensor computed by using density functional perturbation theory (DFPT) with local field effect included (Baroni et al., 2001). Apart from 3d and 4s, the 3s and 3p electrons of Fe are also treated explicitly in the valence region. The Thomas-Fermi wave vectors are 2.57 and 2.56 Å^-1 for pyrite and marcasite FeS₂, respectively. The cut-off energy of plane-wave basis for wave function expansion is chosen as 400 eV, which is sufficient to converge ${E_{\mathrm{g}}}^{\mathrm{\Gamma }}$ of both FeS₂ polymorphs within 2 meV. In terms of k-point mesh, 64 (4 × 4 × 4) and 120 (5 × 4 × 6) k points are sampled in the first Brillouin zones of pyrite and marcasite for the self-consistent calculations, respectively. Using a finer 6 × 6 × 6 sampling for pyrite will change the band gap by less than 0.01 eV, and hence we consider the results well converged with respect to the k-point mesh. The energy convergence criterion is chosen to be 10^–6 eV for the SCF iterations.

3 Results and Discussion

3.1The GW Results

In this part, we present the electronic band structures of pyrite and marcasite FeS₂ computed by the all-electron GW method. In particular, we analyse the effect of high-energy local-orbitals (HLOs) by comparing the results from GW with the standard LAPW and LAPW + HLOs basis.

3.1.1 Convergence of QP Energies with Respect to HLOs Parameters

To achieve a balance between the computational cost and numerical accuracy of the LAPW + HLOs based GW method, we have to decide an optimized HLOs setup for the FeS₂ polymorphs of interest. That is to say, certain convergence with respect to the two HLOs parameters, namely n_LO and Δl_LO, must be achieved for the QP band structures of both polymorphs, while the number of basis functions should be kept as few as possible. To simplify the notation, we denote the setup of HLOs by (n_LO, Δl_LO) so that (1, 1) indicates a set of HLOs with n_LO = 1 and Δl_LO = 1, for example. Since we are most interested in the band gaps (direct and indirect) of the systems, we choose the indirect band gap from the X point to the Γ point, $E_{\mathrm{g}}^{\text{X–}\mathrm{\Gamma }}$, as the descriptor for the band structure, and investigate its dependence on the two HLOs parameters for marcasite.

Before discussing the results, we briefly illustrate the appropriateness of this choice. First of all, $E_{\mathrm{g}}^{\text{X–}\mathrm{\Gamma }}$ is a representative band gap energy for pyrite and marcasite FeS₂. This is because in both phases, the topmost valence state at the X point, X_v, is close to the valence band maximum (VBM) and the bottommost conduction state at the Γ point, Γ_c, is the conduction band minimum (CBM) (that is the case for marcasite given the coarse 2 × 1 × 2 k mesh in the convergence study). Second, either X_v or Γ_c has similar atomic contributions in the two polymorphs, and the effects of HLOs on such states are also similar, as shown in the results for other polymorphs like zinc-blende and wurtzite ZnO (Jiang and Blaha, 2016). Therefore the parameters optimized for marcasite are considered transferable and can be applied to the pyrite polymorph. Last but not least, as we will discuss later, the effects of HLOs on X_v and Γ_c differ significantly, avoiding considerable error cancellation in change of the QP correction to the band gap upon including HLOs.

Figure 2 summarizes the results of convergence study for the G₀W₀@PBE method. $E_{\mathrm{g}}^{\text{X–}\mathrm{\Gamma }}$ (Figure 2A) is about 1.0 eV with the standard LAPW basis (n_LO = 0) and is significantly increased by extending LAPW with HLOs. One can see that the convergence rate of $E_{\mathrm{g}}^{\text{X–}\mathrm{\Gamma }}$ with respect to n_LO differs with different Δl_LO, and is faster for lower Δl_LO. The reverse is also true, i.e. the convergence with respect to Δl_LO is faster when n_LO is smaller. It clearly indicates that the convergence with respect to n_LO and Δl_LO is coupled. Increasing HLOs parameters from (4, 4) to (5, 5) changes $E_{\mathrm{g}}^{\text{X–}\mathrm{\Gamma }}$ by less than 0.05 eV, indicating that HLOs (4,4) is able to deliver an adequate accuracy. Therefore, unless stated otherwise, HLOs (4, 4), amounting to 196 and 145 HLOs for Fe and S atoms, respectively, is considered optimized and will be used in the subsequent GW calculations denoted by LAPW + HLOs. The energy parameters for HLOs (4, 4) can be found in Table 1.

FIGURE 2. Dependence on HLOs parameters of G₀W₀@PBE calculated (A) the indirect band gap between X and Γ, $E_{\mathrm{g}}^{\text{X–}\mathrm{\Gamma }}$, and self-energy corrections to (B) bottommost conduction state at the Γ point, $\mathrm{\Delta }{ϵ}_{\text{c}}^{\mathrm{\Gamma }}$, and (C) topmost valence state at the X point, $\mathrm{\Delta }{ϵ}_{\text{v}}^X$, for marcasite FeS₂.

TABLE 1. Energy parameters (unit: Rydberg) of high-energy local orbitals (HLOs) used in pyrite FeS₂ corresponding to HLOs parameters n_LO = 4, Δl_LO = 4 for Fe and S. Those for marcasite are essentially the same with difference by 0.02 Ry at most in each element and l channel.

It is worth noting that the effect of including HLOs on the QP correction to $E_{\mathrm{g}}^{\text{X–}\mathrm{\Gamma }}$ is different from those on the valence and conduction states. To illustrate this, we show the dependence on n_LO and Δl_LO of the self-energy corrections to Γ_c ($\mathrm{\Delta }{\epsilon }_{\text{c}}^{\mathrm{\Gamma }}$) and X_v ($\mathrm{\Delta }{\epsilon }_{\text{v}}^X$) states in Figures 2B,C, respectively. Both $\mathrm{\Delta }{\epsilon }_{\text{c}}^{\mathrm{\Gamma }}$ and $\mathrm{\Delta }{\epsilon }_{\text{v}}^X$ decrease with increasing n_LO or Δl_LO, but the former converges much faster than the latter, which agrees with the general trend observed previously (Jiang and Blaha, 2016; Zhang and Jiang, 2019). With the standard LAPW basis, G₀W₀@PBE gives $\mathrm{\Delta }{\epsilon }_{\text{c}}^{\mathrm{\Gamma }}=0.48$ eV and $\mathrm{\Delta }{\epsilon }_{\text{v}}^X=0.76$ eV, indicating a negative QP correction to the band gap, which is rarely observed in LDA/GGA-based GW calculations of semiconductors (Jiang and Blaha, 2016). When HLOs (5, 5) are included, $\mathrm{\Delta }{\epsilon }_{\text{c}}^{\mathrm{\Gamma }}$decreases by 0.3 eV, much smaller compared to the decreasing of 1.2 eV in $\mathrm{\Delta }{\epsilon }_{\text{v}}^X$. Such biased effects of including HLOs on valence and conduction band states can be attributed to the difference in atomic characteristics between the states, and will be further discussed in the following sections.

3.1.2 Quasi-Particle Band Gaps

After having obtained the optimized HLOs, we perform the PBE-based GW calculations for pyrite and marcasite FeS₂ with the LAPW + HLOs basis set, and compare with the PBE method and GW with the standard LAPW basis.

The band gaps of pyrite and marcasite FeS₂ calculated by PBE and GW methods are presented in Table 2. The fundamental band gaps are obtained by computing the band energies along the k-point paths indicated in the right panel of Figure 1. In the PBE reference, pyrite and marcasite are predicted to have indirect fundamental band gaps of 0.70 and 0.83 eV, respectively. Our PBE results are consistent with those from previous all-electron LAPW study (Schena et al., 2013) and close to the recently reported optical band gaps obtained from diffuse reflectance spectroscopy (Sánchez et al., 2016). However, our PBE band gap for pyrite is slightly larger than several reported PBE results (Sun et al., 2011; Kolb and Kolpak, 2013; Lazić et al., 2013; Li et al., 2015; Zhang et al., 2018). This can be attributed to the use of different lattice structures in those studies (Eyert et al., 1998; Lazić et al., 2013; Schena et al., 2013) from the current work. Particularly, geometry optimization by PBE (Eyert et al., 1998; Schena et al., 2013) generally gives a longer S-S dimer, which leads to smaller splitting between bonding and anti-bonding S-3pσ orbitals and a consequent shrink in the band gap.

TABLE 2. Fundamental band gap (indicated by “fund.”) and other direct and indirect band gaps (unit: eV) for pyrite and marcasite FeS₂ calculated by PBE and GW methods. Results from previous GW studies and experimental measurements are presented for comparison. To simplify the notation, we use “L” and “L + H” to denote the standard LAPW and LAPW + HLOs basis sets, respectively. PBE is used as the starting point for G₀W₀ and GW₀ calculations unless stated otherwise.

For GW calculations with the standard LAPW basis, the QP fundamental band gaps by G₀W₀@PBE are smaller than the PBE counterparts in both FeS₂ polymorphs. Pyrite FeS₂ is predicted to have a band gap of only 0.06 eV, which is 0.64 eV smaller than that by PBE. The negative QP correction for pyrite band gap has been reported by Schena et al. (2013). The QP fundamental band gap for marcasite predicted by G₀W₀@PBE (LAPW) is also smaller than PBE, while the change (0.26 eV) is less dramatic than that for pyrite. Such negative QP corrections to LDA/GGA band gaps are uncommon in GW studies for closed-shell systems (Klimes et al., 2014; Jiang and Blaha, 2016; van Setten et al., 2017; Zhang and Jiang, 2019) as well as open-shell d/f-electron semiconductors (Jiang, 2018). Switching on self-consistency of the Green’s function by GW₀@PBE further reduces the fundamental band gaps of FeS₂. In particular, pyrite is predicted to be metallic by GW₀@PBE, which disagrees qualitatively with its semiconducting nature in experiment (Ennaoui et al., 1993). For other direct and indirect band gaps, those for Γ → Γ and X → Γ in pyrite and marcasite and M → Γ in marcasite are decreased from PBE to G₀W₀@PBE (LAPW). The decrease is largest for the Γ → Γ gap in the two phases, 0.71 and 0.86 eV for pyrite and marcasite, respectively. On the other hand, the gaps for X → X in pyrite, Γ → T and X → T in marcasite are increased by 0.28, 0.25, and 0.58 eV, respectively. However, it should be noted that the distinction in signs of corrections to the QP gaps in different channels should not be considered as intrinsic for FeS₂. Instead, it is an artifact as a result of the incomplete basis, which we will discuss in details below.

Now we turn to the LAPW + HLOs-based GW calculations. With the G₀W₀@PBE method, including HLOs increases the QP fundamental gap by 0.98 eV for pyrite and 0.58 eV for marcasite. The resulting G₀W₀@PBE band gaps are 1.04 and 1.15 eV for pyrite and marcasite, respectively. In contrast, all band gaps investigated are increased by G₀W₀ with LAPW + HLOs compared to their PBE counterparts. We note that HLOs have distinct effects among band gaps for different channels. Once the HLOs are included, band gaps for channels with the conduction state at the Γ point are increased by about 1 eV. On the other hand, the QP correction to the X → X band gap in pyrite increases by only 0.18 eV. Moreover, the gaps for Γ → T and X → T in marcasite even decrease. With the LAPW + HLOs basis, using GW₀ to switch on partial self-consistency further increases the band gaps, but the change is moderate and no more than 0.1 eV.

As explained at the beginning, the fundamental band gap of FeS₂ has been controversial in the recent decades, partly due to the widely varying experimental values (Ennaoui et al., 1993). In the present study, the GW₀@PBE method with the LAPW + HLOs predicts that pyrite and marcasite have indirect fundamental band gaps of 1.14 and 1.16 eV, respectively. The GW₀ gap of pyrite is slightly larger than the generally accepted experimental value of 0.95 eV (Ennaoui et al., 1993). Furthermore, the fact that the two polymorphs have almost identical band gaps is consistent with the optical measurements by Sánchez et al. (2016), although our predicted band gaps are about 0.3 eV larger. However, it should be noted that one must take exciton binding energy E_B into account for a meaningful comparison between the QP fundamental band gap and experimentally measured optical gap. The difference between the fundamental and optical gaps can be significant when the exciton is localized, i.e. of Frenkel type (Fox, 2010). On the other hand, while it is more straightforward to compare the QP gap with spectral data from direct and inverse PES (Folkerts et al., 1987; Mamiya et al., 1997), the resolutions of available measurements for pyrite FeS₂ are too low to extract a meaningful gap value for comparison. Moreover, to the best knowledge of the authors, no data of combined PES/IPS measurements are available for marcasite. Therefore, further experimental studies are required to determine and verify the band gaps of the FeS₂ polymorphs.

To close this part, we highlight that the present work resolve two issues reported in previous GW studies in terms of QP band structures of FeS₂. First, Schena et al. (2013) performed a G₀W₀@PBE study on pyrite and marcasite FeS₂ with similar HLOs-extended LAPW basis. The fundamental band gap of pyrite was estimated as about 0.3 eV, by which the authors claimed to explain the low V_OC encountered in the pyrite solar cell. However, according to our convergence study, such a small band gap is likely to result from inadequate convergence with respect to HLOs. More specifically, the largest angular momentum of HLOs $l_{\mathrm{m}\mathrm{a}\mathrm{x}}^{(\mathrm{L}\mathrm{O})}$ used in Schena et al. (2013) is 3, i.e. f orbital, while $l_{\mathrm{m}\mathrm{a}\mathrm{x}}^{(\mathrm{L}\mathrm{O})}=2+4=6$ (i orbital) is used in the optimized HLOs of the present work. As a result, the highest energy covered by HLOs in Schena et al. (2013) (800 eV) is much smaller than that used in the present work (about 1800 eV). Second, fully self-consistent GW (scGW) and quasi-particle self-consistent GW (QPscGW or QSGW) calculations have also been carried out to study the band structure of pyrite, and give apparently satisfactory results (Lehner et al., 2012; Kolb and Kolpak, 2013). However, variants of self-consistent GW without taking the vertex function into account tend to overestimate the band gaps of typical semiconductors, as indicated by several works (Shishkin and Kresse, 2007; Deguchi et al., 2016; Cao et al., 2017; Grumet et al., 2018). Thus the error cancellation between the general tendency of overestimating band gaps of semiconductors and the numerical inaccuracy in the LAPW basis or the use of conventional pseudo-potentials could contribute to the apparent agreement between the generally accepted band gap and the self-consistent GW results. Of course, without looking into computational details of previous self-consistent GW calculations, this is just our speculation. Further investigations are needed to fully clarify this issue. We also note that similar LAPW + HLOs calculation has been conducted for the pyrite phase by Ouarab and Boumaour (2017) and gives a band gap (0.97 eV) close to ours.

3.1.3 Quasi-Particle Band Structure

To further illustrate the significance of HLOs in applying the GW methods to FeS₂, we present the QP band structures of pyrite and marcasite FeS₂ calculated from the G₀W₀@PBE method with the LAPW + HLOs basis set, and compare the results to those with the standard LAPW basis.

Figure 3 shows the electronic band structures of the two FeS₂ phases from different methods. Note that the bands are aligned to the CBM at the Γ point for a better view of QP correction to the valence states. With PBE, pyrite (Figure 3A) is found to be an indirect band gap material with the CBM located at the Γ point and the VBM near the X point along Γ–X. The top valence bands within 1 eV below the VBM are dominated by the localized Fe-3d states, also manifested by their flat dispersion. The dispersive bands about 2 eV below the VBM are mainly composed of S-3p states and well separated from the Fe-3d (t_2g) valence bands. In the conduction band region, the lowermost conduction bands are also largely composed of Fe-3d (e_g), except for the states close to the Γ point with predominant S-3p characters. Particularly, the CBM Γ_c state is exclusively formed by the σ anti-bonding overlapping of S-3p orbitals in the S-S dimer (see projected bands in Figure 4A). Valence and conduction bands with strong Fe-3d characters are separated by about 2 eV. For marcasite, an indirect band gap is also observed, with the CBM located at the T point (T_c) and the VBM along Γ–X (Δ_v). Both states at the VBM and CBM of marcasite are of dominant Fe-3d characters (Figure 4D), in contrast to pyrite where CBM is of pure S-3p characters. The wider Fe-3d valence bands overlap with the S-3p bands near about 1.5 eV below the VBM, which indicates stronger covalent bonding between Fe and S in marcasite than in pyrite.

FIGURE 3. Comparison of band structures computed from PBE (black dotted), G₀W₀@PBE with the standard LAPW basis (blue dashed) and G₀W₀@PBE with the LAPW basis extended by optimized HLOs (LAPW + HLOs, red solid) for (A) pyrite and (B) marcasite FeS₂. The conduction band minimum is aligned as the energy zero and marked by the black dash-dotted line.

FIGURE 4. (A) Projected band structure from GW with the LAPW + HLOs basis, (B) self-energy corrections Δɛ to Kohn-Sham states, and (C) the difference between Δɛ with LAPW + HLOs and LAPW basis against the weight of Fe-d characters $w_{n\mathbf{k}}^{\mathrm{F}\mathrm{e}-d}$defined by Eq. 12 for pyrite FeS₂. GW are performed at the G₀W₀ level. The quasi-particle (A) and Kohn-Sham (B) energies are aligned to the corresponding valence band maximum. In (A), projections of states on Fe-d and S-p orbitals are proportional to the diameters of red and blue circles, respectively. (D–F) are the counterparts for the marcasite phase.

Then we compare the QP band structures obtained from G₀W₀@PBE with the LAPW and LAPW + HLOs basis (Figure 3). With the standard LAPW basis, G₀W₀@PBE predicts pyrite almost as a semimetal with a nearly vanishing band gap (Figure 3A). Dispersion of the conduction band around the Γ point and the separation between the Fe-3d and S-3p valence bands are enhanced compared to the PBE reference. For marcasite (Figure 3B), although a noticeable gap (0.57 eV) is predicted by G₀W₀@PBE, the band edges are different from those in PBE: the CBM is located at the Γ point (Γ_c) and the VBM in the middle of the Z–Γ path (Λ_v). The change in the nature of band edges from semi-local functional to GW method is also observed by Schena et al. (2013). Once HLOs are included in the basis set, QP band gaps of both phases are dramatically enlarged. The fundamental gaps of pyrite and marcasite are 1.04 and 1.15 eV, respectively, which are 0.2 ∼ 0.3 eV larger than the optical gaps from absorption spectra (Sánchez et al., 2016). Band edges of marcasite by GW are also recovered to those by PBE. The comparison indicates that both negative QP corrections to band gaps and change of band edges in GW (LAPW) are indeed artifacts due to the inadequate numerical accuracy of the basis set.

To better understand how the HLOs basis functions influence the QP band structures of FeS₂, we scrutinize the QP correction to Kohn-Sham state Δɛ, defined by the difference between the QP energy ɛ_QP and the KS energy ɛ_KS, i.e. Δɛ ≡ɛ_QP − ɛ_KS. For pyrite, with the standard LAPW basis, Δɛ to the CBM is smaller than those to the valence Fe-3d t_2g and conduction Fe-3d e_g states as shown in Figure 4B. Particularly, Δɛ for the VBM is about 0.7 eV greater than that for the CBM. This leads to a up-shift of Fe-3d states with respect to the CBM on a whole. Extending LAPW with HLOs reduces Δɛ for all states, but the reduction in Δɛ to the VBM is more than that to the CBM by about 1.0 eV, resulting in the sign change of the QP correction to the band gap. Similar conclusion can be drawn from Δɛ in the marcasite phase (Figure 4E). With the standard LAPW method, Δɛ to T_c exceeds that to Γ_c by more than 1.1 eV. Consequently, Γ_c drops down below T_c and becomes the CBM, as we have seen in Figure 3B. Upon including HLOs, Δɛ to T_c is reduced more significantly than Δɛ to Γ_c such that T_c recovers the conduction band edge as in PBE.

Such biased effects of HLOs are clearly associated with the atomic characteristics of Kohn-Sham states, as we have demonstrated in the GW calculations of cuprous and silver halides (Zhang and Jiang, 2019). In Figures 4C,F, we plot the difference between Δɛ computed by G₀W₀ with LAPW + HLOs and LAPW against the weight of Fe-d characters of the Kohn-Sham orbitals $\left|{\psi }_{n\mathbf{k}}\rangle \right.⟩$, $w_{n\mathbf{k}}^{\mathrm{F}\mathrm{e}-d}$, defined by

where ${\varphi }_{l=2,m}^{{\text{Fe}}_i}$ represents the pre-defined atomic function centered on the ith Fe atom featuring spherical harmonic function $Y_2^m$. The negative difference implies that including HLOs generally brings down Δɛ. Moreover, the difference is more dramatic for states with larger $w_{n\mathbf{k}}^{\mathrm{F}\mathrm{e}-d}$, indicating that numerical error is more significant for states with stronger Fe-d characters in GW calculations with the incomplete LAPW basis.

3.1.4 GW Density of States

To end this section, we present the GW calculated density of states (DOS) of FeS₂ polymorphs in Figure 5. The results for pyrite FeS₂ are shown in Figure 5A. Due to different definitions of the Fermi level in theoretical results and experimental spectral data, we have shifted the experimental data to match up the highest valence peak near the Fermi level. With this alignment, the overall DOS from G₀W₀ (LAPW + HLOs) agrees well with the energy distribution curves (EDCs) from the PES experiments. The width of the valence Fe-3d band and separation between the Fe-3d and S-3p valence bands are consistent with the UPS experiment by Mamiya et al. (1997) and the XPS experiment by Folkerts et al. (1987). The location of the first peak in the conduction band region is also in good agreement with the BIS data (Folkerts et al., 1987). Interestingly, although G₀W₀ (LAPW) underestimates the band gap severely, the location of the first peak in the conduction region is almost identical to that by G₀W₀ (LAPW + HLOs), probably due to the error cancellation between QP corrections to the valence and conduction Fe-3d bands. However, such fortuitous cancellation does not hold in the valence region as inferred by the too deep S-3p band in the G₀W₀ (LAPW) results.

FIGURE 5. Density of states (DOS) computed from PBE (black dotted), G₀W₀@PBE with the standard LAPW basis (blue dashed) and G₀W₀@PBE with the LAPW + HLOs basis (red solid) for (A) pyrite and (B) marcasite FeS₂. Energy distribution curves (EDCs) of pyrite extracted from photo-electron spectroscopy (PES) are presented for comparison. Each dataset is normalized with respect to its highest peak. Theoretical data are aligned to its valence band maximum as energy zero. The XPS + BIS and UPS + BIS data for pyrite are obtained from Folkerts et al. (1987) and Mamiya et al. (1997), respectively. To take into account the different definitions of the Fermi level in theory and experiment, a rigid shift of 0.60 and 0.40 eV are employed for the EDCs from XPS + BIS and UPS + BIS, respectively, to match the highest valence peaks below the Fermi level.

Figure 5B shows the calculated DOS for marcasite. Regardless of the theoretical method used, the valence Fe-3d band of marcasite has larger width than that of pyrite, indicating a stronger Fe-S interaction in the marcasite phase. In the conduction region, a sharp peak is observed with the G₀W₀ (LAPW) method, while only a plateau is found with G₀W₀ (LAPW + HLOs). However, the sharp peak is actually an artifact of wrongly pushed up Fe-3d conduction bands due to the inaccuracy of the standard LAPW basis as explained above.

3.2 Results From Hybrid Functionals

As mentioned in the introduction, previous theoretical studies found that various hybrid functionals, which are typically able to describe the band gaps of semiconductors quite accurately, performed badly for FeS₂. In this section, we look into this issue and present results by several hybrid schemes including the DSH functional with system-dependent parameters.

3.2.1 Band Gaps by Hybrid Functionals

Band gaps computed by different hybrid functionals are collected in Table 3. The widely used PBE0 and HSE06 functionals have been reported to predict fundamental gaps of pyrite and marcasite FeS₂ larger than 2 eV in the literature (Sun et al., 2011; Choi et al., 2012; Hu et al., 2012; Schena et al., 2013; Liu et al., 2019), which is confirmed by our results. DSH, the hybrid functional with system-tuned parameters, does not improve the prediction over PBE0 and HSE06. This is surprising, given that DSH has been previously shown to outperform several other hybrids in evaluating band structures for wide- and narrow-gap systems (Cui et al., 2018; Liu et al., 2020), including PBE0, HSE06 and the dielectric-dependent hybrid (DDH) functionals (Marques et al., 2011; Skone et al., 2014). SX-PBE screened exchange functional gives band gaps of FeS₂ significantly smaller than the hybrids mentioned above, but the gaps are still larger than those from GW with the LAPW + HLOs basis (Table 2) by about 0.5 eV.

TABLE 3. Fundamental band gap (indicated by “fund.”) and other direct and indirect band gaps (unit: eV) for pyrite and marcasite FeS₂ calculated by different hybrid functionals. Results from other theoretical studies and experimental measurements are presented as comparison.

Considering that the one-shot DSH, i.e. DSH0, may outperform the self-consistent scheme in some transition metal compounds (Cui et al., 2018; Liu et al., 2020), we also employ DSH0 to calculate the two FeS₂ polymorphs. The macroscopic dielectric constant calculated with PBE ${{\epsilon }_{\mathrm{M}}}^{\text{PBE}}$ is 20.6 for pyrite, which agrees well with ${{\epsilon }_{\mathrm{M}}}^{\text{PBE}}=21$ from a previous study (Choi et al., 2012). DSH0 with ${{\epsilon }_{\mathrm{M}}}^{\text{PBE}}$ predicts smaller band gaps than DSH, but the values are still above 2 eV. Meanwhile, DSH0 with experimentally obtained ɛ_M = 10.9 (Husk and Seehra, 1978) gives the pyrite band gap of 2.72 eV. In contrast, a modified HSE functional (MHSE) with HSE06 screening parameter and 10% hybrid ratio, which is roughly equal to the inverse of the experimental dielectric constant, as suggested by Liu et al. (2019), gives band gaps close to the GW₀ (LAPW + HLOs) result. The MHSE results agree with those by Liu et al. (2019) and seem to verify the suggestion by Schena et al. (2013) of using 1/ɛ_M as the hybrid ratio in the HSE-type screened hybrid functional.

3.2.2 Band Structures by Hybrid Functionals

As summarized above, the investigated hybrid functionals except for MHSE fail to give reasonable predictions for the band gaps of pyrite and marcasite FeS₂. In this section, we take a close look at the band structures computed from these methods to understand the failure.

The band structures for pyrite calculated from selected hybrid functionals are shown in the upper panel of Figure 6. With PBE0 and HSE06 (Figures 6A,B), the fundamental band gap is a direct one with both VBM and CBM located at the Γ point. An indirect fundamental gap is obtained by SX-PBE and DSH (Figures 6C,D), but the VBM is different from that in PBE or the GW method (Figure 3A). In addition, compared to the GW (LAPW + HLOs) results, the separation between valence Fe-3d and S-3p bands is reduced and the splitting between the valence and conduction Fe-3d bands is significantly increased by the hybrid functionals. We note that both features can be understood tentatively as a result of increased ligand field strength from the perspective of ligand field theory. This indicates an overestimated interaction between the ligand S-3pσ and Fe-3d orbitals in the selected hybrid functionals than that in PBE. The overestimation is most significant in the DSH method (Figure 6D), where the state of predominant S-3pπ characters along the M–Γ path becomes the VBM and conduction Fe-3d bands are raised beyond 6 eV above the Fermi level.

FIGURE 6. Band structures of pyrite (upper panel) and marcasite (lower panel) FeS₂ computed from different hybrid functional methods. From left to right, the methods used are PBE0 (A,E), HSE06 (B,F), SX-PBE (C,G) and DSH (D,H), respectively. The valence band maximum is aligned as the energy zero and indicated by the black dash-dotted line.

We can observe similar features in marcasite band structures from hybrid functionals, as shown in the lower panel of Figure 6. In the valence band region, the S-3p bands are pushed up relatively to Fe-3d bands compared with PBE and GW. The increase is so significant that the VBM along Γ–X, which is mainly of Fe-3d in PBE and GW, is now of predominant S-3p characters. This also leads to a considerable overlap between the two sets of bands in the energy window 1 ∼ 3 eV below the Fermi level. The conduction bands are also shifted to higher energies. However, the shifts are larger for the conduction Fe-3d bands than for S-3p. For the DSH method as an extreme case, the Fe-3d bands are raised up too high and even separated from the S-3p bands.

The radical failure of DSH invites a close inspection of feasibility of DSH for FeS₂. As a preliminary exploration to the possible cause, we make a direct comparison between the inverse static dielectric function used in the DSH with ${{\epsilon }_{\mathrm{M}}}^{\text{PBE}}$ and that from the RPA calculation with the LAPW + HLOs basis in pyrite FeS₂ as a function of the length of wave vector in the long-range limit, i.e. q → 0. The inverse dielectric function corresponding to DSH reads (Liu et al., 2020)

We note that Eq. 13 differs from the inverse of Eq. 9 because in the derivation of DSH, the exponential function is replaced by erfc [Cui et al. (2018) for more details]. As shown in Figure 7, while DSH overestimates ɛ⁻¹ and underestimates the screening in the short-wavelength region, i.e. near |G| = 0, DSH0 model dielectric function with ${{\epsilon }_{\mathrm{M}}}^{\text{PBE}}$ closely resembles that from RPA calculation, which is similar to the observation by Liu et al. (2020) in transition metal oxides. Hence we consider that the screening effect is reasonably captured in DSH0. Further investigation is needed to understand the cause for the failure of DSH for FeS₂.

FIGURE 7. Inverse dielectric functions of pyrite used in the DSH model and calculated from RPA@PBE with the LAPW + HLOs basis set. Dielectric constant computed by PBE is used in the one-shot DSH0 method.

4 Conclusion

In the present study, we have investigated the electronic band structures of two FeS₂ polymorphs, namely pyrite and marcasite, by using methods in different frameworks. With the all-electron many-body GW method implemented in the LAPW framework, we find that by using GW₀@PBE with the LAPW + HLOs basis, pyrite and marcasite are predicted to have indirect fundamental band gaps of 1.14 and 1.16 eV, respectively. The closeness of band gaps for the two polymorphs agrees with the experimental observation (Sánchez et al., 2016). The pyrite band gap from GW₀@PBE with LAPW + HLOs is very close to the generally accepted experimental value (Ennaoui et al., 1993) and the corresponding density of states also agrees well with energy distribution curves obtained from the photoelectron spectroscopy measurements (Folkerts et al., 1987; Mamiya et al., 1997). In contrast, with the standard LAPW basis, PBE-based G₀W₀ and GW₀ both lead to negative QP correction to the PBE fundamental gap, which is rarely observed in LDA/GGA-based G₀W₀ and GW₀ treatments of semiconductors. The splitting between Fe-3d and S-3p valence bands of pyrite is also significantly overestimated compared to experiment. These artifacts exist not only in calculations with the standard LAPW basis, but also in those with LAPW basis extended by an inadequately converged HLOs (Schena et al., 2013). Therefore in order to eliminate such artifacts, it is instrumental to carefully converge the fundamental band gap with respect to the two controlling parameters, namely n_LO and Δl_LO. We have further studied electronic band structures of FeS₂ polymorphs with different hybrid functionals, including PBE0, HSE06, the screened exchange SX-PBE and the recently developed DSH functional with system-tuned hybridization parameters. We find that all those methods overestimate the band gaps of the two polymorphs by 0.5 ∼ 1.9 eV compared to the results obtained from G₀W₀ (LAPW + HLOs). The overestimation by PBE0 and HSE06 as reported in the literature is reproduced in this work. Furthermore, either self-consistent or one-shot DSH method fails to improve over the conventional fixed-parameter hybrid functionals. By comparing the model dielectric function used in DSH with that from RPA calculation with LAPW + HLOs in pyrite, we point out that the failure of DSH may not be caused by the insufficiency of the dielectric model used and therefore requires further investigation. Our investigations clearly show that accurate prediction of electronic band structures of FeS₂ polymorphs poses a stringent test on the state-of-the-art first-principles approaches, and the GW method based on semi-local density approximation performs well for this difficult system if it is practiced with well-converged numerical accuracy.

Finally, we note that further work in the following aspects can be done to shed more light onto the band gap problem of FeS₂ in terms of GW and hybrid functional calculations. For one thing, it is possible to build the screened Coulomb interaction W using the KS states from the LAPW calculations and calculate the self-energy Σ with G from LAPW + HLOs. One can compare it with GW using LAPW to see whether it is the inaccurate band summation in W or G to blame. For another, replacing the PBE with the hybrid functional as starting point will be worthwhile to evaluate the dependence of G₀W₀/GW₀ results on initial input for FeS₂. Particularly, considering the severe overestimation of FeS₂ band gaps by the hybrid functionals, it is of great interest to see whether G₀W₀/GW₀ can produce a negative QP correction to the gap from hybrid functional calculations such that the experimental gap is approached from above.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

M-YZ performed the calculations, conducted the analysis and wrote the manuscript. HJ revised the manuscript and supervised all the work. All authors listed contributed to the article and approved it for publication.

Funding

National Natural Science Foundation of China (grant numbers: 21873005 and 21911530231).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors acknowledge the support by High-performance Computing Platform of Peking University for the computational resources.

References