First-principles study of phase transition, elastic and thermodynamic properties of ZnSe at high pressure

The structural and elastic properties of ZnSe with B3 and B1 phases under different pressure have been investigated by the first principle method based on density functional theory. The obtained structural parameters of ZnSe in both B3 and B1 structures are in good agreement with the available values. The transition pressure of ZnSe from B3 to B1 was predicted as 14.85 GPa by using the enthalpy–pressure data, which is well in line with experimental result. According to the obtained elastic constants, the elastic properties such as bulk modulus, shear modulus, Young’s modulus, ductile/brittle behavior and elastic anisotropy as a function of pressure for polycrystalline of ZnSe are discussed in details. In the frame work of quasi-harmonic Debye model, the temperature and pressure dependencies of the Debye temperature and heat capacity of ZnSe are obtained and discussed in the wide ranges.

The II-VI semiconductor compounds have attracted extensively attention due to their remarkable physical properties [1][2][3][4] and their numerous potential applications 5,6 . Edwards et al. 7 firstly studied the pressure-induced structural phase transformations of II-VI group compounds. Since then, theoretical and experimental studies on the structural and physical behavior of II-VI materials under high pressure are investigated by several research groups.
ZnSe, a significant member of the II-VI semiconductor compounds, which is a light yellow solid compound with a band gap of 2.70 eV 8 . It is an excellent candidate for fabrication of visual displays and photodetectors, etc [9][10][11][12] . It is an essential prerequisite for material's synthesis and application to investigate the fundamental physical properties of ZnSe. Therefore, a great deal of researches have been done on the optical 13,14 , structural 2,15 , electronic 16,17 , and thermodynamic [18][19][20] properties of ZnSe at ambient pressure. According to the result of Edwards et al. 7 , there may be a phase transition for ZnSe with applied pressure. So several research groups studied on the high-pressure behavior of ZnSe and confirmed that it has some polymorphic structures with the pressure increase. Karzel et al. 21 reported the polycrystalline ZnSe from the B3 phase to the B1 phase happens at 13.0 GPa which encouraged great interest for theorists studying of structural stability under elevated pressure. Smelyansky et al. 22 reported the phase transition pressure data as 15 GPa from full-potential linear augmented plane wave (FP-LAPW). While Biering et al. 23 calculated it as 12.94 GPa using the projector augmented wave (PAW) method. When the phase transition happens with the increasing pressure, the nature of crystal structure would be different. However, it is difficult to obtain the exact value of these properties under high pressure in experimental studies. But the fundamental physical properties at elevated pressure, are extraordinary significances for the condensed matter physics, which will contribute to the understanding of the nature of materials. Therefore, the theoretical study could be a powerful tool to acquaint the ZnSe under elevated pressure, owing to the advance in theoretical methods. However, to the best of my knowledge, there are only a few references investigating the elastic and thermodynamic properties of the B3-type ZnSe. Especially, the behaviors of elastic and thermodynamic properties of the B1-type ZnSe are rarely considered under high pressure. In this work, we have focused on the structural phase transition and the elastic properties as a function of pressure for both B3 and B1 phases by using plane-wave pseudopotential density functional theory (DFT). Meanwhile, some detailed thermodynamic property at elevated pressure and temperature have been calculated through the quasi-harmonic Debye model 24 .

open theoretical Methods
The calculations in present work are performed using the pseudopotentials plane-wave approach in the frame of density functional theory (DFT) as implemented in Cambridge Serial Total Energy Package (CASTEP) code 25,26 . For structural calculations, the Perdew-wang-1991(PW91) 27 formulation of the generalized gradient approximation (GGA) was chosen as the optimum exchange correlation of electrons. The Broyden-Fletcher-Goldfarb-Shannon (BFGS) algorithm 28 , which provides a very efficient method to achieve the geometry with a minimum energy, was applied to relax the crystal structure to reach the ground state. In order to obtain an optimum geometric, the kinetic energy cutoff is set as 500 eV and k point separation in the Brillouin zone of the reciprocal space is 8 × 8 × 8 for both B3 and B1 phases. Pseudo atomic calculations of ZnSe are performed for Zn 3d 10 4s 2 and Se 4s 2 4p 4 . The space group of B3 structure is F43m, the positions of the atoms Zn and Se are (0, 0, 0) and (0.25, 0.25, 0.25). The space group of B1 structure is Fm3m, the Zn and Se atoms are located at (0, 0, 0) and (0.5, 0.5, 0.5), respectively. To achieve some reliable results, the self-consistent convergence of the total energy is less than 5.0×10 −7 eV/atom.

