VIBRATIONAL, OPTICAL AND THERMO-DYNAMICAL BEHAVIOUR OF ZB STRUCTURED BeZ (Z=S, Se AND Te) USING ATK-DFT

The present research is a systematic computational study focused on structural, mechanical, electronic, vibrational, optical and thermo-dynamical properties of zincblende (B3) structured beryllium chalcogenides BeZ (Z=S, Se, Te) compounds using ATK-DFT method using PZ and PBEsol exchange and correlation potentials within the local density approximation (LDA) and the generalized gradient approximation (GGA) respectively and their comparison. The k-point and energy cut-off values were tested and provided convergence in self-consistent calculations. The structural parameters such as lattice constant, bulk modulus, second order elastic constants (C11, C12, C44) and material properties (B, G, Y and σ) for these crystals are computed and discussed. To explain the electronic properties, electronic energy band structure, complex band structures, phonon band structure, phonon density of state and electron density distribution are plotted. The effect of pressure on elastic constant, material properties and phase transitions are also studied, including phase transition from ZB structure to NiAs appearing at 53 GPa, 49 GPa and 33 GPa for BeS, BeSe, and BeTe respectively.


Introduction
The II-VI Semiconductor group material beryllium chalcogenides BeZ (Z=S, Se and Te) are promising materials and have attracted the attention of researchers due to large band gaps (2.7 eV-5.5 eV) [1] and other technological applications in electronic devices; they have remarkable physical, chemical and mechanical properties [2][3][4][5][6][7][8][9][10][11]. Extensive work on these chalcogenides has been done based on optical band gap characterization and phonon properties under various physical conditions of pressure, temperature etc.
These chalcogenides of beryllium, BeZ (Z= S, Se and Te) seem to be exceptionally favourable materials for laser diodes [12] and have lower ionic radius ratio in comparison with other wide-gap chalcogenides, except BeO [2]. At low pressure, beryllium chalcogenides BeZ (Z= S, Se and Te) remain in the fourfold binary zinc-blende (B3) structure. In contrast, available results reported by experimental studies on beryllium chalcogenides show that BeS changes its phase to B8 from B3 at 58 GPa [6]. At the same time, BeSe and BeTe have the transition of the first order from B3 (ZB) to B8 (NiAs) structure at 56±5 GPa and 35±5 GPa [2] of pressure, respectively. These chalcogenides can be easily grown by molecular-beam epitaxy method on several substrates. Owning a higher value of the bulk modulus, beryllium chalcogenides have high values of hardness and stability [13].
Theoretical studies of beryllium chalcogenides have been preferred as compared to the experimental research by most of the researchers due to the highly toxic nature of these materials. In the recent past, numerous experimental and theoretical researchers have reported on beryllium chalcogenides about structural, elastic and electronic properties and illustrated variations in several properties such as electronic, mechanical and phase transition properties . A self-consistent OPW method [25] and the APW method [26] based on first principles were used for the first energy-band calculations of BeS and other beryllium chalcogenides compounds. Further, other theoretical workers have also reported results on structural [27][28], electronic [26][27][28][29], optical, elastic, thermodynamic and lattice dynamical  properties for beryllium chalcogenides with a variation of different physical and chemical properties.
In conclusion, from the preceding, it is clear that numerous studies of beryllium chalcogenides have been done for lattice dynamical, electronic and mechanical properties at different pressures during the past decade, but still, only a few results on the optical constants, elastic constants, electron density, phonon properties are reported till now. The phonon properties of solids play a vigorous role in the extensive studies of several physical properties such as stability of lattice, thermal properties and phase transition have not been reported in a systematized way for beryllium chalcogenides.
In this paper, the author systematically computes and examines the structural, mechanical, electronic, optical, vibrational, and thermo-dynamical properties of B3 structured BeZ materials by very precise ATK-DFT formulations at ground state with LDA and GGA calculation method with the aim of comprehensive results over the experimental and theoretical works on beryllium chalcogenides compounds. A comparative analysis of the presently computed results of this research work with previously reported theoretical and experimental results and found to be in good agreement and favourable in terms of most of the properties which are calculated with ATK-VNL software package. Additionally, the effect of pressure on mechanical properties is also computed and discussed to predict the phase transition behaviour of these chalcogenides. The complete paper is organized into four sections, in section two of this paper, the description of the structural arrangements of the compounds taken for the study and computational methods adopted. The third section is dedicated to the discussion and comparison of the computed and previously reported outcomes, and conclusions based on our results are discussed in section 4.

