BaZrSe$3$⁠: Ab initio study of anion substitution for bandgap tuning in a chalcogenide material

Transition metal perovskite chalcogenide materials have been recently proposed as candidates for solar cells and other optoelectronic applications. Unlike traditional semiconductors, based on monatomic (group IV) or diatomic (III–V or II–VI) systems, the large number of possible perovskites opens up a large potential for tunability, making them promising for various applications. Recently, lead halide perovskites of the form APbX$3$⁠, where A is a monovalent cation and X is a halogen, have caused a boom in the solar cell research community.^1–4 Since the proposal of these materials in 2009, power conversion efficiencies of solar cells built from these materials have risen up to 22.7% in laboratory conditions, comparable to the performance of state-of-the-art commercial devices, such as silicon solar cells. Due to facile synthesis and high abundance of precursors, devices made from perovskite materials may be produced with low cost compared to conventional technologies such as Si or CdTe.^5–7 However, instability, especially under conditions of high humidity, and environmental concerns about Pb have hindered the commercial use and application of these materials.

In efforts to find materials that solve these practical concerns while keeping the advantages that come with perovskite materials, numerous efforts have been made to search for new perovskite materials. Perovskite oxide materials are expected to be more resilient in the presence of moisture. Unfortunately, their bandgaps are too high for most solar cell applications; the perovskite oxide with the lowest known bandgap is BiFeO$3$⁠, with a bandgap of 2.7 eV.⁸ Various studies^8–11 have shown that the bandgaps of these oxides can be tuned by doping and alloying with different oxides, but the high intrinsic disorder of these compounds may cause problems with stability and interfacial matching within a photovoltaic cell. Thus, it is of interest to substitute O with another chalcogenide to reduce the bandgap.

Sun et al.¹² calculated bandgaps of a wide variety of transition metal chalcogenides (TMC) of the form ABX$3$ (A = Ca, Sr, Ba; B = Ti, Zr, Hf; X = S, Se) assuming a set crystal geometry. The high relative abundance of these elements suggests that production of these materials may be possible at low cost. Their results show that the substitution of O with S or Se lowers the bandgap into a more feasible range for solar cell applications. According to their analysis, four Se-containing compounds, namely, CaZrSe$3$⁠, SrZrSe$3$⁠, BaZrSe$3$⁠, and BaHfSe$3$⁠, were found to have bandgaps corresponding to theoretical power conversion efficiencies larger than 28 percent.

Until recently, there has been little experimental characterization of the optoelectronic properties of TMC materials. Niu et al.¹³ synthesized and characterized BaZrS$3$ as well as two phases of SrZrS$3$⁠: the “needlelike” $α$-phase and the “distorted perovskite” $β$-phase. While the bandgaps obtained for these compounds (1.81 eV, 1.53 eV, and 2.53 eV) were not optimal, stability at high temperatures and large absorption coefficients were observed.

In this paper, we investigate the substitution of S with Se in the crystal BaZrS$3$ to form BaZrSe$3$⁠. This is done with the hope that this material will have attractive properties. To examine this material, we use evolutionary methods to find its ground state structure. We then demonstrate the stability of this structure using two methods: calculation of the phonon dispersion curves and ab initio molecular dynamics (AIMD) simulations. We examine the electronic properties of the synthesized TMC materials EuZrSe$3$⁠, $α$-SrZrS$3$⁠, $β$-SrZrS$3$⁠, and BaZrS$3$ and determine the effect of the chalcogenide substitution on the electronic and optical properties of BaZrSe$3$⁠. This will allow us to evaluate its feasibility as a light absorber for solar cells.

Crystal structures of $α$-SrZrS$3$⁠, $β$-SrZrS$3$⁠, and BaZrS$3$ are taken directly from X-ray crystallography measurements of Niu et al.¹³ The crystal structure of EuZrSe$3$ was taken from experimentally determined values from the Inorganic Crystal Structure Database. The crystal structure for BaZrSe$3$ was predicted using evolutionary methods, as implemented in the code USPEX.^14–16 The USPEX code, developed by Oganov et al., uses local structure optimization, real-space representation, and variational operators designed to mimic natural evolution.

