Computational study of structural, elastic, electronic, phonon dispersion relation and thermodynamic properties of orthorhombic CaZrS$_3$ for optoelectronic applications

Chalcogenide perovskites offer superior thermal and aqueous stability as well as a benign elemental composition compared to organic halide perovskites for optoelectronic applications. In this study, the structural, electrical, elastic, phonon dispersion, and thermodynamic features of the orthorhombic phase of chalcogenide perovskite CaZrS$_3$ (space group Pnma) were examined by first principles calculations utilizing the plane wave pseudopotentials (PW-PPs) in generalized gradient approximations (GGA). The ground state properties such as lattice parameters, unit cell volume, bulk modulus, and its derivative were calculated and are in a good agreement with existing findings. The mechanical properties such as bulk modulus, shear modulus, Young's modulus and elastic anisotropy were calculated from the obtained elastic constants. The ratio of bulk modulus to shear modulus confirms that the orthorhombic phase of CaZrS$_3$ is a ductile material. The absence of negative frequencies in phonon dispersion curve and the phonon density of states give an indication that the structure is dynamically stable. Finally, thermodynamic parameters such as free energy, entropy, and heat capacity were calculated with variation in temperature. The estimated findings follow the same pattern as previous efforts.


Introduction
The quest for technological advancement, particularly in the field of semiconductors, plays a significant role in several optoelectronic, photonic, and energy technologies. Among current semiconductors, the prevailing materials like silicon, group III-V, and group II-VI are typically constructed by a fourfold coordinated tetrahedral network of covalent bonds. There have been major successes in developing solar cell semiconductor materials such as Si, GaAs, CuIn Ga 1− Se 2 (CIGS), and lead halide perovskite based materials [1][2][3]. Over the past decades, inorganic-organic/organometal lead halides have been extensively studied since the early 20th century [4,5]. In 1970, for the first time, Weber reported the synthesis and physical properties of −CH 3 NH 3 PbX 3 − (X = Cl, Br, I) organometal lead halide perovskite [4]. During this time, the organometal lead halide perovskites emerged as a candidate and very promising materials for light harvesting in the solar cell as reported in 2009 [5]. The latest area in solar cell materials is the organic-inorganic hybrid lead (Pb) halide perovskites, for which the efficiency reached 22.7% in 2018 [3], starting from 3.8% in 2009 [5].These hybrid perovskites, however, have poor thermal and moisture stability, as well as the presence of lead toxicity [6]. The intrinsic poor long-term stability of CH 3 NH 3 PbI 3 -based perovskite solar cells, as well as the presence of toxicity, has hindered their industrial application and commercialization. Inorganic lead-free perovskites have recently been widely investigated to address these issues. Recently, chalcogenide perovskites have received attention due to their promising photovoltaic and thermoelectric properties, with initial studies conducted on oxide perovskites that have good band gaps for optical absorption [7]. Chalcogenide perovskites assume an ABX 3 configuration with A, and B are elements with a combined valence of 6 (with different valences), while X is typically S or Se. These materials belong to a new class of ionic semiconductors. The band gap of these materials can be systematically tuned in a wide range from ultra-violet to infrared. Due to their predicted strong iconicity, they may exhibit unique physical properties such as being free of deep-level defects, which is beneficial for energy harvesting and other optoelectronic applications [8]. It should be emphasized that oxide perovskites, with a chemical formula ABO 3 , have long been a focus of active research. This family of materials exhibits unusually rich properties ranging from colossal magnetoresistance, ferroelectricity to superconductivity and charge density waves, resulting from the interplay of different degrees of freedom with similar energy scales. The intriguing physics of their chalcogenide counterparts, however, is largely unexplored [8].
In recent years, theoretical calculations based on the density functional theory (DFT) have been used to reveal and predict the structural, mechanical, electrical, optical, and thermal properties of crystal materials. CaZrS 3 , which is the focus of this work, is a family of chalcogenide perovskite crystal material that was considered theoretically for optoelectronic applications [8][9][10][11][12]. The relaxed lattice parameters and band gap of CaZrS 3 were estimated using DFT as implemented in VASP and the Perdew, Burk, and Ernzerhof (GGA-PBE) generalized gradient approximation [13]. In addition, the lattice parameter and the band gap of CaZrS 3 were also calculated using the FPLAPW method with DFT as implemented in WIEN2K, approximating the exchange correlation potential with PBE-GGA, Engel-Vosko (EV) method and Hubbard parameter (GGA+U) [12]. However, to our knowledge, the structural, elastic, electronic, phonon dispersion and thermodynamic properties of CaZrS 3 are not yet well investigated for optoelectronic application. Moreover, investigation of elastic, phonon dispersion and thermodynamic properties of CaZrS 3 using the first principle computational methods remains unexplored. In this paper, the structural, elastic, electronic, and phonon dispersion relation and thermal properties of CaZrS 3 are carefully examined. The electronic properties are calculated by considering PBE-GGA [14] and also DFT with the Hubbard functional (DFT+U) [15] for exchange correlation potential using Quantum ESPRESSO package (QE). In addition, phonon dispersion and the thermodynamic properties of CaZrS 3 are studied using a 1×1×2 (in , , and direction, respectively) supercell containing 40 atoms created in a PHONOPY package [16].