Computational methods
The zinc blende crystal structured beryllium chalcogenides BeZ (Z= S, Se, Te) of cubic space group 216 (F4-3m) fourfold with four metallic Be atoms located at (0,0,0) and four nonmetallic Z (Z= S, Se, Te) atoms are located at (0.25, 0.25, 0.25) as shown in Figure 1  The entire calculations for beryllium chalcogenides were performed using ATK-VNL 2016.4 version software package based on python programming language [31][32][33]. All the three crystals of beryllium chalcogenides (BeS, BeSe and BeTe) were calculated using Atomistix Tool Kit (ATK). A highly precise first principles DFT calculations of two calculators Local Density Approximation (LDA) and generalized-gradient approximation (GGA). DFT formulation [34] with a numerical linear combination of atomic orbitals (LCAO) basis sets at room temperature (300K) using density mesh cutoff of 75 Hartree and sampling of 5×5×5 k points were performed on a Monkhorst-Pack [35] grid with unpolarised spin. These k-point and energy cut-off values were tested and provided convergence in self-consistent calculations. The exchange and correlation potentials were treated within the Local Density Approximation (LDA) of Perdew and Zunger (PZ) [36]; and generalized-gradient approximation (GGA) of Perdew Burke Ernzerhof solids (PBEsol) for the interaction of the electrons with the ion cores.
The band structure of binary crystals BeS, BeSe and BeTe were calculated for Brillouin zone between the symmetry points G, X, W, L, G, K, X, U, W, K, L by both the set of calculators LDA and GGA. The dynamical matrix was solved for 3×3×3 repetitions for all possible symmetry, and a two dimensional (2-D) visualization of the complex band structure is plotted for BeS, BeSe and BeTe based on the method described by Chang and Schulman [37].
The colour-coding of the k-values operates both to the real and the complex portions of band structure and makes it simpler to understand the closeness of complex bands with the real ones. Further, acoustic sum rule was applied for the calculations of phonon band structure for same Brillouin zone between the symmetry points and phonon density of state (PDOS). The electron density of beryllium chalcogenides is calculated by the occupied eigenstates of the Kohn-Sham Hamiltonian.
The elastic constants have been calculated using ATK-DFT by the Universal Linearly Independent Coupling Strain (ULICS) vectors [38] and the dependent material properties, e.g., bulk modulus, shear modulus, Young's modulus and Poison ratio, are also calculated.

