Atypical Pressure Dependent Structural Phonon and Thermodynamic Characteristics of Zinc-Blende BeO

1. Introduction

Novel alkaline earth II-oxides (XO: X = Be, Mg, Zn and Cd) occur in different crystalline structures ranging from the wurtzite (wz) → zinc blende (zb) → rock salt (rs) phases. These II-O semiconductors are considered valuable in a broad range of applications [1,2,3,4,5,6]. Unlike other XOs, the beryllium oxide (BeO), commonly known as “beryllia” has attracted significant attention in recent years due to its fascinating electrical, mechanical and thermal properties. Among other II-Os, the BeO possesses a relatively low density, wide bandgap (WBG) ${\mathrm{E}}_{\mathrm{g}}~10.6\mathrm{e}\mathrm{V},$ high rigidity, high melting point 2570°C, and strong chemical bond$.$ As a ceramic, it exhibits an electrical insulating feature with exceptional thermal conductivity κ, and high temperature stability [7,8,9,10,11]. The unique basic properties of BeO have offered major contributions to materials science, condensed matter physics, and electrical engineering. Invaluable physical traits of the material have encouraged many physicists and chemists to investigate its electronic, elastic, thermal and phonon characteristics at ambient, and high-pressure environments [1,2,3,4,5,6]. In recent years, several engineers have meticulously designed BeO-based key electronic components in developing advanced micro-/nanoelectronics, suitable for high temperature applications. Major devices that have been acquired include the ultraviolet (UV) detectors, light emitting diodes (LEDs), laser diodes (LDs), gas sensors, and solar cells, etc. [1,2,3,4,5,6,7,8,9,10,11] Some of these modules are being incorporated in satellites/aircraft for space exploration, in nuclear reactors as neutron reflector; as well as in medical equipment for radiation detection and thermal management needs [7,8,9,10,11]. In medicine BeO ceramics have contributed to lasers used in DNA sequencers and tissue analysis. As a UV transparent conducting oxide (TCO), BeO has been employed in flat-panel displays. Many innovative BeO-based devices have also indicated significant potentials in managing excess heat generation during the operations of electronic and optoelectronic systems.

In BeO, the prospects of spontaneous polarization by an external electric field E can alter the presence of intrinsic defect structures [7,8,9,10,11]. Although not enough conclusive evidence exists of spontaneous polarization modifying the impurity configurations by E in BeO, it is theoretically possible from the available research on related materials [7,8,9,10,11]. An interplay between E and intrinsic defects is a complex subject and requires further study. Nevertheless, it has been substantiated that change of defect structure can alter both the bulk and surface characteristics of the material. Such changes are known to radically impact on the performance of BeO-based devices [7,8,9,10,11]. As a gate dielectric and oxygen diffusion barrier, BeO is established an extremely assuring high-k dielectric for designing metal–oxide–semiconductor field effect transistors (MOSFETs) [1,2,3,4,5,6]. It is also succeeded for engineering many electronic components for power generation requisites. As a well-known refractory compound, BeO is used as a heat sink and substrate in many high-power electronic devices to effectively dissipate heat and prevent overheating. In nuclear applications, the material has been employed as a neutron reflector and/or moderator for improving thermal management [12,13,14,15]. BeO has served as a dispersion phase fuel matrix (DPFM) in nuclear reactors due to its excellent mechanical, chemical, properties with high radiation resistance and good compatibility with uranium dioxide UO₂. The DPFM has not only improved the overall performance of reactors but also provided safety to its fuel. In band gap engineering BeO is frequently mixed with other iso-structural XOs to prepare ternary X_xBe_1-xO and quaternary X_xY_yBe_1-x-yO alloys. In low-dimensional heterostructures (LDHs) [viz., multi quantum wells (MQWs) and superlattices (SLs)] the integration of BeO-based epifilms has offered unlimited opportunities to develop electro-optical device structures. In recent years, many photonic devices are meticulously blended to accomplish and improve flexible micro- /nano-electronics [1,2,3,4,5,6].

To prepare ultrathin XO-based films, different growth methods have been frequently employed [16,17,18,19,20,21,22,23,24,25,26,27,28]. The epitaxial techniques commonly considered are the chemical vapor deposition (CVD), metal organic (MOCVD) [16,17], laser ablation [18,19,20], molecular beam epitaxy (MBE), pulsed laser deposition, and plasma-assisted PA-MBE [21,22,23,24,25,26,27,28], etc. Characterization efforts are also made [29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51] for analyzing and monitoring epifilms and their fundamental features using reflection high-energy electron diffraction (RHEED) [29,30], Auger electron spectroscopy (AES) [31], He+-Rutherford backscattering spectrometry (RBS) [32], atomic force microscopy (AFM) [33,34], high-resolution X-ray diffraction [34,35] (HR-XRD), cross sectional transmission electron microscopy (XTEM), photoluminescence (PL) [36,37,38,39,40], absorption, Fourier transform infrared (FTIR) [41,42], Raman scattering spectroscopy (RSS) [43,44,45,46,47,48], spectroscopic ellipsometry (SE), [49,50,51], etc. Certainly, these experimental methods have validated different crystal structures and also helped evaluate epifilm thickness, strain, electrical, and optical traits of intrinsic defects. In BeO material, the effect of high-pressure P using a diamond anvil cell is studied for assessing its structural, elastic, electronic, and dynamical properties. Analyses of such properties have validated its importance in high-temperature applications. Moreover, these efforts have facilitated in designing many optoelectronic devices viz., insulated gate bipolar transistors (IGBTs); TCO films for charge carrier transport; UV-LEDs; optoelectronic gas sensors; heterostructure field effect transistors (HFETs); high electron mobility transistors (HEMTs) etc. [1,2,3,4,5,6,7,8,9,10,11] Integration of such devices in flexible transparent nano-electronics/nano-photonics has helped achieving high-resolution medical imaging tools for physicians to diagnose, treat patients by recuperating their illnesses, and diseases [12,13,14,15].

Despite the successful growth [16,17,18,19,20,21,22,23,24,25,26,27,28] and characterization [29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51] of epifilms, there are still many basic issues of BeO-based device structures that remained unresolved. For instance, the structural, electronic, elastic, and vibrational traits of epitaxially grown epilayers have not yet been thoroughly investigated. The study of lattice dynamics of BeO is important for understanding its role in different device performance. Vibrational properties

are crucial for evaluating electron-phonon interactions, transport coefficients, as well as assessing the reliability and functionality of electronic devices. Among other II-Os, the phonons in BeO are suggested as the major heat carriers [7,8,9,10,11]. In designing device structures, the zb phase is preferred due to its unique optical and electronic properties. Specifically, the zb BeO has offered many advantages of lower carrier scattering and higher doping efficiencies as compared to the wz structure. The zb BeO is known to have higher κ value for causing higher thermal properties as compared to the wz BeO [7,8,9,10,11].