Computational methods
In this study, the DFT as implemented in the QE [17] within the generalized gradient approximation (GGA) functional [14] and with the Hubbard correction (DFT+U) [15] was used. The effective Hubbard parameter ( eff ) was calculated iteratively for Zr-orbitals. For this study, a cell with 20 atoms (4-Ca, 4-Zr, and 12-S) in orthorhombic phase for structural, elastic and electronic property calculations was used. The ultra-soft pseudopotentials (US-PP) were used to treat the interaction of the electrons with the ion cores as in [18]. The corresponding valence electrons considered for the calculations are Ca -[Ar]4 2 , Zr -[Kr]4 2 5 2 , and S -[Ne]3 2 3 4 . Crystal structure optimization was done using a plane wave cutoff energy of 60 Ry and the Brillouin zone with a 3 × 3 × 3 Monkhorst-Pack -point grid [19] based on the convergence criteria energy 10 −4 Ry, force 10 −3 Ry/Bohr, and cell pressure 0.5 kbar. Using the optimized structure, the elastic properties were calculated by THERMO_PW package within Quantum ESPRESSO package [17]. In addition, the thermal properties were calculated with the help of PHONOPY package [16]. To study the phonon dispersion relation and thermodynamic properties of CaZrS 3 , a supercell of 1 × 1 × 2 in , , and direction with 40 atoms was created and used for computations.

Crystal structure
Most of the material properties are governed by their crystal structures. A stable crystal structure is the one with the lowest energy arrangement of atoms at a given temperature and pressure [20]. In this study, an orthorhombic phase of CaZrS 3 [space group (62) point group 2ℎ ( )] with the crystallographic structure of GdFeO 3 -type was considered. Each Ca atom is centrally located in the spatial region defined by its neighbouring S and Zr atoms in the lattice, which is composed of distorted ZrS 6 corner-sharing octahedral as visualized in figure 1. The ionic components are Ca +2 , Zr +4 and S −2 ions. The structural stability of perovskite materials in general is determined by the Goldschmidt tolerance factor ( ) [13], where Ca + , Zr − and S − are ionic radii for Ca +2 , Zr +4 and S −2 ions, respectively. The sizes of the ions are known to have a major influence on the structural distortion of perovskites. Thus, materials with a tolerance factor of 0.71 ⩽ ⩽ 0.9 result in a distorted perovskite structure with tilted octahedral, while for 0.9 < < 1.0 the materials have an ideal cubic structure, and when the tolerance factor is much higher ( > 1) or lower ( < 0.71), non-perovskite structures are commonly formed [21,22]. The ionic radii of Ca, Zr and S, and the calculated tolerance factor CaZrS 3 are shown in table 1. Here, since the calculated value falls within a distorted perovskite range, CaZrS 3 stability is defined, and Pnma symmetry is adopted. from the convergence test were utilized. Structural optimization was performed setting the convergence criteria; change in energy, Δ = 1.0×10 −4 Ry and change in force, Δ = 1.0×10 −3 Ry/Å. The optimized equilibrium lattice constants were = 6.57, = 7.06, and = 9.63. These values are in good agreement with the theoretical and experimental values, as summarized in table 2. Moreover, a series of strained lattices were used to calculate the static lattice potential corresponding to total energy. From these results, the equilibrium unit cell volume, bulk modulus, and its pressure derivative can be calculated. A series of total energy calculations as a function of volume can be fitted to an equation of state (EOS) according to Murnaghan [23]: where 0 is an equilibrium bulk modulus that effectively measures the curvature of the energy versus volume curve about the relaxed volume 0 , and ′ 0 is the derivative of the bulk modulus. The calculated values of the bulk modulus, equilibrium unit cell volume and the dimensionless bulk modulus derivative of CaZrS 3 are given in table 2. The calculated values of the unit cell volume and bulk modulus are in a good agreement with the experimental and the previous theoretical values, respectively as shown in table 2.