Structural and mechanical properties
Based on computational methodology described in section 2, at ambient conditions the calculated equilibrium lattice constant (a), three independent second order elastic constant (C11, C12 and C44) and bulk modulus (B0) for beryllium chalcogenides (BeS, BeSe, BeTe) are tabulated in Table 1. The results of the present study are compared with available experimental [23,28], and theoretical [5-7, 9, 11, 14, 27-28, 39-40] findings published so far for these crystals. The equilibrium lattice constant computed based on LDA and GGA are exactly similar for beryllium chalcogenides (BeS, BeSe and BeTe) in the present study (ATK-DFT) because of comparable parameters of a calculator and are in best agreement with available experimental values [23]. The deviations in the values of lattice constants with experimental results for BeS, BeSe and BeTe are 0.01Å (-0.20%), 0.07Å (-0.13%) and 0.08Å (1.42%), respectively. Further, the calculated lattice constants are increasing in BeS→BeSe→BeTe order; this trend may be due to the increase of atomic radii in the same order S (1.84Å)→Se(1.98Å)→Te(2.21Å) [41].
The knowledge of elastic constant may enable one to develop a relationship among mechanical and dynamical properties of a crystal as well as may forecasts the description of the forces acting, stability and stiffness of a solid material. Henceforth, the elastic constant for zinc blende semiconducting beryllium chalcogenides are calculated at zero fermi level keeping homogeneous strain. The calculated values of elastic constant presented in Table 1 satisfy the stability criteria in ZB crystals; C11-C12 > 0, C11 > 0, C44 > 0, C11 + 2C12 > 0 and C12 < B < C11 [28]. It is observed that elastic constants C11, C12 and C44 calculated by LDA and GGA have nearly similar values and the value of C11-C12 is reasonably high for zinc blende structure. The value of C44 is increasing in the order BeSe→ BeS→ BeTe and while the calculated B0= (C11 + 2C12)/3 is decreasing from BeS→ BeSe→ BeTe which confirms the validity of Cohen relation [42]. The calculated results for B0 are in good agreement with the experimental results of Madelung [23] and better than other theoretical studies [7,9,11,14,28].
In the present study, LDA computed elastic constants C11, C12 and C44 are in better agreement with available theoretical results [5,9,11,14,28] as compared to GGA computed elastic constants. Due to the absence of experimental result for elastic constants C11, C12 and C44 of these materials, our results have not been compared so far. It is trusted that due to the toxic nature of these material experimental studies have not been done so far.
The next major focus of the present study is on the material properties of the beryllium chalcogenides (BeS, BeSe, BeTe) and is proposed to be unique. Voigt-Reuss-Hill approach [43] has been applied to evaluate bulk modulus (B), Shear modulus (G), Young's Modulus (Y) along each axis and poison ratio (σ) for each face. The calculated material properties of these chalcogens are listed in Table 2 and compared with other experimental and theoretical results [5,9,11,40]