Results and Discussion
Pressure-induced structural phase transition. For both B3 and B1 structures of ZnSe, a series of different values of primitive cell volume are set to calculate the total energy. The calculated total energies as a function of volume for both structures of ZnSe are displayed in Fig. 1. According to the result shown in Fig. 1, it is clear to see that the ZnSe with B3 structure is a more stable phase. In order to obtain the equilibrium lattice constants a, the bulk modulus B 0 and its pressure derivative B ' 0 , the total energy E vs. volume is fitted to the Birch-Murnaghan equation of states (EOS) 29 . The results are listed in Table 1, which are also compared with some other theoretical and experimental results. The calculated values of lattice parameters are slight overestimated and the bulk modulus are little underestimated corresponding to the experimental data 21 . The overestimation in the lattice parameters and underestimation in the bulk modulus is a common feature with GGA 30,31 . However, the calculated values using GGA for both B3 and B1 phases agree well with the corresponding experimental value 21 and some available theoretical data 2,22,23,32,33 . In order to find out the most stable phase at finite pressure, the free energy of two structures under different pressures should be calculated. According to thermodynamic stability theory, the lower phase free energy corresponds to the more stable phase. It is well known that the free energy defined as G = E + PV-TS. The last term (TS) can be neglected because our density functional calculations are essentially performed at 0 K. So the free energy G reduces to enthalpy relation H (H = E + PV). When the phase transition happens, the enthalpies of two different structures would be the same. The pressure dependence of the two phases ZnSe are illustrated in the inset of Fig. 1. It is obvious from the inset of Fig. 1 that the enthalpy curve of the B3 structure crosses that of the B1 structure around 14.85 GPa, implying that a solid phase transition from B3 to B1 induced by pressure occurs at 14.85 GPa. The obtained phase transition pressure is quite agree with the previous experimental result 21 and some other theoretical calculations 22,23 .  Single-crystal elastic constants and related properties. As is known to all, the elastic constant is very important parameter for the elastic material and can represent the amount of the elasticity of a material. To investigate the elastic constants of ZnSe, the non-volume conserving method is applied. The elastic constants C ijkl can be described as follows 34,35 Where σ ij stands for the applied stress, e kl corresponds to the Eulerian strain tensors, X denotes the coordinate. The fourth-rank tensor C has generally 21 independent components. If the symmetry of the system is account into consideration, the number of independent components reduces. For a cubic structure, there are only three independent elastic constants, namely C 11 , C 12 , and C 44 . The obtained elastic constants for both phases have been given in Table 2 which also contains some experimental and theoretical values for comparison. The discrepancy between our computed results and the experimental values 36 for the elastic constants of B3 structure is acceptable. Meanwhile there is well in line between our obtained results and some previous theoretical results 32 Fig. 2. It is observed that all the elastic constants for both structures in the considered range of pressure increase with increasing pressure showing a monotonic behavior. For both B3 and B1 phases, it is obvious that C 11 vary more under the impact of pressure than others. The elastic constant C 11 is related to the elasticity in length, which is changed with the longitudinal strain, and the C 12 and C 44 are represented the elasticity in shape. Therefore, the pressure has a much more significant influence on elasticity in length. Moreover, the behavior of elastic constants for B3 phase is still in line with the work of Wang et al. 37 by calculation.
Polycrystalline elastic moduli and related properties. For the sake of ensuring the mechanical stability of ZnSe, the elastic constants of ZnSe in the considered range of pressure should satisfy the mechanical stability criteria. The bulk modulus B and shear modulus G are parameter that must be investigated. The reason is that the bulk modulus B is the macroscopic property of a material which reflects the material's resistance to homogeneous compression and the shear modulus G represents the resistance of the material to shear strain. In order to obtain the bulk modulus B and shear modulus G, the Hill model is applied, which takes the arithmetic average of the Voigt 40   www.nature.com/scientificreports www.nature.com/scientificreports/ where B R denotes Reuss bulk modulus, B V represents Voigt bulk modulus, G R is Reuss shear modulus, G V is the Voigt shear modulus, given as: The Young's modulus E can be regarded as an index to measure the difficulty of producing elastic deformation. It can be given by Figure 3 depicts the pressure dependence of the bulk modulus B, shear modulus G, and Young modulus E of ZnSe in both structures. It is obvious from Fig. 3 that all these elastic moduli (B, G and E) for both B3 and B1 phases increased monotonically with the increasing pressure in the considered range of pressure. The effects of the pressure on B and E are larger than that on G, indicating that the pressure can significantly improve the anti-compression ability and stiffness of ZnSe. The above results show that the increasing pressure can enhance the elastic properties of ZnSe.
In order to distinguish the brittle (ductile) behavior of materials, an important empirical relationship, which is the value of B/G, has been proposed by Pugh 43 . The critical threshold value for analyzing the brittle-ductile behavior of materials is approximately 1.75. As the ratio is up to 1.75, the material shows as a ductile manner. Otherwise, it would possess a brittle character. Figure 4 describes the behavior of B/G ratio for both B3 and B1 structures in the investigated pressure range. The obtained B/G ratio for polycrystalline ZnSe are all larger than 1.75, implying that the polycrystal tend to perform in a ductile manner. Figure 4 further illustrates that the increase of hydrostatic pressure tends to enhance the ductile behavior of ZnSe. The effect of pressure on the ductile of B3 ZnSe is larger than that on the ductile of B1 ZnSe.
Elastic anisotropy. The elasticity anisotropy of the material can reflect the difference properties between in two directions perpendicular. So it is necessary for crystal physics and engineering science to investigate the elasticity anisotropy of the material. The elastic anisotropy can be described the universal elastic anisotropy index A U which is developed by Ostoja-Starzewski 44 for crystal with any symmetry. The A U can be written as follows: where B R denotes Reuss bulk modulus, B V represents Voigt bulk modulus, G R corresponds to Reuss shear modulus, G V is the Voigt shear modulus, respectively. For the case of isotropic crystals, the universal elastic anisotropy index is equal to zero. Contrarily, any value deviate from zero implies the degree of single crystal anisotropy. Figure 5 shows the variation of the A U obtained from our studies under surveyed pressure range. It is obvious www.nature.com/scientificreports www.nature.com/scientificreports/ that the universal elastic anisotropy index A U for both phases are larger than zero and significantly increases with increase of pressure, indicating that the elastic anisotropy in both structures would rise rapidly in the investigated pressure range. The universal elastic anisotropy index of ZnSe in B1 phase increases in smaller slope, indicating the impact of pressure on B1 phase is smaller than that on B3 phase. The minimum value of A U in B3 structure is still bigger than the maximum value of A U in B1 structure in the pressure range of 0-30 GPa, revealing that ZnSe in B3 phase is more anisotropy.
To investigate the anisotropic characteristics of ZnSe, a curved surface of a three-dimensional (3D) representation of the elastic anisotropy of the two cubic ZnSe single crystal was further attempted, which can be expressed by the directional dependences of reciprocals of Young's modulus. For a cubic crystal system, the directional dependence of the Young's modulus can be expressed as 45 . 1  4  11  2  4  22  3  4  33  1  2  2  2  12  1  2  3  2  13  3  2  3  2  23  1  2  2  2  66  1  2  3  2  55  2  2  3  2  44 where S ij corresponds to the elastic compliance constants and l 1 , l 2 and l 3 denote the direction cosines. Figure 6 displays the directional Young's modulus of both phases under different hydrostatic pressures. The surface contours of the Young's modulus of both structures become more anisotropic geometry with an increasing pressure, revealing that ZnSe in both B3 and B1 structures tend to become more anisotropic as the hydrostatic pressure increases. The obtained result is in good agreement with the response of the universal elastic anisotropy index A U illustrated in Fig. 5. Figure 6 further depicts that the higher the pressure is, the larger the 3D figures of Young's modulus is. This result agrees well with the trend of E depicted in Fig. 3. As for both structures at 14.85 GPa, the difference in their the 3D directional dependence of Young's modulus is highly apparent. The 3D directional dependence of the Young's modulus of B3 phase show remarkable anisotropic geometry, indicating that B3 ZnSe is more elastically anisotropic than B1 ZnSe. This result is also in agreement with the above results from the pressure-dependent the universal elastic anisotropy index A U shown in Fig. 5.
Thermodynamic properties. The Debye temperature not only reflects the degree of dynamic distortion of the crystal lattice, but also represents the interatomic binding force of the substance. Many physical quantities of the material are related to it, such as elasticity, hardness, melting point and specific heat. In order to investigate the thermodynamic properties of both phases of ZnSe, the Debye temperature versus temperature and hydrostatic  www.nature.com/scientificreports www.nature.com/scientificreports/ pressure of both phases of ZnSe with a periodic boundary condition have been calculated. Meanwhile, the heat capacity as one of the most important thermodynamic properties was also given in this article. To explore the thermodynamic properties of both phases of ZnSe, the quasi-harmonic Debye model 24,46 is carried out. In which the non-equilibrium Gibbs function G* (V; P, T) is defined the following form: vib Where E(V) denotes the total energy per unit cell, V corresponds to the volume of the molecular system, T stands for the temperature of system, P represents the constant pressure condition, A vib is the vibrational term, which can be written as 47,48