Mechanical properties
The elastic constant of crystals gives fundamental information for the study of mechanical characteristics of materials, as they are related to the mechanical properties of the material such as the elastic moduli, Poisson's ratio, and elastic anisotropy factor of materials. To calculate the elastic constants, we applied the non-volume-conserving method. The complete elastic constant tensor was determined from calculations of the stresses induced by small deformations of the equilibrium primitive cell. The elastic constant tensors are given by [25][26][27]; where stands for the Helmholtz free energy, and are the applied stress and Eulerian strain tensors, and stands for the coordinates.

23701-4
In this case, for orthorhombic system, there are nine independent elastic constants that should satisfy the well-known Born stability criteria [28].
The calculated elastic constants (table 3) satisfy the mechanical stability conditions above and > 0, 12 < 11 and 13 < 1/2( 11 + 33 ). Using the Voigt-Reuss-Hill (VRH) average approximation, mechanical parameters such as the Young modulus ( ), Poisson's ratio ( ) and shear modulus ( ) are computed from the calculated elastic constants as [28,29], the elastic constants and calculated mechanical properties from the elastic constants are shown in table 3 and table 4, respectively. Bulk modulus is an essential physical parameter in describing the compressibility of solids under the hydrostatic pressure. Large value of bulk modulus results in higher compressibility of a solid material. The Bulk modulus value (76.36 GPa) calculated from elastic constants is found to be close to the value calculated from the equation of state (81.8 GPa) as in table 2. The shear modulus is another important parameter that can describe the shape change under the shear force. The larger the shear modulus is, the higher is the shape change resistance of the solid material. The calculated value of the shear modulus is 42.73 GPa and it is comparable to the value obtained by [10]. Furthermore, determining the ratio of / is important for understanding the brittle and ductile behavior of materials in the material fabrication. The ductility and brittleness of materials can be determined based on the value of / ratio according to Pugh [30]. The cutoff value is 1.75. When / > 1.75, the material behaves in a ductile manner, otherwise, it exhibits brittle properties. From table 4, our calculated ratio for / value is 1.78 for CaZrS 3 in GdFeO 3 -type phase, showing that CaZrS 3 in this phase is ductile. Another mechanical parameter that provides information about the feature of the bonding forces is Poisson's ratio ( ). In the evaluation of Poisson's ratio, the values 0.25 and 0.5 are the lower and upper limits of the central force, respectively [31]. From table 4, the calculated Poisson's ratio ( ) is 0.26 (which is between 0.25 and 0.5) indicating that the inter-atomic forces are central.
The universal anisotropic index [ = (5 / ) + ( / ) − 6] is a measure to define the elastic anisotropic or isotropic characteristics based on the contributions of both bulk and shear modulus [32]. It is one of the important physical parameters used to study the service life time of materials. The material is isotropic if the value of = 0; otherwise it ( ≠ 0) refers to the anisotropic mechanical properties. Any value smaller or greater than zero represents a higher extent of anisotropy. Based on this, for orthorhombic phase of CaZrS 3 , the calculated value for is 0.35, indicating that the CaZrS 3 was found to be anisotropic.