Electronic, vibrational and optical properties
The calculated electronic energy band structures of BeS, BeSe and BeTe for Brillouin zone between the symmetry points G, X, W, L, G, K, X, U, W, K, L by both the set of calculators LDA and GGA, are plotted in Figure 2. From Figure 2, it is reasonably apparent that beryllium chalcogenide crystals have indirect band gaps (Γ-X) of 3.069 eV (for BeS), 2.937 eV (for BeSe) and 1.919 eV (for BeTe), respectively calculated by LDA calculator, indicating that these three compounds are semiconducting. Further, LDA is preferable over GGA as the computed energy band gaps in Table 3 are compared with available experimental results [44][45][46] and other theoretical results [14,30,40,[47][48]. Our computed values of (Γ-Γ=4.941 eV) band gap of BeSe, (Γ-Γ=3.505 eV), (Γ-X=1.919 eV) band gap of BeTe are fairly in good agreement with experiment results [44][45][46], however the small deviation may be due to van der Waals interactions, which was not taken into account during calculations. This is noteworthy to mention that moving from S to Te (S→Se→Te), the changing of the energy of chalcogens p bands go up due to rise of the lattice parameters which is the normal behaviour found for other II-VI zinc blende semiconducting compounds.   Further, complex band structures provide a comprehensive characterization of the static and dynamical electronic behaviour of semiconducting compounds. It is well explained that the complex band structure offers useful information about conventional band structure and its transport behaviour by studying the imaginary components of wave vectors [37]. A two dimensional (2D) visualization of the complex bandstructure is also plotted for BeS, BeSe and BeTe in Figure 3 and the colour-coding of the k-values (0.0 to 1.0) applies both to the real and the complex parts of the band structure, which makes it simpler to see whether the complex bands are closed to the real ones. The imaginary sections of wave vectors depict exponentially decreasing tendency of wave functions and allowed us to forecast both static and dynamic behaviour of compounds, e.g., electron tunnelling. It is evident from the Figure 3 that BeTe is having more complex band structure in real as well in imaginary parts which suggest the preference of BeTe for a broad range of applications: surfaces and interfaces, electron transportation, magnetic tunnel junctions and laser diodes over BeS and BeSe.

Fig. 3. 2D visualization of the complex bandstructure of (a). BeS; (b). BeSe; (c). BeTe. The colour-coding of the k-values applies both to the real and the complex parts of the bandstructure, easier to see where the complex bands are attached to the real ones.
To characterize the phonon modes and dynamical stability behavior of beryllium chalcogenides, acoustic sum rule was applied for calculation of phonon band structure along major symmetry directions (G, X, W, L, G, K, X, U, W, K, L) of the Brillouin zone (BZ) both the set of results from calculators LDA and GGA are shown in Figure 4. The calculated phonon frequencies at the high-symmetry points Γ, X and L for BeZ (Z= S, Se and Te) compared with previously reported results are listed in Table 4. From Figure 4 and Table 4, it is evident that both the calculators LDA and GGA have similar values of phonon frequencies for BeS and closer values for BeSe and BeTe. Our results are in good agreement with experimental results [49][50] and preferably excellent as compared to other theoretical results [5,7,10,[50][51][52] for the major symmetry points.  From Figure 4 and Table 4, it is predictable that the dispersion of transverse (TO) and longitudinal optical phonon (LO) increases upward in the sequence of BeS→BeSe→BeTe along Γ-X with respect to acoustic and optical phonons and the similar characteristics were confirmed by Srivastava et al. [7]. In the phonon band structure, moving toward Γ to K direction, a unique trend of the longitudinal optical branch (LO) indicates a decrease in a band and increase in dispersion in the sequence of BeS→BeSe→BeTe respectively. Such patterns were never reported for other zinc-blends semiconductors.
To understand the degeneracy and splitting of optical phonon modes, phonon density of state (PDOS) corresponding to phonon dispersion curves explained above are plotted in Figure 5 by both the calculator LDA and GGA. It is clear from Figure 5 that for BeS LDA and GGA peaks are closer to each other while for BeSe and BeTe, some variation in the peaks are observed. It is also obvious from Figure 5  The phonon density of state shows complete closeness to all characteristics of phonon dispersion curves. The phonon dispersion curve of beryllium chalcogenides completely confirms that despite the lower weight of S atoms, it has larger frequencies than Se and Te atoms.

Fig. 5. Phonon Density of State (PDOS) of BeZ (Z=S, Se, Te) by LDA (solid green lines) and GGA (Solid blue lines) exchange co-relations of ATK-DFT.
To recognize the importance of bonding in BeZ Compounds, the electron density distribution calculated by GGA calculator is shown in Figure 6. The electron densities near Be atom is higher than nonmetals Z (Z = S, Se, Te) atom, which indicates guiding drift of charges. The colour patterns indicate the increase and decrease of electron density relatively transfer of electrons between Be and Y atoms. From Figure 6, it is evident that the electron density of BeS ranges from 0.029 to 1.3 Å -3 , BeSe from 0.025 to 0.97 Å -3 and BeTe from 0.025 to 0.57 Å -3 . Since the electronegativity of Te is lower as compared to other chalcogens (S=2.58 > Se=2.55 > Te=2.1) which signifies that the BeTe has stronger bonding capabilities than BeSe and BeS owing to mixed covalent-ionic bonding. The increase in the delocalized region of electrons in the interstitial region, confirms the partial metallic character of beryllium chalcogenides. The calculated values of total fractional polarization show a systematic variation and increase moving from BeTe (1.04 C/m 2 ), BeSe (1.25 C/m 2 ), to BeS (1.36 C/m 2 ). In contrast, effective masses show a different characteristics BeS (1.65) → BeTe(1.57) → BeSe (1.55), maybe due to stronger bonding of BeTe. Additionally, the real and imaginary parts of optical spectrum calculated by GGA for beryllium chalcogenides are shown in Figure 7 and observed that the static dielectric constant (ε͚ ) increases in the order BeS (6.53) → BeSe (7.122) → BeTe (7.515). Also, the refractive indices (not reported here) have the same nature increases in the order BeS→BeSe→BeTe.

Thermo-dynamical properties
To understand the complete effect of phonons over the thermodynamical behaviour, we have studied the impact of phonons over the thermodynamical behaviour of beryllium chalcogenides. Further, understanding of the phonon density of states is quite essential for the determination of thermodynamical functions, which is possible only when the proper estimation of PD curve in all B-zone is available. The Debye temperature variation graphs are plotted in Figure 8 for the BeS, BeSe and BeTe, respectively in temperature ranges of 0-400 K. The computed value of Debye temperature at 0 K is 720 K, 451 K and 332 K for BeS, BeSe and BeTe, respectively, which seems to be in good agreement with earlier reported theoretical publications [5,11]. Since, since no experimental results are reported so far on Debye temperature variation, our results may be a reliable prediction for experimental researchers in future. The Debye temperature variation of beryllium chalcogenides in the present study is quite similar to other semiconducting materials at all temperature ranges. From Figure 8, the value of Poisson's ratio may be correlated in binary beryllium chalcogenides as BeS→ BeSe→ BeTe which shows the increasing value Poisson's ratio and quite similar to previously reported values in mechanical properties.

Effect of pressure on mechanical Properties
To have comprehensive information about the pressure dependency on mechanical properties for BeZ, we have computed the effect of pressure ranging from 0-50 GPa on the second order elastic constant (SOEC) and material properties (B, G, Y) using LDA. To supplement the study further, elastic constant variation and material properties (B, G, Y) with respect to the pressure difference for beryllium chalcogenides are presented in Figures 9 and 10, respectively. It is prominently observed that the values of second order elastic constants (C11, C12, C44) has a linear increase with the increase in the pressure. The elastic constant C44 is quite less dependent on the pressure as compared to C11 maybe since transverse strain originates a transformation in shape at constant volume [33]. The other mechanical properties such as bulk modulus (B), shear modulus (G) and Young's modulus (Y) also increase linearly with increase in pressure. Further, at the pressures above 53GPa, 49GPa and 33GPa in case of BeS, BeSe and BeTe, respectively, the elastic constants are not obeying the stability criteria i.e. 0 < C11-C12, 0 < C11, 0 < C44, 0 < C11 + 2C12 and C11 > B > C12. It may be because of the volatility of the zinc blende (B3) phase.

Conclusions
In the present study, a comparative computational analysis of structural, mechanical, electronic, vibrational and thermodynamical properties of zinc-blende (B3) structured BeZ (Z=S, Se, Te) compounds is performed using ATK-DFT method within LDA and GGA at ambient pressure. The calculated equilibrium lattice constants, elastic constants and bulk modulus of beryllium chalcogenides compounds are in reasonably good agreement with the experimental and previously reported theoretical results. The available experimental and other theoretical results on indirect energy band gaps are compared with our computed results and found to be in good agreement. The indirect band gap values calculated based on LDA are found to be better than those based on GGA. The electronic band structure, complex band structure, phonon dispersion band structure, density of states and electron density distribution shows the semiconducting nature at ambient condition for beryllium chalcogenides. In the present study, three optical and three acoustic modes are spotted for BeS, BeSe and BeTe respectively and the existence of these modes may be justified because of binary atom arrangement in the primitive cell. A remarkable property at ambient condition, i.e., dynamical stability of beryllium chalcogenides in ZB (B3) phase, is spotted during the study.
The phonon density of states (PDOS) shows clear supremacy of Be atom for beryllium chalcogenides materials in the region of high vibrational frequency. The LO-TO splitting on the surfaces of BeZ (Z=S, Se, Te) is in good agreement with previously reported results [5,7]. Precisely, these chalcogens own a significant hybridization of s and p orbitals. The s orbital of Be atom revealing blend of covalent and ionic bonding. Further, electron density distribution reveals that BeS and BeSe are covalent while BeTe is covalent-ionic in chemical bonding nature. The polarization and effective mass is also reported and discussed with valid justification. The real and imaginary portions of the dielectric function are computed, and additionally, static dielectric constants ε͚ and refractive index for these chalcogen materials have also been computed. The thermodynamical computations provide a good description of the Debye Temperature of these materials. The Debye temperature variations have established a reasonably good agreement with the formerly published theoretical results. However, no experimental result is available for comparison of Debye Temperature variations of BeS, BeSe and BeTe. There is a reduction in ΘD value when moving forward as BeS→BeSe→BeTe with the rise in temperature. The pressure dependency on mechanical properties such as second-order elastic constants and material properties (B, G, Y) have been computed for the range 0-50 GPa, using LDA. The computed phase transition pressures from zinc blende to NiAs structure are found to be 53 GPa (BeS), 49 GPa (BeSe) and 33 GPa (BeTe).