To begin the evolutionary search, an initial population of 200 structures with randomly chosen space groups is generated. The total energy is calculated via density functional theory (DFT) using Vienna ab initio Simulation Package (VASP).^17,18 Optimization of these structures occurs in four steps that turn gradually finer, the kinetic energy cutoff for the plane wave expansion in the final step being 400 eV. The optimized structures form the first generation. Subsequent generations containing 40 structures are then formed. Some of these are modifications of the best structures (those with lowest energy) of the previous generation, while others are randomly generated. This process is repeated until a convergence to a clear minimum is achieved.

To verify the zero-temperature stability of the most stable crystal structure of BaZrSe$3$ predicted by USPEX, we constructed a convex hull of energies using the predicted BaZrSe$3$ compound and all compounds containing Ba, Zr, or Se present in the Materials Project database.¹⁹ The energy at each composition on the convex hull is such that it is of lowest energy compared to any other phase or linear combination of phases at that composition. To ensure that the energies of different compositions can be directly compared, the same set of calculation parameters are used to perform a VASP total energy calculation on all the compounds considered for the construction of the convex hull. Note that since this convex hull includes the new hypothetical BaZrSe$3$ structure, a value of zero for the energy of BaZrSe$3$ with respect to the convex hull would indicate stability with respect to the compounds present in the database.

Phonon calculations on the most stable structure of BaZrSe$3$⁠, as predicted by USPEX, were performed using the Quantum-Espresso distribution for materials simulation.²⁰ The initial electronic structure was calculated with a uniform Monkhorst–Pack grid of $k$-points of size 8 $×$ 8 $×$ 8 with fixed electronic occupations. We employ the Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional²¹ with ultrasoft pseudized atomic cores from the GBRV (Garrity-Bennett-Rabe-Vanderbilt) repository.²² We select wavefunction and charge density kinetic energy cutoffs of 680 eV and 5400 eV, respectively. The phonon frequencies were then calculated using density-functional perturbation theory (DFPT)^23–25 on a 4 $×$ 4 $×$ 4 phonon-momentum mesh.

The ab initio molecular dynamics (AIMD) simulations were carried out using the VASP simulation package which implements Born-Oppenheimer AIMD based on DFT. We use projector-augmented wave (PAW)^26,27 pseudopotentials for the description of the electron-ion-core interactions. The approach implemented in VASP is based on exact DFT-based evaluation of the instantaneous electronic ground state at finite temperature (with a free energy as variational quantity) at each MD step. The electron exchange-correlation potential is calculated within generalized gradient approximation with PBE potential (GGA-PBE). The $Γ$ point of the supercell is used to expand the wave functions with kinetic energy cutoff of 500 eV. A time step of 2 fs is used for the integration of the equations of motion, and the electronic and ionic temperatures are controlled via a Nose–Hoover thermostat.^28,29

The simulated supercell is composed of 160 atoms corresponding to a $2×4×1$ supercell of the relaxed predicted structure. The structures are thermalized at 200 K, 300 K, 400 K, and 500 K for about 20 ps under isothermal-isochoric conditions (⁠$NVT$⁠) to obtain well equilibrated hot crystals. Data for the calculation of radial distribution functions are gathered and averaged over the last 4 ps of the simulation.

Electronic band structures and densities of states were determined using the full-potential linearized augmented plane wave method, as implemented in the WIEN2k code.³⁰ Here, space is divided into two regions: an interior region, consisting of nonoverlapping muffin-tin spheres centered on each atom, and an interstitial region, consisting of the space between the spheres. In the interior region, the electronic wave function is atomiclike, expanded in spherical harmonics up to $ℓmax=10$⁠. In the interstitial region, it is expanded using a plane wave basis set with a maximum wave vector of $kmax$⁠, chosen such that $Rmtkmax=9$⁠, where $Rmt$ is the radius of the smallest muffin-tin sphere in the unit cell. For the Ba atom, the $6s$ orbital constitutes the valence states, while the $5s$ and $5p$ orbitals constitute the semicore states. For the Sr atom, the $5s$ orbital constitutes the valence states, while the $4s$ and $4p$ orbitals constitute the semicore states. For the Eu atom, the $4f$ and $6s$ orbitals constitute the valence states, while the $5s$⁠, $5p$⁠, and $5d$ orbitals constitute the semicore states. For the Zr atom, the $4d$ and $5s$ orbitals constitute the valence states, while the $4s$ and $4p$ orbitals constitute the semicore states. For the S atom, the $3s$ and $3p$ orbitals constitute the valence states, and all others are considered core states. For the Se atom, the $3d$⁠, $4s$⁠, and $4p$ orbitals constitute the valence states, and all others are considered core states. The charge density is Fourier-expanded with a maximum wave vector of $14/a0$⁠, where $a0$ is the Bohr radius. The HSE06 hybrid functional, which is known to produce accurate results for the bandgaps of semiconductors, was adopted in our calculations.