Debye temperature
The thermal properties of a solid material are related to two physical parameters: Debye temperature ( D ) and melting temperature , respectively. The Debye temperature is another essential physical term that can be used to characterize solid-state physics phenomena such as lattice vibration, elastic constants, specific heat, and melting point. The magnitude of the Debye temperature is helpful to know the thermal conductivity of solid materials. The higher the value of the Debye temperature is, the higher is its thermal conductivity. The Debye temperature ( D ) of CaZrS 3 can be estimated from the averaged sound velocity, , given by [33], where ℎ is Plank's constant, B is Boltzmann's constant, A is Avogadro's number, is density, is molecular weight, is the number of atoms in a formula unit, is the longitudinal sound velocity, and is the transverse sound velocity. The longitudinal and transverse sound velocities can be obtained from density, shear and bulk modulus of the material as: Moreover, the Debye average sound velocity which represents the maximum frequency of the material is described by D = B / 1/2 . The melting point of a material depends on Debye temperature; a larger Debye temperature of the material shows a higher melting temperature [34]. The melting temperature of a stable phase of CaZrS 3 can be determined based on an elastic constant 11 using [35], The calculated values of longitudinal sound velocity, transverse sound velocity and Debye temperature for CaZrS 3 are given in table 5. It was observed that the calculated values for Debye temperature and melting point are 415.5 K, and 1267.6 K, respectively.

Electronic properties: density of states (DOS) and band structure
The energy band structure and density of states of materials are used to determine the electronic properties of solid materials. The accessible electronic energy levels of solid materials are represented by electronic band structures. In this study, the electronic band structure along the high symmetry direction of the Brillouin zone was estimated using the GGA-PBE functional and the Hubbard correction (GGA+U) for exchange correlation potential, as shown in figure 2. The GGA-PBE functional fails to approximate the exact exchange correlation potential, since the band gap value obtained by approximating the exchange correlation potential with GGA-PBE functional is 1.23 eV (table 6), which appears to be underestimated as compared to experimental band gap values of 1.90 eV [8,24]. Furthermore, the band gap was calculated with the Hubbard correction for on-site interaction, yielding a band gap of 1.88 eV, which is close to the experimental result [8,24].
The DOS is also used to describe how state occupancy behaves at different energy levels. It provides information on both occupied and empty states. The states that are available for occupancy have a high  1.90 DOS at a certain energy level. However, there is no state occupied at DOS equal to zero. In this study, the total and partial densities of states were obtained for the equilibrium states of the phases using GGA-PBE correlation interaction and also with GGA+U as shown in figure 3 and figure 4. The density of states is discontinuous for the width from the top of the valence band to the bottom of the conduction band which is normally refered to the band gap of the system. Moreover, figure 4 shows that the maximum valence band is mainly contributed by S-2 orbitals and the minimum conduction band is mainly dominated by Zr-3 orbitals. On the other hand, Ca-4 orbitals are observed on both maximum valence and minimum conduction bands, and the rest of the orbitals have a small contribution.