1/1/3 s
Where M is the molecular mass per formula unit, the f(σ) is described in detail 48,50 . B s denotes the adiabatic bulk modulus, which is approximated written in the form of 46 Therefore, the non-equilibrium Gibbs function G*(V; P, T) can be minimized with respect to volume V as follows: The thermal EOS V(P, T) can be obtained by solving Eq. (13). Furthermore, the heat capacity C v is given by the following equation 51 .
The Debye temperature θ and heat capacity C v for both phases of ZnSe at different temperatures and pressures are displayed in Figs. 7 and 8, respectively. The pressure dependence of Debye temperature θ of ZnSe at different temperature is shown in Fig. 7(a). The obtained Debye temperature for B3 ZnSe at zero pressure and 300 K, is 244.42 K, which agrees well with the experimental result 19 . The value of Debye temperature for B1 ZnSe at 14.85 GPa and 300 K locates at 376.64 K. To the best of my knowledge, there is no comparable result of the Debye temperature for B1 ZnSe. So the calculated values will provide a basis for future research work. It is obvious from Fig. 7(a) that as the temperature is fixed, the Debye temperature of ZnSe increases almost linearly www.nature.com/scientificreports www.nature.com/scientificreports/ with an increasing pressure. It is also found the Debye temperature of B3 is smaller than that of B1 at 14.85 GPa. Figure 7(b,c) indicate that the trend of Debye temperature at a given pressure would decreases slightly with temperature. It is clear that the Debye temperature is sensitive to the pressure than the temperature.
As displayed in Fig. 8(a), the heat capacity C v at a given temperature would decrease slightly with the increase of pressure. Under high temperature, the heat capacity is almost constant in the pressure range of 0-30 GPa. The decreasing rate of C v is more obvious under low temperature (e.g. 50 K). The heat capacity for B3 phase at zero pressure and 300 K is 48.27 J·mol −1 K −1 , which is reasonable compared with the available values 20. The C v for B1 phase at 14.85 GPa and 300 K is 46.17 J·mol −1 K −1 , which is reported for the first time, and provides a meaning reference for future research work. As displayed in Fig. 8(b,c), the heat capacity increase rapidly at temperature below the Debye temperatures and slowly approaches a linear increase and then the increasing rate of C v is near zero. The acoustic vibrations play a dominant role in the vibrational excitations at temperatures below the Debye temperatures, so C v strictly follows the T 3 -law 24 . The variation of the heat capacity C v with temperature at intermediate temperatures is dominated by the details of vibration of the atoms. Because of the anharmonic approximation of the Debye model, C v of ZnSe for a given pressure increases rapidly from 0 to 300 K. On account of the impact of anharmonic on the heat capacity C v is suppressed at higher temperature, the heat capacity gradually approaches the Dulong-Pettit limit, which is common to all solids at temperatures far above the Debye temperature. It is remarkable that the temperature effect on C v is much greater than the pressure effect on C v of ZnSe.

Summary and conclusion
The structural and elastic properties of ZnSe for both B3 and B1 structures under different pressures are investigated via the first-principles plane-wave pseudopotential method based on density functional theory (DFT). From the usual condition of equal enthalpies, the phase transition of ZnSe from B3 to B1 occurs at the pressure of 14.85 GPa. According to the obtained elastic constants, the pressure dependence of the bulk modulus, shear modulus and Young's modulus of ZnSe are calculated and discussed in detail. The increasing rate of B/G vs. pressure is remarkably, indicating that ZnSe in two phases is a ductile material and the ductility increases with pressure. The obtained anisotropic indexes and the direction dependence of the Young's modulus demonstrate that ZnSe in B3 phase is more anisotropic than B1 phase and the elastic anisotropy of both phases become stronger with an increasing pressure. The thermodynamic properties of ZnSe, such as Debye temperature and heat capacity as a function of the pressure and temperature are successfully investigated by quasi-harmonic Debye modeling.