Inelastic neutron scattering (INS) spectroscopy is considered as one of the most powerful experimental techniques for studying the frequency-wavevector relationships of phonons in solids [47]. Careful analysis of energy lost and/or energy gained by the scattered neutrons in INS has been employed in determining the phonon dispersions${,\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$. Except for zb ZnO, and wz ZnO [47], no experimental results exist for the lattice dynamics of zb BeO material. In the absence of INS measurements, RSS and FTIR spectroscopies are commonly employed [41,42,43,44,45,46,47,48] for assessing the phonon frequencies at high critical points in the Brillouin zone (BZ). In crystals with free charge carriers, Raman scattering has been exploited for comprehending the carrier-phonon interactions and phonon-assisted optical transitions [47,48,49,50,51]. In polar materials, while RSS perceives the inelastic scattering of light by phonons, the FTIR spectroscopy measures the absorption of infrared light by the vibrational modes.

In different solids, P-dependent shifts of phonon frequencies are measured experimentally by using RSS and FTIR spectroscopies. In semiconductors, the analyses of such results have helped gain valuable insights of their bonding characteristics (i.e., the hardening and/or softening of inter-molecular/ intra-molecular interactions). The information of P-dependent phonon shifts in II-VIs, III-Ns and IV-IV materials, has been used earlier for assessing the structural arrangements of atoms for identifying their phase transitions ${\mathrm{P}}_{\mathrm{t}}$. Raman scattering spectroscopy is also exploited in understanding the T-dependent mode Grüneisen parameter$\mathsf{\gamma }\left(\mathrm{T}\right).$ In a solid, the information of $\mathsf{\gamma }\left(\mathrm{T}\right)$ at a given T signifies the sensitivity of P-dependent phonons to its volume change. The P-dependent vibrational modes in different materials have also helped researchers to understand their optical and thermodynamical properties including the specific heat ${\mathrm{C}}_v\left(\mathrm{T}\right)$, Debye temperature ${‌}_{\mathrm{D}}\left(\mathrm{T}\right),$ thermal expansion coefficient (TEC) $\mathsf{\alpha }\left(\mathrm{T}\right)$, etc. [6].

From theoretical standpoints, several ab initio calculations of different wz XOs are performed earlier to understand their structural, electronic, and optical properties [52,53,54,55,56,57,58,59,60,61,62,63,64,65,66]. Typically, the first-principles simulations have employed either the density functional theories (DFT) and/or the density functional perturbation (DFP) methods. Calculations are known by using the full potential linear augmented plane wave (FP-LAPW), frozen phonon (FP), molecular dynamics (MD) and ab initio (AI-MD) schemes [52,53,54,55]. Very few studies exist, however, for materials in the zb phase. Some of the reported results on the zb ZnO material for the electronic and vibrational characteristics using first-principles methods are either inconsistent and/or questionable [52,53]. For instance, Chibueze [52] has recently adopted an ab initio approach to study the mechanical, phonon and electronic properties of zb ZnO. In the framework of DFT, the author [52] has employed a generalized gradient approximation (GGA) by considering a revised Perdew-Burke-Ernzerhof GGA method. In zb ZnO, the author has claimed [52] achieving the degenerate phonon energies at the center of the BZ [i.e., near${\overrightarrow{\mathrm{q}}0,\mathsf{\omega }}_{\mathrm{L}\mathrm{O}\left(‌\right)}={\mathsf{\omega }}_{\mathrm{T}\mathrm{O}\left(‌\right)}~379{\mathrm{c}\mathrm{m}}^{-1}].$ This result for the zb ZnO material is in complete disagreement with the INS experimental as well as first principles theoretical results of phonon dispersions ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ reported by Serrano et al. [48,49] and other LDA simulations [54].

Very few simulations are reported on the phonon dispersions of zb BeO [66]. Absolutely no attempts are made to comprehend the P-dependent ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$and thermodynamical characteristics. High P-phonon studies reflect how the energies of atoms, molecules, or other structural units in solids are arranged for causing modifications in their atomic and/or molecular vibrations. Such studies are required in BeO to decide whether the material is stable or prone to instability. Systematic simulations by exploiting the realistic lattice dynamical models are essential. At ambient and high P, meticulous investigations of zb BeO are crucial to characterize its phonon dispersions ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ and thermodynamical properties [viz., thermal expansion coefficient $\mathsf{\alpha }\left(\mathrm{T}\right)$; specific heat ${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right);$and Debye temperature ${\mathsf{\Theta }}_{\mathrm{D}}\left(\mathrm{T}\right)]$. These traits are expected to play vital roles in evaluating its use in high-T applications especially for thermal management needs in nuclear reactors [7,8,9,10,11].

This paper aims to report the results of a methodical study on the novel zb BeO material by exploiting a rigid-ion model (RIM) [67,68,69,70,71,72,73]. Calculations are performed to understand the structural, phonon and thermodynamical properties of zb BeO at ambient and high P. In Section 2.1, different crystal structures of BeO are described. The salient feature of RIM is outlined in Section 2.2 and 2.2.1–2.2.4 with complete details reported elsewhere [67]. The RIM includes both the short-range and long-range Coulomb interactions. Unique optimization procedures [68] are applied to estimate the interatomic force constants (IFCs) at P = 0 GPa, and P = 20 GPa [75,76,77,78,79]. To achieve the optimized set of IFCs for zb BeO we have incorporated accurate values of phonon frequencies at Γ, X and L critical points in the BZ and the P-coefficients ∂${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$/∂P as input, while employing the elastic constants ${\mathrm{c}}_{\mathrm{i}\mathrm{j}}$, their pressure derivatives ∂${\mathrm{c}}_{\mathrm{i}\mathrm{j}}$/∂P (cf. Sections 3, 3.1-3.2), equilibrium lattice-constants $a_0$ and P-dependent $a/a_0$ as constraints. For zb BeO and using the Murnaghan’s equation of state [74] we have simulated P-dependent volume $V/V_0$ and/or lattice constant ratio $a/a_0$ by considering appropriate values of bulk modulus ${\mathrm{B}}_0$ and its pressure derivative ${\mathrm{B}}_0^{\mathrm{\text{'}}}$ (cf. Sections 3, 3.1-3.2). In zb BeO, the RIM results of phonon dispersions ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right),$ and one-phonon g(ω) density of states (DOS) at P = 0 GPa, and P = 20 GPa are systematically employed in the quasi-harmonic approximation (QHA) to study (cf. Section 3, 3.1-3.3) its thermodynamical properties [e.g., Debye temperature ${‌}_{\mathrm{D}}$(T), heat capacity at constant volume ${\mathrm{C}}_{\mathrm{V}}$(T)]. By incorporating the results of Grüneisen dispersions ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$, we have meticulously achieved the T-dependent Grüneisen parameter γ(T), and linear thermal expansion coefficient α(T). Unlike zb ZnO, the analysis of γ(T) and α(T) for the zb BeO material has revealed atypical behavior (cf. Section 3, 3.1-3.3). The possible explanation for such peculiar trends is ascribed to the combined effect of the short bond length and strong electron localization due to small core size of the electronic charge close to Be in BeO. Results of RIM calculations are compared/contrasted against the existing experimental and first-principles data [80,81,82,83], with concluding remarks presented in Section 4.

2. Theoretical Background

2.1. Structural Properties

Unlike other XOs, the earlier P-dependent measurements for the BeO material have exhibited two main crystal phases. These are (i) the hexagonal wz (B4) $\mathrm{P}{6}_3\mathrm{m}\mathrm{c}\left({\mathrm{C}}_{6\mathrm{v}}^4\right),$and (ii) the sodium chloride or rock salt (rs) (B1) $\mathrm{F}\mathrm{m}\overline{3}\mathrm{m}\left({\mathrm{O}}_{\mathrm{h}}^5\right)$structures. For low P and T (< 2500 K), while BeO occurs in the most stable wz structure it can attain, however, the rs phase at higher P and T [1,2,3,4,5,6]. By using an X-ray synchrotron radiation source, Mori et al. [75] have carried out phase transition studies under high P for BeO. The authors have established that the phase transition ${\mathrm{P}}_{\mathrm{t}}$ appears at $~137\mathrm{G}\mathrm{P}\mathrm{a}$ from wz B4 → rs B1 structure$.$ Structural changes did not evolve for P ≤ 55 GPa. No phase transition${\mathrm{P}}_{\mathrm{t}}$is detected from wz B4 → zb (B3) $\mathrm{F}\overline{4}3\mathrm{m}\left({\mathrm{T}}_{\mathrm{d}}^5\right)$phase [75]. For theorists, the studies of phase transition have faced considerable challenges [84,85,86,87,88,89,90,91,92,93,94]. Unlike experimental studies, the calculations in BeO [84,85,86,87,88,89,90,91,92,93,94] predicted the phase transitions ${\mathrm{P}}_{\mathrm{t}}$ from wz (B4) → zb (B3) → rs (B1). However, the projected values of ${\mathrm{P}}_{\mathrm{t}}$ are perceived varying between 22 GPa $<{\mathrm{P}}_{\mathrm{t}}$ $<$147 GPa. In Figure 1, we have displayed the three possible crystal structures of XOs.

Phase transition

Several DFT simulations are available for the BeO material [84,85,86,87,88,89,90,91,92,93,94] to comprehend its structural phase transitions from B₄ → B₃ → B₁. However, the studies have predicted significantly different values of ${\mathrm{P}}_{\mathrm{t}}$. From these calculations, an average value of ${\mathrm{P}}_{\mathrm{t}}$has been established at nearly ~84 GPa. For BeO, a soft nonlocal pseudopotential method is adopted earlier by Van Camp et al. [93] in the local density function (LDF) approximation. The authors [93] have projected the phase transitions ${\mathrm{P}}_{\mathrm{t}}$from B4 → B3 and B3 → B1 at $~$74 and $~$137 GPa, respectively. Using all electronic and full-potential electronic structure method in the framework of LDF, Boettger et al. [94] have reported the corresponding transitions at ${\mathrm{P}}_{\mathrm{t}}~$63–76 GPa and ${\mathrm{P}}_{\mathrm{t}}~$95 GPa. Later, Park et al. [92] suggested the phase transitions in BeO from B4 → B3 at ${\mathrm{P}}_{\mathrm{t}}~$91 GPa, and B3 → B1 at ${\mathrm{P}}_{\mathrm{t}}~$147 GPa by adopting a first-principles soft nonlocal pseudopotential method in the generalized gradient approximation. In Table 1, we have listed the calculated values of the equilibrium ${\mathrm{P}}_{\mathrm{t}}$obtained for the BeO material between B₄ → B₃ → B₁ by several researchers using different ab initio methods [84,85,86,87,88,89,90,91,92,93,94]. Obviously, the high-pressure status of BeO has been and still is extremely ambiguous [75,84,85,86,87,88,89,90,91,92,93,94].

2.2. Computational Methodology of Lattice Dynamics

In a crystal lattice, the simulated phonon dispersions ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$are used to describe its collective atomic vibrations. In semiconductors, the study of lattice dynamics has played a crucial role for understanding their structural, thermal, optical, and electrical properties. For comprehending ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$in zb materials, two types of theoretical methods are frequently employed. These are: (i) the microscopic or first-principles approaches [52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67], and (ii) the macroscopic techniques [68,69,70,71,72,73,74]. Phonon dispersion studies in solids using microscopic theories rely on their fundamental physical principles which allow detailed understanding of atomic interactions. Macroscopic schemes have, however, employed phenomenological models which require systematic evaluation of the general interatomic force constants. The later methods have often used the simplified representations for numerical evaluation of IFCs by linking them to the material’s fundamental traits rather than examining individual atomic interactions.

In essence, the microscopic techniques provide substantial details for understanding the vibrational characteristics of materials, while the macroscopic methods offer, simplified and efficient perspectives via IFCs to study the lattice dynamical features. For zb BeO, we have exploited a macroscopic RIM (cf. Section 2.2.1) and methodically obtained the optimized sets of its eleven IFCs (cf. Section 2.2.3) at ambient P by incorporating accurate values of phonons (see Table 2), elastic constants ${\mathrm{c}}_{\mathrm{i}\mathrm{j}}$, and lattice constant $a_o.$

2.2.1. Rigid-ion Model

The macroscopic RIM that we have adopted here to study the vibrational characteristics of zb BeO material includes both the short-range (up to 2nd nearest neighbors) and long-range Coulomb interactions [67]. In this model, an atom is identified by using two indices: l and κ. The term l represents the number of unit cells while κ signifies the two types of atoms (i.e., κ = 1 (O) and κ = 2 (Be)). In RIM [67], the polarization of zb crystal lattice is determined by point ions displacements from their equilibrium positions anticipating them to be rigid and non-polarizable. For zb BeO, the phonon dispersions ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ are obtained in the harmonic approximation by solving the crystal Hamiltonian [67]:

$\mathrm{H}=\sum _{\mathrm{l}\mathsf{\alpha }}{\displaystyle \frac{{\mathrm{p}}_{\mathsf{\alpha }}^2\left(\mathrm{l}\right)}{2{\mathrm{M}}_{‌}}}+{‌}_0+{\displaystyle \frac{1}{2}}\sum _{\mathrm{l}\mathsf{\alpha },{\mathrm{l}}^{\mathrm{\text{'}}}{‌}^{\mathrm{\text{'}}}\mathsf{\beta }}{‌}_{\mathsf{\alpha }\mathsf{\beta }}(\mathrm{l},{\mathrm{l}}^{\mathrm{\text{'}}}{‌}^{\mathrm{\text{'}}}){\mathrm{u}}_{\mathsf{\alpha }}\left(\mathrm{l}\right){\mathrm{u}}_{\mathsf{\beta }}\left({\mathrm{l}}^{\mathrm{\text{'}}}{‌}^{\mathrm{\text{'}}}\right)$, (1)

where ${\mathrm{u}}_{\mathsf{\alpha }}\left(\mathrm{l}\right)$ represents the displacement of $\mathsf{\alpha }$-component for the

^th (≡1, 2) type atom from equilibrium in $l$^th unit cell. The term ${\mathrm{p}}_{\mathsf{\alpha }}\left(\mathrm{l}\right)$ signifies the corresponding components of its momentum. For zb crystals, one can write the potential energy ${‌}_{\mathsf{\alpha }\mathsf{\beta }}\left(\mathrm{l},{\mathrm{l}}^{\mathrm{\text{'}}}{‌}^{\mathrm{\text{'}}}\right),$ by splitting it into a short-range repulsive, and long-range attractive Coulomb part as [67]:

${‌}_{\mathsf{\alpha }\mathsf{\beta }}\left(\mathrm{l},{\mathrm{l}}^{\mathrm{\text{'}}}{‌}^{\mathrm{\text{'}}}\right)={‌}_{\mathsf{\alpha }\mathsf{\beta }}^{\mathrm{s}}\left(\mathrm{l},{\mathrm{l}}^{\mathrm{\text{'}}}{‌}^{\mathrm{\text{'}}}\right)+{{\mathrm{Z}}_{‌}{\mathrm{Z}}_{{‌}^{\mathrm{\text{'}}}}}_{\mathsf{\alpha }\mathsf{\beta }}^{\mathrm{C}}(\mathrm{l},{\mathrm{l}}^{\mathrm{\text{'}}}{‌}^{\mathrm{\text{'}}})$ , (2)

with ${\mathrm{Z}}_{‌}$e (≡${\mathrm{Z}}_{\mathrm{e}\mathrm{f}\mathrm{f}}$) being the charge on the ${‌}^{\mathrm{t}\mathrm{h}}(\equiv 1,2)$type ions.

The equations of motion in the harmonic approximation take the following form [67]:

${\mathrm{M}}_{‌}{\ddot{\mathrm{u}}}_{\mathsf{\alpha }}\left(\mathrm{l}\right)=-\sum _{{\mathrm{l}}^{\mathrm{\text{'}}}{‌}^{\mathrm{\text{'}}}\mathsf{\beta }}{‌}_{\mathsf{\alpha }\mathsf{\beta }}(\mathrm{l},{\mathrm{l}}^{\mathrm{\text{'}}}{‌}^{\mathrm{\text{'}}}){\mathrm{u}}_{\mathsf{\beta }}\left({\mathrm{l}}^{\mathrm{\text{'}}}{‌}^{\mathrm{\text{'}}}\right)$, (3)

where the atomic displacement is expressed as a plane wave of the type:

${\mathrm{u}}_{\mathsf{\alpha }}\left(l|\overrightarrow{\mathrm{q}}\mathrm{j}\right)={\displaystyle \frac{1}{\sqrt{{\mathrm{M}}_{‌}}}}{\mathrm{e}}_{\mathsf{\alpha }}\left(\right|\overrightarrow{\mathrm{q}}\mathrm{j}){\mathrm{e}}^{\mathrm{i}[\overrightarrow{\mathrm{q}}\overrightarrow{\mathrm{x}}\left(l\right)-{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)\mathrm{t}]}$; with $\mathsf{\alpha }=$ x, y, z , (4)

with $\mathrm{t}$ being the time; $\overrightarrow{x}\left(\mathrm{l}\right)$ and ${\mathrm{M}}_{‌}$ represent the position and mass of the $\left(l\right)$ atom, respectively. Substituting Equation (4) into Equation (3), it is possible to re-write the equations of motion as [67]:

${\mathsf{\omega }}_{\mathrm{j}}^2\left(\overrightarrow{\mathrm{q}}\right){\mathrm{e}}_{\mathsf{\alpha }}\left(\right|\overrightarrow{\mathrm{q}}\mathrm{j})=\sum _{{‌}^{\mathrm{\text{'}}}\mathsf{\beta }}{\mathrm{D}}_{\mathsf{\alpha }\mathsf{\beta }}^{\mathrm{s}\mathrm{C}}(\mathsf{\kappa }\mathsf{\kappa }\mathrm{\text{'}}\left|\overrightarrow{\mathrm{q}}\right){\mathrm{e}}_{\mathsf{\beta }}\left(\mathsf{\kappa }\mathrm{\text{'}}\right|\overrightarrow{\mathrm{q}}\mathrm{j});\mathsf{\kappa },\mathsf{\kappa }\mathrm{\text{'}}=1,2$ (5)

where ${\mathrm{D}}_{\mathsf{\alpha }\mathsf{\beta }}^{\mathrm{s}\mathrm{C}}\left(\mathsf{\kappa }\mathsf{\kappa }\mathrm{\text{'}}\right|\overrightarrow{\mathrm{q}})\left[\equiv {\mathrm{D}}_{\mathsf{\alpha }\mathsf{\beta }}^{\mathrm{s}}\left(\mathsf{\kappa }\mathsf{\kappa }\mathrm{\text{'}}\right|\overrightarrow{\mathrm{q}})+{\mathrm{D}}_{\mathsf{\alpha }\mathsf{\beta }}^{\mathrm{C}}(\mathsf{\kappa }\mathsf{\kappa }\mathrm{\text{'}}\left|\overrightarrow{\mathrm{q}}\right)\right]$represents the dynamical matrix comprising of the short-range ${\mathrm{D}}_{\mathsf{\alpha }\mathsf{\beta }}^{\mathrm{s}}\left(\mathsf{\kappa }\mathsf{\kappa }\mathrm{\text{'}}\right|\overrightarrow{\mathrm{q}}),$and long range Coulomb ${\mathrm{D}}_{\mathsf{\alpha }\mathsf{\beta }}^{\mathrm{s}}\left(\mathsf{\kappa }\mathsf{\kappa }\mathrm{\text{'}}\right|\overrightarrow{\mathrm{q}})$interactions. For each mode ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ the components of its eigenvectors ${\mathrm{e}}_{\mathsf{\alpha }}\left(\right|\overrightarrow{\mathrm{q}}\mathrm{j})$in Equation (5) satisfy the familiar orthogonality [67]:

$\sum _{\mathsf{\alpha }\mathsf{\kappa }}{\mathrm{e}}_{\mathsf{\alpha }}^{\mathrm{*}}\left(\right|\overrightarrow{\mathrm{q}}\mathrm{j}\left){\mathrm{e}}_{\mathsf{\alpha }}\right(\left|{\overrightarrow{\mathrm{q}}\mathrm{j}}^{\mathrm{\text{'}}}\right)={\mathsf{\delta }}_{\mathrm{j}{\mathrm{j}}^{\mathrm{\text{'}}}}$, (6a)

and closure relations

$\sum _{\mathrm{j}}{\mathrm{e}}_{\mathsf{\alpha }}^{\mathrm{*}}\left(\mathsf{\kappa }\mathrm{\text{'}}\right|\overrightarrow{\mathrm{q}}\mathrm{j}\left){\mathrm{e}}_{\mathsf{\beta }}\right(\left|{\overrightarrow{\mathrm{q}}\mathrm{j}}^{\mathrm{\text{'}}}\right)={\mathsf{\delta }}_{{‌}^{\mathrm{\text{'}}}}{\mathsf{\delta }}_{\mathsf{\alpha }\mathsf{\beta }}$. (6b)

2.2.2. The Quasi Harmonic Approximation

In the harmonic approximation the knowledge of IFCs at ambient P for zb BeO is required to obtain the phonon dispersions${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ (cf. Section 3). The inclusion of higher terms in Equation (1) makes the exact solutions to the equations of motion impossible. Thus, approximations must be made. A QHA is employed to calculate thermodynamical properties of materials for predicting their phase diagram at higher P and T [84]. In this approximation, the simulations are made to study how the atomic vibrations (phonons) of a solid change with temperature, T and volume V. In QHA, the interatomic potential is still terminated at the quadratic term, allowing the definition of normal modes as in the harmonic approximation. But now the IFCs are allowed to alter with the change of interatomic distance $a_0$. This means that the strength of IFCs is modified either by the variation of volume V, or pressure P. Thus, one expects the phonon dispersions ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ and density of states DOS of the material (cf. Section 3, 3.1-3.2) to shift with the change of V (or P), be it due to thermal expansion or by an externally applied stress (pressure, P).

2.2.3. Interatomic Force Constants of zb BeO at P = 0 GPa

Systematic procedures [68] are followed to optimize the eleven RIM parameters at ambient conditions P = 0 GPa for zb BeO. Appropriate values of phonon frequencies at critical points (see 2nd column of Table 2) are used as input while the elastic- ${\mathrm{c}}_{\mathrm{i}\mathrm{j}}$ and lattice- ${\mathrm{a}}_0$constants are employed as constraints. To obtain IFCs (see Table 3) at P = 20 GPa, we have considered the values of ${\partial \mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$/∂P as input from high-P RSS measurements [95] or reliable theoretical methods and used simulated pressure-dependent $\partial {\mathrm{c}}_{\mathrm{i}\mathrm{j}}/\partial \mathrm{P}$ elastic- and $\partial {\mathrm{a}}_0/\partial \mathrm{P}$ lattice constants (cf. Section 3) as constraints. Once the optimized set of RIM interatomic force constants $\mathrm{A},\mathrm{B},{\mathrm{C}}_{\mathsf{\kappa }},{\mathrm{D}}_{\mathsf{\kappa }},{\mathrm{E}}_{\mathsf{\kappa }},{\mathrm{F}}_{\mathsf{\kappa }}$, and effective charge ${\mathrm{Z}}_{\mathrm{e}\mathrm{f}\mathrm{f}}(\equiv {\mathrm{Z}}_{‌}$e) are obtained for the zb BeO material at ambient pressure (P= 0) and high pressure (P = 20 GPa), it is straight forward to simulate${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right),$ one phonon DOS by using Equation (5).

2.3. Thermodynamical properties

Application of isotropic P in solids causes the decrease (increase) of atomic distances (bonding interactions). These changes trigger the variations in their phonon frequencies. The T- and/or P-dependent lattice parameter ${\mathrm{a}}_0$ in a solid is commonly measured by HR-XRD [76,77]. In most semiconductors the effect of P is known to produce stiffening of the longitudinal-, transverse-optical (${\mathsf{\omega }}_{\mathrm{L}\mathrm{O}},{\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}$) phonons, longitudinal acoustic (${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}$) modes while instigating softening in their transverse acoustic (${\mathsf{\omega }}_{\mathrm{T}\mathrm{A}}$) phonon branches. Pressure-induced stiffening/softening influences the phase space for phonon scattering and thus affects their thermal conductivity. Besides thermal conductivity κ$\left(\mathrm{T}\right)$, specific heat ${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right)$, thermal expansion coefficient $\mathsf{\alpha }\left(\mathrm{T}\right),$and Grüneisen parameters $\mathsf{\gamma }\left(\mathrm{T}\right)$, are the other three most important characteristics of the zb semiconductors [76,77]. These features have contributed to the fundamental understanding of lattice anharmonicity.

The impact of P on ${\mathrm{C}}_v\left(\mathrm{T}\right)$ is rather complex. The compression often leads to increase in the vibrational frequencies. One, therefore, expects shift of ${\mathrm{C}}_v\left(\mathrm{T}\right)$ peak to higher T. Debye temperature relates to the maximum vibrational frequency of a solid and can decrease at lower T due to softening of ${\mathsf{\omega }}_{\mathrm{T}\mathrm{A}}$ modes [76,77]. Grüneisen parameter describes the anharmonicity of lattice vibrations and is related to the P dependence of phonon modes. Thermal expansion coefficient describes how a material expands or contracts with the change of T and can be affected by pressure, P. Researchers have typically used a differential scanning calorimetry [76,77] to study experimentally the specific heat ${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right)$ of a material. In this method the crystal is heated at a controlled rate by comparing heat flow to the reference sample. This study allows estimation of ${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right)$ based on the temperature change and heat absorbed by the crystal. Debye temperature ${‌}_{\mathrm{D}}\left(\mathrm{T}\right)$ is assessed from the experimental data of specific heat. Thermal expansion coefficient $\mathsf{\alpha }\left(\mathrm{T}\right)$is generally measured in a three-terminal capacitance dilatometer [76]. Except for the wz BeO, no experimental TEC results

are known for the zb BeO material. Thermodynamic characteristics of zb BeO material can help evaluate the functionality of its use in thermal-management applications [1,2,3,4,5,6,7,8,9,10,11]. In the absence of experimental data, one must study these thermodynamical traits (cf. Section 2.3.1) by using a realistic lattice dynamical model.

2.3.1. Thermodynamical properties

At a constant P, the thermal expansion coefficient $\mathsf{\alpha }\left(\mathrm{T}\right)$ of a material reflects the fractional change of its dimension (either length l, or V) per degree change in the temperature, T. It accounts for the V- and/or P-dependent phonon frequencies. In an isotropic solid $\mathsf{\alpha }\left(\mathrm{T}\right)$ can be calculated in the QHA by minimizing the Helmholtz free energy with respect to V, while considering the temperature dependence of lattice constant and phonon mode frequencies [76,77]:

$\mathsf{\alpha }={\displaystyle \frac{1}{3{\mathrm{B}}_0}}\sum _{\overrightarrow{\mathrm{q}},\mathrm{j}}{\mathrm{C}}_{\mathrm{v}}\left(\overrightarrow{\mathrm{q}},\mathrm{j}\right){\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$, (7)

where ${\mathrm{B}}_0$ is the bulk modulus. In Equation (7), the summation on the right-hand side is performed over the phonon modes$\left(\overrightarrow{\mathrm{q}},\mathrm{j}\right)$. At a constant volume, the vibrational frequencies ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ of a material are assumed independent of T. The relative change of phonon frequencies ${\displaystyle \frac{\mathrm{d}{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)}{\mathrm{d}\mathrm{V}}}$ with volume V (or pressure P, ${\displaystyle \frac{\mathrm{d}{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)}{\mathrm{d}\mathrm{P}}}$) is usually described by the mode Grüneisen dispersions ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right).$ This dimensionless quantity is defined as [76,77]:

${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)=-{\displaystyle \frac{\mathrm{d}\left(\mathrm{l}\mathrm{n}{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)\right)}{\mathrm{d}\mathrm{l}\mathrm{n}\mathrm{V}}}=-{\displaystyle \frac{\mathrm{V}}{{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)}}{\displaystyle \frac{\mathrm{d}{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)}{\mathrm{d}\mathrm{V}}}={\displaystyle \frac{{\mathrm{B}}_0}{{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)}}{\displaystyle \frac{\mathrm{d}{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)}{\mathrm{d}\mathrm{P}}}$ . (8)

The above Equation (8) serves to quantify the vibrational anharmonicity. Positive sign of ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$implies that the phonon frequencies are increasing with the decrease in V (or increase in P). This causes a positive coefficient of thermal expansion and allows for dimensional stability of the material when exposed to T variations. This crucial property of the material is essential for applications requiring high thermal conductivity and electrical insulation. The low positive $\mathsf{\alpha }\left(\mathrm{T}\right)$ in a semiconductor helps minimize stress and potential damage caused by T fluctuation. The relevance of a positive $\mathsf{\alpha }\left(\mathrm{T}\right)$ for zb BeO will be discussed in Sections (3.3-3.4).

Using ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right),$ thermal Grüneisen parameter $\mathsf{\gamma }\left(\mathrm{T}\right)$ can be obtained by using [76,77]:

$\mathsf{\gamma }\left(\mathrm{T}\right)={\displaystyle \frac{\sum _{\overrightarrow{\mathrm{q}},\mathrm{j}}{\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right){\mathrm{C}}_{\mathrm{v}}\left(\overrightarrow{\mathrm{q}},\mathrm{j}\right)}{\sum _{\overrightarrow{\mathrm{q}},\mathrm{j}}{\mathrm{C}}_{\mathrm{v}}\left(\overrightarrow{\mathrm{q}},\mathrm{j}\right)}}$ , (9)

where the contribution of each mode ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ to $\mathsf{\gamma }\left(\mathrm{T}\right)$is weighted by its contribution to the specific heat${\mathrm{C}}_{\mathrm{v}}\left(\overrightarrow{\mathrm{q}},\mathrm{j}\right)$. The denominator in Equation (9) is the heat capacity at constant volume which can take the following form [76,77]:

${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right)={\mathrm{k}}_{\mathrm{B}}\sum _{\overrightarrow{\mathrm{q}},\mathrm{j}}{\left({\displaystyle \frac{\mathrm{ħ}{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)}{2{\mathrm{k}}_{\mathrm{B}}\mathrm{T}}}\right)}^2{\displaystyle \frac{1}{{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}}^2\left(\frac{\mathrm{ħ}{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)}{2{\mathrm{k}}_{\mathrm{B}}\mathrm{T}}\right)}}$ . (10)

Where T is the temperature; ${\mathrm{k}}_{\mathrm{B}}$ and $\mathrm{ħ}$ are the Boltzmann and Planck constants, respectively. It is possible to express ${\mathrm{C}}_{\mathrm{v}}$ via the calculated phonon dispersions ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ and DOS $\mathrm{g}\left(\mathsf{\omega }\right)$. Thus, an equivalent form of Equation (10) can be rewritten as [76,77]:

${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right)={\mathrm{N}\mathrm{r}\mathrm{k}}_{\mathrm{B}}{\int }_0^{\mathrm{\infty }}\mathrm{d}\mathsf{\omega }\mathrm{g}\left(\mathsf{\omega }\right)\left({\displaystyle \frac{\mathrm{ħ}}{k_{B}T}}\right){\displaystyle \frac{\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{ħ}{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)}{{\mathrm{k}}_{\mathrm{B}}\mathrm{T}}\right)}{{\left[\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{ħ}{\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)}{{\mathrm{k}}_{\mathrm{B}}\mathrm{T}}\right)-1\right]}^2}}$ (11)

By using Equations (9-11), the T-dependence on $\mathsf{\gamma }\left(\mathrm{T}\right)$ are evaluated [76,77]:

$\mathsf{\gamma }\left(\mathrm{T}\right)={\displaystyle \frac{\sum _{\overrightarrow{\mathrm{q}},\mathrm{j}}{\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right){\mathrm{C}}_{\mathrm{v}}\left(\overrightarrow{\mathrm{q}},\mathrm{j}\right)}{{\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right)}}$ . (12)

From the Debye’s equation, ${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right)$ can also be expressed as [76,77]:

${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right)={9\mathrm{r}\mathrm{k}}_{\mathrm{B}}{\left({\displaystyle \frac{\mathrm{T}}{{‌}_{\mathrm{D}}\left(\mathrm{T}\right)}}\right)}^3{\int }_0^{{‌}_{\mathrm{D}}\left(\mathrm{T}\right)/\mathrm{T}}{\displaystyle \frac{{\mathrm{x}}^4{\mathrm{e}}^{\mathrm{x}}}{{\left({\mathrm{e}}^{\mathrm{x}}-1\right)}^2}}\mathrm{d}\mathrm{x}$ , (13)

In zb BeO, once the complete phonon dispersions ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right),\mathrm{D}\mathrm{O}\mathrm{S},$ Grüneisen dispersions ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$are obtained for the wave vectors throughout the BZ, thermodynamical characteristics including the Debye temperature ${‌}_{\mathrm{D}}\left(\mathrm{T}\right),$ $\mathsf{\gamma }\left(\mathrm{T}\right)$ or $\mathsf{\alpha }\left(\mathrm{T}\right)$can be easily simulated.

3. Numerical Computations, Results and Discussions

For zb BeO, the RIM phonon dispersions ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right),$ and DOS $\mathrm{g}\left(\mathsf{\omega }\right)$are required for assessing its thermodynamic characteristics. By using the structural and elastic properties of the material we have obtained in Section 2.2.3, the unique set of eleven IFCs at ambient P = 0 GPa. The simulated P-dependent lattice and elastic constants (cf. Section 3.1) of zb BeO are necessary for obtaining IFCs at higher P≠ 0 GPa.

3.1. Pressure Dependent Lattice and Elastic Constabts of zb BeO

Murnaghan equation of state (MEOS) is a simple relationship which links the applied P to the lattice constant a (or volume V) of a material [75]. It can be used to understand how much a solid compress by increasing isotropic P. By using appropriate parameters (see Table 2) of zb BeO we have calculated the P-dependent changes in ${\displaystyle \frac{a}{a_0}}\left[\mathrm{o}\mathrm{r}{\displaystyle \frac{\mathrm{V}}{{\mathrm{V}}_0}}\right]$ [75]:

${\displaystyle \frac{a}{a_0}}={\left[{\displaystyle \frac{B_0^{\text{'}}}{B_0}}P+1\right]}^{-{\displaystyle \frac{1}{3B_0^{\text{'}}}}}$or${\displaystyle \frac{V}{V_0}}={\left[{\displaystyle \frac{B_0^{\text{'}}}{B_0}}P+1\right]}^{-1/B_0^{\text{'}}}$ , (14)

where ${\mathrm{B}}_0$is the bulk modulus and ${\mathrm{B}}_0^{\mathrm{\text{'}}}$ is its pressure derivative$.$The results of P-dependent variations of ${\displaystyle \frac{a}{a_0}}$ and ${\displaystyle \frac{\mathrm{V}}{{\mathrm{V}}_0}}$, are displayed in Figure 2 a) and 2b) respectively for 0 < P < 120 GPa. Modifications in P-dependent elastic constants $c_{ij}$ displayed in Figure 2c) revealed satisfying the mechanical stability conditions in zb BeO.

3.1.1. Interatomic Force Constants of zb BeO at P≠ 0 Gpa

For zb BeO we have optimized the set of eleven RIM parameters (see Section 2.2.3) at ambient P. To achieve IFCs at P ≠ 0, the P-dependent ∂ω/∂P phonons at a few critical points are used as input while ∂a/∂P (see Figure 2 (a)) and ∂c_ij/∂P (see Figure 2 (c)) are employed as constraints. Perusal of Table 3 has clearly indicated that IFCs of zb BeO changed significantly with P. These outcomes are consistent with similar results obtained earlier in many compound semiconductors [69,70,71,72,73,74]. To appreciate the importance of two sets of RIM parameters, we have considered a linear interpolation scheme [69,70]:

${\mathrm{a}}_{\mathrm{i}}\left(\mathrm{P}\ne 0\right)={\mathrm{a}}_{\mathrm{i}}\left(\mathrm{P}=0\right)+\mathrm{P}{\displaystyle \frac{\mathrm{d}{\mathrm{a}}_{\mathrm{i}}}{\mathrm{d}\mathrm{P}}}$ , (15)

where${\mathrm{a}}_{\mathrm{I}}$ (I = 1, 11) represents the values of RIM parameters. By using Equation (15) and IFCs (see Table 3) one can calculate P-dependent force constants at any desired pressure. These RIM force constants have helped us assessing the P-dependent phonon dispersions (see Section 3.1.2) for comprehending different thermodynamical (see Section 3.1.3) characteristics of zb BeO.

3.1.2. Phonon Dispersions and Density of States

For zb BeO, the RIM calculations of phonon dispersions${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ at ambient P = 0 (blue color lines) and P = 20 Gpa (red color lines) are displayed in Figure 3 a) along the high symmetry directions (Γ → X → K → Γ → L → X → W → L) of the BZ. Simulated one-phonon DOS $\mathrm{g}\left(\mathsf{\omega }\right),$ at P = 0 and P = 20 Gpa are also displayed in Figure 3 b). In the absence of experimental data, the RIM phonon values at P = 0 Gpa are compared reasonably well with the ab initio results [66]. Like ambient pressure P (= 0), the RIM results at P = 20 Gpa (red color lines) in zb BeO material for ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ and $\mathrm{g}\left(\mathsf{\omega }\right)$ exhibited identical features (Figures 3 a), 3 b)). Unlike zb ZnO and other II-VI semiconductors, in zb BeO our study has noticed atypical phonon features [95,96,97]. The study projected no appreciable P-dependent changes in the ${\mathsf{\omega }}_{\mathrm{T}\mathrm{A}}$modes while all other phonons are shifted to higher frequencies at P = 20 Gpa. Similar outcomes can be examined from the calculated one phonon DOS (red color lines) revealing no noticeable changes in the ${\mathsf{\omega }}_{\mathrm{T}\mathrm{A}}$modes. Other high frequency phonon branches have indicated upward shifts resulting in a relatively wider minimum phonon gap between 850 cm^-1 – 960 cm^-1.

The study of phonon dispersions ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ for the zb BeO has revealed a large optical phonon splitting ~353 cm^-1${(\mathsf{\omega }}_{\mathrm{L}\mathrm{O}\left(‌\right)}-{\mathsf{\omega }}_{\mathrm{T}\mathrm{O}\left(‌\right)})$at the center of the BZ (i.e., at Γ point). Moving away from Γ point, the ${\mathsf{\omega }}_{\mathrm{L}\mathrm{O}}\left({\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}\right)$modes have exhibited high (low) dispersive behavior along the [Γ - X], [X- Γ], and [Γ - L] directions (see Figure 3 a)). The effect of such behavior is clearly reflected in the calculated DOS $\mathrm{g}\left(\mathsf{\omega }\right)$(see Figure 3 b)). Results of $\mathrm{g}\left(\mathsf{\omega }\right)$disclosed a broad band between 825 cm^-1-1075 cm^-1 with a minimum gap appearing in the frequency range of 744 cm^-1 – 825 cm^-1 at P = 0 Gpa (blue colored arrows), which shifts to higher frequency range of 850 cm^-1 – 960 cm^-1at P = 20 Gpa (red colored arrows).

Different from II-VI materials, another unusual behavior is noticed in zb BeO along the Γ-X and Γ-L directions where the ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}$modes (see Figures 3 a)) exhibit higher frequencies than ${\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}$phonons (i.e., ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}>{\mathsf{\omega }}_{\mathrm{T}\mathrm{O}})$. The higher frequency of ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}$ mode than ${\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}$could potentially influence thermal conductivity by impacting phonon scattering mechanisms. Again, unlike zb BeO, most zb II-VI materials including zb ZnO have shown flatness of ${\mathsf{\omega }}_{\mathrm{T}\mathrm{A}}$ branches along the [Γ - X], and [Γ - L] directions. Accordingly, in $\mathrm{g}\left(\mathsf{\omega }\right)$ our study (see Figures 3 b)) perceived no clear phonon gap between the acoustic and optical phonon branches. Interestingly, identical trends are noted in zb BN and diamond C, for their ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$[96,97,98,99,100]. This provides a strong support to our earlier assertion that the structural, elastic and vibrational characteristics of zb BeO exhibit similarities with super hard BN and C materials [71]. Moreover, in zb BeO the eigenvectors are linked to ${\mathsf{\omega }}_{\mathrm{L}\mathrm{O}\left(\mathrm{X}\right)}$vibrations of lighter Be atom and heavier O oscillation to ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}\left(\mathrm{X}\right)}$ phonon [66]. In Table 4, results of RIM phonon frequencies at high critical points are listed along with linear pressure coefficients ${\mathrm{a}}_{\mathrm{j}}^{\mathrm{P}}$(cm^-1/Gpa) and Grüneisen parameters ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$.

The results of ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ and $\mathrm{g}\left(\mathsf{\omega }\right)$ for zb BeO at ambient and higher P are used next to simulate the thermodynamical characteristics. Major impact of T is reported on different quantities including the Debye temperature ${‌}_{\mathrm{D}}\left(\mathrm{T}\right)$, specific heat ${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right),$Grüneisen dispersions ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right),$and linear thermal expansion coefficient $\mathsf{\alpha }\left(\mathrm{T}\right)$ (see Section 3.1 and 3.2).

(cm^-1/Gpa) and Grüneisen parameters ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$are also reported.

3.1.3. Debye Temperature and Specific Heat

In Figure 4 a) we have displayed our calculated RIM results of zb BeO for the Debye temperature${‌}_{\mathrm{D}}\left(\mathrm{T}\right)$at P = 0 (blue colored line) and P = 20 Gpa (red colored lines) in the T range, 0 ≤ T ≤ 1600 K. Similar calculations are reported in Figure 4 (b) for the specific heat ${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right)$ for T between 0 ≤ T ≤ 1800 K at P = 0 (blue colored line) and P = 20 Gpa (red colored lines), respectively.

At P = 0 GPa, the RIM calculation of ${\mathsf{\Theta }}_{\mathrm{D}}$for zb BeO provided a value of ${\mathsf{\Theta }}_{\mathrm{D}}\left(\mathrm{T}\to 0\right)$ ~ 1390 K at nearly 0 K. It attained a minimum $\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}$ ${\mathsf{\Theta }}_{\mathrm{D}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\left(\mathrm{T}\right)~1150\mathrm{K}$at ~ 124 K and reached a higher value of ${\mathsf{\Theta }}_{\mathrm{D}}\left(297\right)$ ~ 1187 K at RT, respectively. In the absence of P-dependent measurements of ${\mathsf{\Theta }}_{\mathrm{D}}$, the study for zb BeO at P = 20 GPa predicted a slight decrease in ${\mathsf{\Theta }}_{\mathrm{D}}\left(\mathrm{T}\to 0\right)$ ~ 1370 K, while attaining ${\mathsf{\Theta }}_{\mathrm{D}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\left(\mathrm{T}\right)$ ∼ 1177 K at 93 K and achieving ${\mathsf{\Theta }}_{\mathrm{D}}\left(297\right)$~1291 K at RT, respectively. The simulated trends in zb BeO have agreed reasonably well with the behavior noticed in many compound semiconductors [69,70,71,72,73,74]. The rise of ${\mathsf{\Theta }}_{\mathrm{D}}\left(\mathrm{T}\right)$at T > RT

can be attributed to the P-dependent stiffening of high frequency ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}$, ${\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}$, and ${\mathsf{\omega }}_{\mathrm{L}\mathrm{O}}$ phonons. In Table 5, we have summarized our RIM results by comparing them with the limited theoretical and/or experimental data for ${\mathsf{\Theta }}_{\mathrm{D}}\left(0\right)$, ${\mathsf{\Theta }}_{\mathrm{D}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\left(\mathrm{T}\right)$, ${\mathsf{\Theta }}_{\mathrm{D}}\left(297\right),{\mathsf{\Theta }}_{\mathrm{D}}\left(\mathrm{H}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{T}\right)$in (K). Similar comparison of RIM calculations of ${\mathrm{C}}_{\mathrm{V}}$(100), ${\mathrm{C}}_{\mathrm{V}}$(297), ${\mathrm{C}}_{\mathrm{V}}$(High T) in (J/mol-K); and $\mathsf{\alpha }\left(\mathrm{T}\right)$ in (10^-6 K⁻¹) [6,80,81] with experimental data are also reported in Table 5.

3.2. Pressure Dependent Characteristics of zb BeO

By using the diamond anvil cell, many P-dependent Raman scattering studies are performed for comprehending the acoustical and optical phonons in several tetrahedrally coordinated elemental group IV and II-VI, III-V compound semiconductors [69,70,71,72,73,101]. Except for III-Ns and SiC, the existing data on most materials have confirmed [101] observing pressure-induced mode softening ${\mathsf{\gamma }}_{\mathrm{T}\mathrm{A}(\mathrm{X},\mathrm{L})}$of zone edge transverse acoustical${\mathsf{\omega }}_{\mathrm{T}\mathrm{A}(\mathrm{X},\mathrm{L})}$ $\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{s}.$ The P-induced stiffening is perceived, however, in every material for their zone edge longitudinal acoustical ${\mathsf{\gamma }}_{\mathrm{L}\mathrm{A}(\mathrm{X},\mathrm{L})};$ ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}(\mathrm{X},\mathrm{L})}$and high frequency optical ${\mathsf{\gamma }}_{\mathrm{L}\mathrm{O},\mathrm{T}\mathrm{O}(\mathrm{X},\mathrm{L})};{\mathsf{\omega }}_{\mathrm{L}\mathrm{O},\mathrm{T}\mathrm{O}(\mathrm{X},\mathrm{L})}$ phonons. Compression induced softening of ${\mathsf{\omega }}_{\mathrm{T}\mathrm{A}(\mathrm{X},\mathrm{L})}$phonons in II-VI materials and zb ZnO have triggered the negative tension in their bonds due to increased repulsion of electron-charge overlaps. Such mode softening is responsible for attributing negative thermal expansion coefficients $\mathsf{\alpha }\left(\mathrm{T}\right)$in many elemental and II-VI compound semiconductors [101] including zb ZnO [102].

3.2.1. Grüneisen Dispersions

$\mathrm{F}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{z}\mathrm{b}\mathrm{B}\mathrm{e}\mathrm{O}$we have reported in Table 4, our RIM results of the linear pressure coefficients ${\mathrm{a}}_{\mathrm{j}}^{\mathrm{P}}$(≡ $\partial {\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)/\partial \mathrm{P}$) and ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{f}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{B}\mathrm{Z}$. The methodology succinctly outlined in Section 2.3.1 is adopted to simulate $\overrightarrow{\mathrm{q}}$-dependent Grüneisen dispersions ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$for zb BeO. It can be noted that the Grüneisen constants ${\mathsf{\gamma }}_{\mathrm{T}\mathrm{A}(\mathrm{X},\mathrm{L})}$ for the zone edge ${\mathsf{\omega }}_{\mathrm{T}\mathrm{A}(\mathrm{X},\mathrm{L})}\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{s}$ revealed very small negative values (see Table 4) as compared to the large positive results noticed for high frequency longitudinal acoustic ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}$and optical ${\mathsf{\omega }}_{\mathrm{L}\mathrm{O}};{\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}$phonons. For zb BeO, the T-dependent simulations of linear thermal expansion coefficient $\mathsf{\alpha }\left(\mathrm{T}\right)$reported in Section 3.2.2 will be compared/contrasted and discussed with the existing experimental [80,81] and theoretical [82,83] data. In Figure 5, we have displayed our systematically simulated results of the mode Grüneisen dispersions ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ for zb BeO along the high symmetry directions.

The results of Grüneisen dispersions (see Figure 5) and the selected values listed in Table 4 at a few high critical points have attested some important facts that for the zb BeO material, the ${\mathsf{\gamma }}_{\mathrm{T}\mathrm{A}(\mathrm{X},\mathrm{L})}$ of zone edge ${\mathsf{\omega }}_{\mathrm{T}\mathrm{A}\left(\mathrm{X},\mathrm{L}\right)}$phonons exhibit much smaller negative values as compared to the large positive values of ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ for high frequency longitudinal acoustical ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}(\mathrm{X},\mathrm{L})}$and optical ${\mathsf{\omega }}_{\mathrm{L}\mathrm{O},\mathrm{T}\mathrm{O}(\mathrm{X},\mathrm{L})}$ modes. The impact of these$\overrightarrow{\mathrm{q}}$-dependent Grüneisen dispersions ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$on the linear thermal expansion coefficient $\mathsf{\alpha }\left(\mathrm{T}\right)$of zb BeO is reported in Section 3.2.2.

3.2.2. Thermal Expansion Coefficient

In Section 3.3.1, we have described the RIM results of mode Grüneisen dispersions ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ for the zb BeO material. By carefully incorporating values of ${\mathsf{\gamma }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ and ${\mathrm{C}}_{\mathrm{v}}\left(\mathrm{T}\right)$ in Equations (7-9), we have simulated T-dependent linear thermal expansion coefficient $\mathsf{\alpha }\left(\mathrm{T}\right)$. For zb BeO the results are displayed in Figure 6 for 0 ≤ T ≤ 1900 K. Unlike zb ZnO which showed negative thermal expansion [102] at low temperature, the RIM study of zb BeO has exhibited positive $\mathsf{\alpha }\left(\mathrm{T}\right)$ in the entire T range. Theoretical study has provided (green color line) a reasonably good comparison with the existing experimental data of Slack and Bartram [80] (open red square) and Kozlovskii et al. [81] (open inverted black triangle). Obviously, the RIM simulations have established an important fact that the sum of positive high frequency Grüneisen $\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{s}({\mathsf{\gamma }}_{\mathrm{L}\mathrm{O},\mathrm{T}\mathrm{O}(\mathrm{X},\mathrm{L})};$ ${\mathsf{\omega }}_{\mathrm{L}\mathrm{O},\mathrm{T}\mathrm{O}(\mathrm{X},\mathrm{L})})$ in Equation (7) are dominated over the very small negative values of ${\mathsf{\gamma }}_{\mathrm{T}\mathrm{A}(\mathrm{X},\mathrm{L})}$. Thus, at low T our study for zb BeO has offered no negative (see Figure 6) values of $\mathsf{\alpha }\left(\mathrm{T}\right)$.

4. Concluding Remarks

Understanding the electronic, structural, and thermal properties of a solid requires precise knowledge of its vibrational characteristics [77]. In the absence of INS data for zb XOs, theoretical simulations of lattice dynamics

and thermo dynamical properties have played crucial roles in solid state physics, materials science, and electronics. Among the II-Os, the wide bandgap BeO material has offered an exceptional thermal conductivity, chemical stability and high melting temperature. BeO is considered valuable for designing optoelectronics, power electronics, and thermal management systems. To develop LDHs-based advanced flexible micro-, nano electronics, the knowledge of phonon dispersion ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$ and thermal properties [e.g., ${‌}_{\mathrm{D}}\left(\mathrm{T}\right),{\mathrm{C}}_{\mathrm{V}}\left(\mathrm{T}\right),\mathsf{\kappa }\left(\mathrm{T}\right),$ γ(T), α(T)] of zb BeO is crucial. By adopting a realistic RIM, we have reported results of our methodical studies at ambient and high-pressure P to comprehend its structural, phonon, and thermodynamical properties. At ambient P = 0 GPa, the study has provided accurate ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$in good agreement with the first-principles calculations [66]. At higher P = 20 GPa, the predictions are made for the ${\mathsf{\omega }}_{\mathrm{j}}\left(\overrightarrow{\mathrm{q}}\right)$, as well as T-dependent ${‌}_{\mathrm{D}}\left(\mathrm{T}\right)$and ${\mathrm{C}}_{\mathrm{V}}\left(\mathrm{T}\right)$.

Unlike zb ZnO, our RIM study has exhibited atypical phonon and thermodynamical properties for the zb BeO material. In zb BeO, the ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}$mode frequencies along Γ-X and Γ-L directions are higher than the ${\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}$phonons. Vibrational characteristics in zb ZnO are analogous to the conventional II-VI compound semiconductors. Like zb BeO, unusual trends in the phonon characteristics are also noticed in zb BN and diamond C materials [96,97,98,99,100]. These features have provided strong corroboration to our earlier assertion that the structural, elastic and vibrational traits of zb BeO exhibit similarities to the super hard BN and C materials [71]. In zb BeO, phonons are the major heat carriers. In ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}$ modes, the Be and O atoms vibrate in the same direction as the wave propagation. However, in the ${\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}$ phonons, the atoms move perpendicular to the wave propagation. As ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}>{\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}$, one expects ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}$modes to carry more energy than the ${\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}$phonons. This could lead to increased scattering rates and possibly lower the thermal conductivity of zb BeO. The interplay between ${\mathsf{\omega }}_{\mathrm{L}\mathrm{A}}$ and ${\mathsf{\omega }}_{\mathrm{T}\mathrm{O}}$modes and their respective scattering mechanisms is a challenging issue and requires further studies for accurately predicting and understanding thermal transport properties in zb BeO.