The optical properties of BaZrSe$3$ were obtained by using the wave functions from the electronic structure calculation with the HSE06 hybrid functional. In particular, the imaginary part of the dielectric function, $ε″(ω)$⁠, was obtained using the momentum matrix elements,³¹ 

Here, $f$ is the Fermi–Dirac distribution function, $ℏω$ is the incident photon energy, $|kn⟩$ is a crystal wave function with energy eigenvalue $Ekn$⁠, and $p$ is the momentum operator. The real part of the dielectric function, $ε′(ω)$⁠, can be obtained using a Kramers–Kronig transformation. The absorption coefficient $I(ω)$ is then calculated as

As an accurate bandgap is essential for the calculation of optical and electronic properties, we also perform calculations with the $GW0$ method to further check the bandgaps of these materials. $GW0$ calculations were performed using VASP within the random phase approximation. Here, the self-energy is approximated as a product of the Green’s function $G$ and a screened interaction term $W0$⁠. The quasiparticle energies for the Green’s function were obtained via four iterations of a self-consistency loop. The quasiparticle energies within the $GW0$ step were calculated using a total of 1000 bands and a $k$-point mesh of size $4×4×4$⁠.

The crystal structures of BaZrSe$3$ with lowest energy resulting from the evolutionary search are shown in Table I. The most stable structure was found to be orthorhombic with space group (Pnma). Although this shares the same space group as the compound BaZrS$3$⁠, the selenide adopts a needlelike phase, as depicted in Figure 1. The next most stable structure, with a total energy 10 meV/atom higher than the ground state structure, was also orthorhombic, with space group (Pnm2$1$⁠). This is a slightly expanded and distorted form of the structure of BaZrS$3$⁠, which adopts a distorted perovskite phase. Finally, the third most stable structure was monoclinic with space group (P2$1$/m). All of these structures possess inversion symmetry.

The calculated energy of the most stable structure of BaZrSe$3$⁠, as predicted by USPEX, falls on the convex hull, which demonstrates that the compound is stable with respect to the compounds present in the Materials Project database. Figure 2 shows the phase diagram. Here, the blue dot represents BaZrSe$3$⁠, while each green dot represents a stable compound from the Materials Project. These compounds were used in the calculation of the convex hull energy.

The phonon dispersion curves for BaZrSe$3$ are shown in Fig. 3. The lack of imaginary frequencies suggests that this particular structure of this compound is stable at zero temperature.

To investigate the dynamical stability of this crystal at nonzero temperatures, we employ ab initio molecular dynamics (AIMD). With this method, a system can be pushed away from a local minimum by the introduction of temperature, thereby making it possible to explore the stability of the crystal structure at finite temperatures.

To this end, we thermalized the predicted zero-temperature crystal structure of BaZrSe$3$ at different temperatures ranging from 200 to 500 K and analyzed the partial radial distribution function for the Ba-Se and Zr-Se pairs.

To characterize the structural stability of the BaZrSe$3$ hot crystals, we computed the partial radial distribution functions (pRDFs). In Fig. 4, we show the partial $g(r)$ for the Ba-Se (a) and Zr-Se (b) for the BaZrSe$3$ crystal at temperatures ranging from 200 to 500 K. For the hot crystals, the pRDF is averaged over several structures taken over the last 4 ps of the AIMD simulation.

The radial distribution functions of the hot crystals reveal correlation peaks at the same exact locations as the pRDF computed for the 0 K crystal. This indicates that the Ba-Se and Zr-Se shells do not change significantly due to the dynamics of the system at different temperatures. It is worth noting that the smearing of the peaks is due to the averaging over several structures.