Phonon dispersion relation
Phonon vibration plays an essential role in dynamic behaviors and in thermal properties, which are central topics in fundamental issues of materials science. The phonon frequency of crystalline structures is one of the fundamental aspects when considering the phase stability, phase transformations, and thermodynamics of these materials. The phonon density of states ( ) is given by [36,37] where is the number of unit cells in a crystal. Divided by , ( ) is normalized so that the integral over frequency becomes 3 , where is the number of atoms.  Considering the atom specific phonon density of states projected along a unit direction vectorˆ, is defined as [37] ( From the canonical distribution in statistical mechanics for phonons under the harmonic approximation, the energy of the phonon system is given as where , B and ℏ are the temperature, the Boltzmann constant, and the reduced Planck constant, respectively. Here, from statistical mechanics,ˆ= 1/[exp(ℏ / B ) − 1] gives the mean phonon number distribution function. A supercell of 1 × 1 × 2 (in , , and -direction) containing 40 atoms was created to study the phonon dispersion relation for CaZrS 3 , using a PHONOPY package with Quantum ESPRESSO package as implemented in [16]. It is known that a crystal constituent of 40 atoms in bulk system (three dimensions , , coordinates) has 120 degrees of freedom (3 , where is the number of atoms). The phonon dispersion relation for frequency bands and frequency density of states were calculated and displayed as shown in figure 5 and figure 6, respectively. As indicated in figure 5, it was observed that there are three (3) acoustic branches and 117 optical branches (3 − 3) mode of vibrations. Here, also the results showed that CaZrS 3 possesses no imaginary phonon frequency modes. Hence, it is structurally and lattice dynamically stable. This finding agrees with the results of the analysis of the elastic constants and Goldschmidt tolerance factor.

Thermodynamic properties
Thermodynamic properties of materials are one of the foundations of solid-state science and industry. The investigation of these properties is important in order to determine their specific behavior when these materials are subjected to high pressure and temperature. Using the thermodynamic relations, a number of thermal properties, such as constant volume heat capacity , Helmholtz free energy , and entropy , can be computed as functions of temperature as [26,36,38] = ∑︁ Here also, a supercell of of 1 × 1 × 2 (in , , and -direction) containing 40 atoms was used to study the thermodynamic properties of CaZrS 3 , using a PHONOPY package with Quantum ESPRESSO package as implemented in [16]. At finite temperatures ranging from 0 K to 1000 K, thermodynamic parameters such as enthalpy, free energy, entropy, and heat capacity were computed and plotted in figure 7. In this consideration, the volume and temperature are independent variables. From figure 7, one can observe that below 10 K, the values of the entropy and heat capacity are almost zero. The free energy diminishes gradually as the temperature rises, whereas the entropy rises rapidly, following reasonable trends shown in [36,39,40]. As a result, the enthalpy increases linearly with the increment of temperature. The increase in enthalpy at high temperatures leads to a decrease in the free energy which is associated with defects. It is clearly observed that for the temperature below 400 K, the heat capacity increases rapidly, whereas for temperature above 400 K it increases slowly (almost increasing linearly) with temperature and gradually approaches the Dulong-Petit limit (classical limit) owing to the anharmonic approximations of the Debye model as observed in [36,41,42]. The calculated heat capacity graph is also smooth and continuous confirming that there is no phase change occurring in CaZrS 3 up to 1000 K [36].

Conclusion
Chalcogenide perovskites have emerged as a non-toxic and stable photovoltaic material that has similar optoelectronic capabilities to lead halide hybrid perovskites. In this study, the first-principles calculations using the Quantum ESPRESSO software package were used to study the structural, electrical, elastic, phonon dispersion relation, and temperature-dependent thermodynamic characteristics of orthorhombic CaZrS 3 for optoelectronic applications. The THERMO_PW package was used to compute the elastic characteristics of the material using the optimized structure. The PHONOPY package was also used to calculate the phonon dispersion and thermal characteristics. The computed lattice parameters = 6.57, = 7.06, and = 9.63 correspond well to theoretical and experimental results. The elastic constants were used to calculate mechanical parameters of CaZrS 3 , such as the bulk modulus, shear modulus, Young's modulus, and elastic anisotropy. CaZrS 3 was found to be classified as a ductile material according to the determined value of / ratio (1.78). Poisson's ratio obtained the value of 0.26, indicating that interatomic forces are central. The measured global anisotropic index value of 0.35 validates the anisotropic nature of CaZrS 3 . The band gap value of CaZrS 3 was calculated by approximating the exchange correlation potential with GGA/PBE and GGA+U. The band gap value calculated using GGA/PBE is 1.23 eV, which is 35% percent less than the experimental result. However, the computed band gap value using GGA+U is 1.88 eV, which is close to the experimental band gap value. The absence of imaginary (negative frequencies in the figures) frequencies in phonon dispersion curve and the phonon density of states give an indication that the structure is dynamically stable. The temperature dependence of thermodynamic parameters including enthalpy, entropy, free energy, and heat capacity is calculated and analyzed.