As compared to zb ZnO, the zone edge Grüneisen constants ${\mathsf{\gamma }}_{\mathrm{T}\mathrm{A}(\mathrm{X},\mathrm{L})}$ for the ${\mathsf{\omega }}_{\mathrm{T}\mathrm{A}(\mathrm{X},\mathrm{L})}$ $\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{s}$ in zb BeO material have revealed very small negative values. Thus, our systematic calculations in the zb ZnO have exhibited negative $\mathsf{\alpha }\left(\mathrm{T}\right)$ at low T while positive results are predicted for zb BeO material in the entire 0 ≤ T ≤ 1900 K range. In semiconductors, the mode-softening of ${\mathsf{\gamma }}_{\mathrm{T}\mathrm{A}(\mathrm{X},\mathrm{L})}$triggers negative tension in their atomic bonding due to increased repulsion of electron-charge overlaps to cause phase transitions. Both the zb BeO, ZnO and their Be_xZn_1-xO alloys [103] are known to play important roles for managing heat dissipation in high-power electronics, optoelectronic devices, and flexible nano modules operating at varying temperatures which specifically require precise dimensional stability. Controlling and designing novel materials with a specific value of $\mathsf{\alpha }\left(\mathrm{T}\right)$ is essential for managing heat especially in miniaturized high-density devices are needed in the high-speed integrated electronics. By combining zb BeO and ZnO materials with positive and negative thermal expansion, it is quite likely to create LDH-based structures with tailored $\mathsf{\alpha }\left(\mathrm{T}\right)$ even for achieving near-zero value.