Figure 5 shows snapshots of the BaZrSe$3$ crystal at 0, 200, 300, 400, and 500 K along the $y$-axis. From the snapshots we observed that same structural units are preserved regardless of the temperature, that is, each Ba, Zr, and Se atom has the same chemical environment even at the highest temperature.

In Table II, we present our calculations for the bandgap of these TMC materials using different methods. We note that while experimental values reflect the values of the optical gaps, the values reported from the calculations are those of the electronic gaps. However, due to weak excitonic effects in most inorganic semiconductors, the electronic and optical gaps are very close. Thus, when comparing our results to experiment, we ignore this distinction. A low bandgap is essential for a light absorber in solar cells, so we must take caution to verify the accuracy of the methods used to evaluate the bandgap. When using methods like GGA–PBE, the bandgap tends to be underestimated, as observed by the calculations we performed with the known sulfide compounds. There are several approaches that have been used to accurately compute the bandgap, but are they more computationally expensive. The $GW0$ method tended to overestimate the bandgap in our case. The HSE06 hybrid functional gives the results that agree most with experiment, achieving errors of 16%, 6%, and 2%, for $α$-SrZrS$3$⁠, $β$-SrZrS$3$⁠, and BaZrS$3$⁠, respectively. Hence, we adopt the HSE06 hybrid functional for calculating the bandgap of EuZrSe$3$ and BaZrSe$3$⁠, where no experimental measurements are available. The bandgap of EuZrSe$3$ is calculated to be around 0.6–0.7 eV, suggesting that this material is a narrow-gap semiconductor.

In Figs. 6 and 7, we present the electronic band structures and densities of states for BaZrS$3$ and BaZrSe$3$⁠, respectively. Because both sulfur and selenium are significantly less electronegative than oxygen, the bandgaps of both BaZrS$3$ and BaZrSe$3$ are much lower than that of the oxygen analog BaZrO$3$⁠, which has a bandgap of 4.8 eV.³² In both systems, the main contribution to the valence band is derived from the $p$ orbitals on the chalcogen atoms, and the main contribution to the conduction band is derived from the $d$ orbitals on Zr atoms. As indicated in Table II and Figs. 6 and 7, the bandgap has been significantly lowered in BaZrSe$3$ to a predicted value of 1.11 eV, a value that makes it a candidate as a light absorber in solar cells.

The optical absorption of BaZrSe$3$ is shown in Fig. 8. From the plot, the absorption coefficient is seen to quickly exceed $105$ cm$−1$ after surpassing the value of the bandgap at 1.11 eV. This behavior of rapid increase of the absorption coefficient at energies above the bandgap is comparable to that of BaZrS$3$¹³ and indicates potential for facile light harvesting. Additionally, we used the method of Tauc et al.³³ to estimate the bandgap from the optical absorption curve. They show that the relation

holds true in a particular regime. Here, $hν$ is the incident photon energy, $Egap$ is the bandgap energy, $α$ is the absorption coefficient, and A is some constant. For direct allowed transitions, as in the case of BaZrSe$3$⁠, $n$ is taken as $1/2$⁠. Plotting $(αhν)1/n$ against $hν$ and then extrapolating the linear regime to the $x$-intercept give the bandgap. Using this method, we obtain a bandgap of 1.14 eV, compared with the calculated value of 1.11 eV. The excellent agreement between these two numbers is indicative of the validity of Eq. (3) for the band structure of BaZrSe$3$⁠.

Based on our calculations, we demonstrate that the HSE06 hybrid functional performs with reasonable accuracy on experimentally measured TMC materials. We conclude that substitution of S with Se forms a compound that is predicted to be stable based on phonon frequency calculations and ab initio molecular dynamics simulations. Large optical absorption coefficients make BaZrSe$3$ a promising material as a light absorber for solar cell applications. It is hoped that this work will stimulate experimental investigation of the compound BaZrSe$3$⁠.

M.O. and R.J. acknowledge support from the National Science Foundation (NSF) under CREST Grant No. HRD-1547723 and PREM Grant No. DMR-1523588. Q.C. and I.D. acknowledge support from the NSF under CAREER Grant No. DMR-1654625.